#include "mpi.h"
#include "emdata.h"
#include "sirt_Cart.h"
#include "utilcomm_Cart.h"
#include "project3d_Cart.h"
#include "project3d.h"
Include dependency graph for sirt_Cart.cpp:
Go to the source code of this file.
Defines | |
#define | ROW 0 |
#define | COL 1 |
Functions | |
int | recons3d_sirt_mpi_Cart (MPI_Comm comm_2d, MPI_Comm comm_row, MPI_Comm comm_col, EMData **images, float *angleshift, EMData *&xvol, int nangloc, int radius, float lam, int maxit, std::string symmetry, float tol) |
#define COL 1 |
Definition at line 12 of file sirt_Cart.cpp.
#define ROW 0 |
Definition at line 11 of file sirt_Cart.cpp.
int recons3d_sirt_mpi_Cart | ( | MPI_Comm | comm_2d, | |
MPI_Comm | comm_row, | |||
MPI_Comm | comm_col, | |||
EMData ** | images, | |||
float * | angleshift, | |||
EMData *& | xvol, | |||
int | nangloc, | |||
int | radius, | |||
float | lam, | |||
int | maxit, | |||
std::string | symmetry, | |||
float | tol | |||
) |
Definition at line 16 of file sirt_Cart.cpp.
References bckpj3_Cart(), COL, cord, dm, EMDeleteArray(), fwdpj3_Cart(), EMAN::EMData::get_data(), getcb2sph(), getnnz(), ierr, images, make_proj_mat(), myptrstart, nnz, nnzloc, nrays, nraysloc, nx, phi, PI, ptrs, ROW, EMAN::EMData::set_size(), setpart_gr1(), sph2cb(), sphpart(), sqrt(), theta, and EMAN::EMData::to_zero().
Referenced by main().
00022 { 00023 MPI_Status mpistatus; 00024 int ncpus, my2dpid, ierr; 00025 int mycoords[2], dims[2], periods[2]; 00026 int srpid, srcoords[2]; // Send/receive info 00027 double t0; 00028 00029 // Get dims and my coordinates 00030 MPI_Cart_get(comm_2d, 2, dims, periods, mycoords); 00031 MPI_Comm_rank(comm_2d, &my2dpid); //Get my pid in the new 2D topology 00032 ncpus = dims[ROW]*dims[COL]; 00033 00034 int * psize; 00035 int * nbase; 00036 00037 int nangles; 00038 psize = new int[dims[ROW]]; 00039 nbase = new int[dims[ROW]]; 00040 MPI_Allreduce(&nangloc, &nangles, 1, MPI_INT, MPI_SUM, comm_col); 00041 00042 int nsym = 0; 00043 // get image size from first image 00044 int nx = images[0]->get_xsize(); 00045 00046 // make radius as large as possible if the user didn't provide one 00047 if ( radius == -1 ) radius = nx/2 - 1; 00048 00049 Vec3i volsize, origin; 00050 volsize[0] = nx; 00051 volsize[1] = nx; 00052 volsize[2] = nx; 00053 origin[0] = nx/2+1; 00054 origin[1] = nx/2+1; 00055 origin[2] = nx/2+1; 00056 00057 // this is not currently needed, because the stack that gets passed to 00058 // sirt will have its background subtracted already 00059 // ierr = CleanStack(comm, images, nangloc, radius, volsize, origin); 00060 00061 // vector of symmetrized angles 00062 std::vector<float> symangles(3,0.0); 00063 00064 // kluge, make sure if its 1.0 + epsilon it still works; 00065 float old_rnorm = 1.00001; 00066 00067 00068 // Now distribute the volume (in spherical format) among columns of 00069 // processors and use nnz to determine the splitting. 00070 // Note: ptrs and coord are on all processors, but voldata is only on 00071 // Proc 0 00072 00073 int nrays, nnz; 00074 00075 ierr = getnnz(volsize, radius, origin, &nrays, &nnz); 00076 00077 int nnzloc, nraysloc; 00078 int * ptrs = new int[nrays+1]; 00079 int * cord = new int[3*nrays]; 00080 int *nnzpart = new int[dims[COL]]; 00081 int *nnzbase = new int[dims[COL]+1]; 00082 int *ptrstart = new int[dims[COL]+1]; 00083 00084 ierr = getcb2sph(volsize, radius, origin, nnz, ptrs, cord); 00085 nnzloc = setpart_gr1(comm_2d, nnz, nnzpart, nnzbase); 00086 nraysloc = sphpart(comm_2d, nrays, ptrs, nnzbase, ptrstart); 00087 00088 int myptrstart = ptrstart[mycoords[COL]]; 00089 int nnzall[dims[COL]]; 00090 for (int i = 0; i<dims[COL]; i++) 00091 nnzall[i] = ptrs[ptrstart[i+1]] - ptrs[ptrstart[i]]; 00092 00093 nnzloc = nnzall[mycoords[COL]]; 00094 00095 float *bvol_loc = new float[nnzloc]; 00096 float *bvol = new float[nnzloc]; 00097 float *xvol_sphloc = new float[nnzloc]; 00098 float *pxvol_loc = new float[nnzloc]; 00099 float *pxvol = new float[nnzloc]; 00100 float * grad_loc = new float[nnzloc]; 00101 for (int i=0; i< nnzloc; i++){ 00102 xvol_sphloc[i] = 0.0; 00103 bvol[i] = 0.0; 00104 bvol_loc[i] = 0.0; 00105 pxvol_loc[i] = 0.0; 00106 pxvol[i] = 0.0; 00107 grad_loc[i] = 0.0; 00108 } 00109 00110 EMData * current_image; 00111 float phi, theta, psi; 00112 Transform3D RA; 00113 Transform3D Tf; 00114 nsym = Tf.get_nsym(symmetry); 00115 Transform3D::EulerType EULER_SPIDER = Transform3D::SPIDER; 00116 Dict angdict; 00117 //printf("nsym = %d\n", nsym); 00118 int iter = 1; 00119 00120 double rnorm = 0.0, rnorm_loc = 0.0; 00121 double bnorm = 0.0, bnorm_loc = 0.0; 00122 //float * grad = new float[nnz]; 00123 float * image_data; 00124 float * projected_data_loc = new float[nangloc*nx*nx]; 00125 float * projected_data = new float[nangloc*nx*nx]; 00126 00127 float dm[8]; 00128 00129 int restarts = 0; 00130 00131 t0 = MPI_Wtime(); 00132 00133 while (iter <= maxit) { 00134 if ( iter == 1 ) { 00135 if ( restarts == 0 ) { 00136 // only do this if we aren't restarting due to lam being 00137 // too large 00138 for ( int i = 0 ; i < nangloc ; ++i ) { 00139 current_image = images[i]; 00140 image_data = current_image->get_data(); 00141 // retrieve the angles and shifts associated with 00142 // each image from the array angleshift. 00143 phi = angleshift[5*i + 0]; 00144 theta = angleshift[5*i + 1]; 00145 psi = angleshift[5*i + 2]; 00146 00147 // need to change signs here because the input shifts 00148 // are shifts associated with 2-D images. Because 00149 // the projection operator actually shifts the volume 00150 // the signs should be negated here 00151 dm[6] = -angleshift[5*i + 3]; 00152 dm[7] = -angleshift[5*i + 4]; 00153 // make an instance of Transform3D with the angles 00154 RA = Transform3D(EULER_SPIDER, phi, theta, psi); 00155 for ( int ns = 1 ; ns < nsym + 1 ; ++ns ) { 00156 // compose it with each symmetry rotation in turn 00157 // shifts (stored in dm[6] and dm[7] remain fixed 00158 Tf = Tf.get_sym(symmetry, ns) * RA; 00159 angdict = Tf.get_rotation(EULER_SPIDER); 00160 phi = (float) angdict["phi"] * PI/180.0; 00161 theta = (float) angdict["theta"] * PI/180.0; 00162 psi = (float) angdict["psi"] * PI/180.0; 00163 make_proj_mat(phi, theta, psi, dm); 00164 // accumulate the back-projected images in bvol_loc 00165 ierr = bckpj3_Cart(volsize, nraysloc, nnzloc, dm, 00166 origin, radius, ptrs, cord, 00167 myptrstart, image_data, bvol_loc); 00168 } 00169 } 00170 // reduce bvol_loc so each processor has the full volume 00171 //ierr = MPI_Allreduce(bvol_loc, bvol, nnz, MPI_FLOAT, MPI_SUM, 00172 // comm); 00173 // Now an all reduce along the columns 00174 ierr = MPI_Allreduce (bvol_loc, bvol, nnzloc, MPI_FLOAT, 00175 MPI_SUM, comm_col); 00176 00177 } 00178 00179 // calculate the norm of the backprojected volume 00180 bnorm_loc = 0.0; 00181 for ( int j = 0 ; j < nnzloc ; ++j ) { 00182 bnorm_loc += bvol[j] * (double) bvol[j]; 00183 grad_loc[j]= bvol[j]; 00184 } 00185 ierr = MPI_Allreduce (&bnorm_loc, &bnorm, 1, MPI_DOUBLE, MPI_SUM, 00186 comm_row); 00187 00188 bnorm /= nnz; 00189 bnorm = sqrt(bnorm); 00190 // printf("bnorm = %f\n", bnorm); 00191 } else { 00192 00193 // iterate over symmetry 00194 for ( int ns = 1 ; ns < nsym + 1 ; ++ns ) { 00195 // reset local and global 2-D projections to zeros 00196 for (int i=0; i<nangloc*nx*nx; i++){ 00197 projected_data_loc[i] = 0.0; 00198 projected_data[i] = 0.0; 00199 } 00200 00201 // project from 3-D to 2-D 00202 for ( int i = 0 ; i < nangloc ; ++i ) { 00203 // retrieve the angles and shifts from angleshift 00204 RA = Transform3D(EULER_SPIDER, angleshift[5*i + 0], 00205 angleshift[5*i + 1], angleshift[5*i + 2]); 00206 00207 // need to change signs here because the input shifts 00208 // are shifts associated with 2-D images. Because 00209 // the projection operator actually shifts the volume 00210 // the signs should be negated here 00211 dm[6] = -angleshift[5*i + 3]; 00212 dm[7] = -angleshift[5*i + 4]; 00213 00214 // apply symmetry transformation 00215 Tf = Tf.get_sym(symmetry, ns) * RA; 00216 angdict = Tf.get_rotation(EULER_SPIDER); 00217 00218 phi = (float) angdict["phi"] * PI/180.0; 00219 theta = (float) angdict["theta"] * PI/180.0; 00220 psi = (float) angdict["psi"] * PI/180.0; 00221 make_proj_mat(phi, theta, psi, dm); 00222 00223 ierr = fwdpj3_Cart(volsize, nraysloc, nnzloc, dm, 00224 origin, radius, ptrs, cord, myptrstart, 00225 xvol_sphloc, &projected_data_loc[nx*nx*i]); 00226 } 00227 00228 00229 // Now perform global sum on 2-D images along the rows 00230 ierr = MPI_Allreduce(projected_data_loc, projected_data, 00231 nangloc*nx*nx, MPI_FLOAT, MPI_SUM, comm_row); 00232 00233 // backproject from 2-D to 3-D 00234 for ( int i = 0 ; i < nangloc ; ++i ) { 00235 // retrieve the angles and shifts from angleshift 00236 RA = Transform3D(EULER_SPIDER, angleshift[5*i + 0], 00237 angleshift[5*i + 1], angleshift[5*i + 2]); 00238 00239 // need to change signs here because the input shifts 00240 // are shifts associated with 2-D images. Because 00241 // the projection operator actually shifts the volume 00242 // the signs should be negated here 00243 dm[6] = -angleshift[5*i + 3]; 00244 dm[7] = -angleshift[5*i + 4]; 00245 00246 // apply symmetry transformation 00247 Tf = Tf.get_sym(symmetry, ns) * RA; 00248 angdict = Tf.get_rotation(EULER_SPIDER); 00249 00250 // reset the array in which projected data are stored 00251 phi = (float) angdict["phi"] * PI/180.0; 00252 theta = (float) angdict["theta"] * PI/180.0; 00253 psi = (float) angdict["psi"] * PI/180.0; 00254 make_proj_mat(phi, theta, psi, dm); 00255 00256 // accumulate P^TPxvol in pxvol_loc 00257 ierr = bckpj3_Cart(volsize, nraysloc, nnzloc, dm, origin, 00258 radius, ptrs, cord, myptrstart, 00259 &projected_data[nx*nx*i], pxvol_loc); 00260 } 00261 } // endfor (ns=1,nsym) 00262 00263 ierr = MPI_Allreduce(pxvol_loc, pxvol, nnzloc, MPI_FLOAT, MPI_SUM, 00264 comm_col); 00265 00266 for ( int j = 0 ; j < nnzloc ; ++j ) { 00267 grad_loc[j] = bvol[j]; 00268 grad_loc[j] -= pxvol[j]; 00269 } 00270 } 00271 00272 rnorm_loc = 0.0; 00273 for ( int j = 0 ; j < nnzloc ; ++j ) { 00274 rnorm_loc += grad_loc[j]* (double) grad_loc[j]; 00275 } 00276 // printf("I am (%d, %d), with rnorm = %11.3e\n", 00277 // mycoords[ROW], mycoords[COL], rnorm); 00278 ierr = MPI_Allreduce (&rnorm_loc, &rnorm, 1, MPI_DOUBLE, MPI_SUM, 00279 comm_row); 00280 rnorm /= nnz; 00281 rnorm = sqrt(rnorm); 00282 if ( my2dpid == 0 ) 00283 printf("iter = %3d, rnorm / bnorm = %11.3e, rnorm = %11.3e\n", 00284 iter, rnorm / bnorm, rnorm); 00285 // if on the second pass, rnorm is greater than bnorm, 00286 // lam is probably set too high reduce it by a factor of 2 00287 // and start over 00288 // if ( iter == 2 && rnorm / bnorm > old_rnorm ) { 00289 if ( rnorm / bnorm > old_rnorm ) { 00290 // but don't do it more than 20 times 00291 if ( restarts > 20 ) { 00292 if ( my2dpid == 0 ) 00293 printf("Failure to converge, even with lam = %f\n", lam); 00294 break; 00295 } else { 00296 ++restarts; 00297 iter = 1; 00298 lam /= 2.0; 00299 // reset these 00300 // kluge, make sure if its 1.0 + epsilon it still works 00301 old_rnorm = 1.0001; 00302 for ( int j = 0 ; j < nnzloc ; ++j ) { 00303 xvol_sphloc[j] = 0.0; 00304 pxvol_loc[j] = 0.0; 00305 } 00306 if ( my2dpid == 0 ) 00307 printf("reducing lam to %11.3e, restarting\n", lam); 00308 continue; 00309 } // endif (restarts) 00310 } // endif (rnorm/bnorm) 00311 00312 // if changes are sufficiently small, 00313 // or if no further progress is made, terminate 00314 if ( rnorm / bnorm < tol || rnorm / bnorm > old_rnorm ) { 00315 if ( my2dpid == 0 ) 00316 printf("Terminating with rnorm/bnorm = %11.3e, "); 00317 printf("tol = %11.3e, old_rnorm = %11.3e\n", 00318 rnorm/bnorm, tol, old_rnorm); 00319 break; 00320 } 00321 // update the termination threshold 00322 old_rnorm = rnorm / bnorm; 00323 // update the reconstructed volume 00324 00325 for ( int j = 0 ; j < nnzloc ; ++j ) { 00326 xvol_sphloc[j] += lam * grad_loc[j]; 00327 // reset it so it's ready to accumulate for the next iteration 00328 pxvol_loc[j] = 0.0; 00329 } 00330 ++iter; 00331 } 00332 00333 if (my2dpid == 0) printf("Total time in SIRT = %11.3e\n", MPI_Wtime()-t0); 00334 00335 // Bring all parts of the spherical volume back together and turn it 00336 // into cube format. 00337 if (mycoords[ROW] == 0 && mycoords[COL] != 0 ){ 00338 srcoords[ROW] = 0; 00339 srcoords[COL] = 0; 00340 MPI_Cart_rank(comm_2d, srcoords, &srpid); 00341 MPI_Send(xvol_sphloc, nnzloc, MPI_FLOAT, 0, my2dpid, comm_2d); 00342 } 00343 00344 if (mycoords[ROW] == 0 && mycoords[COL] == 0 ){ 00345 xvol->set_size(nx, nx, nx); 00346 xvol->to_zero(); 00347 float * voldata = xvol->get_data(); 00348 float * xvol_sph = new float[nnz]; 00349 00350 for(int i=0; i<nnzloc; i++) 00351 xvol_sph[i] = xvol_sphloc[i]; 00352 00353 for(int i=1; i< dims[COL]; i++){ 00354 srcoords[ROW] = 0; 00355 srcoords[COL] = i; 00356 MPI_Cart_rank(comm_2d, srcoords, &srpid); 00357 MPI_Recv(&xvol_sph[ptrs[ptrstart[i]]-1], nnzall[i], MPI_FLOAT, 00358 srpid, srpid, comm_2d, &mpistatus); 00359 } 00360 00361 // unpack the spherical volume back out into the original EMData object 00362 ierr = sph2cb(xvol_sph, volsize, nrays, radius, nnz, ptrs, cord, voldata); 00363 EMDeleteArray(xvol_sph); 00364 } // endif (mycord[ROW]== 0...) 00365 EMDeleteArray(grad_loc); 00366 EMDeleteArray(pxvol_loc); 00367 EMDeleteArray(bvol_loc); 00368 00369 EMDeleteArray(ptrs); 00370 EMDeleteArray(cord); 00371 //EMDeleteArray(projected_data); 00372 00373 EMDeleteArray(psize); 00374 EMDeleteArray(nbase); 00375 00376 delete [] projected_data_loc; 00377 00378 return 0; // recons3d_sirt_mpi 00379 }