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

Fourier space 3D reconstruction This is a modified version of the normal FourierReconstructor which is aware of the SSNR information stored in individual class-average headers as "ctf_snr_total" and "ctf_wiener_filtered". More...

#include <reconstructor.h>

Inheritance diagram for EMAN::WienerFourierReconstructor:
Inheritance graph
[legend]
Collaboration diagram for EMAN::WienerFourierReconstructor:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 WienerFourierReconstructor ()
 Default constructor calls load_default_settings()
virtual ~WienerFourierReconstructor ()
 Deconstructor calls free_memory()
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 string get_name () const
 Get the unique name of the reconstructor.
virtual string get_desc () const
 Get the one line description of the reconstructor.

Static Public Member Functions

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

Static Public Attributes

static const string NAME = "wiener_fourier"

Protected Member Functions

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.

Private Member Functions

 WienerFourierReconstructor (const WienerFourierReconstructor &that)
 A pixel inserter pointer which inserts pixels into the 3D volume using one of a variety of insertion methods.
WienerFourierReconstructoroperator= (const WienerFourierReconstructor &)
 Disallow assignment.

Detailed Description

Fourier space 3D reconstruction This is a modified version of the normal FourierReconstructor which is aware of the SSNR information stored in individual class-average headers as "ctf_snr_total" and "ctf_wiener_filtered".

It will perform a reconstruction with a nonisotropic Wiener filter applied to the final reconstruction, and will 'undo' the Wiener filter on the individual class-averages if ctf_wiener_filtered is set. This represents something which was not possible to accomplish in EMAN1, and should produce superior results, with proper anisotropic filtering. Still, the filtration makes the assumption that the original SNR estimates were accurate, and that the data will average completely coherently, which is not truly the case. This may produce models which are somewhat underfiltered in the Wiener sense, but since B-factor corrections are not applied in the ctf.auto averager, this effect is likely already more than compensated for.

Definition at line 548 of file reconstructor.h.


Constructor & Destructor Documentation

EMAN::WienerFourierReconstructor::WienerFourierReconstructor ( ) [inline]

Default constructor calls load_default_settings()

Definition at line 554 of file reconstructor.h.

Referenced by NEW().

{};
virtual EMAN::WienerFourierReconstructor::~WienerFourierReconstructor ( ) [inline, virtual]

Deconstructor calls free_memory()

Definition at line 559 of file reconstructor.h.

{ }
EMAN::WienerFourierReconstructor::WienerFourierReconstructor ( const WienerFourierReconstructor that) [private]

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

Disallow copy construction


Member Function Documentation

int WienerFourierReconstructor::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
weightThis is ignored except for it's sign, since the SSNR from the particle header is used instead
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::FourierReconstructor.

Definition at line 1082 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, EMAN::FourierReconstructor::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");

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

//      tmp_data->write_image("dbug.hdf",0);
        
        // Remove the current slice first (not threadsafe, but otherwise performance would be awful)
        if (sub) do_insert_slice_work(slice, *rotation, -weight);

        // 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 WienerFourierReconstructor::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 from EMAN::FourierReconstructor.

Definition at line 1126 of file reconstructor.cpp.

References EMAN::dot(), dt, EMAN::EMData::get_data(), EMAN::Symmetry3D::get_symmetries(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), 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, 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
        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 / 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];
                        }
                }
        }

        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 WienerFourierReconstructor::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 from EMAN::FourierReconstructor.

Definition at line 1027 of file reconstructor.cpp.

