#include "mpi.h"
#include "emdata.h"
#include "projector.h"
#include "alignoptions.h"
Include dependency graph for ali3d_d_mpi.h:
This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Defines | |
#define | PI 3.141592653589793 |
Functions | |
std::vector< int > | Numrinit (int first_ring, int last_ring, int skip=1, std::string mode="F") |
std::vector< float > | ringwe (std::vector< int > numr, std::string mode="F") |
int | Applyws (EMData *circ, std::vector< int > numr, std::vector< float > wr) |
EMData * | recons3d_4nn (std::string stack_name, std::vector< int > list_proj, int npad=4) |
int | ali3d_d (MPI_Comm comm, EMData *&volume, EMData **projdata, EMData **cleandata, float *angleshift, int nloc, AlignOptions &options, char *fname_base) |
#define PI 3.141592653589793 |
Definition at line 11 of file ali3d_d_mpi.h.
Referenced by ali3d_d(), LBD_Cart(), Numrinit(), recons3d_CGLS_mpi_Cart(), recons3d_sirt_mpi(), recons3d_sirt_mpi_Cart(), ringwe(), and unified().
int ali3d_d | ( | MPI_Comm | comm, | |
EMData *& | volume, | |||
EMData ** | projdata, | |||
EMData ** | cleandata, | |||
float * | angleshift, | |||
int | nloc, | |||
AlignOptions & | options, | |||
char * | fname_base | |||
) |
Definition at line 20 of file ali3d_d_mpi.cpp.
References anglelist, Applyws(), EMAN::EMData::calc_fourier_shell_correlation(), cb2sph(), cord, EMDeleteArray(), EMDeletePtr(), even_angles(), fcalc(), EMAN::Reconstructor::finish(), Frngs(), EMAN::EMData::get_attr(), get_attr(), AlignOptions::get_CTF(), EMAN::EMData::get_data(), AlignOptions::get_dtheta(), AlignOptions::get_first_ring(), AlignOptions::get_have_angles(), AlignOptions::get_last_ring(), AlignOptions::get_mask3D(), AlignOptions::get_maxit(), EMAN::EMData::get_ndim(), AlignOptions::get_ri(), AlignOptions::get_rstep(), AlignOptions::get_sirt_lam(), AlignOptions::get_sirt_maxit(), AlignOptions::get_sirt_tol(), AlignOptions::get_snr(), AlignOptions::get_step(), AlignOptions::get_symmetry(), AlignOptions::get_use_sirt(), AlignOptions::get_xrng(), EMAN::EMData::get_xsize(), AlignOptions::get_yrng(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), getnnz(), ierr, infomask(), EMAN::Reconstructor::insert_slice(), mul_img(), multiref_polar_ali_2d(), nnz, nrays, numr, Numrinit(), nx, ny, EMAN::EMData::phase_cog(), phi, PI, Polar2Dm(), EMAN::EMData::process_inplace(), EMAN::EMData::project(), ptrs, recons3d_sirt_mpi(), ringwe(), EMAN::EMData::set_attr(), EMAN::EMData::set_attr_dict(), EMAN::EMData::set_size(), setpart(), EMAN::Reconstructor::setup(), theta, weight, wr, and EMAN::EMData::write_image().
Referenced by main().
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 }
int Applyws | ( | EMData * | circ, | |
std::vector< int > | numr, | |||
std::vector< float > | wr | |||
) |
Definition at line 746 of file ali3d_d_mpi.cpp.
References circ.
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 }
std::vector<int> Numrinit | ( | int | first_ring, | |
int | last_ring, | |||
int | skip = 1 , |
|||
std::string | mode = "F" | |||
) |
Definition at line 690 of file ali3d_d_mpi.cpp.
References EMAN::MAXFFT, numr, and PI.
Referenced by ali3d_d().
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 }
EMData* recons3d_4nn | ( | std::string | stack_name, | |
std::vector< int > | list_proj, | |||
int | npad = 4 | |||
) |
Definition at line 772 of file ali3d_d_mpi.cpp.
References EMDeletePtr(), EMAN::Reconstructor::finish(), EMAN::Reconstructor::insert_slice(), phi, proj, EMAN::Reconstructor::setup(), and theta.
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 }
std::vector<float> ringwe | ( | std::vector< int > | numr, | |
std::string | mode = "F" | |||
) |
Definition at line 725 of file ali3d_d_mpi.cpp.
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 }