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) = denominator*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
03333 else m_volume->nn_ctf(m_wptr, padfft, tsym[isym], mult);
03334
03335 }
03336
03337
03338 return 0;
03339
03340 }
03341
03342 #define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
03343 EMData* nn4_ctfReconstructor::finish(bool)
03344 {
03345 m_volume->set_array_offsets(0, 1, 1);
03346 m_wptr->set_array_offsets(0, 1, 1);
03347 m_volume->symplane0_ctf(m_wptr);
03348
03349 int box = 7;
03350 int vol = box*box*box;
03351 int kc = (box-1)/2;
03352 vector< float > pow_a( 3*kc+1, 1.0 );
03353 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03354 pow_a[3*kc]=0.0;
03355
03356
03357 float max = max3d( kc, pow_a );
03358 float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03359 float osnr = 1.0f/m_snr;
03360
03361
03362 int ix,iy,iz;
03363 for (iz = 1; iz <= m_vnzp; iz++) {
03364 for (iy = 1; iy <= m_vnyp; iy++) {
03365 for (ix = 0; ix <= m_vnxc; ix++) {
03366 if ( (*m_wptr)(ix,iy,iz) > 0.0f) {
03367 int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
03368 int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
03369 float tmp=0.0;
03370 if( m_varsnr ) {
03371 float freq = sqrt( (float)(ix*ix+iyp*iyp+izp*izp) );
03372 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+freq*osnr)*m_sign;
03373 } else {
03374 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+osnr)*m_sign;
03375 }
03376
03377 if( m_weighting == ESTIMATE ) {
03378 int cx = ix;
03379 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03380 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03381 float sum = 0.0;
03382 for( int ii = -kc; ii <= kc; ++ii ) {
03383 int nbrcx = cx + ii;
03384 if( nbrcx >= m_vnxc ) continue;
03385 for( int jj= -kc; jj <= kc; ++jj ) {
03386 int nbrcy = cy + jj;
03387 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03388 for( int kk = -kc; kk <= kc; ++kk ) {
03389 int nbrcz = cz + jj;
03390 if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
03391 if( nbrcx < 0 ) {
03392 nbrcx = -nbrcx;
03393 nbrcy = -nbrcy;
03394 nbrcz = -nbrcz;
03395 }
03396
03397 int nbrix = nbrcx;
03398 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
03399 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
03400 if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0.0 ) {
03401 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
03402 sum = sum + pow_a[c];
03403
03404
03405 }
03406 }
03407 }
03408 }
03409 float wght = 1.0f / ( 1.0f - alpha * sum );
03410
03411
03412
03413
03414
03415
03416
03417
03418
03419 tmp = tmp * wght;
03420 }
03421 (*m_volume)(2*ix,iy,iz) *= tmp;
03422 (*m_volume)(2*ix+1,iy,iz) *= tmp;
03423 }
03424 }
03425 }
03426 }
03427
03428
03429 m_volume->do_ift_inplace();
03430 int npad = m_volume->get_attr("npad");
03431 m_volume->depad();
03432 circumf( m_volume, npad );
03433 m_volume->set_array_offsets( 0, 0, 0 );
03434
03435 return 0;
03436 }
03437 #undef tw
03438
03439
03440
03441
03442
03443
03444 nn4_ctf_rectReconstructor::nn4_ctf_rectReconstructor()
03445 {
03446 m_volume = NULL;
03447 m_wptr = NULL;
03448 }
03449
03450 nn4_ctf_rectReconstructor::nn4_ctf_rectReconstructor( const string& symmetry, int size, int npad, float snr, int sign )
03451 {
03452 setup( symmetry, size, npad, snr, sign );
03453 }
03454
03455 nn4_ctf_rectReconstructor::~nn4_ctf_rectReconstructor()
03456 {
03457
03458
03459
03460
03461
03462 }
03463
03464 void nn4_ctf_rectReconstructor::setup()
03465 {
03466 if( ! params.has_key("sizeprojection") ) throw std::logic_error("Error: projection size is not given");
03467 m_sizeofprojection = params["sizeprojection"];
03468 int npad = params.has_key("npad") ? int(params["npad"]) : 4;
03469
03470 int sign = 1;
03471 string symmetry = params.has_key("symmetry")? params["symmetry"].to_str() : "c1";
03472
03473 float snr = params["snr"];
03474
03475 m_varsnr = params.has_key("varsnr") ? int(params["varsnr"]) : 0;
03476 setup( symmetry, m_sizeofprojection, npad, snr, sign );
03477
03478 }
03479
03480 void nn4_ctf_rectReconstructor::setup( const string& symmetry, int sizeprojection, int npad, float snr, int sign )
03481 {
03482 m_weighting = ESTIMATE;
03483 if( params.has_key("weighting") ) {
03484 int tmp = int( params["weighting"] );
03485 if( tmp==0 ) m_weighting = NONE;
03486 }
03487
03488 m_wghta = 0.2f;
03489 m_wghtb = 0.004f;
03490
03491 m_symmetry = symmetry;
03492 m_npad = npad;
03493 m_sign = sign;
03494 m_nsym = Transform::get_nsym(m_symmetry);
03495
03496 m_snr = snr;
03497 if (params.has_key("sizex")) m_vnx = params["sizex"];
03498 else if (params.has_key("xratio")) {
03499 float temp=params["xratio"];
03500 m_vnx=int(float(sizeprojection)*temp);
03501 } else m_vnx=sizeprojection;
03502
03503 if (params.has_key("sizey")) m_vny = params["sizey"];
03504 else if (params.has_key("yratio")) {
03505 float temp=params["yratio"];
03506 m_vny=int(float(sizeprojection)*temp);
03507 }
03508 else m_vny=sizeprojection;
03509
03510 if (params.has_key("sizez")) m_vnz = params["sizez"];
03511 else if (params.has_key("zratio")) {
03512 float temp=params["zratio"];
03513 m_vnz=int(float(sizeprojection)*temp);
03514 }
03515 else m_vnz=sizeprojection;
03516
03517
03518 m_xratio=float(m_vnx)/float(sizeprojection);
03519 m_yratio=float(m_vny)/float(sizeprojection);
03520 m_zratio=float(m_vnz)/float(sizeprojection);
03521
03522
03523
03524
03525 m_vnxp = m_vnx*npad;
03526 m_vnyp = m_vny*npad;
03527 m_vnzp = m_vnz*npad;
03528
03529 m_vnxc = m_vnxp/2;
03530 m_vnyc = m_vnyp/2;
03531 m_vnzc = m_vnzp/2;
03532
03533 buildFFTVolume();
03534 buildNormVolume();
03535 }
03536
03537 void nn4_ctf_rectReconstructor::buildFFTVolume() {
03538 int offset = 2 - m_vnxp%2;
03539
03540 m_volume = params["fftvol"];
03541
03542 if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
03543 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
03544 m_volume->to_zero();
03545 }
03546
03547 m_volume->set_nxc(m_vnxp/2);
03548 m_volume->set_complex(true);
03549 m_volume->set_ri(true);
03550 m_volume->set_fftpad(true);
03551 m_volume->set_attr("npad", m_npad);
03552 m_volume->set_array_offsets(0,1,1);
03553 }
03554
03555 void nn4_ctf_rectReconstructor::buildNormVolume()
03556 {
03557 m_wptr = params["weight"];
03558
03559 if( m_wptr->get_xsize() != m_vnxc+1 && m_wptr->get_ysize() != m_vnyp && m_wptr->get_zsize() != m_vnzp ) {
03560 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03561 m_wptr->to_zero();
03562 }
03563
03564 m_wptr->set_array_offsets(0,1,1);
03565
03566 }
03567
03568 int nn4_ctf_rectReconstructor::insert_slice(const EMData* const slice, const Transform& t, const float)
03569 {
03570
03571 if (!slice) {
03572 LOGERR("try to insert NULL slice");
03573 return 1;
03574 }
03575
03576 int buffed = slice->get_attr_default( "buffed", 0 );
03577 if( buffed > 0 ) {
03578 int mult = slice->get_attr_default( "mult", 1 );
03579 insert_buffed_slice( slice, mult );
03580 return 0;
03581 }
03582
03583 int padffted= slice->get_attr_default("padffted", 0);
03584
03585 if( padffted==0 && (slice->get_xsize()!=slice->get_ysize()) )
03586 {
03587
03588 LOGERR("Tried to insert a slice that is the wrong size.");
03589 return 1;
03590 }
03591
03592 EMData* padfft = NULL;
03593
03594 if( padffted != 0 ) padfft = new EMData(*slice);
03595 else padfft = padfft_slice( slice, t, m_npad );
03596
03597 int mult= slice->get_attr_default("mult", 1);
03598
03599 Assert( mult > 0 );
03600 insert_padfft_slice( padfft, t, mult );
03601
03602 checked_delete( padfft );
03603
03604 return 0;
03605 }
03606
03607 int nn4_ctf_rectReconstructor::insert_buffed_slice( const EMData* buffed, int mult )
03608 {
03609 const float* bufdata = buffed->get_data();
03610 float* cdata = m_volume->get_data();
03611 float* wdata = m_wptr->get_data();
03612
03613 int npoint = buffed->get_xsize()/4;
03614 for( int i=0; i < npoint; ++i ) {
03615
03616 int pos2 = int( bufdata[4*i] );
03617 int pos1 = pos2 * 2;
03618 cdata[pos1 ] += bufdata[4*i+1]*mult;
03619 cdata[pos1+1] += bufdata[4*i+2]*mult;
03620 wdata[pos2 ] += bufdata[4*i+3]*mult;
03621
03622
03623
03624
03625
03626 }
03627 return 0;
03628 }
03629
03630
03631 int nn4_ctf_rectReconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
03632 {
03633 Assert( padfft != NULL );
03634 float tmp = padfft->get_attr("ctf_applied");
03635 int ctf_applied = (int) tmp;
03636 vector<Transform> tsym = t.get_sym_proj(m_symmetry);
03637 for (unsigned int isym=0; isym < tsym.size(); isym++) {
03638 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);
03639
03640 else m_volume->insert_rect_slice_ctf(m_wptr, padfft, tsym[isym], m_sizeofprojection, m_xratio, m_yratio, m_zratio, m_npad, mult);
03641
03642 }
03643
03644
03645 return 0;
03646
03647 }
03648
03649 #define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
03650 EMData* nn4_ctf_rectReconstructor::finish(bool)
03651 {
03652 m_volume->set_array_offsets(0, 1, 1);
03653 m_wptr->set_array_offsets(0, 1, 1);
03654 m_volume->symplane0_rect(m_wptr);
03655
03656 int box = 7;
03657 int vol = box*box*box;
03658 int kc = (box-1)/2;
03659 vector< float > pow_a( 3*kc+1, 1.0 );
03660 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03661 pow_a[3*kc]=0.0;
03662
03663
03664 float max = max3d( kc, pow_a );
03665 float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03666 float osnr = 1.0f/m_snr;
03667
03668
03669 int ix,iy,iz;
03670 for (iz = 1; iz <= m_vnzp; iz++) {
03671 for (iy = 1; iy <= m_vnyp; iy++) {
03672 for (ix = 0; ix <= m_vnxc; ix++) {
03673 if ( (*m_wptr)(ix,iy,iz) > 0.0f) {
03674 int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
03675 int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
03676 float tmp=0.0;
03677 if( m_varsnr )
03678 {
03679 float freq = sqrt( (float)(ix*ix/(m_xratio*m_xratio)+iyp*iyp/(m_zratio*m_yratio)+izp*izp) );
03680 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+freq*osnr)*m_sign;
03681 }
03682 else
03683 {
03684 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+osnr)*m_sign;
03685 }
03686
03687 if( m_weighting == ESTIMATE ) {
03688 int cx = ix;
03689 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03690 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03691 float sum = 0.0;
03692 for( int ii = -kc; ii <= kc; ++ii ) {
03693 int nbrcx = cx + ii;
03694 if( nbrcx >= m_vnxc ) continue;
03695 for( int jj= -kc; jj <= kc; ++jj ) {
03696 int nbrcy = cy + jj;
03697 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03698 for( int kk = -kc; kk <= kc; ++kk ) {
03699 int nbrcz = cz + jj;
03700 if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
03701 if( nbrcx < 0 ) {
03702 nbrcx = -nbrcx;
03703 nbrcy = -nbrcy;
03704 nbrcz = -nbrcz;
03705 }
03706
03707 int nbrix = nbrcx;
03708 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
03709 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
03710 if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0.0 ) {
03711 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
03712 sum = sum + pow_a[c];
03713
03714
03715 }
03716 }
03717 }
03718 }
03719 float wght = 1.0f / ( 1.0f - alpha * sum );
03720
03721
03722
03723
03724
03725
03726
03727
03728
03729 tmp = tmp * wght;
03730 }
03731 (*m_volume)(2*ix,iy,iz) *= tmp;
03732 (*m_volume)(2*ix+1,iy,iz) *= tmp;
03733 }
03734 }
03735 }
03736 }
03737
03738
03739 m_volume->do_ift_inplace();
03740 int npad = m_volume->get_attr("npad");
03741 m_volume->depad();
03742 circumf_rect( m_volume, npad );
03743 m_volume->set_array_offsets( 0, 0, 0 );
03744 return 0;
03745 }
03746 #undef tw
03747
03748
03749
03750
03751
03752
03753
03754 nnSSNR_ctfReconstructor::nnSSNR_ctfReconstructor()
03755 {
03756 m_volume = NULL;
03757 m_wptr = NULL;
03758 m_wptr2 = NULL;
03759 m_wptr3 = NULL;
03760 }
03761
03762 nnSSNR_ctfReconstructor::nnSSNR_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign)
03763 {
03764 m_volume = NULL;
03765 m_wptr = NULL;
03766 m_wptr2 = NULL;
03767 m_wptr3 = NULL;
03768
03769 setup( symmetry, size, npad, snr, sign );
03770 }
03771
03772 nnSSNR_ctfReconstructor::~nnSSNR_ctfReconstructor()
03773 {
03774
03775
03776
03777
03778
03779
03780 }
03781
03782 void nnSSNR_ctfReconstructor::setup()
03783 {
03784 int size = params["size"];
03785 int npad = params["npad"];
03786 int sign = params["sign"];
03787 float snr = params["snr"];
03788 string symmetry;
03789 if( params.has_key("symmetry") ) symmetry = params["symmetry"].to_str();
03790 else symmetry = "c1";
03791
03792 setup( symmetry, size, npad, snr, sign );
03793 }
03794 void nnSSNR_ctfReconstructor::setup( const string& symmetry, int size, int npad, float snr, int sign )
03795 {
03796
03797 m_weighting = ESTIMATE;
03798 m_wghta = 0.2f;
03799 m_wghtb = 0.004f;
03800 wiener = 1;
03801
03802 m_symmetry = symmetry;
03803 m_npad = npad;
03804 m_nsym = Transform::get_nsym(m_symmetry);
03805
03806 m_sign = sign;
03807 m_snr = snr;
03808
03809 m_vnx = size;
03810 m_vny = size;
03811 m_vnz = size;
03812
03813 m_vnxp = size*npad;
03814 m_vnyp = size*npad;
03815 m_vnzp = size*npad;
03816
03817 m_vnxc = m_vnxp/2;
03818 m_vnyc = m_vnyp/2;
03819 m_vnzc = m_vnzp/2;
03820
03821 buildFFTVolume();
03822 buildNormVolume();
03823 buildNorm2Volume();
03824 buildNorm3Volume();
03825 }
03826
03827 void nnSSNR_ctfReconstructor::buildFFTVolume() {
03828
03829 int offset = 2 - m_vnxp%2;
03830 m_volume = params["fftvol"];
03831
03832 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
03833 m_volume->to_zero();
03834
03835 m_volume->set_fftodd(m_vnxp % 2);
03836
03837 m_volume->set_nxc(m_vnxc);
03838 m_volume->set_complex(true);
03839 m_volume->set_ri(true);
03840 m_volume->set_fftpad(true);
03841 m_volume->set_attr("npad", m_npad);
03842 m_volume->set_array_offsets(0,1,1);
03843 }
03844
03845
03846 void nnSSNR_ctfReconstructor::buildNormVolume()
03847 {
03848 m_wptr = params["weight"];
03849 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03850 m_wptr->to_zero();
03851 m_wptr->set_array_offsets(0,1,1);
03852 }
03853
03854 void nnSSNR_ctfReconstructor::buildNorm2Volume() {
03855
03856 m_wptr2 = params["weight2"];
03857 m_wptr2->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03858 m_wptr2->to_zero();
03859 m_wptr2->set_array_offsets(0,1,1);
03860 }
03861
03862 void nnSSNR_ctfReconstructor::buildNorm3Volume() {
03863
03864 m_wptr3 = params["weight3"];
03865 m_wptr3->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03866 m_wptr3->to_zero();
03867 m_wptr3->set_array_offsets(0,1,1);
03868 }
03869
03870 int nnSSNR_ctfReconstructor::insert_slice(const EMData *const slice, const Transform& t, const float) {
03871
03872 if (!slice) {
03873 LOGERR("try to insert NULL slice");
03874 return 1;
03875 }
03876 int padffted= slice->get_attr_default("padffted", 0);
03877 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx) ) {
03878
03879 LOGERR("Tried to insert a slice that is the wrong size.");
03880 return 1;
03881 }
03882 EMData* padfft = NULL;
03883
03884 if( padffted != 0 ) padfft = new EMData(*slice);
03885 else padfft = padfft_slice( slice, t, m_npad );
03886
03887 int mult= slice->get_attr_default("mult", 1);
03888
03889 Assert( mult > 0 );
03890 insert_padfft_slice( padfft, t, mult );
03891
03892 checked_delete( padfft );
03893 return 0;
03894 }
03895 int nnSSNR_ctfReconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
03896 {
03897 Assert( padfft != NULL );
03898
03899
03900 for (int isym=0; isym < m_nsym; isym++) {
03901 Transform tsym = t.get_sym(m_symmetry, isym);
03902 m_volume->nn_SSNR_ctf(m_wptr, m_wptr2, m_wptr3, padfft, tsym, mult);
03903 }
03904
03905 return 0;
03906 }
03907
03908 #define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
03909 EMData* nnSSNR_ctfReconstructor::finish(bool)
03910 {
03911
03912
03913
03914
03915
03916
03917
03918
03919
03920
03921
03922
03923
03924
03925 int kz, ky;
03926 int box = 7;
03927 int kc = (box-1)/2;
03928 float alpha = 0.0;
03929 float argx, argy, argz;
03930 vector< float > pow_a( 3*kc+1, 1.0 );
03931 float w = params["w"];
03932 float dx2 = 1.0f/float(m_vnxc)/float(m_vnxc);
03933 float dy2 = 1.0f/float(m_vnyc)/float(m_vnyc);
03934 #ifdef _WIN32
03935 float dz2 = 1.0f/_cpp_max(float(m_vnzc),1.0f)/_cpp_max(float(m_vnzc),1.0f);
03936 int inc = Util::round(float(_cpp_max(_cpp_max(m_vnxc,m_vnyc),m_vnzc))/w);
03937 #else
03938 float dz2 = 1.0f/std::max(float(m_vnzc),1.0f)/std::max(float(m_vnzc),1.0f);
03939 int inc = Util::round(float(std::max(std::max(m_vnxc,m_vnyc),m_vnzc))/w);
03940 #endif //_WIN32
03941
03942 EMData* SSNR = params["SSNR"];
03943 SSNR->set_size(inc+1,4,1);
03944
03945
03946
03947
03948 EMData* vol_ssnr = params["vol_ssnr"];
03949 vol_ssnr->set_size(m_vnxp+ 2 - m_vnxp%2, m_vnyp ,m_vnzp);
03950 vol_ssnr->to_zero();
03951 if ( m_vnxp % 2 == 0 ) vol_ssnr->set_fftodd(0);
03952 else vol_ssnr->set_fftodd(1);
03953 vol_ssnr->set_nxc(m_vnxc);
03954 vol_ssnr->set_complex(true);
03955 vol_ssnr->set_ri(true);
03956 vol_ssnr->set_fftpad(false);
03957
03958 float *nom = new float[inc+1];
03959 float *denom = new float[inc+1];
03960 int *ka = new int[inc+1];
03961 int *nn = new int[inc+1];
03962 float wght = 1.f;
03963 for (int i = 0; i <= inc; i++) {
03964 nom[i] = 0.0f;
03965 denom[i] = 0.0f;
03966 nn[i] = 0;
03967 ka[i] = 0;
03968 }
03969 m_volume->symplane2(m_wptr, m_wptr2, m_wptr3);
03970 if ( m_weighting == ESTIMATE ) {
03971 int vol = box*box*box;
03972 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03973 pow_a[3*kc] = 0.0;
03974 float max = max3d( kc, pow_a );
03975 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03976 }
03977 for (int iz = 1; iz <= m_vnzp; iz++) {
03978 if ( iz-1 > m_vnzc ) kz = iz-1-m_vnzp; else kz = iz-1;
03979 argz = float(kz*kz)*dz2;
03980 for (int iy = 1; iy <= m_vnyp; iy++) {
03981 if ( iy-1 > m_vnyc ) ky = iy-1-m_vnyp; else ky = iy-1;
03982 argy = argz + float(ky*ky)*dy2;
03983 for (int ix = 0; ix <= m_vnxc; ix++) {
03984 float Kn = (*m_wptr3)(ix,iy,iz);
03985 argx = std::sqrt(argy + float(ix*ix)*dx2);
03986 int r = Util::round(float(inc)*argx);
03987 if ( r >= 0 && Kn > 4.5f ) {
03988 if ( m_weighting == ESTIMATE ) {
03989 int cx = ix;
03990 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03991 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03992 float sum = 0.0;
03993 for( int ii = -kc; ii <= kc; ++ii ) {
03994 int nbrcx = cx + ii;
03995 if( nbrcx >= m_vnxc ) continue;
03996 for ( int jj= -kc; jj <= kc; ++jj ) {
03997 int nbrcy = cy + jj;
03998 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03999 for( int kk = -kc; kk <= kc; ++kk ) {
04000 int nbrcz = cz + jj;
04001 if ( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
04002 if( nbrcx < 0 ) {
04003 nbrcx = -nbrcx;
04004 nbrcy = -nbrcy;
04005 nbrcz = -nbrcz;
04006 }
04007 int nbrix = nbrcx;
04008 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
04009 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
04010 if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
04011 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
04012 sum = sum + pow_a[c];
04013 }
04014 }
04015 }
04016 }
04017
04018 wght = 1.0f / ( 1.0f - alpha * sum );
04019 }
04020 float nominator = std::norm(m_volume->cmplx(ix,iy,iz))/(*m_wptr)(ix,iy,iz);
04021 float denominator = ((*m_wptr2)(ix,iy,iz)-std::norm(m_volume->cmplx(ix,iy,iz))/(*m_wptr)(ix,iy,iz))/(Kn-1.0f);
04022
04023 if( (ix>0 || (kz>=0 && (ky>=0 || kz!=0)))) {
04024 if ( r <= inc ) {
04025 nom[r] += nominator*wght;
04026 denom[r] += denominator/(*m_wptr)(ix,iy,iz)*wght;
04027 nn[r] += 2;
04028 ka[r] += int(Kn);
04029 }
04030
04031
04032
04033
04034
04035
04036
04037
04038
04039
04040
04041
04042
04043
04044
04045 }
04046 (*vol_ssnr)(2*ix, iy-1, iz-1) = denominator*wght;
04047 }
04048 }
04049 }
04050 }
04051 for (int i = 0; i <= inc; i++) {
04052 (*SSNR)(i,0,0) = nom[i];
04053 (*SSNR)(i,1,0) = denom[i];
04054 (*SSNR)(i,2,0) = static_cast<float>(nn[i]);
04055 (*SSNR)(i,3,0) = static_cast<float>(ka[i]);
04056 }
04057 vol_ssnr->update();
04058
04059 delete[] nom;
04060 delete[] denom;
04061 delete[] nn;
04062 delete[] ka;
04063
04064 return 0;
04065 }
04066 #undef tw
04067
04068
04069
04070
04071 void EMAN::dump_reconstructors()
04072 {
04073 dump_factory < Reconstructor > ();
04074 }
04075
04076 map<string, vector<string> > EMAN::dump_reconstructors_list()
04077 {
04078 return dump_factory_list < Reconstructor > ();
04079 }
04080
04081
04082
04083
04084
04085 using std::ofstream;
04086 using std::ifstream;
04087
04088
04089 newfile_store::newfile_store( const string& filename, int npad, bool ctf )
04090 : m_bin_file( filename + ".bin" ),
04091 m_txt_file( filename + ".txt" )
04092 {
04093 m_npad = npad;
04094 m_ctf = ctf;
04095 }
04096
04097 newfile_store::~newfile_store( )
04098 {
04099 }
04100
04101 void newfile_store::add_image( EMData* emdata, const Transform& tf )
04102 {
04103 if( m_bin_of == NULL ) {
04104 m_bin_of = shared_ptr<ofstream>( new ofstream(m_bin_file.c_str(), std::ios::out|std::ios::binary) );
04105 m_txt_of = shared_ptr<ofstream>( new ofstream(m_txt_file.c_str()) );
04106 }
04107
04108
04109 EMData* padfft = padfft_slice( emdata, tf, m_npad );
04110
04111 int nx = padfft->get_xsize();
04112 int ny = padfft->get_ysize();
04113 int n2 = ny / 2;
04114 int n = ny;
04115
04116 float voltage=0.0f, pixel=0.0f, Cs=0.0f, ampcont=0.0f, bfactor=0.0f, defocus=0.0f;
04117
04118 if( m_ctf ) {
04119 Ctf* ctf = emdata->get_attr( "ctf" );
04120 Dict params = ctf->to_dict();
04121 voltage = params["voltage"];
04122 pixel = params["apix"];
04123 Cs = params["cs"];
04124 ampcont = params["ampcont"];
04125 bfactor = params["bfactor"];
04126 defocus = params["defocus"];
04127 if(ctf) {delete ctf; ctf=0;}
04128 }
04129
04130 vector<point_t> points;
04131 for( int j=-ny/2+1; j <= ny/2; j++ ) {
04132 int jp = (j>=0) ? j+1 : ny+j+1;
04133 for( int i=0; i <= n2; ++i ) {
04134 int r2 = i*i + j*j;
04135 if( (r2<ny*ny/4) && !( (i==0) && (j<0) ) ) {
04136 float ctf;
04137 if( m_ctf ) {
04138 float ak = std::sqrt( r2/float(ny*ny) )/pixel;
04139 ctf = Util::tf( defocus, ak, voltage, Cs, ampcont, bfactor, 1);
04140 } else {
04141 ctf = 1.0;
04142 }
04143
04144 float xnew = i*tf[0][0] + j*tf[1][0];
04145 float ynew = i*tf[0][1] + j*tf[1][1];
04146 float znew = i*tf[0][2] + j*tf[1][2];
04147 std::complex<float> btq;
04148 if (xnew < 0.) {
04149 xnew = -xnew;
04150 ynew = -ynew;
04151 znew = -znew;
04152 btq = conj(padfft->cmplx(i,jp-1));
04153 } else {
04154 btq = padfft->cmplx(i,jp-1);
04155 }
04156
04157 int ixn = int(xnew + 0.5 + n) - n;
04158 int iyn = int(ynew + 0.5 + n) - n;
04159 int izn = int(znew + 0.5 + n) - n;
04160 if ((ixn <= n2) && (iyn >= -n2) && (iyn <= n2) && (izn >= -n2) && (izn <= n2)) {
04161 int ixf, iyf, izf;
04162 if (ixn >= 0) {
04163 int iza, iya;
04164 if (izn >= 0)
04165 iza = izn + 1;
04166 else
04167 iza = n + izn + 1;
04168
04169 if (iyn >= 0)
04170 iya = iyn + 1;
04171 else
04172 iya = n + iyn + 1;
04173
04174 ixf = ixn;
04175 iyf = iya;
04176 izf = iza;
04177 } else {
04178 int izt, iyt;
04179 if (izn > 0)
04180 izt = n - izn + 1;
04181 else
04182 izt = -izn + 1;
04183
04184 if (iyn > 0)
04185 iyt = n - iyn + 1;
04186 else
04187 iyt = -iyn + 1;
04188
04189 ixf = -ixn;
04190 iyf = iyt;
04191 izf = izt;
04192 }
04193
04194
04195 int pos2 = ixf + (iyf-1)*nx/2 + (izf-1)*ny*nx/2;
04196 float ctfv1 = btq.real() * ctf;
04197 float ctfv2 = btq.imag() * ctf;
04198 float ctf2 = ctf*ctf;
04199
04200 point_t p;
04201 p.pos2 = pos2;
04202 p.real = ctfv1;
04203 p.imag = ctfv2;
04204 p.ctf2 = ctf2;
04205
04206 points.push_back( p );
04207 }
04208 }
04209 }
04210 }
04211
04212
04213 int npoint = points.size();
04214 std::istream::off_type offset = (m_offsets.size()==0) ? 0 : m_offsets.back();
04215 offset += npoint*sizeof(point_t);
04216 m_offsets.push_back( offset );
04217
04218 *m_txt_of << m_offsets.back() << std::endl;
04219 m_bin_of->write( (char*)(&points[0]), sizeof(point_t)*npoint );
04220 checked_delete( padfft );
04221 }
04222
04223 void newfile_store::get_image( int id, EMData* buf )
04224 {
04225 if( m_offsets.size()==0 ) {
04226 ifstream is( m_txt_file.c_str() );
04227 std::istream::off_type off;
04228 while( is >> off ) {
04229 m_offsets.push_back( off );
04230 }
04231
04232 m_bin_if = shared_ptr<std::ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
04233 }
04234
04235 Assert( m_bin_if != NULL );
04236
04237 std::istream::off_type offset = (id==0) ? 0 : m_offsets[id-1];
04238 Assert( offset >= 0 );
04239 m_bin_if->seekg( offset, std::ios::beg );
04240
04241
04242 if( m_bin_if->bad() || m_bin_if->fail() || m_bin_if->eof() ) {
04243 std::cout << "bad or fail or eof while fetching id, offset: " << id << " " << offset << std::endl;
04244 throw std::logic_error( "bad happen" );
04245 }
04246
04247 int bufsize = (m_offsets[id] - offset)/sizeof(float);
04248 if( buf->get_xsize() != bufsize ) {
04249 buf->set_size( bufsize, 1, 1 );
04250 }
04251
04252 char* data = (char*)(buf->get_data());
04253 m_bin_if->read( data, sizeof(float)*bufsize );
04254 buf->update();
04255 }
04256
04257 void newfile_store::read( int nprj )
04258 {
04259 if( m_offsets.size()==0 ) {
04260 ifstream is( m_txt_file.c_str() );
04261 std::istream::off_type off;
04262 while( is >> off ) {
04263 m_offsets.push_back( off );
04264 }
04265 }
04266
04267 if( m_bin_if==NULL ) {
04268 m_bin_if = shared_ptr< ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
04269 }
04270
04271
04272 int npoint = m_offsets[0]/sizeof(point_t);
04273 std::ios::off_type prjsize = m_offsets[0];
04274
04275 try {
04276 m_points.resize(nprj * npoint);
04277 }
04278 catch( std::exception& e ) {
04279 std::cout << "Error: " << e.what() << std::endl;
04280 }
04281
04282 int ip = 0;
04283 for( int i=0; i < nprj; ++i ) {
04284 m_bin_if->read( (char*)(&m_points[ip]), prjsize );
04285 if( m_bin_if->bad() || m_bin_if->fail() || m_bin_if->eof() ) {
04286 std::cout << "Error: file hander bad or fail or eof" << std::endl;
04287 return;
04288 }
04289 ip += npoint;
04290 }
04291 }
04292
04293 void newfile_store::add_tovol( EMData* fftvol, EMData* wgtvol, const vector<int>& mults, int pbegin, int pend )
04294 {
04295 float* vdata = fftvol->get_data();
04296 float* wdata = wgtvol->get_data();
04297
04298 int npoint = m_offsets[0]/sizeof(point_t);
04299
04300 Assert( int(m_points.size())== (pend - pbegin)*npoint );
04301
04302 for( int iprj=pbegin; iprj < pend; ++iprj ) {
04303 int m = mults[iprj];
04304 if( m==0 ) continue;
04305
04306 int ipt = (iprj-pbegin)*npoint;
04307 for( int i=0; i < npoint; ++i ) {
04308 int pos2 = m_points[ipt].pos2;
04309 int pos1 = pos2*2;
04310
04311 wdata[pos2] += m_points[ipt].ctf2*m;
04312 vdata[pos1] += m_points[ipt].real*m;
04313 vdata[pos1+1] += m_points[ipt].imag*m;
04314 ++ipt;
04315 }
04316 }
04317 }
04318
04319 void newfile_store::restart()
04320 {
04321 m_bin_if = shared_ptr< ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
04322 }
04323
04324 file_store::file_store(const string& filename, int npad, int write, bool ctf)
04325 : m_bin_file(filename + ".bin"),
04326 m_txt_file(filename + ".txt")
04327 {
04328 m_ctf = ctf;
04329 m_prev = -1;
04330 m_npad = npad;
04331 m_write = write;
04332 }
04333
04334 file_store::~file_store()
04335 {
04336 }
04337
04338 void file_store::add_image( EMData* emdata, const Transform& tf )
04339 {
04340
04341 EMData* padfft = padfft_slice( emdata, tf, m_npad );
04342
04343 float* data = padfft->get_data();
04344
04345 if( m_write && m_bin_ohandle == NULL )
04346 {
04347 m_bin_ohandle = shared_ptr< ofstream >( new ofstream(m_bin_file.c_str(), std::ios::out | std::ios::binary) );
04348 m_txt_ohandle = shared_ptr< ofstream >( new ofstream(m_txt_file.c_str() ) );
04349 if( m_ctf )
04350 *m_txt_ohandle << "Cs pixel voltage ctf_applied amp_contrast defocus ";
04351
04352 *m_txt_ohandle << "phi theta psi" << std::endl;
04353 }
04354
04355 m_x_out = padfft->get_xsize();
04356 m_y_out = padfft->get_ysize();
04357 m_z_out = padfft->get_zsize();
04358 m_totsize = m_x_out*m_y_out*m_z_out;
04359
04360 if( m_ctf )
04361 {
04362 Ctf* ctf = padfft->get_attr( "ctf" );
04363 Dict ctf_params = ctf->to_dict();
04364
04365 m_ctf_applied = padfft->get_attr( "ctf_applied" );
04366
04367 m_Cs = ctf_params["cs"];
04368 m_pixel = ctf_params["apix"];
04369 m_voltage = ctf_params["voltage"];
04370 m_amp_contrast = ctf_params["ampcont"];
04371 m_defocuses.push_back( ctf_params["defocus"] );
04372 if(ctf) {delete ctf; ctf=0;}
04373 }
04374
04375 Dict params = tf.get_rotation( "spider" );
04376 float phi = params.get( "phi" );
04377 float tht = params.get( "theta" );
04378 float psi = params.get( "psi" );
04379
04380
04381 m_phis.push_back( phi );
04382 m_thetas.push_back( tht );
04383 m_psis.push_back( psi );
04384
04385 if( m_write )
04386 {
04387 m_bin_ohandle->write( (char*)data, sizeof(float)*m_totsize );
04388
04389 if( m_ctf )
04390 {
04391 *m_txt_ohandle << m_Cs << " ";
04392 *m_txt_ohandle << m_pixel << " ";
04393 *m_txt_ohandle << m_voltage << " ";
04394 *m_txt_ohandle << m_ctf_applied << " ";
04395 *m_txt_ohandle << m_amp_contrast << " ";
04396 *m_txt_ohandle << m_defocuses.back() << " ";
04397 }
04398 *m_txt_ohandle << m_phis.back() << " ";
04399 *m_txt_ohandle << m_thetas.back() << " ";
04400 *m_txt_ohandle << m_psis.back() << " ";
04401 *m_txt_ohandle << m_x_out << " ";
04402 *m_txt_ohandle << m_y_out << " ";
04403 *m_txt_ohandle << m_z_out << " ";
04404 *m_txt_ohandle << m_totsize << std::endl;
04405 }
04406
04407 checked_delete(padfft);
04408
04409 }
04410
04411 void file_store::get_image( int id, EMData* padfft )
04412 {
04413
04414 if( m_phis.size() == 0 ) {
04415 ifstream m_txt_ifs( m_txt_file.c_str() );
04416
04417 if( !m_txt_ifs )
04418 {
04419 std::cerr << "Error: file " << m_txt_file << " does not exist" << std::endl;
04420 }
04421
04422 string line;
04423 std::getline( m_txt_ifs, line );
04424
04425 float first, defocus, phi, theta, psi;
04426
04427
04428
04429 while( m_txt_ifs >> first ) {
04430
04431 if( m_ctf )
04432 {
04433 m_Cs = first;
04434 m_txt_ifs >> m_pixel >> m_voltage;
04435 m_txt_ifs >> m_ctf_applied >> m_amp_contrast;
04436 m_txt_ifs >> defocus >> phi >> theta >> psi;
04437 m_defocuses.push_back( defocus );
04438 }
04439 else
04440 {
04441 phi = first;
04442 m_txt_ifs >> theta >> psi;
04443 }
04444
04445 m_txt_ifs >> m_x_out >> m_y_out >> m_z_out >> m_totsize;
04446 m_phis.push_back( phi );
04447 m_thetas.push_back( theta );
04448 m_psis.push_back( psi );
04449 }
04450 }
04451
04452 Assert( m_ihandle != NULL );
04453
04454 std::istream::off_type offset = id*sizeof(float)*m_totsize;
04455 Assert( offset >= 0 );
04456
04457 if( offset > 0 )
04458 {
04459 m_ihandle->seekg(offset, std::ios::beg);
04460 }
04461
04462 if( m_ihandle->bad() )
04463 {
04464 std::cout << "bad while fetching id, offset: " << id << " " << offset << std::endl;
04465 throw std::logic_error( "bad happen" );
04466 }
04467
04468 if( m_ihandle->fail() )
04469 {
04470 std::cout << "fail while fetching id, offset, curoff: " << id << " " << offset << std::endl;
04471 throw std::logic_error( "fail happen" );
04472 }
04473
04474 if( m_ihandle->eof() )
04475 {
04476 std::cout << "eof while fetching id, offset: " << id << " " << offset << std::endl;
04477 throw std::logic_error( "eof happen" );
04478 }
04479
04480 if( padfft->get_xsize() != m_x_out ||
04481 padfft->get_ysize() != m_y_out ||
04482 padfft->get_zsize() != m_z_out )
04483 {
04484 padfft->set_size(m_x_out, m_y_out, m_z_out);
04485 }
04486
04487 char* data = (char*)(padfft->get_data());
04488 m_ihandle->read( data, sizeof(float)*m_totsize );
04489 padfft->update();
04490
04491 if( m_ctf )
04492 {
04493 padfft->set_attr( "Cs", m_Cs );
04494 padfft->set_attr( "Pixel_size", m_pixel );
04495 padfft->set_attr( "voltage", m_voltage );
04496 padfft->set_attr( "ctf_applied", m_ctf_applied );
04497 padfft->set_attr( "amp_contrast", m_amp_contrast );
04498 padfft->set_attr( "defocus", m_defocuses[id] );
04499 }
04500
04501 padfft->set_attr( "padffted", 1 );
04502 padfft->set_attr( "phi", m_phis[id] );
04503 padfft->set_attr( "theta", m_thetas[id] );
04504 padfft->set_attr( "psi", m_psis[id] );
04505
04506 }
04507
04508 void file_store::restart( )
04509 {
04510 if( m_ihandle == NULL )
04511 {
04512 m_ihandle = shared_ptr< ifstream >( new ifstream(m_bin_file.c_str(), std::ios::in | std::ios::binary) );
04513 }
04514
04515 if( m_ihandle->bad() || m_ihandle->fail() || m_ihandle->eof() )
04516 {
04517 m_ihandle->open( m_bin_file.c_str(), std::ios::binary );
04518 }
04519
04520 m_ihandle->seekg( 0, std::ios::beg );
04521 }
04522
04523