EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | Private Member Functions
EMAN::FourierReconstructor Class Reference

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:
Inheritance graph
[legend]
Collaboration diagram for EMAN::FourierReconstructor:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 FourierReconstructor ()
 Default constructor calls load_default_settings()
virtual ~FourierReconstructor ()
 Deconstructor calls free_memory()
virtual void setup ()
 Setup the Fourier reconstructor

Exceptions:
InvalidValueExceptionWhen one of the input parameters is invalid.

virtual void setup_seed (EMData *seed, float seed_weight)
 Initialize the reconstructor with a seed volume.
virtual EMDatapreprocess_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.
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.
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.
virtual EMDatafinish (bool doift=true)
 Get the reconstructed volume Normally will return the volume in real-space with the requested size.
virtual void clear ()
 clear the volume and tmp_data for use in Monte Carlo reconstructions
virtual string get_name () const
 Get the unique name of the reconstructor.
virtual string get_desc () const
 Get the one line description of the reconstructor.
virtual TypeDict get_param_types () const
 Get the parameter types of this object.

Static Public Member Functions

static ReconstructorNEW ()
 Factory incorporation uses the pointer of this function.

Static Public Attributes

static const string NAME = "fourier"

Protected Member Functions

virtual void load_default_settings ()
 Load default settings.
virtual void free_memory ()
 Frees the memory owned by this object (but not parent objects) Deletes the FourierPixelInserter3D pointer.
virtual void load_inserter ()
 Load the pixel inserter based on the information in params.
virtual void do_insert_slice_work (const EMData *const input_slice, const Transform &euler, const float weight)
 A function to perform the nuts and bolts of inserting an image slice.
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.
virtual bool pixel_at (const float &xx, const float &yy, const float &zz, float *dt)
 This is a mode-2 pixel extractor.

Protected Attributes

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

Private Member Functions

 FourierReconstructor (const FourierReconstructor &that)
 Disallow copy construction.
FourierReconstructoroperator= (const FourierReconstructor &)
 Disallow assignment.

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 1. - 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

        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 369 of file reconstructor.h.


Constructor & Destructor Documentation

EMAN::FourierReconstructor::FourierReconstructor ( ) [inline]

Default constructor calls load_default_settings()

Definition at line 375 of file reconstructor.h.

References load_default_settings().

Referenced by NEW().

virtual EMAN::FourierReconstructor::~FourierReconstructor ( ) [inline, virtual]

Deconstructor calls free_memory()

Definition at line 380 of file reconstructor.h.

References free_memory().

{ free_memory(); }
EMAN::FourierReconstructor::FourierReconstructor ( const FourierReconstructor that) [private]

Disallow copy construction.


Member Function Documentation

void FourierReconstructor::clear ( ) [virtual]

clear the volume and tmp_data for use in Monte Carlo reconstructions

Reimplemented from EMAN::Reconstructor.

Definition at line 537 of file reconstructor.cpp.

References EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::ReconstructorVolumeData::image, EMAN::ReconstructorVolumeData::tmp_data, EMAN::EMData::to_zero(), and to_zero_cuda().

{
        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();
        
}
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_sliceThe EMData slice to be compared
eulerThe orientation of the slice as a Transform object
weightA weighting factor for this slice, generally the number of particles in a class-average. May be ignored by some reconstructors
subFlag 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:
NullPointerExceptionif the input EMData pointer is null
ImageFormatExceptionif the image is complex as opposed to real

Reimplemented from EMAN::Reconstructor.

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 709 of file reconstructor.cpp.

References EMAN::EMData::copy(), do_compare_slice_work(), do_insert_slice_work(), EMAN::EMData::get_attr(), EMAN::EMData::get_attr_default(), NullPointerException, preprocess_slice(), EMAN::EMData::set_attr(), EMAN::Transform::set_mirror(), EMAN::Transform::set_scale(), and EMAN::Transform::set_trans().

{
        // 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_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;

}
void FourierReconstructor::do_compare_slice_work ( EMData input_slice,
const Transform euler,
float  weight 
) [protected, virtual]

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

Parameters:
input_slicethe slice to insert into the 3D volume
eulera transform storing the slice euler angle

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 755 of file reconstructor.cpp.

References EMAN::Transform::copy_matrix_into_array(), determine_slice_agreement_cuda(), EMAN::dot(), dt, EMAN::EMData::get_data(), EMAN::Symmetry3D::get_symmetries(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::ReconstructorVolumeData::image, EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, pixel_at(), power(), EMAN::EMData::set_attr(), sqrt(), EMAN::ReconstructorVolumeData::subnx, EMAN::ReconstructorVolumeData::subny, EMAN::ReconstructorVolumeData::subnz, EMAN::ReconstructorVolumeData::tmp_data, weight, x, and y.

Referenced by determine_slice_agreement().

{

        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());

        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
        bool use_cpu = true;
        
#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]);
                                        dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
                                        vweight+=dt[2];
                                        power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
                                        power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
                                }
                        }
                        //cout << dot << " " << vweight << " " << power << " " << power2 << endl;
                }
        }
        
        dot/=sqrt(power*power2);                // normalize the dot product
