EMAN::FourierReconstructor Class Reference
[a function or class that is CUDA enabled]

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:
InvalidValueException When 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=1.0)
 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 368 of file reconstructor.h.


Constructor & Destructor Documentation

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

Default constructor calls load_default_settings().

Definition at line 374 of file reconstructor.h.

References load_default_settings().

Referenced by NEW().

00374 { load_default_settings(); }

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

Deconstructor calls free_memory().

Definition at line 379 of file reconstructor.h.

References free_memory().

00379 { 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 535 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().

00536 {
00537         bool zeroimage = true;
00538         bool zerotmpimg = true;
00539         
00540 #ifdef EMAN2_USING_CUDA
00541         if(EMData::usecuda == 1) {
00542                 if(image->getcudarwdata()) {
00543                         to_zero_cuda(image->getcudarwdata(),image->get_xsize(),image->get_ysize(),image->get_zsize());
00544                         zeroimage = false;
00545                 }
00546                 if(tmp_data->getcudarwdata()) {
00547                         //to_zero_cuda(tmp_data->getcudarwdata(),image->get_xsize(),image->get_ysize(),image->get_zsize());
00548                         zerotmpimg = false;
00549                 }
00550         }
00551 #endif
00552 
00553         if(zeroimage) image->to_zero();
00554         if(zerotmpimg) tmp_data->to_zero();
00555         
00556 }

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 707 of file reconstructor.cpp.

References EMAN::EMData::copy(), do_compare_slice_work(), do_insert_slice_work(), 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().

00708 {
00709         // Are these exceptions really necessary? (d.woolford)
00710         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00711 
00712 #ifdef EMAN2_USING_CUDA
00713         if(EMData::usecuda == 1) {
00714                 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda(); //copy slice to cuda using the const version
00715         }
00716 #endif
00717 
00718         Transform * rotation;
00719         rotation = new Transform(arg); // assignment operator
00720 
00721         EMData *slice;
00722         if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00723         else slice = preprocess_slice( input_slice, *rotation);
00724 
00725 
00726         // We must use only the rotational component of the transform, scaling, translation and mirroring
00727         // are not implemented in Fourier space, but are in preprocess_slice
00728         rotation->set_scale(1.0);
00729         rotation->set_mirror(false);
00730         rotation->set_trans(0,0,0);
00731         if (sub) do_insert_slice_work(slice, *rotation, -weight);
00732         // Remove the current slice first (not threadsafe, but otherwise performance would be awful)
00733         
00734         // Compare
00735         do_compare_slice_work(slice, *rotation,weight);
00736 
00737         input_slice->set_attr("reconstruct_norm",slice->get_attr("reconstruct_norm"));
00738         input_slice->set_attr("reconstruct_absqual",slice->get_attr("reconstruct_absqual"));
00739 //      input_slice->set_attr("reconstruct_qual",slice->get_attr("reconstruct_qual"));
00740         input_slice->set_attr("reconstruct_weight",slice->get_attr("reconstruct_weight"));
00741 
00742         // Now put the slice back
00743         if (sub) do_insert_slice_work(slice, *rotation, weight);
00744 
00745         delete rotation;
00746         delete slice;
00747 
00748 //      image->update();
00749         return 0;
00750 
00751 }

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_slice the slice to insert into the 3D volume
euler a transform storing the slice euler angle

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 753 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, x, and y.

Referenced by determine_slice_agreement().

00754 {
00755 
00756         float dt[3];    // This stores the complex and weight from the volume
00757         float dt2[2];   // This stores the local image complex
00758         float *dat = input_slice->get_data();
00759         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00760 
00761         float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
00762         float iny=(float)(input_slice->get_ysize());
00763 
00764         double dot=0;           // summed pixel*weight dot product
00765         double vweight=0;               // sum of weights
00766         double power=0;         // sum of inten*weight from volume
00767         double power2=0;                // sum of inten*weight from image
00768         bool use_cpu = true;
00769         
00770 #ifdef EMAN2_USING_CUDA
00771         if(EMData::usecuda == 1) {
00772                 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda();
00773                 if(!image->getcudarwdata()){
00774                         image->copy_to_cuda();
00775                         tmp_data->copy_to_cuda();
00776                 }
00777                 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00778                         Transform t3d = arg*(*it);
00779                         float * m = new float[12];
00780                         t3d.copy_matrix_into_array(m);
00781                         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);
00782                         dot = stats.x;
00783                         vweight = stats.y;
00784                         power = stats.z;
00785                         power2 = stats.w;
00786                         //cout << "CUDA stats " << stats.x << " " << stats.y << " " << stats.z << " " << stats.w << endl;
00787                         use_cpu = false;
00788                 }
00789         }
00790 #endif
00791         if(use_cpu) {
00792                 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00793                         Transform t3d = arg*(*it);
00794                         for (int y = -iny/2; y < iny/2; y++) {
00795                                 for (int x = 0; x <=  inx/2; x++) {
00796                                         if (x==0 && y==0) continue;             // We don't want to use the Fourier origin
00797 
00798                                         float rx = (float) x/(inx-2);   // coords relative to Nyquist=.5
00799                                         float ry = (float) y/iny;
00800 
00801 //                                      if ((rx * rx + Util::square(ry - max_input_dim "xform.projection"/ 2)) > rl)
00802 //                                      continue;
00803 
00804                                         Vec3f coord(rx,ry,0);
00805                                         coord = coord*t3d; // transpose multiplication
00806                                         float xx = coord[0]; // transformed coordinates in terms of Nyquist
00807                                         float yy = coord[1];
00808                                         float zz = coord[2];
00809 
00810 
00811                                         if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
00812 
00813                                         // Map back to actual pixel coordinates in output volume
00814                                         xx=xx*(nx-2);
00815                                         yy=yy*ny;
00816                                         zz=zz*nz;
00817 
00818 
00819                                         int idx = (int)(x * 2 + inx*(y<0?iny+y:y));
00820                                         dt2[0] = dat[idx];
00821                                         dt2[1] = dat[idx+1];
00822 
00823                                         // value returned indirectly in dt
00824                                         if (!pixel_at(xx,yy,zz,dt) || dt[2]==0) continue;
00825 
00826 //                                      printf("%f\t%f\t%f\t%f\t%f\n",dt[0],dt[1],dt[2],dt2[0],dt2[1]);
00827                                         dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
00828                                         vweight+=dt[2];
00829                                         power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
00830                                         power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
00831                                 }
00832                         }
00833                         //cout << dot << " " << vweight << " " << power << " " << power2 << endl;
00834                 }
00835         }
00836         
00837         dot/=sqrt(power*power2);                // normalize the dot product
00838 //      input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)));
00839         input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)));
00840         input_slice->set_attr("reconstruct_absqual",(float)dot);
00841         float rw=weight<=0?1.0f:1.0f/weight;
00842         input_slice->set_attr("reconstruct_qual",(float)(dot*rw/((rw-1.0)*dot+1.0)));   // here weight is a proxy for SNR
00843         input_slice->set_attr("reconstruct_weight",(float)vweight/(float)(subnx*subny*subnz));
00844 //      printf("** %g\t%g\t%g\t%g ##\n",dot,vweight,power,power2);
00845         //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)));
00846 }

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

