EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes
EMAN::Refine3DAlignerGrid Class Reference

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.align3dThe Transform from the previous course alignment. If unspecified the identity matrix is used
deltaThe angluar distance bewteen points on the spehere used in the grid
rangeThe size of the grid. Measured in angluar distance from the north pole
dotransDo a translational search
searchThe 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
searchxThe maximum length of the detectable translational shift in the x direction- if you supply this parameter you can not supply the maxshift parameters
searchyThe maximum length of the detectable translational shift in the y direction- if you supply this parameter you can not supply the maxshift parameters
searchzThe maximum length of the detectable translational shift in the z direction- if you supply this parameter you can not supply the maxshift parameters
verboseTurn this on to have useful information printed to standard out
Author:
John Flanagan
Date:
Mar 2011

Definition at line 1390 of file aligner.h.


Member Function Documentation

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_imgThe image to be compared.
to_img'this_img" is aligned with 'to_img'.
cmp_nameThe comparison method to compare the two images.
cmp_paramsThe parameter dictionary for comparison method.
Returns:
The aligned image.

Implements EMAN::Aligner.

Definition at line 2284 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, 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.

Referenced by align().

{
        if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
                throw ImageDimensionException("This aligner only works for 3D images");
        }

        Transform* t;
        if (params.has_key("xform.align3d") ) {
                // Unlike the 2d refine aligner, this class doesn't require the starting transform's
                // parameters to form the starting guess. Instead the Transform itself
                // is perturbed carefully (using quaternion rotation) to overcome problems that arise
                // when you use orthogonally-based Euler angles
                t = params["xform.align3d"];
        }else {
                t = new Transform(); // is the identity
        }

        int searchx = 0;
        int searchy = 0;
        int searchz = 0;
        bool dotrans = params.set_default("dotrans",1);
        if (params.has_key("search")) {
                vector<string> check;
                check.push_back("searchx");
                check.push_back("searchy");
                check.push_back("searchz");
                for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
                        if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
                }
                int search  = params["search"];
                searchx = search;
                searchy = search;
                searchz = search;
        } else {
                searchx = params.set_default("searchx",3);
                searchy = params.set_default("searchy",3);
                searchz = params.set_default("searchz",3);
        }       

        float delta = params.set_default("delta",1.0f);
        float range = params.set_default("range",10.0f);
        bool verbose = params.set_default("verbose",false);
        
        bool tomography = (cmp_name == "ccc.tomo") ? 1 : 0;
        EMData * tofft = 0;
        if(dotrans || tomography){
                tofft = to->do_fft();
        }

#ifdef EMAN2_USING_CUDA 
        if(EMData::usecuda == 1) {
                if(!this_img->getcudarodata()) this_img->copy_to_cudaro(); // This is safer
                if(!to->getcudarwdata()) to->copy_to_cuda();
                if(to->getcudarwdata()){if(tofft) tofft->copy_to_cuda();}
        }
#endif

        Dict d;
        d["type"] = "eman"; // d is used in the loop below
        Dict best;
