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