EMAN::RT3DLocalTreeAligner Class Reference

3D rotational and translational alignment using a hierarchical method with gradually decreasing downsampling in Fourier space. More...

#include <aligner.h>

Inheritance diagram for EMAN::RT3DLocalTreeAligner:

Collaboration diagram for EMAN::RT3DLocalTreeAligner:

Public Member Functions
virtual EMData *	align (EMData *this_img, EMData *to_img, const string &cmp_name="sqeuclidean", const Dict &cmp_params=Dict()) const
	See Aligner comments for more details. More...
virtual EMData *	align (EMData *this_img, EMData *to_img) const
	See Aligner comments for more details. More...
virtual vector< Dict >	xform_align_nbest (EMData *this_img, EMData *to_img, const unsigned int nsoln, const string &cmp_name, const Dict &cmp_params) const
	See Aligner comments for more details. More...
virtual string	get_name () const
	Get the Aligner's name. More...
virtual string	get_desc () const
virtual TypeDict	get_param_types () const
Public Member Functions inherited from EMAN::Aligner
virtual	~Aligner ()
virtual Dict	get_params () const
	Get the Aligner parameters in a key/value dictionary. More...
virtual void	set_params (const Dict &new_params)
	Set the Aligner parameters using a key/value dictionary. More...

Static Public Member Functions
static Aligner *	NEW ()

Static Public Attributes
static const string	NAME = "rotate_translate_3d_local_tree"

Private Member Functions
bool	testort (EMData *small_this, EMData *small_to, EMData *small_mask, EMData *small_thissq, vector< float > &sigmathisv, vector< float > &sigmatov, vector< float > &s_score, vector< float > &s_coverage, vector< Transform > &s_xform, int i, Dict &upd, Transform initxf, int maxshift) const

Additional Inherited Members
Protected Attributes inherited from EMAN::Aligner
Dict	params

Detailed Description

3D rotational and translational alignment using a hierarchical method with gradually decreasing downsampling in Fourier space.

In theory, very fast, and without need for a "refine" aligner. Comparator is ignored. Uses an inbuilt comparison. This variant uses something like FLCF (local cross correlation) for translational alignment, and should work better when the reference is masked but doing this slows things down.

Parameters

sym	The symmtery to use as the basis of the spherical sampling
verbose	Turn this on to have useful information printed to standard out

Author

Steve Ludtke

Date

April 2015

Definition at line 1959 of file aligner.h.

Member Function Documentation

align() [1/2]

virtual EMData * EMAN::RT3DLocalTreeAligner::align ( EMData *  this_img, EMData *  to_img  ) const

inlinevirtual

See Aligner comments for more details.

Implements EMAN::Aligner.

Definition at line 1968 of file aligner.h.

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

References align().

align() [2/2]

virtual EMData * EMAN::RT3DLocalTreeAligner::align ( EMData *  this_img, EMData *  to_img, const string &  cmp_name = "sqeuclidean", const Dict &  cmp_params = Dict()  ) const

virtual

See Aligner comments for more details.

Implements EMAN::Aligner.

Referenced by align().

get_desc()

virtual string EMAN::RT3DLocalTreeAligner::get_desc ( ) const

inlinevirtual

Implements EMAN::Aligner.

Definition at line 1983 of file aligner.h.

                        {
                                return "EXPERIMENATAL - 3D rotational and translational alignment using a hierarchical approach in Fourier space. Should be very fast and not require 'refine' alignment. This variant performs a locally normalized CCF for translational alignment and should thus work better when the reference is masked, but will be slower.";
                        }

get_name()

virtual string EMAN::RT3DLocalTreeAligner::get_name ( ) const

inlinevirtual

Get the Aligner's name.

Each Aligner is identified by a unique name.

Returns

The Aligner's name.

Implements EMAN::Aligner.

Definition at line 1978 of file aligner.h.

                        {
                                return NAME;
                        }

References NAME.

get_param_types()

virtual TypeDict EMAN::RT3DLocalTreeAligner::get_param_types ( ) const

inlinevirtual

Implements EMAN::Aligner.

