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

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

Refine alignment. More...

#include <aligner.h>

Inheritance diagram for EMAN::Refine3DAlignerGrid:

Inheritance graph
[legend]
Collaboration diagram for EMAN::Refine3DAlignerGrid:

Collaboration graph
[legend]
List of all members.

Public Member Functions

virtual EMDataalign (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 EMDataalign (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

AlignerNEW ()

Static Public Attributes

const string NAME = "refine_3d_grid"

Detailed Description

Refine alignment.

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

Parameters:
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
Author:
John Flanagan
Date:
Mar 2011

Definition at line 916 of file aligner.h.


Member Function Documentation

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.

Parameters:
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.
Returns:
The aligned image.

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::Dict::end(), 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, CudaPeakInfo::peak, phi, EMAN::EMData::process(), EMAN::EMData::process_inplace(), CudaPeakInfo::px, CudaPeakInfo::py, CudaPeakInfo::pz, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::Transform::set_trans(), sqrt(), and t.

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]
 

Get the Aligner's name.

Each Aligner is identified by a unique name.

Returns:
The Aligner's name.

Implements EMAN::Aligner.

Definition at line 927 of file aligner.h.

00928                         {
00929                                 return NAME;
00930                         }

virtual TypeDict EMAN::Refine3DAlignerGrid::get_param_types  )  const [inline, virtual]
 

Implements EMAN::Aligner.

Definition at line 942 of file aligner.h.

References EMAN::TypeDict::put().

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                         }

Aligner* EMAN::Refine3DAlignerGrid::NEW  )  [inline, static]
 

Definition at line 937 of file aligner.h.

00938                         {
00939                                 return new Refine3DAlignerGrid();
00940                         }


Member Data Documentation

const string Refine3DAlignerGrid::NAME = "refine_3d_grid" [static]
 

Definition at line 74 of file aligner.cpp.


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