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 1012 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 1018 of file aligner.h.

References align().

01019                         {
01020                                 return align(this_img, to_img, "sqeuclidean", Dict());
01021                         }

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 1931 of file aligner.cpp.

References EMAN::EMData::calc_ccf(), EMAN::EMData::calc_max_location_wrap(), calc_max_location_wrap_cuda(), EMAN::Cmp::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.

01933 {
01934         if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
01935                 throw ImageDimensionException("This aligner only works for 3D images");
01936         }
01937 
01938         Transform* t;
01939         if (params.has_key("xform.align3d") ) {
01940                 // Unlike the 2d refine aligner, this class doesn't require the starting transform's
01941                 // parameters to form the starting guess. Instead the Transform itself
01942                 // is perturbed carefully (using quaternion rotation) to overcome problems that arise
01943                 // when you use orthogonally-based Euler angles
01944                 t = params["xform.align3d"];
01945         }else {
01946                 t = new Transform(); // is the identity
01947         }
01948         
01949         int searchx = 0;
01950         int searchy = 0;
01951         int searchz = 0;
01952         
01953         bool dotrans = params.set_default("dotrans",1);
01954         if (params.has_key("search")) {
01955                 vector<string> check;
01956                 check.push_back("searchx");
01957                 check.push_back("searchy");
01958                 check.push_back("searchz");
01959                 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
01960                         if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
01961                 }
01962                 int search  = params["search"];
01963                 searchx = search;
01964                 searchy = search;
01965                 searchz = search;
01966         } else {
01967                 searchx = params.set_default("searchx",3);
01968                 searchy = params.set_default("searchy",3);
01969                 searchz = params.set_default("searchz",3);
01970         }       
01971         
01972         float delta = params.set_default("delta",1.0f);
01973         float range = params.set_default("range",10.0f);
01974         bool verbose = params.set_default("verbose",false);
01975         
01976         bool tomography = (cmp_name == "ccc.tomo") ? 1 : 0;
01977         EMData * tofft = 0;
01978         if(dotrans || tomography){
01979                 tofft = to->do_fft();
01980         }
01981         
01982 #ifdef EMAN2_USING_CUDA 
01983         if(EMData::usecuda == 1) {
01984                 if(!this_img->getcudarodata()) this_img->copy_to_cudaro(); // This is safer
01985                 if(!to->getcudarwdata()) to->copy_to_cuda();
01986                 if(to->getcudarwdata()){if(tofft) tofft->copy_to_cuda();}
01987         }
01988 #endif
01989 
01990         Dict d;
01991         d["type"] = "eman"; // d is used in the loop below
01992         Dict best;
01993         best["score"] = 0.0f;
01994         bool use_cpu = true;
01995         Transform tran = Transform();
01996         Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
01997         for ( float alt = 0; alt < range; alt += delta) {
01998                 // now compute a sane az step size 
01999                 float saz = 360;
02000                 if(alt != 0) saz = delta/sin(alt*M_PI/180.0f); // This gives consistent az step sizes(arc lengths)
02001                 for ( float az = 0; az < 360; az += saz ){
02002                         if (verbose) {
02003                                 cout << "Trying angle alt " << alt << " az " << az << endl;
02004                         }
02005                         // account for any changes in az
02006                         for( float phi = -range-az; phi < range-az; phi += delta ) {
02007                                 d["alt"] = alt;
02008                                 d["phi"] = phi; 
02009                                 d["az"] = az;
02010                                 Transform tr(d);
02011                                 tr = tr*(*t);   // compose transforms, this moves to the pole (aprox)
02012                                 
02013                                 EMData* transformed = this_img->process("xform",Dict("transform",&tr));
02014                                 
02015                                 //need to do things a bit diffrent if we want to compare two tomos
02016                                 float score = 0.0f;
02017                                 if(dotrans || tomography){
02018                                         EMData* ccf = transformed->calc_ccf(tofft);
02019 #ifdef EMAN2_USING_CUDA 
02020                                         if(to->getcudarwdata()){
02021                                                 use_cpu = false;
02022                                                 CudaPeakInfo* data = calc_max_location_wrap_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
02023                                                 tran.set_trans((float)-data->px, (float)-data->py, (float)-data->pz);
02024                                                 //CudaPeakInfoFloat* data = calc_max_location_wrap_intp_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
02025                                                 //tran.set_trans(-data->xintp, -data->yintp, -data->zintp);
02026                                                 tr = tran*tr; // to reflect the fact that we have done a rotation first and THEN a transformation
02027                                                 if (tomography) {
02028                                                         float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
02029                                                         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
02030                                                 } else {
02031                                                         score = -data->peak;
02032                                                 }
02033                                                 delete data;
02034                                         }
02035 #endif
02036                                         if(use_cpu){
02037                                                 if(tomography) ccf->process_inplace("normalize");
02038                                                 //vector<float> fpoint = ccf->calc_max_location_wrap_intp(searchx,searchy,searchz);
02039                                                 //tran.set_trans(-fpoint[0], -fpoint[1], -fpoint[2]);
02040                                                 //score = -fpoint[3];
02041                                                 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
02042                                                 tran.set_trans((float)-point[0], (float)-point[1], (float)-point[2]);
02043                                                 score = -ccf->get_value_at_wrap(point[0], point[1], point[2]);
02044                                                 tr = tran*tr;// to reflect the fact that we have done a rotation first and THEN a transformation
02045                                                 
02046                                         }
02047                                         delete ccf; ccf =0;
02048                                         delete transformed; transformed = 0;// this is to stop a mem leak
02049                                 }
02050 
02051                                 if(!tomography){
02052                                         if(!transformed) transformed = this_img->process("xform",Dict("transform",&tr)); // we are returning a map
02053                                         score = c->cmp(to,transformed);
02054                                         delete transformed; transformed = 0;// this is to stop a mem leak
02055                                 }
02056                                 
02057                                 if(score < float(best["score"])) {
02058                                         best["score"] = score;
02059                                         best["xform.align3d"] = &tr; // I wonder if this will cause a mem leak?
02060                                 }       
02061                         }
02062                 }
02063         }
02064 
02065         if(tofft) {delete tofft; tofft = 0;}
02066         
02067         //make aligned map;
02068         EMData* best_match = this_img->process("xform",Dict("transform", best["xform.align3d"])); // we are returning a map
02069         best_match->set_attr("xform.align3d", best["xform.align3d"]);
02070         best_match->set_attr("score", float(best["score"]));
02071         
02072         //Free up resources (for an expensive opperation like this move data to and from device is a small % of time)
02073         #ifdef EMAN2_USING_CUDA
02074                 to->copy_from_device();
02075                 this_img->ro_free();
02076                 // May move best_match to device?
02077         #endif
02078         
02079         //debug....
02080         //Transform* zz = best_match->get_attr("xform.align3d");
02081         //Vec3f zzz = zz->get_trans();
02082         //cout << "x " << float(zzz[0]) << " y " << float(zzz[1]) << " z " << float(zzz[2]) << endl;
02083         
02084         return best_match;
02085         
02086 }

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