Reimplemented in EMAN::WienerFourierReconstructor.

Definition at line 640 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, x, and y.

Referenced by determine_slice_agreement(), and insert_slice().

00641 {
00642         // Reload the inserter if the mode has changed
00643 //      string mode = (string) params["mode"];
00644 //      if ( mode != inserter->get_name() )     load_inserter();
00645 
00646 //      int y_in = input_slice->get_ysize();
00647 //      int x_in = input_slice->get_xsize();
00648 //      // Adjust the dimensions to account for odd and even ffts
00649 //      if (input_slice->is_fftodd()) x_in -= 1;
00650 //      else x_in -= 2;
00651 
00652         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00653 
00654         float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
00655         float iny=(float)(input_slice->get_ysize());
00656 
00657 #ifdef EMAN2_USING_CUDA
00658         if(EMData::usecuda == 1) {
00659                 if(!image->getcudarwdata()){
00660                         image->copy_to_cuda();
00661                         tmp_data->copy_to_cuda();
00662                 }
00663                 float * m = new float[12];
00664                 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00665                         Transform t3d = arg*(*it);
00666                         t3d.copy_matrix_into_array(m);
00667                         //cout << "using CUDA " << image->getcudarwdata() << endl;
00668                         insert_slice_cuda(m,input_slice->getcudarwdata(),image->getcudarwdata(),tmp_data->getcudarwdata(),inx,iny,image->get_xsize(),image->get_ysize(),image->get_zsize(), weight);
00669                 }
00670                 delete m;
00671                 return;
00672         }
00673 #endif
00674         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00675                 Transform t3d = arg*(*it);
00676                 for (int y = -iny/2; y < iny/2; y++) {
00677                         for (int x = 0; x <=  inx/2; x++) {
00678 
00679                                 float rx = (float) x/(inx-2.0f);        // coords relative to Nyquist=.5
00680                                 float ry = (float) y/iny;
00681 
00682                                 Vec3f coord(rx,ry,0);
00683                                 coord = coord*t3d; // transpose multiplication
00684                                 float xx = coord[0]; // transformed coordinates in terms of Nyquist
00685                                 float yy = coord[1];
00686                                 float zz = coord[2];
00687 
00688                                 // Map back to real pixel coordinates in output volume
00689                                 xx=xx*(nx-2);
00690                                 yy=yy*ny;
00691                                 zz=zz*nz;
00692 
00693 //                              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",
00694 //                                      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);
00695 //                              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",
00696 //                                      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));
00697 
00698                                 //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);
00699 //                              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);
00700 //                              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);
00701                                 inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),weight);
00702                         }
00703                 }
00704         }
00705 }

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 922 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().

