EMAN2
ali3d_d_mpi.cpp
Go to the documentation of this file.
00001 #include "stdlib.h"
00002 #include "mpi.h"
00003 
00004 // include EMAN2
00005 #include "emdata.h"
00006 #include "emassert.h"
00007 #include "projector.h"
00008 #include "transform.h"
00009 #include "reconstructor.h"
00010 
00011 #include "ali3d_d_mpi.h"
00012 #include "utilcomm.h"
00013 #include "alignoptions.h"
00014 #include "sirt.h"
00015 #include "project3d.h"
00016 #include "fgcalc.h"
00017 
00018 using namespace EMAN;
00019 
00020 int ali3d_d( MPI_Comm comm, EMData*& volume, EMData** projdata, 
00021              EMData** cleandata, float *angleshift, int nloc, 
00022              AlignOptions& options, char* fname_base)
00023 {
00024     int mypid, ncpus;
00025     int ierr;
00026     MPI_Status mpistatus;
00027     double timer, total_timer;
00028         
00029     MPI_Comm_rank(comm,&mypid);
00030     MPI_Comm_size(comm,&ncpus);
00031 
00032     int * psize, * nbase, nangloc, nang;
00033     psize = new int[ncpus];
00034     nbase = new int[ncpus];
00035     MPI_Allreduce(&nloc, &nang, 1, MPI_INT, MPI_SUM, comm);
00036     nangloc = setpart(comm, nang, psize, nbase);
00037 
00038     char out_fname[64];
00039     std::ofstream fsc_out;
00040 
00041     int ndim, nx, ny, nz;
00042 
00043     ndim = volume->get_ndim();
00044 
00045     nx = volume->get_xsize();
00046     ny = volume->get_ysize();
00047     nz = volume->get_zsize();
00048 
00049     int first_ring = options.get_first_ring();
00050     int last_ring  = options.get_last_ring();
00051     int rstep      = options.get_rstep();
00052     int ri         = options.get_ri();
00053     float xrng     = options.get_xrng();
00054     float yrng     = options.get_yrng();
00055     float step     = options.get_step();
00056     float dtheta   = options.get_dtheta();
00057     bool CTF       = options.get_CTF();
00058     bool have_angles = options.get_have_angles();
00059     bool use_sirt  = options.get_use_sirt();
00060     std::string symmetry = options.get_symmetry();
00061     int max_iter       = options.get_maxit();
00062 
00063     Dict maskdict;
00064     maskdict["radius"] = last_ring;
00065     maskdict["fill"]   = 1;
00066     
00067     EMData * mask3D = options.get_mask3D();
00068     if (mask3D == NULL) {
00069         mask3D = new EMData();
00070         
00071         mask3D->set_size(nx,ny,nz);
00072         mask3D->process_inplace("testimage.circlesphere", maskdict);
00073     }
00074                         
00075     EMData * mask2D = new EMData();
00076     mask2D->set_size(nx,ny,1);
00077     mask2D->process_inplace("testimage.circlesphere", maskdict);
00078     if (first_ring > 0) {
00079         // subtract off the smaller ring
00080         EMData * inner_mask2D = new EMData();
00081         inner_mask2D->set_size(nx,ny,1);
00082         maskdict["radius"] = first_ring;
00083         inner_mask2D->process_inplace("testimage.circlesphere", maskdict);
00084         *mask2D -= *inner_mask2D;
00085         EMDeletePtr(inner_mask2D);
00086     }   
00087     
00088     EMData ** sirt_images;
00089     float sirt_tol, sirt_lam;
00090     int sirt_maxit;
00091     Dict recons_params;
00092     std::string recons_name;
00093     float * even_odd_angleshift = new float[5*nloc];
00094 
00095 
00096     if ( use_sirt ) {
00097         sirt_images = new EMData*[nloc];
00098         sirt_tol = options.get_sirt_tol();
00099         sirt_lam = options.get_sirt_lam();
00100         sirt_maxit = options.get_sirt_maxit();
00101     } else {
00102         recons_name = "nn4";
00103         recons_params["symmetry"] = symmetry;
00104         recons_params["size"]     = nx;
00105         recons_params["npad"]     = 4; // will this ever be != 4?
00106 
00107         Dict CTF_params;
00108         CTF_params["sign"] = -1.0;
00109         int CTF_applied;
00110         int padffted;
00111         if (CTF) {
00112             // filter each image according to its ctf header info
00113             for ( int i = 0 ; i < nloc ; ++i ) {
00114                 CTF_applied = projdata[i]->get_attr("ctf_applied");
00115                 if ( CTF_applied == 0 ) {
00116                     CTF_params["defocus"]      = projdata[i]->get_attr("defocus");
00117                     CTF_params["Cs"]           = projdata[i]->get_attr("Cs");
00118                     CTF_params["Pixel_size"]   = projdata[i]->get_attr("Pixel_size");
00119                     CTF_params["B_factor"]     = projdata[i]->get_attr("B_factor");
00120                     CTF_params["amp_contrast"] = projdata[i]->get_attr("amp_contrast");
00121                     projdata[i]->process_inplace("filter.CTF_",CTF_params);
00122                 }
00123                 // I think the CTF filter should be applied before the images have the background subtracted,
00124                 // but I'm not really sure.
00125                 CTF_applied = cleandata[i]->get_attr("ctf_applied");
00126                 if ( CTF_applied == 0 ) {
00127                     CTF_params["defocus"]      = cleandata[i]->get_attr("defocus");
00128                     CTF_params["Cs"]           = cleandata[i]->get_attr("Cs");
00129                     CTF_params["Pixel_size"]   = cleandata[i]->get_attr("Pixel_size");
00130                     CTF_params["B_factor"]     = cleandata[i]->get_attr("B_factor");
00131                     CTF_params["amp_contrast"] = cleandata[i]->get_attr("amp_contrast");
00132                     cleandata[i]->process_inplace("filter.CTF_",CTF_params);
00133                 }
00134             }
00135             // set a few more things for the reconstructor to deal with
00136             recons_name = "nn4_ctf";
00137             recons_params["snr"] = options.get_snr();
00138             recons_params["sign"] = 1.0; // confirm this, sign is not defined in the current python code.
00139             try {
00140                 padffted = projdata[0]->get_attr("padffted");
00141             } catch ( std::exception& e ) { // the only exception thrown by get_attr() is NotExistingObjectException()
00142                 padffted = 0;
00143             }
00144             if (padffted == 1) recons_params["size"] = (float) recons_params["size"] / (float) recons_params["npad"];
00145         }
00146     }
00147 
00148     // The below are done in proj_ali_incore
00149     // Perhaps could be moved inside a major loop iterating over 
00150     // different dtheta
00151     // xrng, and xshift settings.
00152     std::string mode("F");
00153 
00154     // The below ensures that dtheta is sufficiently small to make the inverse problem over-determined
00155     int num_ref = 0;
00156     std::vector<float> ref_angles;
00157     int max_dimension = (nx > ny ? nx : ny);
00158     max_dimension = (max_dimension > nz ? max_dimension : nz);
00159     while (true) {
00160         ref_angles = Util::even_angles(dtheta);
00161         num_ref = ref_angles.size() / 3;
00162         if (num_ref >= max_dimension) break;
00163         dtheta -= 1.0;
00164         if (mypid == 0) printf("Problem is underdetermined: decreasing dtheta to %f\n",dtheta);
00165     }
00166     // should this ever not be true?  Is it really needed?
00167     bool phi_eq_minus_psi = true;
00168     if (phi_eq_minus_psi) {
00169         for ( int i = 0 ; i < num_ref ; ++ i ) {
00170             ref_angles[3*i + 2] = fmod(720.0 - ref_angles[3*i], 360.0);
00171         }
00172     }
00173     
00174     int numref; // used after call to multiref_polar_ali_2d, stores the index of the image with max correlation
00175 
00176     Util::mul_img(volume, mask3D);
00177 
00178     std::vector<int> numr = Numrinit(first_ring, last_ring, rstep, mode);
00179     std::vector<float> wr = ringwe(numr, mode);
00180 
00181     std::vector<EMData *> ref_proj_rings;
00182     EMData * ref_proj_ptr;
00183     EMData * cimage_ptr;
00184     Dict volparams;
00185     volparams["angletype"] = "SPIDER";
00186     volparams["radius"] = ri;
00187     Dict proj_params;
00188     proj_params["ctf_applied"] = 0;
00189     std::vector<float> anglelist(3,0.0);
00190 
00191     float sxo, syo;
00192     std::vector<float> best_alignment(6,0.0);
00193 
00194     Transform3D::EulerType EULER_SPIDER = Transform3D::SPIDER;
00195     EMUtil::ImageType WRITE_SPI = EMUtil::IMAGE_SINGLE_SPIDER;
00196     float compphi, comptheta, comppsi;
00197     float angb, sxb, syb, ct;
00198 
00199     float cnx = floor(nx / 2.0) + 1.0;
00200     float cny = floor(ny / 2.0) + 1.0;
00201 
00202     char iter_string[4];
00203 
00204     Dict make_ref_proj_dict;
00205     make_ref_proj_dict["mask"] = mask2D;
00206     make_ref_proj_dict["no_sigma"] = 1;
00207 
00208     float multiref_res, unified_res;
00209     float * volsph, * voldata;
00210     Vec3i volsize;
00211     Vec3i origin;
00212     volsize[0] = nx;
00213     volsize[1] = ny;
00214     volsize[2] = nz;
00215     origin[0] = nx/2 + 1;
00216     origin[1] = ny/2 + 1;
00217     origin[2] = nz/2 + 1;
00218 
00219     int nnz, nrays;
00220 
00221     ierr = getnnz(volsize, ri, origin, &nrays, &nnz);
00222     int * ptrs = new int[nrays+1];
00223     int * cord = new int[3*nrays];
00224     volsph = new float[nnz];
00225     voldata = volume->get_data();
00226     ierr = cb2sph(voldata, volsize, ri, origin, nnz, ptrs, cord, volsph); 
00227     
00228     float angtrs[5];
00229     float * img_data;
00230 
00231     EMData * vol1 = NULL;
00232     EMData * vol2 = NULL;
00233 
00234     EMData * fftvol;
00235     EMData * weight;
00236     Reconstructor* r;
00237     float * fftvol_send;
00238     float * weight_send;
00239     int fftvol_size;
00240     int weight_size;
00241     float * fftvol_recv;
00242     float * weight_recv;
00243 
00244     std::vector<float> stats;
00245     std::vector<float> fsc_result;
00246     int fsc_size;
00247 
00248     float ring_width = 1.0;
00249     float phi, theta, psi;
00250     int filter_cutoff_index;
00251 
00252     // these are settings to compare with the fourier shell correlation 
00253     float filter_coeff_high = 0.9; 
00254     float filter_freq_low = 0.25;
00255     float hcut_default = 0.49;
00256 
00257     // these get passed to the filter
00258     float lcut, hcut; 
00259 
00260     Dict btwl_dict;
00261 
00262     std::vector<float> ph_cog;
00263     Dict cog_dict;
00264 
00265     total_timer = MPI_Wtime();
00266 
00267     float * peak_corr = new float[nloc];
00268     std::ofstream corr_out;
00269 
00270     // These keep track of how many times we take new parameters based on projection matching
00271     int new_angles, old_angles, new_total, old_total;
00272     
00273     for ( int i = 0 ; i < max_iter ; ++i ) { // This loop is as in ali3d_d
00274         new_angles = 0;
00275         old_angles = 0;
00276         // entering proj_ali_incore
00277         // generate reference projections and push pointers to them on the back of ref_proj_rings
00278         if (mypid == 0) std::cout << "Generating " << num_ref << " reference projections..." << std::endl;
00279         for ( int j = 0 ; j < num_ref ; ++j ) { 
00280             for ( int k = 0 ; k < 3 ; ++k) anglelist[k] = ref_angles[3 * j + k];
00281             volparams["anglelist"] = anglelist;
00282             ref_proj_ptr = volume->project("chao", volparams);// This uses a version of fwdpj3() internally
00283             proj_params["phi"] = anglelist[0];
00284             proj_params["theta"] = anglelist[1];
00285             proj_params["psi"] = anglelist[2];
00286             ref_proj_ptr->set_attr_dict(proj_params);
00287             ref_proj_ptr->process_inplace("normalize.mask", make_ref_proj_dict);
00288 
00289             cimage_ptr = Util::Polar2Dm(ref_proj_ptr, cnx, cny, numr, mode);
00290             Util::Frngs(cimage_ptr, numr);
00291             Applyws(cimage_ptr, numr, wr);
00292             ref_proj_rings.push_back(cimage_ptr);
00293 
00294         }
00295         // match the data images to the reference projections
00296         if (mypid == 0) std::cout << "Matching data to reference projections..." << std::endl;
00297         timer = MPI_Wtime();
00298         for ( int j = 0 ; j < nloc ; ++j ) {
00299             if (have_angles) {
00300                 sxo = angleshift[5*j + 3];
00301                 syo = angleshift[5*j + 4];
00302             } else {
00303                 sxo = projdata[j]->get_attr("s2x");
00304                 syo = projdata[j]->get_attr("s2y");
00305             }
00306             best_alignment = Util::multiref_polar_ali_2d(projdata[j], ref_proj_rings, xrng, yrng, step, mode, numr, cnx - sxo, cny - syo);
00307             numref = (int) best_alignment[4];
00308 
00309             peak_corr[j] = best_alignment[5]; // best_alignment[5] = peak
00310             // Here we implement what goes on in compose_transform2(), from utilities.py
00311             Transform3D R1(EULER_SPIDER , 0.0, 0.0, 0.0);
00312             R1.set_posttrans(Vec3f(best_alignment[1],best_alignment[2],0.0)); // best_alignment[1], best_alignment[2] = sxs, sys
00313             R1.set_scale(1.0);
00314 
00315             Transform3D R2(EULER_SPIDER, 0.0, 0.0, -best_alignment[0]); // best_alignment[0] = ang
00316             R2.set_posttrans(Vec3f(0.0,0.0,0.0));
00317             R2.set_scale(1.0);
00318 
00319             Transform3D Rcomp = R2*R1;
00320             Dict compeuler = Rcomp.get_rotation(EULER_SPIDER);
00321 
00322             compphi   = fmod((float) compeuler["phi"]   + 360.0, 360.0);
00323             comptheta = fmod((float) compeuler["theta"] + 360.0, 360.0);
00324             comppsi   = fmod((float) compeuler["psi"]   + 360.0, 360.0);
00325 
00326             angb = fmod((float) compphi + comppsi,  (float) 360.0);  
00327             sxb  = Rcomp.at(0,3);
00328             syb  = Rcomp.at(1,3);
00329             ct   = Rcomp.get_scale();
00330             // end call to compose_transform2()
00331             Dict composition_dict;
00332             // perhaps should do all of these with set_attr(), rather than set_attr_dict()
00333             if (best_alignment[3] != 0.0) { // best_alignment[3] = mirror
00334                 composition_dict["phi"]   = fmod(ref_angles[3 * numref] + 540.0, 360.0);
00335                 composition_dict["theta"] = 180.0 - ref_angles[3 * numref + 1];
00336                 composition_dict["psi"]   = fmod(540.0 - ref_angles[3 * numref + 2] + angb, 360.0);
00337                 composition_dict["s2x"]   = sxb + sxo;
00338                 composition_dict["s2y"]   = syb + syo;
00339             } else {
00340                 composition_dict["phi"]   = ref_angles[3 * numref];
00341                 composition_dict["theta"] = ref_angles[3 * numref + 1];
00342                 composition_dict["psi"]   = fmod( ref_angles[3 * numref + 2] + angb + 360.0, 360.0);
00343                 composition_dict["s2x"]   = sxb + sxo;
00344                 composition_dict["s2y"]   = syb + syo;
00345             }
00346             if (have_angles) {
00347                 // compare residual for this image under old angles (from refinement) and new angles (from projection matching)
00348                 img_data = projdata[j]->get_data();
00349 
00350                 angtrs[0] = (float) composition_dict["phi"]   * (PI/180.0);
00351                 angtrs[1] = (float) composition_dict["theta"] * (PI/180.0);
00352                 angtrs[2] = (float) composition_dict["psi"]   * (PI/180.0);
00353                 angtrs[3] = (float) composition_dict["s2x"]   * -1.0;
00354                 angtrs[4] = (float) composition_dict["s2y"]   * -1.0;
00355                 
00356                 ierr = fcalc(volsph, volsize, 
00357                              nnz, nrays, origin, ri, 
00358                              ptrs, cord, angtrs, 1, 
00359                              img_data, &multiref_res);
00360 
00361                 // set angtrs to the values in angleshift
00362                 angtrs[0] = angleshift[5*j + 0] * (PI/180.0);
00363                 angtrs[1] = angleshift[5*j + 1] * (PI/180.0);
00364                 angtrs[2] = angleshift[5*j + 2] * (PI/180.0);
00365                 angtrs[3] = angleshift[5*j + 3] * -1.0;
00366                 angtrs[4] = angleshift[5*j + 4] * -1.0;
00367 
00368                 ierr = fcalc(volsph, volsize, 
00369                              nnz, nrays, origin, ri, 
00370                              ptrs, cord, angtrs, 1, 
00371                              img_data, &unified_res);
00372 
00373                 // whichever gives better residual gets used
00374                 if ( multiref_res < unified_res ) {
00375                     ++new_angles;
00376                     angleshift[5*j + 0] = composition_dict["phi"];
00377                     angleshift[5*j + 1] = composition_dict["theta"];
00378                     angleshift[5*j + 2] = composition_dict["psi"];
00379                     angleshift[5*j + 3] = composition_dict["s2x"];
00380                     angleshift[5*j + 4] = composition_dict["s2y"];
00381                 } else { // multiref_res >= unified_res
00382                     ++old_angles;
00383                     composition_dict["phi"]   = angleshift[5*j + 0];
00384                     composition_dict["theta"] = angleshift[5*j + 1];
00385                     composition_dict["psi"]   = angleshift[5*j + 2];
00386                     composition_dict["s2x"]   = angleshift[5*j + 3];
00387                     composition_dict["s2y"]   = angleshift[5*j + 4];
00388                 }
00389             } else { // don't have intial angles, and need to set them
00390                 angleshift[5*j + 0] = composition_dict["phi"];
00391                 angleshift[5*j + 1] = composition_dict["theta"];
00392                 angleshift[5*j + 2] = composition_dict["psi"];
00393                 angleshift[5*j + 3] = composition_dict["s2x"];
00394                 angleshift[5*j + 4] = composition_dict["s2y"];
00395             }
00396             // set the header of projdata[j] and cleandata[j] to the best angle and shift values
00397             projdata[j]->set_attr_dict(composition_dict);
00398             cleandata[j]->set_attr_dict(composition_dict);
00399         }
00400 
00401         MPI_Reduce(&new_angles, &new_total, 1, MPI_INT, MPI_SUM, 0, comm);
00402         MPI_Reduce(&old_angles, &old_total, 1, MPI_INT, MPI_SUM, 0, comm);
00403         if (mypid == 0) {
00404           printf("   Wall clock seconds for alignment, iteration %d = %11.3e\n",
00405                  i, MPI_Wtime() - timer);
00406           if (have_angles)
00407               printf("   \nUsing %d sets of parameters from projection matching, %d sets from previous unified refinement\n\n", new_total, old_total);
00408         }
00409 
00410         sprintf(out_fname, "corr_%s_pm%d.dat", fname_base, i);
00411         corr_out.open(out_fname);
00412         if (mypid == 0) {
00413             for ( int corr_index = 0 ; corr_index < nloc ; ++corr_index ) {
00414                 corr_out << std::scientific << peak_corr[corr_index] << std::endl;
00415             }
00416             for ( int k = 1 ; k < ncpus ; ++k ) {
00417                 ierr = MPI_Recv(peak_corr, psize[k], MPI_FLOAT, k, k, comm, &mpistatus);
00418                 for ( int corr_index = 0 ; corr_index < psize[k] ; ++corr_index ) {
00419                     corr_out << std::scientific << peak_corr[corr_index] << std::endl;
00420                 }
00421             }
00422         } else { // mypid != 0 , send my data to master to write out
00423             ierr = MPI_Send(peak_corr, nloc, MPI_FLOAT, 0, mypid, comm);
00424         }
00425         corr_out.close();
00426 
00427         timer = MPI_Wtime();
00428         // leaving proj_ali_incore
00429 
00430         // build two reconstructions
00431         if (mypid == 0) 
00432            printf("   Building even/odd volumes for FSC calculation...\n");
00433         if ( use_sirt ) {
00434 
00435             // construct a volume using even numbered images
00436             vol1 = new EMData();
00437             // lam should be based on number of images, right?
00438             for ( int j = 0 ; j < nloc ; j += 2 ) {
00439                 sirt_images[j/2] = cleandata[j];
00440                 even_odd_angleshift[5*(j/2) + 0] = angleshift[5*j + 0];
00441                 even_odd_angleshift[5*(j/2) + 1] = angleshift[5*j + 1];
00442                 even_odd_angleshift[5*(j/2) + 2] = angleshift[5*j + 2];
00443                 even_odd_angleshift[5*(j/2) + 3] = angleshift[5*j + 3];
00444                 even_odd_angleshift[5*(j/2) + 4] = angleshift[5*j + 4];
00445             }
00446             recons3d_sirt_mpi(comm, sirt_images, even_odd_angleshift, 
00447                               vol1, (nloc+1)/2, ri, sirt_lam, sirt_maxit, 
00448                               symmetry, sirt_tol);
00449 
00450             // construct a volume using odd numbered images
00451             vol2 = new EMData();
00452             for ( int j = 1 ; j < nloc ; j += 2 ) {
00453                 sirt_images[j/2] = cleandata[j];
00454                 even_odd_angleshift[5*(j/2) + 0] = angleshift[5*j + 0];
00455                 even_odd_angleshift[5*(j/2) + 1] = angleshift[5*j + 1];
00456                 even_odd_angleshift[5*(j/2) + 2] = angleshift[5*j + 2];
00457                 even_odd_angleshift[5*(j/2) + 3] = angleshift[5*j + 3];
00458                 even_odd_angleshift[5*(j/2) + 4] = angleshift[5*j + 4];
00459             }
00460             recons3d_sirt_mpi(comm, sirt_images, even_odd_angleshift, 
00461             vol2, nloc/2, ri, sirt_lam, sirt_maxit, symmetry, sirt_tol);
00462         } else {
00463 
00464             // construct a volume using even numbered images
00465             fftvol = new EMData();
00466             weight = new EMData();
00467             recons_params["fftvol"] = fftvol;
00468             recons_params["weight"] = weight;
00469             r = Factory<Reconstructor>::get(recons_name, recons_params);
00470             r->setup();
00471             for ( int j = 0 ; j < nloc ; j += 2 ) { 
00472                 phi   = projdata[j]->get_attr("phi");
00473                 theta = projdata[j]->get_attr("theta");
00474                 psi   = projdata[j]->get_attr("psi");
00475                 r->insert_slice(cleandata[j],
00476                                 Transform3D(EULER_SPIDER,phi,theta,psi));
00477             }
00478             fftvol_size = fftvol->get_xsize() 
00479                         * fftvol->get_ysize() * fftvol->get_zsize();
00480             weight_size = weight->get_xsize() 
00481                         * weight->get_ysize() * weight->get_zsize();
00482 
00483             fftvol_send = fftvol->get_data();
00484             weight_send = weight->get_data();
00485             fftvol_recv = new float[fftvol_size];
00486             weight_recv = new float[weight_size];
00487             MPI_Allreduce(fftvol_send, fftvol_recv, fftvol_size, 
00488                           MPI_FLOAT, MPI_SUM, comm);
00489             MPI_Allreduce(weight_send, weight_recv, weight_size, 
00490                           MPI_FLOAT, MPI_SUM, comm);
00491             for ( int j = 0 ; j < fftvol_size ; j += 1 ) 
00492                 fftvol_send[j] = fftvol_recv[j];
00493             for ( int j = 0 ; j < weight_size ; j += 1 ) 
00494                 weight_send[j] = weight_recv[j];
00495             EMDeleteArray(fftvol_recv);
00496             EMDeleteArray(weight_recv);
00497             vol1 = r->finish();
00498             EMDeletePtr(fftvol);
00499             EMDeletePtr(weight);
00500 
00501             // construct a volume using odd numbered images
00502             fftvol = new EMData();
00503             weight = new EMData();
00504             recons_params["fftvol"] = fftvol;
00505             recons_params["weight"] = weight;
00506             r = Factory<Reconstructor>::get(recons_name, recons_params);
00507             r->setup();
00508             for ( int j = 1 ; j < nloc ; j += 2 ) { 
00509                 phi   = projdata[j]->get_attr("phi");
00510                 theta = projdata[j]->get_attr("theta");
00511                 psi   = projdata[j]->get_attr("psi");
00512                 r->insert_slice(cleandata[j], 
00513                                 Transform3D(EULER_SPIDER,phi,theta,psi));
00514             }
00515 
00516             fftvol_size = fftvol->get_xsize() 
00517                         * fftvol->get_ysize() * fftvol->get_zsize();
00518             weight_size = weight->get_xsize() 
00519                         * weight->get_ysize() * weight->get_zsize();
00520 
00521             fftvol_send = fftvol->get_data();
00522             weight_send = weight->get_data();
00523             fftvol_recv = new float[fftvol_size];
00524             weight_recv = new float[weight_size];
00525             MPI_Allreduce(fftvol_send, fftvol_recv, fftvol_size, 
00526                           MPI_FLOAT, MPI_SUM, comm);
00527             MPI_Allreduce(weight_send, weight_recv, weight_size, 
00528                           MPI_FLOAT, MPI_SUM, comm);
00529 
00530             for ( int j = 0 ; j < fftvol_size ; j += 1 ) 
00531                fftvol_send[j] = fftvol_recv[j];
00532             for ( int j = 0 ; j < weight_size ; j += 1 ) 
00533                weight_send[j] = weight_recv[j];
00534 
00535             EMDeleteArray(fftvol_recv);
00536             EMDeleteArray(weight_recv);
00537             vol2 = r->finish();
00538             EMDeletePtr(fftvol);
00539             EMDeletePtr(weight);
00540         }
00541 
00542         // calculate and subtract the mean from vol1 and vol2, 
00543         // and apply the 3D mask
00544         //
00545         stats = Util::infomask(vol1, mask3D, false);
00546         *vol1 -= stats[0]; // stats[0] = mean
00547         Util::mul_img(vol1, mask3D);
00548         stats = Util::infomask(vol2, mask3D, false);
00549         *vol2 -= stats[0]; // stats[0] = mean
00550         Util::mul_img(vol2, mask3D);
00551 
00552         // calculate fsc
00553         fsc_result = vol1->calc_fourier_shell_correlation(vol2, ring_width);
00554         fsc_size = fsc_result.size() / 3;
00555         sprintf(out_fname, "fsc_%s_pm%d.dat", fname_base, i);
00556         if (mypid == 0) {
00557             fsc_out.open(out_fname);
00558             for ( int j = 0 ; j < fsc_size ; ++j ) {
00559                 // Note the indexing of fsc: the frequencies are in the 
00560                 // first fsc_size entries, then the correlation coeffs 
00561                 // after that
00562                 fsc_out << std::scientific 
00563                         << fsc_result[j] << '\t'
00564                         << fsc_result[j + fsc_size] << '\t'
00565                         << fsc_result[j + 2 * fsc_size] 
00566                         << std::endl;
00567             }
00568             fsc_out.close();
00569         }
00570         EMDeletePtr(vol1);
00571         EMDeletePtr(vol2);
00572 
00573         // and reconstruct the volume from all images
00574         if (mypid == 0) 
00575            std::cout << "Building 3-D reconstruction ... " << std::endl;
00576         if ( use_sirt ) {
00577             volume = new EMData();
00578             recons3d_sirt_mpi(comm, cleandata, angleshift, volume, nloc, 
00579                               ri, sirt_lam, sirt_maxit, symmetry, sirt_tol);
00580         } else {
00581             fftvol = new EMData();
00582             weight = new EMData();
00583             recons_params["fftvol"] = fftvol;
00584             recons_params["weight"] = weight;
00585             r = Factory<Reconstructor>::get(recons_name, recons_params);
00586             r->setup();
00587             for ( int j = 0 ; j < nloc ; ++j ) { 
00588                 phi   = projdata[j]->get_attr("phi");
00589                 theta = projdata[j]->get_attr("theta");
00590                 psi   = projdata[j]->get_attr("psi");
00591                 r->insert_slice(cleandata[j], 
00592                                 Transform3D(EULER_SPIDER,phi,theta,psi));
00593             }
00594 
00595             fftvol_size = fftvol->get_xsize() 
00596                         * fftvol->get_ysize() * fftvol->get_zsize();
00597             weight_size = weight->get_xsize()
00598                         * weight->get_ysize() * weight->get_zsize();
00599 
00600             fftvol_send = fftvol->get_data();
00601             weight_send = weight->get_data();
00602 
00603             fftvol_recv = new float[fftvol_size];
00604             weight_recv = new float[weight_size];
00605 
00606             MPI_Allreduce(fftvol_send, fftvol_recv, fftvol_size, 
00607                           MPI_FLOAT, MPI_SUM, comm);
00608             MPI_Allreduce(weight_send, weight_recv, weight_size, 
00609                           MPI_FLOAT, MPI_SUM, comm);
00610 
00611             for ( int j = 0 ; j < fftvol_size ; ++j ) {
00612                 fftvol_send[j] = fftvol_recv[j];
00613             }
00614             for ( int j = 0 ; j < weight_size ; ++j ) {
00615                 weight_send[j] = weight_recv[j];
00616             }
00617 
00618             EMDeleteArray(fftvol_recv);
00619             EMDeleteArray(weight_recv);
00620             EMDeletePtr(volume); 
00621 
00622             // r->finish() returns EMData *, must free the old one
00623             volume = r->finish();
00624             EMDeletePtr(fftvol);
00625             EMDeletePtr(weight);
00626         }
00627 
00628         sprintf(out_fname, "vol_%s_pm%d.spi", fname_base, i);
00629         if (mypid == 0) volume->write_image(out_fname, 0, WRITE_SPI);
00630         // End of reconstruction
00631 
00632         // and filter it based on fourier shell correlation
00633         filter_cutoff_index = 0;
00634         while (filter_cutoff_index < fsc_size
00635                && fsc_result[filter_cutoff_index + fsc_size] > filter_coeff_high
00636                && fsc_result[filter_cutoff_index] < filter_freq_low) {
00637             ++filter_cutoff_index;
00638         }
00639 
00640         lcut  = fsc_result[filter_cutoff_index * 3];
00641         // hcut = min(f_l + 0.1, f_h_default)
00642         hcut = (lcut + 0.1 < hcut_default ? lcut + 0.1 : hcut_default); 
00643 
00644         btwl_dict["low_cutoff_frequency"]  = lcut;
00645         btwl_dict["high_cutoff_frequency"] = hcut;
00646         volume->set_attr("npad",1);
00647         volume->process_inplace("filter.lowpass.butterworth", btwl_dict);
00648 
00649         // calculate the center of gravity
00650         ph_cog = volume->phase_cog();
00651         cog_dict["x_shift"] = -ph_cog[0];
00652         cog_dict["y_shift"] = -ph_cog[1];
00653         cog_dict["z_shift"] = -ph_cog[2];
00654         // vol = fshift(vol, -cs[0], -cs[1], -cs[2]
00655         volume->process_inplace("filter.shift", cog_dict);
00656 
00657         // and write to disk
00658         sprintf(out_fname, "vol_f_%s_pm%d.spi", fname_base, i);
00659         if (mypid == 0) volume->write_image(out_fname, 0, WRITE_SPI);
00660 
00661         if (mypid == 0) {
00662             printf("   Wall clock seconds for reconstructions and fsc, iteration %d = %11.3e\n",
00663                    i, MPI_Wtime() - timer);
00664         }
00665 
00666         // clean up ref_proj_rings
00667         ref_proj_rings.clear();
00668     }
00669 
00670     if (mypid == 0) {
00671         std::cout << "Total wall clock seconds for " << max_iter << " iterations on " << ncpus << " processors : ";
00672         std::cout << MPI_Wtime() - total_timer << std::endl;
00673     }
00674 
00675     if (options.get_mask3D() == NULL) EMDeletePtr(mask3D); // If it was created by new in this function, otherwise it belongs to someone else.
00676     EMDeletePtr(mask2D);        
00677     EMDeleteArray(peak_corr);
00678     EMDeleteArray(ptrs);
00679     EMDeleteArray(cord);
00680     EMDeleteArray(volsph);
00681     EMDeleteArray(psize);
00682     EMDeleteArray(nbase);
00683     if ( use_sirt) EMDeleteArray(sirt_images);
00684     EMDeleteArray(even_odd_angleshift);
00685 
00686     return 0;
00687         
00688 }
00689 
00690 std::vector<int> Numrinit(int first_ring, int last_ring, int skip, std::string mode)
00691 {
00692         float dpi;
00693         int MAXFFT = 32768;
00694 
00695         if (mode == "f" or mode == "F") {
00696                 dpi = 2*PI;
00697         } else {
00698                 dpi = PI;
00699         }
00700 
00701         std::vector<int> numr;
00702         int lcirc = 1;
00703         int jp, ip, exponent, m;
00704         for ( int k = first_ring ; k < last_ring + 1 ; k += skip ) {
00705             numr.push_back(k);
00706             jp = int(dpi * k + 0.5);
00707             exponent = -1;
00708             m = 1;
00709             while (m <= jp){
00710                 m <<= 1;
00711                 ++exponent;
00712             }
00713             // exponent is the largest int such that 2**exponent <= jp                                                                                                           
00714             ip = 2 << exponent;
00715             if ( (k + skip <= last_ring && jp > ip + ip / 2) || (k + skip > last_ring && jp > ip + ip / 5)) {
00716                 ip = MAXFFT < 2*ip ? MAXFFT : 2*ip;
00717             }
00718             numr.push_back(lcirc);
00719             numr.push_back(ip);
00720             lcirc += ip;
00721         }       
00722         return numr;
00723 }
00724 
00725 std::vector<float> ringwe(std::vector<int> numr, std::string mode)
00726 {
00727         float dpi;
00728 
00729         if (mode == "f" or mode == "F") {
00730                 dpi = 2*PI;
00731         } else {
00732                 dpi = PI;
00733         }
00734         
00735         int nring = numr.size() / 3;
00736         std::vector<float> wr(nring, 0.0);
00737         float maxrin = (float) numr.back();
00738         for ( int i = 0 ; i < nring ; ++i ) {
00739                 // No idea what this line is supposed to mean.
00740                 // I hope the order of operations hold...
00741                 wr[i] = numr[i * 3]*dpi/((float) numr[2+i*3]) * maxrin / ((float) numr[2+i*3]);
00742         }
00743         return wr;
00744 }
00745 
00746 int Applyws(EMData * circ, std::vector<int> numr, std::vector<float> wr)
00747 {
00748         int nring = numr.size() / 3;
00749         int maxrin = numr.back();
00750         int numr3i, numr2i;
00751         int jc;
00752         float w;
00753         for ( int i = 0 ; i < nring ; ++i ) {
00754                 numr3i = numr[2 + i * 3];
00755                 numr2i = numr[1 + i * 3] - 1;
00756                 w = wr[i];
00757                 // why didn't they overload EMData::operator[] ?
00758                 (*circ)(numr2i) *= w;
00759                 if (numr3i == maxrin) {
00760                         (*circ)(numr2i + 1) *= w;
00761                 } else {
00762                         (*circ)(numr2i + 1) *= 0.5 * w;
00763                 }
00764                 for ( int j = 2 ; j < numr3i ; ++j ) {
00765                         jc = j + numr2i;
00766                         (*circ)(jc) *= w;
00767                 }
00768         }
00769         return 0;
00770 }
00771 
00772 EMData * recons3d_4nn(std::string stack_name, std::vector<int> list_proj, int npad)
00773 // This function should NOT be used, since it reads images from a stack, rather than
00774 // an array of EMData pointers.
00775 {
00776         Transform3D::EulerType Ttype = Transform3D::SPIDER; // might need :: instead of .
00777         
00778         EMData * proj = new EMData();
00779         proj->read_image(stack_name, list_proj[0]);
00780         float phi   = proj->get_attr("phi");
00781         float theta = proj->get_attr("theta");
00782         float psi   = proj->get_attr("psi");
00783         int   size  = proj->get_xsize();
00784         int active  = proj->get_attr("active");
00785         // 
00786         Dict params;
00787         params["symmetry"] = "c1";
00788         params["size"]     = size;
00789         params["npad"]     = npad;
00790 
00791         if (size != proj->get_ysize()) {
00792                 std::cerr << "Use a square image!" << std::endl;
00793                 // Throw an exception, the image must be square.
00794         }
00795         Reconstructor* r = Factory<Reconstructor>::get("nn4", params);
00796         r->setup();
00797         if (active) {
00798                 r->insert_slice(proj, Transform3D(Ttype, phi, theta, psi));
00799         }
00800         int nimages = list_proj.size();
00801         for ( int i = 1 ; i < nimages ; ++i ) {
00802                 proj->read_image(stack_name, list_proj[i]);
00803                 active = proj->get_attr("active");
00804                 if (active) {
00805                         phi   = proj->get_attr("phi");
00806                         theta = proj->get_attr("theta");
00807                         psi   = proj->get_attr("psi");
00808                         r->insert_slice(proj, Transform3D(Ttype,phi,theta,psi));
00809                 }
00810         }       
00811         EMDeletePtr(proj);
00812         return r->finish();
00813 }
00814