EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Member Functions | List of all members
EMAN::RT2DTreeAligner Class Reference

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

#include <aligner.h>

Inheritance diagram for EMAN::RT2DTreeAligner:
Inheritance graph
[legend]
Collaboration diagram for EMAN::RT2DTreeAligner:
Collaboration graph
[legend]

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_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) const
 

Additional Inherited Members

- Protected Attributes inherited from EMAN::Aligner
Dict params
 

Detailed Description

2D 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
verboseTurn this on to have useful information printed to standard out
Author
Steve Ludtke
Date
July 2016

Definition at line 1762 of file aligner.h.

Member Function Documentation

◆ align() [1/2]

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

See Aligner comments for more details.

Implements EMAN::Aligner.

Definition at line 1771 of file aligner.h.

1772 {
1773 return align(this_img, to_img, "sqeuclidean", Dict());
1774 }
EMAN::RT2DTreeAligner::align
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::RT2DTreeAligner::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::RT2DTreeAligner::get_desc ( ) const
inlinevirtual

Implements EMAN::Aligner.

Definition at line 1786 of file aligner.h.

1787 {
1788 return "2D rotational and translational alignment using a hierarchical approach in Fourier space. flip options specifies whether this is RT or RTF. No 'refine' alignment should be required +-1 pixel.";
1789 }

◆ get_name()

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

1782 {
1783 return NAME;
1784 }
EMAN::RT2DTreeAligner::NAME
static const string NAME
Definition: aligner.h:1809

References NAME.

◆ get_param_types()

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

Implements EMAN::Aligner.

Definition at line 1797 of file aligner.h.

1798 {
1799 TypeDict d;
1800// d.put("sigmathis", EMObject::FLOAT,"Only Fourier voxels larger than sigma times this value will be considered");
1801// d.put("initxform", EMObject::TRANSFORM,"The Transform storing the starting position. If unspecified the identity matrix is used");
1802 d.put("flip", EMObject::BOOL,"Include flip in the alignment search. Default True");
1803 d.put("verbose", EMObject::INT,"Turn this on to have useful information printed to standard out.");
1804 d.put("maxshift", EMObject::INT,"Maximum acceptable translation. Used only approximately.");
1805 d.put("maxres", EMObject::FLOAT,"Maximum resolution to consider when full sampling is used");
1806 return d;
1807 }

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

◆ NEW()

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

Definition at line 1791 of file aligner.h.

1792 {
1793 return new RT2DTreeAligner();
1794 }

◆ testort()

bool EMAN::RT2DTreeAligner::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 
) const
private

◆ xform_align_nbest()

vector< Dict > RT2DTreeAligner::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 2731 of file aligner.cpp.

