EMAN2
runhybrid.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 "ali3d_unified_mpi.h"
00010 #include "alignoptions.h"
00011 #include "utilparse.h"
00012 #include "utilcomm.h"
00013 
00014 int main(int argc, char *argv[])
00015 {
00016     MPI_Comm comm = MPI_COMM_WORLD;
00017     int ncpus, mypid, ierr;
00018     int nloc; 
00019     double t0;
00020 
00021     MPI_Status mpistatus;
00022     MPI_Init(&argc,&argv);
00023     MPI_Comm_size(comm,&ncpus);
00024     MPI_Comm_rank(comm,&mypid);
00025     printf("mypid = %d, ncpus = %d\n", mypid, ncpus);
00026     char  volfname[100], stackfname[100],voutfname[100],optionsfname[100];
00027     voutfname[0] = 0; // default to empty string
00028     optionsfname[0] = 0;
00029     EMData **expimages;
00030 
00031     // parse the command line and set filenames 
00032     if (argc < 3) {
00033       if (mypid == 0) {
00034           printf("Not enough arguments to the command...\n");
00035           printf("Usage: sxhrefine -data=<imagestack> ");
00036           printf("-model=<initial 3D volume> "); 
00037           printf("-out=<output filename base string> ");
00038           printf("-options=<options filename>\n");
00039       }
00040       ierr = MPI_Finalize();
00041       exit(1);
00042     }
00043     int ia=0;
00044     while (ia < argc) {
00045        if ( !strncmp(argv[ia],"-data",5) ) {
00046           strcpy(stackfname,&argv[ia][6]);
00047        }
00048        else if ( !strncmp(argv[ia],"-model",6) ) {
00049           strcpy(volfname,&argv[ia][7]);
00050        }
00051        else if ( !strncmp(argv[ia],"-out",4) ) {
00052           strcpy(voutfname,&argv[ia][5]);
00053        }
00054        else if ( !strncmp(argv[ia],"-options",8) ) {
00055            strcpy(optionsfname,&argv[ia][9]);
00056        }
00057        ia++;
00058     }
00059 
00060     // read and broadcast the initial model
00061     t0 = MPI_Wtime();
00062     EMData *volume = new EMData();
00063     ierr = ReadVandBcast(comm, volume, volfname);
00064     if (ierr == 0) {
00065         if (mypid == 0) {
00066            printf("Finished reading and replicating volume\n");
00067            printf("I/O time for reading volume = %11.3e\n",
00068                   MPI_Wtime() - t0);
00069         }
00070     }
00071     else {
00072        if (mypid == 0) printf("Failed to read the model volume %s! exit...\n", volfname);
00073        ierr = MPI_Finalize();
00074        exit(1);
00075     }
00076 
00077     // read and distribute a stack of experimental images
00078     t0 = MPI_Wtime();
00079     ierr = ReadStackandDist(comm, &expimages, stackfname, &nloc);
00080     if (ierr == 0) {
00081         if (mypid == 0) {
00082            printf("Finished reading and distributing image stack\n");
00083            printf("I/O time for reading image stack = %11.3e\n",
00084                   MPI_Wtime() - t0);
00085         }
00086     }
00087     else {
00088        if (mypid == 0) 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 
00111     float * angleshift = new float[5*nloc];
00112     
00113 
00114     // these hard-coded numbers need to be removed eventually
00115     int max_refine_cycle = 0;
00116     int max_iter_unified = 10;
00117 
00118     char out_fname[128];
00119 
00120     AlignOptions options(volsize);
00121 
00122 
00123     std::ifstream option_stream;
00124     std::string current_option;
00125     char option_buffer[100];
00126     int int_option;
00127     float float_option;
00128 
00129     EMData * mask3D = NULL;
00130     if ( optionsfname[0] ) {
00131         ParseAlignOptions(comm, options, optionsfname, volsize[0]*volsize[1]*volsize[2], mask3D);
00132     } else {
00133         if ( mypid == 0 ) printf("Using default alignment options\n");
00134     }
00135     options.set_have_angles(false);
00136 
00137     // the below is a bit of a kluge; maxit in the options file indicates total number of refinement cycles, 
00138     // options.maxit is number of projection matching steps within each refinement cycle
00139     max_refine_cycle = options.get_maxit();
00140     options.set_maxit(1);
00141 
00142     try {
00143         for ( int iter = 0; iter < max_refine_cycle ; ++iter ) {
00144             if (mypid == 0) printf("refinement cycle: %d\n", iter+1);
00145             sprintf(out_fname, "%smajor%d", voutfname, iter);
00146 
00147             ali3d_d(comm, volume, expimages, cleanimages, angleshift, nloc, options, out_fname);
00148 
00149             unified(comm, volume, cleanimages, angleshift, nloc, 
00150                     max_iter_unified, out_fname);
00151             options.set_have_angles(true);
00152         }
00153     }
00154     catch (std::exception const& e) {
00155         printf("%s\n", e.what());
00156     }
00157 
00158     EMDeletePtr(volume);
00159     for ( int i = 0 ; i < nloc; ++i ) {
00160         EMDeletePtr(expimages[i]);
00161     }
00162     EMDeleteArray(expimages);
00163     for ( int i = 0 ; i < nloc; ++i ) {
00164         EMDeletePtr(cleanimages[i]);
00165     }
00166     EMDeleteArray(cleanimages);
00167     EMDeleteArray(angleshift);
00168     if (mask3D != NULL) {
00169         EMDeletePtr(mask3D);
00170     }
00171     ierr = MPI_Finalize();
00172     return 0;
00173 }
00174