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

static AlignerNEW ()

Static Public Attributes

static 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 1384 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 1390 of file aligner.h.

References align().

01391                         {
01392                                 return align(this_img, to_img, "sqeuclidean", Dict());
01393                         }

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

02265 {
02266         if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
02267                 throw ImageDimensionException("This aligner only works for 3D images");
02268         }
02269 
02270         Transform* t;
02271         if (params.has_key("xform.align3d") ) {
02272                 // Unlike the 2d refine aligner, this class doesn't require the starting transform's
02273                 // parameters to form the starting guess. Instead the Transform itself
02274                 // is perturbed carefully (using quaternion rotation) to overcome problems that arise
02275                 // when you use orthogonally-based Euler angles
02276                 t = params["xform.align3d"];
02277         }else {
02278                 t = new Transform(); // is the identity
02279         }
02280 
02281         int searchx = 0;
02282         int searchy = 0;
02283         int searchz = 0;
02284         bool dotrans = params.set_default("dotrans",1);
02285         if (params.has_key("search")) {
02286                 vector<string> check;
02287                 check.push_back("searchx");
02288                 check.push_back("searchy");
02289                 check.push_back("searchz");
02290                 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
02291                         if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
02292                 }
02293                 int search  = params["search"];
02294                 searchx = search;
02295                 searchy = search;
02296                 searchz = search;
02297         } else {
02298                 searchx = params.set_default("searchx",3);
02299                 searchy = params.set_default("searchy",3);
02300                 searchz = params.set_default("searchz",3);
02301         }       
02302 
02303         float delta = params.set_default("delta",1.0f);
02304         float range = params.set_default("range",10.0f);
02305         bool verbose = params.set_default("verbose",false);
02306         
02307         bool tomography = (cmp_name == "ccc.tomo") ? 1 : 0;
02308         EMData * tofft = 0;
02309         if(dotrans || tomography){
02310                 tofft = to->do_fft();
02311         }
02312 
02313 #ifdef EMAN2_USING_CUDA 
02314         if(EMData::usecuda == 1) {
02315                 if(!this_img->getcudarodata()) this_img->copy_to_cudaro(); // This is safer
02316                 if(!to->getcudarwdata()) to->copy_to_cuda();
02317                 if(to->getcudarwdata()){if(tofft) tofft->copy_to_cuda();}
02318         }
02319 #endif
02320 
02321         Dict d;
02322         d["type"] = "eman"; // d is used in the loop below
02323         Dict best;
02324 //      best["score"] = numeric_limits<float>::infinity();
02325         best["score"] = 1.0e37;
02326         bool use_cpu = true;
02327         Transform tran = Transform();
02328         Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
02329         
02330         for ( float alt = 0; alt < range; alt += delta) {
02331                 // now compute a sane az step size 
02332                 float saz = 360;
02333                 if(alt != 0) saz = delta/sin(alt*M_PI/180.0f); // This gives consistent az step sizes(arc lengths)
02334                 for ( float az = 0; az < 360; az += saz ){
02335                         if (verbose) {
02336                                 cout << "Trying angle alt " << alt << " az " << az << endl;
02337                         }
02338                         // account for any changes in az
02339                         for( float phi = -range-az; phi < range-az; phi += delta ) {
02340                                 d["alt"] = alt;
02341                                 d["phi"] = phi; 
02342                                 d["az"] = az;
02343                                 Transform tr(d);
02344                                 tr = tr*(*t);   // compose transforms, this moves to the pole (aprox)
02345                                 
02346                                 //EMData* transformed = this_img->process("xform",Dict("transform",&tr));
02347                                 EMData* transformed;
02348                                 transformed = this_img->process("xform",Dict("transform",&tr));
02349                                 
02350                                 //need to do things a bit diffrent if we want to compare two tomos
02351                                 float score = 0.0f;
02352                                 if(dotrans || tomography){
02353                                         EMData* ccf = transformed->calc_ccf(tofft);
02354 #ifdef EMAN2_USING_CUDA 
02355                                         if(EMData::usecuda == 1){
02356                                                 use_cpu = false;
02357                                                 CudaPeakInfo* data = calc_max_location_wrap_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
02358                                                 tran.set_trans((float)-data->px, (float)-data->py, (float)-data->pz);
02359                                                 //CudaPeakInfoFloat* data = calc_max_location_wrap_intp_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
02360                                                 //tran.set_trans(-data->xintp, -data->yintp, -data->zintp);
02361                                                 tr = tran*tr; // to reflect the fact that we have done a rotation first and THEN a transformation
02362                                                 if (tomography) {
02363                                                         float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
02364                                                         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
02365                                                 } else {
02366                                                         score = -data->peak;
02367                                                 }
02368                                                 delete data;
02369                                         }
02370 #endif
02371                                         if(use_cpu){
02372                                                 if(tomography) ccf->process_inplace("normalize");
02373                                                 //vector<float> fpoint = ccf->calc_max_location_wrap_intp(searchx,searchy,searchz);
02374                                                 //tran.set_trans(-fpoint[0], -fpoint[1], -fpoint[2]);
02375                                                 //score = -fpoint[3];
02376                                                 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
02377                                                 tran.set_trans((float)-point[0], (float)-point[1], (float)-point[2]);
02378                                                 score = -ccf->get_value_at_wrap(point[0], point[1], point[2]);
02379                                                 tr = tran*tr;// to reflect the fact that we have done a rotation first and THEN a transformation
02380                                                 
02381                                         }
02382                                         delete ccf; ccf =0;
02383                                         delete transformed; transformed = 0;// this is to stop a mem leak
02384                                 }
02385 
02386                                 if(!tomography){
02387                                         if(!transformed) transformed = this_img->process("xform",Dict("transform",&tr)); // we are returning a map
02388                                         score = c->cmp(to,transformed);
02389                                         delete transformed; transformed = 0;// this is to stop a mem leak
02390                                 }
02391                                 
02392                                 if(score < float(best["score"])) {
02393 //                                      printf("%f\n",score);
02394                                         best["score"] = score;
02395                                         best["xform.align3d"] = &tr; // I wonder if this will cause a mem leak?
02396                                 }       
02397                         }
02398                 }
02399         }
02400 
02401         if(tofft) {delete tofft; tofft = 0;}
02402         if (c != 0) delete c;
02403         
02404         //make aligned map;
02405         EMData* best_match = this_img->process("xform",Dict("transform", best["xform.align3d"])); // we are returning a map
02406         best_match->set_attr("xform.align3d", best["xform.align3d"]);
02407         best_match->set_attr("score", float(best["score"]));
02408         
02409         return best_match;
02410         
02411 }

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

