#include "mpi.h"
#include "emdata.h"
#include "project3d.h"
#include "sirt.h"
#include "utilcomm.h"
Include dependency graph for runsirt.cpp:
Go to the source code of this file.
Defines | |
#define | PI 3.14159265358979 |
Functions | |
int | main (int argc, char **argv) |
|
Definition at line 9 of file runsirt.cpp. |
|
Definition at line 15 of file runsirt.cpp. References CleanStack(), EMAN::EMData::copy(), EMDeleteArray(), EMDeletePtr(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), ierr, nx, ny, ReadAngTrandDist(), ReadStackandDist(), recons3d_sirt_mpi(), EMAN::Vec3i, and EMAN::EMData::write_image(). 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 }
|