Fourier space 3D reconstruction The Fourier reconstructor is designed to work in an iterative fashion, where similarity ("quality") metrics are used to determine if a slice should be inserted into the 3D in each subsequent iteration. More...

#include <reconstructor.h>

Inheritance diagram for EMAN::FourierReconstructor:

[legend]

Collaboration diagram for EMAN::FourierReconstructor:

[legend]

Public Member Functions
	FourierReconstructor ()
	Default constructor calls load_default_settings() More...

virtual	~FourierReconstructor ()
	Deconstructor calls free_memory() More...

virtual void	setup ()
	Setup the Fourier reconstructor. More...

virtual void	setup_seed (EMData *seed, float seed_weight)
	Initialize the reconstructor with a seed volume. More...

virtual void	setup_seedandweights (EMData seed, EMData weight)
	Initialize the reconstructor with a seed volume, as above. More...

virtual EMData *	preprocess_slice (const EMData *const slice, const Transform &t=Transform())
	Preprocess the slice prior to insertion into the 3D volume this Fourier tranforms the slice and make sure all the pixels are in the right position it always returns a copy of the provided slice, so it should be deleted by someone eventually. More...

virtual int	insert_slice (const EMData *const slice, const Transform &euler, const float weight)
	Insert a slice into a 3D volume, in a given orientation. More...

virtual EMData *	projection (const Transform &euler, int ret_fourier)
	Generates a projection by extracting a slice in Fourier space. More...

virtual int	determine_slice_agreement (EMData *slice, const Transform &euler, const float weight=1.0, bool sub=true)
	Compares a slice to the current reconstruction volume and computes a normalization factor and quality. More...

virtual EMData *	finish (bool doift=true)
	Get the reconstructed volume Normally will return the volume in real-space with the requested size. More...

virtual void	clear ()
	clear the volume and tmp_data for use in Monte Carlo reconstructions More...

virtual string	get_name () const
	Get the unique name of the reconstructor. More...

virtual string	get_desc () const
	Get the one line description of the reconstructor. More...

virtual TypeDict	get_param_types () const
	Get the parameter types of this object. More...

Public Member Functions inherited from EMAN::Reconstructor
	Reconstructor ()

virtual	~Reconstructor ()

int	insert_slice (const EMData *const slice, const Transform &euler)

void	print_params () const
	Print the current parameters to std::out. More...

EMObject &	operator[] (const string &key)

Public Member Functions inherited from EMAN::FactoryBase
	FactoryBase ()

virtual	~FactoryBase ()

Dict	get_params () const
	get a copy of the parameters of this class More...

void	set_params (const Dict &new_params)
	Set new parameters. More...

void	set_param (const string key, const EMObject val)

void	insert_params (const Dict &new_params)
	Insert parameters. More...

Dict	copy_relevant_params (const FactoryBase *const that) const

Public Member Functions inherited from EMAN::ReconstructorVolumeData
	ReconstructorVolumeData ()
	Only constructor All member variables are zeroed. More...

virtual	~ReconstructorVolumeData ()
	Destructor safely frees memory. More...

const EMData *	get_emdata ()
	Get the main image pointer, probably redundant (not used) More...

Static Public Member Functions
static Reconstructor *	NEW ()
	Factory incorporation uses the pointer of this function. More...

Static Public Attributes
static const string	NAME = "fourier"

Protected Member Functions
virtual void	load_default_settings ()
	Load default settings. More...

virtual void	free_memory ()
	Frees the memory owned by this object (but not parent objects) Deletes the FourierPixelInserter3D pointer. More...

virtual void	load_inserter ()
	Load the pixel inserter based on the information in params. More...

virtual void	do_insert_slice_work (const EMData *const input_slice, const Transform &euler, const float weight, const bool corners=false)
	A function to perform the nuts and bolts of inserting an image slice. More...

virtual void	do_compare_slice_work (EMData *input_slice, const Transform &euler, float weight)
	A function to perform the nuts and bolts of comparing an image slice. More...

virtual bool	pixel_at (const float &xx, const float &yy, const float &zz, float *dt)
	This is a mode-2 pixel extractor. More...

Protected Member Functions inherited from EMAN::ReconstructorVolumeData
void	free_memory ()
	Free allocated memorys The inherited class may have allocated image of tmp_data In either case you can safely call this function to delete either of those pointers, even if they bdb:refine_03::threed_00are NULL. More...

virtual void	normalize_threed (const bool sqrt_damp=false, const bool wiener=false)
	Normalize on the assumption that image is a Fourier volume and that tmp_data is a volume of weights corresponding in size to this Fourier volume. More...

virtual void	zero_memory ()
	Sends the pixels in tmp_data and image to zero Convenience only. More...

Protected Attributes
FourierPixelInserter3D *	inserter
	A pixel inserter pointer which inserts pixels into the 3D volume using one of a variety of insertion methods. More...

Protected Attributes inherited from EMAN::FactoryBase
Dict	params
	This is the dictionary the stores the parameters of the object. More...

Protected Attributes inherited from EMAN::ReconstructorVolumeData
EMData *	image
	Inheriting class allocates this, probably in setup(). More...

EMData *	tmp_data
	Inheriting class may allocate this, probably in setup() More...

int	nx

int	nx2

int	ny

int	ny2

int	nz

int	nz2

int	subnx

int	subny

int	subnz

int	subx0

int	suby0

int	subz0

Private Member Functions
	FourierReconstructor (const FourierReconstructor &that)
	Disallow copy construction. More...

FourierReconstructor &	operator= (const FourierReconstructor &)
	Disallow assignment. More...

Detailed Description

Fourier space 3D reconstruction The Fourier reconstructor is designed to work in an iterative fashion, where similarity ("quality") metrics are used to determine if a slice should be inserted into the 3D in each subsequent iteration.

The client creates a Fourier reconstructor to insert real images into a 3D volume. The return image is a real space image

This reconstructor is based on EMAN1's Fourier reconstructor with a handful of modifications including

- Fourier ring correlation (FRC) as opposed to the mean phase residual is used to estimate slice quality. The FRC of the slice in the 3D volume is determined - but the slice is removed from the 3D volume before doing this so the score reflects the extent to which the slice agrees with the contribution of the other slices in the 3D volume. The FRC is converted to SNR using the relationship described by Penczek ( Three-dimensional spectral signal to noise ratio for a class of reconstruction algorithms, JSB 2002 138 (24-46) ) FRC = S/(sqrt(S+N1)sqrt(S+N2)) Where N1 is the noise in the slice of the 3D volume and N2 is the noise in the image slice being inserted. We make the assumption that the noise in the 3D volume is 0 (N1=0) to get FRC^2 = SNR/(1+SNR) which gives a spectral SNR plot - we then divide each SNR value by the number of particles in the class average (seeing as SNR should scale linearly with the number of particles) to get the estimated SNR per contributing particle in this class average. If the particles that have been averaged are not homogenous this score should be low etc. The scaled SNR curve is then converted back to a FRC curve and integrated. This integral is the similarity metric, and depends on how far information extends to in Fourier space - typical values range from 0.05 to 0.2, but can vary substantially depending on the data.

