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

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

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

Static Public Attributes

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

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.

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

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

00520 {
00521         bool zeroimage = true;
00522         bool zerotmpimg = true;
00523         
00524 #ifdef EMAN2_USING_CUDA
00525         if(EMData::usecuda == 1) {
00526                 if(image->getcudarwdata()) {
00527                         to_zero_cuda(image->getcudarwdata(),image->get_xsize(),image->get_ysize(),image->get_zsize());
00528                         zeroimage = false;
00529                 }
00530                 if(tmp_data->getcudarwdata()) {
00531                         //to_zero_cuda(tmp_data->getcudarwdata(),image->get_xsize(),image->get_ysize(),image->get_zsize());
00532                         zerotmpimg = false;
00533                 }
00534         }
00535 #endif
00536 
00537         if(zeroimage) image->to_zero();
00538         if(zerotmpimg) tmp_data->to_zero();
00539         
00540 }

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

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

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

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

Referenced by determine_slice_agreement().

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

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 624 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::FourierPixelInserter3D::insert_pixel(), insert_slice_cuda(), inserter, EMAN::Vec3f, x, and y.

Referenced by determine_slice_agreement(), and insert_slice().

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

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

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

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

References inserter.

00334 {
00335         if (image) { delete image; image=0; }
00336         if (tmp_data) { delete tmp_data; tmp_data=0; }
00337         if ( inserter != 0 )
00338         {
00339                 delete inserter;
00340                 inserter = 0;
00341         }
00342 }

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.

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::TypeDict::put().

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

00585 {
00586         // Are these exceptions really necessary? (d.woolford)
00587         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00588 
00589 #ifdef EMAN2_USING_CUDA
00590         if(EMData::usecuda == 1) {
00591                 if(!input_slice->getcudarwdata()) input_slice->copy_to_cuda_keepcpu(); //copy slice to cuda using the const version
00592         }
00593 #endif
00594 
00595         Transform * rotation;
00596 /*      if ( input_slice->has_attr("xform.projection") ) {
00597                 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
00598         } else {*/
00599         rotation = new Transform(arg); // assignment operator
00600 //      }
00601 
00602         EMData *slice;
00603         if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00604         else slice = preprocess_slice( input_slice, *rotation);
00605 
00606 
00607         // We must use only the rotational component of the transform, scaling, translation and mirroring
00608         // are not implemented in Fourier space, but are in preprocess_slice
00609         rotation->set_scale(1.0);
00610         rotation->set_mirror(false);
00611         rotation->set_trans(0,0,0);
00612 
00613         // Finally to the pixel wise slice insertion
00614         //slice->copy_to_cuda();
00615         do_insert_slice_work(slice, *rotation, weight);
00616         
00617         delete rotation; rotation=0;
00618         delete slice;
00619 
00620 //      image->update();
00621         return 0;
00622 }

void FourierReconstructor::load_default_settings  )  [protected, virtual]
 

Load default settings.

Definition at line 326 of file reconstructor.cpp.

References inserter.

00327 {
00328         inserter=0;
00329         image=0;
00330         tmp_data=0;
00331 }

void FourierReconstructor::load_inserter  )  [protected, virtual]
 

Load the pixel inserter based on the information in params.

Definition at line 346 of file reconstructor.cpp.

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

Referenced by setup(), and setup_seed().

00347 {
00348 //      ss
00349 //      string mode = (string)params["mode"];
00350         Dict parms;
00351         parms["data"] = image;
00352         parms["norm"] = tmp_data->get_data();
00353         // These aren't necessary because we deal with them before calling the inserter
00354 //      parms["subnx"] = nx;
00355 //      parms["subny"] = ny;
00356 //      parms["subnx"] = nz;
00357 //      parms["subx0"] = x0;
00358 //      parms["suby0"] = y0;
00359 //      parms["subz0"] = z0;
00360 
00361         if ( inserter != 0 )
00362         {
00363                 delete inserter;
00364         }
00365 
00366         inserter = Factory<FourierPixelInserter3D>::get((string)params["mode"], parms);
00367         inserter->init();
00368 }

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.

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

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

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

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

void FourierReconstructor::setup  )  [virtual]
 

Setup the Fourier reconstructor.

Exceptions:
InvalidValueException When one of the input parameters is invalid

Implements EMAN::Reconstructor.

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

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

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

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


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


The documentation for this class was generated from the following files:
Generated on Tue Jul 12 13:52:55 2011 for EMAN2 by  doxygen 1.3.9.1