EMAN2
reconstructor.cpp
Go to the documentation of this file.
1/*
2 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
3 * Copyright (c) 2000-2006 Baylor College of Medicine
4 *
5 * This software is issued under a joint BSD/GNU license. You may use the
6 * source code in this file under either license. However, note that the
7 * complete EMAN2 and SPARX software packages have some GPL dependencies,
8 * so you are responsible for compliance with the licenses of these packages
9 * if you opt to use BSD licensing. The warranty disclaimer below holds
10 * in either instance.
11 *
12 * This complete copyright notice must be included in any revised version of the
13 * source code. Additional authorship citations may be added, but existing
14 * author citations must be preserved.
15 *
16 * This program is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation; either version 2 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 *
30 * */
31
32#include "reconstructor.h"
34#include "ctf.h"
35#include "emassert.h"
36#include "symmetry.h"
37#include <cstring>
38#include <fstream>
39#include <iomanip>
40
41#include <gsl/gsl_statistics_double.h>
42#include <gsl/gsl_fit.h>
43
44#ifdef EMAN2_USING_CUDA
46#endif
47
48//#define DEBUG_POINT 1
49
50using namespace EMAN;
51using std::complex;
52
53
54#include <iostream>
55using std::cerr;
56using std::endl;
57using std::cout; // for debug
58
59#include <algorithm>
60// find, for_each
61
62#include <iomanip>
63using std::setprecision;
64
65template < typename T > void checked_delete( T*& x )
66{
67 typedef char type_must_be_complete[ sizeof(T)? 1: -1 ];
68 (void) sizeof(type_must_be_complete);
69 delete x;
70 x = NULL;
71}
72
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";
87
89{
90 force_add<FourierReconstructor>();
91 force_add<FourierIterReconstructor>();
92 force_add<FourierReconstructorSimple2D>();
93// force_add(&BaldwinWoolfordReconstructor::NEW);
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>();
105// force_add<XYZReconstructor>();
106}
107
109{
110public:
111
112 static void init( int winsize, const Ctf* ctf ) {
113 Dict params = ctf->to_dict();
114
115 m_winsize = winsize;
116
117 m_voltage = params["voltage"];
118 m_pixel = params["apix"];
119 m_cs = params["cs"];
120 m_ampcont = params["ampcont"];
121 m_bfactor = params["bfactor"];
122 m_defocus = params["defocus"];
123 m_dza = params["dfdiff"];
124 m_azz = params["dfang"];
127 }
128
129 // static float get_ctf( int r2 ,int i, int j) {
130 // float ak = std::sqrt( r2/float(m_winsize2) )/m_pixel;
131 // if(m_dza == 0.0f) return Util::tf( m_defocus, ak, m_voltage, m_cs, m_ampcont, m_bfactor, 1);
132 // else {
133 // float az = atan2(float(j), float(i));
134 // float dzz = m_defocus - m_dza/2.0f*sin(2*(az+m_azz*M_PI/180.0f));
135 // return Util::tf( dzz, ak, m_voltage, m_cs, m_ampcont, m_bfactor, 1);
136 // }
137 // }
138
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);
141 }
142
143private:
144
146 static float m_cs;
147 static float m_voltage;
148 static float m_pixel;
149 static float m_ampcont;
150 static float m_bfactor;
151 static float m_defocus;
152 static float m_dza;
153 static float m_azz;
154};
155
157
161
162
163void FourierReconstructorSimple2D::setup()
164{
165 nx = params.set_default("nx",0);
166
167 if ( nx < 0 ) throw InvalidValueException(nx, "nx must be positive");
168
169 bool is_fftodd = (nx % 2 == 1);
170
171 ny = nx;
172 nx += 2-is_fftodd;
173
174 image = new EMData();
175 image->set_size(nx, ny);
176 image->set_complex(true);
177 image->set_fftodd(is_fftodd);
178 image->set_ri(true);
179
180 tmp_data = new EMData();
181 tmp_data->set_size(nx/2, nx);
182}
183
184int FourierReconstructorSimple2D::insert_slice(const EMData* const slice, const Transform & euler, const float)
185{
186
187 // Are these exceptions really necessary? (d.woolford)
188 if (!slice) throw NullPointerException("EMData pointer (input image) is NULL");
189
190 if ( slice->get_ndim() != 1 ) throw ImageDimensionException("Image dimension must be 1");
191
192 // I could also test to make sure the image is the right dimensions...
193 if (slice->is_complex()) throw ImageFormatException("The image is complex, expecting real");
194
195 EMData* working_slice = slice->process("xform.phaseorigin.tocorner");
196
197 // Fourier transform the slice
198 working_slice->do_fft_inplace();
199
200 float* rdata = image->get_data();
201 float* norm = tmp_data->get_data();
202 float* dat = working_slice->get_data();
203
204 float g[4];
205 int offset[4];
206 float dt[2];
207 offset[0] = 0; offset[1] = 2; offset[2] = nx; offset[3] = nx+2;
208
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++) {
211
212 float rx = (float) x;
213
214 float xx = rx*cos(alt);
215 float yy = rx*sin(alt);
216 float cc = 1.0;
217
218 if (xx < 0) {
219 xx = -xx;
220 yy = -yy;
221 cc = -1.0;
222 }
223
224 yy += ny / 2;
225
226
227 dt[0] = dat[2*x];
228 dt[1] = cc * dat[2*x+1];
229
230 // PHIL IS INTERESTED FROM HERE DOWN
231 int x0 = (int) floor(xx);
232 int y0 = (int) floor(yy);
233
234 int i = 2*x0 + y0*nx;
235
236 float dx = xx - x0;
237 float dy = yy - y0;
238
239 g[0] = Util::agauss(1, dx, dy, 0, EMConsts::I2G);
240 g[1] = Util::agauss(1, 1 - dx, dy, 0, EMConsts::I2G);
241 g[2] = Util::agauss(1, dx, 1 - dy, 0, EMConsts::I2G);
242 g[3] = Util::agauss(1, 1 - dx, 1 - dy, 0, EMConsts::I2G);
243
244 // At the extreme we can only do some much...
245 if ( x0 == nx-2 ) {
246 int k = i + offset[0];
247 rdata[k] += g[0] * dt[0];
248 rdata[k + 1] += g[0] * dt[1];
249 norm[k/2] += g[0];
250
251 k = i + offset[2];
252 rdata[k] += g[2] * dt[0];
253 rdata[k + 1] += g[2] * dt[1];
254 norm[k/2] += g[2];
255 continue;
256
257 }
258 // capture and accommodate for periodic boundary conditions in the x direction
259 if ( x0 > nx-2 ) {
260 int dif = x0 - (nx-2);
261 x0 -= dif;
262 }
263 // At the extreme we can only do some much...
264 if ( y0 == ny -1 ) {
265 int k = i + offset[0];
266 rdata[k] += g[0] * dt[0];
267 rdata[k + 1] += g[0] * dt[1];
268 norm[k/2] += g[0];
269
270 k = i + offset[1];
271 rdata[k] += g[1] * dt[0];
272 rdata[k + 1] += g[1] * dt[1];
273 norm[k/2] += g[1];
274 continue;
275 }
276 // capture and accommodate for periodic boundary conditions in the y direction
277 if ( y0 > ny-1) {
278 int dif = y0 - (ny-1);
279 y0 -= dif;
280 }
281
282 if (x0 >= nx - 2 || y0 >= ny - 1) continue;
283
284
285
286
287 for (int j = 0; j < 4; j++)
288 {
289 int k = i + offset[j];
290 rdata[k] += g[j] * dt[0];
291 rdata[k + 1] += g[j] * dt[1];
292 norm[k/2] += g[j];
293
294 }
295 }
296
297 return 0;
298
299}
300
302{
304
305 image->process_inplace("xform.fourierorigin.tocorner");
306 image->do_ift_inplace();
307 image->depad();
308 image->process_inplace("xform.phaseorigin.tocenter");
309
310 EMData *ret = image;
311 image = 0;
312 return ret;
313}
314
315void ReconstructorVolumeData::normalize_threed(const bool sqrtnorm,const bool wiener)
316// normalizes the 3-D Fourier volume. Also imposes appropriate complex conjugate relationships
317{
318 float* norm = tmp_data->get_data();
319 float* rdata = image->get_data();
320
321 size_t nnx=tmp_data->get_xsize();
322 size_t nnxy=tmp_data->get_ysize()*nnx;
323
324// printf("%d %d %d %d %d %d\n",subnx,subny,subnz,image->get_xsize(),image->get_ysize(),image->get_zsize());
325
326 // FIXME should throw a sensible error
327 if ( 0 == norm ) throw NullPointerException("The normalization volume was null!");
328 if ( 0 == rdata ) throw NullPointerException("The complex reconstruction volume was null!");
329
330#ifdef DEBUG_POINT
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);
334#endif
335
336 // add_complex_at handles complex conjugate addition, but the normalization image doesn't
337 // take this into account, so we need to sum the appropriate values
338 // only works for whole volumes!
339 if (subx0==0 && subnx==nx && suby0==0 && subny==ny && subz0==0 && subnz==nz) {
340// printf("cc gain correction\n");
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;
344 // x=0 plane
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;
350
351 // This is the x=nx-1 plane
352 idx1+=nnx-1;
353 idx2+=nnx-1;
354 snorm=norm[idx1]+norm[idx2];
355 norm[idx1]=norm[idx2]=snorm;
356 }
357 }
358 // special cases not handled elsewhere
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;
364 }
365// else printf("Subregion, no CC plane correction\n");
366
367 // The math is a bit tricky to explain. Wiener filter is basically SNR/(1+SNR)
368 // In this case, data have already been effectively multiplied by SNR (one image at a time),
369 // so without Wiener filter, we just divide by total SNR. With Wiener filter, we just add
370 // 1.0 to total SNR, and we're done --steve
371 float wfilt=0.0;
372 if (wiener) wfilt=1.0;
373
374 // actual normalization
375 for (size_t i = 0; i < (size_t)subnx * subny * subnz; i += 2) {
376 float d = norm[i/2];
377 if (sqrtnorm) d=sqrt(d);
378 if (d == 0) {
379 rdata[i] = 0;
380 rdata[i + 1] = 0;
381 }
382 else {
383// rdata[i]=1.0/d;
384// rdata[i+1]=0.0; // for debugging only
385 rdata[i] /= d+wfilt;
386 rdata[i + 1] /= d+wfilt;
387 }
388 }
389
390
391 // This task should now be handled by use of add_complex_at, which adds both values when appropriate
392
393 // enforce complex conj, only works on subvolumes if the complete conjugate plane is in the volume
394// if (subx0==0 && subnx>1 && subny==ny && subnz==nz) {
395// for (int z=0; z<=nz/2; z++) {
396// for (int y=1; y<=ny; y++) {
397// if (y==0 && z==0) continue;
398// // x=0
399// size_t i=(size_t)(y%ny)*subnx+(size_t)(z%nz)*subnx*subny;
400// size_t i2=(size_t)(ny-y)*subnx+(size_t)((nz-z)%nz)*subnx*subny;
401// float ar=(rdata[i]+rdata[i2])/2.0f;
402// float ai=(rdata[i+1]-rdata[i2+1])/2.0f;
403// rdata[i]=ar;
404// rdata[i2]=ar;
405// rdata[i+1]=ai;
406// rdata[i2+1]=-ai;
407// }
408// }
409// }
410
411// if (subx0+subnx==nx && subnx>1 && subny==ny && subnz==nz) {
412// for (int z=0; z<=nz/2; z++) {
413// for (int y=1; y<=ny; y++) {
414// if (y==0 && z==0) continue;
415// // x=0
416// size_t i=(size_t)(y%ny)*subnx+(size_t)(z%nz)*subnx*subny+subnx-2;
417// size_t i2=(size_t)(ny-y)*subnx+(size_t)((nz-z)%nz)*subnx*subny+subnx-2;
418// float ar=(rdata[i]+rdata[i2])/2.0f;
419// float ai=(rdata[i+1]-rdata[i2+1])/2.0f;
420// rdata[i]=ar;
421// rdata[i2]=ar;
422// rdata[i+1]=ai;
423// rdata[i2+1]=-ai;
424// }
425// }
426// }
427}
428
431 // default setting behavior - does not override if the parameter is already set
432 params.set_default("verbose",(int)0);
433
434 vector<int> size=params["size"];
435
436 nx = size[0]+2;
437 ny = size[1];
438 nz = size[2];
439 nx2=nx/2-2;
440 ny2=ny/2;
441 nz2=nz/2;
442
444 subx0=suby0=subz0=0;
445
446 if (nx%2==1) throw ImageDimensionException("ERROR: FourierIterReconstructor only supports even box sizes");
447
448 // The fourier volume where the reconstruction lives
449 if (image) delete image;
450 image = new EMData(nx,ny,nz);
451 image->set_complex(true);
452 image->set_fftodd(false);
453 image->set_ri(true);
454 image->to_zero();
455
456 // This is the normalization volume
457 if (tmp_data) delete tmp_data;
458 tmp_data = new EMData(nx/2,ny,nz);
459 tmp_data->to_zero();
460
461 // this is the reference volume used to define the interpolation kernel
462 // regular "setup" doesn't have one, in which case we default to a local gaussian kernel
463 // This should only happen during the first iteration
464 if (ref_vol) { delete ref_vol; ref_vol=0; }
465// ref_vol = new EMData(nx,ny,nz);
466// ref_vol->set_complex(true);
467// ref_vol->set_fftodd(false);
468// ref_vol->set_ri(true);
469// ref_vol->to_one();
470
471
472
473 if ( (bool) params["verbose"] )
474 {
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);
478 }
479
480}
481void FourierIterReconstructor::setup_seed(EMData* seed,float seed_weight) {
482 setup(); // ok, a little dumb to initialize ref_vol then immediately delete it
483 if (seed->is_complex()) ref_vol=seed->copy(); // we assume tocorner was called properly and dimensions are correct
484 else {
485 ref_vol=seed->do_fft();
486 ref_vol->process_inplace("xform.phaseorigin.tocorner");
487 }
488}
489
491 setup(); // ok, a little dumb to initialize ref_vol then immediately delete it
492 if (seed->is_complex()) ref_vol=seed->copy(); // we assume tocorner was called properly and dimensions are correct
493 else {
494 ref_vol=seed->do_fft();
495 ref_vol->process_inplace("xform.phaseorigin.tocorner");
496 }
497 printf("Warning: FourierIterReconstructor ignores seed weights");
498}
499
500// warning: this clears any existing seed volume
502 setup();
503}
504
506 // Shift the image pixels so the real space origin is now located at the phase origin (at the bottom left of the image)
507 EMData* return_slice = 0;
508 Transform tmp(t);
509 tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly, this way we only do 2d translation,scaling and mirroring
510
511 if (tmp.is_identity()) return_slice=slice->copy();
512 else return_slice = slice->process("xform",Dict("transform",&tmp));
513
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()));
517
518 return_slice->set_attr("reconstruct_preproc",(int)1);
519 return return_slice;
520}
521
522int FourierIterReconstructor::insert_slice(const EMData* const slice, const Transform & arg, const float weight) {
523 Transform * rotation;
524 rotation = new Transform(arg); // assignment operator
525 // We must use only the rotational component of the transform, scaling, translation and mirroring
526 // are not implemented in Fourier space, but are in preprocess_slice
527 rotation->set_scale(1.0);
528 rotation->set_mirror(false);
529 rotation->set_trans(0,0,0);
530
531// if (slice->get_attr_default("reconstruct_preproc",(int) 0)) throw ImageDimensionException("ERROR: FourierIterReconstructor requires preprocess_slice() to be called in advance");
532
533 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
534
535 float inx=(float)(slice->get_xsize()); // x/y dimensions of the input image
536 float iny=(float)(slice->get_ysize());
537 size_t ny2sq=ny*ny/4;
538
539 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
540 Transform t3d = (*rotation)*(*it);
541 for (int y = -iny/2; y < iny/2; y++) {
542 for (int x = 0; x < inx/2; x++) {
543
544 float rx = (float) x/(inx-2.0f); // coords relative to Nyquist=.5
545 float ry = (float) y/iny;
546 int r=Util::hypot_fast_int(x,y);
547 if (r>iny/2+3 && abs(inx-iny)<3) continue; // no filling in Fourier corners...
548
549 Vec3f coord(rx,ry,0);
550 coord = coord*t3d; // transpose multiplication
551 float xx = coord[0]; // transformed coordinates in terms of Nyquist
552 float yy = coord[1];
553 float zz = coord[2];
554
555 // Map back to real pixel coordinates in output volume
556 xx=xx*(nx-2);
557 yy=yy*ny;
558 zz=zz*nz;
559
560
561 int x0 = (int) floor(xx-2.5);
562 int y0 = (int) floor(yy-2.5);
563 int z0 = (int) floor(zz-2.5);
564
565 std::complex<float>dt=slice->get_complex_at(x,y); // The data value we want to insert
566
567 if (x0<-nx2-4 || y0<-ny2-4 || z0<-nz2-4 || x0>nx2+3 || y0>ny2+3 || z0>nz2+3 ) continue;
568
569 // no error checking on add_complex_fast, so we need to be careful here
570 int x1=x0+5;
571 int y1=y0+5;
572 int z1=z0+5;
573 if (x0<-nx2) x0=-nx2;
574 if (x1>nx2) x1=nx2;
575 if (y0<-ny2) y0=-ny2;
576 if (y1>ny2) y1=ny2;
577 if (z0<-nz2) z0=-nz2;
578 if (z1>nz2) z1=nz2;
579
580 if (ref_vol) {
581 // The value in the reference volume at the center insertion location (nearest voxel)
582// std::complex<float>nodt=ref_vol->get_complex_at(Util::round(xx),Util::round(yy),Util::round(zz));
583
584 // Here we trilinear interpolate the point value from the reference volume
585 float xx0=floor(xx);
586 float yy0=floor(yy);
587 float zz0=floor(zz);
588 std::complex<float>nodt=Util::trilinear_interpolate_complex(
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);
598
599// float h=1.0f/EMConsts::I5G;
600 float h=1.0f/2.0f; // This value was optimized empirically for a specific test case
601 float gg=1.0; //default no radial weight
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);
607 //gg=abs(nodt); // we use this as a weight
608 //if (gg==0) continue;
609 //gg=1.0; // turn off local weight
610 //std::complex<float>updt=dt*ntdt/nodt; // This is the unweighted target value for i,j,k
611 std::complex<float>updt=dt+ntdt-nodt;
612
613 float rl = Util::hypot3sq((float) i - xx, j - yy, k - zz);
614 gg = Util::fast_exp(-rl *h);
615
616 size_t off;
617 off=image->add_complex_at_fast(i,j,k,updt*gg*weight);
618 tmp_data->get_data()[off/2]+=gg*weight;
619 }
620 }
621 }
622 } else {
623 float h=1.0f/((Util::hypot3sq(xx,yy,zz)/4000)+.15); // gaussian kernel is used as a weight not a kernel. We increase radius of integration with resolution. Changed away from this on 11/10/17
624 h=h<0.1?0.1:h; // new formula has h from 0.2 to 5
625
626 float rl, gg;
627 for (int k = z0 ; k <= z1; k++) {
628 for (int j = y0 ; j <= y1; j++) {
629 for (int i = x0; i <= x1; i ++) {
630 rl = Util::hypot3sq((float) i - xx, j - yy, k - zz);
631 gg = Util::fast_exp(-rl *h);
632
633 size_t off;
634 off=image->add_complex_at_fast(i,j,k,dt*gg*weight);
635 //if (off==image->get_size())
636 //printf("%d %d %d\t%g %g %g %g\t%ld\n",i,j,k,dt.real(),dt.imag(),gg,weight,off);
637 tmp_data->get_data()[off/2]+=gg*weight; // This is for interpolation rather than gridding
638 }
639 }
640 }
641 }
642 }
643 }
644 }
645
646
647 delete rotation; rotation=0;
648 return 0;
649}
650
651int FourierIterReconstructor::determine_slice_agreement(EMData* input_slice, const Transform & arg, const float weight,bool sub) { return 0; }
652
655
656 if (doift) {
657 image->do_ift_inplace();
658 image->depad();
659 image->process_inplace("xform.phaseorigin.tocenter");
660 }
661
662 image->update();
663
664 if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) {
665 if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
666 tmp_data->write_image((const char *)params["savenorm"]);
667 }
668
669 //Since we give up the ownership of the pointer to-be-returned, it's caller's responsibility to delete the returned image.
670 //So we wrap this function with return_value_policy< manage_new_object >() in libpyReconstructor2.cpp to hand over ownership to Python.
671 EMData *ret=image;
672 image=0;
673
674 return ret;
675}
676
678 if (image) delete image;
679 if (tmp_data) delete tmp_data;
680 if (ref_vol) delete ref_vol;
682}
683
684
687{
688 inserter=0;
689 image=0;
690 tmp_data=0;
691}
692
694{
695 if (image) { delete image; image=0; }
696 if (tmp_data) { delete tmp_data; tmp_data=0; }
697 if ( inserter != 0 )
698 {
699 delete inserter;
700 inserter = 0;
701 }
702}
703
704#include <sstream>
705
707{
708// ss
709// string mode = (string)params["mode"];
710 Dict parms;
711 parms["data"] = image;
712 parms["norm"] = tmp_data->get_data();
713 // These aren't necessary because we deal with them before calling the inserter
714// parms["subnx"] = nx;
715// parms["subny"] = ny;
716// parms["subnx"] = nz;
717// parms["subx0"] = x0;
718// parms["suby0"] = y0;
719// parms["subz0"] = z0;
720
721 if ( inserter != 0 )
722 {
723 delete inserter;
724 }
725
726 inserter = Factory<FourierPixelInserter3D>::get((string)params["mode"], parms);
727 inserter->init();
728}
729
731{
732 // default setting behavior - does not override if the parameter is already set
733 params.set_default("mode","gauss_2");
734 params.set_default("verbose",(int)0);
735
736 vector<int> size=params["size"];
737
738 nx = size[0];
739 ny = size[1];
740 nz = size[2];
741 nx2=nx/2-1;
742 ny2=ny/2;
743 nz2=nz/2;
744
745#ifdef RECONDEBUG
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;
749#endif
750
751 // Adjust nx if for Fourier transform even odd issues
752 bool is_fftodd = (nx % 2 == 1);
753 // The Fourier transform requires one extra pixel in the x direction,
754 // which is two spaces in memory, one each for the complex and the
755 // real components
756 nx += 2-is_fftodd;
757
758 if (params.has_key("subvolume")) {
759 vector<int> sub=params["subvolume"];
760 subx0=sub[0];
761 suby0=sub[1];
762 subz0=sub[2];
763 subnx=sub[3];
764 subny=sub[4];
765 subnz=sub[5];
766
767 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
768 throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
769
770 }
771 else {
772 subx0=suby0=subz0=0;
773 subnx=nx;
774 subny=ny;
775 subnz=nz;
776 }
777
778
779 // Odd dimension support is here atm, but not above.
780 if (image) delete image;
781 image = new EMData();
782 image->set_size(subnx, subny, subnz);
783 image->set_complex(true);
784 image->set_fftodd(is_fftodd);
785 image->set_ri(true);
786// printf("%d %d %d\n\n",subnx,subny,subnz);
787 image->to_zero();
788
789 if (params.has_key("subvolume")) {
790 image->set_attr("subvolume_x0",subx0);
791 image->set_attr("subvolume_y0",suby0);
792 image->set_attr("subvolume_z0",subz0);
793 image->set_attr("subvolume_full_nx",nx);
794 image->set_attr("subvolume_full_ny",ny);
795 image->set_attr("subvolume_full_nz",nz);
796 }
797
798 if (tmp_data) delete tmp_data;
799 tmp_data = new EMData();
800 tmp_data->set_size(subnx/2, subny, subnz);
801 tmp_data->to_zero();
802 tmp_data->update();
803
805
806#ifdef RECONDEBUG
807 printf("copied\n");
808 inserter->ddata=ddata;
809 inserter->dnorm=dnorm;
810#endif
811
812 if ( (bool) params["verbose"] )
813 {
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);
817 }
818}
819
820void FourierReconstructor::setup_seed(EMData* seed,float seed_weight) {
821 // default setting behavior - does not override if the parameter is already set
822 params.set_default("mode","gauss_2");
823
824 vector<int> size=params["size"];
825
826 nx = size[0];
827 ny = size[1];
828 nz = size[2];
829 nx2=nx/2-1;
830 ny2=ny/2;
831 nz2=nz/2;
832
833
834 // Adjust nx if for Fourier transform even odd issues
835 bool is_fftodd = (nx % 2 == 1);
836 // The Fourier transform requires one extra pixel in the x direction,
837 // which is two spaces in memory, one each for the complex and the
838 // real components
839 nx += 2-is_fftodd;
840
841 if (params.has_key("subvolume")) {
842 vector<int> sub=params["subvolume"];
843 subx0=sub[0];
844 suby0=sub[1];
845 subz0=sub[2];
846 subnx=sub[3];
847 subny=sub[4];
848 subnz=sub[5];
849
850 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
851 throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
852
853 }
854 else {
855 subx0=suby0=subz0=0;
856 subnx=nx;
857 subny=ny;
858 subnz=nz;
859 }
860
861 if (seed->get_xsize()!=subnx || seed->get_ysize()!=subny || seed->get_zsize()!=subnz || !seed->is_complex())
862 throw ImageDimensionException("The dimensions of the seed volume do not match the reconstruction size");
863
864 // Odd dimension support is here atm, but not above.
865 image = seed->copy();
866 image->mult(seed_weight); // The lack of this line was what was causing all of the strange normalization issues
867
868
869
870 if (params.has_key("subvolume")) {
871 image->set_attr("subvolume_x0",subx0);
872 image->set_attr("subvolume_y0",suby0);
873 image->set_attr("subvolume_z0",subz0);
874 image->set_attr("subvolume_full_nx",nx);
875 image->set_attr("subvolume_full_ny",ny);
876 image->set_attr("subvolume_full_nz",nz);
877 }
878
879 if (tmp_data) delete tmp_data;
880 tmp_data = new EMData();
881 tmp_data->set_size(subnx/2, subny, subnz);
882 tmp_data->to_value(seed_weight);
883
884 // The values at x=0 and x=nyquist are added twice in the map normally, but the normalization is only added once
885 // This is compensated in the final normalization. In the seed volume we accomodate this by reducing the seed
886 // alternatively we could double the actual values. Impact is similar, if not 100% identical
887 if (subx0==0) {
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);
891 }
892 }
893 tmp_data->set_value_at(0,subny/2,0,seed_weight);
894 tmp_data->set_value_at(0,0,subnz/2,seed_weight);
895 }
896 if (subnx+subx0==nx) {
897 for (int k=0; k<subnz; k++) {
898 for (int j=0; j<subny; j++) {
899 tmp_data->set_value_at(subnx/2-1,j,k,seed_weight/2.0);
900 }
901 }
902 tmp_data->set_value_at(subnx/2-1,subny/2,0,seed_weight);
903 tmp_data->set_value_at(subnx/2-1,0,subnz/2,seed_weight);
904 }
905
906
907
909
910#ifdef DEBUG_POINT
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);
914#endif
915
916
917 if ( (bool) params["quiet"] == false )
918 {
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;
923 }
924}
925
927 // default setting behavior - does not override if the parameter is already set
928
929 // WARNING - when seed_weight is provided it must already be compensated at x=0 and x=nx-1 (should be 1/2 the actual value)
930 params.set_default("mode","gauss_2");
931
932 vector<int> size=params["size"];
933
934 nx = size[0];
935 ny = size[1];
936 nz = size[2];
937 nx2=nx/2-1;
938 ny2=ny/2;
939 nz2=nz/2;
940
941
942 // Adjust nx if for Fourier transform even odd issues
943 bool is_fftodd = (nx % 2 == 1);
944 // The Fourier transform requires one extra pixel in the x direction,
945 // which is two spaces in memory, one each for the complex and the
946 // real components
947 nx += 2-is_fftodd;
948
949 if (params.has_key("subvolume")) {
950 vector<int> sub=params["subvolume"];
951 subx0=sub[0];
952 suby0=sub[1];
953 subz0=sub[2];
954 subnx=sub[3];
955 subny=sub[4];
956 subnz=sub[5];
957
958 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
959 throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
960
961 }
962 else {
963 subx0=suby0=subz0=0;
964 subnx=nx;
965 subny=ny;
966 subnz=nz;
967 }
968
969 if (seed->get_xsize()!=subnx || seed->get_ysize()!=subny || seed->get_zsize()!=subnz || !seed->is_complex())
970 throw ImageDimensionException("The dimensions of the seed volume do not match the reconstruction size");
971
972 // Odd dimension support is here atm, but not above.
973 image = seed->copy(); // note that if the provided seed has not been 'premultiplied' by the weights, bad things may result!
974
975 if (params.has_key("subvolume")) {
976 image->set_attr("subvolume_x0",subx0);
977 image->set_attr("subvolume_y0",suby0);
978 image->set_attr("subvolume_z0",subz0);
979 image->set_attr("subvolume_full_nx",nx);
980 image->set_attr("subvolume_full_ny",ny);
981 image->set_attr("subvolume_full_nz",nz);
982 }
983
984 if (tmp_data) delete tmp_data;
985 tmp_data=seed_weight->copy();
986
988
989 if ( (bool) params["quiet"] == false )
990 {
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;
995 }
996}
997
998
1000{
1001 bool zeroimage = true;
1002 bool zerotmpimg = true;
1003
1004#ifdef EMAN2_USING_CUDA
1005 if(EMData::usecuda == 1) {
1006 if(image->getcudarwdata()) {
1007 to_zero_cuda(image->getcudarwdata(),image->get_xsize(),image->get_ysize(),image->get_zsize());
1008 zeroimage = false;
1009 }
1010 if(tmp_data->getcudarwdata()) {
1011 //to_zero_cuda(tmp_data->getcudarwdata(),image->get_xsize(),image->get_ysize(),image->get_zsize());
1012 zerotmpimg = false;
1013 }
1014 }
1015#endif
1016
1017 if(zeroimage) image->to_zero();
1018 if(zerotmpimg) tmp_data->to_zero();
1019
1020}
1021
1023{
1024#ifdef EMAN2_USING_CUDA
1025 if(EMData::usecuda == 1) {
1026 if(!slice->getcudarwdata()) slice->copy_to_cuda(); //copy slice to cuda using the const version
1027 }
1028#endif
1029 // Shift the image pixels so the real space origin is now located at the phase origin (at the bottom left of the image)
1030 EMData* return_slice = 0;
1031 Transform tmp(t);
1032 tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly, this way we only do 2d translation,scaling and mirroring
1033
1034 if (tmp.is_identity()) return_slice=slice->copy();
1035 else return_slice = slice->process("xform",Dict("transform",&tmp));
1036
1037 return_slice->process_inplace("xform.phaseorigin.tocorner");
1038
1039// printf("@@@ %d\n",(int)return_slice->get_attr("nx"));
1040 // Fourier transform the slice
1041
1042#ifdef EMAN2_USING_CUDA
1043 if(EMData::usecuda == 1 && return_slice->getcudarwdata()) {
1044 return_slice->do_fft_inplace_cuda(); //a CUDA FFT inplace is quite slow as there is a lot of mem copying.
1045 }else{
1046 return_slice->do_fft_inplace();
1047 }
1048#else
1049 return_slice->do_fft_inplace();
1050#endif
1051
1052// printf("%d\n",(int)return_slice->get_attr("nx"));
1053
1054 return_slice->mult((float)sqrt(1.0f/(return_slice->get_ysize())*return_slice->get_xsize()));
1055
1056 // Shift the Fourier transform so that it's origin is in the center (bottom) of the image.
1057// return_slice->process_inplace("xform.fourierorigin.tocenter");
1058
1059 return_slice->set_attr("reconstruct_preproc",(int)1);
1060 return return_slice;
1061}
1062
1063
1064int FourierReconstructor::insert_slice(const EMData* const input_slice, const Transform & arg, const float oweight)
1065{
1066 // Are these exceptions really necessary? (d.woolford)
1067 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
1068
1069
1070// image->set_attr("apix_x",input_slice->get_attr("apix_x"));
1071// image->set_attr("apix_y",input_slice->get_attr("apix_y"));
1072// image->set_attr("apix_z",input_slice->get_attr("apix_z"));
1073
1074#ifdef EMAN2_USING_CUDA
1075 if(EMData::usecuda == 1) {
1076 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda(); //copy slice to cuda using the const version
1077 }
1078#endif
1079
1080 bool usessnr=params.set_default("usessnr",false);
1081 bool corners=params.set_default("corners",false);
1082 float weight=oweight;
1083 if (usessnr) {
1084 if (input_slice->has_attr("class_ssnr")) weight=-1.0; // negative weight is a flag for using SSNR
1085 else weight=0;
1086 }
1087
1088 if (weight==0) return -1;
1089
1090Transform * rotation;
1091/* if ( input_slice->has_attr("xform.projection") ) {
1092 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
1093 } else {*/
1094 rotation = new Transform(arg); // assignment operator
1095// }
1096
1097 EMData *slice;
1098 if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
1099 else slice = preprocess_slice( input_slice, *rotation);
1100
1101
1102 // We must use only the rotational component of the transform, scaling, translation and mirroring
1103 // are not implemented in Fourier space, but are in preprocess_slice
1104 rotation->set_scale(1.0);
1105 rotation->set_mirror(false);
1106 rotation->set_trans(0,0,0);
1107
1108 // Finally to the pixel wise slice insertion
1109 //slice->copy_to_cuda();
1110// EMData *s2=slice->do_ift();
1111// s2->write_image("is.hdf",-1);
1112 do_insert_slice_work(slice, *rotation, weight, corners);
1113
1114 delete rotation; rotation=0;
1115 delete slice;
1116
1117// image->update();
1118 return 0;
1119}
1120
1121// note that negative weight is a prompt for using SSNR from header
1122void FourierReconstructor::do_insert_slice_work(const EMData* const input_slice, const Transform & arg,const float weight,const bool corners)
1123{
1124 // Reload the inserter if the mode has changed
1125// string mode = (string) params["mode"];
1126// if ( mode != inserter->get_name() ) load_inserter();
1127
1128// int y_in = input_slice->get_ysize();
1129// int x_in = input_slice->get_xsize();
1130// // Adjust the dimensions to account for odd and even ffts
1131// if (input_slice->is_fftodd()) x_in -= 1;
1132// else x_in -= 2;
1133
1134 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
1135
1136 float inx=(float)(input_slice->get_xsize()); // x/y dimensions of the input image
1137 float iny=(float)(input_slice->get_ysize());
1138
1139 if (abs(inx-iny)>2 && weight<0) printf("WARNING: Fourier Reconstruction failure. SSNR flag set with asymmetric dimensions on input image\n");
1140
1141#ifdef EMAN2_USING_CUDA
1142 if(EMData::usecuda == 1) {
1143 if(!image->getcudarwdata()){
1144 image->copy_to_cuda();
1145 tmp_data->copy_to_cuda();
1146 }
1147 float * m = new float[12];
1148 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1149 Transform t3d = arg*(*it);
1151 //cout << "using CUDA " << image->getcudarwdata() << endl;
1152 insert_slice_cuda(m,input_slice->getcudarwdata(),image->getcudarwdata(),tmp_data->getcudarwdata(),inx,iny,image->get_xsize(),image->get_ysize(),image->get_zsize(), weight);
1153 }
1154 delete m;
1155 return;
1156 }
1157#endif
1158 vector<float> ssnr;
1159 float sscale = 1.0f;
1160 if (weight<0) {
1161 ssnr=input_slice->get_attr("class_ssnr");
1162 sscale=2.0*(ssnr.size()-1)/iny;
1163 }
1164
1165 float rweight=weight;
1166 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1167 Transform t3d = arg*(*it);
1168 for (int y = -iny/2; y < iny/2; y++) {
1169 for (int x = 0; x < inx/2; x++) {
1170
1171 float rx = (float) x/(inx-2.0f); // coords relative to Nyquist=.5
1172 float ry = (float) y/iny;
1173 if (corners) {
1174 if (weight<0) {
1175 int r=Util::hypot_fast_int(x,y);
1176 rweight=Util::get_max(0.0f,ssnr[int(r*sscale)]);
1177 }
1178 }
1179 else {
1180 int r=Util::hypot_fast_int(x,y);
1181 if (r>iny/2 && abs(inx-iny)<3) continue; // no filling in Fourier corners...
1182
1183 if (weight<0) rweight=Util::get_max(0.0f,ssnr[int(r*sscale)]);
1184 }
1185// printf("%d\t%f\n",int(r*sscale),rweight);
1186
1187 Vec3f coord(rx,ry,0);
1188 coord = coord*t3d; // transpose multiplication
1189 float xx = coord[0]; // transformed coordinates in terms of Nyquist
1190 float yy = coord[1];
1191 float zz = coord[2];
1192
1193 // Map back to real pixel coordinates in output volume
1194 xx=xx*(nx-2);
1195 yy=yy*ny;
1196 zz=zz*nz;
1197
1198// if (x==10 && y==0) printf("10,0 -> %1.2f,%1.2f,%1.2f\t(%5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f) %1.0f %d\n",
1199// xx,yy,zz,t3d.at(0,0),t3d.at(0,1),t3d.at(0,2),t3d.at(1,0),t3d.at(1,1),t3d.at(1,2),t3d.at(2,0),t3d.at(2,1),t3d.at(2,2),inx,nx);
1200// if (x==0 && y==10 FourierReconstructor:) printf("0,10 -> %1.2f,%1.2f,%1.2f\t(%5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f)\n",
1201// xx,yy,zz,t3d.at(0,0),t3d.at(0,1),t3d.at(0,2),t3d.at(1,0),t3d.at(1,1),t3d.at(1,2),t3d.at(2,0),t3d.at(2,1),t3d.at(2,2));
1202
1203 //printf("%3.1f %3.1f %3.1f\t %1.4f %1.4f\t%1.4f\n",xx,yy,zz,input_slice->get_complex_at(x,y).real(),input_slice->get_complex_at(x,y).imag(),weight);
1204// if (floor(xx)==45 && floor(yy)==45 &&floor(zz)==0) printf("%d. 45 45 0\t %d %d\t %1.4f %1.4f\t%1.4f\n",(int)input_slice->get_attr("n"),x,y,input_slice->get_complex_at(x,y).real(),input_slice->get_complex_at(x,y).imag(),weight);
1205// if (floor(xx)==21 && floor(yy)==21 &&floor(zz)==0) printf("%d. 21 21 0\t %d %d\t %1.4f %1.4f\t%1.4f\n",(int)input_slice->get_attr("n"),x,y,input_slice->get_complex_at(x,y).real(),input_slice->get_complex_at(x,y).imag(),weight);
1206 inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),rweight);
1207 }
1208 }
1209 }
1210}
1211
1212int FourierReconstructor::determine_slice_agreement(EMData* input_slice, const Transform & arg, const float weight,bool sub)
1213{
1214 // Are these exceptions really necessary? (d.woolford)
1215 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
1216
1217#ifdef EMAN2_USING_CUDA
1218 if(EMData::usecuda == 1) {
1219 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda(); //copy slice to cuda using the const version
1220 }
1221#endif
1222
1223 Transform * rotation;
1224 rotation = new Transform(arg); // assignment operator
1225
1226 EMData *slice;
1227 if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
1228 else slice = preprocess_slice( input_slice, *rotation);
1229
1230
1231 // We must use only the rotational component of the transform, scaling, translation and mirroring
1232 // are not implemented in Fourier space, but are in preprocess_slice
1233 rotation->set_scale(1.0);
1234 rotation->set_mirror(false);
1235 rotation->set_trans(0,0,0);
1236 if (sub) do_insert_slice_work(slice, *rotation, -weight);
1237 // Remove the current slice first (not threadsafe, but otherwise performance would be awful)
1238
1239 // Compare
1240 do_compare_slice_work(slice, *rotation,weight);
1241
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"));
1245// input_slice->set_attr("reconstruct_qual",slice->get_attr("reconstruct_qual"));
1246 input_slice->set_attr("reconstruct_weight",slice->get_attr("reconstruct_weight"));
1247
1248 // Now put the slice back
1249 if (sub) do_insert_slice_work(slice, *rotation, weight);
1250
1251 delete rotation;
1252 delete slice;
1253
1254// image->update();
1255 return 0;
1256
1257}
1258
1259void FourierReconstructor::do_compare_slice_work(EMData* input_slice, const Transform & arg,float weight)
1260{
1261
1262 float dt[3]; // This stores the complex and weight from the volume
1263 float dt2[2]; // This stores the local image complex
1264 float *dat = input_slice->get_data();
1265 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
1266
1267 float inx=(float)(input_slice->get_xsize()); // x/y dimensions of the input image
1268 float iny=(float)(input_slice->get_ysize());
1269 float apix=input_slice->get_attr("apix_x");
1270 float rmax=iny*apix/12.0; // radius at 12 A, use as cutoff for low res insertion quality
1271 if (rmax>iny/2.0-1.0) rmax=iny/2.0-1.0; // in case of large A/pix
1272 if (rmax<5.0) rmax=5.0; // can't imagine this will happen
1273
1274 double dotlow=0; // low resolution dot product
1275 double dlpow=0;
1276 double dlpow2=0;
1277
1278 double dot=0; // summed pixel*weight dot product
1279 double vweight=0; // sum of weights
1280 double power=0; // sum of inten*weight from volume
1281 double power2=0; // sum of inten*weight from image
1282// double amp=0,amp2=0; // used for normalization, weighted amplitude averages under specific conditions
1283 bool use_cpu = true;
1284
1285 int nval=0;
1286
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();
1292 tmp_data->copy_to_cuda();
1293 }
1294 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1295 Transform t3d = arg*(*it);
1296 float * m = new float[12];
1298 float4 stats = determine_slice_agreement_cuda(m,input_slice->getcudarwdata(),image->getcudarwdata(),tmp_data->getcudarwdata(),inx,iny,image->get_xsize(),image->get_ysize(),image->get_zsize(), weight);
1299 dot = stats.x;
1300 vweight = stats.y;
1301 power = stats.z;
1302 power2 = stats.w;
1303 //cout << "CUDA stats " << stats.x << " " << stats.y << " " << stats.z << " " << stats.w << endl;
1304 use_cpu = false;
1305 }
1306 }
1307#endif
1308 if(use_cpu) {
1309 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1310 Transform t3d = arg*(*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; // We don't want to use the Fourier origin
1314
1315 float rx = (float) x/(inx-2); // coords relative to Nyquist=.5
1316 float ry = (float) y/iny;
1317
1318// if ((rx * rx + Util::square(ry - max_input_dim "xform.projection"/ 2)) > rl)
1319// continue;
1320
1321 Vec3f coord(rx,ry,0);
1322 coord = coord*t3d; // transpose multiplication
1323 float xx = coord[0]; // transformed coordinates in terms of Nyquist
1324 float yy = coord[1];
1325 float zz = coord[2];
1326
1327
1328 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
1329
1330 // Map back to actual pixel coordinates in output volume
1331 xx=xx*(nx-2);
1332 yy=yy*ny;
1333 zz=zz*nz;
1334
1335
1336 int idx = (int)(x * 2 + inx*(y<0?iny+y:y));
1337 dt2[0] = dat[idx];
1338 dt2[1] = dat[idx+1];
1339
1340 // value returned indirectly in dt
1341 if (!pixel_at(xx,yy,zz,dt) || dt[2]==0) continue;
1342
1343// printf("%f\t%f\t%f\t%f\t%f\n",dt[0],dt[1],dt[2],dt2[0],dt2[1]);
1344 float r=(float)Util::hypot_fast(x,y);
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];
1349 }
1350
1351
1352// if (r<6) continue;
1353 vweight+=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];
1357// printf("%d\t%d\t%f\t%f\t%f\n",x,y,dt[0],dt[1],dt[2]);
1358// if (r>2 && dt[2]>=0.5) {
1359// amp+=hypot(dt[0]*dt[0],dt[1]*dt[1])*dt[2];
1360// amp2+=hypot(dt2[0]*dt2[0],dt2[1]*dt2[1])*dt[2];
1361// }
1362 nval++;
1363 }
1364 }
1365 //cout << dot << " " << vweight << " " << power << " " << power2 << endl;
1366 }
1367 }
1368
1369 dot/=sqrt(power*power2); // normalize the dot product
1370 if (dlpow*dlpow2>0) dotlow/=sqrt(dlpow*dlpow2);
1371// input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)));
1372 input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)));
1373// input_slice->set_attr("reconstruct_norm",(float)(amp2<=0?1.0:amp/amp2));
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))); // here weight is a proxy for SNR
1378 input_slice->set_attr("reconstruct_weight",(float)(vweight/(float)(nval)));
1379// printf("** %g\t%g\t%g\t%g ##\n",dot,vweight,power,power2);
1380 //printf("** %f %f %f ##\n",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)),(float)dot,(float)(dot*weight/((weight-1.0)*dot+1.0)));
1381}
1382
1383EMData* FourierReconstructor::projection(const Transform &euler, int ret_fourier) {
1384
1385 if (subx0!=0 || suby0!=0 || subz0!=0 || subnx!=nx || subny!=ny ||subnz!=nz)
1386 throw ImageDimensionException("ERROR: Reconstructor->projection() does not work with subvolumes");
1387
1388 EMData *ret = new EMData(nx,ny,1);
1389 ret->set_complex(1);
1390 ret->set_ri(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"));
1394
1395 // We must use only the rotational component of the transform, scaling, translation and mirroring
1396 // are not implemented in Fourier space, but are in preprocess_slice
1397 Transform rotation(euler);
1398 rotation.set_scale(1.0);
1399 rotation.set_mirror(false);
1400 rotation.set_trans(0,0,0);
1401
1402 float dt[3]; // This stores the complex and weight from the volume
1403
1404// float apix=input_slice->get_attr("apix_x");
1405// float rmax=iny*apix/12.0; // radius at 12 A, use as cutoff for low res insertion quality
1406
1407 for (int y = -ny/2; y < ny/2; y++) {
1408 for (int x = 0; x <= nx/2; x++) {
1409
1410 float rx = (float) x/(nx-2); // coords relative to Nyquist=.5
1411 float ry = (float) y/ny;
1412
1413
1414 Vec3f coord(rx,ry,0);
1415 coord = coord*rotation; // transpose multiplication
1416 float xx = coord[0]; // transformed coordinates in terms of Nyquist
1417 float yy = coord[1];
1418 float zz = coord[2];
1419
1420
1421 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
1422
1423 // Map back to actual pixel coordinates in output volume
1424 xx=xx*(nx-2);
1425 yy=yy*ny;
1426 zz=zz*nz;
1427
1428 if (!pixel_at(xx,yy,zz,dt) || dt[2]<0.005) continue;
1429// if (x==0) printf("%d %d %g %g %g\n",x,y,dt[0],dt[1],dt[2]);
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])); // division by dt[2] already handled in pixel_at
1432 }
1433 }
1434
1435 Transform translation;
1436 translation.set_trans(euler.get_trans());
1437 ret->process_inplace("xform",Dict("transform",&translation));
1438 if (!ret_fourier) {
1439 ret->process_inplace("xform.phaseorigin.tocenter");
1440 EMData *tmp=ret->do_ift();
1441 delete ret;
1442 ret=tmp;
1443 }
1444 return ret;
1445}
1446
1447
1448bool FourierReconstructor::pixel_at(const float& xx, const float& yy, const float& zz, float *dt)
1449{
1450 int x0 = (int) Util::round(xx);
1451 int y0 = (int) Util::round(yy);
1452 int z0 = (int) Util::round(zz);
1453
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;
1459 dt[2]=norm;
1460// printf("%d %d %d %g %g %g\n",x0,y0,z0,dt[0],dt[1],dt[2]);
1461 return true;
1462
1463 //Trilinear interpolation version caused prominent interpolation artifacts due to high noise levels
1464 // between 1/2 Nyquist and Nyquist. While the nearest neighbor interpolation used above still has
1465 // artifacts, it seems better for this puprose.
1466
1467// int x0 = (int) floor(xx);
1468// int y0 = (int) floor(yy);
1469// int z0 = (int) floor(zz);
1470//
1471// float *rdata=image->get_data();
1472// float *norm=tmp_data->get_data();
1473// float normsum=0,normsum2=0;
1474//
1475// dt[0]=dt[1]=dt[2]=0.0;
1476//
1477// if (nx==subnx) { // normal full reconstruction
1478// if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
1479//
1480//
1481//
1482// int x1=x0+1;
1483// int y1=y0+1;
1484// int z1=z0+1;
1485// if (x0<-nx2) x0=-nx2;
1486// if (x1>nx2) x1=nx2;
1487// if (y0<-ny2) y0=-ny2;
1488// if (y1>ny2) y1=ny2;
1489// if (z0<-nz2) z0=-nz2;
1490// if (z1>nz2) z1=nz2;
1491//
1492// size_t idx=0;
1493// float r, gg;
1494// for (int k = z0 ; k <= z1; k++) {
1495// for (int j = y0 ; j <= y1; j++) {
1496// for (int i = x0; i <= x1; i ++) {
1497// idx=image->get_complex_index_fast(i,j,k);
1498// // r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
1499// // gg = Util::fast_exp(-r / EMConsts::I2G);
1500// gg=(1.0-fabs(i-xx))*(1.0-fabs(j-yy))*(1.0-fabs(k-zz)); // Change from Gaussian to trilinear, 6/7/20
1501//
1502// dt[0]+=gg*rdata[idx];
1503// dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
1504// //dt[2]+=norm[idx/2]*gg;
1505// normsum2+=gg;
1506// normsum+=gg*norm[idx/2];
1507// if (i==0) {
1508// normsum2+=gg;
1509// normsum+=gg*norm[idx/2];
1510// }
1511// }
1512// }
1513// }
1514// if (normsum==0) return false;
1515// dt[0]/=normsum;
1516// dt[1]/=normsum;
1517// dt[2]=normsum/normsum2;
1518// #ifdef DEBUG_POINT
1519// if (z0==23 && y0==0 && x0==113) {
1520// std::complex<float> pv = image->get_complex_at(113,0,23);
1521// float pnv = tmp_data->get_value_at(113,0,23);
1522// printf("pixat: %1.4g\t%1.4g\t%1.4g\t%1.4g\t%1.4g\t%1.4g\n",pv.real()/pnv,pv.imag()/pnv,pnv,dt[0],dt[1],dt[2]);
1523// }
1524// #endif
1525// // printf("%1.2f,%1.2f,%1.2f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\n",xx,yy,zz,dt[0],dt[1],dt[2],rdata[idx],rdata[idx+1]);
1526// return true;
1527// }
1528// else { // for subvolumes, not optimized yet
1529// size_t idx;
1530// float r, gg;
1531// for (int k = z0 ; k <= z0 + 1; k++) {
1532// for (int j = y0 ; j <= y0 + 1; j++) {
1533// for (int i = x0; i <= x0 + 1; i ++) {
1534// r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
1535// idx=image->get_complex_index(i,j,k,subx0,suby0,subz0,nx,ny,nz);
1536// gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
1537//
1538// dt[0]+=gg*rdata[idx];
1539// dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
1540// dt[2]+=norm[idx/2];
1541// normsum+=gg;
1542// }
1543// }
1544// }
1545//
1546// if (normsum==0) return false;
1547// return true;
1548// }
1549}
1550
1551
1553{
1554// float *norm = tmp_data->get_data();
1555// float *rdata = image->get_data();
1556#ifdef EMAN2_USING_CUDA
1557 if(EMData::usecuda == 1 && image->getcudarwdata()){
1558 cout << "copy back from CUDA" << endl;
1559 image->copy_from_device();
1560 tmp_data->copy_from_device();
1561 }
1562#endif
1563
1564 bool sqrtnorm=params.set_default("sqrtnorm",false);
1565 normalize_threed(sqrtnorm);
1566// printf("%f\t%f\t%f\n",tmp_data->get_value_at(67,19,1),image->get_value_at(135,19,1),image->get_value_at(134,19,1));
1567
1568// This compares single precision sum to double precision sum near the origin
1569#ifdef RECONDEBUG
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];
1576 a/=b;
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));
1578 }
1579 }
1580 }
1581#endif
1582
1583// tmp_data->write_image("density.mrc");
1584
1585 // we may as well delete the tmp data now... it saves memory and the calling program might
1586 // need memory after it gets the return volume.
1587 // If this delete didn't happen now, it would happen when the deconstructor was called --david
1588 // no longer a good idea with the new iterative scheme -- steve
1589// if ( tmp_data != 0 )
1590// {
1591// delete tmp_data;
1592// tmp_data = 0;
1593// }
1594
1595/* image->process_inplace("xform.fourierorigin.tocorner");*/
1596
1597 if (doift) {
1598 image->do_ift_inplace();
1599 image->depad();
1600 image->process_inplace("xform.phaseorigin.tocenter");
1601 }
1602 // If the image was padded it should be the original size, as the client would expect
1603 // I blocked the rest, it is almost certainly incorrect PAP 07/31/08
1604 // No, it's not incorrect. You are wrong. You have the meaning of nx mixed up. DSAW 09/23/cd
1605 // This should now be handled in the calling program --steve 11/03/09
1606// bool is_fftodd = (nx % 2 == 1);
1607// if ( (nx-2*(!is_fftodd)) != output_x || ny != output_y || nz != output_z )
1608// {
1609// FloatPoint origin( (nx-output_x)/2, (ny-output_y)/2, (nz-output_z)/2 );
1610// FloatSize region_size( output_x, output_y, output_z);
1611// Region clip_region( origin, region_size );
1612// image->clip_inplace( clip_region );
1613// }
1614
1615 // Should be an "if (verbose)" here or something
1616 //print_stats(quality_scores);
1617
1618 image->update();
1619
1620
1621 if (params.has_key("normout")) {
1622// if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
1623 EMData *normout=(EMData*) params["normout"];
1624 normout->set_data(tmp_data->copy()->get_data());
1625 normout->update();
1626 }
1627
1628 if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) {
1629 if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
1630 tmp_data->write_image((const char *)params["savenorm"]);
1631 }
1632
1633 delete tmp_data;
1634 tmp_data=0;
1635 //Since we give up the ownership of the pointer to-be-returned, it's caller's responsibility to delete the returned image.
1636 //So we wrap this function with return_value_policy< manage_new_object >() in libpyReconstructor2.cpp to hand over ownership to Python.
1637 EMData *ret=image;
1638 image=0;
1639
1640 return ret;
1641}
1642
1643int WienerFourierReconstructor::insert_slice(const EMData* const input_slice, const Transform & arg, const float weight)
1644{
1645 // Are these exceptions really necessary? (d.woolford)
1646 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
1647
1648 Transform * rotation;
1649/* if ( input_slice->has_attr("xform.projection") ) {
1650 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
1651 } else {*/
1652 rotation = new Transform(arg); // assignment operator
1653// }
1654
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.");
1657
1658 EMData *slice;
1659 if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
1660 else slice = preprocess_slice( input_slice, *rotation);
1661
1662
1663 // We must use only the rotational component of the transform, scaling, translation and mirroring
1664 // are not implemented in Fourier space, but are in preprocess_slice
1665 rotation->set_scale(1.0);
1666 rotation->set_mirror(false);
1667 rotation->set_trans(0,0,0);
1668
1669 // Finally to the pixel wise slice insertion
1670 do_insert_slice_work(slice, *rotation, weight);
1671
1672 delete rotation; rotation=0;
1673 delete slice;
1674
1675// image->update();
1676 return 0;
1677}
1678
1679void WienerFourierReconstructor::do_insert_slice_work(const EMData* const input_slice, const Transform & arg,const float inweight)
1680{
1681
1682 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
1683
1684 float inx=(float)(input_slice->get_xsize()); // x/y dimensions of the input image
1685 float iny=(float)(input_slice->get_ysize());
1686
1687 int undo_wiener=(int)input_slice->get_attr_default("ctf_wiener_filtered",0); // indicates whether we need to undo a wiener filter before insertion
1688// if (undo_wiener) throw UnexpectedBehaviorException("wiener_fourier does not yet accept already Wiener filtered class-averages. Suggest using ctf.auto averager for now.");
1689
1690 vector<float> snr=input_slice->get_attr("ctf_snr_total");
1691 float sub=1.0;
1692 if (inweight<0) sub=-1.0;
1693 float weight;
1694
1695 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1696 Transform t3d = arg*(*it);
1697 for (int y = -iny/2; y < iny/2; y++) {
1698 for (int x = 0; x <= inx/2; x++) {
1699
1700 float rx = (float) x/(inx-2.0f); // coords relative to Nyquist=.5
1701 float ry = (float) y/iny;
1702
1703 // This deals with the SNR weight
1704 float rn = (float)hypot(rx,ry);
1705 if (rn>=.5) continue; // no SNR in the corners, and we're going to mask them later anyway
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;
1709 else {
1710 rn-=rni;
1711 weight=Util::linear_interpolate(snr[rni],snr[rni+1],rn);
1712 }
1713// if (weight>500.0) printf("%f %d %d %f %f %d %f\n",weight,x,y,rx,ry,rni);
1714
1715 Vec3f coord(rx,ry,0);
1716 coord = coord*t3d; // transpose multiplication
1717 float xx = coord[0]; // transformed coordinates in terms of Nyquist
1718 float yy = coord[1];
1719 float zz = coord[2];
1720
1721 // Map back to real pixel coordinates in output volume
1722 xx=xx*(nx-2);
1723 yy=yy*ny;
1724 zz=zz*nz;
1725
1726// printf("%f\n",weight);
1727 if (undo_wiener) inserter->insert_pixel(xx,yy,zz,(input_slice->get_complex_at(x,y))*((weight+1.0f)/weight),weight*sub);
1728 else inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),weight*sub);
1729 }
1730 }
1731 }
1732}
1733
1734int WienerFourierReconstructor::determine_slice_agreement(EMData* input_slice, const Transform & arg, const float weight,bool sub)
1735{
1736 // Are these exceptions really necessary? (d.woolford)
1737 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
1738
1739 Transform * rotation;
1740 rotation = new Transform(arg); // assignment operator
1741
1742 EMData *slice;
1743 if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
1744 else slice = preprocess_slice( input_slice, *rotation);
1745
1746
1747 // We must use only the rotational component of the transform, scaling, translation and mirroring
1748 // are not implemented in Fourier space, but are in preprocess_slice
1749 rotation->set_scale(1.0);
1750 rotation->set_mirror(false);
1751 rotation->set_trans(0,0,0);
1752
1753// tmp_data->write_image("dbug.hdf",0);
1754
1755 // Remove the current slice first (not threadsafe, but otherwise performance would be awful)
1756 if (sub) do_insert_slice_work(slice, *rotation, -weight);
1757
1758 // Compare
1759 do_compare_slice_work(slice, *rotation,weight);
1760
1761 input_slice->set_attr("reconstruct_norm",slice->get_attr("reconstruct_norm"));
1762 input_slice->set_attr("reconstruct_absqual",slice->get_attr("reconstruct_absqual"));
1763// input_slice->set_attr("reconstruct_qual",slice->get_attr("reconstruct_qual"));
1764 input_slice->set_attr("reconstruct_weight",slice->get_attr("reconstruct_weight"));
1765
1766 // Now put the slice back
1767 if (sub) do_insert_slice_work(slice, *rotation, weight);
1768
1769
1770 delete rotation;
1771 delete slice;
1772
1773// image->update();
1774 return 0;
1775
1776}
1777
1778void WienerFourierReconstructor::do_compare_slice_work(EMData* input_slice, const Transform & arg,float weight)
1779{
1780
1781 float dt[3]; // This stores the complex and weight from the volume
1782 float dt2[2]; // This stores the local image complex
1783 float *dat = input_slice->get_data();
1784 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
1785
1786 float inx=(float)(input_slice->get_xsize()); // x/y dimensions of the input image
1787 float iny=(float)(input_slice->get_ysize());
1788
1789 double dot=0; // summed pixel*weight dot product
1790 double vweight=0; // sum of weights
1791 double power=0; // sum of inten*weight from volume
1792 double power2=0; // sum of inten*weight from image
1793 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1794 Transform t3d = arg*(*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; // We don't want to use the Fourier origin
1798
1799 float rx = (float) x/(inx-2); // coords relative to Nyquist=.5
1800 float ry = (float) y/iny;
1801
1802// if ((rx * rx + Util::square(ry - max_input_dim / 2)) > rl)
1803// continue;
1804
1805 Vec3f coord(rx,ry,0);
1806 coord = coord*t3d; // transpose multiplication
1807 float xx = coord[0]; // transformed coordinates in terms of Nyquist
1808 float yy = coord[1];
1809 float zz = coord[2];
1810
1811
1812 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
1813
1814 // Map back to actual pixel coordinates in output volume
1815 xx=xx*(nx-2);
1816 yy=yy*ny;
1817 zz=zz*nz;
1818
1819
1820 int idx = (int)(x * 2 + inx*(y<0?iny+y:y));
1821 dt2[0] = dat[idx];
1822 dt2[1] = dat[idx+1];
1823
1824 // value returned indirectly in dt
1825 if (!pixel_at(xx,yy,zz,dt) || dt[2]<=0) continue;
1826
1827// printf("%f\t%f\t%f\t%f\t%f\n",dt[0],dt[1],dt[2],dt2[0],dt2[1]);
1828 dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
1829 vweight+=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];
1832 }
1833 }
1834 }
1835
1836 dot/=sqrt(power*power2); // normalize the dot product
1837// input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)));
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))); // here weight is a proxy for SNR
1842 input_slice->set_attr("reconstruct_weight",(float)vweight/(float)(subnx*subny*subnz));
1843// printf("** %g\t%g\t%g\t%g ##\n",dot,vweight,power,power2);
1844 //printf("** %f %f %f ##\n",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)),(float)dot,(float)(dot*weight/((weight-1.0)*dot+1.0)));
1845}
1846
1847bool WienerFourierReconstructor::pixel_at(const float& xx, const float& yy, const float& zz, float *dt)
1848{
1849 int x0 = (int) floor(xx);
1850 int y0 = (int) floor(yy);
1851 int z0 = (int) floor(zz);
1852
1853 float *rdata=image->get_data();
1854 float *norm=tmp_data->get_data();
1855 float normsum=0,normsum2=0;
1856
1857 dt[0]=dt[1]=dt[2]=0.0;
1858
1859 if (nx==subnx) { // normal full reconstruction
1860 if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
1861
1862 // no error checking on add_complex_fast, so we need to be careful here
1863 int x1=x0+1;
1864 int y1=y0+1;
1865 int z1=z0+1;
1866 if (x0<-nx2) x0=-nx2;
1867 if (x1>nx2) x1=nx2;
1868 if (y0<-ny2) y0=-ny2;
1869 if (y1>ny2) y1=ny2;
1870 if (z0<-nz2) z0=-nz2;
1871 if (z1>nz2) z1=nz2;
1872
1873 size_t idx=0;
1874 float r, gg;
1875 for (int k = z0 ; k <= z1; k++) {
1876 for (int j = y0 ; j <= y1; j++) {
1877 for (int i = x0; i <= x1; i ++) {
1878 r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
1879 idx=image->get_complex_index_fast(i,j,k);
1880 gg = Util::fast_exp(-r / EMConsts::I2G);
1881
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;
1885 normsum2+=gg;
1886 normsum+=gg*norm[idx/2];
1887 }
1888 }
1889 }
1890 if (normsum==0) return false;
1891 dt[0]/=normsum;
1892 dt[1]/=normsum;
1893 dt[2]/=normsum2;
1894// printf("%1.2f,%1.2f,%1.2f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\n",xx,yy,zz,dt[0],dt[1],dt[2],rdata[idx],rdata[idx+1]);
1895 return true;
1896 }
1897 else { // for subvolumes, not optimized yet
1898 size_t idx;
1899 float r, gg;
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 ++) {
1903 r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
1904 idx=image->get_complex_index(i,j,k,subx0,suby0,subz0,nx,ny,nz);
1905 gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
1906
1907 dt[0]+=gg*rdata[idx];
1908 dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
1909 dt[2]+=norm[idx/2];
1910 normsum+=gg;
1911 }
1912 }
1913 }
1914
1915 if (normsum==0) return false;
1916 return true;
1917 }
1918}
1919
1920
1922{
1923
1924 bool sqrtnorm=params.set_default("sqrtnorm",false);
1925 normalize_threed(sqrtnorm,true); // true is the wiener filter
1926
1927 if (doift) {
1928 image->do_ift_inplace();
1929 image->depad();
1930 image->process_inplace("xform.phaseorigin.tocenter");
1931 }
1932
1933 image->update();
1934
1935 if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) {
1936 if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
1937 tmp_data->write_image((const char *)params["savenorm"]);
1938 }
1939
1940 delete tmp_data;
1941 tmp_data=0;
1942 EMData *ret=image;
1943 image=0;
1944
1945 return ret;
1946}
1947
1948/*
1949void BaldwinWoolfordReconstructor::setup()
1950{
1951 //This is a bit of a hack - but for now it suffices
1952 params.set_default("mode","nearest_neighbor");
1953 WienerFourierReconstructor::setup();
1954 // Set up the Baldwin Kernel
1955 int P = (int)((1.0+0.25)*max_input_dim+1);
1956 float r = (float)(max_input_dim+1)/(float)P;
1957 dfreq = 0.2f;
1958 if (W != 0) delete [] W;
1959 int maskwidth = params.set_default("maskwidth",2);
1960 W = Util::getBaldwinGridWeights(maskwidth, (float)P, r,dfreq,0.5f,0.2f);
1961}
1962
1963EMData* BaldwinWoolfordReconstructor::finish(bool doift)
1964{
1965 tmp_data->write_image("density.mrc");
1966 image->process_inplace("xform.fourierorigin.tocorner");
1967 image->do_ift_inplace();
1968 image->depad();
1969 image->process_inplace("xform.phaseorigin.tocenter");
1970
1971 if ( (bool) params.set_default("postmultiply", false) )
1972 {
1973 cout << "POST MULTIPLYING" << endl;
1974 // now undo the Fourier convolution with real space division
1975 float* d = image->get_data();
1976 float N = (float) image->get_xsize()/2.0f;
1977 N *= N;
1978 size_t rnx = image->get_xsize();
1979 size_t rny = image->get_ysize();
1980 size_t rnxy = rnx*rny;
1981 int cx = image->get_xsize()/2;
1982 int cy = image->get_ysize()/2;
1983 int cz = image->get_zsize()/2;
1984 size_t idx;
1985 for (int k = 0; k < image->get_zsize(); ++k ){
1986 for (int j = 0; j < image->get_ysize(); ++j ) {
1987 for (int i =0; i < image->get_xsize(); ++i ) {
1988 float xd = (float)(i-cx); xd *= xd;
1989 float yd = (float)(j-cy); yd *= yd;
1990 float zd = (float)(k-cz); zd *= zd;
1991 float weight = exp((xd+yd+zd)/N);
1992 idx = k*rnxy + j*rnx + i;
1993 d[idx] *= weight;
1994 }
1995 }
1996 }
1997 }
1998 image->update();
1999 return image;
2000}
2001
2002#include <iomanip>
2003
2004// int BaldwinWoolfordReconstructor::insert_slice_weights(const Transform& t3d)
2005// {
2006// bool fftodd = image->is_fftodd();
2007// int rnx = nx-2*!fftodd;
2008//
2009// float y_scale = 1.0, x_scale = 1.0;
2010//
2011// if ( ny != rnx )
2012// {
2013// if ( rnx > ny ) y_scale = (float) rnx / (float) ny;
2014// else x_scale = (float) ny / (float) rnx;
2015// }
2016//
2017// int tnx = tmp_data->get_xsize();
2018// int tny = tmp_data->get_ysize();
2019// int tnz = tmp_data->get_zsize();
2020//
2021// vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
2022// for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
2023// Transform n3d = t3d*(*it);
2024//
2025// for (int y = 0; y < tny; y++) {
2026// for (int x = 0; x < tnx; x++) {
2027//
2028// float rx = (float) x;
2029// float ry = (float) y;
2030//
2031// if ( ny != rnx )
2032// {
2033// if ( rnx > ny ) ry *= y_scale;
2034// else rx *= x_scale;
2035// }
2036// // float xx = rx * n3d[0][0] + (ry - tny/2) * n3d[1][0];
2037// // float yy = rx * n3d[0][1] + (ry - tny/2) * n3d[1][1];
2038// // float zz = rx * n3d[0][2] + (ry - tny/2) * n3d[1][2];
2039//
2040// Vec3f coord(rx,(ry - tny/2),0);
2041// coord = coord*n3d; // transpose multiplication
2042// float xx = coord[0];
2043// float yy = coord[1];
2044// float zz = coord[2];
2045//
2046// if (xx < 0 ){
2047// xx = -xx;
2048// yy = -yy;
2049// zz = -zz;
2050// }
2051//
2052// yy += tny/2;
2053// zz += tnz/2;
2054// insert_density_at(xx,yy,zz);
2055// }
2056// }
2057// }
2058//
2059// return 0;
2060// }
2061
2062void BaldwinWoolfordReconstructor::insert_density_at(const float& x, const float& y, const float& z)
2063{
2064 int xl = Util::fast_floor(x);
2065 int yl = Util::fast_floor(y);
2066 int zl = Util::fast_floor(z);
2067
2068 // w is the windowing width
2069 int w = params.set_default("maskwidth",2);
2070 float wsquared = (float) w*w;
2071 float dw = 1.0f/w;
2072 dw *= dw;
2073
2074 // w minus one - this control the number of
2075 // pixels/voxels to the left of the main pixel
2076 // that will have density
2077 int wmox = w-1;
2078 int wmoy = w-1;
2079 int wmoz = w-1;
2080
2081 // If any coordinate is incedental with a vertex, then
2082 // make sure there is symmetry in density accruing.
2083 // i.e. the window width must be equal in both directions
2084 if ( ((float) xl) == x ) wmox = w;
2085 if ( ((float) yl) == y ) wmoy = w;
2086 if ( ((float) zl) == z ) wmoz = w;
2087
2088 float* d = tmp_data->get_data();
2089 int tnx = tmp_data->get_xsize();
2090 int tny = tmp_data->get_ysize();
2091 int tnz = tmp_data->get_zsize();
2092 size_t tnxy = tnx*tny;
2093
2094 int mode = params.set_default("mode","nearest_neighbor");
2095
2096 for(int k = zl-wmoz; k <= zl+w; ++k ) {
2097 for(int j = yl-wmoy; j <= yl+w; ++j) {
2098 for( int i = xl-wmox; i <= xl+w; ++i) {
2099 float fac = 1.0;
2100 int ic = i, jc = j, kc = k;
2101
2102 // Fourier space is periodic, which is enforced
2103 // by the next 6 if statements. These if statements
2104 // assume that the Fourier DC components is at
2105 // (0,ny/2,nz/2).
2106 if ( i <= 0 ) {
2107
2108 if ( x != 0 && i == 0 ) fac = 1.0;
2109 else if ( x == 0 && i < 0) continue;
2110// if (i < 0 ) ic = -i;
2111 if (i < 0 ) {
2112 continue;
2113 ic = -i;
2114 jc = tny-jc;
2115 kc = tnz-kc;
2116 }
2117 }
2118 if ( ic >= tnx ) ic = 2*tnx-ic-1;
2119 if ( jc < 0 ) jc = tny+jc;
2120 if ( jc >= tny ) jc = jc-tny;
2121 if ( kc < 0 ) kc = tnz+kc;
2122 if ( kc >= tnz ) kc = kc-tnz;
2123// if ( ic >= tnx ) continue;
2124// if ( jc < 0 ) continue;
2125// if ( jc >= tny ) continue;
2126// if ( kc < 0 ) continue;
2127// if ( kc >= tnz ) continue;
2128 // This shouldn't happen
2129 // Debug remove later
2130 if ( ic < 0 ) { cout << "wo 1" << endl; }
2131 if ( ic >= tnx ){ cout << "wo 2" << endl; }
2132 if ( jc < 0 ) { cout << "wo 3" << endl; }
2133 if ( jc >= tny ) { cout << "wo 4" << endl; }
2134 if ( kc < 0 ) { cout << "wo 5" << endl; }
2135 if ( kc >= tnz ) { cout << "wo 6" << endl; }
2136
2137
2138 float zd = (z-(float)k);
2139 float yd = (y-(float)j);
2140 float xd = (x-(float)i);
2141 zd *= zd; yd *= yd; xd *= xd;
2142 float distsquared = xd+yd+zd;
2143 // We enforce a spherical kernel
2144 if ( mode == 1 && distsquared > wsquared ) continue;
2145
2146// float f = fac*exp(-dw*(distsquared));
2147 float f = fac*exp(-2.467f*(distsquared));
2148 // Debug - this error should never occur.
2149 if ( (kc*tnxy+jc*tnx+ic) >= tnxy*tnz ) throw OutofRangeException(0,tnxy*tnz,kc*tnxy+jc*tnx+ic, "in density insertion" );
2150 d[kc*tnxy+jc*tnx+ic] += f;
2151 }
2152 }
2153 }
2154}
2155
2156int BaldwinWoolfordReconstructor::insert_slice(const EMData* const input_slice, const Transform & t, const float weight)
2157{
2158 Transform * rotation;
2159 if ( input_slice->has_attr("xform.projection") ) {
2160 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
2161 } else {
2162 rotation = new Transform(t); // assignment operator
2163 }
2164 Transform tmp(*rotation);
2165 tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly
2166
2167 Vec2f trans = tmp.get_trans_2d();
2168 float scale = tmp.get_scale();
2169 bool mirror = tmp.get_mirror();
2170 EMData* slice = 0;
2171 if (trans[0] != 0 || trans[1] != 0 || scale != 1.0 ) {
2172 slice = input_slice->process("xform",Dict("transform",&tmp));
2173 } else if ( mirror == true ) {
2174 slice = input_slice->process("xform.flip",Dict("axis","x"));
2175 }
2176 if ( slice == 0 ) {
2177 slice = input_slice->process("xform.phaseorigin.tocorner");
2178 } else {
2179 slice->process_inplace("xform.phaseorigin.tocorner");
2180 }
2181
2182 slice->do_fft_inplace();
2183 slice->process_inplace("xform.fourierorigin.tocenter");
2184 float *dat = slice->get_data();
2185 float dt[2];
2186
2187 bool fftodd = image->is_fftodd();
2188 int rnx = nx-2*!fftodd;
2189
2190 float y_scale = 1.0, x_scale = 1.0;
2191
2192 if ( ny != rnx )
2193 {
2194 if ( rnx > ny ) y_scale = (float) rnx / (float) ny;
2195 else x_scale = (float) ny / (float) rnx;
2196 }
2197
2198 int tnx = tmp_data->get_xsize();
2199 int tny = tmp_data->get_ysize();
2200 int tnz = tmp_data->get_zsize();
2201
2202 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
2203// float weight = params.set_default("weight",1.0f);
2204
2205 rotation->set_scale(1.0); rotation->set_mirror(false); rotation->set_trans(0,0,0);
2206 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
2207 Transform t3d = (*rotation)*(*it);
2208
2209 for (int y = 0; y < tny; y++) {
2210 for (int x = 0; x < tnx; x++) {
2211 float rx = (float) x;
2212 float ry = (float) y;
2213
2214 if ( ny != rnx )
2215 {
2216 if ( rnx > ny ) ry *= y_scale;
2217 else rx *= x_scale;
2218 }
2219
2220// float xx = rx * n3d[0][0] + (ry - tny/2) * n3d[1][0];
2221// float yy = rx * n3d[0][1] + (ry - tny/2) * n3d[1][1];
2222// float zz = rx * n3d[0][2] + (ry - tny/2) * n3d[1][2];
2223
2224 Vec3f coord(rx,(ry - tny/2),0);
2225 coord = coord*t3d; // transpose multiplication
2226 float xx = coord[0];
2227 float yy = coord[1];
2228 float zz = coord[2];
2229
2230
2231 float cc = 1;
2232 if (xx < 0 ){
2233 xx = -xx;
2234 yy = -yy;
2235 zz = -zz;
2236 cc = -1;
2237 }
2238
2239 yy += tny/2;
2240 zz += tnz/2;
2241
2242 int idx = x * 2 + y * (slice->get_xsize());
2243 dt[0] = dat[idx];
2244 dt[1] = cc * dat[idx+1];
2245
2246 insert_pixel(xx,yy,zz,dt);
2247 }
2248 }
2249 }
2250
2251 if(rotation) {delete rotation; rotation=0;}
2252 delete slice;
2253
2254 return 0;
2255}
2256
2257void BaldwinWoolfordReconstructor::insert_pixel(const float& x, const float& y, const float& z, const float dt[2])
2258{
2259 int xl = Util::fast_floor(x);
2260 int yl = Util::fast_floor(y);
2261 int zl = Util::fast_floor(z);
2262
2263 // w is the windowing width
2264 int w = params.set_default("maskwidth",2);
2265 float wsquared = (float) w*w;
2266 float dw = 1.0f/w;
2267 dw *= dw;
2268
2269 int wmox = w-1;
2270 int wmoy = w-1;
2271 int wmoz = w-1;
2272
2273 // If any coordinate is incedental with a vertex, then
2274 // make sure there is symmetry in density accruing.
2275 // i.e. the window width must be equal in both directions
2276 if ( ((float) xl) == x ) wmox = w;
2277 if ( ((float) yl) == y ) wmoy = w;
2278 if ( ((float) zl) == z ) wmoz = w;
2279
2280 float* we = tmp_data->get_data();
2281 int tnx = tmp_data->get_xsize();
2282 int tny = tmp_data->get_ysize();
2283 int tnz = tmp_data->get_zsize();
2284 int tnxy = tnx*tny;
2285
2286 int rnx = 2*tnx;
2287 int rnxy = 2*tnxy;
2288
2289 int mode = params.set_default("mode","nearest_neighbor");
2290
2291 float* d = image->get_data();
2292 for(int k = zl-wmoz; k <= zl+w; ++k ) {
2293 for(int j = yl-wmoy; j <= yl+w; ++j) {
2294 for( int i = xl-wmox; i <= xl+w; ++i) {
2295 float fac = 1.0;
2296 int ic = i, jc = j, kc = k;
2297
2298 // Fourier space is periodic, which is enforced
2299 // by the next 6 if statements. These if statements
2300 // assume that the Fourier DC component is at
2301 // (0,ny/2,nz/2).
2302 float negfac=1.0;
2303 if ( i <= 0 ) {
2304 if ( x != 0 && i == 0 ) fac = 1.0;
2305 else if ( x == 0 && i < 0) continue;
2306 if (i < 0 ) {
2307 continue;
2308 ic = -i;
2309 jc = tny-jc;
2310 kc = tnz-kc;
2311 negfac=-1.0;
2312 }
2313 }
2314 if ( ic >= tnx ) ic = 2*tnx-ic-1;
2315 if ( jc < 0 ) jc = tny+jc;
2316 if ( jc >= tny ) jc = jc-tny;
2317 if ( kc < 0 ) kc = tnz+kc;
2318 if ( kc >= tnz ) kc = kc-tnz;
2319// if ( ic >= tnx ) continue;
2320// if ( jc < 0 ) continue;
2321// if ( jc >= tny ) continue;
2322// if ( kc < 0 ) continue;
2323// if ( kc >= tnz ) continue;
2324
2325 float zd = (z-(float)k);
2326 float yd = (y-(float)j);
2327 float xd = (x-(float)i);
2328 zd *= zd; yd *= yd; xd *= xd;
2329 float distsquared = xd+yd+zd;
2330// float f = fac*exp(-dw*(distsquared));
2331 float f = fac*exp(-2.467f*(distsquared));
2332 float weight = f/we[kc*tnxy+jc*tnx+ic];
2333 // debug - this error should never occur
2334 if ( (kc*rnxy+jc*rnx+2*ic+1) >= rnxy*tnz ) throw OutofRangeException(0,rnxy*tnz,kc*rnxy+jc*rnx+2*ic+1, "in pixel insertion" );
2335 size_t k = kc*rnxy+jc*rnx+2*ic;
2336
2337 float factor, dist,residual;
2338 int sizeW,sizeWmid,idx;
2339 switch (mode) {
2340 case 0:
2341 d[k] += weight*f*dt[0];
2342 d[k+1] += negfac*weight*f*dt[1];
2343 cout << "hello" << endl;
2344 break;
2345
2346 case 1:
2347 // We enforce a spherical kernel
2348 if ( distsquared > wsquared ) continue;
2349
2350 sizeW = (int)(1+2*w/dfreq);
2351 sizeWmid = sizeW/2;
2352
2353 dist = sqrtf(distsquared);
2354 idx = (int)(sizeWmid + dist/dfreq);
2355 if (idx >= sizeW) throw InvalidValueException(idx, "idx was greater than or equal to sizeW");
2356 residual = dist/dfreq - (int)(dist/dfreq);
2357 if ( fabs(residual) > 1) throw InvalidValueException(residual, "Residual was too big");
2358
2359 factor = (W[idx]*(1.0f-residual)+W[idx+1]*residual)*weight;
2360
2361 d[k] += dt[0]*factor;
2362 d[k+1] += dt[1]*factor;
2363 break;
2364
2365 default:
2366 throw InvalidValueException(mode, "The mode was unsupported in BaldwinWoolfordReconstructor::insert_pixel");
2367 break;
2368 }
2369 }
2370 }
2371 }
2372}
2373
2374// void BaldwinWoolfordReconstructor::insert_pixel(const float& x, const float& y, const float& z, const float dt[2])
2375// {
2376// int xl = Util::fast_floor(x);
2377// int yl = Util::fast_floor(y);
2378// int zl = Util::fast_floor(z);
2379//
2380// // w is the windowing width
2381// int w = params.set_default("maskwidth",2);
2382// float dw = 1.0/w;
2383// dw *= dw;
2384// // dw = 2;
2385// // cout << w << endl;
2386// // int w = 3;
2387// // w minus one - this control the number of
2388// // pixels/voxels to the left of the main pixel
2389// // that will have density
2390// int wmox = w-1;
2391// int wmoy = w-1;
2392// int wmoz = w-1;
2393//
2394// // If any coordinate is incedental with a vertex, then
2395// // make sure there is symmetry in density accruing.
2396// // i.e. the window width must be equal in both directions
2397// if ( ((float) xl) == x ) wmox = w;
2398// if ( ((float) yl) == y ) wmoy = w;
2399// if ( ((float) zl) == z ) wmoz = w;
2400//
2401// float* d = tmp_data->get_data();
2402// int tnx = tmp_data->get_xsize();
2403// int tny = tmp_data->get_ysize();
2404// int tnz = tmp_data->get_zsize();
2405// int tnxy = tnx*tny;
2406//
2407// float weight = 1.0;
2408// //
2409// for(int k = zl-wmoz; k <= zl+w; ++k ) {
2410// for(int j = yl-wmoy; j <= yl+w; ++j) {
2411// for( int i = xl-wmox; i <= xl+w; ++i) {
2412// float fac = 1.0;
2413// int ic = i, jc = j, kc = k;
2414//
2415// // Fourier space is periodic, which is enforced
2416// // by the next 6 if statements. These if statements
2417// // assume that the Fourier DC components is at
2418// // (0,ny/2,nz/2).
2419// if ( i <= 0 ) {
2420// if ( x != 0 && i == 0 ) fac = 1.0;
2421// else if ( x == 0 && i < 0) continue;
2422// // if (i < 0 ) ic = -i;
2423// if (i < 0 ) {
2424// ic = -i;
2425// jc = tny-jc;
2426// kc = tnz-kc;
2427// }
2428// }
2429// if ( ic >= tnx ) ic = 2*tnx-ic-1;
2430// if ( jc < 0 ) jc = tny+jc;
2431// if ( jc >= tny ) jc = jc-tny;
2432// if ( kc < 0 ) kc = tnz+kc;
2433// if ( kc >= tnz ) kc = kc-tnz;
2434// // This shouldn't happen
2435// // Debug remove later
2436// if ( ic < 0 ) { cout << "wo 1" << endl; }
2437// if ( ic >= tnx ){ cout << "wo 2" << endl; }
2438// if ( jc < 0 ) { cout << "wo 3" << endl; }
2439// if ( jc >= tny ) { cout << "wo 4" << endl; }
2440// if ( kc < 0 ) { cout << "wo 5" << endl; }
2441// if ( kc >= tnz ) { cout << "wo 6" << endl; }
2442//
2443//
2444// float zd = (z-(float)k);
2445// float yd = (y-(float)j);
2446// float xd = (x-(float)i);
2447// zd *= zd; yd *= yd; xd *= xd;
2448// // Debug - this error should never occur.
2449// if ( (kc*tnxy+jc*tnx+ic) >= tnxy*tnz ) throw OutofRangeException(0,tnxy*tnz,kc*tnxy+jc*tnx+ic, "in weight determination insertion" );
2450// // float f = fac*exp(-dw*(xd+yd+zd)*0.5);
2451// float f = exp(-2.467*(xd+yd+zd));
2452// weight += f*(d[kc*tnxy+jc*tnx+ic] - f);
2453// }
2454// }
2455// }
2456// weight = 1.0/weight;
2457// int rnx = 2*tnx;
2458// int rnxy = 2*tnxy;
2459// d = image->get_data();
2460// for(int k = zl-wmoz; k <= zl+w; ++k ) {
2461// for(int j = yl-wmoy; j <= yl+w; ++j) {
2462// for( int i = xl-wmox; i <= xl+w; ++i) {
2463// float fac = 1.0;
2464// int ic = i, jc = j, kc = k;
2465//
2466// // Fourier space is periodic, which is enforced
2467// // by the next 6 if statements. These if statements
2468// // assume that the Fourier DC components is at
2469// // (0,ny/2,nz/2).
2470// float negfac=1.0;
2471// if ( i <= 0 ) {
2472// if ( x != 0 && i == 0 ) fac = 1.0;
2473// else if ( x == 0 && i < 0) continue;
2474// if (i < 0 ) {
2475// continue;
2476// ic = -i;
2477// jc = tny-jc;
2478// kc = tnz-kc;
2479// negfac=-1.0;
2480// }
2481// }
2482// if ( ic >= tnx ) ic = 2*tnx-ic-1;
2483// if ( jc < 0 ) jc = tny+jc;
2484// if ( jc >= tny ) jc = jc-tny;
2485// if ( kc < 0 ) kc = tnz+kc;
2486// if ( kc >= tnz ) kc = kc-tnz;
2487// // This shouldn't happen
2488// // Debug remove later
2489//
2490//
2491// float zd = (z-(float)k);
2492// float yd = (y-(float)j);
2493// float xd = (x-(float)i);
2494// zd *= zd; yd *= yd; xd *= xd;
2495// // float f = fac*exp(-dw*(xd+yd+zd));
2496// float f = exp(-4.934*(xd+yd+zd));
2497// // Debug - this error should never occur.
2498// if ( (kc*rnxy+jc*rnx+2*ic+1) >= rnxy*tnz ) throw OutofRangeException(0,rnxy*tnz,kc*rnxy+jc*rnx+2*ic+1, "in pixel insertion" );
2499//
2500// d[kc*rnxy+jc*rnx+2*ic] += weight*f*dt[0];
2501// d[kc*rnxy+jc*rnx+2*ic+1] += negfac*weight*f*dt[1];
2502// }
2503// }
2504// }
2505// }
2506*/
2507
2509{
2510 vector<int> size=params["size"];
2511 nx = size[0];
2512 ny = size[1];
2513 nz = size[2];
2514 image = new EMData(nx,ny,nz);
2515 image->to_zero();
2516
2517 if (!slices.empty()) {
2518 for (std::vector<EMData *>::iterator it = slices.begin() ; it != slices.end(); ++it) delete *it;
2519 slices.clear();
2520 }
2521 xforms.clear();
2522}
2523
2525{
2526
2527// EMData* return_slice = slice->process("normalize.edgemean");
2528// return_slice->process_inplace("filter.linearfourier");
2529
2530 EMData* return_slice;
2531
2532// return_slice = slice->process("filter.linearfourier");
2533 return_slice = slice->copy();
2534
2535// Transform tmp(t);
2536// tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly
2537// Vec2f trans = tmp.get_trans_2d();
2538// float scale = tmp.get_scale();
2539// bool mirror = tmp.get_mirror();
2540// if (trans[0] != 0 || trans[1] != 0 || scale != 1.0 ) {
2541// return_slice->transform(tmp);
2542// }
2543// if ( mirror == true ) {
2544// return_slice->process_inplace("xform.flip",Dict("axis","x"));
2545// }
2546
2547 return return_slice;
2548}
2549
2550int RealMedianReconstructor::insert_slice(const EMData* const input, const Transform &t, const float weight)
2551{
2552 if (!input) {
2553 LOGERR("try to insert NULL slice");
2554 return 1;
2555 }
2556
2557 if (input->get_xsize() != input->get_ysize() || input->get_xsize() != nx) {
2558 LOGERR("tried to insert image that was not correction dimensions");
2559 return 1;
2560 }
2561
2562 slices.push_back(input->copy());
2563 xforms.push_back(t);
2564
2565 return 0;
2566}
2567
2568// This is just a dummy function. It does not return meaningful numbers
2569int RealMedianReconstructor::determine_slice_agreement(EMData* input_slice, const Transform & arg, const float weight,bool sub)
2570{
2571 // Are these exceptions really necessary? (d.woolford)
2572 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
2573
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);
2577
2578 return 0;
2579
2580}
2581
2583{
2584
2585 int mode = params.set_default("mode",0); // confusing, mode means mode of operation
2586
2587 // We do the actual reconstruction here from all of the slices
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);
2595
2596 // single interpolated value from 2-D
2597 float val=(*it)->sget_value_at_interp(imvec[0]+nx/2,imvec[1]+ny/2);
2598 if (val!=0) vals.push_back(val);
2599
2600 // uses all 4 of the nearest values from the 2-D image
2601// int xx=floor(imvec[0]+nx/2);
2602// int yy=floor(imvec[1]+ny/2);
2603// float val=(*it)->sget_value_at(xx,yy);
2604// if (val!=0) vals.push_back(val);
2605// val=(*it)->sget_value_at(xx+1,yy);
2606// if (val!=0) vals.push_back(val);
2607// val=(*it)->sget_value_at(xx+1,yy+1);
2608// if (val!=0) vals.push_back(val);
2609// val=(*it)->sget_value_at(xx,yy+1);
2610// if (val!=0) vals.push_back(val);
2611 }
2612 if (vals.size()==0) continue;
2613
2614// // debugging, save all values for specific voxels to look at statistics
2615// if ((x==515&&y==461&&z==281) ||
2616// (x==516&&y==462&&z==299) ||
2617// (x==516&&y==461&&z==262) ||
2618// (x==212&&y==409&&z==304)) {
2619// char fsp[80];
2620// sprintf(fsp,"stat_%d_%d_%d.txt",x,y,z);
2621// FILE *out=fopen(fsp,"w");
2622// for (std::vector<float>::iterator v = vals.begin(); v!=vals.end(); ++v) fprintf(out,"%f\n",*v);
2623// fclose(out);
2624// }
2625
2626 switch(mode) {
2627 // median, the namesake of the reconstructor
2628 case 0:
2629 std::sort(vals.begin(),vals.end());
2630 //printf("%d %d %d %d\n",x,y,z,vals.size());
2631 image->set_value_at(x,y,z,vals[vals.size()/2]); // for even sizes, not quite right, and should include possibility of local average
2632 break;
2633
2634 // This is something like a mode. It progressively subdivides the numeric axis in halves
2635 // until reaching a small number of values which are then averaged
2636 case 1:
2637 {
2638 // We only use this method if we have at least 6 possible values for the voxel, otherwise we just average
2639 if (vals.size()<6) {
2640 float val=0;
2641 for (int i=0; i<vals.size(); i++) val+=vals[i];
2642 image->set_value_at(x,y,z,val/vals.size());
2643 break;
2644 }
2645 // Sort the list to simplify the remaining logic
2646 std::sort(vals.begin(),vals.end());
2647 int v0=0,v1=vals.size(),vs=v1/2; // bracket values, v0 first, v1 last+1, vs is first value in second zone
2648
2649 // we continue until we have a zone containing 4 or fewer values
2650 //int count=0;
2651 while (v1-v0>8) {
2652 float cen=(vals[v0]+vals[v1-1])/2.0f; // value halfway between the extremes
2653 for (vs=v0; vs<v1; vs++) { if (vals[vs]>=cen) break; } // find the first value in the upper zone
2654 // shift to the zone containing more values
2655 if (v1-vs>vs-v0) v0=vs;
2656 else v1=vs;
2657 //printf("%d %d %d %d\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",v0,vs,v1,vals.size(),cen,vals[0],vals[v0],vals[vs-1],vals[vs],vals[v1],vals[vals.size()-1]);
2658 //count++;
2659 //if (count>20) exit(1);
2660 }
2661 float val=0;
2662 for (int i=v0; i<v1; i++) val+=vals[i];
2663 image->set_value_at(x,y,z,val/(v1-v0));
2664 break;
2665 }
2666 // median with a bias towards zero (positive and negative), using 1/3 or 2/3 position in the sorted array
2667 case 2:
2668 {
2669 if (vals.size()<6) {
2670 float val=0;
2671 for (int i=0; i<vals.size(); i++) val+=vals[i];
2672 image->set_value_at(x,y,z,val/vals.size());
2673 break;
2674 }
2675 std::sort(vals.begin(),vals.end());
2676 //printf("%d %d %d %d\n",x,y,z,vals.size());
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;
2680 image->set_value_at(x,y,z,v); // for even sizes, not quite right, and should include possibility of local average
2681 break;
2682 }
2683 // bias the median towards the half of the points spanning the narrower range
2684 case 3:
2685 {
2686 if (vals.size()<6) {
2687 float val=0;
2688 for (int i=0; i<vals.size(); i++) val+=vals[i];
2689 image->set_value_at(x,y,z,val/vals.size());
2690 break;
2691 }
2692 std::sort(vals.begin(),vals.end());
2693 //printf("%d %d %d %d\n",x,y,z,vals.size());
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;
2698 image->set_value_at(x,y,z,v); // for even sizes, not quite right, and should include possibility of local average
2699 break;
2700 }
2701 case 4:
2702 {
2703 // too few values, return median
2704 if (vals.size()<10) {
2705 std::sort(vals.begin(),vals.end());
2706 image->set_value_at(x,y,z,vals[vals.size()/2]); // for even sizes, not quite right, and should include possibility of local average
2707 break;
2708 }
2709
2710 //compute stats
2711 int n=vals.size();
2712 float mean=0,sigma=0;
2713 for (int i=0; i<n; i++) { mean+=vals[i]; sigma+=vals[i]*vals[i]; }
2714 mean/=n;
2715 sigma=sqrt(sigma/n-mean);
2716 // low & high refer to angle not value
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;
2719
2720 // If both extreme tilts roughly agree, then we assume that's the right range to use and we average values near that
2721 // by computing the mean of all values within 1 sigma of all values
2722 if (fabs(low-high)<sigma) {
2723 mean=(low+high)*4.0f;
2724 int c=8;
2725 for (int i=4; i<n-4; i++) {
2726 if (fabs(vals[i]-mean/c)<sigma) { mean+=vals[i]; c++; }
2727 }
2728 image->set_value_at(x,y,z,mean/c);
2729 }
2730 // This means the extremes disagree and we have to pick which value is best, we go with min(fabs)
2731 else {
2732 if (fabs(low)<fabs(high)) mean=low*4.0;
2733 else mean=high*4.0;
2734 int c=4;
2735 for (int i=4; i<n-4; i++) {
2736 if (fabs(vals[i]-mean/c)<sigma) { mean+=vals[i]; c++; }
2737 }
2738 image->set_value_at(x,y,z,mean/c);
2739 }
2740 break;
2741 }
2742 }
2743 }
2744 }
2745 }
2746
2747
2748// transform->set_scale(1.0);
2749// transform->set_mirror(false);
2750// transform->set_trans(0,0,0);
2751// transform->invert();
2752
2753// tmp->transform(t.inverse()); // This was incorrect. inverse() was missing
2754// image->add(*tmp);
2755
2756// if(transform) {delete transform; transform=0;}
2757
2758 // apply symmetry if requested
2759 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params["sym"]);
2760 vector<Transform> syms = sym->get_syms();
2761
2762 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
2763
2764// it->printme();
2765 Transform t=*it;
2766 EMData *tmpcopy = image->process("xform",Dict("transform",(EMObject)&t));
2767 image->add(*tmpcopy);
2768 delete tmpcopy;
2769 }
2770
2771 image->mult(1.0f/(float)sym->get_nsym());
2772 delete sym;
2773// if (image->get_xsize()==image->get_ysize() && image->get_ysize()==image->get_zsize()) {
2774// image->process_inplace("mask.sharp",Dict("outer_radius",image->get_xsize()/2-1));
2775// }
2776// else printf("No masking %d %d %d\n",image->get_xsize(),image->get_ysize(),image->get_zsize());
2777
2778 // erase out copy of the projection data
2779 if (!slices.empty()) {
2780 for (std::vector<EMData *>::iterator it = slices.begin() ; it != slices.end(); ++it) delete *it;
2781 slices.clear();
2782 }
2783 xforms.clear();
2784
2785 EMData *ret = image;
2786 if (!doift) {
2787 ret=image->do_fft();
2788 ret->process_inplace("xform.phaseorigin.tocorner");
2789 delete image;
2790 }
2791 image = 0 ;
2792 return ret;
2793}
2794
2795
2797{
2798 image = new EMData();
2799 vector<int> size=params["size"];
2800 nx = size[0];
2801 ny = size[1];
2802 nz = size[2];
2803 image->set_size(nx, ny, nz);
2804}
2805
2807{
2808
2809// EMData* return_slice = slice->process("normalize.edgemean");
2810// return_slice->process_inplace("filter.linearfourier");
2811
2812 EMData* return_slice;
2813
2814 return_slice = slice->process("filter.linearfourier");
2815// return_slice = slice->copy();
2816
2817// Transform tmp(t);
2818// tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly
2819// Vec2f trans = tmp.get_trans_2d();
2820// float scale = tmp.get_scale();
2821// bool mirror = tmp.get_mirror();
2822// if (trans[0] != 0 || trans[1] != 0 || scale != 1.0 ) {
2823// return_slice->transform(tmp);
2824// }
2825// if ( mirror == true ) {
2826// return_slice->process_inplace("xform.flip",Dict("axis","x"));
2827// }
2828
2829 return return_slice;
2830}
2831
2832int BackProjectionReconstructor::insert_slice(const EMData* const input, const Transform &t, const float)
2833{
2834 if (!input) {
2835 LOGERR("try to insert NULL slice");
2836 return 1;
2837 }
2838
2839 if (input->get_xsize() != input->get_ysize() || input->get_xsize() != nx) {
2840 LOGERR("tried to insert image that was not correction dimensions");
2841 return 1;
2842 }
2843
2844// Transform * transform;
2845// if ( input->has_attr("xform.projection") ) {
2846// transform = (Transform*) (input->get_attr("xform.projection")); // assignment operator
2847// } else {
2848// transform = new Transform(t); // assignment operator
2849// }
2850 EMData* slice = preprocess_slice(input, t);
2851
2852 // Clearly weight isn't a useful concept in back-projection without compensating with an exact-filter
2853// float weight = params["weight"];
2854// slice->mult(weight);
2855
2856 EMData *tmp = new EMData();
2857 tmp->set_size(nx, ny, nz);
2858
2859 float *slice_data = slice->get_data();
2860 float *tmp_data = tmp->get_data();
2861
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);
2866 }
2867 tmp->update();
2868
2869// transform->set_scale(1.0);
2870// transform->set_mirror(false);
2871// transform->set_trans(0,0,0);
2872// transform->invert();
2873
2874 tmp->transform(t.inverse()); // This was incorrect. inverse() was missing
2875 image->add(*tmp);
2876
2877// if(transform) {delete transform; transform=0;}
2878 delete tmp;
2879 delete slice;
2880
2881 return 0;
2882}
2883
2884int BackProjectionReconstructor::determine_slice_agreement(EMData* input_slice, const Transform & arg, const float weight,bool sub)
2885{
2886 // Are these exceptions really necessary? (d.woolford)
2887 if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
2888
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);
2892
2893 return 0;
2894
2895}
2896
2898{
2899
2900 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params["sym"]);
2901 vector<Transform> syms = sym->get_syms();
2902
2903 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
2904
2905// it->printme();
2906 Transform t=*it;
2907 EMData *tmpcopy = image->process("xform",Dict("transform",(EMObject)&t));
2908 image->add(*tmpcopy);
2909 delete tmpcopy;
2910 }
2911
2912 image->mult(1.0f/(float)sym->get_nsym());
2913 delete sym;
2914 if (image->get_xsize()==image->get_ysize() && image->get_ysize()==image->get_zsize()) {
2915 image->process_inplace("mask.sharp",Dict("outer_radius",image->get_xsize()/2-1));
2916 }
2917 else printf("No masking %d %d %d\n",image->get_xsize(),image->get_ysize(),image->get_zsize());
2918
2919 EMData *ret = image;
2920 image = 0 ;
2921 return ret;
2922}
2923
2924EMData* EMAN::padfft_slice( const EMData* const slice, const Transform& t, int npad )
2925{
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; // If 2, it means it is a Fourier file.
2932
2933 if( ndim==2 && (nx-extension)!=ny ) {
2934 LOGERR("Input image must be square!");
2935 throw ImageDimensionException("Input image must be square!");
2936 }
2937
2938 EMData* padfftslice = NULL;
2939 if( padffted == 0) {
2940 // process 2D slice or 1D line -- subtract the average outside of the circle, zero-pad, fft extend, and fft
2941 EMData* temp = slice->average_circ_sub();
2942
2943 padfftslice = temp->norm_pad( false, npad );
2944 checked_delete( temp );
2945
2946 padfftslice->do_fft_inplace();
2947 } else {
2948 padfftslice = new EMData(*slice);
2949 }
2950
2951 // shift the projection
2952 Vec2f trans = t.get_trans_2d();
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));
2956
2957 int remove = slice->get_attr_default("remove", 0);
2958 padfftslice->set_attr( "remove", remove );
2959 padfftslice->center_origin_fft();
2960 return padfftslice;
2961}
2962
2963
2964//####################################################################################
2965
2966
2968
2969float max2d( int kc, const vector<float>& pow_a )
2970{
2971 float max = 0.0;
2972 for( int i=-kc; i <= kc; ++i ) {
2973 for( int j=-kc; j <= kc; ++j ) {
2974 if( i==0 && j==0 ) continue;
2975 {
2976 int c = 2*kc+1 - std::abs(i) - std::abs(j);
2977 max = max + pow_a[c];
2978 }
2979 }
2980 }
2981 return max;
2982}
2983
2984float max3d( int kc, const vector<float>& pow_a )
2985{
2986 float max = 0.0;
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;
2991 // if( i!=0 )
2992 {
2993 int c = 3*kc+1 - std::abs(i) - std::abs(j) - std::abs(k);
2994 max = max + pow_a[c];
2995 // max = max + c * c;
2996 // max = max + c;
2997 }
2998 }
2999 }
3000 }
3001 return max;
3002}
3003
3004
3005#define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
3006
3007void circumfnn( EMData* win , int npad)
3008{
3009 float *tw = win->get_data();
3010 // correct for the fall-off of NN interpolation using sinc functions
3011 // mask and subtract circumference average
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);
3017
3018 int IP = ix/2+1;
3019 int JP = iy/2+1;
3020 int KP = iz/2+1;
3021
3022 // sinc functions tabulated for fall-off
3023 float* sincx = new float[IP+1];
3024 float* sincy = new float[JP+1];
3025 float* sincz = new float[KP+1];
3026
3027 sincx[0] = 1.0f;
3028 sincy[0] = 1.0f;
3029 sincz[0] = 1.0f;
3030
3031 float cor;
3032 if( npad == 1 ) cor = 1.0;
3033 else cor = 4.0;
3034
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);
3046 }
3047 }
3048
3049 delete[] sincx;
3050 delete[] sincy;
3051 delete[] sincz;
3052
3053 float TNR = 0.0f;
3054 size_t m = 0;
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) {
3060 TNR += tw(i,j,k);
3061 ++m;
3062 }
3063 }
3064 }
3065 }
3066
3067 TNR /=float(m);
3068
3069
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;
3076 }
3077 }
3078 }
3079
3080}
3081
3082
3083void circumftrl( EMData* win , int npad)
3084{
3085 float *tw = win->get_data();
3086 // correct for the fall-off of tri-linear interpolation using sinc^2 functions
3087 // mask and subtract circumference average
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);
3093
3094 int IP = ix/2+1;
3095 int JP = iy/2+1;
3096 int KP = iz/2+1;
3097
3098 // sinc functions tabulated for fall-off
3099 float* sincx = new float[IP+1];
3100 float* sincy = new float[JP+1];
3101 float* sincz = new float[KP+1];
3102
3103 sincx[0] = 1.0f;
3104 sincy[0] = 1.0f;
3105 sincz[0] = 1.0f;
3106
3107 float cor;
3108 if( npad == 1 ) cor = 1.0;
3109 else cor = 4.0;
3110
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);
3122 }
3123 }
3124
3125 delete[] sincx;
3126 delete[] sincy;
3127 delete[] sincz;
3128
3129 float TNR = 0.0f;
3130 size_t m = 0;
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) {
3136 TNR += tw(i,j,k);
3137 ++m;
3138 }
3139 }
3140 }
3141 }
3142
3143 TNR /=float(m);
3144
3145
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;
3152 }
3153 }
3154 }
3155
3156}
3157
3158
3159
3160//####################################################################################
3161
3162//** nn4 reconstructor
3163
3165{
3166 m_volume = NULL;
3167 m_wptr = NULL;
3168}
3169
3170nn4Reconstructor::nn4Reconstructor( const string& symmetry, int size, int npad )
3171{
3172 setup( symmetry, size, npad );
3174 //print_params();
3175}
3176
3178{
3179 //if( m_delete_volume ) checked_delete(m_volume);
3180
3181 //if( m_delete_weight ) checked_delete( m_wptr );
3182
3183 //checked_delete( m_result );
3184}
3185
3186
3188{
3189
3190 int size = params["size"];
3191 int npad = params["npad"];
3192
3193 string symmetry;
3194 if( params.has_key("symmetry") ) symmetry = params["symmetry"].to_str();
3195 else symmetry = "c1";
3196
3197 if( params.has_key("ndim") ) m_ndim = params["ndim"];
3198 else m_ndim = 3;
3199
3200 if( params.has_key( "snr" ) ) m_osnr = 1.0f/float( params["snr"] );
3201 else m_osnr = 0.0;
3202
3203 setup( symmetry, size, npad );
3204}
3205
3206void nn4Reconstructor::setup( const string& symmetry, int size, int npad )
3207{
3209 m_wghta = 0.2f;
3210
3211 m_symmetry = symmetry;
3212 m_npad = npad;
3214
3215 m_vnx = size;
3216 m_vny = size;
3217 m_vnz = (m_ndim==3) ? size : 1;
3218
3219 m_vnxp = size*npad;
3220 m_vnyp = size*npad;
3221 m_vnzp = (m_ndim==3) ? size*npad : 1;
3222
3223 m_vnxc = m_vnxp/2;
3224 m_vnyc = m_vnyp/2;
3225 m_vnzc = (m_ndim==3) ? m_vnzp/2 : 1;
3226
3229
3230}
3231
3232
3234 int offset = 2 - m_vnxp%2;
3235
3236 m_volume = params["fftvol"];
3237
3238 if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
3239 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
3240 m_volume->to_zero();
3241 }
3242 // ----------------------------------------------------------------
3243 // Added by Zhengfan Yang on 03/15/07
3244 // Original author: please check whether my revision is correct and
3245 // other Reconstructor need similiar revision.
3246 if ( m_vnxp % 2 == 0 ) m_volume->set_fftodd(0);
3247 else m_volume->set_fftodd(1);
3248 // ----------------------------------------------------------------
3249
3250 m_volume->set_nxc(m_vnxp/2);
3251 m_volume->set_complex(true);
3252 m_volume->set_ri(true);
3253 m_volume->set_fftpad(true);
3254 m_volume->set_attr("npad", m_npad);
3255 m_volume->set_array_offsets(0,1,1);
3256}
3257
3259
3260 m_wptr = params["weight"];
3261
3262 if( m_wptr->get_xsize() != m_vnxc+1 &&
3263 m_wptr->get_ysize() != m_vnyp &&
3264 m_wptr->get_zsize() != m_vnzp ) {
3265 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
3266 m_wptr->to_zero();
3267 }
3268 m_wptr->set_array_offsets(0,1,1);
3269}
3270
3271void printImage( const EMData* line )
3272{
3273 Assert( line->get_zsize()==1 );
3274
3275
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) );
3280 printf( "\n" );
3281 }
3282}
3283
3284
3285
3286int nn4Reconstructor::insert_slice(const EMData* const slice, const Transform& t, const float weight) {
3287 // sanity checks
3288 if (!slice) {
3289 LOGERR("try to insert NULL slice");
3290 return 1;
3291 }
3292
3293 int padffted= slice->get_attr_default( "padffted", 0 );
3294 if( m_ndim==3 ) {
3295 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx) ) {
3296 // FIXME: Why doesn't this throw an exception?
3297 LOGERR("Tried to insert a slice that is the wrong size.");
3298 return 1;
3299 }
3300 } else {
3301 Assert( m_ndim==2 );
3302 if( slice->get_ysize() !=1 ) {
3303 LOGERR( "for 2D reconstruction, a line is excepted" );
3304 return 1;
3305 }
3306 }
3307 if( weight > 0.0f ) {
3308
3309 EMData* padfft = padfft_slice( slice, t, m_npad );
3310
3311 if( m_ndim==3 ) {
3312 insert_padfft_slice( padfft, t, weight );
3313 } else {
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 );
3321 if( xnew < 0.0 ) {
3322 xnew *= -1;
3323 ynew *= -1;
3324 btqi *= -1;
3325 }
3326
3327 int ixn = int(xnew+0.5+m_vnxp) - m_vnxp;
3328 int iyn = int(ynew+0.5+m_vnyp) - m_vnyp;
3329
3330 if(iyn < 0 ) iyn += m_vnyp;
3331
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;
3335 }
3336
3337 }
3338 checked_delete( padfft );
3339 }
3340 return 0;
3341}
3342
3343int nn4Reconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, float weight )
3344{
3345 Assert( padfft != NULL );
3346
3347 vector<Transform> tsym = t.get_sym_proj(m_symmetry);
3348 for (unsigned int isym=0; isym < tsym.size(); isym++) m_volume->nn( m_wptr, padfft, tsym[isym], weight);
3349
3350 return 0;
3351}
3352
3353
3355
3356 if( m_ndim == 3 ) {
3357 m_volume->symplane0(m_wptr);
3358 } else {
3359 for( int i=1; i <= m_vnyp; ++i ) {
3360
3361 if( (*m_wptr)(0, i, 1)==0.0 ) {
3362 int j = m_vnyp + 1 - i;
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);
3366 }
3367 }
3368 }
3369
3370
3371 int box = 7;
3372 int kc = (box-1)/2;
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);
3375 pow_a.back()=0.0;
3376
3377 float alpha = 0.0;
3378 if( m_ndim==3) {
3379 int vol = box*box*box;
3380 float max = max3d( kc, pow_a );
3381 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
3382 } else {
3383 int ara = box*box;
3384 float max = max2d( kc, pow_a );
3385 alpha = ( 1.0f - 1.0f/(float)ara ) / max;
3386 }
3387 int ix,iy,iz;
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) {//(*v) should be treated as complex!!
3392 float tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+m_osnr);
3393 if( m_weighting == ESTIMATE && false) { // HERE
3394 int cx = ix;
3395 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
3396 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
3397 float sum = 0.0;
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;
3404
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;
3409 if( nbrcx < 0 ) {
3410 nbrcx = -nbrcx;
3411 nbrcy = -nbrcy;
3412 nbrcz = -nbrcz;
3413 }
3414 int nbrix = nbrcx;
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];
3420 }
3421 }
3422 }
3423 }
3424 float wght = 1.0f / ( 1.0f - alpha * sum );
3425 tmp = tmp * wght;
3426 }
3427//cout<<" mvol "<<ix<<" "<<iy<<" "<<iz<<" "<<(*m_volume)(2*ix,iy,iz)<<" "<<(*m_volume)(2*ix+1,iy,iz)<<" "<<tmp<<" "<<m_osnr<<endl;
3428 (*m_volume)(2*ix,iy,iz) *= tmp;
3429 (*m_volume)(2*ix+1,iy,iz) *= tmp;
3430 }
3431 }
3432 }
3433 }
3434
3435 //if(m_ndim==2) printImage( m_volume );
3436
3437 // back fft
3438 m_volume->do_ift_inplace();
3439
3440 // EMData* win = m_volume->window_center(m_vnx);
3441 int npad = m_volume->get_attr("npad");
3442 m_volume->depad();
3443 circumfnn( m_volume, npad );
3444 m_volume->set_array_offsets( 0, 0, 0 );
3445
3446 return 0;
3447}
3448#undef tw
3449
3450
3452{
3453 m_volume = NULL;
3454 m_wptr = NULL;
3455}
3456
3457nn4_rectReconstructor::nn4_rectReconstructor( const string& symmetry, int size, int npad )
3458{
3459 m_volume = NULL;
3460 m_wptr = NULL;
3461 setup( symmetry, size, npad );
3463 print_params();
3464}
3465
3467{
3468 //if( m_delete_volume ) checked_delete(m_volume);
3469
3470 //if( m_delete_weight ) checked_delete( m_wptr );
3471
3472 //checked_delete( m_result );
3473}
3474
3475
3477{
3478 m_sizeofprojection = params["sizeprojection"];
3479 int npad = params["npad"];
3480 m_count=0;
3481
3482 string symmetry;
3483 if( params.has_key("symmetry") ) symmetry = params["symmetry"].to_str();
3484 else symmetry = "c1";
3485
3486 if( params.has_key("ndim") ) m_ndim = params["ndim"];
3487 else m_ndim = 3;
3488
3489 if( params.has_key( "snr" ) ) m_osnr = 1.0f/float( params["snr"] );
3490 else m_osnr = 0.0;
3491
3492 setup( symmetry, m_sizeofprojection, npad );
3493}
3494
3495void nn4_rectReconstructor::setup( const string& symmetry, int sizeprojection, int npad )
3496{
3498 m_wghta = 0.2f;
3499 m_symmetry = symmetry;
3500 m_npad = npad;
3502
3503 if( params.has_key("sizex") ) m_vnx = params["sizex"];
3504 else if(params.has_key("xratio"))
3505 {
3506 float temp=params["xratio"];
3507 m_vnx=int(float(sizeprojection)*temp);
3508 }
3509 else m_vnx=sizeprojection;
3510
3511 if( params.has_key("sizey") ) m_vny = params["sizey"];
3512 else if (params.has_key("yratio"))
3513 {
3514 float temp=params["yratio"];
3515 m_vny=int(float(sizeprojection)*temp);
3516 }
3517 else m_vny=sizeprojection;
3518
3519 if( params.has_key("sizez") )
3520 m_vnz = params["sizez"];
3521 else
3522 if (params.has_key("zratio"))
3523 {
3524 float temp=params["zratio"];
3525 m_vnz=int(float(sizeprojection)*temp);
3526 }
3527 else
3528 m_vnz = (m_ndim==3) ? sizeprojection : 1;
3529
3530 m_xratio=float(m_vnx)/float(sizeprojection);
3531 m_yratio=float(m_vny)/float(sizeprojection);
3532 m_zratio=float(m_vnz)/float(sizeprojection);
3533
3534 m_vnxp = m_vnx*npad;
3535 m_vnyp = m_vny*npad;
3536 m_vnzp = (m_ndim==3) ? m_vnz*npad : 1;
3537
3538 m_vnxc = m_vnxp/2;
3539 m_vnyc = m_vnyp/2;
3540 m_vnzc = (m_ndim==3) ? m_vnzp/2 : 1;
3541
3544}
3545
3546
3548 int offset = 2 - m_vnxp%2;
3549
3550 m_volume = params["fftvol"];
3551
3552 if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
3553 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
3554 m_volume->to_zero();
3555 }
3556 // ----------------------------------------------------------------
3557 // Added by Zhengfan Yang on 03/15/07
3558 // Original author: please check whether my revision is correct and
3559 // other Reconstructor need similiar revision.
3560 if ( m_vnxp % 2 == 0 ) m_volume->set_fftodd(0);
3561 else m_volume->set_fftodd(1);
3562 // ----------------------------------------------------------------
3563
3564 m_volume->set_nxc(m_vnxp/2);
3565 m_volume->set_complex(true);
3566 m_volume->set_ri(true);
3567 m_volume->set_fftpad(true);
3568 m_volume->set_attr("npad", m_npad);
3569 m_volume->set_array_offsets(0,1,1);
3570}
3571
3573
3574 m_wptr = params["weight"];
3575
3576 if( m_wptr->get_xsize() != m_vnxc+1 &&
3577 m_wptr->get_ysize() != m_vnyp &&
3578 m_wptr->get_zsize() != m_vnzp ) {
3579 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
3580 m_wptr->to_zero();
3581 }
3582
3583 m_wptr->set_array_offsets(0,1,1);
3584}
3585
3586int nn4_rectReconstructor::insert_slice(const EMData* const slice, const Transform& t, const float weight) {
3587 // sanity checks
3588
3589
3590 if (!slice) {
3591 LOGERR("try to insert NULL slice");
3592 return 1;
3593 }
3594
3595 int padffted= slice->get_attr_default( "padffted", 0 );
3596 if( m_ndim==3 ) {
3597 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_sizeofprojection) ) {
3598 // FIXME: Why doesn't this throw an exception?
3599 LOGERR("Tried to insert a slice that is the wrong size.");
3600 return 1;
3601 }
3602 }
3603 if (m_ndim==2) {
3604 if( slice->get_ysize() !=1 ) {
3605 LOGERR( "for 2D reconstruction, a line is excepted" );
3606 return 1;
3607 }
3608 }
3609 if( weight > 0.0f ) {
3610
3611 EMData* padfft = padfft_slice( slice, t, m_npad );
3612
3613 //Assert( mult > 0 );
3614
3615 if( m_ndim==3 ) {
3616 insert_padfft_slice( padfft, t, weight );
3617 } else {
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;
3622 int loop_range;
3623 float temp1,temp2;
3624 temp1=m_xratio*cos(alpha)*float(m_sizeofprojection*m_npad)/2;
3625 temp2=m_yratio*sin(alpha)*float(m_sizeofprojection*m_npad)/2;
3626 ellipse_length=sqrt(temp1*temp1+temp2*temp2);
3627 ellipse_length_int=int(ellipse_length);
3628 ellipse_step=0.5f*(m_sizeofprojection*m_npad)/float(ellipse_length_int);
3629 loop_range=ellipse_length_int;
3630 cos_alpha=temp1/ellipse_length;
3631 sin_alpha=temp2/ellipse_length;
3632 if(m_count%100==0) {
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;
3637 std::cout<<"x and y ratio=="<<m_xratio<<" "<<m_yratio<<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;
3644 }
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 );
3655 if( xnew < 0.0 ) {
3656 xnew *= -1;
3657 ynew *= -1;
3658 btqi *= -1;
3659 }
3660
3661 int ixn = int(xnew+0.5+m_vnxp) - m_vnxp;
3662 int iyn = int(ynew+0.5+m_vnyp) - m_vnyp;
3663
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;
3670 }
3671
3672
3673 }
3674 checked_delete( padfft );
3675 return 0;
3676 } else return 0;
3677}
3678
3679
3680
3681
3683{
3684 Assert( padded != NULL );
3685
3686 vector<Transform> tsym = t.get_sym_proj(m_symmetry);
3687 for (unsigned int isym=0; isym < tsym.size(); isym++)
3688 m_volume->insert_rect_slice(m_wptr, padded, tsym[isym], m_sizeofprojection, m_xratio, m_yratio, m_zratio, m_npad, weight);
3689
3690 return 0;
3691
3692}
3693
3694#define tw(i,j,k) tw[ i-1 + (j-1+(k-1)*iy)*ix ]
3695void circumfnn_rect( EMData* win , int npad)
3696{
3697 float *tw = win->get_data();
3698 // correct for the fall-off
3699 // mask and subtract circumfnnerence average
3700 int ix = win->get_xsize();
3701 int iy = win->get_ysize();
3702 int iz = win->get_zsize();
3703
3704 int IP = ix/2+1;
3705 int JP = iy/2+1;
3706 int KP = iz/2+1;
3707
3708 // sinc functions tabulated for fall-off
3709 float* sincx = new float[IP+1];
3710 float* sincy = new float[JP+1];
3711 float* sincz = new float[KP+1];
3712
3713 sincx[0] = 1.0f;
3714 sincy[0] = 1.0f;
3715 sincz[0] = 1.0f;
3716
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);
3728 }
3729 }
3730
3731 delete[] sincx;
3732 delete[] sincy;
3733 delete[] sincz;
3734
3735
3736
3737 float dxx = 1.0f/float(0.25*ix*ix);
3738 float dyy = 1.0f/float(0.25*iy*iy);
3739
3740
3741
3742 float LR2=(float(ix)/2-1)*(float(ix)/2-1)*dxx;
3743
3744 float TNR = 0.0f;
3745 size_t m = 0;
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) {
3751 TNR += tw(i,j,k);
3752 ++m;
3753 }
3754 }
3755 }
3756 }
3757
3758 TNR /=float(m);
3759
3760
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;
3767 }
3768 }
3769 }
3770
3771}
3772#undef tw
3773
3775{
3776
3777 if( m_ndim==3 ) {
3778 m_volume->symplane0_rect(m_wptr);
3779 } else {
3780 for( int i=1; i <= m_vnyp; ++i ) {
3781
3782 if( (*m_wptr)(0, i, 1)==0.0 ) {
3783 int j = m_vnyp + 1 - i;
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);
3787 }
3788 }
3789 }
3790
3791
3792 int box = 7;
3793 int kc = (box-1)/2;
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);
3796 pow_a.back()=0.0;
3797
3798 float alpha = 0.0;
3799 if( m_ndim==3) {
3800 int vol = box*box*box;
3801 float max = max3d( kc, pow_a );
3802 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
3803 } else {
3804 int ara = box*box;
3805 float max = max2d( kc, pow_a );
3806 alpha = ( 1.0f - 1.0f/(float)ara ) / max;
3807 }
3808
3809 int ix,iy,iz;
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) {//(*v) should be treated as complex!!
3814 float tmp;
3815 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+m_osnr);
3816
3817 if( m_weighting == ESTIMATE ) {
3818 int cx = ix;
3819 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
3820 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
3821 float sum = 0.0;
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;
3828
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;
3833 if( nbrcx < 0 ) {
3834 nbrcx = -nbrcx;
3835 nbrcy = -nbrcy;
3836 nbrcz = -nbrcz;
3837 }
3838 int nbrix = nbrcx;
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];
3844 }
3845 }
3846 }
3847 }
3848 float wght = 1.0f / ( 1.0f - alpha * sum );
3849 tmp = tmp * wght;
3850 }
3851 (*m_volume)(2*ix,iy,iz) *= tmp;
3852 (*m_volume)(2*ix+1,iy,iz) *= tmp;
3853 }
3854 }
3855 }
3856 }
3857
3858 //if(m_ndim==2) printImage( m_volume );
3859
3860 // back fft
3861 m_volume->do_ift_inplace();
3862
3863
3864 int npad = m_volume->get_attr("npad");
3865 m_volume->depad();
3866 circumfnn_rect( m_volume, npad );
3867 m_volume->set_array_offsets( 0, 0, 0 );
3868
3869 return 0;
3870}
3871
3872
3873// Added By Zhengfan Yang on 03/16/07
3874// Beginning of the addition
3875// --------------------------------------------------------------------------------
3876
3878{
3879 m_volume = NULL;
3880 m_wptr = NULL;
3881 m_wptr2 = NULL;
3882}
3883
3884nnSSNR_Reconstructor::nnSSNR_Reconstructor( const string& symmetry, int size, int npad)
3885{
3886 m_volume = NULL;
3887 m_wptr = NULL;
3888 m_wptr2 = NULL;
3889
3890 setup( symmetry, size, npad );
3891}
3892
3894{
3895 //if( m_delete_volume ) checked_delete(m_volume);
3896
3897 //if( m_delete_weight ) checked_delete( m_wptr );
3898
3899 //if( m_delete_weight2 ) checked_delete( m_wptr2 );
3900
3901 //checked_delete( m_result );
3902}
3903
3905{
3906 int size = params["size"];
3907 int npad = params["npad"];
3908
3909 string symmetry;
3910 if( params.has_key("symmetry") ) symmetry = params["symmetry"].to_str();
3911 else symmetry = "c1";
3912
3913 setup( symmetry, size, npad );
3914}
3915
3916void nnSSNR_Reconstructor::setup( const string& symmetry, int size, int npad )
3917{
3918
3920 m_wghta = 0.2f;
3921 m_wghtb = 0.004f;
3922
3923 m_symmetry = symmetry;
3924 m_npad = npad;
3926
3927 m_vnx = size;
3928 m_vny = size;
3929 m_vnz = size;
3930
3931 m_vnxp = size*npad;
3932 m_vnyp = size*npad;
3933 m_vnzp = size*npad;
3934
3935 m_vnxc = m_vnxp/2;
3936 m_vnyc = m_vnyp/2;
3937 m_vnzc = m_vnzp/2;
3938
3942
3943}
3944
3946
3947 m_volume = params["fftvol"];
3948 m_volume->set_size(m_vnxp+ 2 - m_vnxp%2,m_vnyp,m_vnzp);
3949 m_volume->to_zero();
3950 if ( m_vnxp % 2 == 0 ) m_volume->set_fftodd(0);
3951 else m_volume->set_fftodd(1);
3952
3953 m_volume->set_nxc(m_vnxc);
3954 m_volume->set_complex(true);
3955 m_volume->set_ri(true);
3956 m_volume->set_fftpad(true);
3957 m_volume->set_attr("npad", m_npad);
3958 m_volume->set_array_offsets(0,1,1);
3959}
3960
3962
3963 m_wptr = params["weight"];
3964 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
3965 m_wptr->to_zero();
3966 m_wptr->set_array_offsets(0,1,1);
3967}
3968
3970
3971 m_wptr2 = params["weight2"];
3972 m_wptr2->set_size(m_vnxc+1,m_vnyp,m_vnzp);
3973 m_wptr2->to_zero();
3974 m_wptr2->set_array_offsets(0,1,1);
3975}
3976
3977
3978int nnSSNR_Reconstructor::insert_slice(const EMData* const slice, const Transform& t, const float weight) {
3979 // sanity checks
3980 if (!slice) {
3981 LOGERR("try to insert NULL slice");
3982 return 1;
3983 }
3984 if( weight > 0.0f ) {
3985 int padffted=slice->get_attr_default( "padffted", 0 );
3986
3987 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx) ) {
3988 // FIXME: Why doesn't this throw an exception?
3989 LOGERR("Tried to insert a slice that has wrong size.");
3990 return 1;
3991 }
3992
3993 EMData* padfft = padfft_slice( slice, t, m_npad );
3994
3995 insert_padfft_slice( padfft, t, weight );
3996
3997 checked_delete( padfft );
3998 return 0;
3999 } else return 0;
4000}
4001
4002int nnSSNR_Reconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, float weight )
4003{
4004 Assert( padfft != NULL );
4005 // insert slice for all symmetry related positions
4006 vector<Transform> tsym = t.get_sym_proj(m_symmetry);
4007 for (int isym=0; isym < m_nsym; isym++) m_volume->nn_SSNR( m_wptr, m_wptr2, padfft, tsym[isym], weight);
4008
4009 return 0;
4010}
4011
4012
4014{
4015/*
4016 I changed the code on 05/15/2008 so it only returns variance.
4017 Lines commented out are marked by //#
4018 The version prior to the currect changes is r1.190
4019 PAP
4020*/
4021 int kz, ky;
4022 //#int iix, iiy, iiz;
4023 int box = 7;
4024 int kc = (box-1)/2;
4025 float alpha = 0.0;
4026 float argx, argy, argz;
4027 vector< float > pow_a( 3*kc+1, 1.0 );
4028 float w = params["w"];
4029 EMData* SSNR = params["SSNR"];
4030 //#EMData* vol_ssnr = new EMData();
4031 //#vol_ssnr->set_size(m_vnxp, m_vnyp, m_vnzp);
4032 //#vol_ssnr->to_zero();
4033 //# new line follow
4034 EMData* vol_ssnr = params["vol_ssnr"];
4035 vol_ssnr->set_size(m_vnxp+ 2 - m_vnxp%2, m_vnyp ,m_vnzp);
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);
4043 //#
4044
4045 float dx2 = 1.0f/float(m_vnxc)/float(m_vnxc);
4046 float dy2 = 1.0f/float(m_vnyc)/float(m_vnyc);
4047 float dz2 = 1.0f/Util::get_max(float(m_vnzc),1.0f)/Util::get_max(float(m_vnzc),1.0f);
4048 int inc = Util::round(float(Util::get_max(m_vnxc,m_vnyc,m_vnzc))/w);
4049
4050 SSNR->set_size(inc+1,4,1);
4051
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];
4056 float wght = 1.0f;
4057 for (int i = 0; i <= inc; i++) {
4058 nom[i] = 0.0f;
4059 denom[i] = 0.0f;
4060 nn[i] = 0;
4061 ka[i] = 0;
4062 }
4063
4064 m_volume->symplane1(m_wptr, m_wptr2);
4065
4066 if ( m_weighting == ESTIMATE ) {
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);
4069 pow_a[3*kc] = 0.0;
4070 float max = max3d( kc, pow_a );
4071 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
4072 }
4073
4074 for (int iz = 1; iz <= m_vnzp; iz++) {
4075 if ( iz-1 > m_vnzc ) kz = iz-1-m_vnzp; else kz = iz-1;
4076 argz = float(kz*kz)*dz2;
4077 for (int iy = 1; iy <= m_vnyp; iy++) {
4078 if ( iy-1 > m_vnyc ) ky = iy-1-m_vnyp; else ky = iy-1;
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);
4083 int r = Util::round(float(inc)*argx);
4084 if ( r >= 0 && Kn > 4.5f ) {
4085 if ( m_weighting == ESTIMATE ) {
4086 int cx = ix;
4087 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
4088 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
4089
4090 float sum = 0.0;
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;
4100 if( nbrcx < 0 ) {
4101 nbrcx = -nbrcx;
4102 nbrcy = -nbrcy;
4103 nbrcz = -nbrcz;
4104 }
4105 int nbrix = nbrcx;
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];
4111 }
4112 }
4113 }
4114 }
4115 wght = 1.0f / ( 1.0f - alpha * sum );
4116 } // end of ( m_weighting == ESTIMATE )
4117 float nominator = std::norm(m_volume->cmplx(ix,iy,iz)/Kn);
4118 float denominator = ((*m_wptr2)(ix,iy,iz)-nominator)/(Kn-1.0f);
4119 // Skip Friedel related values
4120 if( (ix>0 || (kz>=0 && (ky>=0 || kz!=0)))) {
4121 if ( r <= inc ) {
4122 nom[r] += nominator*wght;
4123 denom[r] += denominator/Kn*wght;
4124 nn[r] += 2;
4125 ka[r] += int(Kn);
4126 }
4127/*
4128 //#float tmp = Util::get_max(nominator/denominator/Kn-1.0f,0.0f);
4129 // Create SSNR as a 3D array (-n/2:n/2+n%2-1)
4130 iix = m_vnxc + ix; iiy = m_vnyc + ky; iiz = m_vnzc + kz;
4131 if( iix >= 0 && iix < m_vnxp && iiy >= 0 && iiy < m_vnyp && iiz >= 0 && iiz < m_vnzp )
4132 (*vol_ssnr)(iix, iiy, iiz) = tmp;
4133 // Friedel part
4134 iix = m_vnxc - ix; iiy = m_vnyc - ky; iiz = m_vnzc - kz;
4135 if( iix >= 0 && iix < m_vnxp && iiy >= 0 && iiy < m_vnyp && iiz >= 0 && iiz < m_vnzp )
4136 (*vol_ssnr)(iix, iiy, iiz) = tmp;
4137*/
4138
4139 }
4140 (*vol_ssnr)(2*ix, iy-1, iz-1) = nominator*wght;//Kn;//denominator*wght;//
4141 //(*vol_ssnr)(2*ix, iy-1, iz-1) = real(m_volume->cmplx(ix,iy,iz))*wght/Kn;
4142 //(*vol_ssnr)(2*ix+1, iy-1, iz-1) = imag(m_volume->cmplx(ix,iy,iz))*wght/Kn;
4143 } // end of Kn>4.5
4144 }
4145 }
4146 }
4147
4148 for (int i = 0; i <= inc; i++) {
4149 (*SSNR)(i,0,0) = nom[i];
4150 (*SSNR)(i,1,0) = denom[i]; // variance
4151 (*SSNR)(i,2,0) = static_cast<float>(nn[i]);
4152 (*SSNR)(i,3,0) = static_cast<float>(ka[i]);
4153 }
4154 vol_ssnr->update();
4155
4156 delete[] nom;
4157 delete[] denom;
4158 delete[] nn;
4159 delete[] ka;
4160
4161 return 0;
4162}
4163
4164// -----------------------------------------------------------------------------------
4165// End of this addition
4166
4167//####################################################################################
4168//** nn4 ctf reconstructor
4169
4171{
4172 m_volume = NULL;
4173 m_wptr = NULL;
4174}
4175
4176nn4_ctfReconstructor::nn4_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign )
4177{
4178 setup( symmetry, size, npad, snr, sign );
4179}
4180
4182{
4183 //if( m_delete_volume ) checked_delete(m_volume);
4184
4185 //if( m_delete_weight ) checked_delete( m_wptr );
4186
4187 //checked_delete( m_result );
4188}
4189
4191{
4192 if( ! params.has_key("size") ) throw std::logic_error("Error: image size is not given");
4193
4194 int size = params["size"];
4195 int npad = params.has_key("npad") ? int(params["npad"]) : 2;
4196 // int sign = params.has_key("sign") ? int(params["sign"]) : 1;
4197 int sign = 1;
4198 string symmetry = params.has_key("symmetry")? params["symmetry"].to_str() : "c1";
4199
4200 float snr = params["snr"];
4201
4202 m_varsnr = params.has_key("varsnr") ? int(params["varsnr"]) : 0;
4203 setup( symmetry, size, npad, snr, sign );
4204
4205}
4206
4207void nn4_ctfReconstructor::setup( const string& symmetry, int size, int npad, float snr, int sign )
4208{
4210 if( params.has_key("weighting") ) {
4211 if( int( params["weighting"])==0 ) m_weighting = NONE;
4212 }
4213
4214
4215
4216 m_wghta = 0.2f;
4217 m_wghtb = 0.004f;
4218
4219 m_symmetry = symmetry;
4220 m_npad = npad;
4221 m_sign = sign;
4223
4224 m_snr = snr;
4225
4226 m_vnx = size;
4227 m_vny = size;
4228 m_vnz = size;
4229
4230 m_vnxp = size*npad;
4231 m_vnyp = size*npad;
4232 m_vnzp = size*npad;
4233
4234 m_vnxc = m_vnxp/2;
4235 m_vnyc = m_vnyp/2;
4236 m_vnzc = m_vnzp/2;
4237
4240}
4241
4243 int offset = 2 - m_vnxp%2;
4244
4245 m_volume = params["fftvol"];
4246
4247 if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
4248 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
4249 m_volume->to_zero();
4250 }
4251
4252 if ( m_vnxp % 2 == 0 ) m_volume->set_fftodd(0);
4253 else m_volume->set_fftodd(1);
4254 m_volume->set_nxc(m_vnxp/2);
4255 m_volume->set_complex(true);
4256 m_volume->set_ri(true);
4257 m_volume->set_fftpad(true);
4258 m_volume->set_attr("npad", m_npad);
4259 m_volume->set_array_offsets(0,1,1);
4260}
4261
4263{
4264 m_wptr = params["weight"];
4265
4266 if( m_wptr->get_xsize() != m_vnxc+1 && m_wptr->get_ysize() != m_vnyp && m_wptr->get_zsize() != m_vnzp ) {
4267 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
4268 m_wptr->to_zero();
4269 }
4270
4271 m_wptr->set_array_offsets(0,1,1);
4272
4273}
4274
4275int nn4_ctfReconstructor::insert_slice(const EMData* const slice, const Transform& t, const float weight)
4276{
4277 // sanity checks
4278 if (!slice) {
4279 LOGERR("try to insert NULL slice");
4280 return 1;
4281 }
4282 if( weight > 0.0f ) {
4283 int buffed = slice->get_attr_default( "buffed", 0 );
4284 if( buffed > 0 ) {
4285 insert_buffed_slice( slice, weight );
4286 return 0;
4287 }
4288
4289 int padffted= slice->get_attr_default("padffted", 0);
4290 if( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx) ) {
4291 // FIXME: Why doesn't this throw an exception?
4292 LOGERR("Tried to insert a slice that is the wrong size.");
4293 return 1;
4294 }
4295
4296 EMData* padfft = padfft_slice( slice, t, m_npad );
4297
4298 float tmp = padfft->get_attr_default("ctf_applied", 0);
4299 int ctf_applied = (int) tmp;
4300
4301 // Generate 2D CTF (EMData object)
4302 //ctf_store_real::init( padfft->get_ysize(), padfft->get_attr( "ctf" ) );
4303 //EMData* ctf2d = ctf_store_real::get_ctf_real(); //This is in 2D projection plane
4304 int winsize = padfft->get_ysize();
4305 Ctf* ctf = padfft->get_attr( "ctf" );
4306 Dict params = ctf->to_dict();
4307 EMData* ctf2d = Util::ctf_img_real(winsize , winsize, 1, params["defocus"], params["apix"],params["voltage"], params["cs"], \
4308 params["ampcont"], params["bfactor"], params["dfdiff"], params["dfang"], 1);
4309 params.clear();
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();
4313
4314 size_t size = (size_t)nx*ny*nz;
4315 if (!ctf_applied) {
4316 for (int i = 0; i < size; ++i) padfft->cmplx(i) *= ctf2d_ptr[i]; // Multiply padfft by CTF
4317 }
4318
4319 for (int i = 0; i < size; ++i) ctf2d_ptr[i] *= ctf2d_ptr[i]; // Square 2D CTF
4320
4321 insert_padfft_slice(padfft, ctf2d, t, weight);
4322
4323 checked_delete( ctf2d );
4324 checked_delete( padfft );
4325
4326 }
4327 return 0;
4328}
4329
4330int nn4_ctfReconstructor::insert_buffed_slice( const EMData* buffed, float weight )
4331{
4332 const float* bufdata = buffed->get_data();
4333 float* cdata = m_volume->get_data();
4334 float* wdata = m_wptr->get_data();
4335
4336 int npoint = buffed->get_xsize()/4;
4337 for( int i=0; i < npoint; ++i ) {
4338
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;
4344/*
4345 std::cout << "pos1, pos2, ctfv1, ctfv2, ctf2: ";
4346 std::cout << pos1 << " " << bufdata[5*i+1] << " " << bufdata[5*i+2] << " ";
4347 std::cout << pos2 << " " << bufdata[5*i+4] << std::endl;
4348 */
4349 }
4350 return 0;
4351}
4352
4353int nn4_ctfReconstructor::insert_padfft_slice( EMData* padfft, EMData* ctf2d2, const Transform& t, float weight)
4354{
4355 Assert( padfft != NULL );
4356
4357 vector<float> abc_list;
4358 int abc_list_len = 0;
4359 if (m_volume->has_attr("smear"))
4360 {
4361 abc_list = m_volume->get_attr("smear");
4362 abc_list_len = abc_list.size();
4363 }
4364
4365 vector<Transform> tsym = t.get_sym_proj(m_symmetry);
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);
4369 else
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]);
4372 }
4373 }
4374 return 0;
4375}
4376
4378{
4379 m_volume->set_array_offsets(0, 1, 1);
4380 m_wptr->set_array_offsets(0, 1, 1);
4381 m_volume->symplane0_ctf(m_wptr);
4382
4383 int box = 7;
4384 int vol = box*box*box;
4385 int kc = (box-1)/2;
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);
4388 pow_a[3*kc]=0.0;
4389
4390
4391 float max = max3d( kc, pow_a );
4392 float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
4393 float osnr = 1.0f/m_snr;
4394
4395 // normalize
4396 int ix,iy,iz;
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) {//(*v) should be treated as complex!!
4401 float tmp=0.0f;
4402 if( m_varsnr ) {
4403 int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
4404 int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
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);//*m_sign;
4407 } else {
4408 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+osnr);//*m_sign;
4409 }
4410
4411 if( m_weighting == ESTIMATE ) {
4412 int cx = ix;
4413 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
4414 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
4415 float sum = 0.0;
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;
4425 if( nbrcx < 0 ) {
4426 nbrcx = -nbrcx;
4427 nbrcy = -nbrcy;
4428 nbrcz = -nbrcz;
4429 }
4430
4431 int nbrix = nbrcx;
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];
4437 // if(ix%20==0 && iy%20==0 && iz%20==0)
4438 // std::cout << boost::format( "%4d %4d %4d %4d %10.3f" ) % nbrix % nbriy % nbriz % c % sum << std::endl;
4439 }
4440 }
4441 }
4442 }
4443 float wght = 1.0f / ( 1.0f - alpha * sum );
4444/*
4445 if(ix%10==0 && iy%10==0)
4446 {
4447 std::cout << boost::format( "%4d %4d %4d " ) % ix % iy %iz;
4448 std::cout << boost::format( "%10.3f %10.3f %10.3f " ) % tmp % wght % sum;
4449 std:: << boost::format( "%10.3f %10.3e " ) % pow_b[r] % alpha;
4450 std::cout << std::endl;
4451 }
4452 */
4453 tmp = tmp * wght;
4454 }
4455 (*m_volume)(2*ix,iy,iz) *= tmp;
4456 (*m_volume)(2*ix+1,iy,iz) *= tmp;
4457 }
4458 }
4459 }
4460 }
4461
4462 // back fft
4463 m_volume->do_ift_inplace();
4464 int npad = m_volume->get_attr("npad");
4465 m_volume->depad();
4466 circumfnn( m_volume, npad );
4467 m_volume->set_array_offsets( 0, 0, 0 );
4468
4469 return 0;
4470}
4471
4472
4473//####################################################################################
4474//** nn4 ctfw reconstructor
4475
4477{
4478 m_volume = NULL;
4479 m_wptr = NULL;
4480}
4481
4482nn4_ctfwReconstructor::nn4_ctfwReconstructor( const string& symmetry, int size, int npad, float snr, int sign, int do_ctf )
4483{
4484 setup( symmetry, size, npad, snr, sign, do_ctf );
4485}
4486
4488{
4490
4491 //if( m_delete_weight ) checked_delete( m_wptr );
4492
4493 //checked_delete( m_result );
4494}
4495
4497{
4498 if( ! params.has_key("size") ) throw std::logic_error("Error: image size is not given");
4499
4500 int size = params["size"];
4501 int npad = params.has_key("npad") ? int(params["npad"]) : 2;
4502 // int sign = params.has_key("sign") ? int(params["sign"]) : 1;
4503 int sign = 1;
4504 string symmetry = params.has_key("symmetry")? params["symmetry"].to_str() : "c1";
4505
4506 float snr = params["snr"];
4507 int do_ctf = params["do_ctf"];
4508
4509 setup( symmetry, size, npad, snr, sign, do_ctf );
4510
4511}
4512
4513void nn4_ctfwReconstructor::setup( const string& symmetry, int size, int npad, float snr, int sign, int do_ctf )
4514{
4516 if( params.has_key("weighting") ) {
4517 if( int( params["weighting"])==0 ) m_weighting = NONE;
4518 }
4519
4520 m_wghta = 0.2f;
4521 m_wghtb = 0.004f;
4522
4523 m_symmetry = symmetry;
4524 m_npad = npad;
4525 m_sign = sign;
4527 m_do_ctf = do_ctf;
4528
4529 m_snr = snr;
4530
4531 m_vnx = size;
4532 m_vny = size;
4533 m_vnz = size;
4534
4535 //m_vnxp = size*npad;
4536 //m_vnyp = size*npad;
4537 //m_vnzp = size*npad;
4538 m_vnxp = size;
4539 m_vnyp = size;
4540 m_vnzp = size;
4541
4542 m_vnxc = m_vnxp/2;
4543 m_vnyc = m_vnyp/2;
4544 m_vnzc = m_vnzp/2;
4545
4548 m_refvol = params["refvol"];
4549
4550}
4551
4553 int offset = 2 - m_vnxp%2;
4554
4555 m_volume = params["fftvol"];
4556
4557 if( m_volume->get_xsize() != m_vnxp+offset && m_volume->get_ysize() != m_vnyp && m_volume->get_zsize() != m_vnzp ) {
4558 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
4559 m_volume->to_zero();
4560 }
4561
4562 if ( m_vnxp % 2 == 0 ) m_volume->set_fftodd(0);
4563 else m_volume->set_fftodd(1);
4564 m_volume->set_nxc(m_vnxp/2);
4565 m_volume->set_complex(true);
4566 m_volume->set_ri(true);
4567 m_volume->set_fftpad(true);
4568 m_volume->set_attr("npad", m_npad);
4569 m_volume->set_array_offsets(0,1,1);
4570}
4571
4573{
4574 m_wptr = params["weight"];
4575
4576 if( m_wptr->get_xsize() != m_vnxc+1 && m_wptr->get_ysize() != m_vnyp && m_wptr->get_zsize() != m_vnzp ) {
4577 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
4578 m_wptr->to_zero();
4579 }
4580
4581 m_wptr->set_array_offsets(0,1,1);
4582
4583}
4584
4585int nn4_ctfwReconstructor::insert_slice(const EMData* const slice, const Transform& t, const float weight)
4586{
4587 // sanity checks
4588 if (!slice) {
4589 LOGERR("try to insert NULL slice");
4590 return 1;
4591 }
4592 if(weight >0.0f) {
4593 /*
4594 int buffed = slice->get_attr_default( "buffed", 0 );
4595 if( buffed > 0 ) {
4596 insert_buffed_slice( slice, weight );
4597 return 0;
4598 }
4599 */
4600
4601 EMData* padfft = padfft_slice( slice, t, m_npad );
4602
4603 EMData* ctf2d = NULL;
4604 if( m_do_ctf == 1 ) {
4605 float tmp = padfft->get_attr_default("ctf_applied", 0);
4606 int ctf_applied = (int) tmp;
4607
4608 // Generate 2D CTF (EMData object)
4609 //ctf_store_real::init( padfft->get_ysize(), padfft->get_attr( "ctf" ) );
4610 int winsize = padfft->get_ysize();
4611 Ctf* ctf = padfft->get_attr( "ctf" );
4612 Dict params = ctf->to_dict();
4613 ctf2d = Util::ctf_img_real(winsize , winsize, 1, params["defocus"], params["apix"],params["voltage"], params["cs"], \
4614 params["ampcont"], params["bfactor"], params["dfdiff"], params["dfang"], 1);
4615 params.clear();
4616 if(ctf) {delete ctf; ctf=0;}
4617 //ctf2d = ctf_store_real::get_ctf_real(); //This is in 2D projection plane
4618 int nx=ctf2d->get_xsize(),ny=ctf2d->get_ysize(),nz=ctf2d->get_zsize();
4619 float *ctf2d_ptr = ctf2d->get_data();
4620
4621 size_t size = (size_t)nx*ny*nz;
4622 if (!ctf_applied) {
4623 for (int i = 0; i < size; ++i) padfft->cmplx(i) *= ctf2d_ptr[i]; // Multiply padfft by CTF
4624 }
4625
4626 for (int i = 0; i < size; ++i) ctf2d_ptr[i] *= ctf2d_ptr[i]; // Squared 2D CTF
4627 } else {
4628 int nx=padfft->get_xsize(),ny=padfft->get_ysize(),nz=padfft->get_zsize();
4629 ctf2d = new EMData();
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;
4634 }
4635
4636 vector<float> bckgnoise;
4637 bckgnoise = slice->get_attr("bckgnoise");
4638
4639 insert_padfft_slice_weighted(padfft, ctf2d, bckgnoise, t, weight);
4640
4641 checked_delete( ctf2d );
4642 checked_delete( padfft );
4643 bckgnoise.clear();
4644
4645 }
4646 return 0;
4647}
4648
4649int nn4_ctfwReconstructor::insert_padfft_slice_weighted( EMData* padfft, EMData* ctf2d2, vector<float> bckgnoise, const Transform& t, float weight )
4650{
4651 Assert( padfft != NULL );
4652
4653 vector<float> abc_list;
4654 int abc_list_len = 0;
4655 if (m_volume->has_attr("smear")) {
4656 abc_list = m_volume->get_attr("smear");
4657 abc_list_len = abc_list.size();
4658 }
4659
4660 vector<Transform> tsym = t.get_sym_proj(m_symmetry);
4661 for (unsigned int isym=0; isym < tsym.size(); isym++) {
4662 if (abc_list_len == 0)
4663 m_volume->nn_ctfw(m_wptr, padfft, ctf2d2, m_npad, bckgnoise, tsym[isym], weight);
4664 else
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]);
4667 }
4668 return 0;
4669}
4670
4671
4673{
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);
4677 if(m_volume->is_fftodd()) m_volume->symplane0_odd(m_wptr);
4678 else m_volume->symplane0_ctf(m_wptr);
4679 bool do_invert = false;
4680 bool refvol_present = (*m_refvol)(0) > 0.0f ;
4681 /*
4682 int box = 7;
4683 int vol = box*box*box;
4684 int kc = (box-1)/2;
4685 vector< float > pow_a( 3*kc+1, 1.0 );
4686 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
4687 pow_a[3*kc]=0.0;
4688
4689
4690 float max = max3d( kc, pow_a );
4691 float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
4692 float osnr = 1.0f/m_snr;
4693 */
4694
4695
4696 int ix,iy,iz;
4697 // refvol carries fsc
4698 int limitres = m_vnyc-1;
4699 if( refvol_present ) { // If fsc is set to zero, it will be straightforward reconstruction with snr = 1
4700 for (ix = 0; ix < m_vnyc; ix++) {
4701 //cout<<" fsc "<< m_vnyc-ix-1 <<" "<<m_vnyc<<" "<<(*m_refvol)(m_vnyc-ix-1)<<endl;
4702 if( (*m_refvol)(m_vnyc-ix-1) == 0.0f ) limitres = m_vnyc-ix-2;
4703 }
4704 } else {
4705//cout<<" limitres "<<limitres<<endl;
4706
4707 vector<float> count(m_vnyc+1, 0.0f);
4708 vector<float> sigma2(m_vnyc+1, 0.0f);
4709
4710 // compute sigma2
4711 for (iz = 1; iz <= m_vnzp; iz++) {
4712 int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
4713 float argz = float(izp*izp);
4714 for (iy = 1; iy <= m_vnyp; iy++) {
4715 int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
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))) { //Skip Friedel related values
4719 float r = std::sqrt(argy + float(ix*ix));
4720 int ir = int(r);
4721 if (ir <= limitres) {
4722 float frac = r - float(ir);
4723 float qres = 1.0f - frac;
4724 float temp = (*m_wptr)(ix,iy,iz);
4725 //cout<<" WEIGHTS "<<jx<<" "<<jy<<" "<<ir<<" "<<temp<<" "<<frac<<endl;
4726 //cout<<" WEIGHTS "<<ix<<" "<<iy-1<<" "<<iz-1<<" "<<temp<<" "<<endl;
4727 sigma2[ir] += temp*qres;
4728 sigma2[ir+1] += temp*frac;
4729 count[ir] += qres;
4730 count[ir+1] += frac;
4731 }
4732 }
4733 }
4734 }
4735 }
4736 for (ix = 0; ix <= m_vnyc; ix++) {
4737 if( count[ix] > 0.0f ) (*m_refvol)(ix) = sigma2[ix]/count[ix];
4738 //cout<<" sigma2 "<< ix <<" "<<sigma2[ix]<<endl;
4739 }
4740 /*
4741 # This part will have to be done on python level.
4742 float fudge = m_refvol->get_attr("fudge");
4743 // now counter will serve to keep fsc-derived stuff
4744 for (ix = 0; ix <= limitres; ix++) count[ix] = fudge * sigma2[ix] * (1.0f - (*m_refvol)(ix))/(*m_refvol)(ix); //fudge?
4745 count[limitres+1] = count[limitres];
4746 //for (ix = 0; ix <= limitres+1; ix++) cout<<" tau2 "<< ix <<" "<<count[ix]<<endl;
4747 */
4748
4749 }
4750
4751
4752 // normalize
4753 float osnr = 0.0f;
4754 for (iz = 1; iz <= m_vnzp; iz++) {
4755 int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
4756 float argz = float(izp*izp);
4757 for (iy = 1; iy <= m_vnyp; iy++) {
4758 int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
4759 float argy = argz + float(iyp*iyp);
4760 for (ix = 0; ix <= m_vnxc; ix++) {
4761 float r = std::sqrt(argy + float(ix*ix));
4762 int ir = int(r);
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);
4769 //if(osnr == 0.0f) osnr = 1.0f/(0.001*(*m_wptr)(ix,iy,iz));
4770 //cout<<" "<<iz<<" "<<iy<<" "<<" "<<ix<<" "<<ir<<" "<<(*m_wptr)(ix,iy,iz)<<" "<<osnr<<" "<<(*m_volume)(2*ix,iy,iz)<<" "<<(*m_volume)(2*ix+1,iy,iz)<<endl;
4771 } else osnr = 0.0f;
4772
4773 float tmp = ((*m_wptr)(ix,iy,iz)+osnr);
4774
4775 if(do_invert){
4776 if(tmp>0.0f) {
4777 //cout<<" mvol "<<ix<<" "<<iy<<" "<<iz<<" "<<(*m_volume)(2*ix,iy,iz)<<" "<<(*m_volume)(2*ix+1,iy,iz)<<" "<<tmp<<" "<<osnr<<endl;
4778 (*m_volume)(2*ix,iy,iz) /= tmp;
4779 (*m_volume)(2*ix+1,iy,iz) /= tmp;
4780 } else {
4781 (*m_volume)(2*ix,iy,iz) = 0.0f;
4782 (*m_volume)(2*ix+1,iy,iz) = 0.0f;
4783 }
4784 } else (*m_wptr)(ix,iy,iz) = tmp;
4785 }
4786 } else {
4787 (*m_volume)(2*ix,iy,iz) = 0.0f;
4788 (*m_volume)(2*ix+1,iy,iz) = 0.0f;
4789 }
4790 }
4791 }
4792 }
4793
4794 if(do_invert) {
4795 m_volume->center_origin_fft();
4796 // back fft
4797 m_volume->do_ift_inplace();
4798 int npad = m_volume->get_attr("npad");
4799 m_volume->depad();
4800 if( compensate ) circumftrl( m_volume, npad );
4801 m_volume->set_array_offsets( 0, 0, 0 );
4802 }
4803
4804 return 0;
4805}
4806
4807
4808
4809#ifdef False
4811{
4812 m_volume->set_array_offsets(0, 1, 1);
4813 m_wptr->set_array_offsets(0, 1, 1);
4814 m_volume->symplane0_ctf(m_wptr);
4815
4816 int box = 7;
4817 int vol = box*box*box;
4818 int kc = (box-1)/2;
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);
4821 pow_a[3*kc]=0.0;
4822
4823
4824 float max = max3d( kc, pow_a );
4825 float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
4826 float osnr = 1.0f/m_snr;
4827
4828 // normalize
4829 int ix,iy,iz;
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) {//(*v) should be treated as complex!!
4834 float tmp=0.0f;
4835 if( m_varsnr ) {
4836 int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
4837 int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
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);//*m_sign;
4840 } else {
4841 tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+osnr);//*m_sign;
4842 }
4843
4844 if( m_weighting == ESTIMATE ) {
4845 int cx = ix;
4846 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
4847 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
4848 float sum = 0.0;
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;
4858 if( nbrcx < 0 ) {
4859 nbrcx = -nbrcx;
4860 nbrcy = -nbrcy;
4861 nbrcz = -nbrcz;
4862 }
4863
4864 int nbrix = nbrcx;
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];
4870 // if(ix%20==0 && iy%20==0 && iz%20==0)
4871 // std::cout << boost::format( "%4d %4d %4d %4d %10.3f" ) % nbrix % nbriy % nbriz % c % sum << std::endl;
4872 }
4873 }
4874 }
4875 }
4876 float wght = 1.0f / ( 1.0f - alpha * sum );
4877/*
4878 if(ix%10==0 && iy%10==0)
4879 {
4880 std::cout << boost::format( "%4d %4d %4d " ) % ix % iy %iz;
4881 std::cout << boost::format( "%10.3f %10.3f %10.3f " ) % tmp % wght % sum;
4882 std:: << boost::format( "%10.3f %10.3e " ) % pow_b[r] % alpha;
4883 std::cout << std::endl;
4884 }
4885 */
4886 tmp = tmp * wght;
4887 }
4888 (*m_volume)(2*ix,iy,iz) *= tmp;
4889 (*m_volume)(2*ix+1,iy,iz) *= tmp;
4890 }
4891 }
4892 }
4893 }
4894
4895 // back fft
4896 m_volume->do_ift_inplace();
4897 int npad = m_volume->get_attr("npad");
4898 m_volume->depad();
4899 circumfnn( m_volume, npad );
4900 m_volume->set_array_offsets( 0, 0, 0 );
4901
4902 return 0;
4903}
4904
4905
4906{
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);
4910
4911 int box = 7;
4912 int vol = box*box*box;
4913 int kc = (box-1)/2;
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);
4916 pow_a[3*kc]=0.0;
4917
4918
4919 float max = max3d( kc, pow_a );
4920 //float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
4921 float osnr = 1.0f/m_snr;
4922
4923 vector<float> sigma2(m_vnyc+1, 0.0f);
4924 vector<float> count(m_vnyc+1, 0.0f);
4925
4926 int ix,iy,iz;
4927 // compute sigma2
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))) { //Skip Friedel related values
4936 float r = std::sqrt(argy + float(ix*ix));
4937 int ir = int(r);
4938 if (ir <= m_vnyc) {
4939 float frac = r - float(ir);
4940 float qres = 1.0f - frac;
4941 float temp = (*m_wptr)(ix,iy,iz);
4942 //cout<<" WEIGHTS "<<jx<<" "<<jy<<" "<<ir<<" "<<temp<<" "<<frac<<endl;
4943 //cout<<" WEIGHTS "<<ix<<" "<<iy-1<<" "<<iz-1<<" "<<temp<<" "<<endl;
4944 sigma2[ir] += temp*qres;
4945 sigma2[ir+1] += temp*frac;
4946 count[ir] += qres;
4947 count[ir+1] += frac;
4948 }
4949 }
4950 }
4951 }
4952 }
4953 for (ix = 0; ix <= m_vnyc+1; ix++) {
4954 if( sigma2[ix] > 0.0f ) sigma2[ix] = count[ix]/sigma2[ix];
4955 //cout<<" 1/sigma2 "<< ix <<" "<<sigma2[ix]<<endl;
4956 }
4957 // now counter will serve to keep fsc-derived stuff
4958 // refvol carries fsc
4959 float fudge = m_refvol->get_attr("fudge");
4960 for (ix = 0; ix <= m_vnyc+1; ix++)
4961 count[ix] = Util::get_max(0.0f, Util::get_min( 0.999f, (*m_refvol)(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]; //fudge?
4965 }
4966//for (ix = 0; ix <= m_vnyc+1; ix++) cout<<" tau2 "<< ix <<" "<<count[ix]<<" m_wptr "<<(*m_wptr)(ix,1,1)<<endl;
4967