Implements EMAN::Aligner.

Definition at line 1028 of file aligner.h.

01029                         {
01030                                 return "Refines a preliminary 3D alignment using a simplex algorithm. Subpixel precision.";
01031                         }

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 1023 of file aligner.h.

01024                         {
01025                                 return NAME;
01026                         }

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

Implements EMAN::Aligner.

Definition at line 1038 of file aligner.h.

References EMAN::TypeDict::put().

01039                         {
01040                                 TypeDict d;
01041                                 d.put("xform.align3d", EMObject::TRANSFORM,"The Transform storing the starting guess. If unspecified the identity matrix is used");
01042                                 d.put("delta", EMObject::FLOAT, "The angular step size. Default is 1." );
01043                                 d.put("range", EMObject::FLOAT, "The angular range size. Default is 10." );
01044                                 d.put("dotrans", EMObject::BOOL,"Do a translational search. Default is True(1)");
01045                                 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.");
01046                                 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.");
01047                                 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.");
01048                                 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");
01049                                 d.put("verbose", EMObject::BOOL,"Turn this on to have useful information printed to standard out.");
01050                                 return d;
01051                         }

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

Definition at line 1033 of file aligner.h.

01034                         {
01035                                 return new Refine3DAlignerGrid();
01036                         }


Member Data Documentation

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

Definition at line 77 of file aligner.cpp.


The documentation for this class was generated from the following files:
Generated on Thu Nov 17 12:45:17 2011 for EMAN2 by  doxygen 1.3.9.1