mpi_classify_CL2D.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_PROJECTIONS
#define _PROG_VQ_PROJECTIONS

#include "core/metadata_db.h"
#include "data/polar.h"
#include "core/histogram.h"
#include "data/numerical_tools.h"
#include "core/xmipp_program.h"

class MpiNode;
template<typename T>
class Image;

class CL2DAssignment
{
public:
    double corr;   // Negative corrCodes indicate invalid particles
    double likelihood; // Only valid if robust criterion
    double shiftx;
    double shifty;
    double psi;
    size_t objId;
    bool flip;

    CL2DAssignment();

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

    void copyAlignment(const CL2DAssignment &alignment);
};

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

class CL2DClass {
public:
    // Projection
    MultidimArray<double> P;
    
    // Update for next iteration
    MultidimArray<double> Pupdate;

    // Polar Fourier transform of the projection at full size
    Polar<std::complex <double> > polarFourierP;

    // Rotational correlation for best_rotation
    MultidimArray<double> rotationalCorr;

    // Plans for the best_rotation
    Polar_fftw_plans *plans;

    // Correlation aux
    CorrelationAux corrAux;

    // Rotational correlation aux
    RotationalCorrelationAux rotAux;

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

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

    // Correlations of the next non-class members
    std::vector<double> nextNonClassCorr;

    // Histogram of the correlations of the current class members
    Histogram1D histClass;

    // Histogram of the correlations of the current non-class members
    Histogram1D histNonClass;

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

    CL2DClass();

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

    ~CL2DClass();

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

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

    inline void updateNonProjection(double corr, bool force=false)
    {
        if (corr>0 || force)
            nextNonClassCorr.push_back(corr);
    }

    void transferUpdate(bool centerReference=true);

    void fitBasic(MultidimArray<double> &I, CL2DAssignment &result,  bool reverse=false);

    void fit(MultidimArray<double> &I, CL2DAssignment &result);

    void lookForNeighbours(const std::vector<CL2DClass *> listP, int K);
private:
    Polar<double> fitBasic_aux;
};

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

class CL2D {
public:
    size_t Nimgs=0;

    MetaDataDb *SF=nullptr;

    std::vector<CL2DClass *> P;

    CL2D() {}

    ~CL2D();

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

    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, bool shareNonCorr);

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

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

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

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

    int cleanEmptyNodes();

    void splitNode(CL2DClass *node,
        CL2DClass *&node1, CL2DClass *&node2,
        std::vector<size_t> &finalAssignment) const;

    void splitFirstNode();
};

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

    FileName fnOut;

    FileName fnODir;

    int Niter;

    int Ncodes0;

    int Ncodes;

    int Nneighbours;

    double PminSize;
    
    bool useCorrelation;

    bool classicalMultiref;
    
    bool classifyAllImages;

    bool classicalSplit;

    int NSplitTrials;

    double maxShift;

    bool normalizeImages;

    bool mirrorImages;

    bool useThresholdMask;

    double threshold;

    bool alignImages;

    ProgClassifyCL2D(int argc, char** argv);

    ~ProgClassifyCL2D();

    ProgClassifyCL2D(const ProgClassifyCL2D &)=delete;
    ProgClassifyCL2D(const ProgClassifyCL2D &&)=delete;
    ProgClassifyCL2D & operator =(const ProgClassifyCL2D &)=delete;
    ProgClassifyCL2D & operator =(const ProgClassifyCL2D &&)=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
    CL2D vq;

    // Mpi node
    MpiNode *node;

    // Maxshift squared
    double maxShift2;

    // Gaussian interpolator
    GaussianInterpolator gaussianInterpolator;

    // Image dimensions
    size_t Ydim, Xdim;

    MultidimArray<int> mask;

    double sigma;
};
#endif