References EMAN::EMData::get_attr(), EMAN::EMData::get_attr_default(), EMAN::EMData::get_complex_at(), EMAN::Symmetry3D::get_symmetries(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::FourierPixelInserter3D::insert_pixel(), EMAN::FourierReconstructor::inserter, EMAN::Util::linear_interpolate(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, sub(), weight, x, and y.

Referenced by determine_slice_agreement(), and insert_slice().

{

        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());
        
        int undo_wiener=(int)input_slice->get_attr_default("ctf_wiener_filtered",0);    // indicates whether we need to undo a wiener filter before insertion
//      if (undo_wiener) throw UnexpectedBehaviorException("wiener_fourier does not yet accept already Wiener filtered class-averages. Suggest using ctf.auto averager for now.");
        
        vector<float> snr=input_slice->get_attr("ctf_snr_total");
        float sub=1.0;
        if (inweight<0) sub=-1.0;
        float 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;

                                // This deals with the SNR weight
                                float rn = (float)hypot(rx,ry);
                                if (rn>=.5) continue;           // no SNR in the corners, and we're going to mask them later anyway
                                rn*=snr.size()*2.0f;
                                int rni=(int)floor(rn);
                                if ((unsigned int)rni>=snr.size()-1) weight=snr[snr.size()-1]*sub;
                                else {
                                        rn-=rni;
                                        weight=Util::linear_interpolate(snr[rni],snr[rni+1],rn);
                                }
//                              if (weight>500.0) printf("%f %d %d %f %f %d %f\n",weight,x,y,rx,ry,rni);
                                
                                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;

//                              printf("%f\n",weight);
                                if (undo_wiener) inserter->insert_pixel(xx,yy,zz,(input_slice->get_complex_at(x,y))*((weight+1.0f)/weight),weight*sub);
                                else inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),weight*sub);
                        }
                }
        }
}
EMData * WienerFourierReconstructor::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::FourierReconstructor.

Definition at line 1269 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().

{

        bool sqrtnorm=params.set_default("sqrtnorm",false);
        normalize_threed(sqrtnorm,true);                // true is the wiener filter

        if (doift) {
                image->do_ift_inplace();
                image->depad();
                image->process_inplace("xform.phaseorigin.tocenter");
        }

        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;
        EMData *ret=image;
        image=0;
        
        return ret;
}
virtual string EMAN::WienerFourierReconstructor::get_desc ( ) const [inline, virtual]

Get the one line description of the reconstructor.

Reimplemented from EMAN::FourierReconstructor.

Definition at line 609 of file reconstructor.h.

                {
                        return "Reconstruction via direct Fourier methods using one of a variety of different kernels, most of which are Gaussian based. This version also incorporates a nonisotropic Wiener filter based on SNR estimates stored in the class-average headers by the ctf.auto averager.";
                }
virtual string EMAN::WienerFourierReconstructor::get_name ( ) const [inline, virtual]

Get the unique name of the reconstructor.

Reimplemented from EMAN::FourierReconstructor.

Definition at line 602 of file reconstructor.h.

References NAME.

                {
                        return NAME;
                }
int WienerFourierReconstructor::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.
weightThis is ignored in this reconstructor, since the SSNR from the particle header is used instead
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::FourierReconstructor.

Definition at line 991 of file reconstructor.cpp.

References EMAN::EMData::copy(), do_insert_slice_work(), EMAN::EMData::get_attr_default(), EMAN::EMData::has_attr(), NotExistingObjectException, NullPointerException, EMAN::FourierReconstructor::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");

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

        if (!input_slice->has_attr("ctf_snr_total")) 
                throw NotExistingObjectException("ctf_snr_total","No SNR information present in class-average. Must use the ctf.auto or ctfw.auto averager.");

        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
        do_insert_slice_work(slice, *rotation, weight);

        delete rotation; rotation=0;
        delete slice;

//      image->update();
        return 0;
}
static Reconstructor* EMAN::WienerFourierReconstructor::NEW ( ) [inline, static]

Factory incorporation uses the pointer of this function.

Returns:
a Reconstructor pointer to a newly allocated WienerFourierReconstructor

Reimplemented from EMAN::FourierReconstructor.

Definition at line 617 of file reconstructor.h.

References WienerFourierReconstructor().

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

Disallow assignment.

bool WienerFourierReconstructor::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 from EMAN::FourierReconstructor.

Definition at line 1195 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;
        }
}

Member Data Documentation

const string WienerFourierReconstructor::NAME = "wiener_fourier" [static]

Reimplemented from EMAN::FourierReconstructor.

Definition at line 622 of file reconstructor.h.

Referenced by get_name().


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