EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Member Functions | List of all members
EMAN::RT3DTreeAligner 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::RT3DTreeAligner:
Inheritance graph
[legend]
Collaboration diagram for EMAN::RT3DTreeAligner:
Collaboration graph
[legend]

Public Member Functions

virtual EMDataalign (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 EMDataalign (EMData *this_img, EMData *to_img) const
 See Aligner comments for more details. More...
 
virtual vector< Dictxform_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 AlignerNEW ()
 

Static Public Attributes

static const string NAME = "rotate_translate_3d_tree"
 

Private Member Functions

bool testort (EMData *small_this, EMData *small_to, 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, EMData *mask) 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.

Parameters
symThe symmtery to use as the basis of the spherical sampling
verboseTurn this on to have useful information printed to standard out
Author
Steve Ludtke
Date
April 2015

Definition at line 1889 of file aligner.h.

Member Function Documentation

◆ align() [1/2]

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

See Aligner comments for more details.

Implements EMAN::Aligner.

Definition at line 1898 of file aligner.h.

1899 {
1900 return align(this_img, to_img, "sqeuclidean", Dict());
1901 }
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.

References align().

◆ align() [2/2]

virtual EMData * EMAN::RT3DTreeAligner::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::RT3DTreeAligner::get_desc ( ) const
inlinevirtual

Implements EMAN::Aligner.

Definition at line 1913 of file aligner.h.

1914 {
1915 return "3D rotational and translational alignment using a hierarchical approach in Fourier space. Should be very fast and not require 'refine' alignment. 'this' should be the fixed reference, aligned to symmetry axes and mask if provided. If masking, provide the mask rather than premasking the volume.";
1916 }

◆ get_name()

virtual string EMAN::RT3DTreeAligner::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 1908 of file aligner.h.

1909 {
1910 return NAME;
1911 }
static const string NAME
Definition: aligner.h:1943

References NAME.

◆ get_param_types()

virtual TypeDict EMAN::RT3DTreeAligner::get_param_types ( ) const
inlinevirtual

Implements EMAN::Aligner.

Definition at line 1923 of file aligner.h.

1924 {
1925 TypeDict d;
1926 d.put("mask", EMObject::EMDATA,"A mask under which to do the alignment. Mask relative to 'this'");
1927 d.put("maxshift", EMObject::INT,"maximum shift allowed");
1928 d.put("sym", EMObject::STRING,"The symmtery to use as the basis of the spherical sampling. Default is c1 (no symmetry)");
1929 d.put("sigmathis", EMObject::FLOAT,"Only Fourier voxels larger than sigma times this value will be considered");
1930 d.put("sigmato", EMObject::FLOAT,"Only Fourier voxels larger than sigma times this value will be considered");
1931 d.put("maxres", EMObject::FLOAT,"maximum resolution to compare");
1932 d.put("minres", EMObject::FLOAT,"minimum resolution to compare");
1933 d.put("maxang", EMObject::FLOAT,"maximum angle from initial rotation.");
1934 d.put("initxform", EMObject::TRANSFORMARRAY,"An array of Transforms storing the starting positions.");
1935
1936// d.put("initxform", EMObject::TRANSFORM,"The Transform storing the starting position. If unspecified the identity matrix is used");
1937 d.put("randphi", EMObject::BOOL,"Ignore phi constraint for refine search");
1938 d.put("rand180", EMObject::BOOL,"Ignore 180 rotation for refine search");
1939 d.put("verbose", EMObject::BOOL,"Turn this on to have useful information printed to standard out.");
1940 return d;
1941 }

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

◆ NEW()

static Aligner * EMAN::RT3DTreeAligner::NEW ( )
inlinestatic

Definition at line 1918 of file aligner.h.

1919 {
1920 return new RT3DTreeAligner();
1921 }

◆ testort()

bool RT3DTreeAligner::testort ( EMData small_this,
EMData small_to,
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,
EMData mask 
) const
private

Definition at line 3809 of file aligner.cpp.

3809 {
3810 Transform t;
3811 Dict aap=s_xform[i].get_params("eman");
3812 aap["tx"]=0;
3813 aap["ty"]=0;
3814 aap["tz"]=0;
3815 for (Dict::const_iterator p=upd.begin(); p!=upd.end(); p++) {
3816 aap[p->first]=(float)aap[p->first]+(float)p->second;
3817 }
3818
3819 t.set_params(aap);
3820 float maxang=params.set_default("maxang",-1.0);
3821 bool randphi=params.set_default("randphi",false);
3822 bool rand180=params.set_default("rand180",false);
3823
3824 if (maxang>0){
3825
3826 Transform tmp=initxf * t.inverse();
3827
3828 if (randphi){
3829 Dict td=t.get_params("eman");
3830 Dict ti=initxf.get_params("eman");
3831 td["phi"]=ti["phi"];
3832 Transform tmp1=Transform(td);
3833 tmp=initxf * tmp1.inverse();
3834 }
3835
3836 float r=tmp.get_params("spin")["omega"];
3837 if (rand180){
3838 r=r>90?(180-r):r;
3839 }
3840 if(r>maxang)
3841 return false;
3842
3843 }
3844
3845 int ny=small_this->get_ysize();
3846 if (params.has_key("initxform")){
3847 // when doing refinement, search around the given position
3848 t.set_trans(initxf.get_trans()*ny/(float)params["boxsize"]);
3849 aap=t.get_params("eman");
3850 }
3851
3852 // rotate in Fourier space then use a CCF to find translation
3853 EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
3854 EMData *ccf=small_to->calc_ccf(stt);
3855// EMData *ccf=small_to->calc_flcf(stt);
3856 IntPoint ml=ccf->calc_max_location_wrap(maxshift, maxshift, maxshift);
3857
3858
3859 aap["tx"]=int(aap["tx"])+(int)ml[0];
3860 aap["ty"]=int(aap["ty"])+(int)ml[1];
3861 aap["tz"]=int(aap["tz"])+(int)ml[2];
3862 t.set_params(aap);
3863 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
3864
3865// float maxres = params.set_default("maxres",0.0f);
3866// float minres = params.set_default("minres",0.0f);
3867
3868 float sim=st2->cmp("fsc.tomo.auto",small_to,Dict("sigmaimg",sigmathisv,"sigmawith",sigmatov, "pmin", 4, "pmax",ny/3));
3869// float sim=st2->cmp("ccc.tomo.thresh",small_to,Dict("sigmaimg",sigmathis,"sigmawith",sigmato));
3870// printf("\nTESTORT %6.1f %6.1f %6.1f\t%4d %4d %4d\t%1.5g\t%1.5g %d (%d)",
3871// 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());
3872
3873 delete ccf;
3874 // If the score is better than before, we update this particular best value
3875 if (sim<s_score[i]) {
3876 s_score[i]=sim;
3877 s_coverage[i]=st2->get_attr("fft_overlap");
3878 s_xform[i]=t;
3879 delete stt;
3880 delete st2;
3881 return true;
3882 }
3883 delete stt;
3884 delete st2;
3885 return false;
3886}
Dict params
Definition: aligner.h:147
Const iterator support for the Dict object This is just a wrapper, everything is inherited from the m...
Definition: emobject.h:674
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
iterator end()
Definition: emobject.cpp:1061
iterator begin()
Definition: emobject.cpp:1045
type set_default(const string &key, type val)
Default setting behavior This can be achieved using a template - d.woolford Jan 2008 (before there wa...
Definition: emobject.h:569
bool has_key(const string &key) const
Ask the Dictionary if it as a particular key.
Definition: emobject.h:511
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
EMObject is a wrapper class for types including int, float, double, etc as defined in ObjectType.
Definition: emobject.h:123
IntPoint defines an integer-coordinate point in a 1D/2D/3D space.
Definition: geometry.h:192
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
Dict get_params(const string &euler_type) const
Get the parameters of the entire transform, using a specific euler convention.
Definition: transform.cpp:479
Vec3f get_trans() const
Get the post trans as a vec3f.
Definition: transform.cpp:1046
void set_params(const Dict &d)
Set the parameters of the entire transform.
Definition: transform.cpp:254
Transform inverse() const
Get the inverse of this transformation matrix.
Definition: transform.cpp:1327
void set_trans(const float &x, const float &y, const float &z=0)
Set the post translation component.
Definition: transform.cpp:1036
EMData * calc_ccf(EMData *with=0, fp_flag fpflag=CIRCULANT, bool center=false)
Calculate Cross-Correlation Function (CCF).
Definition: emdata.cpp:1499

References EMAN::Dict::begin(), EMAN::EMData::calc_ccf(), 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 > RT3DTreeAligner::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 3423 of file aligner.cpp.

3423 {
3424 if (nrsoln == 0) throw InvalidParameterException("ERROR (RT3DTreeAligner): nsoln must be >0"); // What was the user thinking?
3425
3426 int nsoln = nrsoln*2;
3427 if (nrsoln<32) nsoln=64; // we start with at least 32 solutions, but then gradually decrease with increasing scale
3428
3429 float sigmathis = params.set_default("sigmathis",0.01f);
3430 float sigmato = params.set_default("sigmato",0.01f);
3431 int verbose = params.set_default("verbose",0);
3432 float maxres = params.set_default("maxres",-1.0f);
3433 EMData *mask = params.set_default("mask",(EMData *)0);
3434
3435 // !!!!!! IMPORTANT NOTE - we are inverting the order of this and to here to match convention in other aligners, to compensate
3436 // the Transform is inverted before being returned
3437 EMData *base_this;
3438 EMData *base_to;
3439 if (mask) {
3440 if (this_img->is_complex()) {
3441 EMData *tmp = this_img->do_ift();
3442 tmp->process_inplace("xform.phaseorigin.tocenter");
3443 tmp->process_inplace("normalize.mask",Dict("mask",mask,"apply_mask",1));
3444 base_to=tmp->do_fft();
3445 base_to->process_inplace("xform.phaseorigin.tocorner");
3446 delete tmp;
3447 }
3448 else {
3449 EMData *tmp = this_img->process("normalize.mask",Dict("mask",mask,"apply_mask",1));
3450 base_to=tmp->do_fft();
3451 base_to->process_inplace("xform.phaseorigin.tocorner");
3452 delete tmp;
3453 }
3454
3455 if (to->is_complex()) {
3456 base_this=to->copy();
3457 }
3458 else {
3459 base_this=to->do_fft();
3460 base_this->process_inplace("xform.phaseorigin.tocorner");
3461 }
3462 }
3463 else {
3464 if (this_img->is_complex()) base_to=this_img->copy();
3465 else {
3466 base_to=this_img->do_fft();
3467 base_to->process_inplace("xform.phaseorigin.tocorner");
3468 }
3469
3470 if (to->is_complex()) base_this=to->copy();
3471 else {
3472 base_this=to->do_fft();
3473 base_this->process_inplace("xform.phaseorigin.tocorner");
3474 }
3475 }
3476
3477
3478 if (base_this->get_xsize()!=base_this->get_ysize()+2 || base_this->get_ysize()!=base_this->get_zsize()
3479 || 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");
3480
3481// base_this->process_inplace("mask.wedgefill", Dict("thresh_sigma", sigmathis));
3482// base_to->process_inplace("mask.wedgefill", Dict("thresh_sigma", sigmato));
3483
3484
3485 base_this->process_inplace("xform.fourierorigin.tocenter"); // easier to chop out Fourier subvolumes
3486 base_to->process_inplace("xform.fourierorigin.tocenter");
3487
3488 float apix=(float)this_img->get_attr("apix_x");
3489 int ny=this_img->get_ysize();
3490 params["boxsize"]=ny;
3491
3492 int maxny=ny;
3493 if (0)//maxres>0)
3494 maxny=4*int(ny*apix/maxres/2+1);
3495 if (verbose>0)
3496 printf("\n\n*******\nmax resolution %1.2f, box size %d\n", maxres, maxny);
3497
3498
3499 Transform initxf;
3500
3501// int downsample=floor(ny/20); // Minimum shrunken box size is 20^3
3502
3503 vector<float> s_score(nsoln,0.0f);
3504 vector<float> s_coverage(nsoln,0.0f);
3505 vector<float> s_step(nsoln*3,7.5f);
3506 vector<Transform> s_xform(nsoln);
3507 if (verbose>0) printf("%d solutions\n",nsoln);
3508
3509
3510 int curiter=-1;
3511 int sexp_start=4;
3512 if (params.has_key("initxform")){
3513 const vector< Transform > xfs=params["initxform"];
3514 nsoln=xfs.size();
3515 initxf.set_params(xfs[0].get_params("eman"));
3516 for (unsigned int i=0; i<nsoln; i++){
3517 s_xform[i].set_params(xfs[i].get_params("eman"));
3518 }
3519 sexp_start=5;
3520 curiter=0; // bypass the first round of global orientation search
3521// for (int i=0; i<nsoln*3; i++) {
3522// s_step[i]=.5;
3523// }
3524 }
3525 else if (params.has_key("maxang")) { // if we got maxshift or maxang without initxform, we need to set up identity xform
3526 nsoln=1; // Default transform should be identity, so just set number to 1
3527 curiter=0; // skip global search
3528 sexp_start=5;
3529 }
3530
3531 int maxshift00=(int)params.set_default("maxshift",ny/4);
3532 float maxang=params.set_default("maxang",-1.0);
3533
3534// float dstep[3] = {7.5,7.5,7.5}; // we take steps for each of the 3 angles, may be positive or negative
3535 string axname[] = {"az","alt","phi"};
3536 int lastss=24;
3537 // We start with 32^3, 64^3 ...
3538 for (int sexp=sexp_start; sexp<10; sexp++) {
3539 curiter++;
3540// for (int sexp=4; sexp<5; sexp++) {
3541 int ss=pow(2.0,sexp);
3542 if (ss>maxny) ss=maxny;
3543 if (ss<24) ss=24; // 16 may be too small, but 32 takes too long...
3544 else if (ss<48) ss=48; // 16 may be too small, but 32 takes too long...
3545 if (verbose>0) printf("\nSize %d\n",ss);
3546
3547 int maxshift=maxshift00*ss/ny;
3548 if (maxshift00<0) maxshift=-1;
3549
3550 //ss=good_size(ny/ds);
3551 EMData *small_this=base_this->get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
3552 EMData *small_to= base_to-> get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
3553 small_this->process_inplace("xform.fourierorigin.tocorner");
3554 small_this->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
3555 small_to->process_inplace("xform.fourierorigin.tocorner");
3556 small_to->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
3557
3558 if (ss<maxny){
3559 // skip lp filter at full sampling seems to help..
3560 small_this->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.33f));
3561 small_to->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.33f));
3562 }
3563
3564 // these are cached for speed in the comparator
3565 vector<float>sigmathisv=small_this->calc_radial_dist(ss/2,0,1,4);
3566 vector<float>sigmatov=small_to->calc_radial_dist(ss/2,0,1,4);
3567 for (int i=0; i<ss/2; i++) {
3568 sigmathisv[i]*=sigmathisv[i]*sigmathis;
3569 sigmatov[i]*=sigmatov[i]*sigmato;
3570 }
3571
3572 // This is a solid estimate for very complete searching, 2.5 is a bit arbitrary
3573 // make sure the altitude step hits 90 degrees, not absolutely necessary for this, but can't hurt
3574// float astep = 89.999/floor(pi/(2.5*2.0*atan(2.0/ss)));
3575 float astep = 89.999/floor(pi/(1.5*2.0*atan(2.0/ss)));
3576
3577 // This is drawn from single particle analysis testing, which in that case insures that enough sampling to
3578 // reasonably fill Fourier space is achieved, but doesn't perfectly apply to SPT
3579// 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
3580
3581 // This insures we make at least one real effort at each level
3582 for (int i=0; i<s_score.size(); i++) {
3583 s_score[i]=1.0e24; // reset the scores since the different scales will not match
3584 if (fabs(s_step[i*3+0])<astep/4.0) s_step[i*3+0]*=2.0;
3585 if (fabs(s_step[i*3+1])<astep/4.0) s_step[i*3+1]*=2.0;
3586 if (fabs(s_step[i*3+2])<astep/4.0) s_step[i*3+2]*=2.0;
3587 }
3588
3589 // This is for the first loop, we do a full search in a heavily downsampled space
3590 if (curiter==0){ //s_coverage[0]==0.0f) {
3591 // Genrate points on a sphere in an asymmetric unit
3592 if (verbose>1) printf("stage 1 - ang step %1.2f\n",astep);
3593 Dict d;
3594 d["inc_mirror"] = true;
3595 d["delta"] = astep;
3596 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
3597 // We don't generate for phi, since this can produce a very large number of orientations
3598 vector<Transform> transforms = sym->gen_orientations((string)params.set_default("orientgen","eman"),d);
3599 if (verbose>0) printf("%d orientations to test (%lu)\n",(int)(transforms.size()*(360.0/astep)),transforms.size());
3600 if (transforms.size()<30) continue; // for very high symmetries we will go up to 32 instead of 24
3601
3602 // We iterate over all orientations in an asym triangle (alt & az) then deal with phi ourselves
3603// for (std::vector<Transform>::iterator t = transforms.begin(); t!=transforms.end(); ++t) { // iterator form was causing all sorts of problems
3604 for (unsigned int it=0; it<transforms.size(); it++) {
3605
3606 if (verbose>2) {
3607 printf(" %d/%lu \r",it,transforms.size());
3608 fflush(stdout);
3609 }
3610 for (float phi=0; phi<360.0; phi+=astep) {
3611 Transform t = transforms[it];
3612 Dict aap=t.get_params("eman");
3613 aap["phi"]=phi;
3614 aap["tx"]=0;
3615 aap["ty"]=0;
3616 aap["tz"]=0;
3617 t.set_params(aap);
3618 t.invert();
3619 aap=t.get_params("eman");
3620
3621 // somewhat strangely, rotations are actually much more expensive than FFTs, so we use a CCF for translation
3622 EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
3623 EMData *ccf=small_to->calc_ccf(stt);
3624 IntPoint ml=ccf->calc_max_location_wrap();
3625
3626 aap["tx"]=(int)ml[0];
3627 aap["ty"]=(int)ml[1];
3628 aap["tz"]=(int)ml[2];
3629 t.set_params(aap);
3630 delete stt;
3631 delete ccf;
3632 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
3633
3634// float sim=stt->cmp("ccc.tomo.thresh",small_to,Dict("sigmaimg",sigmathis,"sigmawith",sigmato));
3635 float sim=stt->cmp("ccc.tomo.thresh",small_to);
3636
3637// float sim=stt->cmp("fsc.tomo.auto",small_to,Dict("sigmaimg",sigmathisv,"sigmawith",sigmatov));
3638// float sim=stt->cmp("fsc.tomo.auto",small_to);
3639
3640 // 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
3641 // If we find an existing 'best' angle within range, then we either replace it or skip
3642 int worst=-1;
3643 float worstv=1.0e20;
3644 for (int i=0; i<nsoln; i++) {
3645 if (s_coverage[i]==0.0) continue; // hasn't been set yet
3646 Transform tdif=s_xform[i].inverse();
3647 tdif=tdif*t;
3648 float adif=tdif.get_rotation("spin")["omega"];
3649 if (adif<astep*2.5) {
3650 worst=i;
3651// printf("= %1.3f\n",adif);
3652 }
3653 }
3654
3655 // if we weren't close to an existing angle, then we find the lowest current score and use that
3656 if (worst==-1) {
3657 // First we find the worst solution in the list of possible best solutions, or the first
3658 // solution which is currently "empty"
3659 for (int i=0; i<nsoln; i++) {
3660 if (s_coverage[i]==0.0) { worst=i; break; }
3661 if (s_score[i]<worstv) {worst=i; worstv=s_score[i];}
3662 }
3663 }
3664
3665 // If the current solution is better than the 'worst' of the previous solutions, then we
3666 // displace it. Note that there is no sorting performed here
3667 if (sim<s_score[worst]) {
3668 s_score[worst]=sim;
3669 s_coverage[worst]=stt->get_attr("fft_overlap");
3670 s_xform[worst]=t;
3671 //printf("%f\t%f\t%d\n",s_score[worst],s_coverage[worst],worst);
3672 }
3673 delete stt;
3674 }
3675 }
3676 if (verbose>2) printf("\n");
3677// for (int i=0; i<nsoln; i++) {
3678// Dict d=s_xform[i].get_params("eman");
3679// printf("%d, %f, %f, %f, %f\n",i,(float)d["alt"], (float)d["az"], (float)d["phi"], s_score[i]);
3680// }
3681
3682 }
3683 // Once we have our initial list of best locations, we just refine each possibility individually
3684 else {
3685 // We generate a search pattern around each existing solution
3686 float res=float(ny)/float(ss)*apix*2;
3687 if (verbose>1) printf("stage 2 - maxres %1.2f, astep %1.2f\n",res, astep);
3688 if (ss>=maxny && params.has_key("initxform")){
3689 // last iteration, add back the initial transform so we do not do worse than the last run
3690 s_xform[nsoln]=initxf;
3691 nsoln+=1;
3692 }
3693
3694 for (int i=0; i<nsoln; i++) {
3695
3696 if (verbose>2) {
3697 printf(" %d\t%d\r",i,nsoln);
3698 fflush(stdout);
3699 }
3700 // We work an axis at a time until we get where we want to be. Somewhat like a simplex
3701 int changed=1;
3702 Dict upd;
3703 testort(small_this,small_to,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift,mask);
3704 while (changed) {
3705 changed=0;
3706 for (int axis=0; axis<3; axis++) {
3707 if (fabs(s_step[i*3+axis])<astep/4.0) continue; // skip axes where we already have enough precision on this axis
3708 upd[axname[axis]]=s_step[i*3+axis];
3709 // when moving az, we move phi in the opposite direction by the same amount since the two are singular at alt=0
3710 // phi continues to move independently. I believe this should produce a more monotonic energy surface
3711 if (axis==0) upd[axname[2]]=-s_step[i*3+axis];
3712
3713 int r=testort(small_this,small_to,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift,mask);
3714
3715 // If we fail, we reverse direction with a slightly smaller step and try that
3716 // Whether this fails or not, we move on to the next axis
3717 if (r) changed=1;
3718 else {
3719 s_step[i*3+axis]*=-0.75;
3720 upd[axname[axis]]=s_step[i*3+axis];
3721 r=testort(small_this,small_to,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift,mask);
3722 if (r) changed=1;
3723 }
3724 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);
3725 }
3726 if (verbose>3) {
3727 Dict aap=s_xform[i].get_params("eman");
3728 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]);
3729 }
3730
3731 if (!changed) {
3732// for (int j=0; j<3; j++) s_step[i*3+j]*-0.75;
3733 changed=1;
3734 }
3735 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;
3736 }
3737
3738 // Ouch, exhaustive (local) search
3739// for (int daz=-1; daz<=1; daz++) {
3740// for (int dalt=-1; dalt<=1; dalt++) {
3741// for (int dphi=-1; dphi<=1; dphi++) {
3742// Dict upd;
3743// upd["az"]=daz*astep;
3744// upd["alt"]=dalt*astep;
3745// upd["phi"]=dphi*astep;
3746// int r=testort(small_this,small_to,s_score,s_coverage,s_xform,i,upd);
3747// }
3748// }
3749// }
3750 }
3751 }
3752 // lazy earlier in defining s_ vectors, so lazy here too and inefficiently sorting
3753 // We are sorting inside the outermost loop so we can decrease the number of solutions
3754 // before we get to the finest precision
3755 for (unsigned int i=0; i<nsoln-1; i++) {
3756 for (unsigned int j=i+1; j<nsoln; j++) {
3757 if (s_score[i]>s_score[j]) {
3758 float t=s_score[i]; s_score[i]=s_score[j]; s_score[j]=t;
3759 t=s_coverage[i]; s_coverage[i]=s_coverage[j]; s_coverage[j]=t;
3760 Transform tt=s_xform[i]; s_xform[i]=s_xform[j]; s_xform[j]=tt;
3761 }
3762 }
3763 }
3764
3765 // At each level of sampling we (potentially) decrease the number of answers we check in detail
3766 // assuming we are gradually homing in on the best solution
3767 nsoln/=2;
3768 if (nsoln<nrsoln) nsoln=nrsoln;
3769
3770
3771 delete small_this;
3772 delete small_to;
3773
3774 lastss=ss;
3775 if (ss>=maxny && curiter>0) break;
3776 }
3777
3778 delete base_this;
3779 delete base_to;
3780
3781 // note the translations are wrong when the alignment stops before ss==ny
3782 if (lastss<ny){
3783 for (unsigned int i = 0; i < nrsoln; ++i ) {
3784 Transform tt=s_xform[i];
3785 Vec3f v=tt.get_trans();
3786 v=v*ny/lastss;
3787 tt.set_trans(v);
3788 s_xform[i]=tt;
3789
3790 }
3791 }
3792
3793 // initialize results
3794 vector<Dict> solns;
3795 for (unsigned int i = 0; i < nrsoln; ++i ) {
3796 Dict d;
3797 d["score"] = s_score[i];
3798 d["coverage"] = s_coverage[i];
3799 s_xform[i].invert(); // this is because we inverted the order of the input images above to match convention
3800 d["xform.align3d"] = &s_xform[i];
3801 solns.push_back(d);
3802 }
3803
3804 return solns;
3805}
virtual Dict get_params() const
Get the Aligner parameters in a key/value dictionary.
Definition: aligner.h:112
size_t size() const
Ask the Dictionary for its size.
Definition: emobject.h:519
EMData * get_clip(const Region &area, const float fill=0) const
Get an inclusive clip.
Definition: emdata.cpp:592
vector< float > calc_radial_dist(int n, float x0, float dx, int inten)
calculates radial distribution.
Definition: emdata.cpp:2781
static T * get(const string &instance_name)
Definition: emobject.h:781
bool testort(EMData *small_this, EMData *small_to, 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, EMData *mask) const
Definition: aligner.cpp:3809
Region defines a 2D or 3D rectangular region specified by its origin coordinates and all edges' sizes...
Definition: geometry.h:497
Symmetry3D - A base class for 3D Symmetry objects.
Definition: symmetry.h:57
vector< Transform > gen_orientations(const string &generatorname="eman", const Dict &parms=Dict())
Ask the Symmetry3D object to generate a set of orientations in its asymmetric unit using an Orientati...
Definition: symmetry.cpp:167
Dict get_rotation(const string &euler_type="eman") const
Get a rotation in any Euler format.
Definition: transform.cpp:829
void invert()
Get the inverse of this transformation matrix.
Definition: transform.cpp:1293
#define InvalidParameterException(desc)
Definition: exception.h:361
#define InvalidCallException(desc)
Definition: exception.h:348

References EMAN::EMData::calc_ccf(), 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 RT3DTreeAligner::NAME = "rotate_translate_3d_tree"
static

Definition at line 1943 of file aligner.h.

Referenced by get_name().


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