Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

runsirt.cpp

Go to the documentation of this file.
00001 #include "mpi.h"
00002 
00003 #include "emdata.h"
00004 
00005 #include "project3d.h"
00006 #include "sirt.h"
00007 #include "utilcomm.h"
00008 
00009 #define PI 3.14159265358979
00010 using namespace EMAN;
00011 
00012 // Reconstruct the 3-D density of a macromolecule from
00013 // a stack of 2-D images using the SIRT algorithm
00014 
00015 int main(int argc, char ** argv)
00016 {
00017    MPI_Comm comm = MPI_COMM_WORLD;
00018    int ncpus, mypid, ierr, mpierr=0;
00019    int nloc, maxit = 0; 
00020    float lam = 0.0, tol = 0.0;
00021    char symstr[10]="none";
00022    std::string symmetry;
00023    double t0;
00024    FILE *fp;
00025 
00026    MPI_Status mpistatus;
00027    MPI_Init(&argc,&argv);
00028    MPI_Comm_size(comm,&ncpus);
00029    MPI_Comm_rank(comm,&mypid);
00030    printf("mypid = %d, ncpus = %d\n", mypid, ncpus);
00031 
00032    char  stackfname[100],voutfname[100], paramfname[100];
00033    EMData **expimages;
00034 
00035    // parse the command line and set filenames  
00036    if (argc < 3) {
00037      if (mypid == 0) {
00038          printf("Not enough arguments to the command...\n");
00039          printf("Usage: runsirt -data=<imagestack> ");
00040          printf("-param=<paramter file that contains angles & shifts> "); 
00041          printf("-out=<output filename base string> ");
00042          printf("-maxit=<max number of iterations> ");
00043          printf("-lam=<lambda> ");
00044          printf("-tol=<convergence tolerance> ");
00045          printf("-sym=<symmtry> \n");
00046      }
00047      ierr = MPI_Finalize();
00048      exit(1);
00049    }
00050    int ia=1;
00051    while (ia < argc) {
00052       if ( !strncmp(argv[ia],"-data",5) ) {
00053          strcpy(stackfname,&argv[ia][6]);
00054       }
00055       else if ( !strncmp(argv[ia],"-param",6) ) {
00056          strcpy(paramfname,&argv[ia][7]);
00057       }
00058       else if ( !strncmp(argv[ia],"-out",4) ) {
00059          strcpy(voutfname,&argv[ia][5]);
00060       }
00061       else if ( !strncmp(argv[ia],"-maxit",6) ) {
00062          maxit = atoi(&argv[ia][7]);
00063       }
00064       else if ( !strncmp(argv[ia],"-lam",4) ) {
00065          lam = atof(&argv[ia][5]);
00066       }
00067       else if ( !strncmp(argv[ia],"-tol",4) ) {
00068          tol = atof(&argv[ia][5]);
00069       }
00070       else if ( !strncmp(argv[ia],"-sym",4) ) {
00071          strcpy(symstr, &argv[ia][5]);
00072       }
00073       else {
00074          if (mypid ==0) printf("invalid option: %s\n", argv[ia]);
00075          ierr = MPI_Finalize();
00076          exit(1);
00077       }
00078       ia++;
00079    }
00080 
00081    // read and distribute a stack of experimental images
00082    t0 = MPI_Wtime();
00083    ierr = ReadStackandDist(comm, &expimages, stackfname, &nloc);
00084    if (ierr == 0) {
00085         if (mypid == 0) {
00086            printf("I/O time for reading image stack = %11.3e\n",
00087                   MPI_Wtime() - t0);
00088         }
00089    }
00090    else {
00091       if (mypid == 0) 
00092          printf("Failed to read the image stack %s! exit...\n",stackfname);
00093       ierr = MPI_Finalize();
00094       exit(1);
00095    }
00096 
00097    // make a copy of the images for removing the background; 
00098    // this stack will be used for reconstruction
00099    EMData** cleanimages = new EMData*[nloc];
00100    for ( int i = 0 ; i < nloc ; ++i) {
00101         cleanimages[i] = expimages[i]->copy();
00102    }
00103 
00104    int nx = cleanimages[0]->get_xsize();
00105    int ny = cleanimages[0]->get_ysize();
00106    int nz = nx;
00107 
00108    Vec3i volsize;
00109    Vec3i origin;
00110    volsize[0] = nx;
00111    volsize[1] = ny;
00112    volsize[2] = nz;
00113    origin[0] = nx/2 + 1;
00114    origin[1] = ny/2 + 1;
00115    origin[2] = nz/2 + 1;
00116    int ri = nx/2 - 2;
00117 
00118    ierr = CleanStack(comm, cleanimages, nloc, ri, volsize, origin);
00119    if (mypid == 0) printf("image size: nx = %d, ny = %d, ri = %d\n",
00120                           nx, ny, ri); 
00121 
00122 
00123    // read angle and shift data and distribute
00124    float * angleshift = new float[5*nloc];
00125    ierr = ReadAngTrandDist(comm, angleshift, paramfname, nloc);
00126    if (ierr!=0) { 
00127       mpierr = MPI_Finalize();
00128       return 1;
00129    }
00130 
00131    // Use xvol to hold reconstructed volume
00132    EMData * xvol = new EMData();
00133 
00134    // set SIRT parameters
00135    if (maxit == 0) maxit = 10;
00136    if (lam == 0.0) lam = 5.0e-6;
00137    if (tol == 0.0) tol = 1.0e-3;
00138    if ( !strncmp(symstr,"none",4) ) {
00139       symmetry = string("c1");
00140    }
00141    else {
00142       symmetry = string(symstr);
00143    }
00144    // write out all options 
00145    if (mypid == 0) {
00146       printf("SIRT options used:\n");
00147       printf("maxit = %d\n", maxit);
00148       printf("lam   = %11.3e\n", lam);
00149       printf("tol   = %11.3e\n", tol);
00150       printf("sym   = %s\n", symmetry.c_str()); 
00151    }
00152 
00153    // call SIRT to reconstruct
00154    recons3d_sirt_mpi(comm, cleanimages, angleshift, xvol, nloc, ri, 
00155                      lam, maxit, symmetry, tol);
00156 
00157    if ( mypid == 0 ) {
00158       printf("Done with SIRT\n");
00159    }
00160 
00161    // write the reconstructed volume to disk
00162    EMUtil::ImageType WRITE_SPI = EMUtil::IMAGE_SINGLE_SPIDER;
00163    if ( mypid == 0 ) {
00164         xvol->write_image(voutfname, 0, WRITE_SPI);
00165    }
00166 
00167    // cleanup
00168    for ( int i = 0 ; i < nloc; ++i ) {
00169        EMDeletePtr(expimages[i]);
00170    }
00171    EMDeleteArray(expimages);
00172    for ( int i = 0 ; i < nloc; ++i ) {
00173        EMDeletePtr(cleanimages[i]);
00174    }
00175    EMDeleteArray(cleanimages);
00176 
00177    EMDeletePtr(xvol);
00178    EMDeleteArray(angleshift);
00179 
00180    ierr = MPI_Finalize();
00181 
00182    return 0; // main
00183 }

Generated on Tue Jun 11 13:46:18 2013 for EMAN2 by  doxygen 1.3.9.1