2731 {
2732 if (nrsoln == 0) throw InvalidParameterException("ERROR (RT2DTreeAligner): nsoln must be >0"); // What was the user thinking?
2733
2734 int nsoln = nrsoln*2;
2735 if (nrsoln<8) nsoln=8; // we start with at least n solutions, but then gradually decrease with increasing scale
2736
2737 // !!!!!! IMPORTANT NOTE - we are inverting the order of 'this' and 'to' here to match convention in other aligners, to compensate
2738 // the Transform is inverted before being returned
2739 EMData *base_this;
2740 EMData *base_to;
2741 if (this_img->is_complex()) base_this=this_img->copy();
2742 else {
2743 base_this=this_img->do_fft();
2744 base_this->process_inplace("xform.phaseorigin.tocorner"); // This was originally .tocorner, taken from RT3DTree, but I think that's probably wrong...
2745 }
2746
2747 if (to->is_complex()) base_to=to->copy();
2748 else {
2749 base_to=to->do_fft();
2750 base_to->process_inplace("xform.phaseorigin.tocorner");
2751 }
2752
2753 int verbose = params.set_default("verbose",0);
2754 int maxshift = params.set_default("maxshift",-1);
2755 int flip = params.set_default("flip",1);
2756 float maxres = params.set_default("maxres",-1.0f);
2757 if (maxres<0.1) maxres=0.1;
2758
2759 if (base_this->get_xsize()!=base_this->get_ysize()+2 || base_to->get_xsize()!=base_to->get_ysize()+2 ) throw InvalidCallException("ERROR (RT2DTreeAligner): requires cubic images with even numbered box sizes");
2760
2761 base_this->process_inplace("xform.fourierorigin.tocenter"); // easier to chop out Fourier subvolumes
2762 base_to->process_inplace("xform.fourierorigin.tocenter");
2763
2764 float apix=(float)this_img->get_attr("apix_x");
2765 int ny=this_img->get_ysize();
2766 if (maxshift<0) maxshift=ny*3/8;
2767 float maxs = ny*apix/maxres;
2768
2769// int downsample=floor(ny/20); // Minimum shrunken box size is 20^3
2770
2771 vector<float> s_score(nsoln,0.0f);
2772 vector<Transform> s_xform(nsoln);
2773// vector<float> s_step(nsoln,7.5);
2774 if (verbose>0) printf("%d solutions\n",nsoln);
2775
2776 // We start with 32^3, 64^3 ...
2777 for (int sexp=5; sexp<10; sexp++) {
2778 int ss=pow(2.0,sexp);
2779 float rescale=2.0;
2780 if (ss>ny) {
2781 rescale=(float)ny/pow(2.0,sexp-1);
2782 ss=ny;
2783 }
2784 if (verbose>0) printf("\nSize %d\n",ss);
2785
2786 //ss=good_size(ny/ds);
2787 // Clearly these regions will be messed up on x=nyquist edge, but we are zeroing that pixel during rotation anyway
2788 EMData *small_this=base_this->get_clip(Region(0,(ny-ss)/2,ss+2,ss));
2789 EMData *small_to= base_to-> get_clip(Region(0,(ny-ss)/2,ss+2,ss));
2790 small_this->process_inplace("xform.fourierorigin.tocorner"); // after clipping back to canonical form
2791 float cut2=0.33*ny;
2792 if (maxs<cut2) cut2=maxs;
2793 small_this->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",3));
2794 small_this->process_inplace("filter.lowpass.gauss",Dict("cutoff_pixels",cut2));
2795 small_to->process_inplace("xform.fourierorigin.tocorner");
2796 small_to->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",3));
2797 small_to->process_inplace("filter.lowpass.gauss",Dict("cutoff_pixels",cut2));
2798
2799 // This is a solid estimate for very complete searching
2800 float astep = 360.0/int(M_PI/atan(1.25/ss)); // Decent angular step in degrees
2801
2802 // This insures we make at least one real effort at each level
2803 for (int i=0; i<nsoln; i++) {
2804 s_score[i]=-1.0e24; // reset the scores since the different scales will not match
2805// if (fabs(s_step[i])<astep/4.0) s_step[i]*=2.0;
2806 }
2807
2808
2809 // This is for the first loop, we do a full search in a (normally) heavily downsampled space
2810 if (sexp==5) {
2811 if (verbose>1) printf("stage 1 - ang step %1.2f, filter %1.3f\n",astep,cut2/(apix*ny));
2812 vector<Transform> transforms;
2813 for (float a=astep/2.0; a<359.9; a+=astep) {
2814 transforms.push_back(Transform(Dict("type","2d","alpha",a,"mirror",0)));
2815 if (flip) transforms.push_back(Transform(Dict("type","2d","alpha",a,"mirror",1)));
2816 }
2817 if (verbose>0) printf("%lu orientations to test\n",transforms.size());
2818
2819 // We iterate over all orientations
2820 for (unsigned int it=0; it<transforms.size(); it++) {
2821 Transform t = transforms[it];
2822
2823 // somewhat strangely, rotations are actually much more expensive than FFTs, so we use a CCF for translation
2824 EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
2825 EMData *ccf=small_to->calc_ccf(stt);
2826
2827 int lmaxs=maxshift/ss;
2828 IntPoint ml(0,0,0);
2829 float sim=ccf->get_attr("minimum");
2830 int mir=t.get_mirror()?-1:1;
2831 for (int y=-lmaxs; y<=lmaxs; y++) {
2832 for (int x=-lmaxs; x<=lmaxs; x++) {
2833 float v=ccf->get_value_at_wrap(x,y);
2834 if (v>sim) {
2835 sim=v;
2836 ml[0]=x;
2837 ml[1]=y;
2838 ml[2]=0;
2839 }
2840 }
2841 }
2842
2843// IntPoint ml=ccf->calc_max_location_wrap();
2844
2845// Dict aap=t.get_params("2d");
2846// aap["tx"]=(int)ml[0];
2847// aap["ty"]=(int)ml[1];
2848// aap["tz"]=(int)ml[2];
2849
2850 t.set_params(Dict("tx",(int)ml[0],"ty",(int)ml[1],"tz",(int)ml[2]));
2851 delete stt;
2852 delete ccf;
2853// 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
2854//
2855// sim=stt->cmp("ccc",small_to);
2856
2857 // could have used a priority queue, but would have required more infrastructure
2858 int worst=0;
2859 // First we find the worst solution in the list of possible best solutions, or the first
2860 // solution which is currently "empty"
2861//printf("a\n");
2862 for (int i=0; i<nsoln; i++) {
2863 if (s_score[i]==-1.0e24) { worst=i; break; }
2864 if (s_score[i]<s_score[worst]) worst=i;
2865 }
2866
2867 // If the current solution is better than the 'worst' of the previous solutions, then we
2868 // displace it. Note that there is no sorting performed here
2869 if (sim>s_score[worst]) {
2870 s_score[worst]=sim;
2871 s_xform[worst]=t;
2872 //printf("%f\t%f\t%d\n",s_score[worst],s_coverage[worst],worst);
2873 }
2874 }
2875
2876 if (verbose>2) {
2877 for (int i=0; i<nsoln; i++) {
2878 Dict aap=s_xform[i].get_params("2d");
2879 printf("%d) %d\t%1.1f\t%1.1f\t%1.2f\t%1.1f\n",i,(int)aap["mirror"],(float)aap["tx"],(float)aap["ty"],(float)aap["alpha"],s_score[i]);
2880 }
2881 }
2882
2883 }
2884 // Once we have our initial list of best locations, we just refine each possibility individually
2885 else {
2886 // We generate a search pattern around each existing solution
2887 if (verbose>1) printf("stage 2 (%1.2f, %1.2f, %d) filter: %1.3f\n",astep,rescale,nsoln,cut2/(apix*ny));
2888 for (int i=0; i<nsoln; i++) {
2889
2890 Dict aap=s_xform[i].get_params("2d");
2891 aap["tx"]=(float)aap["tx"]*rescale; // compensate for scaling steps, usually 2, but not in the last step
2892 aap["ty"]=(float)aap["ty"]*rescale;
2893 float alpha=aap["alpha"];
2894 for (float a=alpha-astep; a<=alpha+astep; a+=astep) {
2895 Transform t=Transform(Dict("type","2d","alpha",a,"mirror",aap["mirror"]));
2896 EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
2897 EMData *ccf=small_to->calc_ccf(stt);
2898// ccf->process_inplace("xform.phaseorigin.tocenter");
2899 int bx,by;
2900 float ba;
2901 int mir=t.get_mirror()?-1:1;
2902 int tx=aap["tx"];
2903 int ty=aap["ty"];
2904 for (int y=-4+ty; y<=4+ty; y++) {
2905 for (int x=-4+tx; x<=4+tx; x++) {
2906 float v=ccf->get_value_at_wrap(x,y);
2907 if (v>s_score[i]) {
2908 t.set_trans(x,y);
2909 s_xform[i]=t;
2910 s_score[i]=v;
2911 }
2912 }
2913 }
2914 delete stt;
2915 delete ccf;
2916 }
2917
2918 if (verbose>2) {
2919 aap=s_xform[i].get_params("2d");
2920 printf("%d) %d\t%1.1f\t%1.1f\t%1.2f\t%1.1f\n",i,(int)aap["mirror"],(float)aap["tx"],(float)aap["ty"],(float)aap["alpha"],s_score[i]);
2921 }
2922 }
2923 }
2924 // lazy earlier in defining s_ vectors, so lazy here too and inefficiently sorting
2925 // We are sorting inside the outermost loop so we can decrease the number of solutions
2926 // before we get to the finest precision
2927 for (unsigned int i=0; i<nsoln-1; i++) {
2928 for (unsigned int j=i+1; j<nsoln; j++) {
2929 if (s_score[i]<s_score[j]) {
2930 float t=s_score[i]; s_score[i]=s_score[j]; s_score[j]=t;
2931 Transform tt=s_xform[i]; s_xform[i]=s_xform[j]; s_xform[j]=tt;
2932 }
2933 }
2934 }
2935
2936 // At each level of sampling we (potentially) decrease the number of answers we check in detail
2937 // assuming we are gradually homing in on the best solution
2938 nsoln/=2;
2939 if (nsoln<nrsoln) nsoln=nrsoln;
2940
2941
2942 delete small_this;
2943 delete small_to;
2944 if (ss==ny) break;
2945 }
2946
2947 delete base_this;
2948 delete base_to;
2949
2950
2951 // initialize results
2952 vector<Dict> solns;
2953 for (unsigned int i = 0; i < nrsoln; ++i ) {
2954 Dict d;
2955 d["score"] = -s_score[i];
2956// s_xform[i].invert(); // this is because we inverted the order of the input images above to match convention
2957 if (s_xform[i].get_mirror()) {
2958 Vec3f t=s_xform[i].get_trans();
2959 t[0]*=-1.0;
2960 s_xform[i].set_trans(t);
2961 }
2962 d["xform.align2d"] = &s_xform[i];
2963 solns.push_back(d);
2964 }
2965 if (verbose>1) {
2966 Dict aap=s_xform[0].get_params("2d");
2967 printf("Final: %d\t%1.1f\t%1.1f\t%1.2f\n",(int)aap["mirror"],(float)aap["tx"],(float)aap["ty"],(float)aap["alpha"]);
2968 }
2969
2970 return solns;
2971}
EMAN::Aligner::params
Dict params
Definition: aligner.h:147
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
EMAN::Dict::set_default
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
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
EMAN::EMData::get_clip
EMData * get_clip(const Region &area, const float fill=0) const
Get an inclusive clip.
Definition: emdata.cpp:592
EMAN::EMObject
EMObject is a wrapper class for types including int, float, double, etc as defined in ObjectType.
Definition: emobject.h:123
EMAN::IntPoint
IntPoint defines an integer-coordinate point in a 1D/2D/3D space.
Definition: geometry.h:192
EMAN::Region
Region defines a 2D or 3D rectangular region specified by its origin coordinates and all edges' sizes...
Definition: geometry.h:497
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
EMAN::Transform::get_mirror
bool get_mirror() const
Query whether x_mirroring is occurring.
Definition: transform.cpp:1250
EMAN::Transform::get_trans
Vec3f get_trans() const
Get the post trans as a vec3f.
Definition: transform.cpp:1046
EMAN::Transform::set_params
void set_params(const Dict &d)
Set the parameters of the entire transform.
Definition: transform.cpp:254
EMAN::Transform::set_trans
void set_trans(const float &x, const float &y, const float &z=0)
Set the post translation component.
Definition: transform.cpp:1036
InvalidParameterException
#define InvalidParameterException(desc)
Definition: exception.h:361
InvalidCallException
#define InvalidCallException(desc)
Definition: exception.h:348
EMAN::EMData::calc_ccf
EMData * calc_ccf(EMData *with=0, fp_flag fpflag=CIRCULANT, bool center=false)
Calculate Cross-Correlation Function (CCF).
Definition: emdata.cpp:1499
y
#define y(i, j)
Definition: projector.cpp:1516
x
#define x(i)
Definition: projector.cpp:1517

References EMAN::EMData::calc_ccf(), EMAN::EMData::get_clip(), EMAN::Transform::get_mirror(), InvalidCallException, InvalidParameterException, EMAN::Aligner::params, EMAN::Dict::set_default(), EMAN::Transform::set_params(), EMAN::Transform::set_trans(), x, and y.

Member Data Documentation

◆ NAME

const string RT2DTreeAligner::NAME = "rotate_translate_tree"
static

Definition at line 1809 of file aligner.h.

Referenced by get_name().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.2