2 - Uses half of the memory used by EMAN1's equivalent reconstruction algorithm

Fourier reconstructor usage

Reconstructor* r = Factory<Reconstructor>::get("fourier", params);
r->setup();
for k in 0:num_iterations-1
        // First do a round of slice quality (metric) determination - only possible if a 3D volume has
        // already been generated (k>0)
        if ( k  > 0 )
                // Determine the agreement of the slices with the previous reconstructed volume (in memory)
                for i in 0:num_slices-1
                        r->determine_slice_agreement(image[i], image[i].euler_orientation);
 
        // Insert the slices into the 3D volume
        // Will decide not to insert the slice if the its "quality" is not good enough
        for i in 0:num_slices-1
                int failure = r->insert_slice(image[i], image[i].euler_orientation);
                if ( failure ) cout << "Slice was not inserted due to poor quality" << endl;
 
// Get the resulting volume
EMData* result = r->finish();
result->write_image("threed.mrc");

Definition at line 388 of file reconstructor.h.

Constructor & Destructor Documentation

◆ FourierReconstructor() [1/2]

EMAN::FourierReconstructor::FourierReconstructor ( )

inline

Default constructor calls load_default_settings()

Definition at line 394 of file reconstructor.h.

394{ load_default_settings(); }

EMAN::FourierReconstructor::load_default_settings

virtual void load_default_settings()

Load default settings.

Definition: reconstructor.cpp:686

References load_default_settings().

Referenced by NEW().

◆ ~FourierReconstructor()

virtual EMAN::FourierReconstructor::~FourierReconstructor ( )

inlinevirtual

Deconstructor calls free_memory()

Definition at line 399 of file reconstructor.h.

399{ free_memory(); }

EMAN::FourierReconstructor::free_memory

virtual void free_memory()

Frees the memory owned by this object (but not parent objects) Deletes the FourierPixelInserter3D poi...

Definition: reconstructor.cpp:693

References free_memory().

◆ FourierReconstructor() [2/2]

EMAN::FourierReconstructor::FourierReconstructor ( const FourierReconstructor & that )

private

Disallow copy construction.

Member Function Documentation

◆ clear()

void FourierReconstructor::clear ( )

virtual

clear the volume and tmp_data for use in Monte Carlo reconstructions

Reimplemented from EMAN::Reconstructor.

Definition at line 999 of file reconstructor.cpp.

{
        bool zeroimage = true;
        bool zerotmpimg = true;
        
#ifdef EMAN2_USING_CUDA
        if(EMData::usecuda == 1) {
                if(image->getcudarwdata()) {
                        to_zero_cuda(image->getcudarwdata(),image->get_xsize(),image->get_ysize(),image->get_zsize());
                        zeroimage = false;
                }
                if(tmp_data->getcudarwdata()) {
                        //to_zero_cuda(tmp_data->getcudarwdata(),image->get_xsize(),image->get_ysize(),image->get_zsize());
                        zerotmpimg = false;
                }
        }
#endif
 
        if(zeroimage) image->to_zero();
        if(zerotmpimg) tmp_data->to_zero();
        
}

References EMAN::ReconstructorVolumeData::image, EMAN::ReconstructorVolumeData::tmp_data, and to_zero_cuda().

◆ determine_slice_agreement()

int FourierReconstructor::determine_slice_agreement	(	EMData *	slice,
		const Transform &	euler,
		const float	weight = `1.0`,
		bool	sub = `true`
	)

virtual

Compares a slice to the current reconstruction volume and computes a normalization factor and quality.

Normalization and quality are returned via attributes set in the passed slice. You may freely mix calls to determine_slice_agreement with calls to insert_slice, but note that determine_slice_agreement can only use information from slices that have already been inserted. Attributes set in the slice are: reconstruct_norm the relative normalization factor which should be applied before inserting the slice reconstruct_qual a scaled quality factor (larger better) for this slice as compared to the existing reconstruction reconstruct_absqual the absolute (not scaled based on weight) quality factor comparing this slice to the existing reconstruction reconstruct_weight the summed weights from all voxels intersecting with the inserted slice, larger -> more overlap with other slices

Parameters

input_slice	The EMData slice to be compared
euler	The orientation of the slice as a Transform object
weight	A weighting factor for this slice, generally the number of particles in a class-average. May be ignored by some reconstructors
sub	Flag indicating whether to subtract the slice from the volume before comparing. May be ignored by some reconstructors

Returns: 0 if OK. 1 if error.

Exceptions

NullPointerException	if the input EMData pointer is null
ImageFormatException	if the image is complex as opposed to real

Reimplemented from EMAN::Reconstructor.

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 1212 of file reconstructor.cpp.

{
        // Are these exceptions really necessary? (d.woolford)
        if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
 
#ifdef EMAN2_USING_CUDA
        if(EMData::usecuda == 1) {
                if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda(); //copy slice to cuda using the const version
        }
#endif
 
        Transform * rotation;
        rotation = new Transform(arg); // assignment operator
 
        EMData *slice;
        if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
        else slice = preprocess_slice( input_slice, *rotation);
 
 
        // We must use only the rotational component of the transform, scaling, translation and mirroring
        // are not implemented in Fourier space, but are in preprocess_slice
        rotation->set_scale(1.0);
        rotation->set_mirror(false);
        rotation->set_trans(0,0,0);
        if (sub) do_insert_slice_work(slice, *rotation, -weight);
        // Remove the current slice first (not threadsafe, but otherwise performance would be awful)
        
        // Compare
        do_compare_slice_work(slice, *rotation,weight);
 
        input_slice->set_attr("reconstruct_norm",slice->get_attr("reconstruct_norm"));
        input_slice->set_attr("reconstruct_absqual",slice->get_attr("reconstruct_absqual"));
        input_slice->set_attr("reconstruct_absqual_lowres",slice->get_attr("reconstruct_absqual_lowres"));
//      input_slice->set_attr("reconstruct_qual",slice->get_attr("reconstruct_qual"));
        input_slice->set_attr("reconstruct_weight",slice->get_attr("reconstruct_weight"));
 
        // Now put the slice back
        if (sub) do_insert_slice_work(slice, *rotation, weight);
 
        delete rotation;
        delete slice;
 
//      image->update();
        return 0;
 
}

References do_compare_slice_work(), do_insert_slice_work(), NullPointerException, preprocess_slice(), EMAN::Transform::set_mirror(), EMAN::Transform::set_scale(), EMAN::Transform::set_trans(), and sub().

◆ do_compare_slice_work()

void FourierReconstructor::do_compare_slice_work	(	EMData *	input_slice,
		const Transform &	euler,
		float	weight
	)

protectedvirtual

A function to perform the nuts and bolts of comparing an image slice.

Parameters

input_slice	the slice to insert into the 3D volume
euler	a transform storing the slice euler angle

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 1259 of file reconstructor.cpp.

