41#include <gsl/gsl_statistics_double.h>
42#include <gsl/gsl_fit.h>
44#ifdef EMAN2_USING_CUDA
63using std::setprecision;
67 typedef char type_must_be_complete[
sizeof(T)? 1: -1 ];
68 (void)
sizeof(type_must_be_complete);
73const string FourierReconstructor::NAME =
"fourier";
74const string FourierIterReconstructor::NAME =
"fourier_iter";
75const string FourierReconstructorSimple2D::NAME =
"fouriersimple2D";
76const string WienerFourierReconstructor::NAME =
"wiener_fourier";
77const string BackProjectionReconstructor::NAME =
"back_projection";
78const string RealMedianReconstructor::NAME =
"real_median";
79const string nn4Reconstructor::NAME =
"nn4";
80const string nn4_rectReconstructor::NAME =
"nn4_rect";
81const string nnSSNR_Reconstructor::NAME =
"nnSSNR";
82const string nn4_ctfReconstructor::NAME =
"nn4_ctf";
83const string nn4_ctfwReconstructor::NAME =
"nn4_ctfw";
84const string nn4_ctfwsReconstructor::NAME =
"nn4_ctfws";
85const string nn4_ctf_rectReconstructor::NAME =
"nn4_ctf_rect";
86const string nnSSNR_ctfReconstructor::NAME =
"nnSSNR_ctf";
90 force_add<FourierReconstructor>();
91 force_add<FourierIterReconstructor>();
92 force_add<FourierReconstructorSimple2D>();
94 force_add<WienerFourierReconstructor>();
95 force_add<RealMedianReconstructor>();
96 force_add<BackProjectionReconstructor>();
97 force_add<nn4Reconstructor>();
98 force_add<nn4_rectReconstructor>();
99 force_add<nnSSNR_Reconstructor>();
100 force_add<nn4_ctfReconstructor>();
101 force_add<nn4_ctfwReconstructor>();
102 force_add<nn4_ctfwsReconstructor>();
103 force_add<nn4_ctf_rectReconstructor>();
104 force_add<nnSSNR_ctfReconstructor>();
112 static void init(
int winsize,
const Ctf* ctf ) {
123 m_dza = params[
"dfdiff"];
124 m_azz = params[
"dfang"];
140 return Util::ctf_img_real(
m_winsize,
m_winsize, 1,
m_defocus,
m_pixel,
m_voltage,
m_cs,
m_ampcont,
m_bfactor,
m_dza,
m_azz, 1);
163void FourierReconstructorSimple2D::setup()
176 image->set_complex(
true);
195 EMData* working_slice = slice->process(
"xform.phaseorigin.tocorner");
198 working_slice->do_fft_inplace();
202 float* dat = working_slice->get_data();
207 offset[0] = 0; offset[1] = 2; offset[2] =
nx; offset[3] =
nx+2;
209 float alt = -((float)(euler.
get_rotation(
"2d"))[
"alpha"])*M_PI/180.0f;
210 for (
int x = 0;
x < working_slice->get_xsize() / 2;
x++) {
212 float rx = (float)
x;
214 float xx = rx*cos(alt);
215 float yy = rx*sin(alt);
228 dt[1] = cc * dat[2*
x+1];
231 int x0 = (int) floor(xx);
232 int y0 = (int) floor(yy);
234 int i = 2*x0 + y0*
nx;
246 int k = i + offset[0];
247 rdata[k] += g[0] * dt[0];
248 rdata[k + 1] += g[0] * dt[1];
252 rdata[k] += g[2] * dt[0];
253 rdata[k + 1] += g[2] * dt[1];
260 int dif = x0 - (
nx-2);
265 int k = i + offset[0];
266 rdata[k] += g[0] * dt[0];
267 rdata[k + 1] += g[0] * dt[1];
271 rdata[k] += g[1] * dt[0];
272 rdata[k + 1] += g[1] * dt[1];
278 int dif = y0 - (
ny-1);
282 if (x0 >=
nx - 2 || y0 >=
ny - 1)
continue;
287 for (
int j = 0; j < 4; j++)
289 int k = i + offset[j];
290 rdata[k] += g[j] * dt[0];
291 rdata[k + 1] += g[j] * dt[1];
305 image->process_inplace(
"xform.fourierorigin.tocorner");
306 image->do_ift_inplace();
308 image->process_inplace(
"xform.phaseorigin.tocenter");
322 size_t nnxy=
tmp_data->get_ysize()*nnx;
331 std::complex<float> pv =
image->get_complex_at(113,0,23);
332 float pnv =
tmp_data->get_value_at(113,0,23);
333 printf(
"norm: %1.4g\t%1.4g\t%1.4g\n",pv.real()/pnv,pv.imag()/pnv,pnv);
341 for (
int z=-
nz/2; z<
nz/2; z++) {
342 for (
int y=0;
y<=
ny/2;
y++) {
343 if (z<=0 && (
y==0||
y==
ny/2))
continue;
345 size_t idx1=(
y==0?0:
ny-
y)*nnx+(z<=0?-z:
nz-z)*nnxy;
346 size_t idx2=
y*nnx+(z<0?
nz+z:z)*nnxy;
347 if (idx1==idx2)
continue;
348 float snorm=norm[idx1]+norm[idx2];
349 norm[idx1]=norm[idx2]=snorm;
354 snorm=norm[idx1]+norm[idx2];
355 norm[idx1]=norm[idx2]=snorm;
359 norm[0 + 0*nnx+
nz/2*nnxy]*=2;
360 norm[0 +
ny/2*nnx+ 0*nnxy]*=2;
361 norm[0 +
ny/2*nnx+
nz/2*nnxy]*=2;
362 norm[
nx/2-1+ 0*nnx+
nz/2*nnxy]*=2;
363 norm[
nx/2-1+
ny/2*nnx+ 0*nnxy]*=2;
372 if (wiener) wfilt=1.0;
377 if (sqrtnorm) d=
sqrt(d);
386 rdata[i + 1] /= d+wfilt;
434 vector<int> size=
params[
"size"];
451 image->set_complex(
true);
452 image->set_fftodd(
false);
473 if ( (
bool)
params[
"verbose"] )
475 cout <<
"3D Fourier dimensions are " <<
nx <<
" " <<
ny <<
" " <<
nz << endl;
476 cout <<
"3D Fourier subvolume is " <<
subnx <<
" " <<
subny <<
" " <<
subnz << endl;
477 printf (
"You will require approximately %1.3g GB of memory to reconstruct this volume\n",((
float)
subnx*
subny*
subnz*
sizeof(
float)*2.5)/1000000000.0);
483 if (seed->is_complex())
ref_vol=seed->copy();
486 ref_vol->process_inplace(
"xform.phaseorigin.tocorner");
492 if (seed->is_complex())
ref_vol=seed->copy();
495 ref_vol->process_inplace(
"xform.phaseorigin.tocorner");
497 printf(
"Warning: FourierIterReconstructor ignores seed weights");
512 else return_slice = slice->process(
"xform",
Dict(
"transform",&tmp));
514 return_slice->process_inplace(
"xform.phaseorigin.tocorner");
515 return_slice->do_fft_inplace();
516 return_slice->mult((
float)
sqrt(1.0f/(return_slice->get_ysize())*return_slice->get_xsize()));
518 return_slice->set_attr(
"reconstruct_preproc",(
int)1);
535 float inx=(float)(slice->get_xsize());
536 float iny=(float)(slice->get_ysize());
537 size_t ny2sq=
ny*
ny/4;
539 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
541 for (
int y = -iny/2;
y < iny/2;
y++) {
542 for (
int x = 0;
x < inx/2;
x++) {
544 float rx = (float)
x/(inx-2.0f);
545 float ry = (float)
y/iny;
547 if (r>iny/2+3 && abs(inx-iny)<3)
continue;
549 Vec3f coord(rx,ry,0);
561 int x0 = (int) floor(xx-2.5);
562 int y0 = (int) floor(yy-2.5);
563 int z0 = (int) floor(zz-2.5);
565 std::complex<float>dt=slice->get_complex_at(
x,
y);
567 if (x0<-
nx2-4 || y0<-
ny2-4 || z0<-nz2-4 || x0>
nx2+3 || y0>
ny2+3 || z0>
nz2+3 )
continue;
589 ref_vol->get_complex_at(xx0,yy0,zz0),
590 ref_vol->get_complex_at(xx0+1,yy0,zz0),
591 ref_vol->get_complex_at(xx0,yy0+1,zz0),
592 ref_vol->get_complex_at(xx0+1,yy0+1,zz0),
593 ref_vol->get_complex_at(xx0,yy0,zz0+1),
594 ref_vol->get_complex_at(xx0+1,yy0,zz0+1),
595 ref_vol->get_complex_at(xx0,yy0+1,zz0+1),
596 ref_vol->get_complex_at(xx0+1,yy0+1,zz0+1),
597 xx-xx0,yy-yy0,zz-zz0);
602 for (
int k = z0 ; k <= z1; k++) {
603 for (
int j = y0 ; j <= y1; j++) {
604 for (
int i = x0; i <= x1; i ++) {
605 if ((
size_t)i*i+j*j+k*k>ny2sq)
continue;
606 std::complex<float>ntdt=
ref_vol->get_complex_at(i,j,k);
611 std::complex<float>updt=dt+ntdt-nodt;
617 off=
image->add_complex_at_fast(i,j,k,updt*gg*weight);
618 tmp_data->get_data()[off/2]+=gg*weight;
627 for (
int k = z0 ; k <= z1; k++) {
628 for (
int j = y0 ; j <= y1; j++) {
629 for (
int i = x0; i <= x1; i ++) {
634 off=
image->add_complex_at_fast(i,j,k,dt*gg*weight);
637 tmp_data->get_data()[off/2]+=gg*weight;
647 delete rotation; rotation=0;
657 image->do_ift_inplace();
659 image->process_inplace(
"xform.phaseorigin.tocenter");
665 if (
tmp_data->get_ysize()%2==0 &&
tmp_data->get_zsize()%2==0)
tmp_data->process_inplace(
"xform.fourierorigin.tocenter");
711 parms[
"data"] =
image;
712 parms[
"norm"] =
tmp_data->get_data();
736 vector<int> size=
params[
"size"];
746 ddata=(
double *)malloc(
sizeof(
double)*250);
747 dnorm=(
double *)malloc(
sizeof(
double)*250);
748 for (
int i=0; i<250; i++) ddata[i]=dnorm[i]=0.0;
783 image->set_complex(
true);
793 image->set_attr(
"subvolume_full_nx",
nx);
794 image->set_attr(
"subvolume_full_ny",
ny);
795 image->set_attr(
"subvolume_full_nz",
nz);
812 if ( (
bool)
params[
"verbose"] )
814 cout <<
"3D Fourier dimensions are " <<
nx <<
" " <<
ny <<
" " <<
nz << endl;
815 cout <<
"3D Fourier subvolume is " <<
subnx <<
" " <<
subny <<
" " <<
subnz << endl;
816 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);
824 vector<int> size=
params[
"size"];
861 if (seed->get_xsize()!=
subnx || seed->get_ysize()!=
subny || seed->get_zsize()!=
subnz || !seed->is_complex())
865 image = seed->copy();
866 image->mult(seed_weight);
874 image->set_attr(
"subvolume_full_nx",
nx);
875 image->set_attr(
"subvolume_full_ny",
ny);
876 image->set_attr(
"subvolume_full_nz",
nz);
888 for (
int k=0; k<
subnz; k++) {
889 for (
int j=0; j<
subny; j++) {
890 tmp_data->set_value_at(0,j,k,seed_weight/2.0);
897 for (
int k=0; k<
subnz; k++) {
898 for (
int j=0; j<
subny; j++) {
911 std::complex<float> pv =
image->get_complex_at(113,0,23);
912 float pnv =
tmp_data->get_value_at(113,0,23);
913 printf(
"seed: %1.4g\t%1.4g\t%1.4g\n",pv.real()/pnv,pv.imag()/pnv,pnv);
917 if ( (
bool)
params[
"quiet"] ==
false )
919 cout <<
"Seeded direct Fourier inversion";
920 cout <<
"3D Fourier dimensions are " <<
nx <<
" " <<
ny <<
" " <<
nz << endl;
921 cout <<
"3D Fourier subvolume is " <<
subnx <<
" " <<
subny <<
" " <<
subnz << endl;
922 cout <<
"You will require approximately " << setprecision(3) << (
subnx*
subny*
subnz*
sizeof(float)*1.5)/1000000000.0 <<
"GB of memory to reconstruct this volume" << endl;
932 vector<int> size=
params[
"size"];
969 if (seed->get_xsize()!=
subnx || seed->get_ysize()!=
subny || seed->get_zsize()!=
subnz || !seed->is_complex())
973 image = seed->copy();
979 image->set_attr(
"subvolume_full_nx",
nx);
980 image->set_attr(
"subvolume_full_ny",
ny);
981 image->set_attr(
"subvolume_full_nz",
nz);
989 if ( (
bool)
params[
"quiet"] ==
false )
991 cout <<
"Seeded direct Fourier inversion";
992 cout <<
"3D Fourier dimensions are " <<
nx <<
" " <<
ny <<
" " <<
nz << endl;
993 cout <<
"3D Fourier subvolume is " <<
subnx <<
" " <<
subny <<
" " <<
subnz << endl;
994 cout <<
"You will require approximately " << setprecision(3) << (
subnx*
subny*
subnz*
sizeof(float)*1.5)/1000000000.0 <<
"GB of memory to reconstruct this volume" << endl;
1001 bool zeroimage =
true;
1002 bool zerotmpimg =
true;
1004#ifdef EMAN2_USING_CUDA
1005 if(EMData::usecuda == 1) {
1006 if(
image->getcudarwdata()) {
1017 if(zeroimage)
image->to_zero();
1018 if(zerotmpimg)
tmp_data->to_zero();
1024#ifdef EMAN2_USING_CUDA
1025 if(EMData::usecuda == 1) {
1026 if(!slice->getcudarwdata()) slice->copy_to_cuda();
1030 EMData* return_slice = 0;
1034 if (tmp.
is_identity()) return_slice=slice->copy();
1035 else return_slice = slice->process(
"xform",
Dict(
"transform",&tmp));
1037 return_slice->process_inplace(
"xform.phaseorigin.tocorner");
1042#ifdef EMAN2_USING_CUDA
1043 if(EMData::usecuda == 1 && return_slice->getcudarwdata()) {
1044 return_slice->do_fft_inplace_cuda();
1046 return_slice->do_fft_inplace();
1049 return_slice->do_fft_inplace();
1054 return_slice->mult((
float)
sqrt(1.0f/(return_slice->get_ysize())*return_slice->get_xsize()));
1059 return_slice->set_attr(
"reconstruct_preproc",(
int)1);
1060 return return_slice;
1074#ifdef EMAN2_USING_CUDA
1075 if(EMData::usecuda == 1) {
1076 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda();
1082 float weight=oweight;
1084 if (input_slice->has_attr(
"class_ssnr")) weight=-1.0;
1088 if (weight==0)
return -1;
1098 if (input_slice->get_attr_default(
"reconstruct_preproc",(
int) 0)) slice=input_slice->copy();
1114 delete rotation; rotation=0;
1136 float inx=(float)(input_slice->get_xsize());
1137 float iny=(float)(input_slice->get_ysize());
1139 if (abs(inx-iny)>2 && weight<0) printf(
"WARNING: Fourier Reconstruction failure. SSNR flag set with asymmetric dimensions on input image\n");
1141#ifdef EMAN2_USING_CUDA
1142 if(EMData::usecuda == 1) {
1143 if(!
image->getcudarwdata()){
1144 image->copy_to_cuda();
1147 float * m =
new float[12];
1148 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1159 float sscale = 1.0f;
1161 ssnr=input_slice->get_attr(
"class_ssnr");
1162 sscale=2.0*(ssnr.size()-1)/iny;
1165 float rweight=weight;
1166 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1168 for (
int y = -iny/2;
y < iny/2;
y++) {
1169 for (
int x = 0;
x < inx/2;
x++) {
1171 float rx = (float)
x/(inx-2.0f);
1172 float ry = (float)
y/iny;
1181 if (r>iny/2 && abs(inx-iny)<3)
continue;
1183 if (weight<0) rweight=
Util::get_max(0.0f,ssnr[
int(r*sscale)]);
1187 Vec3f coord(rx,ry,0);
1189 float xx = coord[0];
1190 float yy = coord[1];
1191 float zz = coord[2];
1217#ifdef EMAN2_USING_CUDA
1218 if(EMData::usecuda == 1) {
1219 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda();
1227 if (input_slice->get_attr_default(
"reconstruct_preproc",(
int) 0)) slice=input_slice->copy();
1242 input_slice->set_attr(
"reconstruct_norm",slice->get_attr(
"reconstruct_norm"));
1243 input_slice->set_attr(
"reconstruct_absqual",slice->get_attr(
"reconstruct_absqual"));
1244 input_slice->set_attr(
"reconstruct_absqual_lowres",slice->get_attr(
"reconstruct_absqual_lowres"));
1246 input_slice->set_attr(
"reconstruct_weight",slice->get_attr(
"reconstruct_weight"));
1264 float *dat = input_slice->get_data();
1267 float inx=(float)(input_slice->get_xsize());
1268 float iny=(float)(input_slice->get_ysize());
1269 float apix=input_slice->get_attr(
"apix_x");
1270 float rmax=iny*apix/12.0;
1271 if (rmax>iny/2.0-1.0) rmax=iny/2.0-1.0;
1272 if (rmax<5.0) rmax=5.0;
1283 bool use_cpu =
true;
1287#ifdef EMAN2_USING_CUDA
1288 if(EMData::usecuda == 1) {
1289 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda();
1290 if(!
image->getcudarwdata()){
1291 image->copy_to_cuda();
1294 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1296 float * m =
new float[12];
1309 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1311 for (
int y = -iny/2;
y < iny/2;
y++) {
1312 for (
int x = 0;
x <= inx/2;
x++) {
1313 if (
x==0 &&
y==0)
continue;
1315 float rx = (float)
x/(inx-2);
1316 float ry = (float)
y/iny;
1321 Vec3f coord(rx,ry,0);
1323 float xx = coord[0];
1324 float yy = coord[1];
1325 float zz = coord[2];
1328 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5)
continue;
1336 int idx = (int)(
x * 2 + inx*(
y<0?iny+
y:
y));
1338 dt2[1] = dat[idx+1];
1341 if (!
pixel_at(xx,yy,zz,dt) || dt[2]==0)
continue;
1345 if (r>3 && r<rmax) {
1346 dotlow+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
1347 dlpow+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
1348 dlpow2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
1354 dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
1355 power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
1356 power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
1370 if (dlpow*dlpow2>0) dotlow/=
sqrt(dlpow*dlpow2);
1372 input_slice->set_attr(
"reconstruct_norm",(
float)(power2<=0?1.0:
sqrt(
power/power2)));
1374 input_slice->set_attr(
"reconstruct_absqual",(
float)
dot);
1375 input_slice->set_attr(
"reconstruct_absqual_lowres",(
float)dotlow);
1376 float rw=weight<=0?1.0f:1.0f/weight;
1377 input_slice->set_attr(
"reconstruct_qual",(
float)(
dot*rw/((rw-1.0)*
dot+1.0)));
1378 input_slice->set_attr(
"reconstruct_weight",(
float)(vweight/(
float)(nval)));
1389 ret->set_complex(1);
1391 ret->set_attr(
"apix_x",
image->get_attr(
"apix_x"));
1392 ret->set_attr(
"apix_y",
image->get_attr(
"apix_y"));
1393 ret->set_attr(
"apix_z",
image->get_attr(
"apix_z"));
1407 for (
int y = -
ny/2;
y <
ny/2;
y++) {
1408 for (
int x = 0;
x <=
nx/2;
x++) {
1410 float rx = (float)
x/(
nx-2);
1411 float ry = (float)
y/
ny;
1414 Vec3f coord(rx,ry,0);
1415 coord = coord*rotation;
1416 float xx = coord[0];
1417 float yy = coord[1];
1418 float zz = coord[2];
1421 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5)
continue;
1428 if (!
pixel_at(xx,yy,zz,dt) || dt[2]<0.005)
continue;
1430 if ((
x==0 ||
x==
nx/2) && (
y!=0 &&
y!=-
ny/2)) { dt[0]/=2; dt[1]/=2; }
1431 ret->set_complex_at(
x,
y,std::complex<float>(dt[0],dt[1]));
1437 ret->process_inplace(
"xform",
Dict(
"transform",&translation));
1439 ret->process_inplace(
"xform.phaseorigin.tocenter");
1440 EMData *tmp=ret->do_ift();
1454 std::complex<float> val=
image->get_complex_at(x0,y0,z0);
1455 size_t idx=
image->get_complex_index_fast(x0,y0,z0)/2;
1456 float norm=
tmp_data->get_value_at_index(idx);
1457 dt[0]=val.real()/norm;
1458 dt[1]=val.imag()/norm;
1556#ifdef EMAN2_USING_CUDA
1557 if(EMData::usecuda == 1 &&
image->getcudarwdata()){
1558 cout <<
"copy back from CUDA" << endl;
1559 image->copy_from_device();
1570 for (
int k=0; k<5; k++) {
1571 for (
int j=0; j<5; j++) {
1572 for (
int i=0; i<5; i++) {
1573 int idx=i*2+j*10+k*50;
1574 std::complex <double> a(ddata[idx],ddata[idx+1]);
1575 double b=dnorm[idx];
1577 printf(
"%d %d %d %1.4lg\t%1.4g %1.4lg\t%1.4g\n",i,j,k,a.real(),
image->get_value_at(i*2,j,k),a.imag(),
image->get_value_at(i*2+1,j,k));
1598 image->do_ift_inplace();
1600 image->process_inplace(
"xform.phaseorigin.tocenter");
1624 normout->set_data(
tmp_data->copy()->get_data());
1629 if (
tmp_data->get_ysize()%2==0 &&
tmp_data->get_zsize()%2==0)
tmp_data->process_inplace(
"xform.fourierorigin.tocenter");
1655 if (!input_slice->has_attr(
"ctf_snr_total"))
1656 throw NotExistingObjectException(
"ctf_snr_total",
"No SNR information present in class-average. Must use the ctf.auto or ctfw.auto averager.");
1659 if (input_slice->get_attr_default(
"reconstruct_preproc",(
int) 0)) slice=input_slice->copy();
1672 delete rotation; rotation=0;
1684 float inx=(float)(input_slice->get_xsize());
1685 float iny=(float)(input_slice->get_ysize());
1687 int undo_wiener=(int)input_slice->get_attr_default(
"ctf_wiener_filtered",0);
1690 vector<float> snr=input_slice->get_attr(
"ctf_snr_total");
1692 if (inweight<0)
sub=-1.0;
1695 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1697 for (
int y = -iny/2;
y < iny/2;
y++) {
1698 for (
int x = 0;
x <= inx/2;
x++) {
1700 float rx = (float)
x/(inx-2.0f);
1701 float ry = (float)
y/iny;
1704 float rn = (float)hypot(rx,ry);
1705 if (rn>=.5)
continue;
1706 rn*=snr.size()*2.0f;
1707 int rni=(int)floor(rn);
1708 if ((
unsigned int)rni>=snr.size()-1) weight=snr[snr.size()-1]*
sub;
1715 Vec3f coord(rx,ry,0);
1717 float xx = coord[0];
1718 float yy = coord[1];
1719 float zz = coord[2];
1743 if (input_slice->get_attr_default(
"reconstruct_preproc",(
int) 0)) slice=input_slice->copy();
1761 input_slice->set_attr(
"reconstruct_norm",slice->get_attr(
"reconstruct_norm"));
1762 input_slice->set_attr(
"reconstruct_absqual",slice->get_attr(
"reconstruct_absqual"));
1764 input_slice->set_attr(
"reconstruct_weight",slice->get_attr(
"reconstruct_weight"));
1783 float *dat = input_slice->get_data();
1786 float inx=(float)(input_slice->get_xsize());
1787 float iny=(float)(input_slice->get_ysize());
1793 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1795 for (
int y = -iny/2;
y < iny/2;
y++) {
1796 for (
int x = 0;
x <= inx/2;
x++) {
1797 if (
x==0 &&
y==0)
continue;
1799 float rx = (float)
x/(inx-2);
1800 float ry = (float)
y/iny;
1805 Vec3f coord(rx,ry,0);
1807 float xx = coord[0];
1808 float yy = coord[1];
1809 float zz = coord[2];
1812 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5)
continue;
1820 int idx = (int)(
x * 2 + inx*(
y<0?iny+
y:
y));
1822 dt2[1] = dat[idx+1];
1825 if (!
pixel_at(xx,yy,zz,dt) || dt[2]<=0)
continue;
1828 dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
1830 power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
1831 power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
1838 input_slice->set_attr(
"reconstruct_norm",(
float)(power2<=0?1.0:
sqrt(
power/power2)));
1839 input_slice->set_attr(
"reconstruct_absqual",(
float)
dot);
1840 float rw=weight<=0?1.0f:1.0f/weight;
1841 input_slice->set_attr(
"reconstruct_qual",(
float)(
dot*rw/((rw-1.0)*
dot+1.0)));
1842 input_slice->set_attr(
"reconstruct_weight",(
float)vweight/(
float)(
subnx*
subny*
subnz));
1849 int x0 = (int) floor(xx);
1850 int y0 = (int) floor(yy);
1851 int z0 = (int) floor(zz);
1855 float normsum=0,normsum2=0;
1857 dt[0]=dt[1]=dt[2]=0.0;
1860 if (x0<-
nx2-1 || y0<-
ny2-1 || z0<-nz2-1 || x0>
nx2 || y0>
ny2 || z0>
nz2 )
return false;
1875 for (
int k = z0 ; k <= z1; k++) {
1876 for (
int j = y0 ; j <= y1; j++) {
1877 for (
int i = x0; i <= x1; i ++) {
1879 idx=
image->get_complex_index_fast(i,j,k);
1882 dt[0]+=gg*
rdata[idx];
1883 dt[1]+=(i<0?-1.0f:1.0f)*gg*
rdata[idx+1];
1884 dt[2]+=norm[idx/2]*gg;
1886 normsum+=gg*norm[idx/2];
1890 if (normsum==0)
return false;
1900 for (
int k = z0 ; k <= z0 + 1; k++) {
1901 for (
int j = y0 ; j <= y0 + 1; j++) {
1902 for (
int i = x0; i <= x0 + 1; i ++) {
1907 dt[0]+=gg*
rdata[idx];
1908 dt[1]+=(i<0?-1.0f:1.0f)*gg*
rdata[idx+1];
1915 if (normsum==0)
return false;
1928 image->do_ift_inplace();
1930 image->process_inplace(
"xform.phaseorigin.tocenter");
1936 if (
tmp_data->get_ysize()%2==0 &&
tmp_data->get_zsize()%2==0)
tmp_data->process_inplace(
"xform.fourierorigin.tocenter");
2510 vector<int> size=
params[
"size"];
2518 for (std::vector<EMData *>::iterator it =
slices.begin() ; it !=
slices.end(); ++it)
delete *it;
2533 return_slice = slice->copy();
2547 return return_slice;
2553 LOGERR(
"try to insert NULL slice");
2557 if (input->get_xsize() != input->get_ysize() || input->get_xsize() !=
nx) {
2558 LOGERR(
"tried to insert image that was not correction dimensions");
2562 slices.push_back(input->copy());
2574 input_slice->set_attr(
"reconstruct_norm",1.0f);
2575 input_slice->set_attr(
"reconstruct_absqual",1.0f);
2576 input_slice->set_attr(
"reconstruct_weight",1.0f);
2588 for (
int z=0; z<
nz; z++) {
2589 for (
int y=0;
y<
ny;
y++) {
2590 for (
int x=0;
x<
nx;
x++) {
2591 std::vector<float> vals;
2592 std::vector<Transform>::iterator itt =
xforms.begin();
2593 for (std::vector<EMData *>::iterator it =
slices.begin() ; it !=
slices.end(); ++it,++itt) {
2594 Vec3f imvec = itt->transform(
x-
nx/2.0f,
y-
ny/2.0f,z-
nz/2.0f);
2597 float val=(*it)->sget_value_at_interp(imvec[0]+
nx/2,imvec[1]+
ny/2);
2598 if (val!=0) vals.push_back(val);
2612 if (vals.size()==0)
continue;
2629 std::sort(vals.begin(),vals.end());
2631 image->set_value_at(
x,
y,z,vals[vals.size()/2]);
2639 if (vals.size()<6) {
2641 for (
int i=0; i<vals.size(); i++) val+=vals[i];
2642 image->set_value_at(
x,
y,z,val/vals.size());
2646 std::sort(vals.begin(),vals.end());
2647 int v0=0,v1=vals.size(),vs=v1/2;
2652 float cen=(vals[
v0]+vals[v1-1])/2.0f;
2653 for (vs=
v0; vs<v1; vs++) {
if (vals[vs]>=cen)
break; }
2655 if (v1-vs>vs-
v0)
v0=vs;
2662 for (
int i=
v0; i<v1; i++) val+=vals[i];
2669 if (vals.size()<6) {
2671 for (
int i=0; i<vals.size(); i++) val+=vals[i];
2672 image->set_value_at(
x,
y,z,val/vals.size());
2675 std::sort(vals.begin(),vals.end());
2677 int p1=vals.size()/3-1;
2678 int p2=vals.size()*2/3;
2679 float v=(fabs(vals[p1])<fabs(vals[p2]))?(vals[p1]+vals[p1-1]+vals[p1+1])/3.0f:(vals[p2]+vals[p2-1]+vals[p2+1])/3.0f;
2686 if (vals.size()<6) {
2688 for (
int i=0; i<vals.size(); i++) val+=vals[i];
2689 image->set_value_at(
x,
y,z,val/vals.size());
2692 std::sort(vals.begin(),vals.end());
2694 int p1=vals.size()/3-1;
2695 int p2=vals.size()/2;
2696 int p3=vals.size()*2/3;
2697 float v=(vals[p2]-vals[0]<vals[vals.size()-1]-vals[p2])?(vals[p1]+vals[p1-1]+vals[p1+1])/3.0f:(vals[p3]+vals[p3-1]+vals[p3+1])/3.0f;
2704 if (vals.size()<10) {
2705 std::sort(vals.begin(),vals.end());
2706 image->set_value_at(
x,
y,z,vals[vals.size()/2]);
2712 float mean=0,sigma=0;
2713 for (
int i=0; i<n; i++) { mean+=vals[i]; sigma+=vals[i]*vals[i]; }
2715 sigma=
sqrt(sigma/n-mean);
2717 float low=(vals[0]+vals[1]+vals[2]+vals[3])/4.0f;
2718 float high=(vals[n-1]+vals[n-2]+vals[n-3]+vals[n-4])/4.0f;
2722 if (fabs(low-high)<sigma) {
2723 mean=(low+high)*4.0f;
2725 for (
int i=4; i<n-4; i++) {
2726 if (fabs(vals[i]-mean/c)<sigma) { mean+=vals[i]; c++; }
2728 image->set_value_at(
x,
y,z,mean/c);
2732 if (fabs(low)<fabs(high)) mean=low*4.0;
2735 for (
int i=4; i<n-4; i++) {
2736 if (fabs(vals[i]-mean/c)<sigma) { mean+=vals[i]; c++; }
2738 image->set_value_at(
x,
y,z,mean/c);
2760 vector<Transform> syms = sym->
get_syms();
2762 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
2767 image->add(*tmpcopy);
2780 for (std::vector<EMData *>::iterator it =
slices.begin() ; it !=
slices.end(); ++it)
delete *it;
2787 ret=
image->do_fft();
2788 ret->process_inplace(
"xform.phaseorigin.tocorner");
2799 vector<int> size=
params[
"size"];
2814 return_slice = slice->process(
"filter.linearfourier");
2829 return return_slice;
2835 LOGERR(
"try to insert NULL slice");
2839 if (input->get_xsize() != input->get_ysize() || input->get_xsize() !=
nx) {
2840 LOGERR(
"tried to insert image that was not correction dimensions");
2857 tmp->set_size(
nx,
ny,
nz);
2859 float *slice_data = slice->get_data();
2862 size_t nxy =
nx *
ny;
2863 size_t nxy_size = nxy *
sizeof(float);;
2864 for (
int i = 0; i <
nz; ++i) {
2865 memcpy(&
tmp_data[nxy * i], slice_data, nxy_size);
2889 input_slice->set_attr(
"reconstruct_norm",1.0f);
2890 input_slice->set_attr(
"reconstruct_absqual",1.0f);
2891 input_slice->set_attr(
"reconstruct_weight",1.0f);
2901 vector<Transform> syms = sym->
get_syms();
2903 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
2908 image->add(*tmpcopy);
2915 image->process_inplace(
"mask.sharp",
Dict(
"outer_radius",
image->get_xsize()/2-1));
2917 else printf(
"No masking %d %d %d\n",
image->get_xsize(),
image->get_ysize(),
image->get_zsize());
2926 int nx = slice->get_xsize();
2927 int ny = slice->get_ysize();
2928 int padffted= slice->get_attr_default(
"padffted", 0);
2929 int ndim = (ny==1) ? 1 : 2;
2930 int iext = slice->is_fftodd();
2931 int extension = (2-iext)*padffted;
2933 if( ndim==2 && (nx-extension)!=ny ) {
2934 LOGERR(
"Input image must be square!");
2938 EMData* padfftslice = NULL;
2939 if( padffted == 0) {
2941 EMData* temp = slice->average_circ_sub();
2943 padfftslice = temp->norm_pad(
false, npad );
2946 padfftslice->do_fft_inplace();
2948 padfftslice =
new EMData(*slice);
2953 float sx = -trans[0];
2954 float sy = -trans[1];
2955 if(sx != 0.0f || sy != 0.0) padfftslice->process_inplace(
"filter.shift",
Dict(
"x_shift", sx,
"y_shift", sy,
"z_shift", 0.0f));
2957 int remove = slice->get_attr_default(
"remove", 0);
2958 padfftslice->set_attr(
"remove", remove );
2959 padfftslice->center_origin_fft();
2969float max2d(
int kc,
const vector<float>& pow_a )
2972 for(
int i=-kc; i <= kc; ++i ) {
2973 for(
int j=-kc; j <= kc; ++j ) {
2974 if( i==0 && j==0 )
continue;
2976 int c = 2*kc+1 - std::abs(i) - std::abs(j);
2977 max = max + pow_a[c];
2984float max3d(
int kc,
const vector<float>& pow_a )
2987 for(
int i=-kc; i <= kc; ++i ) {
2988 for(
int j=-kc; j <= kc; ++j ) {
2989 for(
int k=-kc; k <= kc; ++k ) {
2990 if( i==0 && j==0 && k==0 )
continue;
2993 int c = 3*kc+1 - std::abs(i) - std::abs(j) - std::abs(k);
2994 max = max + pow_a[c];
3005#define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
3009 float *
tw = win->get_data();
3012 int ix = win->get_xsize();
3013 int iy = win->get_ysize();
3014 int iz = win->get_zsize();
3015 int L2 = (ix/2)*(ix/2);
3016 int L2P = (ix/2-1)*(ix/2-1);
3023 float* sincx =
new float[IP+1];
3024 float* sincy =
new float[JP+1];
3025 float* sincz =
new float[KP+1];
3032 if( npad == 1 ) cor = 1.0;
3035 float cdf = M_PI/(cor*ix);
3036 for (
int i = 1; i <= IP; ++i) sincx[i] = sin(i*cdf)/(i*cdf);
3037 cdf = M_PI/(cor*iy);
3038 for (
int i = 1; i <= JP; ++i) sincy[i] = sin(i*cdf)/(i*cdf);
3039 cdf = M_PI/(cor*iz);
3040 for (
int i = 1; i <= KP; ++i) sincz[i] = sin(i*cdf)/(i*cdf);
3041 for (
int k = 1; k <= iz; ++k) {
3042 int kkp = abs(k-KP);
3043 for (
int j = 1; j <= iy; ++j) {
3044 cdf = sincy[abs(j- JP)]*sincz[kkp];
3045 for (
int i = 1; i <= ix; ++i)
tw(i,j,k) /= (sincx[abs(i-IP)]*cdf);
3055 for (
int k = 1; k <= iz; ++k) {
3056 for (
int j = 1; j <= iy; ++j) {
3057 for (
int i = 1; i <= ix; ++i) {
3058 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
3059 if (LR >= (
size_t)L2P && LR<=(size_t)L2) {
3070 for (
int k = 1; k <= iz; ++k) {
3071 for (
int j = 1; j <= iy; ++j) {
3072 for (
int i = 1; i <= ix; ++i) {
3073 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
3074 if (LR<=(
size_t)L2)
tw(i,j,k) -= TNR;
3075 else tw(i,j,k) = 0.0f;
3085 float *
tw = win->get_data();
3088 int ix = win->get_xsize();
3089 int iy = win->get_ysize();
3090 int iz = win->get_zsize();
3091 int L2 = (ix/2)*(ix/2);
3092 int L2P = (ix/2-1)*(ix/2-1);
3099 float* sincx =
new float[IP+1];
3100 float* sincy =
new float[JP+1];
3101 float* sincz =
new float[KP+1];
3108 if( npad == 1 ) cor = 1.0;
3111 float cdf = M_PI/(cor*ix);
3112 for (
int i = 1; i <= IP; ++i) sincx[i] = pow(sin(i*cdf)/(i*cdf),2);
3113 cdf = M_PI/(cor*iy);
3114 for (
int i = 1; i <= JP; ++i) sincy[i] = pow(sin(i*cdf)/(i*cdf),2);
3115 cdf = M_PI/(cor*iz);
3116 for (
int i = 1; i <= KP; ++i) sincz[i] = pow(sin(i*cdf)/(i*cdf),2);
3117 for (
int k = 1; k <= iz; ++k) {
3118 int kkp = abs(k-KP);
3119 for (
int j = 1; j <= iy; ++j) {
3120 cdf = sincy[abs(j- JP)]*sincz[kkp];
3121 for (
int i = 1; i <= ix; ++i)
tw(i,j,k) /= (sincx[abs(i-IP)]*cdf);
3131 for (
int k = 1; k <= iz; ++k) {
3132 for (
int j = 1; j <= iy; ++j) {
3133 for (
int i = 1; i <= ix; ++i) {
3134 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
3135 if (LR >= (
size_t)L2P && LR<=(size_t)L2) {
3146 for (
int k = 1; k <= iz; ++k) {
3147 for (
int j = 1; j <= iy; ++j) {
3148 for (
int i = 1; i <= ix; ++i) {
3149 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
3150 if (LR<=(
size_t)L2)
tw(i,j,k) -= TNR;
3151 else tw(i,j,k) = 0.0f;
3172 setup( symmetry, size, npad );
3190 int size =
params[
"size"];
3191 int npad =
params[
"npad"];
3195 else symmetry =
"c1";
3203 setup( symmetry, size, npad );
3234 int offset = 2 -
m_vnxp%2;
3255 m_volume->set_array_offsets(0,1,1);
3268 m_wptr->set_array_offsets(0,1,1);
3273 Assert( line->get_zsize()==1 );
3276 int nx = line->get_xsize();
3277 int ny = line->get_ysize();
3278 for(
int j=0; j < ny; ++j ) {
3279 for(
int i=0; i < nx; ++i ) printf(
"%10.3f ", line->get_value_at(i,j) );
3289 LOGERR(
"try to insert NULL slice");
3293 int padffted= slice->get_attr_default(
"padffted", 0 );
3295 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=
m_vnx) ) {
3297 LOGERR(
"Tried to insert a slice that is the wrong size.");
3302 if( slice->get_ysize() !=1 ) {
3303 LOGERR(
"for 2D reconstruction, a line is excepted" );
3307 if( weight > 0.0f ) {
3314 float alpha = padfft->get_attr(
"alpha" );
3315 alpha = alpha/180.0f*M_PI;
3316 for(
int i=0; i <
m_vnxc+1; ++i ) {
3317 float xnew = i*cos(alpha);
3318 float ynew = -i*sin(alpha);
3319 float btqr = padfft->get_value_at( 2*i, 0, 0 );
3320 float btqi = padfft->get_value_at( 2*i+1, 0, 0 );
3330 if(iyn < 0 ) iyn +=
m_vnyp;
3332 (*m_volume)( 2*ixn, iyn+1, 1 ) += btqr * weight;
3333 (*m_volume)( 2*ixn+1, iyn+1, 1 ) += btqi * weight;
3334 (*m_wptr)(ixn,iyn+1, 1) += weight;
3345 Assert( padfft != NULL );
3348 for (
unsigned int isym=0; isym < tsym.size(); isym++)
m_volume->nn(
m_wptr, padfft, tsym[isym], weight);
3359 for(
int i=1; i <=
m_vnyp; ++i ) {
3361 if( (*
m_wptr)(0, i, 1)==0.0 ) {
3363 (*m_wptr)(0, i, 1) = (*
m_wptr)(0, j, 1);
3364 (*m_volume)(0, i, 1) = (*
m_volume)(0, j, 1);
3365 (*m_volume)(1, i, 1) = (*
m_volume)(1, j, 1);
3373 vector< float > pow_a(
m_ndim*kc+1, 1.0 );
3374 for(
unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(
m_wghta);
3379 int vol = box*box*box;
3380 float max =
max3d( kc, pow_a );
3381 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
3384 float max =
max2d( kc, pow_a );
3385 alpha = ( 1.0f - 1.0f/(float)ara ) / max;
3388 for (iz = 1; iz <=
m_vnzp; iz++) {
3389 for (iy = 1; iy <=
m_vnyp; iy++) {
3390 for (ix = 0; ix <=
m_vnxc; ix++) {
3391 if ( (*
m_wptr)(ix,iy,iz) > 0) {
3392 float tmp = (-2*((ix+iy+iz)%2)+1)/((*
m_wptr)(ix,iy,iz)+
m_osnr);
3398 for(
int ii = -kc; ii <= kc; ++ii ) {
3399 int nbrcx = cx + ii;
3400 if( nbrcx >=
m_vnxc )
continue;
3401 for(
int jj= -kc; jj <= kc; ++jj ) {
3402 int nbrcy = cy + jj;
3403 if( nbrcy <= -m_vnyc || nbrcy >=
m_vnyc )
continue;
3405 int kcz = (
m_ndim==3) ? kc : 0;
3406 for(
int kk = -kcz; kk <= kcz; ++kk ) {
3407 int nbrcz = cz + kk;
3408 if( nbrcz <= -m_vnyc || nbrcz >=
m_vnyc )
continue;
3415 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 +
m_vnyp;
3416 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 +
m_vnzp;
3417 if( (*
m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
3418 int c =
m_ndim*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
3419 sum = sum + pow_a[c];
3424 float wght = 1.0f / ( 1.0f - alpha * sum );
3428 (*m_volume)(2*ix,iy,iz) *= tmp;
3429 (*m_volume)(2*ix+1,iy,iz) *= tmp;
3441 int npad =
m_volume->get_attr(
"npad");
3444 m_volume->set_array_offsets( 0, 0, 0 );
3461 setup( symmetry, size, npad );
3479 int npad =
params[
"npad"];
3484 else symmetry =
"c1";
3506 float temp=
params[
"xratio"];
3507 m_vnx=int(
float(sizeprojection)*temp);
3509 else m_vnx=sizeprojection;
3514 float temp=
params[
"yratio"];
3515 m_vny=int(
float(sizeprojection)*temp);
3517 else m_vny=sizeprojection;
3524 float temp=
params[
"zratio"];
3525 m_vnz=int(
float(sizeprojection)*temp);
3548 int offset = 2 -
m_vnxp%2;
3569 m_volume->set_array_offsets(0,1,1);
3583 m_wptr->set_array_offsets(0,1,1);
3591 LOGERR(
"try to insert NULL slice");
3595 int padffted= slice->get_attr_default(
"padffted", 0 );
3597 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=
m_sizeofprojection) ) {
3599 LOGERR(
"Tried to insert a slice that is the wrong size.");
3604 if( slice->get_ysize() !=1 ) {
3605 LOGERR(
"for 2D reconstruction, a line is excepted" );
3609 if( weight > 0.0f ) {
3618 float ellipse_length,ellipse_step,cos_alpha,sin_alpha;
3619 int ellipse_length_int;
3620 float alpha = padfft->get_attr(
"alpha" );
3621 alpha = alpha/180.0f*M_PI;
3626 ellipse_length=
sqrt(temp1*temp1+temp2*temp2);
3627 ellipse_length_int=int(ellipse_length);
3629 loop_range=ellipse_length_int;
3630 cos_alpha=temp1/ellipse_length;
3631 sin_alpha=temp2/ellipse_length;
3633 std::cout<<
"#############################################################"<<std::endl;
3634 std::cout<<
"line insert start=="<<
m_count<<std::endl;
3635 std::cout<<
"ellipse length=="<<ellipse_length_int<<
"ellips step=="<<ellipse_step<<std::endl;
3636 std::cout<<
"loop_range"<<loop_range<<std::endl;
3638 std::cout<<
"cos sin of alpha=="<<cos(alpha)<<
" "<<sin(alpha)<<std::endl;
3639 std::cout<<
"cos sin of alpha_new==="<<cos_alpha<<sin_alpha<<std::endl;
3640 std::cout<<
"alpah dig==="<<cos_alpha<<sin_alpha<<std::endl;
3641 std::cout<<
"prjection maximum==="<<loop_range*ellipse_step<<
"ideal maximum"<<
m_sizeofprojection*
m_npad/2<<std::endl;
3642 std::cout<<
"x_size=="<<
m_volume->get_xsize()<<
"y_size=="<<
m_volume->get_ysize()<<std::endl;
3643 std::cout<<
"#############################################################"<<std::endl;
3645 for(
int i=0; i <=loop_range; ++i ) {
3646 float xnew = i*cos_alpha;
3647 float ynew = -i*sin_alpha;
3648 if(
m_count%100==0&&i==loop_range)
3649 std::cout<<
"x_new=="<<xnew<<
"Y_new=="<<ynew<<std::endl;
3650 float btqr=0,btqi=0;
3651 float xprj=i*ellipse_step;
3652 float t=xprj-int(xprj);
3653 btqr = (1-t)*padfft->get_value_at( 2*
int(xprj), 0, 0 )+t*padfft->get_value_at( 2*(1+
int(xprj)), 0, 0 );
3654 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 );
3664 if(iyn < 0 ) iyn +=
m_vnyp;
3665 if(
m_count%100==0&&i==loop_range)
3666 std::cout<<
"xnn=="<<ixn<<
"ynn=="<<iyn<<std::endl;
3667 (*m_volume)( 2*ixn, iyn+1, 1 ) += btqr * weight;
3668 (*m_volume)( 2*ixn+1, iyn+1, 1 ) += btqi * weight;
3669 (*m_wptr)(ixn,iyn+1, 1) += weight;
3684 Assert( padded != NULL );
3687 for (
unsigned int isym=0; isym < tsym.size(); isym++)
3694#define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
3697 float *
tw = win->get_data();
3700 int ix = win->get_xsize();
3701 int iy = win->get_ysize();
3702 int iz = win->get_zsize();
3709 float* sincx =
new float[IP+1];
3710 float* sincy =
new float[JP+1];
3711 float* sincz =
new float[KP+1];
3717 float cdf = M_PI/float(npad*2*ix);
3718 for (
int i = 1; i <= IP; ++i) sincx[i] = sin(i*cdf)/(i*cdf);
3719 cdf = M_PI/float(npad*2*iy);
3720 for (
int i = 1; i <= JP; ++i) sincy[i] = sin(i*cdf)/(i*cdf);
3721 cdf = M_PI/float(npad*2*iz);
3722 for (
int i = 1; i <= KP; ++i) sincz[i] = sin(i*cdf)/(i*cdf);
3723 for (
int k = 1; k <= iz; ++k) {
3724 int kkp = abs(k-KP);
3725 for (
int j = 1; j <= iy; ++j) {
3726 cdf = sincy[abs(j- JP)]*sincz[kkp];
3727 for (
int i = 1; i <= ix; ++i)
tw(i,j,k) /= (sincx[abs(i-IP)]*cdf);
3737 float dxx = 1.0f/float(0.25*ix*ix);
3738 float dyy = 1.0f/float(0.25*iy*iy);
3742 float LR2=(float(ix)/2-1)*(
float(ix)/2-1)*dxx;
3746 for (
int k = 1; k <= iz; ++k) {
3747 for (
int j = 1; j <= iy; ++j) {
3748 for (
int i = 1; i <= ix; ++i) {
3749 float LR = (j-JP)*(j-JP)*dyy+(i-IP)*(i-IP)*dxx;
3750 if (LR<=1.0f && LR >= LR2) {
3761 for (
int k = 1; k <= iz; ++k) {
3762 for (
int j = 1; j <= iy; ++j) {
3763 for (
int i = 1; i <= ix; ++i) {
3764 float LR = (j-JP)*(j-JP)*dyy+(i-IP)*(i-IP)*dxx;
3765 if (LR<=1.0f)
tw(i,j,k)=
tw(i,j,k)-TNR;
3766 else tw(i,j,k) = 0.0f;
3780 for(
int i=1; i <=
m_vnyp; ++i ) {
3782 if( (*
m_wptr)(0, i, 1)==0.0 ) {
3784 (*m_wptr)(0, i, 1) = (*
m_wptr)(0, j, 1);
3785 (*m_volume)(0, i, 1) = (*
m_volume)(0, j, 1);
3786 (*m_volume)(1, i, 1) = (*
m_volume)(1, j, 1);
3794 vector< float > pow_a(
m_ndim*kc+1, 1.0 );
3795 for(
unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(
m_wghta);
3800 int vol = box*box*box;
3801 float max =
max3d( kc, pow_a );
3802 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
3805 float max =
max2d( kc, pow_a );
3806 alpha = ( 1.0f - 1.0f/(float)ara ) / max;
3810 for (iz = 1; iz <=
m_vnzp; iz++) {
3811 for (iy = 1; iy <=
m_vnyp; iy++) {
3812 for (ix = 0; ix <=
m_vnxc; ix++) {
3813 if ( (*
m_wptr)(ix,iy,iz) > 0) {
3815 tmp = (-2*((ix+iy+iz)%2)+1)/((*
m_wptr)(ix,iy,iz)+
m_osnr);
3822 for(
int ii = -kc; ii <= kc; ++ii ) {
3823 int nbrcx = cx + ii;
3824 if( nbrcx >=
m_vnxc )
continue;
3825 for(
int jj= -kc; jj <= kc; ++jj ) {
3826 int nbrcy = cy + jj;
3827 if( nbrcy <= -m_vnyc || nbrcy >=
m_vnyc )
continue;
3829 int kcz = (
m_ndim==3) ? kc : 0;
3830 for(
int kk = -kcz; kk <= kcz; ++kk ) {
3831 int nbrcz = cz + kk;
3832 if( nbrcz <= -m_vnyc || nbrcz >=
m_vnyc )
continue;
3839 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 +
m_vnyp;
3840 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 +
m_vnzp;
3841 if( (*
m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
3842 int c =
m_ndim*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
3843 sum = sum + pow_a[c];
3848 float wght = 1.0f / ( 1.0f - alpha * sum );
3851 (*m_volume)(2*ix,iy,iz) *= tmp;
3852 (*m_volume)(2*ix+1,iy,iz) *= tmp;
3864 int npad =
m_volume->get_attr(
"npad");
3867 m_volume->set_array_offsets( 0, 0, 0 );
3890 setup( symmetry, size, npad );
3906 int size =
params[
"size"];
3907 int npad =
params[
"npad"];
3911 else symmetry =
"c1";
3913 setup( symmetry, size, npad );
3958 m_volume->set_array_offsets(0,1,1);
3966 m_wptr->set_array_offsets(0,1,1);
3974 m_wptr2->set_array_offsets(0,1,1);
3981 LOGERR(
"try to insert NULL slice");
3984 if( weight > 0.0f ) {
3985 int padffted=slice->get_attr_default(
"padffted", 0 );
3987 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=
m_vnx) ) {
3989 LOGERR(
"Tried to insert a slice that has wrong size.");
4004 Assert( padfft != NULL );
4026 float argx, argy, argz;
4027 vector< float > pow_a( 3*kc+1, 1.0 );
4036 vol_ssnr->to_zero();
4037 if (
m_vnxp % 2 == 0 ) vol_ssnr->set_fftodd(0);
4038 else vol_ssnr->set_fftodd(1);
4039 vol_ssnr->set_nxc(
m_vnxc);
4040 vol_ssnr->set_complex(
true);
4041 vol_ssnr->set_ri(
true);
4042 vol_ssnr->set_fftpad(
false);
4050 SSNR->set_size(inc+1,4,1);
4052 float *nom =
new float[inc+1];
4053 float *denom =
new float[inc+1];
4054 int *nn =
new int[inc+1];
4055 int *ka =
new int[inc+1];
4057 for (
int i = 0; i <= inc; i++) {
4067 int vol = box*box*box;
4068 for(
unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(
m_wghta);
4070 float max =
max3d( kc, pow_a );
4071 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
4074 for (
int iz = 1; iz <=
m_vnzp; iz++) {
4076 argz = float(kz*kz)*dz2;
4077 for (
int iy = 1; iy <=
m_vnyp; iy++) {
4079 argy = argz + float(ky*ky)*dy2;
4080 for (
int ix = 0; ix <=
m_vnxc; ix++) {
4081 float Kn = (*m_wptr)(ix,iy,iz);
4082 argx =
std::sqrt(argy +
float(ix*ix)*dx2);
4084 if ( r >= 0 && Kn > 4.5f ) {
4091 for(
int ii = -kc; ii <= kc; ++ii ) {
4092 int nbrcx = cx + ii;
4093 if( nbrcx >=
m_vnxc )
continue;
4094 for (
int jj= -kc; jj <= kc; ++jj ) {
4095 int nbrcy = cy + jj;
4096 if( nbrcy <= -m_vnyc || nbrcy >=
m_vnyc )
continue;
4097 for(
int kk = -kc; kk <= kc; ++kk ) {
4098 int nbrcz = cz + jj;
4099 if ( nbrcz <= -m_vnyc || nbrcz >=
m_vnyc )
continue;
4106 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 +
m_vnyp;
4107 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 +
m_vnzp;
4108 if( (*
m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
4109 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
4110 sum = sum + pow_a[c];
4115 wght = 1.0f / ( 1.0f - alpha * sum );
4117 float nominator = std::norm(
m_volume->cmplx(ix,iy,iz)/Kn);
4118 float denominator = ((*m_wptr2)(ix,iy,iz)-nominator)/(Kn-1.0f);
4120 if( (ix>0 || (kz>=0 && (ky>=0 || kz!=0)))) {
4122 nom[r] += nominator*wght;
4123 denom[r] += denominator/Kn*wght;
4140 (*vol_ssnr)(2*ix, iy-1, iz-1) = nominator*wght;
4148 for (
int i = 0; i <= inc; i++) {
4149 (*SSNR)(i,0,0) = nom[i];
4150 (*SSNR)(i,1,0) = denom[i];
4151 (*SSNR)(i,2,0) =
static_cast<float>(nn[i]);
4152 (*SSNR)(i,3,0) =
static_cast<float>(ka[i]);
4178 setup( symmetry, size, npad, snr, sign );
4192 if( !
params.
has_key(
"size") )
throw std::logic_error(
"Error: image size is not given");
4194 int size =
params[
"size"];
4200 float snr =
params[
"snr"];
4203 setup( symmetry, size, npad, snr, sign );
4243 int offset = 2 -
m_vnxp%2;
4259 m_volume->set_array_offsets(0,1,1);
4271 m_wptr->set_array_offsets(0,1,1);
4279 LOGERR(
"try to insert NULL slice");
4282 if( weight > 0.0f ) {
4283 int buffed = slice->get_attr_default(
"buffed", 0 );
4289 int padffted= slice->get_attr_default(
"padffted", 0);
4290 if( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=
m_vnx) ) {
4292 LOGERR(
"Tried to insert a slice that is the wrong size.");
4298 float tmp = padfft->get_attr_default(
"ctf_applied", 0);
4299 int ctf_applied = (int) tmp;
4304 int winsize = padfft->get_ysize();
4305 Ctf* ctf = padfft->get_attr(
"ctf" );
4310 if(ctf) {
delete ctf; ctf=0;}
4311 int nx=ctf2d->get_xsize(),ny=ctf2d->get_ysize(),nz=ctf2d->get_zsize();
4312 float *ctf2d_ptr = ctf2d->get_data();
4314 size_t size = (size_t)nx*ny*nz;
4316 for (
int i = 0; i < size; ++i) padfft->cmplx(i) *= ctf2d_ptr[i];
4319 for (
int i = 0; i < size; ++i) ctf2d_ptr[i] *= ctf2d_ptr[i];
4332 const float* bufdata = buffed->get_data();
4333 float* cdata =
m_volume->get_data();
4334 float* wdata =
m_wptr->get_data();
4336 int npoint = buffed->get_xsize()/4;
4337 for(
int i=0; i < npoint; ++i ) {
4339 int pos2 = int( bufdata[4*i] );
4340 int pos1 = pos2 * 2;
4341 cdata[pos1 ] += bufdata[4*i+1]*weight;
4342 cdata[pos1+1] += bufdata[4*i+2]*weight;
4343 wdata[pos2 ] += bufdata[4*i+3]*weight;
4355 Assert( padfft != NULL );
4357 vector<float> abc_list;
4358 int abc_list_len = 0;
4361 abc_list =
m_volume->get_attr(
"smear");
4362 abc_list_len = abc_list.size();
4366 for (
unsigned int isym=0; isym < tsym.size(); isym++) {
4367 if (abc_list_len == 0)
4368 m_volume->nn_ctf_exists(
m_wptr, padfft, ctf2d2, tsym[isym], weight);
4370 for (
int i = 0; i < abc_list_len; i += 4) {
4371 m_volume->nn_ctf_exists(
m_wptr, padfft, ctf2d2, tsym[isym] *
Transform(
Dict(
"type",
"SPIDER",
"phi", abc_list[i],
"theta", abc_list[i+1],
"psi", abc_list[i+2])), weight * abc_list[i+3]);
4379 m_volume->set_array_offsets(0, 1, 1);
4380 m_wptr->set_array_offsets(0, 1, 1);
4384 int vol = box*box*box;
4386 vector< float > pow_a( 3*kc+1, 1.0 );
4387 for(
unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(
m_wghta);
4391 float max =
max3d( kc, pow_a );
4392 float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
4393 float osnr = 1.0f/
m_snr;
4397 for (iz = 1; iz <=
m_vnzp; iz++) {
4398 for (iy = 1; iy <=
m_vnyp; iy++) {
4399 for (ix = 0; ix <=
m_vnxc; ix++) {
4400 if ( (*
m_wptr)(ix,iy,iz) > 0.0f) {
4405 float freq =
sqrt( (
float)(ix*ix+iyp*iyp+izp*izp) );
4406 tmp = (-2*((ix+iy+iz)%2)+1)/((*
m_wptr)(ix,iy,iz)+freq*osnr);
4408 tmp = (-2*((ix+iy+iz)%2)+1)/((*
m_wptr)(ix,iy,iz)+osnr);
4416 for(
int ii = -kc; ii <= kc; ++ii ) {
4417 int nbrcx = cx + ii;
4418 if( nbrcx >=
m_vnxc )
continue;
4419 for(
int jj= -kc; jj <= kc; ++jj ) {
4420 int nbrcy = cy + jj;
4421 if( nbrcy <= -m_vnyc || nbrcy >=
m_vnyc )
continue;
4422 for(
int kk = -kc; kk <= kc; ++kk ) {
4423 int nbrcz = cz + jj;
4424 if( nbrcz <= -m_vnyc || nbrcz >=
m_vnyc )
continue;
4432 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 +
m_vnyp;
4433 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 +
m_vnzp;
4434 if( (*
m_wptr)( nbrix, nbriy, nbriz ) == 0.0 ) {
4435 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
4436 sum = sum + pow_a[c];
4443 float wght = 1.0f / ( 1.0f - alpha * sum );
4455 (*m_volume)(2*ix,iy,iz) *= tmp;
4456 (*m_volume)(2*ix+1,iy,iz) *= tmp;
4464 int npad =
m_volume->get_attr(
"npad");
4467 m_volume->set_array_offsets( 0, 0, 0 );
4484 setup( symmetry, size, npad, snr, sign, do_ctf );
4498 if( !
params.
has_key(
"size") )
throw std::logic_error(
"Error: image size is not given");
4500 int size =
params[
"size"];
4506 float snr =
params[
"snr"];
4507 int do_ctf =
params[
"do_ctf"];
4509 setup( symmetry, size, npad, snr, sign, do_ctf );
4553 int offset = 2 -
m_vnxp%2;
4569 m_volume->set_array_offsets(0,1,1);
4581 m_wptr->set_array_offsets(0,1,1);
4589 LOGERR(
"try to insert NULL slice");
4605 float tmp = padfft->get_attr_default(
"ctf_applied", 0);
4606 int ctf_applied = (int) tmp;
4610 int winsize = padfft->get_ysize();
4611 Ctf* ctf = padfft->get_attr(
"ctf" );
4613 ctf2d = Util::ctf_img_real(winsize , winsize, 1,
params[
"defocus"],
params[
"apix"],
params[
"voltage"],
params[
"cs"], \
4616 if(ctf) {
delete ctf; ctf=0;}
4618 int nx=ctf2d->get_xsize(),ny=ctf2d->get_ysize(),nz=ctf2d->get_zsize();
4619 float *ctf2d_ptr = ctf2d->get_data();
4621 size_t size = (size_t)nx*ny*nz;
4623 for (
int i = 0; i < size; ++i) padfft->cmplx(i) *= ctf2d_ptr[i];
4626 for (
int i = 0; i < size; ++i) ctf2d_ptr[i] *= ctf2d_ptr[i];
4628 int nx=padfft->get_xsize(),ny=padfft->get_ysize(),nz=padfft->get_zsize();
4630 ctf2d->set_size(nx/2,ny,nz);
4631 float *ctf2d_ptr = ctf2d->get_data();
4632 size_t size = (size_t)nx*ny*nz/2;
4633 for (
int i = 0; i < size; ++i) ctf2d_ptr[i] = 1.0;
4636 vector<float> bckgnoise;
4637 bckgnoise = slice->get_attr(
"bckgnoise");
4651 Assert( padfft != NULL );
4653 vector<float> abc_list;
4654 int abc_list_len = 0;
4656 abc_list =
m_volume->get_attr(
"smear");
4657 abc_list_len = abc_list.size();
4661 for (
unsigned int isym=0; isym < tsym.size(); isym++) {
4662 if (abc_list_len == 0)
4665 for (
int i = 0; i < abc_list_len; i += 4)
4666 m_volume->nn_ctfw(
m_wptr, padfft, ctf2d2,
m_npad, bckgnoise, tsym[isym] *
Transform(
Dict(
"type",
"SPIDER",
"phi", abc_list[i],
"theta", abc_list[i+1],
"psi", abc_list[i+2])), weight * abc_list[i+3]);
4674 m_volume->set_array_offsets(0, 1, 1);
4675 m_wptr->set_array_offsets(0, 1, 1);
4676 m_refvol->set_array_offsets(0, 1, 1);
4679 bool do_invert =
false;
4680 bool refvol_present = (*m_refvol)(0) > 0.0f ;
4699 if( refvol_present ) {
4700 for (ix = 0; ix <
m_vnyc; ix++) {
4707 vector<float> count(
m_vnyc+1, 0.0f);
4708 vector<float> sigma2(
m_vnyc+1, 0.0f);
4711 for (iz = 1; iz <=
m_vnzp; iz++) {
4713 float argz = float(izp*izp);
4714 for (iy = 1; iy <=
m_vnyp; iy++) {
4716 float argy = argz + float(iyp*iyp);
4717 for (ix = 0; ix <=
m_vnxc; ix++) {
4718 if(ix>0 || (izp>=0 && (iyp>=0 || izp!=0))) {
4719 float r =
std::sqrt(argy +
float(ix*ix));
4721 if (ir <= limitres) {
4722 float frac = r - float(ir);
4723 float qres = 1.0f - frac;
4724 float temp = (*m_wptr)(ix,iy,iz);
4727 sigma2[ir] += temp*qres;
4728 sigma2[ir+1] += temp*frac;
4730 count[ir+1] += frac;
4736 for (ix = 0; ix <=
m_vnyc; ix++) {
4737 if( count[ix] > 0.0f ) (*m_refvol)(ix) = sigma2[ix]/count[ix];
4754 for (iz = 1; iz <=
m_vnzp; iz++) {
4756 float argz = float(izp*izp);
4757 for (iy = 1; iy <=
m_vnyp; iy++) {
4759 float argy = argz + float(iyp*iyp);
4760 for (ix = 0; ix <=
m_vnxc; ix++) {
4761 float r =
std::sqrt(argy +
float(ix*ix));
4763 if (ir <= limitres) {
4764 if ( (*
m_wptr)(ix,iy,iz) > 0.0f) {
4765 if( refvol_present) {
4766 float frac = r - float(ir);
4767 float qres = 1.0f - frac;
4768 osnr = qres*(*m_refvol)(ir) + frac*(*
m_refvol)(ir+1);
4773 float tmp = ((*m_wptr)(ix,iy,iz)+osnr);
4778 (*m_volume)(2*ix,iy,iz) /= tmp;
4779 (*m_volume)(2*ix+1,iy,iz) /= tmp;
4781 (*m_volume)(2*ix,iy,iz) = 0.0f;
4782 (*m_volume)(2*ix+1,iy,iz) = 0.0f;
4784 }
else (*
m_wptr)(ix,iy,iz) = tmp;
4787 (*m_volume)(2*ix,iy,iz) = 0.0f;
4788 (*m_volume)(2*ix+1,iy,iz) = 0.0f;
4798 int npad =
m_volume->get_attr(
"npad");
4801 m_volume->set_array_offsets( 0, 0, 0 );
4812 m_volume->set_array_offsets(0, 1, 1);
4813 m_wptr->set_array_offsets(0, 1, 1);
4817 int vol = box*box*box;
4819 vector< float > pow_a( 3*kc+1, 1.0 );
4820 for(
unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(
m_wghta);
4824 float max =
max3d( kc, pow_a );
4825 float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
4826 float osnr = 1.0f/
m_snr;
4830 for (iz = 1; iz <=
m_vnzp; iz++) {
4831 for (iy = 1; iy <=
m_vnyp; iy++) {
4832 for (ix = 0; ix <=
m_vnxc; ix++) {
4833 if ( (*
m_wptr)(ix,iy,iz) > 0.0f) {
4838 float freq =
sqrt( (
float)(ix*ix+iyp*iyp+izp*izp) );
4839 tmp = (-2*((ix+iy+iz)%2)+1)/((*
m_wptr)(ix,iy,iz)+freq*osnr);
4841 tmp = (-2*((ix+iy+iz)%2)+1)/((*
m_wptr)(ix,iy,iz)+osnr);
4849 for(
int ii = -kc; ii <= kc; ++ii ) {
4850 int nbrcx = cx + ii;
4851 if( nbrcx >=
m_vnxc )
continue;
4852 for(
int jj= -kc; jj <= kc; ++jj ) {
4853 int nbrcy = cy + jj;
4854 if( nbrcy <= -m_vnyc || nbrcy >=
m_vnyc )
continue;
4855 for(
int kk = -kc; kk <= kc; ++kk ) {
4856 int nbrcz = cz + jj;
4857 if( nbrcz <= -m_vnyc || nbrcz >=
m_vnyc )
continue;
4865 int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 +
m_vnyp;
4866 int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 +
m_vnzp;
4867 if( (*
m_wptr)( nbrix, nbriy, nbriz ) == 0.0 ) {
4868 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
4869 sum = sum + pow_a[c];
4876 float wght = 1.0f / ( 1.0f - alpha * sum );
4888 (*m_volume)(2*ix,iy,iz) *= tmp;
4889 (*m_volume)(2*ix+1,iy,iz) *= tmp;
4897 int npad =
m_volume->get_attr(
"npad");
4900 m_volume->set_array_offsets( 0, 0, 0 );
4907 m_volume->set_array_offsets(0, 1, 1);
4908 m_wptr->set_array_offsets(0, 1, 1);
4909 m_volume->symplane0_ctf(m_wptr);
4912 int vol = box*box*box;
4914 vector< float > pow_a( 3*kc+1, 1.0 );
4915 for(
unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
4919 float max =
max3d( kc, pow_a );
4921 float osnr = 1.0f/m_snr;
4923 vector<float> sigma2(m_vnyc+1, 0.0f);
4924 vector<float> count(m_vnyc+1, 0.0f);
4928 for (iz = 1; iz <= m_vnzp; iz++) {
4929 int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
4930 float argz = float(izp*izp);
4931 for (iy = 1; iy <= m_vnyp; iy++) {
4932 int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
4933 float argy = argz + float(iyp*iyp);
4934 for (ix = 0; ix <= m_vnxc; ix++) {
4935 if(ix>0 || (izp>=0 && (iyp>=0 || izp!=0))) {
4936 float r =
std::sqrt(argy +
float(ix*ix));
4939 float frac = r - float(ir);
4940 float qres = 1.0f - frac;
4941 float temp = (*m_wptr)(ix,iy,iz);
4944 sigma2[ir] += temp*qres;
4945 sigma2[ir+1] += temp*frac;
4947 count[ir+1] += frac;
4953 for (ix = 0; ix <= m_vnyc+1; ix++) {
4954 if( sigma2[ix] > 0.0f ) sigma2[ix] = count[ix]/sigma2[ix];
4959 float fudge = m_refvol->get_attr(
"fudge");
4960 for (ix = 0; ix <= m_vnyc+1; ix++)
4962 for (ix = 0; ix <= m_vnyc+1; ix++) count[ix] = count[ix]/(1.0f - count[ix]) * sigma2[ix];
4963 for (ix = 0; ix <= m_vnyc+1; ix++) {
4964 if ( count[ix] >0.0f) count[ix] = fudge/count[ix];