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