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

References align().

01325                         {
01326                                 return align(this_img, to_img, "sqeuclidean", Dict());
01327                         }

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

02111 {
02112         if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
02113                 throw ImageDimensionException("This aligner only works for 3D images");
02114         }
02115 
02116         Transform* t;
02117         if (params.has_key("xform.align3d") ) {
02118                 // Unlike the 2d refine aligner, this class doesn't require the starting transform's
02119                 // parameters to form the starting guess. Instead the Transform itself
02120                 // is perturbed carefully (using quaternion rotation) to overcome problems that arise
02121                 // when you use orthogonally-based Euler angles
02122                 t = params["xform.align3d"];
02123         }else {
02124                 t = new Transform(); // is the identity
02125         }
02126 
02127         int searchx = 0;
02128         int searchy = 0;
02129         int searchz = 0;
02130         bool dotrans = params.set_default("dotrans",1);
02131         if (params.has_key("search")) {
02132                 vector<string> check;
02133                 check.push_back("searchx");
02134                 check.push_back("searchy");
02135                 check.push_back("searchz");
02136                 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
02137                         if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
02138                 }
02139                 int search  = params["search"];
02140                 searchx = search;
02141                 searchy = search;
02142                 searchz = search;
02143         } else {
02144                 searchx = params.set_default("searchx",3);
02145                 searchy = params.set_default("searchy",3);
02146                 searchz = params.set_default("searchz",3);
02147         }       
02148 
02149         float delta = params.set_default("delta",1.0f);
02150         float range = params.set_default("range",10.0f);
02151         bool verbose = params.set_default("verbose",false);
02152         bool fsrotate = params.set_default("fsrotate",0);
02153         
02154         bool tomography = (cmp_name == "ccc.tomo") ? 1 : 0;
02155         EMData * tofft = 0;
02156         if(dotrans || tomography){
02157                 tofft = to->do_fft();
02158         }
02159         EMData * this_imgfft = 0;
02160         if(fsrotate){
02161                 this_imgfft = this_img->do_fft();
02162         }
02163 
02164 #ifdef EMAN2_USING_CUDA 
02165         if(EMData::usecuda == 1) {
02166                 if(!this_img->getcudarodata()) this_img->copy_to_cudaro(); // This is safer
02167                 if(!to->getcudarwdata()) to->copy_to_cuda();
02168                 if(to->getcudarwdata()){if(tofft) tofft->copy_to_cuda();}
02169         }
02170 #endif
02171 
02172         Dict d;
02173         d["type"] = "eman"; // d is used in the loop below
02174         Dict best;
02175         best["score"] = numeric_limits<float>::infinity();
02176         bool use_cpu = true;
02177         Transform tran = Transform();
02178         Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
02179         
02180         for ( float alt = 0; alt < range; alt += delta) {
02181                 // now compute a sane az step size 
02182                 float saz = 360;
02183                 if(alt != 0) saz = delta/sin(alt*M_PI/180.0f); // This gives consistent az step sizes(arc lengths)
02184                 for ( float az = 0; az < 360; az += saz ){
02185                         if (verbose) {
02186                                 cout << "Trying angle alt " << alt << " az " << az << endl;
02187                         }
02188                         // account for any changes in az
02189                         for( float phi = -range-az; phi < range-az; phi += delta ) {
02190                                 d["alt"] = alt;
02191                                 d["phi"] = phi; 
02192                                 d["az"] = az;
02193                                 Transform tr(d);
02194                                 tr = tr*(*t);   // compose transforms, this moves to the pole (aprox)
02195                                 
02196                                 //EMData* transformed = this_img->process("xform",Dict("transform",&tr));
02197                                 EMData* transformed;
02198                                 if(fsrotate){
02199                                         transformed = this_imgfft->process("rotateinfs",Dict("transform",&tr));
02200                                 } else {
02201                                         transformed = this_img->process("xform",Dict("transform",&tr));
02202                                 }
02203                                 
02204                                 //need to do things a bit diffrent if we want to compare two tomos
02205                                 float score = 0.0f;
02206                                 if(dotrans || tomography){
02207                                         EMData* ccf = transformed->calc_ccf(tofft);
02208 #ifdef EMAN2_USING_CUDA 
02209                                         if(EMData::usecuda == 1){
02210                                                 use_cpu = false;
02211                                                 CudaPeakInfo* data = calc_max_location_wrap_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
02212                                                 tran.set_trans((float)-data->px, (float)-data->py, (float)-data->pz);
02213                                                 //CudaPeakInfoFloat* data = calc_max_location_wrap_intp_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
02214                                                 //tran.set_trans(-data->xintp, -data->yintp, -data->zintp);
02215                                                 tr = tran*tr; // to reflect the fact that we have done a rotation first and THEN a transformation
02216                                                 if (tomography) {
02217                                                         float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
02218                                                         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
02219                                                 } else {
02220                                                         score = -data->peak;
02221                                                 }
02222                                                 delete data;
02223                                         }
02224 #endif
02225                                         if(use_cpu){
02226                                                 if(tomography) ccf->process_inplace("normalize");
02227                                                 //vector<float> fpoint = ccf->calc_max_location_wrap_intp(searchx,searchy,searchz);
02228                                                 //tran.set_trans(-fpoint[0], -fpoint[1], -fpoint[2]);
02229                                                 //score = -fpoint[3];
02230                                                 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
02231                                                 tran.set_trans((float)-point[0], (float)-point[1], (float)-point[2]);
02232                                                 score = -ccf->get_value_at_wrap(point[0], point[1], point[2]);
02233                                                 tr = tran*tr;// to reflect the fact that we have done a rotation first and THEN a transformation
02234                                                 
02235                                         }
02236                                         delete ccf; ccf =0;
02237                                         delete transformed; transformed = 0;// this is to stop a mem leak
02238                                 }
02239 
02240                                 if(!tomography){
02241                                         if(!transformed) transformed = this_img->process("xform",Dict("transform",&tr)); // we are returning a map
02242                                         score = c->cmp(to,transformed);
02243                                         delete transformed; transformed = 0;// this is to stop a mem leak
02244                                 }
02245                                 
02246                                 if(score < float(best["score"])) {
02247                                         best["score"] = score;
02248                                         best["xform.align3d"] = &tr; // I wonder if this will cause a mem leak?
02249                                 }       
02250                         }
02251                 }
02252         }
02253 
02254         if(tofft) {delete tofft; tofft = 0;}
02255         if(fsrotate) {delete this_imgfft;}
02256         if (c != 0) delete c;
02257         
02258         //make aligned map;
02259         EMData* best_match = this_img->process("xform",Dict("transform", best["xform.align3d"])); // we are returning a map
02260         best_match->set_attr("xform.align3d", best["xform.align3d"]);
02261         best_match->set_attr("score", float(best["score"]));
02262         
02263         return best_match;
02264         
02265 }

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