00923 {
00924 //      float *norm = tmp_data->get_data();
00925 //      float *rdata = image->get_data();
00926 #ifdef EMAN2_USING_CUDA
00927         if(EMData::usecuda == 1 && image->getcudarwdata()){
00928                 cout << "copy back from CUDA" << endl;
00929                 image->copy_from_device();
00930                 tmp_data->copy_from_device();
00931         }
00932 #endif
00933         
00934         bool sqrtnorm=params.set_default("sqrtnorm",false);
00935         normalize_threed(sqrtnorm);
00936         
00937 //      tmp_data->write_image("density.mrc");
00938 
00939         // we may as well delete the tmp data now... it saves memory and the calling program might
00940         // need memory after it gets the return volume.
00941         // If this delete didn't happen now, it would happen when the deconstructor was called --david
00942         // no longer a good idea with the new iterative scheme -- steve
00943 //      if ( tmp_data != 0 )
00944 //      {
00945 //              delete tmp_data;
00946 //              tmp_data = 0;
00947 //      }
00948 
00949 /*      image->process_inplace("xform.fourierorigin.tocorner");*/
00950 
00951         if (doift) {
00952                 image->do_ift_inplace();
00953                 image->depad();
00954                 image->process_inplace("xform.phaseorigin.tocenter");
00955         }
00956         // If the image was padded it should be the original size, as the client would expect
00957         //  I blocked the rest, it is almost certainly incorrect  PAP 07/31/08
00958         // No, it's not incorrect. You are wrong. You have the meaning of nx mixed up. DSAW 09/23/cd
00959         // This should now be handled in the calling program --steve 11/03/09
00960 //      bool is_fftodd = (nx % 2 == 1);
00961 //      if ( (nx-2*(!is_fftodd)) != output_x || ny != output_y || nz != output_z )
00962 //      {
00963 //              FloatPoint origin( (nx-output_x)/2, (ny-output_y)/2, (nz-output_z)/2 );
00964 //              FloatSize region_size( output_x, output_y, output_z);
00965 //              Region clip_region( origin, region_size );
00966 //              image->clip_inplace( clip_region );
00967 //      }
00968 
00969         // Should be an "if (verbose)" here or something
00970         //print_stats(quality_scores);
00971 
00972         image->update();
00973         
00974         if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) {
00975                 if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
00976                 tmp_data->write_image((const char *)params["savenorm"]);
00977         }
00978 
00979         delete tmp_data;
00980         tmp_data=0;
00981         //Since we give up the ownership of the pointer to-be-returned, it's caller's responsibility to delete the returned image.
00982         //So we wrap this function with return_value_policy< manage_new_object >() in libpyReconstructor2.cpp to hand over ownership to Python.
00983         EMData *ret=image;
00984         image=0;
00985         
00986         return ret;
00987 }

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 349 of file reconstructor.cpp.

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

Referenced by ~FourierReconstructor().

00350 {
00351         if (image) { delete image; image=0; }
00352         if (tmp_data) { delete tmp_data; tmp_data=0; }
00353         if ( inserter != 0 )
00354         {
00355                 delete inserter;
00356                 inserter = 0;
00357         }
00358 }

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

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

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

References NAME.

00451                 {
00452                         return NAME;
00453                 }

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

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

