tomo_align_tilt_series.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:    Carlos Oscar            coss@cnb.csic.es (2008)
 *
 * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 *
 * 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
 *
 *  All comments concerning this program package may be sent to the
 *  e-mail address 'xmipp@cnb.csic.es'
 ***************************************************************************/
#ifndef _PROG_ANGULAR_ASSIGN_FOR_TILTSERIES
#define _PROG_ANGULAR_ASSIGN_FOR_TILTSERIES

#include "core/matrix1d.h"
#include "core/matrix2d.h"
#include "core/metadata_vec.h"
#include "core/multidim_array.h"
#include "core/xmipp_program.h"


double computeAffineTransformation(const MultidimArray<unsigned char> &I1,
                                   const MultidimArray<unsigned char> &I2, int maxShift, int maxIterDE,
                                   const FileName &fn_affine,
                                   Matrix2D<double> &A12, Matrix2D<double> &A21, bool show,
                                   double thresholdAffine, bool globalAffine, bool isMirror,
                                   bool checkRotation);

class Landmark
{
public:
    double x;
    double y;
    int    imgIdx;
    Landmark& operator=(const Landmark &rhs)
    {
        x=rhs.x;
        y=rhs.y;
        imgIdx=rhs.imgIdx;
        return *this;
    }
    int operator<(const Landmark &rhs) const
    {
        return imgIdx<rhs.imgIdx;
    }
};

typedef std::vector<Landmark> LandmarkChain;

/* Forward prototype */
class Alignment;

class ProgTomographAlignment: public XmippProgram
{
public:
    FileName fnSel;

    FileName fnSelOrig;

    FileName fnRoot;

    int numThreads;

    bool globalAffine;

    bool useCriticalPoints;

    size_t Ncritical;

    size_t seqLength;

    size_t blindSeqLength;

    int gridSamples;

    double maxShiftPercentage;

    int maxIterDE;

    double psiMax;

    int maxStep;

    double deltaRot;

    double localSize;

    bool optimizeTiltAngle;

    bool isCapillar;

    bool dontNormalize;

    bool difficult;

    double corrThreshold;

    bool showAffine;

    double thresholdAffine;

    double identifyOutliersZ;

    bool doNotIdentifyOutliers;

    int pyramidLevel;

    int lastStep;

    // iteration counter as a progress measure
    int iteration;

    ~ProgTomographAlignment();

    void defineParams();

    void readParams();

    void show();

    void produceSideInfo();

    void computeAffineTransformations(bool globalAffineToUse);

    void identifyOutliers(bool mark);

    void produceInformationFromLandmarks();

    bool refineLandmark(int ii, int jj, const Matrix1D<double> &rii,
                        Matrix1D<double> &rjj, double &maxCorr, bool tryFourier) const;


    bool refineLandmark(const MultidimArray<double> &pieceii, int jj,
                        Matrix1D<double> &rjj, double actualCorrThreshold,
                        bool reversed, double &maxCorr) const;

    bool refineChain(LandmarkChain &chain, double &corrChain);

    void generateLandmarkSet();

    void writeLandmarkSet(const FileName &fnLandmark) const;

    void writeTransformations(const FileName &fnTransformations) const;

    void readLandmarkSet(const FileName &fnLandmark);

    void removeOutlierLandmarks(const Alignment &alignment);

    void alignImages(const Alignment &alignment);

    void run();

public:
    // MetaData with the input images
    MetaDataVec SF;

    // MetaData with the input images
    MetaDataVec SForig;

    // List of image pointers
    std::vector < MultidimArray<unsigned char> *> img;

    // List of mask pointers
    std::vector < MultidimArray<unsigned char> *> maskImg;

    // Index of the image closest to 0 degrees in tilt
    int iMinTilt;

    // List of tilt angles
    std::vector < double > tiltList;

    // List of names
    std::vector < std::string > name_list;

    // Average forward patch correlation
    Matrix1D<double> avgForwardPatchCorr;

    // Average backward patch correlation
    Matrix1D<double> avgBackwardPatchCorr;

    // Outlier
    Matrix1D<int> isOutlier;

    // Total number of images
    int Nimg;

    // Element [i][j] is the transformation of coordinates of i into
    // coordinates of j
    std::vector< std::vector< Matrix2D<double> > > affineTransformations;

    // List of affine costs
    std::vector < double > correlationList;

    // Landmark matrix (X component)
    Matrix2D<double> allLandmarksX;

    // Landmark matrix (Y component)
    Matrix2D<double> allLandmarksY;

    // Set of all landmarks seen in image i
    std::vector< std::vector<int> > Vseti;

    // Set of all images in which the landmark j is seen
    std::vector< std::vector<int> > Vsetj;

