00001 #include "mpi.h"
00002
00003 #include "emdata.h"
00004
00005 #include "project3d.h"
00006 #include "sirt.h"
00007 #include "utilcomm.h"
00008
00009 #define PI 3.14159265358979
00010 using namespace EMAN;
00011
00012
00013
00014
00015 int main(int argc, char ** argv)
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
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
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
00098
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
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
00132 EMData * xvol = new EMData();
00133
00134
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
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
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
00162 EMUtil::ImageType WRITE_SPI = EMUtil::IMAGE_SINGLE_SPIDER;
00163 if ( mypid == 0 ) {
00164 xvol->write_image(voutfname, 0, WRITE_SPI);
00165 }
00166
00167
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;
00183 }