00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "reconstructor.h"
00037 #include "plugins/reconstructor_template.h"
00038 #include "ctf.h"
00039 #include "emassert.h"
00040 #include "symmetry.h"
00041 #include <cstring>
00042 #include <fstream>
00043 #include <iomanip>
00044 #include <boost/bind.hpp>
00045 #include <boost/format.hpp>
00046
00047 #include <gsl/gsl_statistics_double.h>
00048 #include <gsl/gsl_fit.h>
00049
00050 #ifdef EMAN2_USING_CUDA
00051 #include "cuda/cuda_reconstructor.h"
00052 #endif
00053
00054 using namespace EMAN;
00055 using std::complex;
00056
00057
00058 #include <iostream>
00059 using std::cerr;
00060 using std::endl;
00061 using std::cout;
00062
00063 #include <algorithm>
00064
00065
00066 #include <iomanip>
00067 using std::setprecision;
00068
00069 template < typename T > void checked_delete( T*& x )
00070 {
00071 typedef char type_must_be_complete[ sizeof(T)? 1: -1 ];
00072 (void) sizeof(type_must_be_complete);
00073 delete x;
00074 x = NULL;
00075 }
00076
00077 const string FourierReconstructor::NAME = "fourier";
00078 const string FourierReconstructorSimple2D::NAME = "fouriersimple2D";
00079 const string WienerFourierReconstructor::NAME = "wiener_fourier";
00080 const string BackProjectionReconstructor::NAME = "back_projection";
00081 const string nn4Reconstructor::NAME = "nn4";
00082 const string nn4_rectReconstructor::NAME = "nn4_rect";
00083 const string nnSSNR_Reconstructor::NAME = "nnSSNR";
00084 const string nn4_ctfReconstructor::NAME = "nn4_ctf";
00085 const string nn4_ctf_rectReconstructor::NAME = "nn4_ctf_rect";
00086 const string nnSSNR_ctfReconstructor::NAME = "nnSSNR_ctf";
00087
00088 template <> Factory < Reconstructor >::Factory()
00089 {
00090 force_add<FourierReconstructor>();
00091 force_add<FourierReconstructorSimple2D>();
00092
00093 force_add<WienerFourierReconstructor>();
00094 force_add<BackProjectionReconstructor>();
00095 force_add<nn4Reconstructor>();
00096 force_add<nn4_rectReconstructor>();
00097 force_add<nnSSNR_Reconstructor>();
00098 force_add<nn4_ctfReconstructor>();
00099 force_add<nn4_ctf_rectReconstructor>();
00100 force_add<nnSSNR_ctfReconstructor>();
00101
00102 }
00103
00104 void FourierReconstructorSimple2D::setup()
00105 {
00106 nx = params.set_default("nx",0);
00107
00108 if ( nx < 0 ) throw InvalidValueException(nx, "nx must be positive");
00109
00110 bool is_fftodd = (nx % 2 == 1);
00111
00112 ny = nx;
00113 nx += 2-is_fftodd;
00114
00115 image = new EMData();
00116 image->set_size(nx, ny);
00117 image->set_complex(true);
00118 image->set_fftodd(is_fftodd);
00119 image->set_ri(true);
00120
00121 tmp_data = new EMData();
00122 tmp_data->set_size(nx/2, nx);
00123 }
00124
00125 int FourierReconstructorSimple2D::insert_slice(const EMData* const slice, const Transform & euler, const float)
00126 {
00127
00128
00129 if (!slice) throw NullPointerException("EMData pointer (input image) is NULL");
00130
00131 if ( slice->get_ndim() != 1 ) throw ImageDimensionException("Image dimension must be 1");
00132
00133
00134 if (slice->is_complex()) throw ImageFormatException("The image is complex, expecting real");
00135
00136 EMData* working_slice = slice->process("xform.phaseorigin.tocorner");
00137
00138
00139 working_slice->do_fft_inplace();
00140
00141 float* rdata = image->get_data();
00142 float* norm = tmp_data->get_data();
00143 float* dat = working_slice->get_data();
00144
00145 float g[4];
00146 int offset[4];
00147 float dt[2];
00148 offset[0] = 0; offset[1] = 2; offset[2] = nx; offset[3] = nx+2;
00149
00150 float alt = -((float)(euler.get_rotation("2d"))["alpha"])*M_PI/180.0f;
00151 for (int x = 0; x < working_slice->get_xsize() / 2; x++) {
00152
00153 float rx = (float) x;
00154
00155 float xx = rx*cos(alt);
00156 float yy = rx*sin(alt);
00157 float cc = 1.0;
00158
00159 if (xx < 0) {
00160 xx = -xx;
00161 yy = -yy;
00162 cc = -1.0;
00163 }
00164
00165 yy += ny / 2;
00166
00167
00168 dt[0] = dat[2*x];
00169 dt[1] = cc * dat[2*x+1];
00170
00171
00172 int x0 = (int) floor(xx);
00173 int y0 = (int) floor(yy);
00174
00175 int i = 2*x0 + y0*nx;
00176
00177 float dx = xx - x0;
00178 float dy = yy - y0;
00179
00180 g[0] = Util::agauss(1, dx, dy, 0, EMConsts::I2G);
00181 g[1] = Util::agauss(1, 1 - dx, dy, 0, EMConsts::I2G);
00182 g[2] = Util::agauss(1, dx, 1 - dy, 0, EMConsts::I2G);
00183 g[3] = Util::agauss(1, 1 - dx, 1 - dy, 0, EMConsts::I2G);
00184
00185
00186 if ( x0 == nx-2 ) {
00187 int k = i + offset[0];
00188 rdata[k] += g[0] * dt[0];
00189 rdata[k + 1] += g[0] * dt[1];
00190 norm[k/2] += g[0];
00191
00192 k = i + offset[2];
00193 rdata[k] += g[2] * dt[0];
00194 rdata[k + 1] += g[2] * dt[1];
00195 norm[k/2] += g[2];
00196 continue;
00197
00198 }
00199
00200 if ( x0 > nx-2 ) {
00201 int dif = x0 - (nx-2);
00202 x0 -= dif;
00203 }
00204
00205 if ( y0 == ny -1 ) {
00206 int k = i + offset[0];
00207 rdata[k] += g[0] * dt[0];
00208 rdata[k + 1] += g[0] * dt[1];
00209 norm[k/2] += g[0];
00210
00211 k = i + offset[1];
00212 rdata[k] += g[1] * dt[0];
00213 rdata[k + 1] += g[1] * dt[1];
00214 norm[k/2] += g[1];
00215 continue;
00216 }
00217
00218 if ( y0 > ny-1) {
00219 int dif = y0 - (ny-1);
00220 y0 -= dif;
00221 }
00222
00223 if (x0 >= nx - 2 || y0 >= ny - 1) continue;
00224
00225
00226
00227
00228 for (int j = 0; j < 4; j++)
00229 {
00230 int k = i + offset[j];
00231 rdata[k] += g[j] * dt[0];
00232 rdata[k + 1] += g[j] * dt[1];
00233 norm[k/2] += g[j];
00234
00235 }
00236 }
00237
00238 return 0;
00239
00240 }
00241
00242 EMData *FourierReconstructorSimple2D::finish(bool)
00243 {
00244 normalize_threed();
00245
00246 image->process_inplace("xform.fourierorigin.tocorner");
00247 image->do_ift_inplace();
00248 image->depad();
00249 image->process_inplace("xform.phaseorigin.tocenter");
00250
00251 EMData *ret = image;
00252 image = 0;
00253 return ret;
00254 }
00255
00256 void ReconstructorVolumeData::normalize_threed(const bool sqrt_damp,const bool wiener)
00257
00258 {
00259 float* norm = tmp_data->get_data();
00260 float* rdata = image->get_data();
00261
00262
00263
00264
00265 if ( 0 == norm ) throw NullPointerException("The normalization volume was null!");
00266 if ( 0 == rdata ) throw NullPointerException("The complex reconstruction volume was null!");
00267
00268
00269
00270
00271
00272 float wfilt=0.0;
00273 if (wiener) wfilt=1.0;
00274
00275 for (size_t i = 0; i < (size_t)subnx * subny * subnz; i += 2) {
00276 float d = norm[i/2];
00277 if (sqrt_damp) d*=sqrt(d);
00278 if (d == 0) {
00279 rdata[i] = 0;
00280 rdata[i + 1] = 0;
00281 }
00282 else {
00283
00284
00285 rdata[i] /= d+wfilt;
00286 rdata[i + 1] /= d+wfilt;
00287 }
00288 }
00289
00290
00291 if (subx0==0 && subnx>1 && subny==ny && subnz==nz) {
00292 for (int z=0; z<=nz/2; z++) {
00293 for (int y=1; y<=ny; y++) {
00294 if (y==0 && z==0) continue;
00295
00296 size_t i=(size_t)(y%ny)*subnx+(size_t)(z%nz)*subnx*subny;
00297 size_t i2=(size_t)(ny-y)*subnx+(size_t)((nz-z)%nz)*subnx*subny;
00298 float ar=(rdata[i]+rdata[i2])/2.0f;
00299 float ai=(rdata[i+1]-rdata[i2+1])/2.0f;
00300 rdata[i]=ar;
00301 rdata[i2]=ar;
00302 rdata[i+1]=ai;
00303 rdata[i2+1]=-ai;
00304 }
00305 }
00306 }
00307
00308 if (subx0+subnx==nx && subnx>1 && subny==ny && subnz==nz) {
00309 for (int z=0; z<=nz/2; z++) {
00310 for (int y=1; y<=ny; y++) {
00311 if (y==0 && z==0) continue;
00312
00313 size_t i=(size_t)(y%ny)*subnx+(size_t)(z%nz)*subnx*subny+subnx-2;
00314 size_t i2=(size_t)(ny-y)*subnx+(size_t)((nz-z)%nz)*subnx*subny+subnx-2;
00315 float ar=(rdata[i]+rdata[i2])/2.0f;
00316 float ai=(rdata[i+1]-rdata[i2+1])/2.0f;
00317 rdata[i]=ar;
00318 rdata[i2]=ar;
00319 rdata[i+1]=ai;
00320 rdata[i2+1]=-ai;
00321 }
00322 }
00323 }
00324 }
00325
00326 void FourierReconstructor::load_default_settings()
00327 {
00328 inserter=0;
00329 image=0;
00330 tmp_data=0;
00331 }
00332
00333 void FourierReconstructor::free_memory()
00334 {
00335 if (image) { delete image; image=0; }
00336 if (tmp_data) { delete tmp_data; tmp_data=0; }
00337 if ( inserter != 0 )
00338 {
00339 delete inserter;
00340 inserter = 0;
00341 }
00342 }
00343
00344 #include <sstream>
00345
00346 void FourierReconstructor::load_inserter()
00347 {
00348
00349
00350 Dict parms;
00351 parms["data"] = image;
00352 parms["norm"] = tmp_data->get_data();
00353
00354
00355
00356
00357
00358
00359
00360
00361 if ( inserter != 0 )
00362 {
00363 delete inserter;
00364 }
00365
00366 inserter = Factory<FourierPixelInserter3D>::get((string)params["mode"], parms);
00367 inserter->init();
00368 }
00369
00370 void FourierReconstructor::setup()
00371 {
00372
00373 params.set_default("mode","gauss_2");
00374
00375 vector<int> size=params["size"];
00376
00377 nx = size[0];
00378 ny = size[1];
00379 nz = size[2];
00380 nx2=nx/2-1;
00381 ny2=ny/2;
00382 nz2=nz/2;
00383
00384
00385
00386 bool is_fftodd = (nx % 2 == 1);
00387
00388
00389
00390 nx += 2-is_fftodd;
00391
00392 if (params.has_key("subvolume")) {
00393 vector<int> sub=params["subvolume"];
00394 subx0=sub[0];
00395 suby0=sub[1];
00396 subz0=sub[2];
00397 subnx=sub[3];
00398 subny=sub[4];
00399 subnz=sub[5];
00400
00401 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
00402 throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
00403
00404 }
00405 else {
00406 subx0=suby0=subz0=0;
00407 subnx=nx;
00408 subny=ny;
00409 subnz=nz;
00410 }
00411
00412
00413
00414 if (image) delete image;
00415 image = new EMData();
00416 image->set_size(subnx, subny, subnz);
00417 image->set_complex(true);
00418 image->set_fftodd(is_fftodd);
00419 image->set_ri(true);
00420
00421 image->to_zero();
00422
00423 if (params.has_key("subvolume")) {
00424 image->set_attr("subvolume_x0",subx0);
00425 image->set_attr("subvolume_y0",suby0);
00426 image->set_attr("subvolume_z0",subz0);
00427 image->set_attr("subvolume_full_nx",nx);
00428 image->set_attr("subvolume_full_ny",ny);
00429 image->set_attr("subvolume_full_nz",nz);
00430 }
00431
00432 if (tmp_data) delete tmp_data;
00433 tmp_data = new EMData();
00434 tmp_data->set_size(subnx/2, subny, subnz);
00435 tmp_data->to_zero();
00436 tmp_data->update();
00437
00438 load_inserter();
00439
00440 if ( (bool) params["quiet"] == false )
00441 {
00442 cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
00443 cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
00444 printf ("You will require approximately %1.3g GB of memory to reconstruct this volume\n",((float)subnx*subny*subnz*sizeof(float)*1.5)/1000000000.0);
00445 }
00446 }
00447
00448 void FourierReconstructor::setup_seed(EMData* seed,float seed_weight) {
00449
00450 params.set_default("mode","gauss_2");
00451
00452 vector<int> size=params["size"];
00453
00454 nx = size[0];
00455 ny = size[1];
00456 nz = size[2];
00457 nx2=nx/2-1;
00458 ny2=ny/2;
00459 nz2=nz/2;
00460
00461
00462
00463 bool is_fftodd = (nx % 2 == 1);
00464
00465
00466
00467 nx += 2-is_fftodd;
00468
00469 if (params.has_key("subvolume")) {
00470 vector<int> sub=params["subvolume"];
00471 subx0=sub[0];
00472 suby0=sub[1];
00473 subz0=sub[2];
00474 subnx=sub[3];
00475 subny=sub[4];
00476 subnz=sub[5];
00477
00478 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
00479 throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
00480
00481 }
00482 else {
00483 subx0=suby0=subz0=0;
00484 subnx=nx;
00485 subny=ny;
00486 subnz=nz;
00487 }
00488
00489 if (seed->get_xsize()!=subnx || seed->get_ysize()!=subny || seed->get_zsize()!=subnz || !seed->is_complex())
00490 throw ImageDimensionException("The dimensions of the seed volume do not match the reconstruction size");
00491
00492
00493 image = seed;
00494 if (params.has_key("subvolume")) {
00495 image->set_attr("subvolume_x0",subx0);
00496 image->set_attr("subvolume_y0",suby0);
00497 image->set_attr("subvolume_z0",subz0);
00498 image->set_attr("subvolume_full_nx",nx);
00499 image->set_attr("subvolume_full_ny",ny);
00500 image->set_attr("subvolume_full_nz",nz);
00501 }
00502
00503 if (tmp_data) delete tmp_data;
00504 tmp_data = new EMData();
00505 tmp_data->set_size(subnx/2, subny, subnz);
00506 tmp_data->to_value(seed_weight);
00507
00508 load_inserter();
00509
00510 if ( (bool) params["quiet"] == false )
00511 {
00512 cout << "Seeded direct Fourier inversion";
00513 cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
00514 cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
00515 cout << "You will require approximately " << setprecision(3) << (subnx*subny*subnz*sizeof(float)*1.5)/1000000000.0 << "GB of memory to reconstruct this volume" << endl;
00516 }
00517 }
00518
00519 void FourierReconstructor::clear()
00520 {
00521 bool zeroimage = true;
00522 bool zerotmpimg = true;
00523
00524 #ifdef EMAN2_USING_CUDA
00525 if(EMData::usecuda == 1) {
00526 if(image->getcudarwdata()) {
00527 to_zero_cuda(image->getcudarwdata(),image->get_xsize(),image->get_ysize(),image->get_zsize());
00528 zeroimage = false;
00529 }
00530 if(tmp_data->getcudarwdata()) {
00531
00532 zerotmpimg = false;
00533 }
00534 }
00535 #endif
00536
00537 if(zeroimage) image->to_zero();
00538 if(zerotmpimg) tmp_data->to_zero();
00539
00540 }
00541
00542 EMData* FourierReconstructor::preprocess_slice( const EMData* const slice, const Transform& t )
00543 {
00544 #ifdef EMAN2_USING_CUDA
00545 if(EMData::usecuda == 1) {
00546 if(!slice->getcudarwdata()) slice->copy_to_cuda_keepcpu();
00547 }
00548 #endif
00549
00550 EMData* return_slice = 0;
00551 Transform tmp(t);
00552 tmp.set_rotation(Dict("type","eman"));
00553
00554 if (tmp.is_identity()) return_slice=slice->copy();
00555 else return_slice = slice->process("xform",Dict("transform",&tmp));
00556
00557 return_slice->process_inplace("xform.phaseorigin.tocorner");
00558
00559
00560
00561
00562 #ifdef EMAN2_USING_CUDA
00563 if(return_slice->getcudarwdata()) {
00564 return_slice->do_fft_inplace_cuda();
00565 }else{
00566 return_slice->do_fft_inplace();
00567 }
00568 #else
00569 return_slice->do_fft_inplace();
00570 #endif
00571
00572
00573
00574 return_slice->mult((float)sqrt(1.0f/(return_slice->get_ysize())*return_slice->get_xsize()));
00575
00576
00577
00578
00579 return_slice->set_attr("reconstruct_preproc",(int)1);
00580 return return_slice;
00581 }
00582
00583
00584 int FourierReconstructor::insert_slice(const EMData* const input_slice, const Transform & arg, const float weight)
00585 {
00586
00587 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00588
00589 #ifdef EMAN2_USING_CUDA
00590 if(EMData::usecuda == 1) {
00591 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda_keepcpu();
00592 }
00593 #endif
00594
00595 Transform * rotation;
00596
00597
00598
00599 rotation = new Transform(arg);
00600
00601
00602 EMData *slice;
00603 if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00604 else slice = preprocess_slice( input_slice, *rotation);
00605
00606
00607
00608
00609 rotation->set_scale(1.0);
00610 rotation->set_mirror(false);
00611 rotation->set_trans(0,0,0);
00612
00613
00614
00615 do_insert_slice_work(slice, *rotation, weight);
00616
00617 delete rotation; rotation=0;
00618 delete slice;
00619
00620
00621 return 0;
00622 }
00623
00624 void FourierReconstructor::do_insert_slice_work(const EMData* const input_slice, const Transform & arg,const float weight)
00625 {
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00637
00638 float inx=(float)(input_slice->get_xsize());
00639 float iny=(float)(input_slice->get_ysize());
00640
00641 #ifdef EMAN2_USING_CUDA
00642 if(EMData::usecuda == 1) {
00643 if(!image->getcudarwdata()){
00644 image->copy_to_cuda();
00645 tmp_data->copy_to_cuda();
00646 }
00647 float * m = new float[12];
00648 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00649 Transform t3d = arg*(*it);
00650 t3d.copy_matrix_into_array(m);
00651
00652 insert_slice_cuda(m,input_slice->getcudarwdata(),image->getcudarwdata(),tmp_data->getcudarwdata(),inx,iny,image->get_xsize(),image->get_ysize(),image->get_zsize(), weight);
00653 }
00654 delete m;
00655 return;
00656 }
00657 #endif
00658 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00659 Transform t3d = arg*(*it);
00660 for (int y = -iny/2; y < iny/2; y++) {
00661 for (int x = 0; x <= inx/2; x++) {
00662
00663 float rx = (float) x/(inx-2.0f);
00664 float ry = (float) y/iny;
00665
00666 Vec3f coord(rx,ry,0);
00667 coord = coord*t3d;
00668 float xx = coord[0];
00669 float yy = coord[1];
00670 float zz = coord[2];
00671
00672
00673 xx=xx*(nx-2);
00674 yy=yy*ny;
00675 zz=zz*nz;
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685 inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),weight);
00686 }
00687 }
00688 }
00689 }
00690
00691 int FourierReconstructor::determine_slice_agreement(EMData* input_slice, const Transform & arg, const float weight,bool sub)
00692 {
00693
00694 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00695
00696 #ifdef EMAN2_USING_CUDA
00697 if(EMData::usecuda == 1) {
00698 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda();
00699 }
00700 #endif
00701
00702 Transform * rotation;
00703 rotation = new Transform(arg);
00704
00705 EMData *slice;
00706 if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00707 else slice = preprocess_slice( input_slice, *rotation);
00708
00709
00710
00711
00712 rotation->set_scale(1.0);
00713 rotation->set_mirror(false);
00714 rotation->set_trans(0,0,0);
00715 if (sub) do_insert_slice_work(slice, *rotation, -weight);
00716
00717
00718
00719 do_compare_slice_work(slice, *rotation,weight);
00720
00721 input_slice->set_attr("reconstruct_norm",slice->get_attr("reconstruct_norm"));
00722 input_slice->set_attr("reconstruct_absqual",slice->get_attr("reconstruct_absqual"));
00723
00724 input_slice->set_attr("reconstruct_weight",slice->get_attr("reconstruct_weight"));
00725
00726
00727 if (sub) do_insert_slice_work(slice, *rotation, weight);
00728
00729 delete rotation;
00730 delete slice;
00731
00732
00733 return 0;
00734
00735 }
00736
00737 void FourierReconstructor::do_compare_slice_work(EMData* input_slice, const Transform & arg,float weight)
00738 {
00739
00740 float dt[3];
00741 float dt2[2];
00742 float *dat = input_slice->get_data();
00743 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00744
00745 float inx=(float)(input_slice->get_xsize());
00746 float iny=(float)(input_slice->get_ysize());
00747
00748 double dot=0;
00749 double vweight=0;
00750 double power=0;
00751 double power2=0;
00752 bool use_cpu = true;
00753
00754 #ifdef EMAN2_USING_CUDA
00755 if(EMData::usecuda == 1) {
00756 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda();
00757 if(!image->getcudarwdata()){
00758 image->copy_to_cuda();
00759 tmp_data->copy_to_cuda();
00760 }
00761 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00762 Transform t3d = arg*(*it);
00763 float * m = new float[12];
00764 t3d.copy_matrix_into_array(m);
00765 float4 stats = determine_slice_agreement_cuda(m,input_slice->getcudarwdata(),image->getcudarwdata(),tmp_data->getcudarwdata(),inx,iny,image->get_xsize(),image->get_ysize(),image->get_zsize(), weight);
00766 dot = stats.x;
00767 vweight = stats.y;
00768 power = stats.z;
00769 power2 = stats.w;
00770
00771 use_cpu = false;
00772 }
00773 }
00774 #endif
00775 if(use_cpu) {
00776 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00777 Transform t3d = arg*(*it);
00778 for (int y = -iny/2; y < iny/2; y++) {
00779 for (int x = 0; x <= inx/2; x++) {
00780 if (x==0 && y==0) continue;
00781
00782 float rx = (float) x/(inx-2);
00783 float ry = (float) y/iny;
00784
00785
00786
00787
00788 Vec3f coord(rx,ry,0);
00789 coord = coord*t3d;
00790 float xx = coord[0];
00791 float yy = coord[1];
00792 float zz = coord[2];
00793
00794
00795 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
00796
00797
00798 xx=xx*(nx-2);
00799 yy=yy*ny;
00800 zz=zz*nz;
00801
00802
00803 int idx = (int)(x * 2 + inx*(y<0?iny+y:y));
00804 dt2[0] = dat[idx];
00805 dt2[1] = dat[idx+1];
00806
00807
00808 if (!pixel_at(xx,yy,zz,dt) || dt[2]==0) continue;
00809
00810
00811 dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
00812 vweight+=dt[2];
00813 power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
00814 power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
00815 }
00816 }
00817
00818 }
00819 }
00820
00821 dot/=sqrt(power*power2);
00822
00823 input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)));
00824 input_slice->set_attr("reconstruct_absqual",(float)dot);
00825 float rw=weight<=0?1.0f:1.0f/weight;
00826 input_slice->set_attr("reconstruct_qual",(float)(dot*rw/((rw-1.0)*dot+1.0)));
00827 input_slice->set_attr("reconstruct_weight",(float)vweight/(float)(subnx*subny*subnz));
00828
00829
00830 }
00831
00832 bool FourierReconstructor::pixel_at(const float& xx, const float& yy, const float& zz, float *dt)
00833 {
00834 int x0 = (int) floor(xx);
00835 int y0 = (int) floor(yy);
00836 int z0 = (int) floor(zz);
00837
00838 float *rdata=image->get_data();
00839 float *norm=tmp_data->get_data();
00840 float normsum=0,normsum2=0;
00841
00842 dt[0]=dt[1]=dt[2]=0.0;
00843
00844 if (nx==subnx) {
00845 if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
00846
00847
00848 int x1=x0+1;
00849 int y1=y0+1;
00850 int z1=z0+1;
00851 if (x0<-nx2) x0=-nx2;
00852 if (x1>nx2) x1=nx2;
00853 if (y0<-ny2) y0=-ny2;
00854 if (y1>ny2) y1=ny2;
00855 if (z0<-nz2) z0=-nz2;
00856 if (z1>nz2) z1=nz2;
00857
00858 size_t idx=0;
00859 float r, gg;
00860 for (int k = z0 ; k <= z1; k++) {
00861 for (int j = y0 ; j <= y1; j++) {
00862 for (int i = x0; i <= x1; i ++) {
00863 r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00864 idx=image->get_complex_index_fast(i,j,k);
00865 gg = Util::fast_exp(-r / EMConsts::I2G);
00866
00867 dt[0]+=gg*rdata[idx];
00868 dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
00869 dt[2]+=norm[idx/2]*gg;
00870 normsum2+=gg;
00871 normsum+=gg*norm[idx/2];
00872 }
00873 }
00874 }
00875 if (normsum==0) return false;
00876 dt[0]/=normsum;
00877 dt[1]/=normsum;
00878 dt[2]/=normsum2;
00879
00880 return true;
00881 }
00882 else {
00883 size_t idx;
00884 float r, gg;
00885 for (int k = z0 ; k <= z0 + 1; k++) {
00886 for (int j = y0 ; j <= y0 + 1; j++) {
00887 for (int i = x0; i <= x0 + 1; i ++) {
00888 r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00889 idx=image->get_complex_index(i,j,k,subx0,suby0,subz0,nx,ny,nz);
00890 gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
00891
00892 dt[0]+=gg*rdata[idx];
00893 dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
00894 dt[2]+=norm[idx/2];
00895 normsum+=gg;
00896 }
00897 }
00898 }
00899
00900 if (normsum==0) return false;
00901 return true;
00902 }
00903 }
00904
00905
00906 EMData *FourierReconstructor::finish(bool doift)
00907 {
00908
00909
00910 #ifdef EMAN2_USING_CUDA
00911 if(image->getcudarwdata()){
00912 cout << "copy back from CUDA" << endl;
00913 image->copy_from_device();
00914 tmp_data->copy_from_device();
00915 }
00916 #endif
00917
00918 bool sqrtnorm=params.set_default("sqrtnorm",false);
00919 normalize_threed(sqrtnorm);
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935 if (doift) {
00936 image->do_ift_inplace();
00937 image->depad();
00938 image->process_inplace("xform.phaseorigin.tocenter");
00939 }
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956 image->update();
00957
00958 if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) {
00959 if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
00960 tmp_data->write_image((const char *)params["savenorm"]);
00961 }
00962
00963 delete tmp_data;
00964 tmp_data=0;
00965
00966
00967 EMData *ret=image;
00968 image=0;
00969
00970 return ret;
00971 }
00972
00973 int WienerFourierReconstructor::insert_slice(const EMData* const input_slice, const Transform & arg, const float weight)
00974 {
00975
00976 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00977
00978 Transform * rotation;
00979
00980
00981
00982 rotation = new Transform(arg);
00983
00984
00985 if (!input_slice->has_attr("ctf_snr_total"))
00986 throw NotExistingObjectException("ctf_snr_total","No SNR information present in class-average. Must use the ctf.auto or ctfw.auto averager.");
00987
00988 EMData *slice;
00989 if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00990 else slice = preprocess_slice( input_slice, *rotation);
00991
00992
00993
00994
00995 rotation->set_scale(1.0);
00996 rotation->set_mirror(false);
00997 rotation->set_trans(0,0,0);
00998
00999
01000 do_insert_slice_work(slice, *rotation, weight);
01001
01002 delete rotation; rotation=0;
01003 delete slice;
01004
01005
01006 return 0;
01007 }
01008
01009 void WienerFourierReconstructor::do_insert_slice_work(const EMData* const input_slice, const Transform & arg,const float inweight)
01010 {
01011
01012 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
01013
01014 float inx=(float)(input_slice->get_xsize());
01015 float iny=(float)(input_slice->get_ysize());
01016
01017 int undo_wiener=(int)input_slice->get_attr_default("ctf_wiener_filtered",0);
01018
01019
01020 vector<float> snr=input_slice->get_attr("ctf_snr_total");
01021 float sub=1.0;
01022 if (inweight<0) sub=-1.0;
01023 float weight;
01024
01025 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01026 Transform t3d = arg*(*it);
01027 for (int y = -iny/2; y < iny/2; y++) {
01028 for (int x = 0; x <= inx/2; x++) {
01029
01030 float rx = (float) x/(inx-2.0f);
01031 float ry = (float) y/iny;
01032
01033
01034 float rn = (float)hypot(rx,ry);
01035 if (rn>=.5) continue;
01036 rn*=snr.size()*2.0f;
01037 int rni=(int)floor(rn);
01038 if ((unsigned int)rni>=snr.size()-1) weight=snr[snr.size()-1]*sub;
01039 else {
01040 rn-=rni;
01041 weight=Util::linear_interpolate(snr[rni],snr[rni+1],rn);
01042 }
01043
01044
01045 Vec3f coord(rx,ry,0);
01046 coord = coord*t3d;
01047 float xx = coord[0];
01048 float yy = coord[1];
01049 float zz = coord[2];
01050
01051
01052 xx=xx*(nx-2);
01053 yy=yy*ny;
01054 zz=zz*nz;
01055
01056
01057 if (undo_wiener) inserter->insert_pixel(xx,yy,zz,(input_slice->get_complex_at(x,y))*((weight+1.0f)/weight),weight*sub);
01058 else inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),weight*sub);
01059 }
01060 }
01061 }
01062 }
01063
01064 int WienerFourierReconstructor::determine_slice_agreement(EMData* input_slice, const Transform & arg, const float weight,bool sub)
01065 {
01066
01067 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
01068
01069 Transform * rotation;
01070 rotation = new Transform(arg);
01071
01072 EMData *slice;
01073 if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
01074 else slice = preprocess_slice( input_slice, *rotation);
01075
01076
01077
01078
01079 rotation->set_scale(1.0);
01080 rotation->set_mirror(false);
01081 rotation->set_trans(0,0,0);
01082
01083
01084
01085
01086 if (sub) do_insert_slice_work(slice, *rotation, -weight);
01087
01088
01089 do_compare_slice_work(slice, *rotation,weight);
01090
01091 input_slice->set_attr("reconstruct_norm",slice->get_attr("reconstruct_norm"));
01092 input_slice->set_attr("reconstruct_absqual",slice->get_attr("reconstruct_absqual"));
01093
01094 input_slice->set_attr("reconstruct_weight",slice->get_attr("reconstruct_weight"));
01095
01096
01097 if (sub) do_insert_slice_work(slice, *rotation, weight);
01098
01099
01100 delete rotation;
01101 delete slice;
01102
01103
01104 return 0;
01105
01106 }
01107
01108 void WienerFourierReconstructor::do_compare_slice_work(EMData* input_slice, const Transform & arg,float weight)
01109 {
01110
01111 float dt[3];
01112 float dt2[2];
01113 float *dat = input_slice->get_data();
01114 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
01115
01116 float inx=(float)(input_slice->get_xsize());
01117 float iny=(float)(input_slice->get_ysize());
01118
01119 double dot=0;
01120 double vweight=0;
01121 double power=0;
01122 double power2=0;
01123 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01124 Transform t3d = arg*(*it);
01125 for (int y = -iny/2; y < iny/2; y++) {
01126 for (int x = 0; x <= inx/2; x++) {
01127 if (x==0 && y==0) continue;
01128
01129 float rx = (float) x/(inx-2);
01130 float ry = (float) y/iny;
01131
01132
01133
01134
01135 Vec3f coord(rx,ry,0);
01136 coord = coord*t3d;
01137 float xx = coord[0];
01138 float yy = coord[1];
01139 float zz = coord[2];
01140
01141
01142 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
01143
01144
01145 xx=xx*(nx-2);
01146 yy=yy*ny;
01147 zz=zz*nz;
01148
01149
01150 int idx = (int)(x * 2 + inx*(y<0?iny+y:y));
01151 dt2[0] = dat[idx];
01152 dt2[1] = dat[idx+1];
01153
01154
01155 if (!pixel_at(xx,yy,zz,dt) || dt[2]<=0) continue;
01156
01157
01158 dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
01159 vweight+=dt[2];
01160 power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
01161 power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
01162 }
01163 }
01164 }
01165
01166 dot/=sqrt(power*power2);
01167
01168 input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)));
01169 input_slice->set_attr("reconstruct_absqual",(float)dot);
01170 float rw=weight<=0?1.0f:1.0f/weight;
01171 input_slice->set_attr("reconstruct_qual",(float)(dot*rw/((rw-1.0)*dot+1.0)));
01172 input_slice->set_attr("reconstruct_weight",(float)vweight/(float)(subnx*subny*subnz));
01173
01174
01175 }
01176
01177 bool WienerFourierReconstructor::pixel_at(const float& xx, const float& yy, const float& zz, float *dt)
01178 {
01179 int x0 = (int) floor(xx);
01180 int y0 = (int) floor(yy);
01181 int z0 = (int) floor(zz);
01182
01183 float *rdata=image->get_data();
01184 float *norm=tmp_data->get_data();
01185 float normsum=0,normsum2=0;
01186
01187 dt[0]=dt[1]=dt[2]=0.0;
01188
01189 if (nx==subnx) {
01190 if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
01191
01192
01193 int x1=x0+1;
01194 int y1=y0+1;
01195 int z1=z0+1;
01196 if (x0<-nx2) x0=-nx2;
01197 if (x1>nx2) x1=nx2;
01198 if (y0<-ny2) y0=-ny2;
01199 if (y1>ny2) y1=ny2;
01200 if (z0<-nz2) z0=-nz2;
01201 if (z1>nz2) z1=nz2;
01202
01203 size_t idx=0;
01204 float r, gg;
01205 for (int k = z0 ; k <= z1; k++) {
01206 for (int j = y0 ; j <= y1; j++) {
01207 for (int i = x0; i <= x1; i ++) {
01208 r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
01209 idx=image->get_complex_index_fast(i,j,k);
01210 gg = Util::fast_exp(-r / EMConsts::I2G);
01211
01212 dt[0]+=gg*rdata[idx];
01213 dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
01214 dt[2]+=norm[idx/2]*gg;
01215 normsum2+=gg;
01216 normsum+=gg*norm[idx/2];
01217 }
01218 }
01219 }
01220 if (normsum==0) return false;
01221 dt[0]/=normsum;
01222 dt[1]/=normsum;
01223 dt[2]/=normsum2;
01224
01225 return true;
01226 }
01227 else {
01228 size_t idx;
01229 float r, gg;
01230 for (int k = z0 ; k <= z0 + 1; k++) {
01231 for (int j = y0 ; j <= y0 + 1; j++) {
01232 for (int i = x0; i <= x0 + 1; i ++) {
01233 r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
01234 idx=image->get_complex_index(i,j,k,subx0,suby0,subz0,nx,ny,nz);
01235 gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
01236
01237 dt[0]+=gg*rdata[idx];
01238 dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
01239 dt[2]+=norm[idx/2];
01240 normsum+=gg;
01241 }
01242 }
01243 }
01244
01245 if (normsum==0) return false;
01246 return true;
01247 }
01248 }
01249
01250
01251 EMData *WienerFourierReconstructor::finish(bool doift)
01252 {
01253
01254 bool sqrtnorm=params.set_default("sqrtnorm",false);
01255 normalize_threed(sqrtnorm,true);
01256
01257 if (doift) {
01258 image->do_ift_inplace();
01259 image->depad();
01260 image->process_inplace("xform.phaseorigin.tocenter");
01261 }
01262
01263 image->update();
01264
01265 if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) {
01266 if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
01267 tmp_data->write_image((const char *)params["savenorm"]);
01268 }
01269
01270 delete tmp_data;
01271 tmp_data=0;
01272 EMData *ret=image;
01273 image=0;
01274
01275 return ret;
01276 }
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839 void BackProjectionReconstructor::setup()
01840 {
01841 int size = params["size"];
01842 image = new EMData();
01843 nx = size;
01844 ny = size;
01845 if ( (int) params["zsample"] != 0 ) nz = params["zsample"];
01846 else nz = size;
01847 image->set_size(nx, ny, nz);
01848 }
01849
01850 EMData* BackProjectionReconstructor::preprocess_slice(const EMData* const slice, const Transform& t)
01851 {
01852
01853 EMData* return_slice = slice->process("normalize.edgemean");
01854 return_slice->process_inplace("filter.linearfourier");
01855
01856 Transform tmp(t);
01857 tmp.set_rotation(Dict("type","eman"));
01858 Vec2f trans = tmp.get_trans_2d();
01859 float scale = tmp.get_scale();
01860 bool mirror = tmp.get_mirror();
01861 if (trans[0] != 0 || trans[1] != 0 || scale != 1.0 ) {
01862 return_slice->transform(tmp);
01863 } else if ( mirror == true ) {
01864 return_slice = slice->process("xform.flip",Dict("axis","x"));
01865 }
01866
01867 return return_slice;
01868 }
01869
01870 int BackProjectionReconstructor::insert_slice(const EMData* const input, const Transform &t, const float)
01871 {
01872 if (!input) {
01873 LOGERR("try to insert NULL slice");
01874 return 1;
01875 }
01876
01877 if (input->get_xsize() != input->get_ysize() || input->get_xsize() != nx) {
01878 LOGERR("tried to insert image that was not correction dimensions");
01879 return 1;
01880 }
01881
01882 Transform * transform;
01883 if ( input->has_attr("xform.projection") ) {
01884 transform = (Transform*) (input->get_attr("xform.projection"));
01885 } else {
01886 transform = new Transform(t);
01887 }
01888 EMData* slice = preprocess_slice(input, t);
01889
01890 float weight = params["weight"];
01891 slice->mult(weight);
01892
01893 EMData *tmp = new EMData();
01894 tmp->set_size(nx, ny, nz);
01895
01896 float *slice_data = slice->get_data();
01897 float *tmp_data = tmp->get_data();
01898
01899 size_t nxy = nx * ny;
01900 size_t nxy_size = nxy * sizeof(float);;
01901 for (int i = 0; i < nz; ++i) {
01902 memcpy(&tmp_data[nxy * i], slice_data, nxy_size);
01903 }
01904
01905 transform->set_scale(1.0);
01906 transform->set_mirror(false);
01907 transform->set_trans(0,0,0);
01908 transform->invert();
01909
01910 tmp->transform(*transform);
01911 image->add(*tmp);
01912
01913 if(transform) {delete transform; transform=0;}
01914 delete tmp;
01915 delete slice;
01916
01917 return 0;
01918 }
01919
01920 EMData *BackProjectionReconstructor::finish(bool)
01921 {
01922
01923 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params["sym"]);
01924 vector<Transform> syms = sym->get_syms();
01925
01926 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01927
01928 EMData tmpcopy(*image);
01929 tmpcopy.transform(*it);
01930 image->add(tmpcopy);
01931 }
01932
01933 image->mult(1.0f/(float)sym->get_nsym());
01934 delete sym;
01935
01936 EMData *ret = image;
01937 image = 0 ;
01938 return ret;
01939 }
01940
01941 EMData* EMAN::padfft_slice( const EMData* const slice, const Transform& t, int npad )
01942 {
01943 int nx = slice->get_xsize();
01944 int ny = slice->get_ysize();
01945 int ndim = (ny==1) ? 1 : 2;
01946
01947 if( ndim==2 && nx!=ny )
01948 {
01949
01950 throw std::runtime_error("Tried to padfft a 2D slice which is not square.");
01951 }
01952
01953
01954 EMData* temp = slice->average_circ_sub();
01955
01956 Assert( temp != NULL );
01957 EMData* zeropadded = temp->norm_pad( false, npad );
01958 Assert( zeropadded != NULL );
01959 checked_delete( temp );
01960
01961 zeropadded->do_fft_inplace();
01962 EMData* padfftslice = zeropadded;
01963
01964
01965 Vec2f trans = t.get_trans_2d();
01966 float sx = -trans[0];
01967 float sy = -trans[1];
01968 if(sx != 0.0f || sy != 0.0)
01969 padfftslice->process_inplace("filter.shift", Dict("x_shift", sx, "y_shift", sy, "z_shift", 0.0f));
01970
01971 int remove = slice->get_attr_default("remove", 0);
01972 padfftslice->set_attr( "remove", remove );
01973
01974
01975
01976 padfftslice->center_origin_fft();
01977 return padfftslice;
01978 }
01979
01980 nn4Reconstructor::nn4Reconstructor()
01981 {
01982 m_volume = NULL;
01983 m_wptr = NULL;
01984 }
01985
01986 nn4Reconstructor::nn4Reconstructor( const string& symmetry, int size, int npad )
01987 {
01988 m_volume = NULL;
01989 m_wptr = NULL;
01990 setup( symmetry, size, npad );
01991 load_default_settings();
01992 print_params();
01993 }
01994
01995 nn4Reconstructor::~nn4Reconstructor()
01996 {
01997
01998
01999
02000
02001
02002 }
02003
02004 enum weighting_method { NONE, ESTIMATE, VORONOI };
02005
02006 float max2d( int kc, const vector<float>& pow_a )
02007 {
02008 float max = 0.0;
02009 for( int i=-kc; i <= kc; ++i ) {
02010 for( int j=-kc; j <= kc; ++j ) {
02011 if( i==0 && j==0 ) continue;
02012 {
02013 int c = 2*kc+1 - std::abs(i) - std::abs(j);
02014 max = max + pow_a[c];
02015 }
02016 }
02017 }
02018 return max;
02019 }
02020
02021 float max3d( int kc, const vector<float>& pow_a )
02022 {
02023 float max = 0.0;
02024 for( int i=-kc; i <= kc; ++i ) {
02025 for( int j=-kc; j <= kc; ++j ) {
02026 for( int k=-kc; k <= kc; ++k ) {
02027 if( i==0 && j==0 && k==0 ) continue;
02028
02029 {
02030 int c = 3*kc+1 - std::abs(i) - std::abs(j) - std::abs(k);
02031 max = max + pow_a[c];
02032
02033
02034 }
02035 }
02036 }
02037 }
02038 return max;
02039 }
02040
02041
02042 void nn4Reconstructor::setup()
02043 {
02044 int size = params["size"];
02045 int npad = params["npad"];
02046
02047
02048 string symmetry;
02049 if( params.has_key("symmetry") ) symmetry = params["symmetry"].to_str();
02050 else symmetry = "c1";
02051
02052 if( params.has_key("ndim") ) m_ndim = params["ndim"];
02053 else m_ndim = 3;
02054
02055 if( params.has_key( "snr" ) ) m_osnr = 1.0f/float( params["snr"] );
02056 else m_osnr = 0.0;
02057
02058 setup( symmetry, size, npad );
02059 }
02060
02061 void nn4Reconstructor::setup( const string& symmetry, int size, int npad )
02062 {
02063 m_weighting = ESTIMATE;
02064 m_wghta = 0.2f;
02065
02066 m_symmetry = symmetry;
02067 m_npad = npad;
02068 m_nsym = Transform::get_nsym(m_symmetry);
02069
02070 m_vnx = size;
02071 m_vny = size;
02072 m_vnz = (m_ndim==3) ? size : 1;
02073
02074 m_vnxp = size*npad;
02075 m_vnyp = size*npad;
02076 m_vnzp = (m_ndim==3) ? size*npad : 1;
02077
02078 m_vnxc = m_vnxp/2;
02079 m_vnyc = m_vnyp/2;
02080 m_vnzc = (m_ndim==3) ? m_vnzp/2 : 1;
02081
02082 buildFFTVolume();
02083 buildNormVolume();
02084
02085 }
02086
02087
02088 void nn4Reconstructor::buildFFTVolume() {
02089 int offset = 2 - m_vnxp%2;
02090
02091 m_volume = params["fftvol"];
02092
02093 if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
02094 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
02095 m_volume->to_zero();
02096 }
02097
02098
02099
02100
02101 if ( m_vnxp % 2 == 0 ) m_volume->set_fftodd(0);
02102 else m_volume->set_fftodd(1);
02103
02104
02105 m_volume->set_nxc(m_vnxp/2);
02106 m_volume->set_complex(true);
02107 m_volume->set_ri(true);
02108 m_volume->set_fftpad(true);
02109 m_volume->set_attr("npad", m_npad);
02110 m_volume->set_array_offsets(0,1,1);
02111 }
02112
02113 void nn4Reconstructor::buildNormVolume() {
02114
02115 m_wptr = params["weight"];
02116
02117 if( m_wptr->get_xsize() != m_vnxc+1 &&
02118 m_wptr->get_ysize() != m_vnyp &&
02119 m_wptr->get_zsize() != m_vnzp ) {
02120 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02121 m_wptr->to_zero();
02122 }
02123 m_wptr->set_array_offsets(0,1,1);
02124 }
02125
02126 void printImage( const EMData* line )
02127 {
02128 Assert( line->get_zsize()==1 );
02129
02130
02131 int nx = line->get_xsize();
02132 int ny = line->get_ysize();
02133 for( int j=0; j < ny; ++j ) {
02134 for( int i=0; i < nx; ++i ) printf( "%10.3f ", line->get_value_at(i,j) );
02135 printf( "\n" );
02136 }
02137 }
02138
02139
02140
02141 int nn4Reconstructor::insert_slice(const EMData* const slice, const Transform& t, const float) {
02142
02143 if (!slice) {
02144 LOGERR("try to insert NULL slice");
02145 return 1;
02146 }
02147
02148 int padffted= slice->get_attr_default( "padffted", 0 );
02149 if( m_ndim==3 ) {
02150 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx) ) {
02151
02152 LOGERR("Tried to insert a slice that is the wrong size.");
02153 return 1;
02154 }
02155 } else {
02156 Assert( m_ndim==2 );
02157 if( slice->get_ysize() !=1 ) {
02158 LOGERR( "for 2D reconstruction, a line is excepted" );
02159 return 1;
02160 }
02161 }
02162
02163 EMData* padfft = NULL;
02164
02165 if( padffted != 0 ) padfft = new EMData(*slice);
02166 else padfft = padfft_slice( slice, t, m_npad );
02167
02168 int mult= slice->get_attr_default( "mult", 1 );
02169 Assert( mult > 0 );
02170
02171 if( m_ndim==3 ) {
02172 insert_padfft_slice( padfft, t, mult );
02173 } else {
02174 float alpha = padfft->get_attr( "alpha" );
02175 alpha = alpha/180.0f*M_PI;
02176 for(int i=0; i < m_vnxc+1; ++i ) {
02177 float xnew = i*cos(alpha);
02178 float ynew = -i*sin(alpha);
02179 float btqr = padfft->get_value_at( 2*i, 0, 0 );
02180 float btqi = padfft->get_value_at( 2*i+1, 0, 0 );
02181 if( xnew < 0.0 ) {
02182 xnew *= -1;
02183 ynew *= -1;
02184 btqi *= -1;
02185 }
02186
02187 int ixn = int(xnew+0.5+m_vnxp) - m_vnxp;
02188 int iyn = int(ynew+0.5+m_vnyp) - m_vnyp;
02189
02190 if(iyn < 0 ) iyn += m_vnyp;
02191
02192 (*m_volume)( 2*ixn, iyn+1, 1 ) += btqr *float(mult);
02193 (*m_volume)( 2*ixn+1, iyn+1, 1 ) += btqi * float(mult);
02194 (*m_wptr)(ixn,iyn+1, 1) += float(mult);
02195 }
02196
02197 }
02198 checked_delete( padfft );
02199 return 0;
02200 }
02201
02202 int nn4Reconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
02203 {
02204 Assert( padfft != NULL );
02205
02206 vector<Transform> tsym = t.get_sym_proj(m_symmetry);
02207 for (unsigned int isym=0; isym < tsym.size(); isym++) {
02208 m_volume->nn( m_wptr, padfft, tsym[isym], mult);
02209 }
02210
02211
02212
02213
02214
02215 return 0;
02216 }
02217
02218
02219 #define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
02220
02221 void circumf( EMData* win , int npad)
02222 {
02223 float *tw = win->get_data();
02224
02225
02226 int ix = win->get_xsize();
02227 int iy = win->get_ysize();
02228 int iz = win->get_zsize();
02229 int L2 = (ix/2)*(ix/2);
02230 int L2P = (ix/2-1)*(ix/2-1);
02231
02232 int IP = ix/2+1;
02233 int JP = iy/2+1;
02234 int KP = iz/2+1;
02235
02236
02237 float* sincx = new float[IP+1];
02238 float* sincy = new float[JP+1];
02239 float* sincz = new float[KP+1];
02240
02241 sincx[0] = 1.0f;
02242 sincy[0] = 1.0f;
02243 sincz[0] = 1.0f;
02244
02245 float cdf = M_PI/float(npad*2*ix);
02246 for (int i = 1; i <= IP; ++i) sincx[i] = sin(i*cdf)/(i*cdf);
02247 cdf = M_PI/float(npad*2*iy);
02248 for (int i = 1; i <= JP; ++i) sincy[i] = sin(i*cdf)/(i*cdf);
02249 cdf = M_PI/float(npad*2*iz);
02250 for (int i = 1; i <= KP; ++i) sincz[i] = sin(i*cdf)/(i*cdf);
02251 for (int k = 1; k <= iz; ++k) {
02252 int kkp = abs(k-KP);
02253 for (int j = 1; j <= iy; ++j) {
02254 cdf = sincy[abs(j- JP)]*sincz[kkp];
02255 for (int i = 1; i <= ix; ++i) tw(i,j,k) /= (sincx[abs(i-IP)]*cdf);
02256 }
02257 }
02258
02259 delete[] sincx;
02260 delete[] sincy;
02261 delete[] sincz;
02262
02263 float TNR = 0.0f;
02264 size_t m = 0;
02265 for (int k = 1; k <= iz; ++k) {
02266 for (int j = 1; j <= iy; ++j) {
02267 for (int i = 1; i <= ix; ++i) {
02268 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
02269 if (LR >= (size_t)L2P && LR<=(size_t)L2) {
02270 TNR += tw(i,j,k);
02271 ++m;
02272 }
02273 }
02274 }
02275 }
02276
02277 TNR /=float(m);
02278
02279
02280 for (int k = 1; k <= iz; ++k) {
02281 for (int j = 1; j <= iy; ++j) {
02282 for (int i = 1; i <= ix; ++i) {
02283 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
02284 if (LR<=(size_t)L2) tw(i,j,k) -= TNR;
02285 else tw(i,j,k) = 0.0f;
02286
02287 }
02288 }
02289 }
02290
02291 }
02292
02293
02294 EMData* nn4Reconstructor::finish(bool) {
02295
02296 if( m_ndim==3 ) {
02297 m_volume->symplane0(m_wptr);
02298 } else {
02299 for( int i=1; i <= m_vnyp; ++i ) {
02300
02301 if( (*m_wptr)(0, i, 1)==0.0 ) {
02302 int j = m_vnyp + 1 - i;
02303 (*m_wptr)(0, i, 1) = (*m_wptr)(0, j, 1);
02304 (*m_volume)(0, i, 1) = (*m_volume)(0, j, 1);
02305 (*m_volume)(1, i, 1) = (*m_volume)(1, j, 1);
02306 }
02307 }
02308 }
02309
02310
02311 int box = 7;
02312 int kc = (box-1)/2;
02313 vector< float > pow_a( m_ndim*kc+1, 1.0 );
02314 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
02315 pow_a.back()=0.0;
02316
02317 float alpha = 0.0;
02318 if( m_ndim==3) {
02319 int vol = box*box*box;
02320 float max = max3d( kc, pow_a );
02321 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
02322 } else {
02323 int ara = box*box;
02324 float max = max2d( kc, pow_a );
02325 alpha = ( 1.0f - 1.0f/(float)ara ) / max;
02326 }
02327
02328 int ix,iy,iz;
02329 for (iz = 1; iz <= m_vnzp; iz++) {
02330 for (iy = 1; iy <= m_vnyp; iy++) {
02331 for (ix = 0; ix <= m_vnxc; ix++) {
02332 if ( (*m_wptr)(ix,iy,iz) > 0) {
02333 float tmp;
02334 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+m_osnr);
02335
02336 if( m_weighting == ESTIMATE ) {
02337 int cx = ix;
02338 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
02339 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
02340 float sum = 0.0;
02341 for( int ii = -kc; ii <= kc; ++ii ) {
02342 int nbrcx = cx + ii;
02343 if( nbrcx >= m_vnxc ) continue;
02344 for( int jj= -kc; jj <= kc; ++jj ) {
02345 int nbrcy = cy + jj;
02346 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
02347
02348 int kcz = (m_ndim==3) ? kc : 0;
02349 for( int kk = -kcz; kk <= kcz; ++kk ) {
02350 int nbrcz = cz + kk;
02351 if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
02352 if( nbrcx < 0 ) {
02353 nbrcx = -nbrcx;
02354 nbrcy = -nbrcy;
02355 nbrcz = -nbrcz;
02356 }
02357 int nbrix = nbrcx;
02358 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
02359 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
02360 if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
02361 int c = m_ndim*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
02362 sum = sum + pow_a[c];
02363 }
02364 }
02365 }
02366 }
02367 float wght = 1.0f / ( 1.0f - alpha * sum );
02368 tmp = tmp * wght;
02369 }
02370 (*m_volume)(2*ix,iy,iz) *= tmp;
02371 (*m_volume)(2*ix+1,iy,iz) *= tmp;
02372 }
02373 }
02374 }
02375 }
02376
02377
02378
02379
02380 m_volume->do_ift_inplace();
02381
02382
02383 int npad = m_volume->get_attr("npad");
02384 m_volume->depad();
02385 circumf( m_volume, npad );
02386 m_volume->set_array_offsets( 0, 0, 0 );
02387
02388 return 0;
02389 }
02390 #undef tw
02391
02392
02393 nn4_rectReconstructor::nn4_rectReconstructor()
02394 {
02395 m_volume = NULL;
02396 m_wptr = NULL;
02397 }
02398
02399 nn4_rectReconstructor::nn4_rectReconstructor( const string& symmetry, int size, int npad )
02400 {
02401 m_volume = NULL;
02402 m_wptr = NULL;
02403 setup( symmetry, size, npad );
02404 load_default_settings();
02405 print_params();
02406 }
02407
02408 nn4_rectReconstructor::~nn4_rectReconstructor()
02409 {
02410
02411
02412
02413
02414
02415 }
02416
02417
02418 void nn4_rectReconstructor::setup()
02419 {
02420 m_sizeofprojection = params["sizeprojection"];
02421 int npad = params["npad"];
02422 m_count=0;
02423
02424 string symmetry;
02425 if( params.has_key("symmetry") ) symmetry = params["symmetry"].to_str();
02426 else symmetry = "c1";
02427
02428 if( params.has_key("ndim") ) m_ndim = params["ndim"];
02429 else m_ndim = 3;
02430
02431 if( params.has_key( "snr" ) ) m_osnr = 1.0f/float( params["snr"] );
02432 else m_osnr = 0.0;
02433
02434 setup( symmetry, m_sizeofprojection, npad );
02435 }
02436
02437 void nn4_rectReconstructor::setup( const string& symmetry, int sizeprojection, int npad )
02438 {
02439 m_weighting = ESTIMATE;
02440 m_wghta = 0.2f;
02441 m_symmetry = symmetry;
02442 m_npad = npad;
02443 m_nsym = Transform::get_nsym(m_symmetry);
02444
02445 if( params.has_key("sizex") ) m_vnx = params["sizex"];
02446 else if(params.has_key("xratio"))
02447 {
02448 float temp=params["xratio"];
02449 m_vnx=int(float(sizeprojection)*temp);
02450 }
02451 else m_vnx=sizeprojection;
02452
02453 if( params.has_key("sizey") ) m_vny = params["sizey"];
02454 else if (params.has_key("yratio"))
02455 {
02456 float temp=params["yratio"];
02457 m_vny=int(float(sizeprojection)*temp);
02458 }
02459 else m_vny=sizeprojection;
02460
02461 if( params.has_key("sizez") )
02462 m_vnz = params["sizez"];
02463 else
02464 m_vnz = (m_ndim==3) ? sizeprojection : 1;
02465
02466 m_xratio=float(m_vnx)/float(sizeprojection);
02467 m_yratio=float(m_vny)/float(sizeprojection);
02468
02469 m_vnxp = m_vnx*npad;
02470 m_vnyp = m_vny*npad;
02471 m_vnzp = (m_ndim==3) ? m_vnz*npad : 1;
02472
02473 m_vnxc = m_vnxp/2;
02474 m_vnyc = m_vnyp/2;
02475 m_vnzc = (m_ndim==3) ? m_vnzp/2 : 1;
02476
02477 buildFFTVolume();
02478 buildNormVolume();
02479 }
02480
02481
02482 void nn4_rectReconstructor::buildFFTVolume() {
02483 int offset = 2 - m_vnxp%2;
02484
02485 m_volume = params["fftvol"];
02486
02487 if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
02488 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
02489 m_volume->to_zero();
02490 }
02491
02492
02493
02494
02495 if ( m_vnxp % 2 == 0 ) m_volume->set_fftodd(0);
02496 else m_volume->set_fftodd(1);
02497
02498
02499 m_volume->set_nxc(m_vnxp/2);
02500 m_volume->set_complex(true);
02501 m_volume->set_ri(true);
02502 m_volume->set_fftpad(true);
02503 m_volume->set_attr("npad", m_npad);
02504 m_volume->set_array_offsets(0,1,1);
02505 }
02506
02507 void nn4_rectReconstructor::buildNormVolume() {
02508
02509 m_wptr = params["weight"];
02510
02511 if( m_wptr->get_xsize() != m_vnxc+1 &&
02512 m_wptr->get_ysize() != m_vnyp &&
02513 m_wptr->get_zsize() != m_vnzp ) {
02514 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02515 m_wptr->to_zero();
02516 }
02517
02518 m_wptr->set_array_offsets(0,1,1);
02519 }
02520
02521 int nn4_rectReconstructor::insert_slice(const EMData* const slice, const Transform& t, const float) {
02522
02523
02524
02525
02526 if (!slice) {
02527 LOGERR("try to insert NULL slice");
02528 return 1;
02529 }
02530
02531 int padffted= slice->get_attr_default( "padffted", 0 );
02532 if( m_ndim==3 ) {
02533 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_sizeofprojection) ) {
02534
02535
02536 LOGERR("Tried to insert a slice that is the wrong size.");
02537 return 1;
02538 }
02539 }
02540 if (m_ndim==2) {
02541 Assert( m_ndim==2 );
02542 if( slice->get_ysize() !=1 ) {
02543 LOGERR( "for 2D reconstruction, a line is excepted" );
02544 return 1;
02545 }
02546 }
02547
02548 EMData* padfft = NULL;
02549
02550 if( padffted != 0 ) padfft = new EMData(*slice);
02551 else padfft = padfft_slice( slice, t, m_npad );
02552
02553 int mult= slice->get_attr_default( "mult", 1 );
02554 Assert( mult > 0 );
02555
02556 if( m_ndim==3 ) {
02557 insert_padfft_slice( padfft, t, mult );
02558
02559 } else {
02560 float ellipse_length,ellipse_step,cos_alpha,sin_alpha;
02561 int ellipse_length_int;
02562 float alpha = padfft->get_attr( "alpha" );
02563 alpha = alpha/180.0f*M_PI;
02564 int loop_range;
02565 float temp1,temp2;
02566 temp1=m_xratio*cos(alpha)*float(m_sizeofprojection*m_npad)/2;
02567 temp2=m_yratio*sin(alpha)*float(m_sizeofprojection*m_npad)/2;
02568 ellipse_length=sqrt(temp1*temp1+temp2*temp2);
02569 ellipse_length_int=int(ellipse_length);
02570 ellipse_step=0.5f*(m_sizeofprojection*m_npad)/float(ellipse_length_int);
02571 loop_range=ellipse_length_int;
02572 cos_alpha=temp1/ellipse_length;
02573 sin_alpha=temp2/ellipse_length;
02574 if(m_count%100==0)
02575 {
02576 std::cout<<"#############################################################"<<std::endl;
02577 std::cout<<"line insert start=="<<m_count<<std::endl;
02578 std::cout<<"ellipse lenth=="<<ellipse_length_int<<"ellips step=="<<ellipse_step<<std::endl;
02579 std::cout<<"loop_range"<<loop_range<<std::endl;
02580 std::cout<<"x and y ratio=="<<m_xratio<<" "<<m_yratio<<std::endl;
02581 std::cout<<"cos sin of alpha=="<<cos(alpha)<<" "<<sin(alpha)<<std::endl;
02582 std::cout<<"cos sin of alpha_new==="<<cos_alpha<<sin_alpha<<std::endl;
02583 std::cout<<"alpah dig==="<<cos_alpha<<sin_alpha<<std::endl;
02584 std::cout<<"prjection maximum==="<<loop_range*ellipse_step<<"ideal maximum"<<m_sizeofprojection*m_npad/2<<std::endl;
02585 std::cout<<"x_size=="<<m_volume->get_xsize()<<"y_size=="<<m_volume->get_ysize()<<std::endl;
02586 std::cout<<"#############################################################"<<std::endl;
02587
02588
02589
02590
02591 }
02592 for(int i=0; i <=loop_range; ++i ) {
02593 float xnew = i*cos_alpha;
02594 float ynew = -i*sin_alpha;
02595 if(m_count%100==0&&i==loop_range)
02596 std::cout<<"x_new=="<<xnew<<"Y_new=="<<ynew<<std::endl;
02597 float btqr=0,btqi=0;
02598 float xprj=i*ellipse_step;
02599 float t=xprj-int(xprj);
02600 btqr = (1-t)*padfft->get_value_at( 2*int(xprj), 0, 0 )+t*padfft->get_value_at( 2*(1+int(xprj)), 0, 0 );
02601 btqi = (1-t)*padfft->get_value_at( 2*int(xprj)+1, 0, 0 )+t*padfft->get_value_at( 2*(1+int(xprj))+1, 0, 0 );
02602 if( xnew < 0.0 ) {
02603 xnew *= -1;
02604 ynew *= -1;
02605 btqi *= -1;
02606 }
02607
02608 int ixn = int(xnew+0.5+m_vnxp) - m_vnxp;
02609 int iyn = int(ynew+0.5+m_vnyp) - m_vnyp;
02610
02611 if(iyn < 0 ) iyn += m_vnyp;
02612 if(m_count%100==0&&i==loop_range)
02613 std::cout<<"xnn=="<<ixn<<"ynn=="<<iyn<<std::endl;
02614 (*m_volume)( 2*ixn, iyn+1, 1 ) += btqr *float(mult);
02615 (*m_volume)( 2*ixn+1, iyn+1, 1 ) += btqi * float(mult);
02616 (*m_wptr)(ixn,iyn+1, 1) += float(mult);
02617 }
02618
02619
02620 }
02621 checked_delete( padfft );
02622 return 0;
02623 }
02624
02625
02626
02627
02628 int nn4_rectReconstructor::insert_padfft_slice( EMData* padded, const Transform& t, int mult )
02629 {
02630 Assert( padded != NULL );
02631
02632 vector<Transform> tsym = t.get_sym_proj(m_symmetry);
02633 for (unsigned int isym=0; isym < tsym.size(); isym++) {
02634 m_volume->insert_rect_slice(m_wptr, padded, tsym[isym], m_sizeofprojection, m_xratio, m_yratio, m_npad, mult);
02635 }
02636
02637
02638 return 0;
02639
02640 }
02641
02642 #define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
02643 void circumf_rect( EMData* win , int npad)
02644 {
02645 float *tw = win->get_data();
02646
02647
02648 int ix = win->get_xsize();
02649 int iy = win->get_ysize();
02650 int iz = win->get_zsize();
02651
02652 int IP = ix/2+1;
02653 int JP = iy/2+1;
02654 int KP = iz/2+1;
02655
02656
02657 float* sincx = new float[IP+1];
02658 float* sincy = new float[JP+1];
02659 float* sincz = new float[KP+1];
02660
02661 sincx[0] = 1.0f;
02662 sincy[0] = 1.0f;
02663 sincz[0] = 1.0f;
02664
02665 float cdf = M_PI/float(npad*2*ix);
02666 for (int i = 1; i <= IP; ++i) sincx[i] = sin(i*cdf)/(i*cdf);
02667 cdf = M_PI/float(npad*2*iy);
02668 for (int i = 1; i <= JP; ++i) sincy[i] = sin(i*cdf)/(i*cdf);
02669 cdf = M_PI/float(npad*2*iz);
02670 for (int i = 1; i <= KP; ++i) sincz[i] = sin(i*cdf)/(i*cdf);
02671 for (int k = 1; k <= iz; ++k) {
02672 int kkp = abs(k-KP);
02673 for (int j = 1; j <= iy; ++j) {
02674 cdf = sincy[abs(j- JP)]*sincz[kkp];
02675 for (int i = 1; i <= ix; ++i) tw(i,j,k) /= (sincx[abs(i-IP)]*cdf);
02676 }
02677 }
02678
02679 delete[] sincx;
02680 delete[] sincy;
02681 delete[] sincz;
02682
02683
02684
02685 float dxx = 1.0f/float(0.25*ix*ix);
02686 float dyy = 1.0f/float(0.25*iy*iy);
02687
02688
02689
02690 float LR2=(float(ix)/2-1)*(float(ix)/2-1)*dxx;
02691
02692 float TNR = 0.0f;
02693 size_t m = 0;
02694 for (int k = 1; k <= iz; ++k) {
02695 for (int j = 1; j <= iy; ++j) {
02696 for (int i = 1; i <= ix; ++i) {
02697 float LR = (j-JP)*(j-JP)*dyy+(i-IP)*(i-IP)*dxx;
02698 if (LR<=1.0f && LR >= LR2) {
02699 TNR += tw(i,j,k);
02700 ++m;
02701 }
02702 }
02703 }
02704 }
02705
02706 TNR /=float(m);
02707
02708
02709 for (int k = 1; k <= iz; ++k) {
02710 for (int j = 1; j <= iy; ++j) {
02711 for (int i = 1; i <= ix; ++i) {
02712 float LR = (j-JP)*(j-JP)*dyy+(i-IP)*(i-IP)*dxx;
02713 if (LR<=1.0f) tw(i,j,k)=tw(i,j,k)-TNR;
02714 else tw(i,j,k) = 0.0f;
02715 }
02716 }
02717 }
02718
02719 }
02720 #undef tw
02721
02722 EMData* nn4_rectReconstructor::finish(bool)
02723 {
02724
02725 if( m_ndim==3 ) {
02726 m_volume->symplane0_rect(m_wptr);
02727 } else {
02728 for( int i=1; i <= m_vnyp; ++i ) {
02729
02730 if( (*m_wptr)(0, i, 1)==0.0 ) {
02731 int j = m_vnyp + 1 - i;
02732 (*m_wptr)(0, i, 1) = (*m_wptr)(0, j, 1);
02733 (*m_volume)(0, i, 1) = (*m_volume)(0, j, 1);
02734 (*m_volume)(1, i, 1) = (*m_volume)(1, j, 1);
02735 }
02736 }
02737 }
02738
02739
02740 int box = 7;
02741 int kc = (box-1)/2;
02742 vector< float > pow_a( m_ndim*kc+1, 1.0 );
02743 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
02744 pow_a.back()=0.0;
02745
02746 float alpha = 0.0;
02747 if( m_ndim==3) {
02748 int vol = box*box*box;
02749 float max = max3d( kc, pow_a );
02750 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
02751 } else {
02752 int ara = box*box;
02753 float max = max2d( kc, pow_a );
02754 alpha = ( 1.0f - 1.0f/(float)ara ) / max;
02755 }
02756
02757 int ix,iy,iz;
02758 for (iz = 1; iz <= m_vnzp; iz++) {
02759 for (iy = 1; iy <= m_vnyp; iy++) {
02760 for (ix = 0; ix <= m_vnxc; ix++) {
02761 if ( (*m_wptr)(ix,iy,iz) > 0) {
02762 float tmp;
02763 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+m_osnr);
02764
02765 if( m_weighting == ESTIMATE ) {
02766 int cx = ix;
02767 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
02768 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
02769 float sum = 0.0;
02770 for( int ii = -kc; ii <= kc; ++ii ) {
02771 int nbrcx = cx + ii;
02772 if( nbrcx >= m_vnxc ) continue;
02773 for( int jj= -kc; jj <= kc; ++jj ) {
02774 int nbrcy = cy + jj;
02775 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
02776
02777 int kcz = (m_ndim==3) ? kc : 0;
02778 for( int kk = -kcz; kk <= kcz; ++kk ) {
02779 int nbrcz = cz + kk;
02780 if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
02781 if( nbrcx < 0 ) {
02782 nbrcx = -nbrcx;
02783 nbrcy = -nbrcy;
02784 nbrcz = -nbrcz;
02785 }
02786 int nbrix = nbrcx;
02787 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
02788 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
02789 if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
02790 int c = m_ndim*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
02791 sum = sum + pow_a[c];
02792 }
02793 }
02794 }
02795 }
02796 float wght = 1.0f / ( 1.0f - alpha * sum );
02797 tmp = tmp * wght;
02798 }
02799 (*m_volume)(2*ix,iy,iz) *= tmp;
02800 (*m_volume)(2*ix+1,iy,iz) *= tmp;
02801 }
02802 }
02803 }
02804 }
02805
02806
02807
02808
02809 m_volume->do_ift_inplace();
02810
02811
02812 int npad = m_volume->get_attr("npad");
02813 m_volume->depad();
02814 circumf_rect( m_volume, npad );
02815
02816 m_volume->set_array_offsets( 0, 0, 0 );
02817
02818 return 0;
02819 }
02820
02821
02822
02823
02824
02825
02826 nnSSNR_Reconstructor::nnSSNR_Reconstructor()
02827 {
02828 m_volume = NULL;
02829 m_wptr = NULL;
02830 m_wptr2 = NULL;
02831 }
02832
02833 nnSSNR_Reconstructor::nnSSNR_Reconstructor( const string& symmetry, int size, int npad)
02834 {
02835 m_volume = NULL;
02836 m_wptr = NULL;
02837 m_wptr2 = NULL;
02838
02839 setup( symmetry, size, npad );
02840 }
02841
02842 nnSSNR_Reconstructor::~nnSSNR_Reconstructor()
02843 {
02844
02845
02846
02847
02848
02849
02850
02851 }
02852
02853 void nnSSNR_Reconstructor::setup()
02854 {
02855 int size = params["size"];
02856 int npad = params["npad"];
02857
02858 string symmetry;
02859 if( params.has_key("symmetry") ) symmetry = params["symmetry"].to_str();
02860 else symmetry = "c1";
02861
02862 setup( symmetry, size, npad );
02863 }
02864
02865 void nnSSNR_Reconstructor::setup( const string& symmetry, int size, int npad )
02866 {
02867
02868 m_weighting = ESTIMATE;
02869 m_wghta = 0.2f;
02870 m_wghtb = 0.004f;
02871
02872 m_symmetry = symmetry;
02873 m_npad = npad;
02874 m_nsym = Transform::get_nsym(m_symmetry);
02875
02876 m_vnx = size;
02877 m_vny = size;
02878 m_vnz = size;
02879
02880 m_vnxp = size*npad;
02881 m_vnyp = size*npad;
02882 m_vnzp = size*npad;
02883
02884 m_vnxc = m_vnxp/2;
02885 m_vnyc = m_vnyp/2;
02886 m_vnzc = m_vnzp/2;
02887
02888 buildFFTVolume();
02889 buildNormVolume();
02890 buildNorm2Volume();
02891
02892 }
02893
02894 void nnSSNR_Reconstructor::buildFFTVolume() {
02895
02896 m_volume = params["fftvol"];
02897 m_volume->set_size(m_vnxp+ 2 - m_vnxp%2,m_vnyp,m_vnzp);
02898 m_volume->to_zero();
02899 if ( m_vnxp % 2 == 0 ) m_volume->set_fftodd(0);
02900 else m_volume->set_fftodd(1);
02901
02902 m_volume->set_nxc(m_vnxc);
02903 m_volume->set_complex(true);
02904 m_volume->set_ri(true);
02905 m_volume->set_fftpad(true);
02906 m_volume->set_attr("npad", m_npad);
02907 m_volume->set_array_offsets(0,1,1);
02908 }
02909
02910 void nnSSNR_Reconstructor::buildNormVolume() {
02911
02912 m_wptr = params["weight"];
02913 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02914 m_wptr->to_zero();
02915 m_wptr->set_array_offsets(0,1,1);
02916 }
02917
02918 void nnSSNR_Reconstructor::buildNorm2Volume() {
02919
02920 m_wptr2 = params["weight2"];
02921 m_wptr2->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02922 m_wptr2->to_zero();
02923 m_wptr2->set_array_offsets(0,1,1);
02924 }
02925
02926
02927 int nnSSNR_Reconstructor::insert_slice(const EMData* const slice, const Transform& t, const float) {
02928
02929 if (!slice) {
02930 LOGERR("try to insert NULL slice");
02931 return 1;
02932 }
02933
02934 int padffted=slice->get_attr_default( "padffted", 0 );
02935
02936 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx) ) {
02937
02938 LOGERR("Tried to insert a slice that has wrong size.");
02939 return 1;
02940 }
02941
02942 EMData* padfft = NULL;
02943
02944 if( padffted != 0 ) padfft = new EMData(*slice);
02945 else padfft = padfft_slice( slice, t, m_npad );
02946
02947 int mult = slice->get_attr_default("mult", 1);
02948
02949 Assert( mult > 0 );
02950 insert_padfft_slice( padfft, t, mult );
02951
02952 checked_delete( padfft );
02953 return 0;
02954 }
02955
02956 int nnSSNR_Reconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
02957 {
02958 Assert( padfft != NULL );
02959
02960 for (int isym=0; isym < m_nsym; isym++) {
02961 Transform tsym = t.get_sym(m_symmetry, isym);
02962 m_volume->nn_SSNR( m_wptr, m_wptr2, padfft, tsym, mult);
02963 }
02964 return 0;
02965 }
02966
02967
02968 #define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
02969 EMData* nnSSNR_Reconstructor::finish(bool)
02970 {
02971
02972
02973
02974
02975
02976
02977 int kz, ky;
02978
02979 int box = 7;
02980 int kc = (box-1)/2;
02981 float alpha = 0.0;
02982 float argx, argy, argz;
02983 vector< float > pow_a( 3*kc+1, 1.0 );
02984 float w = params["w"];
02985 EMData* SSNR = params["SSNR"];
02986
02987
02988
02989
02990 EMData* vol_ssnr = params["vol_ssnr"];
02991 vol_ssnr->set_size(m_vnxp+ 2 - m_vnxp%2, m_vnyp ,m_vnzp);
02992 vol_ssnr->to_zero();
02993 if ( m_vnxp % 2 == 0 ) vol_ssnr->set_fftodd(0);
02994 else vol_ssnr->set_fftodd(1);
02995 vol_ssnr->set_nxc(m_vnxc);
02996 vol_ssnr->set_complex(true);
02997 vol_ssnr->set_ri(true);
02998 vol_ssnr->set_fftpad(false);
02999
03000
03001 float dx2 = 1.0f/float(m_vnxc)/float(m_vnxc);
03002 float dy2 = 1.0f/float(m_vnyc)/float(m_vnyc);
03003 #ifdef _WIN32
03004 float dz2 = 1.0f/_cpp_max(float(m_vnzc),1.0f)/_cpp_max(float(m_vnzc),1.0f);
03005 int inc = Util::round(float(_cpp_max(_cpp_max(m_vnxc,m_vnyc),m_vnzc))/w);
03006 #else
03007 float dz2 = 1.0f/std::max(float(m_vnzc),1.0f)/std::max(float(m_vnzc),1.0f);
03008 int inc = Util::round(float(std::max(std::max(m_vnxc,m_vnyc),m_vnzc))/w);
03009 #endif //_WIN32
03010 SSNR->set_size(inc+1,4,1);
03011
03012 float *nom = new float[inc+1];
03013 float *denom = new float[inc+1];
03014 int *nn = new int[inc+1];
03015 int *ka = new int[inc+1];
03016 float wght = 1.0f;
03017 for (int i = 0; i <= inc; i++) {
03018 nom[i] = 0.0f;
03019 denom[i] = 0.0f;
03020 nn[i] = 0;
03021 ka[i] = 0;
03022 }
03023
03024 m_volume->symplane1(m_wptr, m_wptr2);
03025
03026 if ( m_weighting == ESTIMATE ) {
03027 int vol = box*box*box;
03028 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03029 pow_a[3*kc] = 0.0;
03030 float max = max3d( kc, pow_a );
03031 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03032 }
03033
03034 for (int iz = 1; iz <= m_vnzp; iz++) {
03035 if ( iz-1 > m_vnzc ) kz = iz-1-m_vnzp; else kz = iz-1;
03036 argz = float(kz*kz)*dz2;
03037 for (int iy = 1; iy <= m_vnyp; iy++) {
03038 if ( iy-1 > m_vnyc ) ky = iy-1-m_vnyp; else ky = iy-1;
03039 argy = argz + float(ky*ky)*dy2;
03040 for (int ix = 0; ix <= m_vnxc; ix++) {
03041 float Kn = (*m_wptr)(ix,iy,iz);
03042 argx = std::sqrt(argy + float(ix*ix)*dx2);
03043 int r = Util::round(float(inc)*argx);
03044 if ( r >= 0 && Kn > 4.5f ) {
03045 if ( m_weighting == ESTIMATE ) {
03046 int cx = ix;
03047 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03048 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03049
03050 float sum = 0.0;
03051 for( int ii = -kc; ii <= kc; ++ii ) {
03052 int nbrcx = cx + ii;
03053 if( nbrcx >= m_vnxc ) continue;
03054 for ( int jj= -kc; jj <= kc; ++jj ) {
03055 int nbrcy = cy + jj;
03056 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03057 for( int kk = -kc; kk <= kc; ++kk ) {
03058 int nbrcz = cz + jj;
03059 if ( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
03060 if( nbrcx < 0 ) {
03061 nbrcx = -nbrcx;
03062 nbrcy = -nbrcy;
03063 nbrcz = -nbrcz;
03064 }
03065 int nbrix = nbrcx;
03066 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
03067 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
03068 if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
03069 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
03070 sum = sum + pow_a[c];
03071 }
03072 }
03073 }
03074 }
03075 wght = 1.0f / ( 1.0f - alpha * sum );
03076 }
03077 float nominator = std::norm(m_volume->cmplx(ix,iy,iz))/Kn;
03078 float denominator = ((*m_wptr2)(ix,iy,iz)-nominator)/(Kn-1.0f);
03079
03080 if( (ix>0 || (kz>=0 && (ky>=0 || kz!=0)))) {
03081 if ( r <= inc ) {
03082 nom[r] += nominator*wght;
03083 denom[r] += denominator/Kn*wght;
03084 nn[r] += 2;
03085 ka[r] += int(Kn);
03086 }
03087
03088
03089
03090
03091
03092
03093
03094
03095
03096
03097
03098
03099
03100
03101
03102
03103 }
03104 (*vol_ssnr)(2*ix, iy-1, iz-1) = denominator*wght;
03105
03106
03107 }
03108 }
03109 }
03110 }
03111
03112 for (int i = 0; i <= inc; i++) {
03113 (*SSNR)(i,0,0) = nom[i];
03114 (*SSNR)(i,1,0) = denom[i];
03115 (*SSNR)(i,2,0) = static_cast<float>(nn[i]);
03116 (*SSNR)(i,3,0) = static_cast<float>(ka[i]);
03117 }
03118 vol_ssnr->update();
03119
03120 delete[] nom;
03121 delete[] denom;
03122 delete[] nn;
03123 delete[] ka;
03124
03125 return 0;
03126 }
03127 #undef tw
03128
03129
03130
03131
03132
03133
03134
03135 nn4_ctfReconstructor::nn4_ctfReconstructor()
03136 {
03137 m_volume = NULL;
03138 m_wptr = NULL;
03139 }
03140
03141 nn4_ctfReconstructor::nn4_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign )
03142 {
03143 setup( symmetry, size, npad, snr, sign );
03144 }
03145
03146 nn4_ctfReconstructor::~nn4_ctfReconstructor()
03147 {
03148
03149
03150
03151
03152
03153 }
03154
03155 void nn4_ctfReconstructor::setup()
03156 {
03157 if( ! params.has_key("size") ) throw std::logic_error("Error: image size is not given");
03158
03159 int size = params["size"];
03160 int npad = params.has_key("npad") ? int(params["npad"]) : 4;
03161
03162 int sign = 1;
03163 string symmetry = params.has_key("symmetry")? params["symmetry"].to_str() : "c1";
03164
03165 float snr = params["snr"];
03166
03167 m_varsnr = params.has_key("varsnr") ? int(params["varsnr"]) : 0;
03168 setup( symmetry, size, npad, snr, sign );
03169
03170 }
03171
03172 void nn4_ctfReconstructor::setup( const string& symmetry, int size, int npad, float snr, int sign )
03173 {
03174 m_weighting = ESTIMATE;
03175 if( params.has_key("weighting") ) {
03176 int tmp = int( params["weighting"] );
03177 if( tmp==0 ) m_weighting = NONE;
03178 }
03179
03180
03181
03182 m_wghta = 0.2f;
03183 m_wghtb = 0.004f;
03184
03185 m_symmetry = symmetry;
03186 m_npad = npad;
03187 m_sign = sign;
03188 m_nsym = Transform::get_nsym(m_symmetry);
03189
03190 m_snr = snr;
03191
03192 m_vnx = size;
03193 m_vny = size;
03194 m_vnz = size;
03195
03196 m_vnxp = size*npad;
03197 m_vnyp = size*npad;
03198 m_vnzp = size*npad;
03199
03200 m_vnxc = m_vnxp/2;
03201 m_vnyc = m_vnyp/2;
03202 m_vnzc = m_vnzp/2;
03203
03204 buildFFTVolume();
03205 buildNormVolume();
03206 }
03207
03208 void nn4_ctfReconstructor::buildFFTVolume() {
03209 int offset = 2 - m_vnxp%2;
03210
03211 m_volume = params["fftvol"];
03212
03213 if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
03214 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
03215 m_volume->to_zero();
03216 }
03217
03218 m_volume->set_nxc(m_vnxp/2);
03219 m_volume->set_complex(true);
03220 m_volume->set_ri(true);
03221 m_volume->set_fftpad(true);
03222 m_volume->set_attr("npad", m_npad);
03223 m_volume->set_array_offsets(0,1,1);
03224 }
03225
03226 void nn4_ctfReconstructor::buildNormVolume()
03227 {
03228 m_wptr = params["weight"];
03229
03230 if( m_wptr->get_xsize() != m_vnxc+1 && m_wptr->get_ysize() != m_vnyp && m_wptr->get_zsize() != m_vnzp ) {
03231 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03232 m_wptr->to_zero();
03233 }
03234
03235 m_wptr->set_array_offsets(0,1,1);
03236
03237 }
03238
03239 int nn4_ctfReconstructor::insert_slice(const EMData* const slice, const Transform& t, const float)
03240 {
03241
03242 if (!slice) {
03243 LOGERR("try to insert NULL slice");
03244 return 1;
03245 }
03246
03247 int buffed = slice->get_attr_default( "buffed", 0 );
03248 if( buffed > 0 ) {
03249 int mult = slice->get_attr_default( "mult", 1 );
03250 insert_buffed_slice( slice, mult );
03251 return 0;
03252 }
03253
03254 int padffted= slice->get_attr_default("padffted", 0);
03255 if( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx) )
03256 {
03257
03258 LOGERR("Tried to insert a slice that is the wrong size.");
03259 return 1;
03260 }
03261
03262 EMData* padfft = NULL;
03263
03264 if( padffted != 0 ) padfft = new EMData(*slice);
03265 else padfft = padfft_slice( slice, t, m_npad );
03266
03267 int mult= slice->get_attr_default("mult", 1);
03268
03269 Assert( mult > 0 );
03270 insert_padfft_slice( padfft, t, mult );
03271
03272 checked_delete( padfft );
03273
03274 return 0;
03275 }
03276
03277 int nn4_ctfReconstructor::insert_buffed_slice( const EMData* buffed, int mult )
03278 {
03279 const float* bufdata = buffed->get_data();
03280 float* cdata = m_volume->get_data();
03281 float* wdata = m_wptr->get_data();
03282
03283 int npoint = buffed->get_xsize()/4;
03284 for( int i=0; i < npoint; ++i ) {
03285
03286 int pos2 = int( bufdata[4*i] );
03287 int pos1 = pos2 * 2;
03288 cdata[pos1 ] += bufdata[4*i+1]*mult;
03289 cdata[pos1+1] += bufdata[4*i+2]*mult;
03290 wdata[pos2 ] += bufdata[4*i+3]*mult;
03291
03292
03293
03294
03295
03296 }
03297 return 0;
03298 }
03299
03300 int nn4_ctfReconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
03301 {
03302 Assert( padfft != NULL );
03303 float tmp = padfft->get_attr("ctf_applied");
03304 int ctf_applied = (int) tmp;
03305
03306 vector<Transform> tsym = t.get_sym_proj(m_symmetry);
03307 for (unsigned int isym=0; isym < tsym.size(); isym++) {
03308 if(ctf_applied) m_volume->nn_ctf_applied(m_wptr, padfft, tsym[isym], mult);
03309
03310 else m_volume->nn_ctf(m_wptr, padfft, tsym[isym], mult);
03311
03312 }
03313
03314
03315 return 0;
03316
03317 }
03318
03319 #define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
03320 EMData* nn4_ctfReconstructor::finish(bool)
03321 {
03322 m_volume->set_array_offsets(0, 1, 1);
03323 m_wptr->set_array_offsets(0, 1, 1);
03324 m_volume->symplane0_ctf(m_wptr);
03325
03326 int box = 7;
03327 int vol = box*box*box;
03328 int kc = (box-1)/2;
03329 vector< float > pow_a( 3*kc+1, 1.0 );
03330 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03331 pow_a[3*kc]=0.0;
03332
03333
03334 float max = max3d( kc, pow_a );
03335 float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03336 float osnr = 1.0f/m_snr;
03337
03338
03339 int ix,iy,iz;
03340 for (iz = 1; iz <= m_vnzp; iz++) {
03341 for (iy = 1; iy <= m_vnyp; iy++) {
03342 for (ix = 0; ix <= m_vnxc; ix++) {
03343 if ( (*m_wptr)(ix,iy,iz) > 0.0f) {
03344 int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
03345 int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
03346 float tmp=0.0;
03347 if( m_varsnr )
03348 {
03349 float freq = sqrt( (float)(ix*ix+iyp*iyp+izp*izp) );
03350 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+freq*osnr)*m_sign;
03351 }
03352 else
03353 {
03354 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+osnr)*m_sign;
03355 }
03356
03357 if( m_weighting == ESTIMATE ) {
03358 int cx = ix;
03359 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03360 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03361 float sum = 0.0;
03362 for( int ii = -kc; ii <= kc; ++ii ) {
03363 int nbrcx = cx + ii;
03364 if( nbrcx >= m_vnxc ) continue;
03365 for( int jj= -kc; jj <= kc; ++jj ) {
03366 int nbrcy = cy + jj;
03367 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03368 for( int kk = -kc; kk <= kc; ++kk ) {
03369 int nbrcz = cz + jj;
03370 if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
03371 if( nbrcx < 0 ) {
03372 nbrcx = -nbrcx;
03373 nbrcy = -nbrcy;
03374 nbrcz = -nbrcz;
03375 }
03376
03377 int nbrix = nbrcx;
03378 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
03379 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
03380 if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0.0 ) {
03381 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
03382 sum = sum + pow_a[c];
03383
03384
03385 }
03386 }
03387 }
03388 }
03389 float wght = 1.0f / ( 1.0f - alpha * sum );
03390
03391
03392
03393
03394
03395
03396
03397
03398
03399 tmp = tmp * wght;
03400 }
03401 (*m_volume)(2*ix,iy,iz) *= tmp;
03402 (*m_volume)(2*ix+1,iy,iz) *= tmp;
03403 }
03404 }
03405 }
03406 }
03407
03408
03409 m_volume->do_ift_inplace();
03410 int npad = m_volume->get_attr("npad");
03411 m_volume->depad();
03412 circumf( m_volume, npad );
03413 m_volume->set_array_offsets( 0, 0, 0 );
03414
03415 return 0;
03416 }
03417 #undef tw
03418
03419
03420
03421
03422
03423
03424 nn4_ctf_rectReconstructor::nn4_ctf_rectReconstructor()
03425 {
03426 m_volume = NULL;
03427 m_wptr = NULL;
03428 }
03429
03430 nn4_ctf_rectReconstructor::nn4_ctf_rectReconstructor( const string& symmetry, int size, int npad, float snr, int sign )
03431 {
03432 setup( symmetry, size, npad, snr, sign );
03433 }
03434
03435 nn4_ctf_rectReconstructor::~nn4_ctf_rectReconstructor()
03436 {
03437
03438
03439
03440
03441
03442 }
03443
03444 void nn4_ctf_rectReconstructor::setup()
03445 {
03446 if( ! params.has_key("sizeprojection") ) throw std::logic_error("Error: projection size is not given");
03447 m_sizeofprojection = params["sizeprojection"];
03448 int npad = params.has_key("npad") ? int(params["npad"]) : 4;
03449
03450 int sign = 1;
03451 string symmetry = params.has_key("symmetry")? params["symmetry"].to_str() : "c1";
03452
03453 float snr = params["snr"];
03454
03455 m_varsnr = params.has_key("varsnr") ? int(params["varsnr"]) : 0;
03456 setup( symmetry, m_sizeofprojection, npad, snr, sign );
03457
03458 }
03459
03460 void nn4_ctf_rectReconstructor::setup( const string& symmetry, int sizeprojection, int npad, float snr, int sign )
03461 {
03462 m_weighting = ESTIMATE;
03463 if( params.has_key("weighting") ) {
03464 int tmp = int( params["weighting"] );
03465 if( tmp==0 ) m_weighting = NONE;
03466 }
03467
03468 m_wghta = 0.2f;
03469 m_wghtb = 0.004f;
03470
03471 m_symmetry = symmetry;
03472 m_npad = npad;
03473 m_sign = sign;
03474 m_nsym = Transform::get_nsym(m_symmetry);
03475
03476 m_snr = snr;
03477 if (params.has_key("sizex")) m_vnx = params["sizex"];
03478 else if (params.has_key("xratio")) {
03479 float temp=params["xratio"];
03480 m_vnx=int(float(sizeprojection)*temp);
03481 } else m_vnx=sizeprojection;
03482
03483 if (params.has_key("sizey")) m_vny = params["sizey"];
03484 else if (params.has_key("yratio")) {
03485 float temp=params["yratio"];
03486 m_vny=int(float(sizeprojection)*temp);
03487 }
03488 else m_vny=sizeprojection;
03489
03490 if( params.has_key("sizez") )
03491 m_vnz = params["sizez"];
03492 else
03493 m_vnz =sizeprojection;
03494
03495 m_xratio=float(m_vnx)/float(sizeprojection);
03496 m_yratio=float(m_vny)/float(sizeprojection);
03497
03498
03499
03500
03501 m_vnxp = m_vnx*npad;
03502 m_vnyp = m_vny*npad;
03503 m_vnzp = m_vnz*npad;
03504
03505 m_vnxc = m_vnxp/2;
03506 m_vnyc = m_vnyp/2;
03507 m_vnzc = m_vnzp/2;
03508
03509 buildFFTVolume();
03510 buildNormVolume();
03511 }
03512
03513 void nn4_ctf_rectReconstructor::buildFFTVolume() {
03514 int offset = 2 - m_vnxp%2;
03515
03516 m_volume = params["fftvol"];
03517
03518 if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
03519 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
03520 m_volume->to_zero();
03521 }
03522
03523 m_volume->set_nxc(m_vnxp/2);
03524 m_volume->set_complex(true);
03525 m_volume->set_ri(true);
03526 m_volume->set_fftpad(true);
03527 m_volume->set_attr("npad", m_npad);
03528 m_volume->set_array_offsets(0,1,1);
03529 }
03530
03531 void nn4_ctf_rectReconstructor::buildNormVolume()
03532 {
03533 m_wptr = params["weight"];
03534
03535 if( m_wptr->get_xsize() != m_vnxc+1 && m_wptr->get_ysize() != m_vnyp && m_wptr->get_zsize() != m_vnzp ) {
03536 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03537 m_wptr->to_zero();
03538 }
03539
03540 m_wptr->set_array_offsets(0,1,1);
03541
03542 }
03543
03544 int nn4_ctf_rectReconstructor::insert_slice(const EMData* const slice, const Transform& t, const float)
03545 {
03546
03547 if (!slice) {
03548 LOGERR("try to insert NULL slice");
03549 return 1;
03550 }
03551
03552 int buffed = slice->get_attr_default( "buffed", 0 );
03553 if( buffed > 0 ) {
03554 int mult = slice->get_attr_default( "mult", 1 );
03555 insert_buffed_slice( slice, mult );
03556 return 0;
03557 }
03558
03559 int padffted= slice->get_attr_default("padffted", 0);
03560
03561 if( padffted==0 && (slice->get_xsize()!=slice->get_ysize()) )
03562 {
03563
03564 LOGERR("Tried to insert a slice that is the wrong size.");
03565 return 1;
03566 }
03567
03568 EMData* padfft = NULL;
03569
03570 if( padffted != 0 ) padfft = new EMData(*slice);
03571 else padfft = padfft_slice( slice, t, m_npad );
03572
03573 int mult= slice->get_attr_default("mult", 1);
03574
03575 Assert( mult > 0 );
03576 insert_padfft_slice( padfft, t, mult );
03577
03578 checked_delete( padfft );
03579
03580 return 0;
03581 }
03582
03583 int nn4_ctf_rectReconstructor::insert_buffed_slice( const EMData* buffed, int mult )
03584 {
03585 const float* bufdata = buffed->get_data();
03586 float* cdata = m_volume->get_data();
03587 float* wdata = m_wptr->get_data();
03588
03589 int npoint = buffed->get_xsize()/4;
03590 for( int i=0; i < npoint; ++i ) {
03591
03592 int pos2 = int( bufdata[4*i] );
03593 int pos1 = pos2 * 2;
03594 cdata[pos1 ] += bufdata[4*i+1]*mult;
03595 cdata[pos1+1] += bufdata[4*i+2]*mult;
03596 wdata[pos2 ] += bufdata[4*i+3]*mult;
03597
03598
03599
03600
03601
03602 }
03603 return 0;
03604 }
03605
03606
03607 int nn4_ctf_rectReconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
03608 {
03609 Assert( padfft != NULL );
03610 float tmp = padfft->get_attr("ctf_applied");
03611 int ctf_applied = (int) tmp;
03612 vector<Transform> tsym = t.get_sym_proj(m_symmetry);
03613 for (unsigned int isym=0; isym < tsym.size(); isym++) {
03614 if(ctf_applied) m_volume->insert_rect_slice_ctf_applied(m_wptr, padfft, tsym[isym], m_sizeofprojection, m_xratio,m_yratio, m_npad, mult);
03615
03616 else m_volume->insert_rect_slice_ctf(m_wptr, padfft, tsym[isym], m_sizeofprojection, m_xratio, m_yratio, m_npad, mult);
03617
03618 }
03619
03620
03621 return 0;
03622
03623 }
03624
03625 #define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
03626 EMData* nn4_ctf_rectReconstructor::finish(bool)
03627 {
03628 m_volume->set_array_offsets(0, 1, 1);
03629 m_wptr->set_array_offsets(0, 1, 1);
03630 m_volume->symplane0_rect(m_wptr);
03631
03632 int box = 7;
03633 int vol = box*box*box;
03634 int kc = (box-1)/2;
03635 vector< float > pow_a( 3*kc+1, 1.0 );
03636 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03637 pow_a[3*kc]=0.0;
03638
03639
03640 float max = max3d( kc, pow_a );
03641 float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03642 float osnr = 1.0f/m_snr;
03643
03644
03645 int ix,iy,iz;
03646 for (iz = 1; iz <= m_vnzp; iz++) {
03647 for (iy = 1; iy <= m_vnyp; iy++) {
03648 for (ix = 0; ix <= m_vnxc; ix++) {
03649 if ( (*m_wptr)(ix,iy,iz) > 0.0f) {
03650 int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
03651 int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
03652 float tmp=0.0;
03653 if( m_varsnr )
03654 {
03655 float freq = sqrt( (float)(ix*ix/(m_xratio*m_xratio)+iyp*iyp/(m_zratio*m_yratio)+izp*izp) );
03656 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+freq*osnr)*m_sign;
03657 }
03658 else
03659 {
03660 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+osnr)*m_sign;
03661 }
03662
03663 if( m_weighting == ESTIMATE ) {
03664 int cx = ix;
03665 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03666 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03667 float sum = 0.0;
03668 for( int ii = -kc; ii <= kc; ++ii ) {
03669 int nbrcx = cx + ii;
03670 if( nbrcx >= m_vnxc ) continue;
03671 for( int jj= -kc; jj <= kc; ++jj ) {
03672 int nbrcy = cy + jj;
03673 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03674 for( int kk = -kc; kk <= kc; ++kk ) {
03675 int nbrcz = cz + jj;
03676 if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
03677 if( nbrcx < 0 ) {
03678 nbrcx = -nbrcx;
03679 nbrcy = -nbrcy;
03680 nbrcz = -nbrcz;
03681 }
03682
03683 int nbrix = nbrcx;
03684 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
03685 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
03686 if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0.0 ) {
03687 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
03688 sum = sum + pow_a[c];
03689
03690
03691 }
03692 }
03693 }
03694 }
03695 float wght = 1.0f / ( 1.0f - alpha * sum );
03696
03697
03698
03699
03700
03701
03702
03703
03704
03705 tmp = tmp * wght;
03706 }
03707 (*m_volume)(2*ix,iy,iz) *= tmp;
03708 (*m_volume)(2*ix+1,iy,iz) *= tmp;
03709 }
03710 }
03711 }
03712 }
03713
03714
03715 m_volume->do_ift_inplace();
03716 int npad = m_volume->get_attr("npad");
03717 m_volume->depad();
03718 circumf_rect( m_volume, npad );
03719 m_volume->set_array_offsets( 0, 0, 0 );
03720 return 0;
03721 }
03722 #undef tw
03723
03724
03725
03726
03727
03728
03729
03730 nnSSNR_ctfReconstructor::nnSSNR_ctfReconstructor()
03731 {
03732 m_volume = NULL;
03733 m_wptr = NULL;
03734 m_wptr2 = NULL;
03735 m_wptr3 = NULL;
03736 }
03737
03738 nnSSNR_ctfReconstructor::nnSSNR_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign)
03739 {
03740 m_volume = NULL;
03741 m_wptr = NULL;
03742 m_wptr2 = NULL;
03743 m_wptr3 = NULL;
03744
03745 setup( symmetry, size, npad, snr, sign );
03746 }
03747
03748 nnSSNR_ctfReconstructor::~nnSSNR_ctfReconstructor()
03749 {
03750
03751
03752
03753
03754
03755
03756 }
03757
03758 void nnSSNR_ctfReconstructor::setup()
03759 {
03760 int size = params["size"];
03761 int npad = params["npad"];
03762 int sign = params["sign"];
03763 float snr = params["snr"];
03764 string symmetry;
03765 if( params.has_key("symmetry") ) symmetry = params["symmetry"].to_str();
03766 else symmetry = "c1";
03767
03768 setup( symmetry, size, npad, snr, sign );
03769 }
03770 void nnSSNR_ctfReconstructor::setup( const string& symmetry, int size, int npad, float snr, int sign )
03771 {
03772
03773 m_weighting = ESTIMATE;
03774 m_wghta = 0.2f;
03775 m_wghtb = 0.004f;
03776 wiener = 1;
03777
03778 m_symmetry = symmetry;
03779 m_npad = npad;
03780 m_nsym = Transform::get_nsym(m_symmetry);
03781
03782 m_sign = sign;
03783 m_snr = snr;
03784
03785 m_vnx = size;
03786 m_vny = size;
03787 m_vnz = size;
03788
03789 m_vnxp = size*npad;
03790 m_vnyp = size*npad;
03791 m_vnzp = size*npad;
03792
03793 m_vnxc = m_vnxp/2;
03794 m_vnyc = m_vnyp/2;
03795 m_vnzc = m_vnzp/2;
03796
03797 buildFFTVolume();
03798 buildNormVolume();
03799 buildNorm2Volume();
03800 buildNorm3Volume();
03801 }
03802
03803 void nnSSNR_ctfReconstructor::buildFFTVolume() {
03804
03805 int offset = 2 - m_vnxp%2;
03806 m_volume = params["fftvol"];
03807
03808 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
03809 m_volume->to_zero();
03810
03811 m_volume->set_fftodd(m_vnxp % 2);
03812
03813 m_volume->set_nxc(m_vnxc);
03814 m_volume->set_complex(true);
03815 m_volume->set_ri(true);
03816 m_volume->set_fftpad(true);
03817 m_volume->set_attr("npad", m_npad);
03818 m_volume->set_array_offsets(0,1,1);
03819 }
03820
03821
03822 void nnSSNR_ctfReconstructor::buildNormVolume()
03823 {
03824 m_wptr = params["weight"];
03825 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03826 m_wptr->to_zero();
03827 m_wptr->set_array_offsets(0,1,1);
03828 }
03829
03830 void nnSSNR_ctfReconstructor::buildNorm2Volume() {
03831
03832 m_wptr2 = params["weight2"];
03833 m_wptr2->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03834 m_wptr2->to_zero();
03835 m_wptr2->set_array_offsets(0,1,1);
03836 }
03837
03838 void nnSSNR_ctfReconstructor::buildNorm3Volume() {
03839
03840 m_wptr3 = params["weight3"];
03841 m_wptr3->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03842 m_wptr3->to_zero();
03843 m_wptr3->set_array_offsets(0,1,1);
03844 }
03845
03846 int nnSSNR_ctfReconstructor::insert_slice(const EMData *const slice, const Transform& t, const float) {
03847
03848 if (!slice) {
03849 LOGERR("try to insert NULL slice");
03850 return 1;
03851 }
03852 int padffted= slice->get_attr_default("padffted", 0);
03853 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx) ) {
03854
03855 LOGERR("Tried to insert a slice that is the wrong size.");
03856 return 1;
03857 }
03858 EMData* padfft = NULL;
03859
03860 if( padffted != 0 ) padfft = new EMData(*slice);
03861 else padfft = padfft_slice( slice, t, m_npad );
03862
03863 int mult= slice->get_attr_default("mult", 1);
03864
03865 Assert( mult > 0 );
03866 insert_padfft_slice( padfft, t, mult );
03867
03868 checked_delete( padfft );
03869 return 0;
03870 }
03871 int nnSSNR_ctfReconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
03872 {
03873 Assert( padfft != NULL );
03874
03875
03876 for (int isym=0; isym < m_nsym; isym++) {
03877 Transform tsym = t.get_sym(m_symmetry, isym);
03878 m_volume->nn_SSNR_ctf(m_wptr, m_wptr2, m_wptr3, padfft, tsym, mult);
03879 }
03880
03881 return 0;
03882 }
03883
03884 #define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
03885 EMData* nnSSNR_ctfReconstructor::finish(bool)
03886 {
03887
03888
03889
03890
03891
03892
03893
03894
03895
03896
03897
03898
03899
03900
03901 int kz, ky;
03902 int box = 7;
03903 int kc = (box-1)/2;
03904 float alpha = 0.0;
03905 float argx, argy, argz;
03906 vector< float > pow_a( 3*kc+1, 1.0 );
03907 float w = params["w"];
03908 float dx2 = 1.0f/float(m_vnxc)/float(m_vnxc);
03909 float dy2 = 1.0f/float(m_vnyc)/float(m_vnyc);
03910 #ifdef _WIN32
03911 float dz2 = 1.0f/_cpp_max(float(m_vnzc),1.0f)/_cpp_max(float(m_vnzc),1.0f);
03912 int inc = Util::round(float(_cpp_max(_cpp_max(m_vnxc,m_vnyc),m_vnzc))/w);
03913 #else
03914 float dz2 = 1.0f/std::max(float(m_vnzc),1.0f)/std::max(float(m_vnzc),1.0f);
03915 int inc = Util::round(float(std::max(std::max(m_vnxc,m_vnyc),m_vnzc))/w);
03916 #endif //_WIN32
03917
03918 EMData* SSNR = params["SSNR"];
03919 SSNR->set_size(inc+1,4,1);
03920
03921
03922
03923
03924 EMData* vol_ssnr = params["vol_ssnr"];
03925 vol_ssnr->set_size(m_vnxp+ 2 - m_vnxp%2, m_vnyp ,m_vnzp);
03926 vol_ssnr->to_zero();
03927 if ( m_vnxp % 2 == 0 ) vol_ssnr->set_fftodd(0);
03928 else vol_ssnr->set_fftodd(1);
03929 vol_ssnr->set_nxc(m_vnxc);
03930 vol_ssnr->set_complex(true);
03931 vol_ssnr->set_ri(true);
03932 vol_ssnr->set_fftpad(false);
03933
03934 float *nom = new float[inc+1];
03935 float *denom = new float[inc+1];
03936 int *ka = new int[inc+1];
03937 int *nn = new int[inc+1];
03938 float wght = 1.f;
03939 for (int i = 0; i <= inc; i++) {
03940 nom[i] = 0.0f;
03941 denom[i] = 0.0f;
03942 nn[i] = 0;
03943 ka[i] = 0;
03944 }
03945 m_volume->symplane2(m_wptr, m_wptr2, m_wptr3);
03946 if ( m_weighting == ESTIMATE ) {
03947 int vol = box*box*box;
03948 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03949 pow_a[3*kc] = 0.0;
03950 float max = max3d( kc, pow_a );
03951 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03952 }
03953 for (int iz = 1; iz <= m_vnzp; iz++) {
03954 if ( iz-1 > m_vnzc ) kz = iz-1-m_vnzp; else kz = iz-1;
03955 argz = float(kz*kz)*dz2;
03956 for (int iy = 1; iy <= m_vnyp; iy++) {
03957 if ( iy-1 > m_vnyc ) ky = iy-1-m_vnyp; else ky = iy-1;
03958 argy = argz + float(ky*ky)*dy2;
03959 for (int ix = 0; ix <= m_vnxc; ix++) {
03960 float Kn = (*m_wptr3)(ix,iy,iz);
03961 argx = std::sqrt(argy + float(ix*ix)*dx2);
03962 int r = Util::round(float(inc)*argx);
03963 if ( r >= 0 && Kn > 4.5f ) {
03964 if ( m_weighting == ESTIMATE ) {
03965 int cx = ix;
03966 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03967 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03968 float sum = 0.0;
03969 for( int ii = -kc; ii <= kc; ++ii ) {
03970 int nbrcx = cx + ii;
03971 if( nbrcx >= m_vnxc ) continue;
03972 for ( int jj= -kc; jj <= kc; ++jj ) {
03973 int nbrcy = cy + jj;
03974 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03975 for( int kk = -kc; kk <= kc; ++kk ) {
03976 int nbrcz = cz + jj;
03977 if ( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
03978 if( nbrcx < 0 ) {
03979 nbrcx = -nbrcx;
03980 nbrcy = -nbrcy;
03981 nbrcz = -nbrcz;
03982 }
03983 int nbrix = nbrcx;
03984 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
03985 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
03986 if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
03987 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
03988 sum = sum + pow_a[c];
03989 }
03990 }
03991 }
03992 }
03993
03994 wght = 1.0f / ( 1.0f - alpha * sum );
03995 }
03996 float nominator = std::norm(m_volume->cmplx(ix,iy,iz))/(*m_wptr)(ix,iy,iz);
03997 float denominator = ((*m_wptr2)(ix,iy,iz)-std::norm(m_volume->cmplx(ix,iy,iz))/(*m_wptr)(ix,iy,iz))/(Kn-1.0f);
03998
03999 if( (ix>0 || (kz>=0 && (ky>=0 || kz!=0)))) {
04000 if ( r <= inc ) {
04001 nom[r] += nominator*wght;
04002 denom[r] += denominator/(*m_wptr)(ix,iy,iz)*wght;
04003 nn[r] += 2;
04004 ka[r] += int(Kn);
04005 }
04006
04007
04008
04009
04010
04011
04012
04013
04014
04015
04016
04017
04018
04019
04020
04021 }
04022 (*vol_ssnr)(2*ix, iy-1, iz-1) = denominator*wght;
04023 }
04024 }
04025 }
04026 }
04027 for (int i = 0; i <= inc; i++) {
04028 (*SSNR)(i,0,0) = nom[i];
04029 (*SSNR)(i,1,0) = denom[i];
04030 (*SSNR)(i,2,0) = static_cast<float>(nn[i]);
04031 (*SSNR)(i,3,0) = static_cast<float>(ka[i]);
04032 }
04033 vol_ssnr->update();
04034
04035 delete[] nom;
04036 delete[] denom;
04037 delete[] nn;
04038 delete[] ka;
04039
04040 return 0;
04041 }
04042 #undef tw
04043
04044
04045
04046
04047 void EMAN::dump_reconstructors()
04048 {
04049 dump_factory < Reconstructor > ();
04050 }
04051
04052 map<string, vector<string> > EMAN::dump_reconstructors_list()
04053 {
04054 return dump_factory_list < Reconstructor > ();
04055 }
04056
04057
04058
04059
04060
04061 using std::ofstream;
04062 using std::ifstream;
04063
04064
04065 newfile_store::newfile_store( const string& filename, int npad, bool ctf )
04066 : m_bin_file( filename + ".bin" ),
04067 m_txt_file( filename + ".txt" )
04068 {
04069 m_npad = npad;
04070 m_ctf = ctf;
04071 }
04072
04073 newfile_store::~newfile_store( )
04074 {
04075 }
04076
04077 void newfile_store::add_image( EMData* emdata, const Transform& tf )
04078 {
04079 if( m_bin_of == NULL ) {
04080 m_bin_of = shared_ptr<ofstream>( new ofstream(m_bin_file.c_str(), std::ios::out|std::ios::binary) );
04081 m_txt_of = shared_ptr<ofstream>( new ofstream(m_txt_file.c_str()) );
04082 }
04083
04084
04085 EMData* padfft = padfft_slice( emdata, tf, m_npad );
04086
04087 int nx = padfft->get_xsize();
04088 int ny = padfft->get_ysize();
04089 int n2 = ny / 2;
04090 int n = ny;
04091
04092 float voltage=0.0f, pixel=0.0f, Cs=0.0f, ampcont=0.0f, bfactor=0.0f, defocus=0.0f;
04093
04094 if( m_ctf ) {
04095 Ctf* ctf = emdata->get_attr( "ctf" );
04096 Dict params = ctf->to_dict();
04097 voltage = params["voltage"];
04098 pixel = params["apix"];
04099 Cs = params["cs"];
04100 ampcont = params["ampcont"];
04101 bfactor = params["bfactor"];
04102 defocus = params["defocus"];
04103 if(ctf) {delete ctf; ctf=0;}
04104 }
04105
04106 vector<point_t> points;
04107 for( int j=-ny/2+1; j <= ny/2; j++ ) {
04108 int jp = (j>=0) ? j+1 : ny+j+1;
04109 for( int i=0; i <= n2; ++i ) {
04110 int r2 = i*i + j*j;
04111 if( (r2<ny*ny/4) && !( (i==0) && (j<0) ) ) {
04112 float ctf;
04113 if( m_ctf ) {
04114 float ak = std::sqrt( r2/float(ny*ny) )/pixel;
04115 ctf = Util::tf( defocus, ak, voltage, Cs, ampcont, bfactor, 1);
04116 } else {
04117 ctf = 1.0;
04118 }
04119
04120 float xnew = i*tf[0][0] + j*tf[1][0];
04121 float ynew = i*tf[0][1] + j*tf[1][1];
04122 float znew = i*tf[0][2] + j*tf[1][2];
04123 std::complex<float> btq;
04124 if (xnew < 0.) {
04125 xnew = -xnew;
04126 ynew = -ynew;
04127 znew = -znew;
04128 btq = conj(padfft->cmplx(i,jp-1));
04129 } else {
04130 btq = padfft->cmplx(i,jp-1);
04131 }
04132
04133 int ixn = int(xnew + 0.5 + n) - n;
04134 int iyn = int(ynew + 0.5 + n) - n;
04135 int izn = int(znew + 0.5 + n) - n;
04136 if ((ixn <= n2) && (iyn >= -n2) && (iyn <= n2) && (izn >= -n2) && (izn <= n2)) {
04137 int ixf, iyf, izf;
04138 if (ixn >= 0) {
04139 int iza, iya;
04140 if (izn >= 0)
04141 iza = izn + 1;
04142 else
04143 iza = n + izn + 1;
04144
04145 if (iyn >= 0)
04146 iya = iyn + 1;
04147 else
04148 iya = n + iyn + 1;
04149
04150 ixf = ixn;
04151 iyf = iya;
04152 izf = iza;
04153 } else {
04154 int izt, iyt;
04155 if (izn > 0)
04156 izt = n - izn + 1;
04157 else
04158 izt = -izn + 1;
04159
04160 if (iyn > 0)
04161 iyt = n - iyn + 1;
04162 else
04163 iyt = -iyn + 1;
04164
04165 ixf = -ixn;
04166 iyf = iyt;
04167 izf = izt;
04168 }
04169
04170
04171 int pos2 = ixf + (iyf-1)*nx/2 + (izf-1)*ny*nx/2;
04172 float ctfv1 = btq.real() * ctf;
04173 float ctfv2 = btq.imag() * ctf;
04174 float ctf2 = ctf*ctf;
04175
04176 point_t p;
04177 p.pos2 = pos2;
04178 p.real = ctfv1;
04179 p.imag = ctfv2;
04180 p.ctf2 = ctf2;
04181
04182 points.push_back( p );
04183 }
04184 }
04185 }
04186 }
04187
04188
04189 int npoint = points.size();
04190 std::istream::off_type offset = (m_offsets.size()==0) ? 0 : m_offsets.back();
04191 offset += npoint*sizeof(point_t);
04192 m_offsets.push_back( offset );
04193
04194 *m_txt_of << m_offsets.back() << std::endl;
04195 m_bin_of->write( (char*)(&points[0]), sizeof(point_t)*npoint );
04196 checked_delete( padfft );
04197 }
04198
04199 void newfile_store::get_image( int id, EMData* buf )
04200 {
04201 if( m_offsets.size()==0 ) {
04202 ifstream is( m_txt_file.c_str() );
04203 std::istream::off_type off;
04204 while( is >> off ) {
04205 m_offsets.push_back( off );
04206 }
04207
04208 m_bin_if = shared_ptr<std::ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
04209 }
04210
04211 Assert( m_bin_if != NULL );
04212
04213 std::istream::off_type offset = (id==0) ? 0 : m_offsets[id-1];
04214 Assert( offset >= 0 );
04215 m_bin_if->seekg( offset, std::ios::beg );
04216
04217
04218 if( m_bin_if->bad() || m_bin_if->fail() || m_bin_if->eof() ) {
04219 std::cout << "bad or fail or eof while fetching id, offset: " << id << " " << offset << std::endl;
04220 throw std::logic_error( "bad happen" );
04221 }
04222
04223 int bufsize = (m_offsets[id] - offset)/sizeof(float);
04224 if( buf->get_xsize() != bufsize ) {
04225 buf->set_size( bufsize, 1, 1 );
04226 }
04227
04228 char* data = (char*)(buf->get_data());
04229 m_bin_if->read( data, sizeof(float)*bufsize );
04230 buf->update();
04231 }
04232
04233 void newfile_store::read( int nprj )
04234 {
04235 if( m_offsets.size()==0 ) {
04236 ifstream is( m_txt_file.c_str() );
04237 std::istream::off_type off;
04238 while( is >> off ) {
04239 m_offsets.push_back( off );
04240 }
04241 }
04242
04243 if( m_bin_if==NULL ) {
04244 m_bin_if = shared_ptr< ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
04245 }
04246
04247
04248 int npoint = m_offsets[0]/sizeof(point_t);
04249 std::ios::off_type prjsize = m_offsets[0];
04250
04251 try {
04252 m_points.resize(nprj * npoint);
04253 }
04254 catch( std::exception& e ) {
04255 std::cout << "Error: " << e.what() << std::endl;
04256 }
04257
04258 int ip = 0;
04259 for( int i=0; i < nprj; ++i ) {
04260 m_bin_if->read( (char*)(&m_points[ip]), prjsize );
04261 if( m_bin_if->bad() || m_bin_if->fail() || m_bin_if->eof() ) {
04262 std::cout << "Error: file hander bad or fail or eof" << std::endl;
04263 return;
04264 }
04265 ip += npoint;
04266 }
04267 }
04268
04269 void newfile_store::add_tovol( EMData* fftvol, EMData* wgtvol, const vector<int>& mults, int pbegin, int pend )
04270 {
04271 float* vdata = fftvol->get_data();
04272 float* wdata = wgtvol->get_data();
04273
04274 int npoint = m_offsets[0]/sizeof(point_t);
04275
04276 Assert( int(m_points.size())== (pend - pbegin)*npoint );
04277
04278 for( int iprj=pbegin; iprj < pend; ++iprj ) {
04279 int m = mults[iprj];
04280 if( m==0 ) continue;
04281
04282 int ipt = (iprj-pbegin)*npoint;
04283 for( int i=0; i < npoint; ++i ) {
04284 int pos2 = m_points[ipt].pos2;
04285 int pos1 = pos2*2;
04286
04287 wdata[pos2] += m_points[ipt].ctf2*m;
04288 vdata[pos1] += m_points[ipt].real*m;
04289 vdata[pos1+1] += m_points[ipt].imag*m;
04290 ++ipt;
04291 }
04292 }
04293 }
04294
04295 void newfile_store::restart()
04296 {
04297 m_bin_if = shared_ptr< ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
04298 }
04299
04300 file_store::file_store(const string& filename, int npad, int write, bool ctf)
04301 : m_bin_file(filename + ".bin"),
04302 m_txt_file(filename + ".txt")
04303 {
04304 m_ctf = ctf;
04305 m_prev = -1;
04306 m_npad = npad;
04307 m_write = write;
04308 }
04309
04310 file_store::~file_store()
04311 {
04312 }
04313
04314 void file_store::add_image( EMData* emdata, const Transform& tf )
04315 {
04316
04317 EMData* padfft = padfft_slice( emdata, tf, m_npad );
04318
04319 float* data = padfft->get_data();
04320
04321 if( m_write && m_bin_ohandle == NULL )
04322 {
04323 m_bin_ohandle = shared_ptr< ofstream >( new ofstream(m_bin_file.c_str(), std::ios::out | std::ios::binary) );
04324 m_txt_ohandle = shared_ptr< ofstream >( new ofstream(m_txt_file.c_str() ) );
04325 if( m_ctf )
04326 *m_txt_ohandle << "Cs pixel voltage ctf_applied amp_contrast defocus ";
04327
04328 *m_txt_ohandle << "phi theta psi" << std::endl;
04329 }
04330
04331 m_x_out = padfft->get_xsize();
04332 m_y_out = padfft->get_ysize();
04333 m_z_out = padfft->get_zsize();
04334 m_totsize = m_x_out*m_y_out*m_z_out;
04335
04336 if( m_ctf )
04337 {
04338 Ctf* ctf = padfft->get_attr( "ctf" );
04339 Dict ctf_params = ctf->to_dict();
04340
04341 m_ctf_applied = padfft->get_attr( "ctf_applied" );
04342
04343 m_Cs = ctf_params["cs"];
04344 m_pixel = ctf_params["apix"];
04345 m_voltage = ctf_params["voltage"];
04346 m_amp_contrast = ctf_params["ampcont"];
04347 m_defocuses.push_back( ctf_params["defocus"] );
04348 if(ctf) {delete ctf; ctf=0;}
04349 }
04350
04351 Dict params = tf.get_rotation( "spider" );
04352 float phi = params.get( "phi" );
04353 float tht = params.get( "theta" );
04354 float psi = params.get( "psi" );
04355
04356
04357 m_phis.push_back( phi );
04358 m_thetas.push_back( tht );
04359 m_psis.push_back( psi );
04360
04361 if( m_write )
04362 {
04363 m_bin_ohandle->write( (char*)data, sizeof(float)*m_totsize );
04364
04365 if( m_ctf )
04366 {
04367 *m_txt_ohandle << m_Cs << " ";
04368 *m_txt_ohandle << m_pixel << " ";
04369 *m_txt_ohandle << m_voltage << " ";
04370 *m_txt_ohandle << m_ctf_applied << " ";
04371 *m_txt_ohandle << m_amp_contrast << " ";
04372 *m_txt_ohandle << m_defocuses.back() << " ";
04373 }
04374 *m_txt_ohandle << m_phis.back() << " ";
04375 *m_txt_ohandle << m_thetas.back() << " ";
04376 *m_txt_ohandle << m_psis.back() << " ";
04377 *m_txt_ohandle << m_x_out << " ";
04378 *m_txt_ohandle << m_y_out << " ";
04379 *m_txt_ohandle << m_z_out << " ";
04380 *m_txt_ohandle << m_totsize << std::endl;
04381 }
04382
04383 checked_delete(padfft);
04384
04385 }
04386
04387 void file_store::get_image( int id, EMData* padfft )
04388 {
04389
04390 if( m_phis.size() == 0 ) {
04391 ifstream m_txt_ifs( m_txt_file.c_str() );
04392
04393 if( !m_txt_ifs )
04394 {
04395 std::cerr << "Error: file " << m_txt_file << " does not exist" << std::endl;
04396 }
04397
04398 string line;
04399 std::getline( m_txt_ifs, line );
04400
04401 float first, defocus, phi, theta, psi;
04402
04403
04404
04405 while( m_txt_ifs >> first ) {
04406
04407 if( m_ctf )
04408 {
04409 m_Cs = first;
04410 m_txt_ifs >> m_pixel >> m_voltage;
04411 m_txt_ifs >> m_ctf_applied >> m_amp_contrast;
04412 m_txt_ifs >> defocus >> phi >> theta >> psi;
04413 m_defocuses.push_back( defocus );
04414 }
04415 else
04416 {
04417 phi = first;
04418 m_txt_ifs >> theta >> psi;
04419 }
04420
04421 m_txt_ifs >> m_x_out >> m_y_out >> m_z_out >> m_totsize;
04422 m_phis.push_back( phi );
04423 m_thetas.push_back( theta );
04424 m_psis.push_back( psi );
04425 }
04426 }
04427
04428 Assert( m_ihandle != NULL );
04429
04430 std::istream::off_type offset = id*sizeof(float)*m_totsize;
04431 Assert( offset >= 0 );
04432
04433 if( offset > 0 )
04434 {
04435 m_ihandle->seekg(offset, std::ios::beg);
04436 }
04437
04438 if( m_ihandle->bad() )
04439 {
04440 std::cout << "bad while fetching id, offset: " << id << " " << offset << std::endl;
04441 throw std::logic_error( "bad happen" );
04442 }
04443
04444 if( m_ihandle->fail() )
04445 {
04446 std::cout << "fail while fetching id, offset, curoff: " << id << " " << offset << std::endl;
04447 throw std::logic_error( "fail happen" );
04448 }
04449
04450 if( m_ihandle->eof() )
04451 {
04452 std::cout << "eof while fetching id, offset: " << id << " " << offset << std::endl;
04453 throw std::logic_error( "eof happen" );
04454 }
04455
04456 if( padfft->get_xsize() != m_x_out ||
04457 padfft->get_ysize() != m_y_out ||
04458 padfft->get_zsize() != m_z_out )
04459 {
04460 padfft->set_size(m_x_out, m_y_out, m_z_out);
04461 }
04462
04463 char* data = (char*)(padfft->get_data());
04464 m_ihandle->read( data, sizeof(float)*m_totsize );
04465 padfft->update();
04466
04467 if( m_ctf )
04468 {
04469 padfft->set_attr( "Cs", m_Cs );
04470 padfft->set_attr( "Pixel_size", m_pixel );
04471 padfft->set_attr( "voltage", m_voltage );
04472 padfft->set_attr( "ctf_applied", m_ctf_applied );
04473 padfft->set_attr( "amp_contrast", m_amp_contrast );
04474 padfft->set_attr( "defocus", m_defocuses[id] );
04475 }
04476
04477 padfft->set_attr( "padffted", 1 );
04478 padfft->set_attr( "phi", m_phis[id] );
04479 padfft->set_attr( "theta", m_thetas[id] );
04480 padfft->set_attr( "psi", m_psis[id] );
04481
04482 }
04483
04484 void file_store::restart( )
04485 {
04486 if( m_ihandle == NULL )
04487 {
04488 m_ihandle = shared_ptr< ifstream >( new ifstream(m_bin_file.c_str(), std::ios::in | std::ios::binary) );
04489 }
04490
04491 if( m_ihandle->bad() || m_ihandle->fail() || m_ihandle->eof() )
04492 {
04493 m_ihandle->open( m_bin_file.c_str(), std::ios::binary );
04494 }
04495
04496 m_ihandle->seekg( 0, std::ios::beg );
04497 }
04498
04499