EMAN2
runali3d_d.cpp
Go to the documentation of this file.
00001 #include "mpi.h"
00002 #include "stdlib.h"
00003 
00004 // include EMAN2
00005 #include "emdata.h"
00006 #include "emassert.h"
00007 
00008 #include "ali3d_d_mpi.h"
00009 #include "alignoptions.h"
00010 #include "utilcomm.h"
00011 
00012 int main(int argc, char *argv[])
00013 {
00014     MPI_Comm comm = MPI_COMM_WORLD;
00015     int ncpus, mypid, ierr;
00016     int nloc; 
00017     double t0;
00018 
00019     MPI_Status mpistatus;
00020     MPI_Init(&argc,&argv);
00021     MPI_Comm_size(comm,&ncpus);
00022     MPI_Comm_rank(comm,&mypid);
00023     printf("mypid = %d, ncpus = %d\n", mypid, ncpus);
00024     char  volfname[100], optionsfname[100], stackfname[100],voutfname[100];
00025     voutfname[0] = 0; // default to empty string
00026     optionsfname[0] = 0;
00027     EMData **expimages;
00028 
00029     // parse the command line and set filenames 
00030     if (argc < 3) {
00031       if (mypid == 0) {
00032           printf("Not enough arguments to the command...\n");
00033           printf("Usage: runali3d_d -data=<imagestack> ");
00034           printf("-model=<initial 3D volume> ");
00035           printf("-out=<output filename base string> "); 
00036           printf("-options=<options filename>\n");
00037       }
00038       ierr = MPI_Finalize();
00039       exit(1);
00040     }
00041     int ia=0;
00042     while (ia < argc) {
00043        if ( !strncmp(argv[ia],"-data",5) ) {
00044           strcpy(stackfname,&argv[ia][6]);
00045        }
00046        else if ( !strncmp(argv[ia],"-model",6) ) {
00047           strcpy(volfname,&argv[ia][7]);
00048        }
00049        else if ( !strncmp(argv[ia],"-out",4) ) {
00050           strcpy(voutfname,&argv[ia][5]);
00051        }
00052        else if ( !strncmp(argv[ia],"-options",8) ) {
00053            strcpy(optionsfname,&argv[ia][9]);
00054        }
00055        ia++;
00056     }
00057 
00058     // read and broadcast the initial model
00059     t0 = MPI_Wtime();
00060     EMData *volume = new EMData();
00061     ierr = ReadVandBcast(comm, volume, volfname);
00062     if (ierr == 0) {
00063        if (mypid == 0) {
00064           printf("Finished reading and replicating volume\n");
00065           printf("I/O time for reading volume = %11.3e\n",
00066                  MPI_Wtime() - t0);
00067        }
00068     }
00069     else {
00070        if (mypid == 0) 
00071           printf("Failed to read the model volume %s! exit...\n", volfname);
00072        ierr = MPI_Finalize();
00073        exit(1);
00074     }
00075     
00076     // read and distribute a stack of experimental images
00077     t0 = MPI_Wtime();
00078     ierr = ReadStackandDist(comm, &expimages, stackfname, &nloc);
00079     if (ierr == 0) {
00080         if (mypid == 0) {
00081            printf("Finished reading and distributing image stack\n");
00082            printf("I/O time for reading image stack = %11.3e\n",
00083                   MPI_Wtime() - t0);
00084         }
00085     }
00086     else {
00087        if (mypid == 0) 
00088           printf("Failed to read the image stack %s! exit...\n",stackfname);
00089        ierr = MPI_Finalize();
00090        exit(1);
00091     }
00092 
00093     Vec3i volsize;
00094     Vec3i origin;
00095     volsize[0] = volume->get_xsize();
00096     volsize[1] = volume->get_ysize();
00097     volsize[2] = volume->get_zsize();
00098     origin[0] = volume->get_xsize()/2 + 1;
00099     origin[1] = volume->get_ysize()/2 + 1;
00100     origin[2] = volume->get_zsize()/2 + 1;
00101     int    ri = volume->get_xsize()/2 - 2;
00102 
00103     // make a copy of the images for removing the background; this stack will be used for reconstruction
00104     EMData** cleanimages = new EMData*[nloc];
00105     for ( int i = 0 ; i < nloc ; ++i) {
00106         cleanimages[i] = expimages[i]->copy();
00107     }
00108     ierr = CleanStack(comm, cleanimages, nloc, ri, volsize, origin);
00109 
00110     float *angleshift = new float[5*nloc];
00111     int maxiter = 1; // need to be able to set this, either on command line or in AlignOptions object
00112     
00113     EMData * mask3D = NULL;
00114     AlignOptions options(volsize);
00115     if ( optionsfname[0] ){
00116         ParseAlignOptions(comm, options, optionsfname, volsize[0]*volsize[1]*volsize[2], mask3D);
00117     } else {
00118         if ( mypid == 0 ) printf("Using default alignment options\n");
00119     }
00120 
00121     options.set_have_angles(false);
00122 
00123     try {
00124        ali3d_d(comm, volume, expimages, cleanimages, angleshift, nloc, options, voutfname);
00125     }
00126     catch (std::exception const& e) {
00127         printf("%s\n", e.what());
00128     }
00129     EMDeletePtr(volume);
00130     for ( int i = 0 ; i < nloc; ++i ) {
00131         EMDeletePtr(expimages[i]);
00132     }
00133     EMDeleteArray(expimages);
00134     for ( int i = 0 ; i < nloc; ++i ) {
00135         EMDeletePtr(cleanimages[i]);
00136     }
00137     EMDeleteArray(cleanimages);
00138     EMDeleteArray(angleshift);
00139     if ( mask3D != NULL ) {
00140         EMDeletePtr(mask3D);
00141     }
00142     ierr = MPI_Finalize();
00143     return 0;
00144 }
00145