{
 
        float dt[3];    // This stores the complex and weight from the volume
        float dt2[2];   // This stores the local image complex
        float *dat = input_slice->get_data();
        vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
 
        float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
        float iny=(float)(input_slice->get_ysize());
        float apix=input_slice->get_attr("apix_x");
        float rmax=iny*apix/12.0;               // radius at 12 A, use as cutoff for low res insertion quality
        if (rmax>iny/2.0-1.0) rmax=iny/2.0-1.0;         // in case of large A/pix
        if (rmax<5.0) rmax=5.0;         // can't imagine this will happen
        
        double dotlow=0;        // low resolution dot product
        double dlpow=0;
        double dlpow2=0;
        
        double dot=0;           // summed pixel*weight dot product
        double vweight=0;               // sum of weights
        double power=0;         // sum of inten*weight from volume
        double power2=0;                // sum of inten*weight from image
//      double amp=0,amp2=0;    // used for normalization, weighted amplitude averages under specific conditions
        bool use_cpu = true;
        
        int nval=0;
 
#ifdef EMAN2_USING_CUDA
        if(EMData::usecuda == 1) {
                if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda();
                if(!image->getcudarwdata()){
                        image->copy_to_cuda();
                        tmp_data->copy_to_cuda();
                }
                for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
                        Transform t3d = arg*(*it);
                        float * m = new float[12];
                        t3d.copy_matrix_into_array(m);
                        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);
                        dot = stats.x;
                        vweight = stats.y;
                        power = stats.z;
                        power2 = stats.w;
                        //cout << "CUDA stats " << stats.x << " " << stats.y << " " << stats.z << " " << stats.w << endl;
                        use_cpu = false;
                }
        }
#endif
        if(use_cpu) {
                for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
                        Transform t3d = arg*(*it);
                        for (int y = -iny/2; y < iny/2; y++) {
                                for (int x = 0; x <=  inx/2; x++) {
                                        if (x==0 && y==0) continue;             // We don't want to use the Fourier origin
 
                                        float rx = (float) x/(inx-2);   // coords relative to Nyquist=.5
                                        float ry = (float) y/iny;
 
//                                      if ((rx * rx + Util::square(ry - max_input_dim "xform.projection"/ 2)) > rl)
//                                      continue;
 
                                        Vec3f coord(rx,ry,0);
                                        coord = coord*t3d; // transpose multiplication
                                        float xx = coord[0]; // transformed coordinates in terms of Nyquist
                                        float yy = coord[1];
                                        float zz = coord[2];
 
 
                                        if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
 
                                        // Map back to actual pixel coordinates in output volume
                                        xx=xx*(nx-2);
                                        yy=yy*ny;
                                        zz=zz*nz;
 
 
                                        int idx = (int)(x * 2 + inx*(y<0?iny+y:y));
                                        dt2[0] = dat[idx];
                                        dt2[1] = dat[idx+1];
 
                                        // value returned indirectly in dt
                                        if (!pixel_at(xx,yy,zz,dt) || dt[2]==0) continue;
 
//                                      printf("%f\t%f\t%f\t%f\t%f\n",dt[0],dt[1],dt[2],dt2[0],dt2[1]);
                                        float r=(float)Util::hypot_fast(x,y);
                                        if (r>3 && r<rmax) {
                                                dotlow+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
                                                dlpow+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
                                                dlpow2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
                                        }
                                        
                                        
//                                      if (r<6) continue; 
                                        vweight+=dt[2];
                                        dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
                                        power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
                                        power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
//                                      printf("%d\t%d\t%f\t%f\t%f\n",x,y,dt[0],dt[1],dt[2]);
//                                      if (r>2 && dt[2]>=0.5) {
//                                              amp+=hypot(dt[0]*dt[0],dt[1]*dt[1])*dt[2];
//                                              amp2+=hypot(dt2[0]*dt2[0],dt2[1]*dt2[1])*dt[2];
//                                      }
                                        nval++;
                                }
                        }
                        //cout << dot << " " << vweight << " " << power << " " << power2 << endl;
                }
        }
        
        dot/=sqrt(power*power2);                // normalize the dot product
        if (dlpow*dlpow2>0) dotlow/=sqrt(dlpow*dlpow2);
//      input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)));
        input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)));
//      input_slice->set_attr("reconstruct_norm",(float)(amp2<=0?1.0:amp/amp2));
        input_slice->set_attr("reconstruct_absqual",(float)dot);
        input_slice->set_attr("reconstruct_absqual_lowres",(float)dotlow);
        float rw=weight<=0?1.0f:1.0f/weight;
        input_slice->set_attr("reconstruct_qual",(float)(dot*rw/((rw-1.0)*dot+1.0)));   // here weight is a proxy for SNR
        input_slice->set_attr("reconstruct_weight",(float)(vweight/(float)(nval)));
//      printf("** %g\t%g\t%g\t%g ##\n",dot,vweight,power,power2);
        //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)));
}

References EMAN::Transform::copy_matrix_into_array(), determine_slice_agreement_cuda(), EMAN::dot(), EMAN::Symmetry3D::get_symmetries(), EMAN::Util::hypot_fast(), EMAN::ReconstructorVolumeData::image, EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, pixel_at(), power(), sqrt(), EMAN::ReconstructorVolumeData::tmp_data, x, and y.

Referenced by determine_slice_agreement().

◆ do_insert_slice_work()

void FourierReconstructor::do_insert_slice_work	(	const EMData *const	input_slice,
		const Transform &	euler,
		const float	weight,
		const bool	corners = `false`
	)

protectedvirtual

A function to perform the nuts and bolts of inserting an image slice.

Parameters

input_slice	the slice to insert into the 3D volume
euler	a transform storing the slice euler angle
weight	weighting factor for this slice (usually number of particles in a class-average)

Definition at line 1122 of file reconstructor.cpp.

{
        // Reload the inserter if the mode has changed
//      string mode = (string) params["mode"];
//      if ( mode != inserter->get_name() )     load_inserter();
 
//      int y_in = input_slice->get_ysize();
//      int x_in = input_slice->get_xsize();
//      // Adjust the dimensions to account for odd and even ffts
//      if (input_slice->is_fftodd()) x_in -= 1;
//      else x_in -= 2;
 
        vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
 
        float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
        float iny=(float)(input_slice->get_ysize());
        
        if (abs(inx-iny)>2 && weight<0) printf("WARNING: Fourier Reconstruction failure. SSNR flag set with asymmetric dimensions on input image\n");
        
#ifdef EMAN2_USING_CUDA
        if(EMData::usecuda == 1) {
                if(!image->getcudarwdata()){
                        image->copy_to_cuda();
                        tmp_data->copy_to_cuda();
                }
                float * m = new float[12];
                for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
                        Transform t3d = arg*(*it);
                        t3d.copy_matrix_into_array(m);
                        //cout << "using CUDA " << image->getcudarwdata() << endl;
                        insert_slice_cuda(m,input_slice->getcudarwdata(),image->getcudarwdata(),tmp_data->getcudarwdata(),inx,iny,image->get_xsize(),image->get_ysize(),image->get_zsize(), weight);
                }
                delete m;
                return;
        }