Definition at line 1993 of file aligner.h.

                        {
                                TypeDict d;
                                d.put("maxshift", EMObject::INT,"maximum shift allowed");
                                d.put("sym", EMObject::STRING,"The symmtery to use as the basis of the spherical sampling. Default is c1 (no symmetry)");
                                d.put("sigmathis", EMObject::FLOAT,"Only Fourier voxels larger than sigma times this value will be considered");
                                d.put("sigmato", EMObject::FLOAT,"Only Fourier voxels larger than sigma times this value will be considered");
                                d.put("maxres", EMObject::FLOAT,"maximum resolution to compare");
                                d.put("minres", EMObject::FLOAT,"minimum resolution to compare");
                                d.put("maxang", EMObject::FLOAT,"maximum angle from initial rotation.");
                                d.put("initxform", EMObject::TRANSFORMARRAY,"An array of Transforms storing the starting positions.");
 
//                              d.put("initxform", EMObject::TRANSFORM,"The Transform storing the starting position. If unspecified the identity matrix is used");
                                d.put("randphi", EMObject::BOOL,"Ignore phi constraint for refine search");
                                d.put("rand180", EMObject::BOOL,"Ignore 180 rotation for refine search");
                                d.put("verbose", EMObject::BOOL,"Turn this on to have useful information printed to standard out.");
                                return d;
                        }

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

NEW()

static Aligner * EMAN::RT3DLocalTreeAligner::NEW ( )

inlinestatic

Definition at line 1988 of file aligner.h.

                        {
                                return new RT3DLocalTreeAligner();
                        }

testort()

bool RT3DLocalTreeAligner::testort ( EMData *  small_this, EMData *  small_to, EMData *  small_mask, EMData *  small_thissq, vector< float > &  sigmathisv, vector< float > &  sigmatov, vector< float > &  s_score, vector< float > &  s_coverage, vector< Transform > &  s_xform, int  i, Dict &  upd, Transform  initxf, int  maxshift  ) const

private

Definition at line 4296 of file aligner.cpp.

                                                                                                                                                                                                                                                                                                   {
        Transform t;
        Dict aap=s_xform[i].get_params("eman");
        aap["tx"]=0;
        aap["ty"]=0;
        aap["tz"]=0;
        for (Dict::const_iterator p=upd.begin(); p!=upd.end(); p++) {
                aap[p->first]=(float)aap[p->first]+(float)p->second;
        }
 
        t.set_params(aap);
        float maxang=params.set_default("maxang",-1.0);
        bool randphi=params.set_default("randphi",false);
        bool rand180=params.set_default("rand180",false);
        
        if (maxang>0){
                
                Transform tmp=initxf * t.inverse();
                
                if (randphi){
                        Dict td=t.get_params("eman");
                        Dict ti=initxf.get_params("eman");
                        td["phi"]=ti["phi"];
                        Transform tmp1=Transform(td);
                        tmp=initxf * tmp1.inverse();
                }               
                
                float r=tmp.get_params("spin")["omega"];
                if (rand180){
                        r=r>90?(180-r):r;
                }
                if(r>maxang)
                        return false;
                
        }
        
        int ny=small_this->get_ysize();
        if (params.has_key("initxform")){
                // when doing refinement, search around the given position
                t.set_trans(initxf.get_trans()*ny/(float)params["boxsize"]);
                aap=t.get_params("eman");
        }
 
        // rotate in Fourier space then use a CCF to find translation
        EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
        EMData *sttsq=small_thissq->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
//      EMData *ccf=small_to->calc_ccf(stt);
        EMData *ccf=small_to->calc_ccf_masked(stt,sttsq,small_mask);
        IntPoint ml=ccf->calc_max_location_wrap(maxshift, maxshift, maxshift);
 
 
        aap["tx"]=int(aap["tx"])+(int)ml[0];
        aap["ty"]=int(aap["ty"])+(int)ml[1];
        aap["tz"]=int(aap["tz"])+(int)ml[2];
        t.set_params(aap);
        EMData *st2=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));        // we have to do 1 slow transform here now that we have the translation
        
