00001 #include "stdlib.h"
00002 #include "mpi.h"
00003
00004
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
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;
00106
00107 Dict CTF_params;
00108 CTF_params["sign"] = -1.0;
00109 int CTF_applied;
00110 int padffted;
00111 if (CTF) {
00112
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
00124
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
00136 recons_name = "nn4_ctf";
00137 recons_params["snr"] = options.get_snr();
00138 recons_params["sign"] = 1.0;
00139 try {
00140 padffted = projdata[0]->get_attr("padffted");
00141 } catch ( std::exception& e ) {
00142 padffted = 0;
00143 }
00144 if (padffted == 1) recons_params["size"] = (float) recons_params["size"] / (float) recons_params["npad"];
00145 }
00146 }
00147
00148
00149
00150
00151
00152 std::string mode("F");
00153
00154
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
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;
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
00253 float filter_coeff_high = 0.9;
00254 float filter_freq_low = 0.25;
00255 float hcut_default = 0.49;
00256
00257
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
00271 int new_angles, old_angles, new_total, old_total;
00272
00273 for ( int i = 0 ; i < max_iter ; ++i ) {
00274 new_angles = 0;
00275 old_angles = 0;
00276
00277
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);
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
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];
00310
00311 Transform3D R1(EULER_SPIDER , 0.0, 0.0, 0.0);
00312 R1.set_posttrans(Vec3f(best_alignment[1],best_alignment[2],0.0));
00313 R1.set_scale(1.0);
00314
00315 Transform3D R2(EULER_SPIDER, 0.0, 0.0, -best_alignment[0]);
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
00331 Dict composition_dict;
00332
00333 if (best_alignment[3] != 0.0) {
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
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
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
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 {
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 {
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
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 {
00423 ierr = MPI_Send(peak_corr, nloc, MPI_FLOAT, 0, mypid, comm);
00424 }
00425 corr_out.close();
00426
00427 timer = MPI_Wtime();
00428
00429
00430
00431 if (mypid == 0)
00432 printf(" Building even/odd volumes for FSC calculation...\n");
00433 if ( use_sirt ) {
00434
00435
00436 vol1 = new EMData();
00437
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
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
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
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
00543
00544
00545 stats = Util::infomask(vol1, mask3D, false);
00546 *vol1 -= stats[0];
00547 Util::mul_img(vol1, mask3D);
00548 stats = Util::infomask(vol2, mask3D, false);
00549 *vol2 -= stats[0];
00550 Util::mul_img(vol2, mask3D);
00551
00552
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
00560
00561
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
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
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
00631
00632
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
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
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
00655 volume->process_inplace("filter.shift", cog_dict);
00656
00657
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
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);
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
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
00740
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
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
00774
00775 {
00776 Transform3D::EulerType Ttype = Transform3D::SPIDER;
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
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