#endif
        vector<float> ssnr;
        float sscale = 1.0f;
        if (weight<0) {
                ssnr=input_slice->get_attr("class_ssnr");
                sscale=2.0*(ssnr.size()-1)/iny;
        }
        
        float rweight=weight;
        for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
                Transform t3d = arg*(*it);
                for (int y = -iny/2; y < iny/2; y++) {
                        for (int x = 0; x < inx/2; x++) {
 
                                float rx = (float) x/(inx-2.0f);        // coords relative to Nyquist=.5
                                float ry = (float) y/iny;
                                if (corners) {
                                        if (weight<0) {
                                                int r=Util::hypot_fast_int(x,y);
                                                rweight=Util::get_max(0.0f,ssnr[int(r*sscale)]);
                                        }
                                }
                                else {
                                        int r=Util::hypot_fast_int(x,y);
                                        if (r>iny/2 && abs(inx-iny)<3) continue;        // no filling in Fourier corners...
 
                                        if (weight<0) rweight=Util::get_max(0.0f,ssnr[int(r*sscale)]);
                                }
//                              printf("%d\t%f\n",int(r*sscale),rweight);
 
                                Vec3f coord(rx,ry,0);
                                coord = coord*t3d; // transpose multiplication
                                float xx = coord[0]; // transformed coordinates in terms of Nyquist
                                float yy = coord[1];
                                float zz = coord[2];
 
                                // Map back to real pixel coordinates in output volume
                                xx=xx*(nx-2);
                                yy=yy*ny;
                                zz=zz*nz;
                                
//                              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",
//                                      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);
//                              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",
//                                      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));
 
                                //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);
//                              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);
//                              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);
                                inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),rweight);
                        }
                }
        }
}

References EMAN::Transform::copy_matrix_into_array(), EMAN::Util::get_max(), EMAN::Symmetry3D::get_symmetries(), EMAN::Util::hypot_fast_int(), EMAN::ReconstructorVolumeData::image, EMAN::FourierPixelInserter3D::insert_pixel(), insert_slice_cuda(), inserter, EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, EMAN::ReconstructorVolumeData::tmp_data, x, and y.

Referenced by determine_slice_agreement(), and insert_slice().

◆ finish()

EMData * FourierReconstructor::finish ( bool doift = true )

virtual

Get the reconstructed volume Normally will return the volume in real-space with the requested size.

The calling application is responsible for removing any padding.

Parameters

doift A flag indicating whether the returned object should be guaranteed to be in real-space (true) or should be left in whatever space the reconstructor generated

Returns: The real space reconstructed volume

Reimplemented from EMAN::Reconstructor.

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 1552 of file reconstructor.cpp.

{
//      float *norm = tmp_data->get_data();
//      float *rdata = image->get_data();
#ifdef EMAN2_USING_CUDA
        if(EMData::usecuda == 1 && image->getcudarwdata()){
                cout << "copy back from CUDA" << endl;
                image->copy_from_device();
                tmp_data->copy_from_device();
        }
#endif
        
        bool sqrtnorm=params.set_default("sqrtnorm",false);
        normalize_threed(sqrtnorm);
//      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));
        
// This compares single precision sum to double precision sum near the origin
#ifdef RECONDEBUG
        for (int k=0; k<5; k++) {
                for (int j=0; j<5; j++) {
                        for (int i=0; i<5; i++) {
                                int idx=i*2+j*10+k*50;
                                std::complex <double> a(ddata[idx],ddata[idx+1]);
                                double b=dnorm[idx];
                                a/=b;
                                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));
                        }
                }
        }
#endif
        
//      tmp_data->write_image("density.mrc");
 
        // we may as well delete the tmp data now... it saves memory and the calling program might
        // need memory after it gets the return volume.
        // If this delete didn't happen now, it would happen when the deconstructor was called --david
        // no longer a good idea with the new iterative scheme -- steve
//      if ( tmp_data != 0 )
//      {
//              delete tmp_data;
//              tmp_data = 0;
//      }
 
/*      image->process_inplace("xform.fourierorigin.tocorner");*/
 
        if (doift) {
                image->do_ift_inplace();
                image->depad();
                image->process_inplace("xform.phaseorigin.tocenter");
        }
        // If the image was padded it should be the original size, as the client would expect
        //  I blocked the rest, it is almost certainly incorrect  PAP 07/31/08
        // No, it's not incorrect. You are wrong. You have the meaning of nx mixed up. DSAW 09/23/cd
        // This should now be handled in the calling program --steve 11/03/09
//      bool is_fftodd = (nx % 2 == 1);
//      if ( (nx-2*(!is_fftodd)) != output_x || ny != output_y || nz != output_z )
//      {
//              FloatPoint origin( (nx-output_x)/2, (ny-output_y)/2, (nz-output_z)/2 );
//              FloatSize region_size( output_x, output_y, output_z);
//              Region clip_region( origin, region_size );
//              image->clip_inplace( clip_region );
//      }
 
        // Should be an "if (verbose)" here or something
        //print_stats(quality_scores);
 
        image->update();
        
        
        if (params.has_key("normout")) {
//              if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
                EMData *normout=(EMData*) params["normout"];
                normout->set_data(tmp_data->copy()->get_data());
                normout->update();
        }
        
        if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) {
                if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
                tmp_data->write_image((const char *)params["savenorm"]);
        }
 
        delete tmp_data;
        tmp_data=0;
        //Since we give up the ownership of the pointer to-be-returned, it's caller's responsibility to delete the returned image.
        //So we wrap this function with return_value_policy< manage_new_object >() in libpyReconstructor2.cpp to hand over ownership to Python.
        EMData *ret=image;
        image=0;
        
        return ret;
}

References EMAN::Dict::has_key(), EMAN::ReconstructorVolumeData::image, EMAN::ReconstructorVolumeData::normalize_threed(), EMAN::FactoryBase::params, EMAN::Dict::set_default(), and EMAN::ReconstructorVolumeData::tmp_data.

◆ free_memory()

void FourierReconstructor::free_memory ( )

protectedvirtual

Frees the memory owned by this object (but not parent objects) Deletes the FourierPixelInserter3D pointer.

Definition at line 693 of file reconstructor.cpp.

{
        if (image) { delete image; image=0; }
        if (tmp_data) { delete tmp_data; tmp_data=0; }
        if ( inserter != 0 )
        {
                delete inserter;
                inserter = 0;
        }
}

References EMAN::ReconstructorVolumeData::image, inserter, and EMAN::ReconstructorVolumeData::tmp_data.

Referenced by ~FourierReconstructor().

◆ get_desc()

virtual string EMAN::FourierReconstructor::get_desc ( ) const

inlinevirtual

Get the one line description of the reconstructor.

Implements EMAN::FactoryBase.

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 485 of file reconstructor.h.

                {
                        return "Reconstruction via direct Fourier methods using one of a variety of different kernels, most of which are Gaussian based";
                }

◆ get_name()

virtual string EMAN::FourierReconstructor::get_name ( ) const

inlinevirtual

Get the unique name of the reconstructor.

Implements EMAN::FactoryBase.

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 478 of file reconstructor.h.

                {
                        return NAME;
                }

References NAME.

