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