//      float maxres = params.set_default("maxres",0.0f);
//      float minres = params.set_default("minres",0.0f);
 
        float sim=st2->cmp("fsc.tomo.auto",small_to,Dict("sigmaimg",sigmathisv,"sigmawith",sigmatov, "pmin", 4, "pmax",ny/3));
//      float sim=st2->cmp("ccc.tomo.thresh",small_to,Dict("sigmaimg",sigmathis,"sigmawith",sigmato));
//      printf("\nTESTORT %6.1f  %6.1f  %6.1f\t%4d %4d %4d\t%1.5g\t%1.5g %d (%d)",
//              float(aap["az"]),float(aap["alt"]),float(aap["phi"]),int(aap["tx"]),int(aap["ty"]),int(aap["tz"]),sim,s_score[i],int(sim<s_score[i]),ccf->get_ysize());
 
        delete ccf;
        // If the score is better than before, we update this particular best value
        if (sim<s_score[i]) {
                s_score[i]=sim;
                s_coverage[i]=st2->get_attr("fft_overlap");
                s_xform[i]=t;
                delete stt;
        delete sttsq;
                delete st2;
                return true;
        }
        delete stt;
        delete sttsq;
        delete st2;
        return false;
}

References EMAN::Dict::begin(), EMAN::EMData::calc_ccf_masked(), EMAN::Dict::end(), EMAN::Transform::get_params(), EMAN::Transform::get_trans(), EMAN::Dict::has_key(), EMAN::Transform::inverse(), EMAN::Aligner::params, EMAN::Dict::set_default(), EMAN::Transform::set_params(), and EMAN::Transform::set_trans().

Referenced by xform_align_nbest().

xform_align_nbest()

vector< Dict > RT3DLocalTreeAligner::xform_align_nbest ( EMData *  this_img, EMData *  to_img, const unsigned int  nsoln, const string &  cmp_name, const Dict &  cmp_params  ) const

virtual

See Aligner comments for more details.

Reimplemented from EMAN::Aligner.

Definition at line 3908 of file aligner.cpp.

                                                                                                                                                                     {
        if (nrsoln == 0) throw InvalidParameterException("ERROR (RT3DTreeAligner): nsoln must be >0"); // What was the user thinking?
 
        int nsoln = nrsoln*2;
        if (nrsoln<32) nsoln=64;                // we start with at least 32 solutions, but then gradually decrease with increasing scale
 
        // !!!!!! IMPORTANT NOTE - we are inverting the order of this and to here to match convention in other aligners, to compensate
        // the Transform is inverted before being returned
        EMData *base_this;
        EMData *base_to;
        EMData *base_thissq;
        EMData *base_mask;
        
        if (this_img->is_complex()) {
                base_to=this_img->copy();
                EMData *tmp = base_to->do_ift();
                tmp->process_inplace("threshold.notzero");
                base_mask=tmp->do_fft();
                delete tmp;
        }
        else {
                base_to=this_img->do_fft();
                base_to->process_inplace("xform.phaseorigin.tocorner");
                EMData *tmp = this_img->process("threshold.notzero");
                base_mask=tmp->do_fft();
                delete tmp;
        }
 
        if (to->is_complex()) {
                base_this=to->copy();
                EMData *tmp = base_this->do_ift();
                tmp->process_inplace("math.squared");
                base_thissq = tmp->do_fft();
                delete tmp;
        }
        else {
                base_this=to->do_fft();
                base_this->process_inplace("xform.phaseorigin.tocorner");
                EMData *tmp = base_this->process("math.squared");
                base_thissq = tmp->do_fft();
                delete tmp;
        }
 
        float sigmathis = params.set_default("sigmathis",0.01f);
        float sigmato = params.set_default("sigmato",0.01f);
        int verbose = params.set_default("verbose",0);
        float maxres = params.set_default("maxres",-1.0f);
 
        if (base_this->get_xsize()!=base_this->get_ysize()+2 || base_this->get_ysize()!=base_this->get_zsize()
                || base_to->get_xsize()!=base_to->get_ysize()+2 || base_to->get_ysize()!=base_to->get_zsize()) throw InvalidCallException("ERROR (RT3DTreeAligner): requires cubic images with even numbered box sizes");
 
//      base_this->process_inplace("mask.wedgefill", Dict("thresh_sigma", sigmathis));
//      base_to->process_inplace("mask.wedgefill", Dict("thresh_sigma", sigmato));
        
        
        base_this->process_inplace("xform.fourierorigin.tocenter");             // easier to chop out Fourier subvolumes
        base_to->process_inplace("xform.fourierorigin.tocenter");
        base_thissq->process_inplace("xform.fourierorigin.tocenter");
        base_mask->process_inplace("xform.fourierorigin.tocenter");
 
        float apix=(float)this_img->get_attr("apix_x");
        int ny=this_img->get_ysize();
        params["boxsize"]=ny;
        int maxshift00=(int)params.set_default("maxshift",ny/4);
 
        int maxny=ny;
        if (0)//maxres>0)
                maxny=4*int(ny*apix/maxres/2+1);
                if (verbose>0)
                        printf("\n\n*******\nmax resolution %1.2f, box size %d\n", maxres, maxny);
        
        float maxang=params.set_default("maxang",-1.0);
        Transform initxf;
        
//      int downsample=floor(ny/20);            // Minimum shrunken box size is 20^3
 
        vector<float> s_score(nsoln,0.0f);
        vector<float> s_coverage(nsoln,0.0f);
        vector<float> s_step(nsoln*3,7.5f);
        vector<Transform> s_xform(nsoln);
        if (verbose>0) printf("%d solutions\n",nsoln);
        
        
        int curiter=-1;
        int sexp_start=4;
        if (params.has_key("initxform")){
                const vector< Transform > xfs=params["initxform"];
                nsoln=xfs.size();
                initxf.set_params(xfs[0].get_params("eman"));
                for (unsigned int i=0; i<nsoln; i++){
                        s_xform[i].set_params(xfs[i].get_params("eman"));
                }
                sexp_start=5;
                curiter=0;
//              for (int i=0; i<nsoln*3; i++) {
//                      s_step[i]=.5;
//              }
        }
 
 
//      float dstep[3] = {7.5,7.5,7.5};         // we take  steps for each of the 3 angles, may be positive or negative
        string axname[] = {"az","alt","phi"};
        int lastss=24;
        // We start with 32^3, 64^3 ...
        for (int sexp=sexp_start; sexp<10; sexp++) {
                curiter++;
//      for (int sexp=4; sexp<5; sexp++) {
                int ss=pow(2.0,sexp);
                if (ss>maxny) ss=maxny;
                if (ss<24) ss=24;               // 16 may be too small, but 32 takes too long...
                else if (ss<48) ss=48;          // 16 may be too small, but 32 takes too long...
                if (verbose>0) printf("\nSize %d\n",ss);
 
                int maxshift=maxshift00*ss/ny;
                if (maxshift00<0) maxshift=-1;
                
                //ss=good_size(ny/ds);
                EMData *small_this=base_this->get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
                EMData *small_to=  base_to->  get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
                EMData *small_thissq=base_thissq->get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
                EMData *small_mask=  base_mask->  get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
                small_this->process_inplace("xform.fourierorigin.tocorner");
                small_this->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
                small_to->process_inplace("xform.fourierorigin.tocorner");
                small_to->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
                small_mask->process_inplace("xform.fourierorigin.tocorner");
                small_mask->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
                small_thissq->process_inplace("xform.fourierorigin.tocorner");
                small_thissq->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
                
                if (ss<maxny){
                        // skip lp filter at full sampling seems to help..
                        small_this->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.33f));
                        small_to->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.33f));
                }
 
                // these are cached for speed in the comparator
                vector<float>sigmathisv=small_this->calc_radial_dist(ss/2,0,1,4);
                vector<float>sigmatov=small_to->calc_radial_dist(ss/2,0,1,4);
                for (int i=0; i<ss/2; i++) {
                        sigmathisv[i]*=sigmathisv[i]*sigmathis;
                        sigmatov[i]*=sigmatov[i]*sigmato;
                }
 
                // This is a solid estimate for very complete searching, 2.5 is a bit arbitrary
                // make sure the altitude step hits 90 degrees, not absolutely necessary for this, but can't hurt
