Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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:

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

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

Static Public Attributes

const string NAME = "fourier"

Protected Member Functions

void load_default_settings ()
 Load default settings.
void free_memory ()
 Frees the memory owned by this object (but not parent objects) Deletes the FourierPixelInserter3D pointer.
void load_inserter ()
 Load the pixel inserter based on the information in params.
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.
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.
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

Definition at line 363 of file reconstructor.h.


Constructor & Destructor Documentation

EMAN::FourierReconstructor::FourierReconstructor  )  [inline]
 

Default constructor calls load_default_settings().

Definition at line 369 of file reconstructor.h.

00369 { load_default_settings(); }

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

Deconstructor calls free_memory().

Definition at line 374 of file reconstructor.h.

References free_memory().

00374 { free_memory(); }

EMAN::FourierReconstructor::FourierReconstructor const FourierReconstructor that  )  [private]
 

Disallow copy construction.


Member Function Documentation

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.

Definition at line 606 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(), EMAN::Transform::set_trans(), and weight.

00607 {
00608         // Are these exceptions really necessary? (d.woolford)
00609         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00610 
00611         Transform * rotation;
00612         rotation = new Transform(arg); // assignment operator
00613 
00614         EMData *slice;
00615         if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00616         else slice = preprocess_slice( input_slice, *rotation);
00617 
00618 
00619         // We must use only the rotational component of the transform, scaling, translation and mirroring
00620         // are not implemented in Fourier space, but are in preprocess_slice
00621         rotation->set_scale(1.0);
00622         rotation->set_mirror(false);
00623         rotation->set_trans(0,0,0);
00624 
00625         // Remove the current slice first (not threadsafe, but otherwise performance would be awful)
00626         if (sub) do_insert_slice_work(slice, *rotation, -weight);
00627 
00628         // Compare
00629         do_compare_slice_work(slice, *rotation,weight);
00630 
00631         input_slice->set_attr("reconstruct_norm",slice->get_attr("reconstruct_norm"));
00632         input_slice->set_attr("reconstruct_absqual",slice->get_attr("reconstruct_absqual"));
00633 //      input_slice->set_attr("reconstruct_qual",slice->get_attr("reconstruct_qual"));
00634         input_slice->set_attr("reconstruct_weight",slice->get_attr("reconstruct_weight"));
00635 
00636         // Now put the slice back
00637         if (sub) do_insert_slice_work(slice, *rotation, weight);
00638 
00639 
00640         delete rotation;
00641         delete slice;
00642 
00643 //      image->update();
00644         return 0;
00645 
00646 }

void FourierReconstructor::do_compare_slice_work EMData input_slice,
const Transform euler,
float  weight
[protected]
 

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 transform3D storing the slice euler angle

Definition at line 648 of file reconstructor.cpp.

References dot(), dt, EMAN::EMData::get_data(), EMAN::Symmetry3D::get_symmetries(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), pixel_at(), power(), EMAN::EMData::set_attr(), sqrt(), EMAN::Vec3f, weight, x, and y.

Referenced by determine_slice_agreement().