Implements EMAN::Aligner.

Definition at line 1400 of file aligner.h.

01401                         {
01402                                 return "Refines a preliminary 3D alignment using a simplex algorithm. Subpixel precision.";
01403                         }

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

References NAME.

01396                         {
01397                                 return NAME;
01398                         }

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

Implements EMAN::Aligner.

Definition at line 1410 of file aligner.h.

References EMAN::EMObject::BOOL, EMAN::EMObject::FLOAT, EMAN::EMObject::INT, EMAN::TypeDict::put(), and EMAN::EMObject::TRANSFORM.

01411                         {
01412                                 TypeDict d;
01413                                 d.put("xform.align3d", EMObject::TRANSFORM,"The Transform storing the starting guess. If unspecified the identity matrix is used");
01414                                 d.put("delta", EMObject::FLOAT, "The angular step size. Default is 1." );
01415                                 d.put("range", EMObject::FLOAT, "The angular range size. Default is 10." );
01416                                 d.put("dotrans", EMObject::BOOL,"Do a translational search. Default is True(1)");
01417                                 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.");
01418                                 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.");
01419                                 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.");
01420                                 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");
01421                                 d.put("verbose", EMObject::BOOL,"Turn this on to have useful information printed to standard out.");
01422                                 return d;
01423                         }

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

Definition at line 1405 of file aligner.h.

01406                         {
01407                                 return new Refine3DAlignerGrid();
01408                         }


Member Data Documentation

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

Definition at line 1425 of file aligner.h.

Referenced by get_name().


The documentation for this class was generated from the following files:
Generated on Fri Aug 10 16:34:48 2012 for EMAN2 by  doxygen 1.4.7