//              float astep = 89.999/floor(pi/(2.5*2.0*atan(2.0/ss)));
                float astep = 89.999/floor(pi/(1.5*2.0*atan(2.0/ss)));
 
                // This is drawn from single particle analysis testing, which in that case insures that enough sampling to
                // reasonably fill Fourier space is achieved, but doesn't perfectly apply to SPT
//              float astep = (float)(89.99/ceil(90.0*9.0/(8.0*sqrt((float)(4300.0/ss)))));     // 8 is (3+speed) from SPA with speed=5
 
                // This insures we make at least one real effort at each level
                for (int i=0; i<s_score.size(); i++) {
                        s_score[i]=1.0e24;      // reset the scores since the different scales will not match
                        if (fabs(s_step[i*3+0])<astep/4.0) s_step[i*3+0]*=2.0;
                        if (fabs(s_step[i*3+1])<astep/4.0) s_step[i*3+1]*=2.0;
                        if (fabs(s_step[i*3+2])<astep/4.0) s_step[i*3+2]*=2.0;
                }
 
                // This is for the first loop, we do a full search in a heavily downsampled space
                if (curiter==0){ //s_coverage[0]==0.0f) {
                        // Genrate points on a sphere in an asymmetric unit
                        if (verbose>1) printf("stage 1 - ang step %1.2f\n",astep);
                        Dict d;
                        d["inc_mirror"] = true;
                        d["delta"] = astep;
                        Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
                        // We don't generate for phi, since this can produce a very large number of orientations
                        vector<Transform> transforms = sym->gen_orientations((string)params.set_default("orientgen","eman"),d);
                        if (verbose>0) printf("%d orientations to test (%lu)\n",(int)(transforms.size()*(360.0/astep)),transforms.size());
                        if (transforms.size()<30) continue; // for very high symmetries we will go up to 32 instead of 24
 
                        // We iterate over all orientations in an asym triangle (alt & az) then deal with phi ourselves
//                      for (std::vector<Transform>::iterator t = transforms.begin(); t!=transforms.end(); ++t) {    // iterator form was causing all sorts of problems
                        for (unsigned int it=0; it<transforms.size(); it++) {
                                
                                if (verbose>2) {
                                        printf("  %d/%lu \r",it,transforms.size());
                                        fflush(stdout);
                                }
                                for (float phi=0; phi<360.0; phi+=astep) {
                                        Transform t = transforms[it];
                                        Dict aap=t.get_params("eman");
                                        aap["phi"]=phi;
                                        aap["tx"]=0;
                                        aap["ty"]=0;
                                        aap["tz"]=0;
                                        t.set_params(aap);
                                        t.invert();
                                        aap=t.get_params("eman");
 
                                        // somewhat strangely, rotations are actually much more expensive than FFTs, so we use a CCF for translation
                                        EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
                                        EMData *sttsq=small_thissq->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
//                                      EMData *ccf=small_to->calc_ccf(stt);
                                        EMData *ccf=small_to->calc_ccf_masked(stt,sttsq,small_mask);
                                        IntPoint ml=ccf->calc_max_location_wrap();
 
                                        aap["tx"]=(int)ml[0];
                                        aap["ty"]=(int)ml[1];
                                        aap["tz"]=(int)ml[2];
                                        t.set_params(aap);
                                        delete stt;
                                        delete sttsq;
                                        delete ccf;
                                        stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));        // we have to do 1 slow transform here now that we have the translation
 
//                                      float sim=stt->cmp("ccc.tomo.thresh",small_to,Dict("sigmaimg",sigmathis,"sigmawith",sigmato));
                                        float sim=stt->cmp("ccc.tomo.thresh",small_to);
                                        
//                                      float sim=stt->cmp("fsc.tomo.auto",small_to,Dict("sigmaimg",sigmathisv,"sigmawith",sigmatov));
//                                      float sim=stt->cmp("fsc.tomo.auto",small_to);
 
                                        // We want to make sure our starting points are somewhat separated from each other, so we replace any angles too close to an existing angle
                                        // If we find an existing 'best' angle within range, then we either replace it or skip
                                        int worst=-1;
                                        float worstv=1.0e20;
                                        for (int i=0; i<nsoln; i++) {
                                                if (s_coverage[i]==0.0) continue;       // hasn't been set yet
                                                Transform tdif=s_xform[i].inverse();
                                                tdif=tdif*t;
                                                float adif=tdif.get_rotation("spin")["omega"];
                                                if (adif<astep*2.5) {
                                                        worst=i;
//                                                      printf("= %1.3f\n",adif);
                                                }
                                        }
 
                                        // if we weren't close to an existing angle, then we find the lowest current score and use that
                                        if (worst==-1) {
                                                // First we find the worst solution in the list of possible best solutions, or the first
                                                // solution which is currently "empty"
                                                for (int i=0; i<nsoln; i++) {
                                                        if (s_coverage[i]==0.0) { worst=i; break; }
                                                        if (s_score[i]<worstv) {worst=i; worstv=s_score[i];}
                                                }
                                        }
 
                                        // If the current solution is better than the 'worst' of the previous solutions, then we
                                        // displace it. Note that there is no sorting performed here
                                        if (sim<s_score[worst]) {
                                                s_score[worst]=sim;
                                                s_coverage[worst]=stt->get_attr("fft_overlap");
                                                s_xform[worst]=t;
                                                //printf("%f\t%f\t%d\n",s_score[worst],s_coverage[worst],worst);
                                        }
                                        delete stt;
                                }
                        }
                        if (verbose>2) printf("\n");