00649 {
00650 
00651         float dt[3];    // This stores the complex and weight from the volume
00652         float dt2[2];   // This stores the local image complex
00653         float *dat = input_slice->get_data();
00654         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00655 
00656         float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
00657         float iny=(float)(input_slice->get_ysize());
00658 
00659         double dot=0;           // summed pixel*weight dot product
00660         double vweight=0;               // sum of weights
00661         double power=0;         // sum of inten*weight from volume
00662         double power2=0;                // sum of inten*weight from image
00663         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00664                 Transform t3d = arg*(*it);
00665                 for (int y = -iny/2; y < iny/2; y++) {
00666                         for (int x = 0; x <=  inx/2; x++) {
00667                                 if (x==0 && y==0) continue;             // We don't want to use the Fourier origin
00668 
00669                                 float rx = (float) x/(inx-2);   // coords relative to Nyquist=.5
00670                                 float ry = (float) y/iny;
00671 
00672 //                              if ((rx * rx + Util::square(ry - max_input_dim / 2)) > rl)
00673 //                                      continue;
00674 
00675                                 Vec3f coord(rx,ry,0);
00676                                 coord = coord*t3d; // transpose multiplication
00677                                 float xx = coord[0]; // transformed coordinates in terms of Nyquist
00678                                 float yy = coord[1];
00679                                 float zz = coord[2];
00680 
00681 
00682                                 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
00683 
00684                                 // Map back to actual pixel coordinates in output volume
00685                                 xx=xx*(nx-2);
00686                                 yy=yy*ny;
00687                                 zz=zz*nz;
00688 
00689 
00690                                 int idx = (int)(x * 2 + inx*(y<0?iny+y:y));
00691                                 dt2[0] = dat[idx];
00692                                 dt2[1] = dat[idx+1];
00693 
00694                                 // value returned indirectly in dt
00695                                 if (!pixel_at(xx,yy,zz,dt) || dt[2]==0) continue;
00696 
00697 //                              printf("%f\t%f\t%f\t%f\t%f\n",dt[0],dt[1],dt[2],dt2[0],dt2[1]);
00698                                 dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
00699                                 vweight+=dt[2];
00700                                 power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
00701                                 power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
00702                         }
00703                 }
00704         }
00705 
00706         dot/=sqrt(power*power2);                // normalize the dot product
00707         input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)));
00708         input_slice->set_attr("reconstruct_absqual",(float)dot);
00709         float rw=weight<=0?1.0f:1.0f/weight;
00710         input_slice->set_attr("reconstruct_qual",(float)(dot*rw/((rw-1.0)*dot+1.0)));   // here weight is a proxy for SNR
00711         input_slice->set_attr("reconstruct_weight",(float)vweight/(float)(subnx*subny*subnz));
00712 //      printf("** %g\t%g\t%g\t%g ##\n",dot,vweight,power,power2);
00713         //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)));
00714 }

void FourierReconstructor::do_insert_slice_work const EMData *const   input_slice,
const Transform euler,
const float  weight
[protected]
 

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 transform3D storing the slice euler angle
weight weighting factor for this slice (usually number of particles in a class-average)

Definition at line 556 of file reconstructor.cpp.

References EMAN::EMData::get_complex_at(), EMAN::Symmetry3D::get_symmetries(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::FourierPixelInserter3D::insert_pixel(), inserter, EMAN::Vec3f, x, and y.

Referenced by determine_slice_agreement(), and insert_slice().

00557 {
00558         // Reload the inserter if the mode has changed
00559 //      string mode = (string) params["mode"];
00560 //      if ( mode != inserter->get_name() )     load_inserter();
00561 
00562 //      int y_in = input_slice->get_ysize();
00563 //      int x_in = input_slice->get_xsize();
00564 //      // Adjust the dimensions to account for odd and even ffts
00565 //      if (input_slice->is_fftodd()) x_in -= 1;
00566 //      else x_in -= 2;
00567 
00568         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00569 
00570         float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
00571         float iny=(float)(input_slice->get_ysize());
00572 
00573         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00574                 Transform t3d = arg*(*it);
00575                 for (int y = -iny/2; y < iny/2; y++) {
00576                         for (int x = 0; x <=  inx/2; x++) {
00577 
00578                                 float rx = (float) x/(inx-2.0f);        // coords relative to Nyquist=.5
00579                                 float ry = (float) y/iny;
00580 
00581                                 Vec3f coord(rx,ry,0);
00582                                 coord = coord*t3d; // transpose multiplication
00583                                 float xx = coord[0]; // transformed coordinates in terms of Nyquist
00584                                 float yy = coord[1];
00585                                 float zz = coord[2];
00586 
00587                                 // Map back to real pixel coordinates in output volume
00588                                 xx=xx*(nx-2);
00589                                 yy=yy*ny;
00590                                 zz=zz*nz;
00591 
00592 //                              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",
00593 //                                      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);
00594 //                              if (x==0 && y==10) 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",
00595 //                                      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));
00596 
00597                                 //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);
00598 //                              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);
00599 //                              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);
00600                                 inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),weight);
00601                         }
00602                 }
00603         }
00604 }

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.

Definition at line 787 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::normalize_threed(), EMAN::EMData::process_inplace(), EMAN::Dict::set_default(), EMAN::EMData::update(), and EMAN::EMData::write_image().

