#include <aligner.h>
Inheritance diagram for EMAN::Refine3DAlignerGrid:
Public Member Functions | |
virtual EMData * | align (EMData *this_img, EMData *to_img, const string &cmp_name="sqeuclidean", const Dict &cmp_params=Dict()) const |
To align 'this_img' with another image passed in through its parameters. | |
virtual EMData * | align (EMData *this_img, EMData *to_img) const |
virtual string | get_name () const |
Get the Aligner's name. | |
virtual string | get_desc () const |
virtual TypeDict | get_param_types () const |
Static Public Member Functions | |
static Aligner * | NEW () |
Static Public Attributes | |
static const string | NAME = "refine_3d_grid" |
Refines a preliminary 3D alignment using a sampling grid. This is a port from tomohunter, but the az sampling scheme is altered cuch that the points on the sphere are equidistant (Improves speed several hundered times). The distance between the points on the sphere is 'delta' and the range(distance from the pole, 0,0,0 position) is given as 'range'. IN general this refinement scheme is a bit slower than the Quaternion Simplex aligner, but perfroms better in the presence of noise(according to current dogma).
xform.align3d | The Transform from the previous course alignment. If unspecified the identity matrix is used | |
delta | The angluar distance bewteen points on the spehere used in the grid | |
range | The size of the grid. Measured in angluar distance from the north pole | |
dotrans | Do a translational search | |
search | The maximum length of the detectable translational shift - if you supply this parameter you can not supply the maxshiftx, maxshifty or maxshiftz parameters. Each approach is mutually exclusive | |
searchx | The maximum length of the detectable translational shift in the x direction- if you supply this parameter you can not supply the maxshift parameters | |
searchy | The maximum length of the detectable translational shift in the y direction- if you supply this parameter you can not supply the maxshift parameters | |
searchz | The maximum length of the detectable translational shift in the z direction- if you supply this parameter you can not supply the maxshift parameters | |
verbose | Turn this on to have useful information printed to standard out |
Definition at line 916 of file aligner.h.
virtual EMData* EMAN::Refine3DAlignerGrid::align | ( | EMData * | this_img, | |
EMData * | to_img | |||
) | const [inline, virtual] |
Implements EMAN::Aligner.
Definition at line 922 of file aligner.h.
References align().
00923 { 00924 return align(this_img, to_img, "sqeuclidean", Dict()); 00925 }
EMData * Refine3DAlignerGrid::align | ( | EMData * | this_img, | |
EMData * | to_img, | |||
const string & | cmp_name = "sqeuclidean" , |
|||
const Dict & | cmp_params = Dict() | |||
) | const [virtual] |
To align 'this_img' with another image passed in through its parameters.
The alignment uses a user-given comparison method to compare the two images. If none is given, a default one is used.
this_img | The image to be compared. | |
to_img | 'this_img" is aligned with 'to_img'. | |
cmp_name | The comparison method to compare the two images. | |
cmp_params | The parameter dictionary for comparison method. |
Implements EMAN::Aligner.
Definition at line 1734 of file aligner.cpp.
References EMAN::EMData::calc_ccf(), EMAN::EMData::calc_max_location_wrap(), calc_max_location_wrap_cuda(), EMAN::EMData::cmp(), data, EMAN::EMData::do_fft(), EMAN::EMData::get_ndim(), get_stats_cuda(), EMAN::EMData::get_value_at_wrap(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Dict::has_key(), ImageDimensionException, InvalidParameterException, EMAN::Aligner::params, phi, EMAN::EMData::process(), EMAN::EMData::process_inplace(), EMAN::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::Transform::set_trans(), sqrt(), and t.
Referenced by align().
01736 { 01737 if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) { 01738 throw ImageDimensionException("This aligner only works for 3D images"); 01739 } 01740 01741 Transform* t; 01742 if (params.has_key("xform.align3d") ) { 01743 // Unlike the 2d refine aligner, this class doesn't require the starting transform's 01744 // parameters to form the starting guess. Instead the Transform itself 01745 // is perturbed carefully (using quaternion rotation) to overcome problems that arise 01746 // when you use orthogonally-based Euler angles 01747 t = params["xform.align3d"]; 01748 }else { 01749 t = new Transform(); // is the identity 01750 } 01751 01752 int searchx = 0; 01753 int searchy = 0; 01754 int searchz = 0; 01755 01756 bool dotrans = params.set_default("dotrans",1); 01757 if (params.has_key("search")) { 01758 vector<string> check; 01759 check.push_back("searchx"); 01760 check.push_back("searchy"); 01761 check.push_back("searchz"); 01762 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) { 01763 if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters"); 01764 } 01765 int search = params["search"]; 01766 searchx = search; 01767 searchy = search; 01768 searchz = search; 01769 } else { 01770 searchx = params.set_default("searchx",3); 01771 searchy = params.set_default("searchy",3); 01772 searchz = params.set_default("searchz",3); 01773 } 01774 01775 float delta = params.set_default("delta",1.0f); 01776 float range = params.set_default("range",10.0f); 01777 bool verbose = params.set_default("verbose",false); 01778 01779 bool tomography = (cmp_name == "ccc.tomo") ? 1 : 0; 01780 EMData * tofft = 0; 01781 if(dotrans || tomography){ 01782 tofft = to->do_fft(); 01783 } 01784 01785 #ifdef EMAN2_USING_CUDA 01786 if(EMData::usecuda == 1) { 01787 if(!this_img->getcudarodata()) this_img->copy_to_cudaro(); // This is safer 01788 if(!to->getcudarwdata()) to->copy_to_cuda(); 01789 if(to->getcudarwdata()){if(tofft) tofft->copy_to_cuda();} 01790 } 01791 #endif 01792 01793 Dict d; 01794 d["type"] = "eman"; // d is used in the loop below 01795 Dict best; 01796 best["score"] = 0.0f; 01797 bool use_cpu = true; 01798 Transform tran = Transform(); 01799 for ( float alt = 0; alt < range; alt += delta) { 01800 // now compute a sane az step size 01801 float saz = 360; 01802 if(alt != 0) saz = delta/sin(alt*M_PI/180.0f); // This gives consistent az step sizes(arc lengths) 01803 for ( float az = 0; az < 360; az += saz ){ 01804 if (verbose) { 01805 cout << "Trying angle alt " << alt << " az " << az << endl; 01806 } 01807 // account for any changes in az 01808 for( float phi = -range-az; phi < range-az; phi += delta ) { 01809 d["alt"] = alt; 01810 d["phi"] = phi; 01811 d["az"] = az; 01812 Transform tr(d); 01813 tr = tr*(*t); // compose transforms, this moves to the pole (aprox) 01814 01815 EMData* transformed = this_img->process("xform",Dict("transform",&tr)); 01816 01817 //need to do things a bit diffrent if we want to compare two tomos 01818 float score = 0.0f; 01819 if(dotrans || tomography){ 01820 EMData* ccf = transformed->calc_ccf(tofft); 01821 #ifdef EMAN2_USING_CUDA 01822 if(to->getcudarwdata()){ 01823 use_cpu = false; 01824 CudaPeakInfo* data = calc_max_location_wrap_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz); 01825 tran.set_trans((float)-data->px, (float)-data->py, (float)-data->pz); 01826 //CudaPeakInfoFloat* data = calc_max_location_wrap_intp_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz); 01827 //tran.set_trans(-data->xintp, -data->yintp, -data->zintp); 01828 tr = tran*tr; // to reflect the fact that we have done a rotation first and THEN a transformation 01829 if (tomography) { 01830 float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize()); 01831 score = -(data->peak - stats.x)/sqrt(stats.y); // Normalize, this is better than calling the norm processor since we only need to normalize one point 01832 } else { 01833 score = -data->peak; 01834 } 01835 delete data; 01836 } 01837 #endif 01838 if(use_cpu){ 01839 if(tomography) ccf->process_inplace("normalize"); 01840 //vector<float> fpoint = ccf->calc_max_location_wrap_intp(searchx,searchy,searchz); 01841 //tran.set_trans(-fpoint[0], -fpoint[1], -fpoint[2]); 01842 //score = -fpoint[3]; 01843 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz); 01844 tran.set_trans((float)-point[0], (float)-point[1], (float)-point[2]); 01845 score = -ccf->get_value_at_wrap(point[0], point[1], point[2]); 01846 tr = tran*tr;// to reflect the fact that we have done a rotation first and THEN a transformation 01847 01848 } 01849 delete ccf; ccf =0; 01850 delete transformed; transformed = 0;// this is to stop a mem leak 01851 } 01852 01853 if(!tomography){ 01854 if(!transformed) transformed = this_img->process("xform",Dict("transform",&tr)); // we are returning a map 01855 score = transformed->cmp(cmp_name,to,cmp_params); //this is not very efficient as it creates a new cmp object for each iteration 01856 delete transformed; transformed = 0;// this is to stop a mem leak 01857 } 01858 01859 if(score < float(best["score"])) { 01860 best["score"] = score; 01861 best["xform.align3d"] = &tr; // I wonder if this will cause a mem leak? 01862 } 01863 } 01864 } 01865 } 01866 01867 if(tofft) {delete tofft; tofft = 0;} 01868 01869 //make aligned map; 01870 EMData* best_match = this_img->process("xform",Dict("transform", best["xform.align3d"])); // we are returning a map 01871 best_match->set_attr("xform.align3d", best["xform.align3d"]); 01872 best_match->set_attr("score", float(best["score"])); 01873 01874 //Free up resources (for an expensive opperation like this move data to and from device is a small % of time) 01875 #ifdef EMAN2_USING_CUDA 01876 to->copy_from_device(); 01877 this_img->ro_free(); 01878 // May move best_match to device? 01879 #endif 01880 01881 //debug.... 01882 //Transform* zz = best_match->get_attr("xform.align3d"); 01883 //Vec3f zzz = zz->get_trans(); 01884 //cout << "x " << float(zzz[0]) << " y " << float(zzz[1]) << " z " << float(zzz[2]) << endl; 01885 01886 return best_match; 01887 01888 }
virtual string EMAN::Refine3DAlignerGrid::get_desc | ( | ) | const [inline, virtual] |
Implements EMAN::Aligner.
Definition at line 932 of file aligner.h.
00933 { 00934 return "Refines a preliminary 3D alignment using a simplex algorithm. Subpixel precision."; 00935 }
virtual string EMAN::Refine3DAlignerGrid::get_name | ( | ) | const [inline, virtual] |
virtual TypeDict EMAN::Refine3DAlignerGrid::get_param_types | ( | ) | const [inline, virtual] |
Implements EMAN::Aligner.
Definition at line 942 of file aligner.h.
References EMAN::EMObject::BOOL, EMAN::EMObject::FLOAT, EMAN::EMObject::INT, EMAN::TypeDict::put(), and EMAN::EMObject::TRANSFORM.
00943 { 00944 TypeDict d; 00945 d.put("xform.align3d", EMObject::TRANSFORM,"The Transform storing the starting guess. If unspecified the identity matrix is used"); 00946 d.put("delta", EMObject::FLOAT, "The angular step size. Default is 1." ); 00947 d.put("range", EMObject::FLOAT, "The angular range size. Default is 10." ); 00948 d.put("dotrans", EMObject::BOOL,"Do a translational search. Default is True(1)"); 00949 d.put("search", EMObject::INT,"The maximum length of the detectable translational shift - if you supply this parameter you can not supply the maxshiftx, maxshifty or maxshiftz parameters. Each approach is mutually exclusive."); 00950 d.put("searchx", EMObject::INT,"The maximum length of the detectable translational shift in the x direction- if you supply this parameter you can not supply the maxshift parameters. Default is 3."); 00951 d.put("searchy", EMObject::INT,"The maximum length of the detectable translational shift in the y direction- if you supply this parameter you can not supply the maxshift parameters. Default is 3."); 00952 d.put("searchz", EMObject::INT,"The maximum length of the detectable translational shift in the z direction- if you supply this parameter you can not supply the maxshift parameters. Default is 3"); 00953 d.put("maxshift", EMObject::INT,"Maximum translation in pixels in any direction. If the solution yields a shift beyond this value in any direction, then the refinement is judged a failure and the original alignment is used as the solution."); 00954 d.put("verbose", EMObject::BOOL,"Turn this on to have useful information printed to standard out."); 00955 return d; 00956 }
static Aligner* EMAN::Refine3DAlignerGrid::NEW | ( | ) | [inline, static] |
const string Refine3DAlignerGrid::NAME = "refine_3d_grid" [static] |