forward_zernike_subtomos.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:    Carlos Oscar Sanchez Sorzano (coss@cnb.csic.es)
 *             David Herreros Calero (dherreros@cnb.csic.es)
 *
 * 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_FORWARD_ZERNIKE_SUBTOMOS
#define _PROG_FORWARD_ZERNIKE_SUBTOMOS

#include "core/xmipp_metadata_program.h"
#include "core/rerunable_program.h"
#include "core/matrix1d.h"
#include <data/blobs.h>
#include "core/xmipp_image.h"
#include "data/fourier_filter.h"
#include "data/fourier_projection.h"


class ProgForwardZernikeSubtomos: public XmippMetadataProgram, public Rerunable
{
public:
    FileName fnVolR;
    FileName fnMaskR;
    FileName fnOutDir;
    int L1, L2;
    Matrix1D<int> vL1, vN, vL2, vM;
    double maxShift;
    double maxAngularChange;
    double maxResol;
    double Ts;
    int Rmax;
    // Optimize alignment
    bool optimizeAlignment;
    //Optimize deformation
    bool optimizeDeformation;
    //Optimize defocus
    bool optimizeDefocus;
    // Ignore CTF
    bool ignoreCTF;
    //Radius optimization
    bool optimizeRadius;
    // Phase Flipped
    bool phaseFlipped;
    // Regularization weight
    double lambda;
    // Maximum radius for the deformation
    int RmaxDef;
    // Number of images (execution mode) (single image == 1, pair == 2, or triplet == 3)
    int num_images;
    // Number alignment parameters to minimize
    int algn_params;
    // Number CTF parameters to minimize
    int ctf_params;
    Matrix1D<double> p;
    int flagEnabled;
    bool useCTF;
    // Tilt angles
    double t1, t2;

public:
    bool resume;
    // 2D and 3D masks in real space
    MultidimArray<int> mask3D, V_mask;
    // Volume size
    size_t Xdim;
    // Images Filename
    FileName fnImage;
    // Input image
    Image<double> V, Vdeformed;
    Image<double> I;
    Image<double> Ifiltered;
    Image<double> Ifilteredp;
    // Theoretical subtomos
    Image<double> P;
    // Filter
    FourierFilter filter;
    FourierFilter filterMW;
    // Transformation matrix
    Matrix2D<double> A;
    // Original angles
    double old_rot, old_tilt, old_psi, deltaRot, deltaTilt, deltaPsi;
    // Original shift
    double old_shiftX, old_shiftY, old_shiftZ, deltaX, deltaY, deltaZ;
    // Original flip
    bool old_flip;
    // CTF Check
    bool hasCTF;
    // Original defocus
    double old_defocusU, old_defocusV, old_defocusAngle, deltaDefocusU, deltaDefocusV, deltaDefocusAngle;
    // Current defoci
    double currentDefocusU, currentDefocusV, currentAngle;
    // CTF filter
    FourierFilter FilterCTF;
    // Vector Size
    int vecSize;
    // Vector containing the degree of the spherical harmonics
    Matrix1D<double> clnm, prev_clnm;
    //Copy of Optimizer steps
    Matrix1D<double> steps_cp;
    //Total Deformation, sumV, sumVd
    double totalDeformation, prior_deformation;
    int sumV;
    // Show optimization
    bool showOptimization;
    // Correlation
    double correlation;
    // Loop step
    int loop_step;
    // Blob
    struct blobtype blob;
    double blob_r;

    double sigma4;
    // Gaussian projection table
    Matrix1D<double> gaussianProjectionTable;

    // Gaussian projection2 table
    Matrix1D<double> gaussianProjectionTable2;

    // Deformation field and positions
    MultidimArray<double> vpos, df;
    std::vector<size_t> idx_z_clnm;
    std::vector<double> z_clnm_diff;
    Matrix2D<double> R;

public:
    enum class Direction { ROTATE, UNROTATE };
    
    ProgForwardZernikeSubtomos();

    ~ProgForwardZernikeSubtomos();

    void readParams();

    void show();

    void defineParams();

    void preProcess();

    void processImage(const FileName &fnImg, const FileName &fnImgOut, const MDRow &rowIn, MDRow &rowOut);

    void numCoefficients(int l1, int l2, int &vecSize);

    // void minimizepos(Matrix1D<double> &vectpos);
    void minimizepos(int L1, int l2, Matrix1D<double> &steps);

    void fillVectorTerms(int l1, int l2, Matrix1D<int> &vL1, Matrix1D<int> &vN, 
                         Matrix1D<int> &vL2, Matrix1D<int> &vM);

    void deformVol(MultidimArray<double> &mVD, const MultidimArray<double> &mV, double &def,
                   double rot, double tilt, double psi);

    void updateCTFImage(double defocusU, double defocusV, double angle);

    double transformImageSph(double *pclnm);

    template<Direction DIRECTION>
    void rotateCoefficients();

    //AJ new
    virtual void finishProcessing();

    virtual void writeImageParameters(MDRow &row);
    //END AJ

    virtual void checkPoint();

    // void removePixels();

    Matrix1D<double> weightsInterpolation3D(double x, double y, double z);
    
    void splattingAtPos(std::array<double, 3> r, double weight, MultidimArray<double> &mP, const MultidimArray<double> &mV);

    void rotatePositions(double rot, double tilt, double psi);

    void preComputeDF();

protected:
    void createWorkFiles() { return Rerunable::createWorkFiles(resume, getInputMd()); }

private:
    std::vector<MDLabel> getLabelsForEmpty() override
    {
        std::vector<MDLabel> labels = getInputMd()->getActiveLabels();
        labels.push_back(MDL_SPH_DEFORMATION);
        labels.push_back(MDL_SPH_COEFFICIENTS);
        if (!(std::find(labels.begin(), labels.end(), MDL_ENABLED) != labels.end()))
            labels.push_back(MDL_ENABLED);
        // labels.push_back(MDL_COST);
        return labels;
    }
};
#endif