00788 {
00789 //      float *norm = tmp_data->get_data();
00790 //      float *rdata = image->get_data();
00791 
00792         bool sqrtnorm=params.set_default("sqrtnorm",false);
00793         normalize_threed(sqrtnorm);
00794 
00795 //      tmp_data->write_image("density.mrc");
00796 
00797         // we may as well delete the tmp data now... it saves memory and the calling program might
00798         // need memory after it gets the return volume.
00799         // If this delete didn't happen now, it would happen when the deconstructor was called --david
00800         // no longer a good idea with the new iterative scheme -- steve
00801 //      if ( tmp_data != 0 )
00802 //      {
00803 //              delete tmp_data;
00804 //              tmp_data = 0;
00805 //      }
00806 
00807 /*      image->process_inplace("xform.fourierorigin.tocorner");*/
00808 
00809         if (doift) {
00810                 image->do_ift_inplace();
00811                 image->depad();
00812                 image->process_inplace("xform.phaseorigin.tocenter");
00813         }
00814         // If the image was padded it should be the original size, as the client would expect
00815         //  I blocked the rest, it is almost certainly incorrect  PAP 07/31/08
00816         // No, it's not incorrect. You are wrong. You have the meaning of nx mixed up. DSAW 09/23/cd
00817         // This should now be handled in the calling program --steve 11/03/09
00818 //      bool is_fftodd = (nx % 2 == 1);
00819 //      if ( (nx-2*(!is_fftodd)) != output_x || ny != output_y || nz != output_z )
00820 //      {
00821 //              FloatPoint origin( (nx-output_x)/2, (ny-output_y)/2, (nz-output_z)/2 );
00822 //              FloatSize region_size( output_x, output_y, output_z);
00823 //              Region clip_region( origin, region_size );
00824 //              image->clip_inplace( clip_region );
00825 //      }
00826 
00827         // Should be an "if (verbose)" here or something
00828         //print_stats(quality_scores);
00829 
00830         image->update();
00831         
00832         if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) {
00833                 if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
00834                 tmp_data->write_image((const char *)params["savenorm"]);
00835         }
00836 
00837         delete tmp_data;
00838         tmp_data=0;
00839         EMData *ret=image;
00840         image=0;
00841         
00842         return ret;
00843 }

void FourierReconstructor::free_memory  )  [protected]
 

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

Reimplemented from EMAN::ReconstructorVolumeData.

Definition at line 309 of file reconstructor.cpp.

References inserter.

00310 {
00311         if (image) { delete image; image=0; }
00312         if (tmp_data) { delete tmp_data; tmp_data=0; }
00313         if ( inserter != 0 )
00314         {
00315                 delete inserter;
00316                 inserter = 0;
00317         }
00318 }

virtual string EMAN::FourierReconstructor::get_desc  )  const [inline, virtual]
 

Get the one line description of the reconstructor.

Implements EMAN::FactoryBase.

Definition at line 448 of file reconstructor.h.

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

virtual string EMAN::FourierReconstructor::get_name  )  const [inline, virtual]
 

Get the unique name of the reconstructor.

Implements EMAN::FactoryBase.

Definition at line 441 of file reconstructor.h.

00442                 {
00443                         return NAME;
00444                 }

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

References EMAN::TypeDict::put().

00465                 {
00466                         TypeDict d;
00467                         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.");
00468                         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");
00469                         d.put("mode", EMObject::STRING, "Optional. Fourier pixel insertion mode name. gauss_2 is the default.");
00470                         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.");
00471                         d.put("verbose", EMObject::BOOL, "Optional. Toggles writing useful information to standard out. Default is false.");
00472                         d.put("subvolume",EMObject::INTARRAY, "Optional. (xorigin,yorigin,zorigin,xsize,ysize,zsize) all in Fourier pixels. Useful for parallelism.");
00473                         d.put("savenorm",EMObject::STRING, "Debug. Will cause the normalization volume to be written directly to the specified file when finish() is called.");
00474                         return d;
00475                 }

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.

Definition at line 523 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(), EMAN::Transform::set_trans(), and weight.

