00001 #include "mpi.h"
00002 #include "stdlib.h"
00003
00004
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;
00028 optionsfname[0] = 0;
00029 EMData **expimages;
00030
00031
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
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
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
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
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
00138
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