◆ get_param_types()

virtual TypeDict EMAN::FourierReconstructor::get_param_types ( ) const

inlinevirtual

Get the parameter types of this object.

Returns: a TypeDict detailing all of the acceptable (and necessary) parameters

Implements EMAN::FactoryBase.

Definition at line 501 of file reconstructor.h.

                {
                        TypeDict d;
                        d.put("size", EMObject::INTARRAY, "Required. The dimensions of the real-space output volume, including any padding (must be handled by the calling application). Assumed that apix x/y/z identical.");
                        d.put("sym", EMObject::STRING, "Required. The symmetry of the reconstructed volume, c?, d?, oct, tet, icos, h?. Default is c1, ie - an asymmetric object");
                        d.put("mode", EMObject::STRING, "Optional. Fourier pixel insertion mode name (nearest_neighbor, gauss_2, gauss_3, gauss_5, gauss_5_slow, gypergeom_5, experimental) gauss_2 is the default.");
                        d.put("corners", EMObject::BOOL, "Optional. If not set, then reconstruction will cover a spherical volume with a radius of nx/2. If set, the full Fourier volume will be reconstructed, but will be ~2x slower.");
                        d.put("sqrtnorm", EMObject::BOOL, "Optional. When normalizing, additionally divides by the sqrt of the normalization factor to damp exaggerated features. Is this justifyable ? No idea (yet). Default is false.");
                        d.put("usessnr", EMObject::BOOL, "Optional. Looks for and uses the class_ssnr header parameter from each slice to weight each voxel during insertion to the reconstruction.");
                        d.put("verbose", EMObject::BOOL, "Optional. Toggles writing useful information to standard out. Default is false.");
                        d.put("quiet", EMObject::BOOL, "Optional. If false, print verbose information.");
                        d.put("subvolume",EMObject::INTARRAY, "Optional. (xorigin,yorigin,zorigin,xsize,ysize,zsize) all in Fourier pixels. Useful for parallelism.");
                        d.put("savenorm",EMObject::STRING, "Debug. Will cause the normalization volume to be written directly to the specified file when finish() is called.");
                        
                        d.put("normout",EMObject::EMDATA, "Will write the normalization volume to the given EMData object file when finish() is called.");
                        return d;
                }

References EMAN::EMObject::BOOL, EMAN::EMObject::EMDATA, EMAN::EMObject::INTARRAY, EMAN::TypeDict::put(), and EMAN::EMObject::STRING.

◆ insert_slice()

int FourierReconstructor::insert_slice	(	const EMData *const	slice,
		const Transform &	euler,
		const float	weight
	)

virtual

Insert a slice into a 3D volume, in a given orientation.

Returns: 0 if successful, 1 otherwise

Parameters

slice	the image slice to be inserted into the 3D volume
euler	Euler angle of this image slice.
weight	A weighting factor for this slice, generally the number of particles in a class-average. May be ignored by some reconstructors

Returns: 0 if OK. 1 if error.

Exceptions

NullPointerException	if the input EMData pointer is null
ImageFormatException	if the image is complex as opposed to real

Reimplemented from EMAN::Reconstructor.

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 1064 of file reconstructor.cpp.

{
        // Are these exceptions really necessary? (d.woolford)
        if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
 
        
//      image->set_attr("apix_x",input_slice->get_attr("apix_x"));
//      image->set_attr("apix_y",input_slice->get_attr("apix_y"));
//      image->set_attr("apix_z",input_slice->get_attr("apix_z"));
 
#ifdef EMAN2_USING_CUDA
        if(EMData::usecuda == 1) {
                if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda(); //copy slice to cuda using the const version
        }
#endif
 
        bool usessnr=params.set_default("usessnr",false);
        bool corners=params.set_default("corners",false);
        float weight=oweight;
        if (usessnr) {
                if (input_slice->has_attr("class_ssnr")) weight=-1.0;   // negative weight is a flag for using SSNR
                else weight=0;
        }
 
        if (weight==0) return -1;
        
Transform * rotation;
/*      if ( input_slice->has_attr("xform.projection") ) {
                rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
        } else {*/
        rotation = new Transform(arg); // assignment operator
//      }
 
        EMData *slice;
        if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
        else slice = preprocess_slice( input_slice, *rotation);
 
 
        // We must use only the rotational component of the transform, scaling, translation and mirroring
        // are not implemented in Fourier space, but are in preprocess_slice
        rotation->set_scale(1.0);
        rotation->set_mirror(false);
        rotation->set_trans(0,0,0);
 
        // Finally to the pixel wise slice insertion
        //slice->copy_to_cuda();
//      EMData *s2=slice->do_ift();
//      s2->write_image("is.hdf",-1);
        do_insert_slice_work(slice, *rotation, weight, corners);
        
        delete rotation; rotation=0;
        delete slice;
 
//      image->update();
        return 0;
}

References do_insert_slice_work(), NullPointerException, EMAN::FactoryBase::params, preprocess_slice(), EMAN::Dict::set_default(), EMAN::Transform::set_mirror(), EMAN::Transform::set_scale(), and EMAN::Transform::set_trans().

◆ load_default_settings()

void FourierReconstructor::load_default_settings ( )

protectedvirtual

Load default settings.

Definition at line 686 of file reconstructor.cpp.

{
        inserter=0;
        image=0;
        tmp_data=0;
}

References EMAN::ReconstructorVolumeData::image, inserter, and EMAN::ReconstructorVolumeData::tmp_data.

Referenced by FourierReconstructor().

◆ load_inserter()

void FourierReconstructor::load_inserter ( )

protectedvirtual

Load the pixel inserter based on the information in params.

Definition at line 706 of file reconstructor.cpp.

{
//      ss
//      string mode = (string)params["mode"];
        Dict parms;
        parms["data"] = image;
        parms["norm"] = tmp_data->get_data();
        // These aren't necessary because we deal with them before calling the inserter
//      parms["subnx"] = nx;
//      parms["subny"] = ny;
//      parms["subnx"] = nz;
//      parms["subx0"] = x0;
//      parms["suby0"] = y0;
//      parms["subz0"] = z0;
 
        if ( inserter != 0 )
        {
                delete inserter;
        }
 
        inserter = Factory<FourierPixelInserter3D>::get((string)params["mode"], parms);
        inserter->init();
}

References EMAN::Factory< T >::get(), EMAN::ReconstructorVolumeData::image, EMAN::FourierPixelInserter3D::init(), inserter, EMAN::FactoryBase::params, and EMAN::ReconstructorVolumeData::tmp_data.

Referenced by setup(), setup_seed(), and setup_seedandweights().

◆ NEW()

static Reconstructor * EMAN::FourierReconstructor::NEW ( )

inlinestatic

Factory incorporation uses the pointer of this function.

Returns: a Reconstructor pointer to a newly allocated FourierReconstructor

Definition at line 493 of file reconstructor.h.

                {
                        return new FourierReconstructor();
                }

References FourierReconstructor().

◆ operator=()

FourierReconstructor & EMAN::FourierReconstructor::operator= ( const FourierReconstructor & )

private