    // Average of the landmarks in a given image
    std::vector< Matrix1D<double> > barpi;

    // Number of landmarks in image i
    Matrix1D<int> ni;

    // Best exhaustive alignmnet
    Alignment* bestPreviousAlignment;

    // Show refinement
    bool showRefinement;
};

class Alignment
{
public:
    const ProgTomographAlignment *prm;
    std::vector< Matrix1D<double> > di;
    std::vector< Matrix1D<double> > rj;
    Matrix1D<double> psi;
    double rot;
    double tilt;
    Matrix1D<double> raxis;

public:
    Alignment(const ProgTomographAlignment *_prm)
    {
        prm=_prm;
        Nimg=MAT_XSIZE(_prm->allLandmarksX);
        Nlandmark=MAT_YSIZE(_prm->allLandmarksX);
        clear();
    }

    void clear()
    {
        psi.initZeros(Nimg);
        rot=90;
        tilt=90;
        raxis.initZeros(3);

        Ai.clear();
        Ait.clear();
        Aip.clear();
        Aipt.clear();
        di.clear();
        diaxis.clear();
        B1i.clear();
        B2i.clear();
        barri.clear();
        Matrix2D<double> dummy23(2,3);
        Matrix2D<double> dummy32(3,2);
        Matrix2D<double> dummy33(3,3);
        Matrix2D<double> dummy22(2,2);
        Matrix1D<double> dummy2(2);
        Matrix1D<double> dummy3(3);
        for (int i=0; i<Nimg; i++)
        {
            Ai.push_back(dummy23);
            Ait.push_back(dummy32);
            Aip.push_back(dummy23);
            Aipt.push_back(dummy32);
            di.push_back(dummy2);
            diaxis.push_back(dummy2);
            B1i.push_back(dummy33);
            B2i.push_back(dummy33);
            barri.push_back(dummy3);
        }
        allLandmarksPredictedX.initZeros(Nlandmark,Nimg);
        allLandmarksPredictedY.initZeros(Nlandmark,Nimg);
        errorLandmark.initZeros(Nlandmark);

        rj.clear();
        for (int j=0; j<Nlandmark; j++)
            rj.push_back(dummy3);
    }

    Alignment & operator= (const Alignment &op)
    {
        if (this!=&op)
        {
            prm=op.prm;
            Nimg=op.Nimg;
            Nlandmark=op.Nlandmark;
            psi=op.psi;
            rot=op.rot;
            tilt=op.tilt;
            Ai=op.Ai;
            Ait=op.Ait;
            Aip=op.Aip;
            Aipt=op.Aipt;
            diaxis=op.diaxis;
            B1i=op.B1i;
            B2i=op.B2i;
            raxis=op.raxis;
            di=op.di;
            rj=op.rj;
            barri=op.barri;
            allLandmarksPredictedX=op.allLandmarksPredictedX;
            allLandmarksPredictedY=op.allLandmarksPredictedY;
            errorLandmark=op.errorLandmark;
        }
        return *this;
    }

    double optimizeGivenAxisDirection();

    void computeGeometryDependentOfAxis();

    void computeGeometryDependentOfRotation();

    double computeError() const;

    void computeErrorForLandmarks();

    void updateModel();

    friend std::ostream& operator<<(std::ostream &out,
                                    Alignment &alignment);

public:
    // Number of images
    int Nimg;

    // Number of landmarks
    int Nlandmark;

    // Set of Ai matrices
    std::vector< Matrix2D<double> > Ai;

    // Set of Ait matrices
    std::vector< Matrix2D<double> > Ait;

    // Set of Ai prime matrices
    std::vector< Matrix2D<double> > Aip;

    // Set of Ai prime transpose matrices
    std::vector< Matrix2D<double> > Aipt;

    // Set of bar ri
    std::vector< Matrix1D<double> > barri;

    // Set of shifts due to raxis
    std::vector< Matrix1D<double> > diaxis;

    // Set of B1i matrices
    std::vector< Matrix2D<double> > B1i;

    // Set of B2i matrices
    std::vector< Matrix2D<double> > B2i;

    // Matrix used for updating raxis
    Matrix2D<double> Binvraxis;

    // Set of predicted landmarks (component X)
    Matrix2D<double> allLandmarksPredictedX;

    // Set of predicted landmarks (component Y)
    Matrix2D<double> allLandmarksPredictedY;

    // Set of errors associated to each landmark
    MultidimArray<double> errorLandmark;
};

double computeAffineTransformation(const MultidimArray<double> &I1,
                                   const MultidimArray<double> &I2, int maxShift, int maxIterDE,
                                   Matrix2D<double> &A12, Matrix2D<double> &A21, bool show,
                                   double thresholdAffine, bool globalAffine, bool isMirror,
                                   int pyramidLevel);
#endif