#include "mpi.h"
#include "emdata.h"
#include "sirt_Cart.h"
#include "utilcomm_Cart.h"
#include "project3d.h"
#include "HyBR_Cart.h"
Include dependency graph for runsirt_Cart.cpp:
Go to the source code of this file.
Defines | |
#define | PI 3.14159265358979 |
#define | ROW 0 |
#define | COL 1 |
Functions | |
int | main (int argc, char **argv) |
#define COL 1 |
Definition at line 12 of file runsirt_Cart.cpp.
Referenced by LBD_Cart(), main(), ReadAngTrandDist_Cart(), ReadStackandDist_Cart(), recons3d_CGLS_mpi_Cart(), recons3d_HyBR_mpi_Cart(), recons3d_sirt_mpi_Cart(), setpart_gc1(), setpart_gr1(), and sphpart().
#define PI 3.14159265358979 |
Definition at line 10 of file runsirt_Cart.cpp.
#define ROW 0 |
Definition at line 11 of file runsirt_Cart.cpp.
Referenced by LBD_Cart(), main(), ReadAngTrandDist_Cart(), ReadStackandDist_Cart(), recons3d_CGLS_mpi_Cart(), recons3d_HyBR_mpi_Cart(), recons3d_sirt_mpi_Cart(), setpart_gc1(), setpart_gr1(), and sphpart().
int main | ( | int | argc, | |
char ** | argv | |||
) |
Definition at line 21 of file runsirt_Cart.cpp.
References CleanStack_Cart(), COL, copy(), EMDeleteArray(), EMDeletePtr(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), ierr, nx, ny, ReadAngTrandDist_Cart(), ReadStackandDist_Cart(), recons3d_HyBR_mpi_Cart(), recons3d_sirt_mpi_Cart(), ROW, and EMAN::EMData::write_image().
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 // parse the command line and set filenames 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]); // Row dimension of the topology 00075 } 00076 else if ( !strncmp(argv[ia],"-coldim",7) ) { 00077 dims[COL] = atoi(&argv[ia][8]); // Column dimension of the topology 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 // Set up the Cartesian virtual topology: comm_2d 00109 periods[ROW] = periods[COL] = 1; // Set the periods for wrap-around 00110 MPI_Cart_create(comm, 2, dims, periods, 1, &comm_2d); 00111 MPI_Comm_rank(comm_2d, &my2dpid); //Get my pid in the new 2D topology 00112 MPI_Cart_coords(comm_2d, my2dpid, 2, mycoords); // Get my coordinates 00113 00114 // printf("MPI_2d: mypid = %d, my2dpid = %d, mycoords = [%d, %d] \n", 00115 //mypid, my2dpid, mycoords[ROW], mycoords[COL]); 00116 00117 /* Create the row-based sub-topology */ 00118 keep_dims[ROW] = 0; 00119 keep_dims[COL] = 1; 00120 MPI_Cart_sub(comm_2d, keep_dims, &comm_row); 00121 00122 /* Create the column-based sub-topology */ 00123 keep_dims[ROW] = 1; 00124 keep_dims[COL] = 0; 00125 MPI_Cart_sub(comm_2d, keep_dims, &comm_col); 00126 00127 // read and distribute a stack of experimental images along row processors 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 // make a copy of the images for removing the background; 00146 // this stack will be used for reconstruction 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 // read angle and shift data and distribute along first column 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 // Use xvol to hold reconstructed volume 00179 EMData * xvol = new EMData(); 00180 00181 // set SIRT parameters 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 // write out all options 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; //use SIRT unless method=1 00204 00205 00206 // call SIRT to reconstruct 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 // write the reconstructed volume to disk 00223 EMUtil::ImageType WRITE_SPI = EMUtil::IMAGE_SINGLE_SPIDER; 00224 if ( my2dpid == 0 ) { 00225 xvol->write_image(voutfname, 0, WRITE_SPI); 00226 } 00227 00228 // cleanup 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; // main 00244 }