00001 #include "mpi.h"
00002
00003 #include "emdata.h"
00004
00005 #include "sirt_Cart.h"
00006 #include "utilcomm_Cart.h"
00007 #include "project3d.h"
00008 #include "HyBR_Cart.h"
00009
00010 #define PI 3.14159265358979
00011 #define ROW 0
00012 #define COL 1
00013
00014 using namespace EMAN;
00015
00016
00017
00018
00019
00020
00021 int main(int argc, char ** argv)
00022 {
00023 MPI_Comm comm = MPI_COMM_WORLD;
00024 int ncpus, mypid, ierr, mpierr=0;
00025 int nloc, maxit=0, method=0;
00026 float lam = 0.0, tol = 0.0;
00027 char symstr[10]="none";
00028 std::string symmetry;
00029 double t0;
00030 FILE *fp;
00031 MPI_Comm comm_2d, comm_row, comm_col;
00032
00033 MPI_Status mpistatus;
00034 MPI_Init(&argc,&argv);
00035 MPI_Comm_size(comm,&ncpus);
00036 MPI_Comm_rank(comm,&mypid);
00037
00038 int dims[2], periods[2], my2dpid, mycoords[2];
00039 int srpid, srcoords[2], keep_dims[2];
00040
00041 char stackfname[100],voutfname[100], paramfname[100];
00042 EMData **expimages;
00043
00044
00045 if (argc < 5) {
00046 if (mypid == 0) {
00047 printf("Not enough arguments to the command...\n");
00048 printf("Usage: runsirt -data=<imagestack> ");
00049 printf("-param=<parameter file that contains angles and shifts> ");
00050 printf("-out=<output filename base string> ");
00051 printf("-maxit=<maximum number of iterations> ");
00052 printf("-lam=<lambda> ");
00053 printf("-tol=<convergence tolerance> ");
00054 printf("-sym=<symmtry> \n");
00055 printf("-rowdim=<row dimension of Cartesian topology> ");
00056 printf("-coldim=<column dimension of Cartesian topology> ");
00057 printf("-method=<method to use: 0 SIRT (default), 1 Hybrid>\n");
00058 }
00059 ierr = MPI_Finalize();
00060 exit(1);
00061 }
00062 int ia=1;
00063 while (ia < argc) {
00064 if ( !strncmp(argv[ia],"-data",5) ) {
00065 strcpy(stackfname,&argv[ia][6]);
00066 }
00067 else if ( !strncmp(argv[ia],"-param",6) ) {
00068 strcpy(paramfname,&argv[ia][7]);
00069 }
00070 else if ( !strncmp(argv[ia],"-out",4) ) {
00071 strcpy(voutfname,&argv[ia][5]);
00072 }
00073 else if ( !strncmp(argv[ia],"-rowdim",7) ) {
00074 dims[ROW] = atoi(&argv[ia][8]);
00075 }
00076 else if ( !strncmp(argv[ia],"-coldim",7) ) {
00077 dims[COL] = atoi(&argv[ia][8]);
00078 }
00079 else if ( !strncmp(argv[ia],"-maxit",6) ) {
00080 maxit = atoi(&argv[ia][7]);
00081 }
00082 else if ( !strncmp(argv[ia],"-lam",4) ) {
00083 lam = atof(&argv[ia][5]);
00084 }
00085 else if ( !strncmp(argv[ia],"-tol",4) ) {
00086 tol = atof(&argv[ia][5]);
00087 }
00088 else if ( !strncmp(argv[ia],"-sym",4) ) {
00089 strcpy(symstr, &argv[ia][5]);
00090 }
00091 else if ( !strncmp(argv[ia],"-method",7) ) {
00092 method = atoi(&argv[ia][8]);
00093 }
00094 else {
00095 if (mypid ==0) printf("invalid option: %s\n", argv[ia]);
00096 ierr = MPI_Finalize();
00097 exit(1);
00098 }
00099 ia++;
00100 }
00101
00102 if (dims[ROW]*dims[COL] != ncpus){
00103 printf("ERROR: rowdim*coldim != ncpus\n");
00104 ierr = MPI_Finalize();
00105 return -1;
00106 }
00107
00108
00109 periods[ROW] = periods[COL] = 1;
00110 MPI_Cart_create(comm, 2, dims, periods, 1, &comm_2d);
00111 MPI_Comm_rank(comm_2d, &my2dpid);
00112 MPI_Cart_coords(comm_2d, my2dpid, 2, mycoords);
00113
00114
00115
00116
00117
00118 keep_dims[ROW] = 0;
00119 keep_dims[COL] = 1;
00120 MPI_Cart_sub(comm_2d, keep_dims, &comm_row);
00121
00122
00123 keep_dims[ROW] = 1;
00124 keep_dims[COL] = 0;
00125 MPI_Cart_sub(comm_2d, keep_dims, &comm_col);
00126
00127
00128 t0 = MPI_Wtime();
00129 ierr = ReadStackandDist_Cart(comm_2d, &expimages, stackfname, &nloc);
00130 if (ierr == 0) {
00131 if (mypid == 0) {
00132 printf("Finished reading and distributing image stack onto ");
00133 printf("Cartesian topology\n");
00134 printf("I/O time for reading image stack = %11.3e\n",
00135 MPI_Wtime() - t0);
00136 }
00137 }
00138 else {
00139 if (mypid == 0)
00140 printf("Failed to read the image stack %s! exit...\n",stackfname);
00141 ierr = MPI_Finalize();
00142 exit(1);
00143 }
00144
00145
00146
00147 EMData** cleanimages = new EMData*[nloc];
00148 for ( int i = 0 ; i < nloc ; ++i) {
00149 cleanimages[i] = expimages[i]->copy();
00150 }
00151
00152 int nx = cleanimages[0]->get_xsize();
00153 int ny = cleanimages[0]->get_ysize();
00154 int nz = nx;
00155
00156 Vec3i volsize;
00157 Vec3i origin;
00158 volsize[0] = nx;
00159 volsize[1] = ny;
00160 volsize[2] = nz;
00161 origin[0] = nx/2 + 1;
00162 origin[1] = ny/2 + 1;
00163 origin[2] = nz/2 + 1;
00164 int ri = nx/2 - 2;
00165
00166 ierr = CleanStack_Cart(comm_col, cleanimages, nloc, ri, volsize, origin);
00167
00168
00169 float * angleshift = new float[5*nloc];
00170
00171 ierr = ReadAngTrandDist_Cart(comm_2d, comm_row, dims, angleshift,
00172 paramfname, nloc);
00173 if (ierr!=0) {
00174 mpierr = MPI_Finalize();
00175 return 1;
00176 }
00177
00178
00179 EMData * xvol = new EMData();
00180
00181
00182
00183 if (maxit==0) maxit = 10;
00184 if (lam == 0.0) lam = 5.0e-6;
00185 if (tol == 0.0) tol = 1.0e-3;
00186 if ( !strncmp(symstr,"none",4) ) {
00187 symmetry = string("c1");
00188 }
00189 else {
00190 symmetry = string(symstr);
00191 }
00192
00193 if (mypid == 0) {
00194 printf("SIRT options used:\n");
00195 printf("maxit = %d\n", maxit);
00196 printf("lam = %11.3e\n", lam);
00197 printf("tol = %11.3e\n", tol);
00198 printf("sym = %s\n", symmetry.c_str());
00199 printf("rowdim = %d\n", dims[ROW]);
00200 printf("coldim = %d\n", dims[COL]);
00201 }
00202
00203 if (method!=1) method = 0;
00204
00205
00206
00207 t0 = MPI_Wtime();
00208 if (method == 0) {
00209 recons3d_sirt_mpi_Cart(comm_2d, comm_row, comm_col, cleanimages,
00210 angleshift, xvol, nloc, ri, lam, maxit,
00211 symmetry, tol);
00212 }
00213 else {
00214 int insolve = 1;
00215 recons3d_HyBR_mpi_Cart(comm_2d, comm_row, comm_col, cleanimages,
00216 angleshift, xvol, nloc, ri, maxit, symmetry, insolve);
00217 }
00218
00219 if ( my2dpid == 0 )
00220 printf("Done with reconstruction\n");
00221
00222
00223 EMUtil::ImageType WRITE_SPI = EMUtil::IMAGE_SINGLE_SPIDER;
00224 if ( my2dpid == 0 ) {
00225 xvol->write_image(voutfname, 0, WRITE_SPI);
00226 }
00227
00228
00229 for ( int i = 0 ; i < nloc; ++i ) {
00230 EMDeletePtr(expimages[i]);
00231 }
00232 EMDeleteArray(expimages);
00233 for ( int i = 0 ; i < nloc; ++i ) {
00234 EMDeletePtr(cleanimages[i]);
00235 }
00236 EMDeleteArray(cleanimages);
00237
00238 EMDeletePtr(xvol);
00239 EMDeleteArray(angleshift);
00240
00241 ierr = MPI_Finalize();
00242
00243 return 0;
00244 }