Disallow assignment.

◆ pixel_at()

bool FourierReconstructor::pixel_at	(	const float &	xx,
		const float &	yy,
		const float &	zz,
		float *	dt
	)

protectedvirtual

This is a mode-2 pixel extractor.

Parameters

xx,yy,zz	voxel coordinates (need not be integers)
dt	float pointer with 3 floats allocated for returned complex value and weight sum

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 1448 of file reconstructor.cpp.

{
        int x0 = (int) Util::round(xx);
        int y0 = (int) Util::round(yy);
        int z0 = (int) Util::round(zz);
 
        std::complex<float> val=image->get_complex_at(x0,y0,z0);
        size_t idx=image->get_complex_index_fast(x0,y0,z0)/2;
        float norm=tmp_data->get_value_at_index(idx);
        dt[0]=val.real()/norm;
        dt[1]=val.imag()/norm;
        dt[2]=norm;
//      printf("%d %d %d    %g %g %g\n",x0,y0,z0,dt[0],dt[1],dt[2]);
        return true;
        
        //Trilinear interpolation version caused prominent interpolation artifacts due to high noise levels
        // between 1/2 Nyquist and Nyquist. While the nearest neighbor interpolation used above still has
        // artifacts, it seems better for this puprose.
        
//      int x0 = (int) floor(xx);
//      int y0 = (int) floor(yy);
//      int z0 = (int) floor(zz);
//      
//      float *rdata=image->get_data();
//      float *norm=tmp_data->get_data();
//      float normsum=0,normsum2=0;
// 
//      dt[0]=dt[1]=dt[2]=0.0;
// 
//      if (nx==subnx) {                        // normal full reconstruction
//              if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
// 
//              
//              
//              int x1=x0+1;
//              int y1=y0+1;
//              int z1=z0+1;
//              if (x0<-nx2) x0=-nx2;
//              if (x1>nx2) x1=nx2;
//              if (y0<-ny2) y0=-ny2;
//              if (y1>ny2) y1=ny2;
//              if (z0<-nz2) z0=-nz2;
//              if (z1>nz2) z1=nz2;
//              
//              size_t idx=0;
//              float r, gg;
//              for (int k = z0 ; k <= z1; k++) {
//                      for (int j = y0 ; j <= y1; j++) {
//                              for (int i = x0; i <= x1; i ++) {
//                                      idx=image->get_complex_index_fast(i,j,k);
// //                                   r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
// //                                   gg = Util::fast_exp(-r / EMConsts::I2G);
//                                      gg=(1.0-fabs(i-xx))*(1.0-fabs(j-yy))*(1.0-fabs(k-zz));  // Change from Gaussian to trilinear, 6/7/20
//                                      
//                                      dt[0]+=gg*rdata[idx];
//                                      dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
//                                      //dt[2]+=norm[idx/2]*gg;
//                                      normsum2+=gg;
//                                      normsum+=gg*norm[idx/2];
//                                      if (i==0) {
//                                              normsum2+=gg;
//                                              normsum+=gg*norm[idx/2];                                
//                                      }
//                              }
//                      }
//              }
//              if (normsum==0) return false;
//              dt[0]/=normsum;
//              dt[1]/=normsum;
//              dt[2]=normsum/normsum2;
// #ifdef DEBUG_POINT
//              if (z0==23 && y0==0 && x0==113) {
//                      std::complex<float> pv = image->get_complex_at(113,0,23);
//                      float pnv = tmp_data->get_value_at(113,0,23);
//                      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]);
//              }
// #endif
// //           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]);
//              return true;
//      } 
//      else {                                  // for subvolumes, not optimized yet
//              size_t idx;
//              float r, gg;
//              for (int k = z0 ; k <= z0 + 1; k++) {
//                      for (int j = y0 ; j <= y0 + 1; j++) {
//                              for (int i = x0; i <= x0 + 1; i ++) {
//                                      r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
//                                      idx=image->get_complex_index(i,j,k,subx0,suby0,subz0,nx,ny,nz);
//                                      gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
//                                      
//                                      dt[0]+=gg*rdata[idx];
//                                      dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
//                                      dt[2]+=norm[idx/2];
//                                      normsum+=gg;
//                              }
//                      }
//              }
//              
//              if (normsum==0)  return false;
//              return true;
//      }
}

References EMAN::ReconstructorVolumeData::image, EMAN::Util::round(), and EMAN::ReconstructorVolumeData::tmp_data.

Referenced by do_compare_slice_work(), and projection().

◆ preprocess_slice()

EMData * FourierReconstructor::preprocess_slice	(	const EMData *const	slice,
		const Transform &	t = `Transform()`
	)

virtual

Preprocess the slice prior to insertion into the 3D volume this Fourier tranforms the slice and make sure all the pixels are in the right position it always returns a copy of the provided slice, so it should be deleted by someone eventually.

Returns: the processed slice

Parameters

slice	the slice to be prepocessed
t	transform to be used for insertion

Exceptions

InvalidValueException when the specified padding value is less than the size of the images

Reimplemented from EMAN::Reconstructor.

Definition at line 1022 of file reconstructor.cpp.

{
#ifdef EMAN2_USING_CUDA
        if(EMData::usecuda == 1) {
                if(!slice->getcudarwdata()) slice->copy_to_cuda(); //copy slice to cuda using the const version
        }
#endif
        // Shift the image pixels so the real space origin is now located at the phase origin (at the bottom left of the image)
        EMData* return_slice = 0;
        Transform tmp(t);
        tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly, this way we only do 2d translation,scaling and mirroring
 
        if (tmp.is_identity()) return_slice=slice->copy();
        else return_slice = slice->process("xform",Dict("transform",&tmp));
 
        return_slice->process_inplace("xform.phaseorigin.tocorner");
 
//      printf("@@@ %d\n",(int)return_slice->get_attr("nx"));
        // Fourier transform the slice
        
#ifdef EMAN2_USING_CUDA
        if(EMData::usecuda == 1 && return_slice->getcudarwdata()) {
                return_slice->do_fft_inplace_cuda(); //a CUDA FFT inplace is quite slow as there is a lot of mem copying.
        }else{
                return_slice->do_fft_inplace();
        }
#else
        return_slice->do_fft_inplace();
#endif
 
//      printf("%d\n",(int)return_slice->get_attr("nx"));
 
        return_slice->mult((float)sqrt(1.0f/(return_slice->get_ysize())*return_slice->get_xsize()));
 
        // Shift the Fourier transform so that it's origin is in the center (bottom) of the image.
//      return_slice->process_inplace("xform.fourierorigin.tocenter");
 
        return_slice->set_attr("reconstruct_preproc",(int)1);
        return return_slice;
}

References EMAN::Transform::is_identity(), EMAN::Transform::set_rotation(), and sqrt().

Referenced by determine_slice_agreement(), EMAN::WienerFourierReconstructor::determine_slice_agreement(), insert_slice(), and EMAN::WienerFourierReconstructor::insert_slice().

◆ projection()

EMData * FourierReconstructor::projection	(	const Transform &	euler,
		int	ret_fourier
	)

