EMAN2
runsirt_Cart.cpp
Go to the documentation of this file.
00001 #include "mpi.h"
00002 
00003 #include "emdata.h"
00004 
00005 #include "sirt_Cart.h"
00006 #include "utilcomm_Cart.h"
00007 #include "project3d.h"
00008 #include "HyBR_Cart.h"
00009 
00010 #define PI 3.14159265358979
00011 #define ROW 0
00012 #define COL 1
00013 
00014 using namespace EMAN;
00015 
00016 // Reconstruct the 3-D density of a macromolecule from
00017 // a stack of 2-D images using the SIRT algorithm
00018 //
00019 // MPI Parallel Implementation using Cartesian virtual topology
00020 
00021 int main(int argc, char ** argv)
00022 {
00023    MPI_Comm comm = MPI_COMM_WORLD;
00024    int ncpus, mypid, ierr, mpierr=0;
00025    int nloc, maxit=0, method=0; 
00026    float lam = 0.0, tol = 0.0;
00027    char symstr[10]="none";
00028    std::string symmetry;
00029    double t0;
00030    FILE *fp;
00031    MPI_Comm comm_2d, comm_row, comm_col;
00032 
00033    MPI_Status mpistatus;
00034    MPI_Init(&argc,&argv);
00035    MPI_Comm_size(comm,&ncpus);
00036    MPI_Comm_rank(comm,&mypid);
00037 
00038    int dims[2], periods[2], my2dpid, mycoords[2];
00039    int srpid, srcoords[2], keep_dims[2];
00040 
00041    char  stackfname[100],voutfname[100], paramfname[100];
00042    EMData **expimages;
00043 
00044    // parse the command line and set filenames  
00045    if (argc < 5) {
00046      if (mypid == 0) {
00047          printf("Not enough arguments to the command...\n");
00048          printf("Usage: runsirt -data=<imagestack> ");
00049          printf("-param=<parameter file that contains angles and shifts> "); 
00050          printf("-out=<output filename base string> ");
00051          printf("-maxit=<maximum number of iterations> ");
00052          printf("-lam=<lambda> ");
00053          printf("-tol=<convergence tolerance> ");
00054          printf("-sym=<symmtry> \n");
00055          printf("-rowdim=<row dimension of Cartesian topology> ");
00056          printf("-coldim=<column dimension of Cartesian topology> ");
00057          printf("-method=<method to use: 0 SIRT (default), 1 Hybrid>\n");
00058      }
00059      ierr = MPI_Finalize();
00060      exit(1);
00061    }
00062    int ia=1;
00063    while (ia < argc) {
00064       if ( !strncmp(argv[ia],"-data",5) ) {
00065          strcpy(stackfname,&argv[ia][6]);
00066       }
00067       else if ( !strncmp(argv[ia],"-param",6) ) {
00068          strcpy(paramfname,&argv[ia][7]);
00069       }
00070       else if ( !strncmp(argv[ia],"-out",4) ) {
00071          strcpy(voutfname,&argv[ia][5]);
00072       }
00073       else if ( !strncmp(argv[ia],"-rowdim",7) ) {
00074          dims[ROW] = atoi(&argv[ia][8]); // Row dimension of the topology
00075       }
00076       else if ( !strncmp(argv[ia],"-coldim",7) ) {
00077          dims[COL] = atoi(&argv[ia][8]); // Column dimension of the topology
00078       }
00079       else if ( !strncmp(argv[ia],"-maxit",6) ) {
00080          maxit = atoi(&argv[ia][7]);
00081       }
00082       else if ( !strncmp(argv[ia],"-lam",4) ) {
00083          lam = atof(&argv[ia][5]);
00084       }
00085       else if ( !strncmp(argv[ia],"-tol",4) ) {
00086          tol = atof(&argv[ia][5]);
00087       }
00088       else if ( !strncmp(argv[ia],"-sym",4) ) {
00089          strcpy(symstr, &argv[ia][5]);
00090       }
00091       else if ( !strncmp(argv[ia],"-method",7) ) {
00092          method = atoi(&argv[ia][8]);
00093       }
00094       else {
00095          if (mypid ==0) printf("invalid option: %s\n", argv[ia]);
00096          ierr = MPI_Finalize();
00097          exit(1);
00098       }
00099       ia++;
00100    }
00101 
00102   if (dims[ROW]*dims[COL] != ncpus){
00103         printf("ERROR: rowdim*coldim != ncpus\n");
00104         ierr = MPI_Finalize();
00105         return -1;
00106   }
00107 
00108   // Set up the Cartesian virtual topology: comm_2d
00109   periods[ROW] = periods[COL] = 1; // Set the periods for wrap-around
00110   MPI_Cart_create(comm, 2, dims, periods, 1, &comm_2d);
00111   MPI_Comm_rank(comm_2d, &my2dpid); //Get my pid in the new 2D topology
00112   MPI_Cart_coords(comm_2d, my2dpid, 2, mycoords); // Get my coordinates
00113   
00114   // printf("MPI_2d: mypid = %d, my2dpid = %d, mycoords = [%d, %d] \n", 
00115   //mypid, my2dpid, mycoords[ROW], mycoords[COL]);
00116 
00117   /* Create the row-based sub-topology */ 
00118   keep_dims[ROW] = 0; 
00119   keep_dims[COL] = 1; 
00120   MPI_Cart_sub(comm_2d, keep_dims, &comm_row); 
00121 
00122   /* Create the column-based sub-topology */ 
00123   keep_dims[ROW] = 1; 
00124   keep_dims[COL] = 0; 
00125   MPI_Cart_sub(comm_2d, keep_dims, &comm_col); 
00126 
00127    // read and distribute a stack of experimental images along row processors
00128    t0 = MPI_Wtime();
00129    ierr = ReadStackandDist_Cart(comm_2d, &expimages, stackfname, &nloc);
00130    if (ierr == 0) {
00131         if (mypid == 0) {
00132            printf("Finished reading and distributing image stack onto ");
00133            printf("Cartesian topology\n");
00134            printf("I/O time for reading image stack = %11.3e\n",
00135                   MPI_Wtime() - t0);
00136         }
00137    }
00138    else {
00139       if (mypid == 0) 
00140          printf("Failed to read the image stack %s! exit...\n",stackfname);
00141       ierr = MPI_Finalize();
00142       exit(1);
00143    }
00144 
00145    // make a copy of the images for removing the background; 
00146    // this stack will be used for reconstruction
00147    EMData** cleanimages = new EMData*[nloc];
00148    for ( int i = 0 ; i < nloc ; ++i) {
00149         cleanimages[i] = expimages[i]->copy();
00150    }
00151 
00152    int nx = cleanimages[0]->get_xsize();
00153    int ny = cleanimages[0]->get_ysize();
00154    int nz = nx;
00155 
00156    Vec3i volsize;
00157    Vec3i origin;
00158    volsize[0] = nx;
00159    volsize[1] = ny;
00160    volsize[2] = nz;
00161    origin[0] = nx/2 + 1;
00162    origin[1] = ny/2 + 1;
00163    origin[2] = nz/2 + 1;
00164    int ri = nx/2 - 2;
00165 
00166    ierr = CleanStack_Cart(comm_col, cleanimages, nloc, ri, volsize, origin);
00167 
00168    // read angle and shift data and distribute along first column
00169    float * angleshift = new float[5*nloc];
00170 
00171    ierr = ReadAngTrandDist_Cart(comm_2d, comm_row, dims, angleshift, 
00172                                 paramfname, nloc);
00173    if (ierr!=0) { 
00174       mpierr = MPI_Finalize();
00175       return 1;
00176    }
00177 
00178    // Use xvol to hold reconstructed volume
00179    EMData * xvol = new EMData();
00180 
00181    // set SIRT parameters
00182   
00183    if (maxit==0) maxit = 10;
00184    if (lam == 0.0) lam = 5.0e-6;
00185    if (tol == 0.0) tol = 1.0e-3;
00186    if ( !strncmp(symstr,"none",4) ) {
00187       symmetry = string("c1");
00188    }
00189    else {
00190       symmetry = string(symstr);
00191    }
00192    // write out all options 
00193    if (mypid == 0) {
00194       printf("SIRT options used:\n");
00195       printf("maxit  = %d\n", maxit);
00196       printf("lam    = %11.3e\n", lam);
00197       printf("tol    = %11.3e\n", tol);
00198       printf("sym    = %s\n", symmetry.c_str()); 
00199       printf("rowdim = %d\n", dims[ROW]); 
00200       printf("coldim = %d\n", dims[COL]); 
00201    }
00202 
00203    if (method!=1) method = 0; //use SIRT unless method=1
00204 
00205 
00206    // call SIRT to reconstruct
00207    t0 = MPI_Wtime();
00208    if (method == 0) {
00209       recons3d_sirt_mpi_Cart(comm_2d, comm_row, comm_col, cleanimages, 
00210                              angleshift, xvol, nloc, ri, lam, maxit, 
00211                              symmetry, tol);
00212    }
00213    else {
00214       int insolve = 1;
00215       recons3d_HyBR_mpi_Cart(comm_2d, comm_row, comm_col, cleanimages, 
00216                              angleshift, xvol, nloc, ri, maxit, symmetry, insolve);
00217    }
00218 
00219    if ( my2dpid == 0 ) 
00220        printf("Done with reconstruction\n");
00221 
00222    // write the reconstructed volume to disk
00223    EMUtil::ImageType WRITE_SPI = EMUtil::IMAGE_SINGLE_SPIDER;
00224    if ( my2dpid == 0 ) {
00225         xvol->write_image(voutfname, 0, WRITE_SPI);
00226    }
00227 
00228    // cleanup
00229    for ( int i = 0 ; i < nloc; ++i ) {
00230        EMDeletePtr(expimages[i]);
00231    }
00232    EMDeleteArray(expimages);
00233    for ( int i = 0 ; i < nloc; ++i ) {
00234        EMDeletePtr(cleanimages[i]);
00235    }
00236    EMDeleteArray(cleanimages);
00237 
00238    EMDeletePtr(xvol);
00239    EMDeleteArray(angleshift);
00240 
00241    ierr = MPI_Finalize();
00242 
00243    return 0; // main
00244 }