mpi_classify_CLTomo.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:    Carlos Oscar            coss@cnb.csic.es (2009)
 *
 * 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_VQ_VOLUMES
#define _PROG_VQ_VOLUMES

#include <interface/frm.h>  // must be included first as it defines _POSIX_C_SOURCE

#include <parallel/xmipp_mpi.h>
#include <core/metadata_db.h>
#include <core/metadata_extension.h>
#include <data/filters.h>
#include <data/polar.h>
#include <core/xmipp_fftw.h>
#include <core/histogram.h>
#include <data/numerical_tools.h>
#include <core/xmipp_program.h>
#include <core/symmetries.h>
#include <vector>
#include "data/mask.h"

class CL3DAssignment
{
public:
    double score;   // Negative corrCodes indicate invalid particles
    double shiftx;
    double shifty;
    double shiftz;
    double psi;
    double tilt;
    double rot;
    size_t objId;

    CL3DAssignment();

    void readAlignment(const Matrix2D<double> &M);

    void copyAlignment(const CL3DAssignment &alignment);
};

std::ostream & operator << (std::ostream &out, const CL3DAssignment& assigned);

class CL3DClass {
public:
    // Centroid
    MultidimArray<double> P, Paux, bgMask;
    
    // Centroid in Fourier space
    MultidimArray< std::complex<double> > Pfourier;

    // Fourier transformer
    FourierTransformer transformer;

    // Auxiliary Fourier image
    MultidimArray< std::complex<double> > Ifourier;

    // Auxiliary Fourier image magnitude
    MultidimArray<double> IfourierMag, IfourierMagSorted;

    // Fourier mask for the experimental image
    MultidimArray<int> IfourierMask, IfourierMaskFRM;

    // Mask for experimental image as a python object
    PyObject *pyIfourierMaskFRM=nullptr;

    // Update for next iteration
    MultidimArray< std::complex<double> > Pupdate;

    // Update for next iteration in real space
    //MultidimArray<double> PupdateReal;

    // Sum weights
    //double weightSum;

    // Update for next iteration
    MultidimArray< double > PupdateMask;

    // Auxiliary for alignment
    MultidimArray<double> Iaux;

    // List of images assigned
    std::vector<CL3DAssignment> currentListImg;

    // List of images assigned
    std::vector<CL3DAssignment> nextListImg;

    // List of neighbour indexes
    std::vector<int> neighboursIdx;

public:
    CL3DClass();

    ~CL3DClass() {}

    CL3DClass(const CL3DClass &other);
    CL3DClass(const CL3DClass &&)=delete;

    CL3DClass & operator=(const CL3DClass &other);
    CL3DClass & operator=(const CL3DClass &&)=delete;

    void updateProjection(MultidimArray<double> &I, const CL3DAssignment &assigned, bool force=false);

    void transferUpdate();

    void constructFourierMask(MultidimArray<double> &I);

    void constructFourierMaskFRM();

    void fitBasic(MultidimArray<double> &I, CL3DAssignment &result);

    void lookForNeighbours(const std::vector<CL3DClass *> listP, int K);
};

struct SDescendingClusterSort
{
     bool operator()(CL3DClass* const& rpStart, CL3DClass* const& rpEnd)
     {
          return rpStart->currentListImg.size() > rpEnd->currentListImg.size();
     }
};

class CL3D {
public:
    size_t Nimgs;

    MetaDataDb *SF;

    std::vector<CL3DClass *> P;
    
public:
    void readImage(Image<double> &I, size_t objId, bool applyGeo) const;

    void initialize(MetaDataDb &_SF,
                    std::vector< MultidimArray<double> > &_codes0);
    
    void shareAssignments(bool shareAssignment, bool shareUpdates);

    void shareSplitAssignments(Matrix1D<int> &assignment, CL3DClass *node1, CL3DClass *node2) const;

    void write(const FileName &fnRoot, int level) const;

    void lookNode(MultidimArray<double> &I, int oldnode,
                  int &newnode, CL3DAssignment &bestAssignment);
    
    void transferUpdates();

    void run(const FileName &fnOut, int level);

    int cleanEmptyNodes();

    void splitNode(CL3DClass *node,
        CL3DClass *&node1, CL3DClass *&node2,
        std::vector<size_t> &finalAssignment,
        bool iterate=true) const;

    void splitFirstNode();
};

class ProgClassifyCL3D: public XmippProgram {
public:
    FileName fnSel;
    
    FileName fnCodes0;

    FileName fnOut;

    FileName fnSym;

    int Niter;

    int Ncodes0;

    int Ncodes;

    int Nneighbours;

    double PminSize;
    
    double maxFreq;

    double sparsity;

    double DWTsparsity;

    bool classifyAllImages;

    bool randomizeStartingOrientation;

    double maxShiftZ;

    double maxShiftY;

    double maxShiftX;


    double maxRot;

    double maxTilt;

    double maxPsi;

    Mask mask;

    bool dontAlign;

    bool generateAlignedVolumes;

    // Symmetry List
    SymList SL;

    ProgClassifyCL3D(int argc, char** argv);
    ProgClassifyCL3D(const ProgClassifyCL3D&)=delete;
    ProgClassifyCL3D(const ProgClassifyCL3D&&)=delete;

    ~ProgClassifyCL3D();

    ProgClassifyCL3D & operator=(const ProgClassifyCL3D &)=delete;
    ProgClassifyCL3D & operator=(const ProgClassifyCL3D &&)=delete;

    void readParams();
    
    void show() const;
    
    void defineParams();
    
    void produceSideInfo();
    
    void run();
public:
    // Selfile with all the input images
    MetaDataDb SF;
    
    // Object Ids
    std::vector<size_t> objId;

    // Structure for the classes
    CL3D vq;

    // Mpi node
    MpiNode *node;

    // Gaussian interpolator
    GaussianInterpolator gaussianInterpolator;

    // Image dimensions
    size_t Zdim, Ydim, Xdim;

    PyObject * frmFunc;

    PyObject * wedgeClass;

    double maxShift;

    MultidimArray<unsigned char> maxFreqMask;
};
#endif