00474                 {
00475                         TypeDict d;
00476                         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.");
00477                         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");
00478                         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.");
00479                         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.");
00480                         d.put("verbose", EMObject::BOOL, "Optional. Toggles writing useful information to standard out. Default is false.");
00481                         d.put("quiet", EMObject::BOOL, "Optional. If false, print verbose information.");
00482                         d.put("subvolume",EMObject::INTARRAY, "Optional. (xorigin,yorigin,zorigin,xsize,ysize,zsize) all in Fourier pixels. Useful for parallelism.");
00483                         d.put("savenorm",EMObject::STRING, "Debug. Will cause the normalization volume to be written directly to the specified file when finish() is called.");
00484                         return d;
00485                 }

int FourierReconstructor::insert_slice ( const EMData *const   slice,
const Transform euler,
const float  weight = 1.0 
) [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 600 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().

00601 {
00602         // Are these exceptions really necessary? (d.woolford)
00603         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00604 
00605 #ifdef EMAN2_USING_CUDA
00606         if(EMData::usecuda == 1) {
00607                 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda(); //copy slice to cuda using the const version
00608         }
00609 #endif
00610 
00611         Transform * rotation;
00612 /*      if ( input_slice->has_attr("xform.projection") ) {
00613                 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
00614         } else {*/
00615         rotation = new Transform(arg); // assignment operator
00616 //      }
00617 
00618         EMData *slice;
00619         if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00620         else slice = preprocess_slice( input_slice, *rotation);
00621 
00622 
00623         // We must use only the rotational component of the transform, scaling, translation and mirroring
00624         // are not implemented in Fourier space, but are in preprocess_slice
00625         rotation->set_scale(1.0);
00626         rotation->set_mirror(false);
00627         rotation->set_trans(0,0,0);
00628 
00629         // Finally to the pixel wise slice insertion
00630         //slice->copy_to_cuda();
00631         do_insert_slice_work(slice, *rotation, weight);
00632         
00633         delete rotation; rotation=0;
00634         delete slice;
00635 
00636 //      image->update();
00637         return 0;
00638 }

void FourierReconstructor::load_default_settings (  )  [protected, virtual]

Load default settings.

Definition at line 342 of file reconstructor.cpp.

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

Referenced by FourierReconstructor().

00343 {
00344         inserter=0;
00345         image=0;
00346         tmp_data=0;
00347 }

void FourierReconstructor::load_inserter (  )  [protected, virtual]

Load the pixel inserter based on the information in params.

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

00363 {
00364 //      ss
00365 //      string mode = (string)params["mode"];
00366         Dict parms;
00367         parms["data"] = image;
00368         parms["norm"] = tmp_data->get_data();
00369         // These aren't necessary because we deal with them before calling the inserter
00370 //      parms["subnx"] = nx;
00371 //      parms["subny"] = ny;
00372 //      parms["subnx"] = nz;
00373 //      parms["subx0"] = x0;
00374 //      parms["suby0"] = y0;
00375 //      parms["subz0"] = z0;
00376 
00377         if ( inserter != 0 )
00378         {
00379                 delete inserter;
00380         }
00381 
00382         inserter = Factory<FourierPixelInserter3D>::get((string)params["mode"], parms);
00383         inserter->init();
00384 }

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

References FourierReconstructor().

00466                 {
00467                         return new FourierReconstructor();
00468                 }

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,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 848 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().

00849 {
00850         int x0 = (int) floor(xx);
00851         int y0 = (int) floor(yy);
00852         int z0 = (int) floor(zz);
00853         
00854         float *rdata=image->get_data();
00855         float *norm=tmp_data->get_data();
00856         float normsum=0,normsum2=0;
00857 
00858         dt[0]=dt[1]=dt[2]=0.0;
00859 
00860         if (nx==subnx) {                        // normal full reconstruction
00861                 if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
00862 
00863                 // no error checking on add_complex_fast, so we need to be careful here
00864                 int x1=x0+1;
00865                 int y1=y0+1;
00866                 int z1=z0+1;
00867                 if (x0<-nx2) x0=-nx2;
00868                 if (x1>nx2) x1=nx2;
00869                 if (y0<-ny2) y0=-ny2;
00870                 if (y1>ny2) y1=ny2;
00871                 if (z0<-nz2) z0=-nz2;
00872                 if (z1>nz2) z1=nz2;
00873                 
00874                 size_t idx=0;
00875                 float r, gg;
00876                 for (int k = z0 ; k <= z1; k++) {
00877                         for (int j = y0 ; j <= y1; j++) {
00878                                 for (int i = x0; i <= x1; i ++) {
00879                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00880                                         idx=image->get_complex_index_fast(i,j,k);
00881                                         gg = Util::fast_exp(-r / EMConsts::I2G);
00882                                         
00883                                         dt[0]+=gg*rdata[idx];
00884                                         dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
00885                                         dt[2]+=norm[idx/2]*gg;
00886                                         normsum2+=gg;
00887                                         normsum+=gg*norm[idx/2];                                
00888                                 }
00889                         }
00890                 }
00891                 if (normsum==0) return false;
00892                 dt[0]/=normsum;
00893                 dt[1]/=normsum;
00894                 dt[2]/=normsum2;
00895 //              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]);
00896                 return true;
00897         } 
00898         else {                                  // for subvolumes, not optimized yet
00899                 size_t idx;
00900                 float r, gg;
00901                 for (int k = z0 ; k <= z0 + 1; k++) {
00902                         for (int j = y0 ; j <= y0 + 1; j++) {
00903                                 for (int i = x0; i <= x0 + 1; i ++) {
00904                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00905                                         idx=image->get_complex_index(i,j,k,subx0,suby0,subz0,nx,ny,nz);
00906                                         gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
00907                                         
00908                                         dt[0]+=gg*rdata[idx];
00909                                         dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
00910                                         dt[2]+=norm[idx/2];
00911                                         normsum+=gg;                            
00912                                 }
00913                         }
00914                 }
00915                 
00916                 if (normsum==0)  return false;
00917                 return true;
00918         }
00919 }

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 558 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(), sqrt(), and t.

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

00559 {
00560 #ifdef EMAN2_USING_CUDA
00561         if(EMData::usecuda == 1) {
00562                 if(!slice->getcudarwdata()) slice->copy_to_cuda(); //copy slice to cuda using the const version
00563         }
00564 #endif
00565         // Shift the image pixels so the real space origin is now located at the phase origin (at the bottom left of the image)
00566         EMData* return_slice = 0;
00567         Transform tmp(t);
00568         tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly, this way we only do 2d translation,scaling and mirroring
00569 
00570         if (tmp.is_identity()) return_slice=slice->copy();
00571         else return_slice = slice->process("xform",Dict("transform",&tmp));
00572 
00573         return_slice->process_inplace("xform.phaseorigin.tocorner");
00574 
00575 //      printf("@@@ %d\n",(int)return_slice->get_attr("nx"));
00576         // Fourier transform the slice
00577         
00578 #ifdef EMAN2_USING_CUDA
00579         if(EMData::usecuda == 1 && return_slice->getcudarwdata()) {
00580                 return_slice->do_fft_inplace_cuda(); //a CUDA FFT inplace is quite slow as there is a lot of mem copying.
00581         }else{
00582                 return_slice->do_fft_inplace();
00583         }
00584 #else
00585         return_slice->do_fft_inplace();
00586 #endif
00587 
00588 //      printf("%d\n",(int)return_slice->get_attr("nx"));
00589 
00590         return_slice->mult((float)sqrt(1.0f/(return_slice->get_ysize())*return_slice->get_xsize()));
00591 
00592         // Shift the Fourier transform so that it's origin is in the center (bottom) of the image.
00593 //      return_slice->process_inplace("xform.fourierorigin.tocenter");
00594 
00595         return_slice->set_attr("reconstruct_preproc",(int)1);
00596         return return_slice;
00597 }

void FourierReconstructor::setup (  )  [virtual]

Setup the Fourier reconstructor

Exceptions:
InvalidValueException When one of the input parameters is invalid.

Implements EMAN::Reconstructor.

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

00387 {
00388         // default setting behavior - does not override if the parameter is already set
00389         params.set_default("mode","gauss_2");
00390 
00391         vector<int> size=params["size"];
00392 
00393         nx = size[0];
00394         ny = size[1];
00395         nz = size[2];
00396         nx2=nx/2-1;
00397         ny2=ny/2;
00398         nz2=nz/2;
00399 
00400 
00401         // Adjust nx if for Fourier transform even odd issues
00402         bool is_fftodd = (nx % 2 == 1);
00403         // The Fourier transform requires one extra pixel in the x direction,
00404         // which is two spaces in memory, one each for the complex and the
00405         // real components
00406         nx += 2-is_fftodd;
00407 
00408         if (params.has_key("subvolume")) {
00409                 vector<int> sub=params["subvolume"];
00410                 subx0=sub[0];
00411                 suby0=sub[1];
00412                 subz0=sub[2];
00413                 subnx=sub[3];
00414                 subny=sub[4];
00415                 subnz=sub[5];
00416 
00417                 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
00418                         throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
00419 
00420         }
00421         else {
00422                 subx0=suby0=subz0=0;
00423                 subnx=nx;
00424                 subny=ny;
00425                 subnz=nz;
00426         }
00427 
00428 
00429         // Odd dimension support is here atm, but not above.
00430         if (image) delete image;
00431         image = new EMData();
00432         image->set_size(subnx, subny, subnz);
00433         image->set_complex(true);
00434         image->set_fftodd(is_fftodd);
00435         image->set_ri(true);
00436 //      printf("%d %d %d\n\n",subnx,subny,subnz);
00437         image->to_zero();
00438 
00439         if (params.has_key("subvolume")) {
00440                 image->set_attr("subvolume_x0",subx0);
00441                 image->set_attr("subvolume_y0",suby0);
00442                 image->set_attr("subvolume_z0",subz0);
00443                 image->set_attr("subvolume_full_nx",nx);
00444                 image->set_attr("subvolume_full_ny",ny);
00445                 image->set_attr("subvolume_full_nz",nz);
00446         }
00447         
00448         if (tmp_data) delete tmp_data;
00449         tmp_data = new EMData();
00450         tmp_data->set_size(subnx/2, subny, subnz);
00451         tmp_data->to_zero();
00452         tmp_data->update();
00453 
00454         load_inserter();
00455 
00456         if ( (bool) params["quiet"] == false )
00457         {
00458                 cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
00459                 cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
00460                 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);
00461         }
00462 }

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 464 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().

