00001 #include "mpi.h"
00002 #include "stdlib.h"
00003
00004
00005 #include "emdata.h"
00006
00007 #include "ali3d_d_mpi.h"
00008 #include "ali3d_unified_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, mpierr=0;
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], paramfname[100], stackfname[100],voutfname[100];
00025 EMData **expimages;
00026 int maxiter=0;
00027
00028
00029 if (argc < 4) {
00030 if (mypid == 0) {
00031 printf("Not enough arguments to the command...\n");
00032 printf("Usage: rununified -data=<imagestack> ");
00033 printf("-model=<initial 3D volume filename> ");
00034 printf("-param=<initial angles&shifts> ");
00035 printf("-out=<output volume filename>\n");
00036 printf("[-maxit=<maximum number of iterations>]\n");
00037 printf("[-sym=<symmtry type>]\n");
00038 printf("[-CTF]\n");
00039 }
00040 ierr = MPI_Finalize();
00041 exit(1);
00042 }
00043
00044 int ia=1;
00045 while (ia < argc) {
00046 if ( !strncmp(argv[ia],"-data",5) ) {
00047 strcpy(stackfname,&argv[ia][6]);
00048 if (mypid == 0)
00049 printf("image stack = %s\n", stackfname);
00050 }
00051 else if ( !strncmp(argv[ia],"-model",6) ) {
00052 strcpy(volfname,&argv[ia][7]);
00053 if (mypid == 0)
00054 printf("initial model = %s\n", volfname);
00055 }
00056 else if ( !strncmp(argv[ia],"-param",6) ) {
00057 strcpy(paramfname,&argv[ia][7]);
00058 if (mypid == 0)
00059 printf("initial parameters = %s\n", paramfname);
00060 }
00061 else if ( !strncmp(argv[ia],"-maxit",6) ) {
00062 maxiter = atoi(&(argv[ia][7]));
00063 if (mypid == 0)
00064 printf("maxiter = %d\n", maxiter);
00065 }
00066 else if ( !strncmp(argv[ia],"-out",4) ) {
00067 strcpy(voutfname,&argv[ia][5]);
00068 if (mypid == 0)
00069 printf("output volume file = %s\n", voutfname);
00070 }
00071 ia++;
00072 }
00073
00074
00075 t0 = MPI_Wtime();
00076 EMData *volume = new EMData();
00077 ierr = ReadVandBcast(comm, volume, volfname);
00078 if (ierr == 0) {
00079 if (mypid == 0) {
00080 printf("Finished reading and replicating volume\n");
00081 printf("I/O time for reading volume = %11.3e\n",
00082 MPI_Wtime() - t0);
00083 }
00084 }
00085 else {
00086 if (mypid == 0) printf("Failed to read the model volume %s! exit...\n", volfname);
00087 ierr = MPI_Finalize();
00088 exit(1);
00089 }
00090
00091
00092 t0 = MPI_Wtime();
00093 ierr = ReadStackandDist(comm, &expimages, stackfname,&nloc);
00094 if (ierr == 0) {
00095 if (mypid == 0) {
00096 printf("Finished reading and distributing image stack\n");
00097 printf("I/O time for reading image stack = %11.3e\n",
00098 MPI_Wtime() - t0);
00099 }
00100 }
00101 else {
00102 if (mypid == 0) printf("Failed to read the image stack %s! exit...\n", stackfname);
00103 ierr = MPI_Finalize();
00104 exit(1);
00105 }
00106
00107 Vec3i volsize;
00108 Vec3i origin;
00109 volsize[0] = volume->get_xsize();
00110 volsize[1] = volume->get_ysize();
00111 volsize[2] = volume->get_zsize();
00112 origin[0] = volume->get_xsize()/2 + 1;
00113 origin[1] = volume->get_ysize()/2 + 1;
00114 origin[2] = volume->get_zsize()/2 + 1;
00115 int ri = volume->get_xsize()/2 - 2;
00116 ierr = CleanStack(comm, expimages, nloc, ri, volsize, origin);
00117
00118
00119 float * angleshift = new float[5*nloc];
00120 ierr = ReadAngTrandDist(comm, angleshift, paramfname, nloc);
00121 if (ierr!=0) {
00122 mpierr = MPI_Finalize();
00123 return 1;
00124 }
00125
00126 if (maxiter <= 0) maxiter = 10;
00127 try {
00128 unified(comm, volume, expimages, angleshift, nloc,
00129 maxiter, voutfname);
00130 }
00131 catch (std::exception const& e) {
00132 printf("%s\n", e.what());
00133 }
00134
00135 #ifdef DEBUG
00136 for (int i = 0; i < nloc; i++)
00137 printf("%11.3e %11.3e %11.3e %11.3e %11.3e\n",
00138 angleshift[5*i+0],
00139 angleshift[5*i+1],
00140 angleshift[5*i+2],
00141 angleshift[5*i+3],
00142 angleshift[5*i+4]);
00143 #endif
00144
00145 EMDeletePtr(volume);
00146 for ( int i = 0 ; i < nloc; ++i ) {
00147 EMDeletePtr(expimages[i]);
00148 }
00149 EMDeleteArray(expimages);
00150 EMDeleteArray(angleshift);
00151 ierr = MPI_Finalize();
00152 return 0;
00153 }
00154