#include <reconstructor.h>
Inheritance diagram for EMAN::FourierReconstructor:
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 EMData * | preprocess_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 EMData * | finish (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 | ||||
static Reconstructor * | NEW () | |||
Factory incorporation uses the pointer of this function. | ||||
Static Public Attributes | ||||
static 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 | ||||
FourierPixelInserter3D * | inserter | |||
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. | ||||
FourierReconstructor & | operator= (const FourierReconstructor &) | |||
Disallow assignment. |
The client creates a Fourier reconstructor to insert real images into a 3D volume. The return image is a real space image
This reconstructor is based on EMAN1's Fourier reconstructor with a handful of modifications including 1. - Fourier ring correlation (FRC) as opposed to the mean phase residual is used to estimate slice quality. The FRC of the slice in the 3D volume is determined - but the slice is removed from the 3D volume before doing this so the score reflects the extent to which the slice agrees with the contribution of the other slices in the 3D volume. The FRC is converted to SNR using the relationship described by Penczek ( Three-dimensional spectral signal to noise ratio for a class of reconstruction algorithms, JSB 2002 138 (24-46) ) FRC = S/(sqrt(S+N1)sqrt(S+N2)) Where N1 is the noise in the slice of the 3D volume and N2 is the noise in the image slice being inserted. We make the assumption that the noise in the 3D volume is 0 (N1=0) to get FRC^2 = SNR/(1+SNR) which gives a spectral SNR plot - we then divide each SNR value by the number of particles in the class average (seeing as SNR should scale linearly with the number of particles) to get the estimated SNR per contributing particle in this class average. If the particles that have been averaged are not homogenous this score should be low etc. The scaled SNR curve is then converted back to a FRC curve and integrated. This integral is the similarity metric, and depends on how far information extends to in Fourier space - typical values range from 0.05 to 0.2, but can vary substantially depending on the data.
2 - Uses half of the memory used by EMAN1's equivalent reconstruction algorithm
Reconstructor* r = Factory<Reconstructor>::get("fourier", params); r->setup(); for k in 0:num_iterations-1 // First do a round of slice quality (metric) determination - only possible if a 3D volume has // already been generated (k>0) if ( k > 0 ) // Determine the agreement of the slices with the previous reconstructed volume (in memory) for i in 0:num_slices-1 r->determine_slice_agreement(image[i], image[i].euler_orientation); // Insert the slices into the 3D volume // Will decide not to insert the slice if the its "quality" is not good enough for i in 0:num_slices-1 int failure = r->insert_slice(image[i], image[i].euler_orientation); if ( failure ) cout << "Slice was not inserted due to poor quality" << endl; // Get the resulting volume EMData* result = r->finish(); result->write_image("threed.mrc");
Definition at line 363 of file reconstructor.h.
EMAN::FourierReconstructor::FourierReconstructor | ( | ) | [inline] |
Default constructor calls load_default_settings().
Definition at line 369 of file reconstructor.h.
References load_default_settings().
Referenced by NEW().
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.
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
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 |
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_default(), NullPointerException, preprocess_slice(), EMAN::EMData::set_attr(), EMAN::Transform::set_mirror(), EMAN::Transform::set_scale(), and EMAN::Transform::set_trans().
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.
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 EMAN::dot(), dt, EMAN::EMData::get_data(), EMAN::Symmetry3D::get_symmetries(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, pixel_at(), power(), EMAN::EMData::set_attr(), sqrt(), EMAN::ReconstructorVolumeData::subnx, EMAN::ReconstructorVolumeData::subny, EMAN::ReconstructorVolumeData::subnz, 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.
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::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, 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.
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 |
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::image, EMAN::ReconstructorVolumeData::normalize_threed(), EMAN::FactoryBase::params, EMAN::EMData::process_inplace(), EMAN::Dict::set_default(), EMAN::ReconstructorVolumeData::tmp_data, EMAN::EMData::update(), and EMAN::EMData::write_image().
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 EMAN::ReconstructorVolumeData::image, inserter, and EMAN::ReconstructorVolumeData::tmp_data.
Referenced by ~FourierReconstructor().
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.
References NAME.
00442 { 00443 return NAME; 00444 }
virtual TypeDict EMAN::FourierReconstructor::get_param_types | ( | ) | const [inline, virtual] |
Get the parameter types of this object.
Implements EMAN::FactoryBase.
Definition at line 464 of file reconstructor.h.
References EMAN::EMObject::BOOL, EMAN::EMObject::INTARRAY, EMAN::TypeDict::put(), and EMAN::EMObject::STRING.
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.
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 |
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(), and EMAN::Transform::set_trans().
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.
References EMAN::ReconstructorVolumeData::image, inserter, and EMAN::ReconstructorVolumeData::tmp_data.
Referenced by FourierReconstructor().
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::EMData::get_data(), EMAN::ReconstructorVolumeData::image, EMAN::FourierPixelInserter3D::init(), inserter, EMAN::FactoryBase::params, and EMAN::ReconstructorVolumeData::tmp_data.
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 }
static Reconstructor* EMAN::FourierReconstructor::NEW | ( | ) | [inline, static] |
Factory incorporation uses the pointer of this function.
Definition at line 456 of file reconstructor.h.
References FourierReconstructor().
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.
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 EMAN::Util::fast_exp(), EMAN::EMData::get_complex_index(), EMAN::EMData::get_complex_index_fast(), EMAN::EMData::get_data(), EMAN::Util::hypot3sq(), EMAN::EMConsts::I2G, EMAN::ReconstructorVolumeData::image, norm(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::nx2, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::ny2, EMAN::ReconstructorVolumeData::nz, EMAN::ReconstructorVolumeData::nz2, rdata, EMAN::ReconstructorVolumeData::subnx, EMAN::ReconstructorVolumeData::subx0, EMAN::ReconstructorVolumeData::suby0, EMAN::ReconstructorVolumeData::subz0, and EMAN::ReconstructorVolumeData::tmp_data.
Referenced by do_compare_slice_work().
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.
slice | the slice to be prepocessed | |
t | transform to be used for insertion |
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
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(), EMAN::ReconstructorVolumeData::image, ImageDimensionException, is_fftodd(), load_inserter(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::nx2, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::ny2, EMAN::ReconstructorVolumeData::nz, EMAN::ReconstructorVolumeData::nz2, EMAN::FactoryBase::params, EMAN::EMData::set_attr(), EMAN::EMData::set_complex(), EMAN::Dict::set_default(), EMAN::EMData::set_fftodd(), EMAN::EMData::set_ri(), EMAN::EMData::set_size(), sub(), EMAN::ReconstructorVolumeData::subnx, EMAN::ReconstructorVolumeData::subny, EMAN::ReconstructorVolumeData::subnz, EMAN::ReconstructorVolumeData::subx0, EMAN::ReconstructorVolumeData::suby0, EMAN::ReconstructorVolumeData::subz0, EMAN::ReconstructorVolumeData::tmp_data, EMAN::EMData::to_zero(), and EMAN::EMData::update().
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.
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(), EMAN::ReconstructorVolumeData::image, ImageDimensionException, EMAN::EMData::is_complex(), is_fftodd(), load_inserter(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::nx2, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::ny2, EMAN::ReconstructorVolumeData::nz, EMAN::ReconstructorVolumeData::nz2, EMAN::FactoryBase::params, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::EMData::set_size(), sub(), EMAN::ReconstructorVolumeData::subnx, EMAN::ReconstructorVolumeData::subny, EMAN::ReconstructorVolumeData::subnz, EMAN::ReconstructorVolumeData::subx0, EMAN::ReconstructorVolumeData::suby0, EMAN::ReconstructorVolumeData::subz0, EMAN::ReconstructorVolumeData::tmp_data, and EMAN::EMData::to_value().
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 }
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(), load_default_settings(), and load_inserter().
const string FourierReconstructor::NAME = "fourier" [static] |