virtual

Generates a projection by extracting a slice in Fourier space.

Parameters

euler	The orientation of the slice as a Transform object
ret_fourier	If true, will return the Fourier transform of the projection

Reimplemented from EMAN::Reconstructor.

Definition at line 1383 of file reconstructor.cpp.

                                                                                {
        
        if (subx0!=0 || suby0!=0 || subz0!=0 || subnx!=nx || subny!=ny ||subnz!=nz) 
                throw ImageDimensionException("ERROR: Reconstructor->projection() does not work with subvolumes");
        
        EMData *ret = new EMData(nx,ny,1);
        ret->set_complex(1);
        ret->set_ri(1);
        ret->set_attr("apix_x",image->get_attr("apix_x"));
        ret->set_attr("apix_y",image->get_attr("apix_y"));
        ret->set_attr("apix_z",image->get_attr("apix_z"));
        
        // We must use only the rotational component of the transform, scaling, translation and mirroring
        // are not implemented in Fourier space, but are in preprocess_slice
        Transform rotation(euler);
        rotation.set_scale(1.0);
        rotation.set_mirror(false);
        rotation.set_trans(0,0,0);
        
        float dt[3];    // This stores the complex and weight from the volume
 
//      float apix=input_slice->get_attr("apix_x");
//      float rmax=iny*apix/12.0;               // radius at 12 A, use as cutoff for low res insertion quality
        
        for (int y = -ny/2; y < ny/2; y++) {
                for (int x = 0; x <=  nx/2; x++) {
 
                        float rx = (float) x/(nx-2);    // coords relative to Nyquist=.5
                        float ry = (float) y/ny;
 
 
                        Vec3f coord(rx,ry,0);
                        coord = coord*rotation; // transpose multiplication
                        float xx = coord[0]; // transformed coordinates in terms of Nyquist
                        float yy = coord[1];
                        float zz = coord[2];
 
 
                        if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
 
                        // Map back to actual pixel coordinates in output volume
                        xx=xx*(nx-2);
                        yy=yy*ny;
                        zz=zz*nz;
 
                        if (!pixel_at(xx,yy,zz,dt) || dt[2]<0.005) continue;
//                      if (x==0) printf("%d %d %g %g %g\n",x,y,dt[0],dt[1],dt[2]);
                        if ((x==0 ||x==nx/2) && (y!=0 && y!=-ny/2)) { dt[0]/=2; dt[1]/=2; }
                        ret->set_complex_at(x,y,std::complex<float>(dt[0],dt[1]));      // division by dt[2] already handled in pixel_at
                }
        }
        
        Transform translation;
        translation.set_trans(euler.get_trans());
        ret->process_inplace("xform",Dict("transform",&translation));
        if (!ret_fourier) {
                ret->process_inplace("xform.phaseorigin.tocenter");
                EMData *tmp=ret->do_ift();
                delete ret;
                ret=tmp;
        }
        return ret;
}

◆ setup()

void FourierReconstructor::setup ( )

virtual

Setup the Fourier reconstructor.

Exceptions

InvalidValueException When one of the input parameters is invalid

Implements EMAN::Reconstructor.

Definition at line 730 of file reconstructor.cpp.

{
        // default setting behavior - does not override if the parameter is already set
        params.set_default("mode","gauss_2");
        params.set_default("verbose",(int)0);
        
        vector<int> size=params["size"];
 
        nx = size[0];
        ny = size[1];
        nz = size[2];
        nx2=nx/2-1;
        ny2=ny/2;
        nz2=nz/2;
 
#ifdef RECONDEBUG
        ddata=(double *)malloc(sizeof(double)*250);
        dnorm=(double *)malloc(sizeof(double)*250);
        for (int i=0; i<250; i++) ddata[i]=dnorm[i]=0.0;
#endif
 
        // Adjust nx if for Fourier transform even odd issues
        bool is_fftodd = (nx % 2 == 1);
        // The Fourier transform requires one extra pixel in the x direction,
        // which is two spaces in memory, one each for the complex and the
        // real components
        nx += 2-is_fftodd;
 
        if (params.has_key("subvolume")) {
                vector<int> sub=params["subvolume"];
                subx0=sub[0];
                suby0=sub[1];
                subz0=sub[2];
                subnx=sub[3];
                subny=sub[4];
                subnz=sub[5];
 
                if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
                        throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
 
        }
        else {
                subx0=suby0=subz0=0;
                subnx=nx;
                subny=ny;
                subnz=nz;
        }
 
 
        // Odd dimension support is here atm, but not above.
        if (image) delete image;
        image = new EMData();
        image->set_size(subnx, subny, subnz);
        image->set_complex(true);
        image->set_fftodd(is_fftodd);
        image->set_ri(true);
//      printf("%d %d %d\n\n",subnx,subny,subnz);
        image->to_zero();
 
        if (params.has_key("subvolume")) {
                image->set_attr("subvolume_x0",subx0);
                image->set_attr("subvolume_y0",suby0);
                image->set_attr("subvolume_z0",subz0);
                image->set_attr("subvolume_full_nx",nx);
                image->set_attr("subvolume_full_ny",ny);
                image->set_attr("subvolume_full_nz",nz);
        }
        
        if (tmp_data) delete tmp_data;
        tmp_data = new EMData();
        tmp_data->set_size(subnx/2, subny, subnz);
        tmp_data->to_zero();
        tmp_data->update();
 
        load_inserter();
 
#ifdef RECONDEBUG
        printf("copied\n");
        inserter->ddata=ddata;
        inserter->dnorm=dnorm;
#endif
        
        if ( (bool) params["verbose"] )
        {
                cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
                cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
                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);
        }
}

◆ setup_seed()

void FourierReconstructor::setup_seed	(	EMData *	seed,
		float	seed_weight
	)

virtual

Initialize the reconstructor with a seed volume.

This can be used to provide some 'default' value when there is missing data in Fourier space. The passed 'seed' must be of the appropriate padded size, must be in Fourier space, and the same EMData* object will be returned by finish(), meaning the Reconstructor is implicitly taking ownership of the object. However, in Python, this means the seed may be passed in without copying, as the same EMData will be coming back out at the end. The seed_weight determines how 'strong' the seed volume should be as compared to other inserted slices in Fourier space.

Exceptions

InvalidValueException When one of the input parameters is invalid

Reimplemented from EMAN::Reconstructor.