00524 {
00525         // Are these exceptions really necessary? (d.woolford)
00526         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00527 
00528         Transform * rotation;
00529 /*      if ( input_slice->has_attr("xform.projection") ) {
00530                 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
00531         } else {*/
00532         rotation = new Transform(arg); // assignment operator
00533 //      }
00534 
00535         EMData *slice;
00536         if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00537         else slice = preprocess_slice( input_slice, *rotation);
00538 
00539 
00540         // We must use only the rotational component of the transform, scaling, translation and mirroring
00541         // are not implemented in Fourier space, but are in preprocess_slice
00542         rotation->set_scale(1.0);
00543         rotation->set_mirror(false);
00544         rotation->set_trans(0,0,0);
00545 
00546         // Finally to the pixel wise slice insertion
00547         do_insert_slice_work(slice, *rotation, weight);
00548 
00549         delete rotation; rotation=0;
00550         delete slice;
00551 
00552 //      image->update();
00553         return 0;
00554 }

void EMAN::FourierReconstructor::load_default_settings  )  [inline, protected]
 

Load default settings.

Definition at line 483 of file reconstructor.h.

00484                 {
00485                         inserter=0;
00486                         image=0;
00487                         tmp_data=0;
00488                 }

void FourierReconstructor::load_inserter  )  [protected]
 

Load the pixel inserter based on the information in params.

Definition at line 322 of file reconstructor.cpp.

References EMAN::Factory< T >::get(), EMAN::EMData::get_data(), EMAN::FourierPixelInserter3D::init(), and inserter.

Referenced by setup(), and setup_seed().

00323 {
00324 //      ss
00325 //      string mode = (string)params["mode"];
00326         Dict parms;
00327         parms["data"] = image;
00328         parms["norm"] = tmp_data->get_data();
00329         // These aren't necessary because we deal with them before calling the inserter
00330 //      parms["subnx"] = nx;
00331 //      parms["subny"] = ny;
00332 //      parms["subnx"] = nz;
00333 //      parms["subx0"] = x0;
00334 //      parms["suby0"] = y0;
00335 //      parms["subz0"] = z0;
00336 
00337         if ( inserter != 0 )
00338         {
00339                 delete inserter;
00340         }
00341 
00342         inserter = Factory<FourierPixelInserter3D>::get((string)params["mode"], parms);
00343         inserter->init();
00344 }

Reconstructor* EMAN::FourierReconstructor::NEW  )  [inline, static]
 

Factory incorporation uses the pointer of this function.

Returns:
a Reconstructor pointer to a newly allocated FourierReconstructor

Definition at line 456 of file reconstructor.h.

00457                 {
00458                         return new FourierReconstructor();
00459                 }

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]
 

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

Definition at line 716 of file reconstructor.cpp.

References dt, EMAN::Util::fast_exp(), EMAN::EMData::get_complex_index(), EMAN::EMData::get_complex_index_fast(), EMAN::EMData::get_data(), EMAN::Util::hypot3sq(), norm(), and rdata.

Referenced by do_compare_slice_work().

00717 {
00718         int x0 = (int) floor(xx);
00719         int y0 = (int) floor(yy);
00720         int z0 = (int) floor(zz);
00721         
00722         float *rdata=image->get_data();
00723         float *norm=tmp_data->get_data();
00724         float normsum=0;
00725 
00726         dt[0]=dt[1]=dt[2]=0.0;
00727 
00728         if (nx==subnx) {                        // normal full reconstruction
00729                 if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
00730 
00731                 // no error checking on add_complex_fast, so we need to be careful here
00732                 int x1=x0+1;
00733                 int y1=y0+1;
00734                 int z1=z0+1;
00735                 if (x0<-nx2) x0=-nx2;
00736                 if (x1>nx2) x1=nx2;
00737                 if (y0<-ny2) y0=-ny2;
00738                 if (y1>ny2) y1=ny2;
00739                 if (z0<-nz2) z0=-nz2;
00740                 if (z1>nz2) z1=nz2;
00741                 
00742                 size_t idx;
00743                 float r, gg;
00744                 for (int k = z0 ; k <= z1; k++) {
00745                         for (int j = y0 ; j <= y1; j++) {
00746                                 for (int i = x0; i <= x1; i ++) {
00747                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00748                                         idx=image->get_complex_index_fast(i,j,k);
00749                                         gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
00750                                         
00751                                         dt[0]+=gg*rdata[idx];
00752                                         dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
00753                                         dt[2]+=norm[idx/2];
00754                                         normsum+=gg;                            
00755                                 }
00756                         }
00757                 }
00758                 dt[0]/=normsum;
00759                 dt[1]/=normsum;
00760 //              printf("%1.2f,%1.2f,%1.2f\t%1.3f\t%1.3f\t%1.3f\n",xx,yy,zz,dt[0],dt[1],dt[2]);
00761                 return true;
00762         } 
00763         else {                                  // for subvolumes, not optimized yet
00764                 size_t idx;
00765                 float r, gg;
00766                 for (int k = z0 ; k <= z0 + 1; k++) {
00767                         for (int j = y0 ; j <= y0 + 1; j++) {
00768                                 for (int i = x0; i <= x0 + 1; i ++) {
00769                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00770                                         idx=image->get_complex_index(i,j,k,subx0,suby0,subz0,nx,ny,nz);
00771                                         gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
00772                                         
00773                                         dt[0]+=gg*rdata[idx];
00774                                         dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
00775                                         dt[2]+=norm[idx/2];
00776                                         normsum+=gg;                            
00777                                 }
00778                         }
00779                 }
00780                 
00781                 if (normsum==0)  return false;
00782                 return true;
00783         }
00784 }

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 495 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 determine_slice_agreement(), and insert_slice().