//                      for (int i=0; i<nsoln; i++) {
//                              Dict d=s_xform[i].get_params("eman");
//                              printf("%d, %f, %f, %f, %f\n",i,(float)d["alt"], (float)d["az"], (float)d["phi"], s_score[i]);
//                      }
 
                }
                // Once we have our initial list of best locations, we just refine each possibility individually
                else {
                        // We generate a search pattern around each existing solution
                        float res=float(ny)/float(ss)*apix*2;
                        if (verbose>1) printf("stage 2 - maxres %1.2f, astep %1.2f\n",res, astep);
                        if (ss>=maxny && params.has_key("initxform")){
                                // last iteration, add back the initial transform so we do not do worse than the last run
                                s_xform[nsoln]=initxf;
                                nsoln+=1;
                        }
                        
                        for (int i=0; i<nsoln; i++) {
 
                                if (verbose>2) {
                                        printf("  %d\t%d\r",i,nsoln);
                                        fflush(stdout);
                                }
                                // We work an axis at a time until we get where we want to be. Somewhat like a simplex
                                int changed=1;
                                Dict upd;
                                testort(small_this,small_to,small_mask,small_thissq,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift);
                                while (changed) {
                                        changed=0;
                                        for (int axis=0; axis<3; axis++) {
                                                if (fabs(s_step[i*3+axis])<astep/4.0) continue;         // skip axes where we already have enough precision on this axis
                                                upd[axname[axis]]=s_step[i*3+axis];
                                                // when moving az, we move phi in the opposite direction by the same amount since the two are singular at alt=0
                                                // phi continues to move independently. I believe this should produce a more monotonic energy surface
                                                if (axis==0) upd[axname[2]]=-s_step[i*3+axis];
 
                                                int r=testort(small_this,small_to,small_mask,small_thissq,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift);
 
                                                // If we fail, we reverse direction with a slightly smaller step and try that
                                                // Whether this fails or not, we move on to the next axis
                                                if (r) changed=1;
                                                else {
                                                        s_step[i*3+axis]*=-0.75;
                                                        upd[axname[axis]]=s_step[i*3+axis];
                                                        r=testort(small_this,small_to,small_mask,small_thissq,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift);
                                                        if (r) changed=1;
                                                }
                                                if (verbose>4) printf("\nX %1.3f\t%1.3f\t%1.3f\t%d\t",s_step[i*3],s_step[i*3+1],s_step[i*3+2],changed);
                                        }
                                        if (verbose>3) {
                                                        Dict aap=s_xform[i].get_params("eman");
                                                        printf("\n%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t(%1.3f)",s_step[i*3],s_step[i*3+1],s_step[i*3+2],float(aap["az"]),float(aap["alt"]),float(aap["phi"]),s_score[i]);
                                        }
 
                                        if (!changed) {
//                                              for (int j=0; j<3; j++) s_step[i*3+j]*-0.75;
                                                changed=1;
                                        }
                                        if (fabs(s_step[i*3])<astep/4 && fabs(s_step[i*3+1])<astep/4 && fabs(s_step[i*3+2])<astep/4) changed=0;
                                }
 
                                // Ouch, exhaustive (local) search
//                              for (int daz=-1; daz<=1; daz++) {
//                                      for (int dalt=-1; dalt<=1; dalt++) {
//                                              for (int dphi=-1; dphi<=1; dphi++) {
//                                                      Dict upd;
//                                                      upd["az"]=daz*astep;
//                                                      upd["alt"]=dalt*astep;
//                                                      upd["phi"]=dphi*astep;
//                                                      int r=testort(small_this,small_to,s_score,s_coverage,s_xform,i,upd);
//                                              }
//                                      }
//                              }
                        }
                }
                // lazy earlier in defining s_ vectors, so lazy here too and inefficiently sorting
                // We are sorting inside the outermost loop so we can decrease the number of solutions
                // before we get to the finest precision
                for (unsigned int i=0; i<nsoln-1; i++) {
                        for (unsigned int j=i+1; j<nsoln; j++) {
                                if (s_score[i]>s_score[j]) {
                                        float t=s_score[i]; s_score[i]=s_score[j]; s_score[j]=t;
                                        t=s_coverage[i]; s_coverage[i]=s_coverage[j]; s_coverage[j]=t;
                                        Transform tt=s_xform[i]; s_xform[i]=s_xform[j]; s_xform[j]=tt;
                                }
                        }
                }
 
                // At each level of sampling we (potentially) decrease the number of answers we check in detail
                // assuming we are gradually homing in on the best solution
                nsoln/=2;
                if (nsoln<nrsoln) nsoln=nrsoln;
 
 
                delete small_this;
                delete small_to;
                delete small_thissq;
                delete small_mask;
                
                lastss=ss;
                if (ss>=maxny && curiter>0) break;
        }
 
        delete base_this;
        delete base_to;
        delete base_thissq;
        delete base_mask;
        
        // note the translations are wrong when the alignment stops before ss==ny
        if (lastss<ny){
                for (unsigned int i = 0; i < nrsoln; ++i ) {
                        Transform tt=s_xform[i];
                        Vec3f v=tt.get_trans();
                        v=v*ny/lastss;
                        tt.set_trans(v);
                        s_xform[i]=tt;
                        
                }
        }
 
        // initialize results
        vector<Dict> solns;
        for (unsigned int i = 0; i < nrsoln; ++i ) {
                Dict d;
                d["score"] = s_score[i];
                d["coverage"] = s_coverage[i];
                s_xform[i].invert();    // this is because we inverted the order of the input images above to match convention
                d["xform.align3d"] = &s_xform[i];
                solns.push_back(d);
        }
 
        return solns;
}

References EMAN::EMData::calc_ccf_masked(), EMAN::EMData::calc_radial_dist(), EMAN::Symmetry3D::gen_orientations(), EMAN::Factory< T >::get(), EMAN::EMData::get_clip(), EMAN::Aligner::get_params(), EMAN::Transform::get_params(), EMAN::Transform::get_rotation(), EMAN::Transform::get_trans(), EMAN::Dict::has_key(), InvalidCallException, InvalidParameterException, EMAN::Transform::inverse(), EMAN::Transform::invert(), EMAN::Aligner::params, EMAN::Dict::set_default(), EMAN::Transform::set_params(), EMAN::Transform::set_trans(), EMAN::Dict::size(), and testort().

Member Data Documentation

NAME

const string RT3DLocalTreeAligner::NAME = "rotate_translate_3d_local_tree"

static

Definition at line 2012 of file aligner.h.

Referenced by get_name().

---

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

• aligner.h
• aligner.cpp