00464                                                                     {
00465         // default setting behavior - does not override if the parameter is already set
00466         params.set_default("mode","gauss_2");
00467 
00468         vector<int> size=params["size"];
00469 
00470         nx = size[0];
00471         ny = size[1];
00472         nz = size[2];
00473         nx2=nx/2-1;
00474         ny2=ny/2;
00475         nz2=nz/2;
00476 
00477 
00478         // Adjust nx if for Fourier transform even odd issues
00479         bool is_fftodd = (nx % 2 == 1);
00480         // The Fourier transform requires one extra pixel in the x direction,
00481         // which is two spaces in memory, one each for the complex and the
00482         // real components
00483         nx += 2-is_fftodd;
00484 
00485         if (params.has_key("subvolume")) {
00486                 vector<int> sub=params["subvolume"];
00487                 subx0=sub[0];
00488                 suby0=sub[1];
00489                 subz0=sub[2];
00490                 subnx=sub[3];
00491                 subny=sub[4];
00492                 subnz=sub[5];
00493 
00494                 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
00495                         throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
00496 
00497         }
00498         else {
00499                 subx0=suby0=subz0=0;
00500                 subnx=nx;
00501                 subny=ny;
00502                 subnz=nz;
00503         }
00504 
00505         if (seed->get_xsize()!=subnx || seed->get_ysize()!=subny || seed->get_zsize()!=subnz || !seed->is_complex())
00506                 throw ImageDimensionException("The dimensions of the seed volume do not match the reconstruction size");
00507 
00508         // Odd dimension support is here atm, but not above.
00509         image = seed;
00510         if (params.has_key("subvolume")) {
00511                 image->set_attr("subvolume_x0",subx0);
00512                 image->set_attr("subvolume_y0",suby0);
00513                 image->set_attr("subvolume_z0",subz0);
00514                 image->set_attr("subvolume_full_nx",nx);
00515                 image->set_attr("subvolume_full_ny",ny);
00516                 image->set_attr("subvolume_full_nz",nz);
00517         }
00518 
00519         if (tmp_data) delete tmp_data;
00520         tmp_data = new EMData();
00521         tmp_data->set_size(subnx/2, subny, subnz);
00522         tmp_data->to_value(seed_weight);
00523 
00524         load_inserter();
00525 
00526         if ( (bool) params["quiet"] == false )
00527         {
00528                 cout << "Seeded direct Fourier inversion";
00529                 cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
00530                 cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
00531                 cout << "You will require approximately " << setprecision(3) << (subnx*subny*subnz*sizeof(float)*1.5)/1000000000.0 << "GB of memory to reconstruct this volume" << endl;
00532         }
00533 }


Member Data Documentation

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

Referenced by get_name().


The documentation for this class was generated from the following files:
Generated on Tue Jun 11 12:44:58 2013 for EMAN2 by  doxygen 1.4.7