00496 {
00497         // Shift the image pixels so the real space origin is now located at the phase origin (at the bottom left of the image)
00498         EMData* return_slice = 0;
00499         Transform tmp(t);
00500         tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly, this way we only do 2d translation,scaling and mirroring
00501 
00502         if (tmp.is_identity()) return_slice=slice->copy();
00503         else return_slice = slice->process("xform",Dict("transform",&tmp));
00504 
00505         return_slice->process_inplace("xform.phaseorigin.tocorner");
00506 
00507 //      printf("@@@ %d\n",(int)return_slice->get_attr("nx"));
00508         // Fourier transform the slice
00509         return_slice->do_fft_inplace();
00510 
00511 //      printf("%d\n",(int)return_slice->get_attr("nx"));
00512 
00513         return_slice->mult((float)sqrt(1.0f/(return_slice->get_ysize())*return_slice->get_xsize()));
00514 
00515         // Shift the Fourier transform so that it's origin is in the center (bottom) of the image.
00516 //      return_slice->process_inplace("xform.fourierorigin.tocenter");
00517 
00518         return_slice->set_attr("reconstruct_preproc",(int)1);
00519         return return_slice;
00520 }

void FourierReconstructor::setup  )  [virtual]
 

Setup the Fourier reconstructor.

Exceptions:
InvalidValueException When one of the input parameters is invalid

Implements EMAN::Reconstructor.

Definition at line 346 of file reconstructor.cpp.

References EMAN::Dict::has_key(), ImageDimensionException, is_fftodd(), load_inserter(), 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::EMData::to_zero(), and EMAN::EMData::update().

00347 {
00348         // default setting behavior - does not override if the parameter is already set
00349         params.set_default("mode","gauss_2");
00350 
00351         vector<int> size=params["size"];
00352 
00353         nx = size[0];
00354         ny = size[1];
00355         nz = size[2];
00356         nx2=nx/2-1;
00357         ny2=ny/2;
00358         nz2=nz/2;
00359 
00360 
00361         // Adjust nx if for Fourier transform even odd issues
00362         bool is_fftodd = (nx % 2 == 1);
00363         // The Fourier transform requires one extra pixel in the x direction,
00364         // which is two spaces in memory, one each for the complex and the
00365         // real components
00366         nx += 2-is_fftodd;
00367 
00368         if (params.has_key("subvolume")) {
00369                 vector<int> sub=params["subvolume"];
00370                 subx0=sub[0];
00371                 suby0=sub[1];
00372                 subz0=sub[2];
00373                 subnx=sub[3];
00374                 subny=sub[4];
00375                 subnz=sub[5];
00376 
00377                 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
00378                         throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
00379 
00380         }
00381         else {
00382                 subx0=suby0=subz0=0;
00383                 subnx=nx;
00384                 subny=ny;
00385                 subnz=nz;
00386         }
00387 
00388 
00389         // Odd dimension support is here atm, but not above.
00390         if (image) delete image;
00391         image = new EMData();
00392         image->set_size(subnx, subny, subnz);
00393         image->set_complex(true);
00394         image->set_fftodd(is_fftodd);
00395         image->set_ri(true);
00396         image->to_zero();
00397 
00398         if (params.has_key("subvolume")) {
00399                 image->set_attr("subvolume_x0",subx0);
00400                 image->set_attr("subvolume_y0",suby0);
00401                 image->set_attr("subvolume_z0",subz0);
00402                 image->set_attr("subvolume_full_nx",nx);
00403                 image->set_attr("subvolume_full_ny",ny);
00404                 image->set_attr("subvolume_full_nz",nz);
00405         }
00406         
00407         if (tmp_data) delete tmp_data;
00408         tmp_data = new EMData();
00409         tmp_data->set_size(subnx/2, subny, subnz);
00410         tmp_data->to_zero();
00411         tmp_data->update();
00412 
00413         load_inserter();
00414 
00415         if ( (bool) params["quiet"] == false )
00416         {
00417                 cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
00418                 cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
00419                 cout << "You will require approximately " << setprecision(3) << (subnx*subny*subnz*sizeof(float)*1.5)/1000000000.0 << "GB of memory to reconstruct this volume" << endl;
00420         }
00421 }

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

References EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Dict::has_key(), ImageDimensionException, EMAN::EMData::is_complex(), is_fftodd(), load_inserter(), EMAN::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::EMData::set_size(), sub(), and EMAN::EMData::to_value().

00423                                                                     {
00424         // default setting behavior - does not override if the parameter is already set
00425         params.set_default("mode","gauss_2");
00426 
00427         vector<int> size=params["size"];
00428 
00429         nx = size[0];
00430         ny = size[1];
00431         nz = size[2];
00432         nx2=nx/2-1;
00433         ny2=ny/2;
00434         nz2=nz/2;
00435 
00436 
00437         // Adjust nx if for Fourier transform even odd issues
00438         bool is_fftodd = (nx % 2 == 1);
00439         // The Fourier transform requires one extra pixel in the x direction,
00440         // which is two spaces in memory, one each for the complex and the
00441         // real components
00442         nx += 2-is_fftodd;
00443 
00444         if (params.has_key("subvolume")) {
00445                 vector<int> sub=params["subvolume"];
00446                 subx0=sub[0];
00447                 suby0=sub[1];
00448                 subz0=sub[2];
00449                 subnx=sub[3];
00450                 subny=sub[4];
00451                 subnz=sub[5];
00452 
00453                 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
00454                         throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
00455 
00456         }
00457         else {
00458                 subx0=suby0=subz0=0;
00459                 subnx=nx;
00460                 subny=ny;
00461                 subnz=nz;
00462         }
00463 
00464         if (seed->get_xsize()!=subnx || seed->get_ysize()!=subny || seed->get_zsize()!=subnz || !seed->is_complex())
00465                 throw ImageDimensionException("The dimensions of the seed volume do not match the reconstruction size");
00466 
00467         // Odd dimension support is here atm, but not above.
00468         image = seed;
00469         if (params.has_key("subvolume")) {
00470                 image->set_attr("subvolume_x0",subx0);
00471                 image->set_attr("subvolume_y0",suby0);
00472                 image->set_attr("subvolume_z0",subz0);
00473                 image->set_attr("subvolume_full_nx",nx);
00474                 image->set_attr("subvolume_full_ny",ny);
00475                 image->set_attr("subvolume_full_nz",nz);
00476         }
00477 
00478         if (tmp_data) delete tmp_data;
00479         tmp_data = new EMData();
00480         tmp_data->set_size(subnx/2, subny, subnz);
00481         tmp_data->to_value(seed_weight);
00482 
00483         load_inserter();
00484 
00485         if ( (bool) params["quiet"] == false )
00486         {
00487                 cout << "Seeded direct Fourier inversion";
00488                 cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
00489                 cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
00490                 cout << "You will require approximately " << setprecision(3) << (subnx*subny*subnz*sizeof(float)*1.5)/1000000000.0 << "GB of memory to reconstruct this volume" << endl;
00491         }
00492 }


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

Referenced by do_insert_slice_work(), free_memory(), and load_inserter().

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

Definition at line 73 of file reconstructor.cpp.


The documentation for this class was generated from the following files:
Generated on Fri Apr 30 15:39:31 2010 for EMAN2 by  doxygen 1.3.9.1