//      best["score"] = numeric_limits<float>::infinity();
        best["score"] = 1.0e37;
        bool use_cpu = true;
        Transform tran = Transform();
        Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
        
        for ( float alt = 0; alt < range; alt += delta) {
                // now compute a sane az step size 
                float saz = 360;
                if(alt != 0) saz = delta/sin(alt*M_PI/180.0f); // This gives consistent az step sizes(arc lengths)
                for ( float az = 0; az < 360; az += saz ){
                        if (verbose) {
                                cout << "Trying angle alt " << alt << " az " << az << endl;
                        }
                        // account for any changes in az
                        for( float phi = -range-az; phi < range-az; phi += delta ) {
                                d["alt"] = alt;
                                d["phi"] = phi; 
                                d["az"] = az;
                                Transform tr(d);
                                tr = tr*(*t);   // compose transforms, this moves to the pole (aprox)
                                
                                //EMData* transformed = this_img->process("xform",Dict("transform",&tr));
                                EMData* transformed;
                                transformed = this_img->process("xform",Dict("transform",&tr));
                                
                                //need to do things a bit diffrent if we want to compare two tomos
                                float score = 0.0f;
                                if(dotrans || tomography){
                                        EMData* ccf = transformed->calc_ccf(tofft);
#ifdef EMAN2_USING_CUDA 
                                        if(EMData::usecuda == 1){
                                                use_cpu = false;
                                                CudaPeakInfo* data = calc_max_location_wrap_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
                                                tran.set_trans((float)-data->px, (float)-data->py, (float)-data->pz);
                                                //CudaPeakInfoFloat* data = calc_max_location_wrap_intp_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
                                                //tran.set_trans(-data->xintp, -data->yintp, -data->zintp);
                                                tr = tran*tr; // to reflect the fact that we have done a rotation first and THEN a transformation
                                                if (tomography) {
                                                        float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
                                                        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
                                                } else {
                                                        score = -data->peak;
                                                }
                                                delete data;
                                        }
#endif
                                        if(use_cpu){
                                                if(tomography) ccf->process_inplace("normalize");
                                                //vector<float> fpoint = ccf->calc_max_location_wrap_intp(searchx,searchy,searchz);
                                                //tran.set_trans(-fpoint[0], -fpoint[1], -fpoint[2]);
                                                //score = -fpoint[3];
                                                IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
                                                tran.set_trans((float)-point[0], (float)-point[1], (float)-point[2]);
                                                score = -ccf->get_value_at_wrap(point[0], point[1], point[2]);
                                                tr = tran*tr;// to reflect the fact that we have done a rotation first and THEN a transformation
                                                
                                        }
                                        delete ccf; ccf =0;
                                        delete transformed; transformed = 0;// this is to stop a mem leak
                                }

                                if(!tomography){
                                        if(!transformed) transformed = this_img->process("xform",Dict("transform",&tr)); // we are returning a map
                                        score = c->cmp(to,transformed);
                                        delete transformed; transformed = 0;// this is to stop a mem leak
                                }
                                
                                if(score < float(best["score"])) {
//                                      printf("%f\n",score);
                                        best["score"] = score;
                                        best["xform.align3d"] = &tr; // I wonder if this will cause a mem leak?
                                }       
                        }
                }
        }

        if(tofft) {delete tofft; tofft = 0;}
        if (c != 0) delete c;
        
        //make aligned map;
        EMData* best_match = this_img->process("xform",Dict("transform", best["xform.align3d"])); // we are returning a map
        best_match->set_attr("xform.align3d", best["xform.align3d"]);
        best_match->set_attr("score", float(best["score"]));
        
        return best_match;
        
}
virtual EMData* EMAN::Refine3DAlignerGrid::align ( EMData this_img,
EMData to_img 
) const [inline, virtual]

Implements EMAN::Aligner.

Definition at line 1396 of file aligner.h.

References align().

                        {
                                return align(this_img, to_img, "sqeuclidean", Dict());
                        }
virtual string EMAN::Refine3DAlignerGrid::get_desc ( ) const [inline, virtual]

Implements EMAN::Aligner.

Definition at line 1406 of file aligner.h.

                        {
                                return "Refines a preliminary 3D alignment using a simplex algorithm. Subpixel precision.";
                        }
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 1401 of file aligner.h.

References NAME.

                        {
                                return NAME;
                        }
virtual TypeDict EMAN::Refine3DAlignerGrid::get_param_types ( ) const [inline, virtual]

Implements EMAN::Aligner.

Definition at line 1416 of file aligner.h.

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

                        {
                                TypeDict d;
                                d.put("xform.align3d", EMObject::TRANSFORM,"The Transform storing the starting guess. If unspecified the identity matrix is used");
                                d.put("delta", EMObject::FLOAT, "The angular step size. Default is 1." );
                                d.put("range", EMObject::FLOAT, "The angular range size. Default is 10." );
                                d.put("dotrans", EMObject::BOOL,"Do a translational search. Default is True(1)");
                                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.");
                                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.");
                                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.");
                                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");
                                d.put("verbose", EMObject::BOOL,"Turn this on to have useful information printed to standard out.");
                                return d;
                        }
static Aligner* EMAN::Refine3DAlignerGrid::NEW ( ) [inline, static]

Definition at line 1411 of file aligner.h.

                        {
                                return new Refine3DAlignerGrid();
                        }

Member Data Documentation

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

Definition at line 1431 of file aligner.h.

Referenced by get_name().


The documentation for this class was generated from the following files: