aligner.h

Go to the documentation of this file.

/*
 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
 * Copyright (c) 2000-2006 Baylor College of Medicine
 *
 * This software is issued under a joint BSD/GNU license. You may use the
 * source code in this file under either license. However, note that the
 * complete EMAN2 and SPARX software packages have some GPL dependencies,
 * so you are responsible for compliance with the licenses of these packages
 * if you opt to use BSD licensing. The warranty disclaimer below holds
 * in either instance.
 *
 * This complete copyright notice must be included in any revised version of the
 * source code. Additional authorship citations may be added, but existing
 * author citations must be preserved.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 * */

#ifndef eman__aligner_h__
#define eman__aligner_h__ 1


#include "emobject.h"


namespace EMAN
{
        class EMData;
        class Cmp;

        class Aligner
        {
          public:
                virtual ~ Aligner()
                {
                }

                virtual EMData *align(EMData * this_img, EMData * to_img) const = 0;

                virtual EMData *align(EMData * this_img, EMData * to_img,
                                                          const string & cmp_name, const Dict& cmp_params) const = 0;

                virtual string get_name() const = 0;

                virtual string get_desc() const = 0;
                virtual Dict get_params() const
                {
                        return params;
                }

                virtual void set_params(const Dict & new_params)
                {
                        params = new_params;
                }

                virtual TypeDict get_param_types() const = 0;

                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 {
                        vector<Dict> solns;
                        return solns;
                }

          protected:
                mutable Dict params;

//              /** Get a Transform pointer that is currently stored in the image header corresponding to the given key.
//               * If non existant then initialize a new Transform pointer, set it as the image attribute, and return it.
//               * The point being that if you have to access the xform.align2d or the xform.align3d
//               * attributes from multiple aligners then you always want to get the same one.
//               * @param key the alignment key such as "xform.align2d" or "xform.align3d"
//               * @param image the image from which the Transform pointer will be extracted (and potentially assigned to if non existant)
//               * @return the Transform pointer that is currently stored in the image header corresponding to the given key.
//               */
//              static Transform* get_align_attr(const string& key, EMData* const image );
//
//              static Transform* get_set_align_attr(const string& key, EMData* const to_image, const EMData* const from_image  );
        };

        class TranslationalAligner:public Aligner
        {
          public:
                virtual EMData * align(EMData * this_img, EMData * to_img,
                                                const string & cmp_name="dot", const Dict& cmp_params = Dict()) const;

                virtual EMData * align(EMData * this_img, EMData * to_img) const
                {
                        return align(this_img, to_img, "", Dict());
                }

                virtual string get_name() const
                {
                        return "translational";
                }

                virtual string get_desc() const
                {
                        return "Translational 2D and 3D alignment by cross-correlation";
                }

                static Aligner *NEW()
                {
                        return new TranslationalAligner();
                }

                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("intonly", EMObject::INT,"Integer pixel translations only");
                        d.put("maxshift", EMObject::INT,"Maximum translation in pixels");
                        d.put("masked", EMObject::INT,"Treat zero pixels in 'this' as a mask for normalization (default false)");
                        d.put("nozero", EMObject::INT,"Zero translation not permitted (useful for CCD images)");
                        return d;
                }
        };

        class RotationalAligner:public Aligner
        {
          public:
                virtual EMData * align(EMData * this_img, EMData * to_img,
                                                const string & cmp_name = "dot", const Dict& cmp_params = Dict()) const;

                virtual EMData * align(EMData * this_img, EMData * to_img) const
                {
                        return align(this_img, to_img, "", Dict());
                }

                virtual string get_name() const
                {
                        return "rotational";
                }

                virtual string get_desc() const
                {
                        return "Performs rotational alignment,works accurately if the image is precentered, normally called internally in combination with translational and flip alignment";
                }

                static Aligner *NEW()
                {
                        return new RotationalAligner();
                }

                static EMData * align_180_ambiguous(EMData * this_img, EMData * to_img, int rfp_mode = 0);

                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("rfp_mode", EMObject::INT,"Either 0,1 or 2. A temporary flag for testing the rotational foot print. O is the original eman1 way. 1 is just using calc_ccf without padding. 2 is using calc_mutual_correlation without padding.");
                        return d;
                }
        };

        class RotatePrecenterAligner:public Aligner
        {
          public:
                virtual EMData * align(EMData * this_img, EMData * to_img,
                                                const string & cmp_name = "dot", const Dict& cmp_params = Dict()) const;

                virtual EMData * align(EMData * this_img, EMData * to_img) const
                {
                        return align(this_img, to_img, "", Dict());
                }

                virtual string get_name() const
                {
                        return "rotate_precenter";
                }

                virtual string get_desc() const
                {
                        return "Performs rotational alignment and works accurately if the image is precentered";
                }

                static Aligner *NEW()
                {
                        return new RotatePrecenterAligner();
                }

                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        return d;
                }
        };

        class RotateTranslateAligner:public Aligner
        {
          public:
                  virtual EMData * align(EMData * this_img, EMData * to_img,
                                           const string & cmp_name="dot", const Dict& cmp_params = Dict()) const;

                virtual EMData * align(EMData * this_img, EMData * to_img) const
                {
                        return align(this_img, to_img, "sqeuclidean", Dict());
                }

                virtual string get_name() const
                {
                        return "rotate_translate";
                }

                virtual string get_desc() const
                {
                        return "Performs rotational alignment and follows this with translational alignment.";
                }

                static Aligner *NEW()
                {
                        return new RotateTranslateAligner();
                }

                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        //d.put("usedot", EMObject::INT);
                        d.put("maxshift", EMObject::INT, "Maximum translation in pixels");
                        d.put("nozero", EMObject::INT,"Zero translation not permitted (useful for CCD images)");
                        d.put("rfp_mode", EMObject::INT,"Either 0,1 or 2. A temporary flag for testing the rotational foot print");
                        return d;
                }
        };

        class RotateTranslateBestAligner:public Aligner
        {
          public:
                virtual EMData * align(EMData * this_img, EMData * to_img,
                                           const string & cmp_name="dot", const Dict& cmp_params = Dict()) const;

                virtual EMData * align(EMData * this_img, EMData * to_img) const
                {
                        return align(this_img, to_img, "frc", Dict());
                }

                virtual string get_name() const
                {
                        return "rotate_translate_best";
                }

                virtual string get_desc() const
                {
                        return "Full 2D alignment using 'Rotational' and 'Translational', also incorporates 2D 'Refine' alignments.";
                }

                static Aligner *NEW()
                {
                        return new RotateTranslateBestAligner();
                }

                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("maxshift", EMObject::INT, "Maximum translation in pixels");
                        d.put("snr", EMObject::FLOATARRAY, "signal to noise ratio array");
                        return d;
                }


        };

        class RotateFlipAligner:public Aligner
        {
          public:
                virtual EMData * align(EMData * this_img, EMData * to_img,
                                           const string & cmp_name="dot", const Dict& cmp_params = Dict()) const;
                virtual EMData * align(EMData * this_img, EMData * to_img) const
                {
                        return align(this_img, to_img, "", Dict());
                }
                virtual string get_name() const
                {
                        return "rotate_flip";
                }

                virtual string get_desc() const
                {
                        return "Performs two rotational alignments, one using the original image and one using the hand-flipped image. Decides which alignment is better using a comparitor and returns it";
                }

                static Aligner *NEW()
                {
                        return new RotateFlipAligner();
                }

                virtual TypeDict get_param_types() const
                {
                        return static_get_param_types();
                }

                static TypeDict static_get_param_types() {
                        TypeDict d;

                        d.put("imask", EMObject::INT);
                        d.put("rfp_mode", EMObject::INT,"Either 0,1 or 2. A temporary flag for testing the rotational foot print");
                        return d;
                }

        };

        class RotateTranslateFlipAligner:public Aligner
        {
          public:
                virtual EMData * align(EMData * this_img, EMData * to_img,
                                           const string & cmp_name="dot", const Dict& cmp_params = Dict()) const;
                virtual EMData * align(EMData * this_img, EMData * to_img) const
                {
                        return align(this_img, to_img, "sqeuclidean", Dict());
                }

                virtual string get_name() const
                {
                        return "rotate_translate_flip";
                }

                virtual string get_desc() const
                {
                        return " Does two 'rotate_translate' alignments, one to accommodate for possible handedness change. Decided which alignment is better using a comparitor and returns the aligned image as the solution";
                }

                static Aligner *NEW()
                {
                        return new RotateTranslateFlipAligner();
                }

                virtual TypeDict get_param_types() const
                {
                        return static_get_param_types();
                }

                static TypeDict static_get_param_types() {
                        TypeDict d;

                        d.put("flip", EMObject::EMDATA);
                        d.put("usedot", EMObject::INT);
                        d.put("maxshift", EMObject::INT, "Maximum translation in pixels");
                        d.put("rfp_mode", EMObject::INT,"Either 0,1 or 2. A temporary flag for testing the rotational foot print");
                        return d;
                }
        };

        class RTFExhaustiveAligner:public Aligner
        {
          public:
                virtual EMData * align(EMData * this_img, EMData * to_img,
                                           const string & cmp_name="dot", const Dict& cmp_params = Dict()) const;
                virtual EMData * align(EMData * this_img, EMData * to_img) const
                {
                        return align(this_img, to_img, "sqeuclidean", Dict());
                }

                virtual string get_name() const
                {
                        return "rtf_exhaustive";
                }

                virtual string get_desc() const
                {
                        return "Experimental full 2D alignment with handedness check using semi-exhaustive search (not necessarily better than RTFBest)";
                }

                static Aligner *NEW()
                {
                        return new RTFExhaustiveAligner();
                }

                virtual TypeDict get_param_types() const
                {
                        TypeDict d;

                        d.put("flip", EMObject::EMDATA);
                        d.put("maxshift", EMObject::INT, "Maximum translation in pixels");
                        return d;
                }
        };

        class RTFSlowExhaustiveAligner:public Aligner
        {
          public:
                virtual EMData * align(EMData * this_img, EMData * to_img,
                                                const string & cmp_name, const Dict& cmp_params) const;
                virtual EMData * align(EMData * this_img, EMData * to_img) const
                {
                        return align(this_img, to_img, "sqeuclidean", Dict());
                }
                virtual string get_name() const
                {
                        return "rtf_slow_exhaustive";
                }

                virtual string get_desc() const
                {
                        return "Experimental full 2D alignment with handedness check using more exhaustive search (not necessarily better than RTFBest)";
                }

                static Aligner *NEW()
                {
                        return new RTFSlowExhaustiveAligner();
                }

                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("flip", EMObject::EMDATA,"Optional. This is the flipped version of the images that is being aligned. If specified it will be used for the handedness check, if not a flipped copy of the image will be made");
                        d.put("maxshift", EMObject::INT,"The maximum length of the detectable translational shift");
                        d.put("transtep", EMObject::FLOAT,"The translation step to take when honing the alignment, which occurs after coarse alignment");
                        d.put("angstep", EMObject::FLOAT,"The angular step (in degrees) to take in the exhaustive search for the solution angle. Typically very small i.e. 3 or smaller.");
                        return d;
                }
        };

        class RefineAligner:public Aligner
        {
          public:
                virtual EMData * align(EMData * this_img, EMData * to_img,
                                           const string & cmp_name="dot", const Dict& cmp_params = Dict()) const;

                virtual EMData * align(EMData * this_img, EMData * to_img) const
                {
                        return align(this_img, to_img, "sqeuclidean", Dict());
                }

                virtual string get_name() const
                {
                        return "refine";
                }

                virtual string get_desc() const
                {
                        return "Refines a preliminary 2D alignment using a simplex algorithm. Subpixel precision.";
                }

                static Aligner *NEW()
                {
                        return new RefineAligner();
                }

                virtual TypeDict get_param_types() const
                {
                        TypeDict d;

                        d.put("mode", EMObject::INT, "Currently unused");
                        d.put("xform.align2d", EMObject::TRANSFORM, "The Transform storing the starting guess. If unspecified the identity matrix is used");
                        d.put("stepx", EMObject::FLOAT, "The x increment used to create the starting simplex. Default is 1");
                        d.put("stepy", EMObject::FLOAT, "The y increment used to create the starting simplex. Default is 1");
                        d.put("stepaz", EMObject::FLOAT, "The rotational increment used to create the starting simplex. Default is 5");
                        d.put("precision", EMObject::FLOAT, "The precision which, if achieved, can stop the iterative refinement before reaching the maximum iterations. Default is 0.04.");
                        d.put("maxiter", EMObject::INT,"The maximum number of iterations that can be performed by the Simplex minimizer");
                        d.put("maxshift", EMObject::INT,"Maximum translation in pixels in any direction. If the solution yields a shift beyond this value in any direction, then the refinement is judged a failure and the original alignment is used as the solution.");
                        return d;
                }
        };


        class Refine3DAligner:public Aligner
        {
                public:
                        virtual EMData * align(EMData * this_img, EMData * to_img,
                                                   const string & cmp_name="dot", const Dict& cmp_params = Dict()) const;

                        virtual EMData * align(EMData * this_img, EMData * to_img) const
                        {
                                return align(this_img, to_img, "sqeuclidean", Dict());
                        }

                        virtual string get_name() const
                        {
                                return "refine.3d";
                        }

                        virtual string get_desc() const
                        {
                                return "Refines a preliminary 3D alignment using a simplex algorithm. Subpixel precision.";
                        }

                        static Aligner *NEW()
                        {
                                return new Refine3DAligner();
                        }

                        virtual TypeDict get_param_types() const
                        {
                                TypeDict d;
                                d.put("xform.align3d", EMObject::TRANSFORM,"The Transform storing the starting guess. If unspecified the identity matrix is used");
                                d.put("stepx", EMObject::FLOAT, "The x increment used to create the starting simplex. Default is 1");
                                d.put("stepy", EMObject::FLOAT,"The y increment used to create the starting simplex. Default is 1");
                                d.put("stepz", EMObject::FLOAT, "The z increment used to create the starting simplex. Default is 1." );
                                d.put("stepphi", EMObject::FLOAT, "The phi incremenent used to creat the starting simplex. This is the increment applied to the inplane rotation. Default is 5." );
                                d.put("stepdelta", EMObject::FLOAT,"The angular increment which represents a good initial step along the sphere, thinking in terms of quaternions. Default is 5.");
                                d.put("precision", EMObject::FLOAT, "The precision which, if achieved, can stop the iterative refinement before reaching the maximum iterations. Default is 0.04." );
                                d.put("maxiter", EMObject::INT, "The maximum number of iterations that can be performed by the Simplex minimizer. Default is 60.");
                                d.put("maxshift", EMObject::INT,"Maximum translation in pixels in any direction. If the solution yields a shift beyond this value in any direction, then the refinement is judged a failure and the original alignment is used as the solution.");
                                return d;
                        }
        };

        class RT3DGridAligner:public Aligner
        {
                public:
                        virtual EMData * align(EMData * this_img, EMData * to_img,
                                                   const string & cmp_name, const Dict& cmp_params) const;
                        virtual EMData * align(EMData * this_img, EMData * to_img) const
                        {
                                return align(this_img, to_img, "dot.tomo", Dict());
                        }


                        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;

                        virtual string get_name() const
                        {
                                return "rt.3d.grid";
                        }

                        virtual string get_desc() const
                        {
                                return "3D rotational and translational alignment using specified ranges and maximum shifts";
                        }

                        static Aligner *NEW()
                        {
                                return new RT3DGridAligner();
                        }

                        virtual TypeDict get_param_types() const
                        {
                                TypeDict d;
                                d.put("daz", EMObject::FLOAT,"The angle increment in the azimuth direction. Default is 10");
                                d.put("raz", EMObject::FLOAT,"The range of angles to sample in the azimuth direction. Default is 360.");
                                d.put("dphi", EMObject::FLOAT,"The angle increment in the phi direction. Default is 10.");
                                d.put("rphi", EMObject::FLOAT,"The range of angles to sample in the phi direction. Default is 180.");
                                d.put("dalt", EMObject::FLOAT,"The angle increment in the altitude direction. Default is 10.");
                                d.put("ralt", EMObject::FLOAT,"The range of angles to sample in the altitude direction. Default is 180.");
                                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;
                        }
        };

        class RT3DSphereAligner:public Aligner
        {
                public:
                        virtual EMData * align(EMData * this_img, EMData * to_img,
                                                                   const string & cmp_name, const Dict& cmp_params) const;
                        virtual EMData * align(EMData * this_img, EMData * to_img) const
                        {
                                return align(this_img, to_img, "sqeuclidean", Dict());
                        }


                        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;

                        virtual string get_name() const
                        {
                                return "rt.3d.sphere";
                        }

                        virtual string get_desc() const
                        {
                                return "3D rotational and translational alignment using spherical sampling. Can reduce the search space if symmetry is supplied";
                        }

                        static Aligner *NEW()
                        {
                                return new RT3DSphereAligner();
                        }

                        virtual TypeDict get_param_types() const
                        {
                                TypeDict d;
                                d.put("sym", EMObject::STRING,"The symmtery to use as the basis of the spherical sampling. Default is c1 (asymmetry).");
                                d.put("orientgen", EMObject::STRING,"Advanced. The orientation generation strategy. Default is eman");
                                d.put("delta", EMObject::FLOAT,"Angle the separates points on the sphere. This is exclusive of the \'n\' paramater. Default is 10");
                                d.put("n", EMObject::INT,"An alternative to the delta argument, this is the number of points you want generated on the sphere. Default is OFF");
                                d.put("dphi", EMObject::FLOAT,"The angle increment in the phi direction. Default is 10.");
                                d.put("rphi", EMObject::FLOAT,"The range of angles to sample in the phi direction. Default is 180.");
                                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;
                        }
        };


        class CUDA_Aligner
        {
          public:
                CUDA_Aligner();

                void finish();

                void setup(int nima, int nx, int ny, int ring_length, int nring, int ou, float step, int kx, int ky, bool ctf);

                void insert_image(EMData *image, int num);

                void filter_stack(vector<float> ctf_params, int id);
                
                void sum_oe(vector<float> ctf_params, vector<float> ali_params, EMData* ave1, EMData *ave2, int id);

                vector<float> alignment_2d(EMData *ref_image, vector<float> sx, vector<float> sy, int id, int silent);

          private:
                float *image_stack, *image_stack_filtered;
                float *ccf;
                int NIMA, NX, NY, RING_LENGTH, NRING, OU, KX, KY;
                bool CTF;
                float STEP;
        };


        template <> Factory < Aligner >::Factory();

        void dump_aligners();
        map<string, vector<string> > dump_aligners_list();
}

#endif