Definition at line 820 of file reconstructor.cpp.

                                                                    {
        // default setting behavior - does not override if the parameter is already set
        params.set_default("mode","gauss_2");
 
        vector<int> size=params["size"];
 
        nx = size[0];
        ny = size[1];
        nz = size[2];
        nx2=nx/2-1;
        ny2=ny/2;
        nz2=nz/2;
 
 
        // Adjust nx if for Fourier transform even odd issues
        bool is_fftodd = (nx % 2 == 1);
        // The Fourier transform requires one extra pixel in the x direction,
        // which is two spaces in memory, one each for the complex and the
        // real components
        nx += 2-is_fftodd;
 
        if (params.has_key("subvolume")) {
                vector<int> sub=params["subvolume"];
                subx0=sub[0];
                suby0=sub[1];
                subz0=sub[2];
                subnx=sub[3];
                subny=sub[4];
                subnz=sub[5];
 
                if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
                        throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
 
        }
        else {
                subx0=suby0=subz0=0;
                subnx=nx;
                subny=ny;
                subnz=nz;
        }
 
        if (seed->get_xsize()!=subnx || seed->get_ysize()!=subny || seed->get_zsize()!=subnz || !seed->is_complex())
                throw ImageDimensionException("The dimensions of the seed volume do not match the reconstruction size");
 
        // Odd dimension support is here atm, but not above.
        image = seed->copy();
        image->mult(seed_weight);               // The lack of this line was what was causing all of the strange normalization issues
 
        
 
        if (params.has_key("subvolume")) {
                image->set_attr("subvolume_x0",subx0);
                image->set_attr("subvolume_y0",suby0);
                image->set_attr("subvolume_z0",subz0);
                image->set_attr("subvolume_full_nx",nx);
                image->set_attr("subvolume_full_ny",ny);
                image->set_attr("subvolume_full_nz",nz);
        }
 
        if (tmp_data) delete tmp_data;
        tmp_data = new EMData();
        tmp_data->set_size(subnx/2, subny, subnz);
        tmp_data->to_value(seed_weight);
        
        // The values at x=0 and x=nyquist are added twice in the map normally, but the normalization is only added once
        // This is compensated in the final normalization. In the seed volume we accomodate this by reducing the seed
        // alternatively we could double the actual values. Impact is similar, if not 100% identical
        if (subx0==0) {
                for (int k=0; k<subnz; k++) {
                        for (int j=0; j<subny; j++) {
                                tmp_data->set_value_at(0,j,k,seed_weight/2.0);
                        }
                }
                tmp_data->set_value_at(0,subny/2,0,seed_weight);
                tmp_data->set_value_at(0,0,subnz/2,seed_weight);
        }
        if (subnx+subx0==nx) {
                for (int k=0; k<subnz; k++) {
                        for (int j=0; j<subny; j++) {
                                tmp_data->set_value_at(subnx/2-1,j,k,seed_weight/2.0);
                        }
                }
                tmp_data->set_value_at(subnx/2-1,subny/2,0,seed_weight);
                tmp_data->set_value_at(subnx/2-1,0,subnz/2,seed_weight);
        }
        
 
        
        load_inserter();
        
#ifdef DEBUG_POINT
        std::complex<float> pv = image->get_complex_at(113,0,23);
        float pnv = tmp_data->get_value_at(113,0,23);
        printf("seed: %1.4g\t%1.4g\t%1.4g\n",pv.real()/pnv,pv.imag()/pnv,pnv);
#endif
 
 
        if ( (bool) params["quiet"] == false )
        {
                cout << "Seeded direct Fourier inversion";
                cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
                cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
                cout << "You will require approximately " << setprecision(3) << (subnx*subny*subnz*sizeof(float)*1.5)/1000000000.0 << "GB of memory to reconstruct this volume" << endl;
        }
}

◆ setup_seedandweights()

void FourierReconstructor::setup_seedandweights	(	EMData *	seed,
		EMData *	weight
	)

virtual

Initialize the reconstructor with a seed volume, as above.

In this case the initial weight map is also provided explicitly, rather than a single weight value.

Reimplemented from EMAN::Reconstructor.

Definition at line 926 of file reconstructor.cpp.

                                                                                {
        // default setting behavior - does not override if the parameter is already set
        
        // 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)
        params.set_default("mode","gauss_2");
 
        vector<int> size=params["size"];
 
        nx = size[0];
        ny = size[1];
        nz = size[2];
        nx2=nx/2-1;
        ny2=ny/2;
        nz2=nz/2;
 
 
        // Adjust nx if for Fourier transform even odd issues
        bool is_fftodd = (nx % 2 == 1);
        // The Fourier transform requires one extra pixel in the x direction,
        // which is two spaces in memory, one each for the complex and the
        // real components
        nx += 2-is_fftodd;
 
        if (params.has_key("subvolume")) {
                vector<int> sub=params["subvolume"];
                subx0=sub[0];
                suby0=sub[1];
                subz0=sub[2];
                subnx=sub[3];
                subny=sub[4];
                subnz=sub[5];
 
                if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
                        throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
 
        }
        else {
                subx0=suby0=subz0=0;
                subnx=nx;
                subny=ny;
                subnz=nz;
        }
 
        if (seed->get_xsize()!=subnx || seed->get_ysize()!=subny || seed->get_zsize()!=subnz || !seed->is_complex())
                throw ImageDimensionException("The dimensions of the seed volume do not match the reconstruction size");
 
        // Odd dimension support is here atm, but not above.
        image = seed->copy();   // note that if the provided seed has not been 'premultiplied' by the weights, bad things may result!
        
        if (params.has_key("subvolume")) {
                image->set_attr("subvolume_x0",subx0);
                image->set_attr("subvolume_y0",suby0);
                image->set_attr("subvolume_z0",subz0);
                image->set_attr("subvolume_full_nx",nx);
                image->set_attr("subvolume_full_ny",ny);
                image->set_attr("subvolume_full_nz",nz);
        }
 
        if (tmp_data) delete tmp_data;
        tmp_data=seed_weight->copy();
 
        load_inserter();
 
        if ( (bool) params["quiet"] == false )
        {
                cout << "Seeded direct Fourier inversion";
                cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
                cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
                cout << "You will require approximately " << setprecision(3) << (subnx*subny*subnz*sizeof(float)*1.5)/1000000000.0 << "GB of memory to reconstruct this volume" << endl;
        }
}

Member Data Documentation

◆ inserter

FourierPixelInserter3D* EMAN::FourierReconstructor::inserter

protected

A pixel inserter pointer which inserts pixels into the 3D volume using one of a variety of insertion methods.

Definition at line 555 of file reconstructor.h.

Referenced by EMAN::WienerFourierReconstructor::do_insert_slice_work(), do_insert_slice_work(), free_memory(), load_default_settings(), load_inserter(), and setup().

◆ NAME

const string FourierReconstructor::NAME = "fourier"

static

Definition at line 519 of file reconstructor.h.

Referenced by get_name().

The documentation for this class was generated from the following files:

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Protected Attributes

Private Member Functions

Detailed Description

Constructor & Destructor Documentation

◆ FourierReconstructor() [1/2]

◆ ~FourierReconstructor()

◆ FourierReconstructor() [2/2]

Member Function Documentation

◆ clear()

◆ determine_slice_agreement()

◆ do_compare_slice_work()

◆ do_insert_slice_work()

◆ finish()

◆ free_memory()

◆ get_desc()

◆ get_name()

◆ get_param_types()

◆ insert_slice()

◆ load_default_settings()

◆ load_inserter()

◆ NEW()

◆ operator=()

◆ pixel_at()

◆ preprocess_slice()

◆ projection()

◆ setup()

◆ setup_seed()

◆ setup_seedandweights()

Member Data Documentation

◆ inserter

◆ NAME