Implements EMAN::Aligner.

Definition at line 1334 of file aligner.h.

01335                         {
01336                                 return "Refines a preliminary 3D alignment using a simplex algorithm. Subpixel precision.";
01337                         }

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

References NAME.

01330                         {
01331                                 return NAME;
01332                         }

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

Implements EMAN::Aligner.

Definition at line 1344 of file aligner.h.

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

01345                         {
01346                                 TypeDict d;
01347                                 d.put("xform.align3d", EMObject::TRANSFORM,"The Transform storing the starting guess. If unspecified the identity matrix is used");
01348                                 d.put("delta", EMObject::FLOAT, "The angular step size. Default is 1." );
01349                                 d.put("range", EMObject::FLOAT, "The angular range size. Default is 10." );
01350                                 d.put("dotrans", EMObject::BOOL,"Do a translational search. Default is True(1)");
01351                                 d.put("fsrotate", EMObject::BOOL,"Do rotations in Fourier space");
01352                                 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.");
01353                                 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.");
01354                                 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.");
01355                                 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");
01356                                 d.put("verbose", EMObject::BOOL,"Turn this on to have useful information printed to standard out.");
01357                                 return d;
01358                         }

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

Definition at line 1339 of file aligner.h.

01340                         {
01341                                 return new Refine3DAlignerGrid();
01342                         }


Member Data Documentation

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

Definition at line 1360 of file aligner.h.

Referenced by get_name().


The documentation for this class was generated from the following files:
Generated on Thu May 3 10:08:49 2012 for EMAN2 by  doxygen 1.4.7