//      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_absqual",(float)dot);
        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)(subnx*subny*subnz));
//      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)));
}
void FourierReconstructor::do_insert_slice_work ( const EMData *const  input_slice,
const Transform euler,
const float  weight 
) [protected, virtual]

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

Parameters:
input_slicethe slice to insert into the 3D volume
eulera transform storing the slice euler angle
weightweighting factor for this slice (usually number of particles in a class-average)

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 642 of file reconstructor.cpp.

References EMAN::Transform::copy_matrix_into_array(), EMAN::EMData::get_complex_at(), EMAN::Symmetry3D::get_symmetries(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), 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, weight, x, and y.

Referenced by determine_slice_agreement(), and insert_slice().

{
        // 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());

#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
        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;

                                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),weight);
                        }
                }
        }
}
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:
doiftA 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 924 of file reconstructor.cpp.

References EMAN::EMData::depad(), EMAN::EMData::do_ift_inplace(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Dict::has_key(), EMAN::ReconstructorVolumeData::image, EMAN::ReconstructorVolumeData::normalize_threed(), EMAN::FactoryBase::params, EMAN::EMData::process_inplace(), EMAN::Dict::set_default(), EMAN::ReconstructorVolumeData::tmp_data, EMAN::EMData::update(), and EMAN::EMData::write_image().

{
//      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);
        
//      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("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;
}
void FourierReconstructor::free_memory ( ) [protected, virtual]

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

Reimplemented from EMAN::ReconstructorVolumeData.

Definition at line 351 of file reconstructor.cpp.

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

Referenced by ~FourierReconstructor().

{
        if (image) { delete image; image=0; }
        if (tmp_data) { delete tmp_data; tmp_data=0; }
        if ( inserter != 0 )
        {
                delete inserter;
                inserter = 0;
        }
}
virtual string EMAN::FourierReconstructor::get_desc ( ) const [inline, virtual]

Get the one line description of the reconstructor.

Implements EMAN::FactoryBase.

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 458 of file reconstructor.h.

                {
                        return "Reconstruction via direct Fourier methods using one of a variety of different kernels, most of which are Gaussian based";
                }
virtual string EMAN::FourierReconstructor::get_name ( ) const [inline, virtual]

Get the unique name of the reconstructor.

Implements EMAN::FactoryBase.

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 451 of file reconstructor.h.

References NAME.

                {
                        return NAME;
                }
virtual TypeDict EMAN::FourierReconstructor::get_param_types ( ) const [inline, virtual]

Get the parameter types of this object.

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

Implements EMAN::FactoryBase.

Definition at line 474 of file reconstructor.h.

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

                {
                        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, "Optional. 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("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("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.");
                        return d;
                }
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:
slicethe image slice to be inserted into the 3D volume
eulerEuler angle of this image slice.
weightA 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:
NullPointerExceptionif the input EMData pointer is null
ImageFormatExceptionif the image is complex as opposed to real

Reimplemented from EMAN::Reconstructor.

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 602 of file reconstructor.cpp.

References EMAN::EMData::copy(), do_insert_slice_work(), EMAN::EMData::get_attr_default(), NullPointerException, preprocess_slice(), EMAN::Transform::set_mirror(), EMAN::Transform::set_scale(), and EMAN::Transform::set_trans().

{
        // 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;
/*      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();
        do_insert_slice_work(slice, *rotation, weight);
        
        delete rotation; rotation=0;
        delete slice;

//      image->update();
        return 0;
}
void FourierReconstructor::load_default_settings ( ) [protected, virtual]

Load default settings.

Definition at line 344 of file reconstructor.cpp.

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

Referenced by FourierReconstructor().

{
        inserter=0;
        image=0;
        tmp_data=0;
}
void FourierReconstructor::load_inserter ( ) [protected, virtual]

Load the pixel inserter based on the information in params.

Definition at line 364 of file reconstructor.cpp.

References EMAN::EMData::get_data(), EMAN::ReconstructorVolumeData::image, EMAN::FourierPixelInserter3D::init(), inserter, EMAN::FactoryBase::params, and EMAN::ReconstructorVolumeData::tmp_data.

Referenced by setup(), and setup_seed().

{
//      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();
}
static Reconstructor* EMAN::FourierReconstructor::NEW ( ) [inline, static]

Factory incorporation uses the pointer of this function.

Returns:
a Reconstructor pointer to a newly allocated FourierReconstructor

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 466 of file reconstructor.h.

References FourierReconstructor().

                {
                        return new FourierReconstructor();
                }
FourierReconstructor& EMAN::FourierReconstructor::operator= ( const FourierReconstructor ) [private]

Disallow assignment.

bool FourierReconstructor::pixel_at ( const float &  xx,
const float &  yy,
const float &  zz,
float *  dt 
) [protected, virtual]

This is a mode-2 pixel extractor.

Parameters:
xx,yy,zzvoxel coordinates (need not be integers)
dtfloat pointer with 3 floats allocated for returned complex value and weight sum

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 850 of file reconstructor.cpp.

References EMAN::Util::fast_exp(), EMAN::EMData::get_complex_index(), EMAN::EMData::get_complex_index_fast(), EMAN::EMData::get_data(), EMAN::Util::hypot3sq(), EMAN::EMConsts::I2G, EMAN::ReconstructorVolumeData::image, norm(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::nx2, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::ny2, EMAN::ReconstructorVolumeData::nz, EMAN::ReconstructorVolumeData::nz2, rdata, EMAN::ReconstructorVolumeData::subnx, EMAN::ReconstructorVolumeData::subx0, EMAN::ReconstructorVolumeData::suby0, EMAN::ReconstructorVolumeData::subz0, and EMAN::ReconstructorVolumeData::tmp_data.

Referenced by do_compare_slice_work().

{
        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;

                // no error checking on add_complex_fast, so we need to be careful here
                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 ++) {
                                        r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
                                        idx=image->get_complex_index_fast(i,j,k);
                                        gg = Util::fast_exp(-r / EMConsts::I2G);
                                        
                                        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 (normsum==0) return false;
                dt[0]/=normsum;
                dt[1]/=normsum;
                dt[2]/=normsum2;
//              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;
        }
}
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:
slicethe slice to be prepocessed
ttransform to be used for insertion
Exceptions:
InvalidValueExceptionwhen the specified padding value is less than the size of the images

Reimplemented from EMAN::Reconstructor.

Definition at line 560 of file reconstructor.cpp.

References EMAN::EMData::copy(), EMAN::EMData::do_fft_inplace(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::Transform::is_identity(), EMAN::EMData::mult(), EMAN::EMData::process(), EMAN::EMData::process_inplace(), EMAN::EMData::set_attr(), EMAN::Transform::set_rotation(), and sqrt().

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

{
#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;
}
void FourierReconstructor::setup ( ) [virtual]

Setup the Fourier reconstructor

Exceptions:
InvalidValueExceptionWhen one of the input parameters is invalid.

Implements EMAN::Reconstructor.

Definition at line 388 of file reconstructor.cpp.

References EMAN::Dict::has_key(), EMAN::ReconstructorVolumeData::image, ImageDimensionException, is_fftodd(), load_inserter(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::nx2, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::ny2, EMAN::ReconstructorVolumeData::nz, EMAN::ReconstructorVolumeData::nz2, EMAN::FactoryBase::params, EMAN::EMData::set_attr(), EMAN::EMData::set_complex(), EMAN::Dict::set_default(), EMAN::EMData::set_fftodd(), EMAN::EMData::set_ri(), EMAN::EMData::set_size(), sub(), EMAN::ReconstructorVolumeData::subnx, EMAN::ReconstructorVolumeData::subny, EMAN::ReconstructorVolumeData::subnz, EMAN::ReconstructorVolumeData::subx0, EMAN::ReconstructorVolumeData::suby0, EMAN::ReconstructorVolumeData::subz0, EMAN::ReconstructorVolumeData::tmp_data, EMAN::EMData::to_zero(), and EMAN::EMData::update().

{
        // 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;
        }


        // 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();

        if ( (bool) params["quiet"] == false )
        {
                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);
        }
}
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:
InvalidValueExceptionWhen one of the input parameters is invalid

Reimplemented from EMAN::Reconstructor.

Definition at line 466 of file reconstructor.cpp.

References EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Dict::has_key(), EMAN::ReconstructorVolumeData::image, ImageDimensionException, EMAN::EMData::is_complex(), is_fftodd(), load_inserter(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::nx2, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::ny2, EMAN::ReconstructorVolumeData::nz, EMAN::ReconstructorVolumeData::nz2, EMAN::FactoryBase::params, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::EMData::set_size(), sub(), EMAN::ReconstructorVolumeData::subnx, EMAN::ReconstructorVolumeData::subny, EMAN::ReconstructorVolumeData::subnz, EMAN::ReconstructorVolumeData::subx0, EMAN::ReconstructorVolumeData::suby0, EMAN::ReconstructorVolumeData::subz0, EMAN::ReconstructorVolumeData::tmp_data, and EMAN::EMData::to_value().

                                                                    {
        // 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;
        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);

        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

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

Definition at line 524 of file reconstructor.h.

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

const string FourierReconstructor::NAME = "fourier" [static]

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 488 of file reconstructor.h.

Referenced by get_name().


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