micrograph.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors: Carlos Oscar (coss@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 _CORE_XMIPPMICROGRAPH_H
#define _CORE_XMIPPMICROGRAPH_H

/* ************************************************************************* */
/* INCLUDES                                                                  */
/* ************************************************************************* */
#include <vector>
#include "filters.h"
#include <core/xmipp_funcs.h>
#include <core/multidim_array.h>
#include <core/xmipp_image.h>
#include <core/xmipp_fftw.h>
#include <core/metadata_vec.h>
#include <core/xmipp_error.h>

/* ************************************************************************* */
/* FORWARD DEFINITIONS                                                       */
/* ************************************************************************* */
// This forward definitions are needed for defining operators functions that
// use other clases type

/* ************************************************************************* */
/* MICROGRAPHY                                                               */
/* ************************************************************************* */


struct Particle_coords
{
    int label;
    int X;
    int Y;
    bool valid;
    double cost;
    double scoreVar;
    double scoreGini;
};


class Micrograph
{
public:
    struct Point
    {
        double x;
        double y;
    };
private:
    /* This image will contain a single particle from the micrograph,
       this is done to avoid asking/freeing memory all time. */
    Image<double>            single_particle;
    
    FileName                 fn_coords;
    FileName                 fn_micrograph;
    MDRowVec                 ctfRow; //MetaData row with CTF params
    FileName                 fn_inf;
    int                      X_window_size;
    int                      Y_window_size;
    int                      datatype;
    int                      swapbyte;
    int                      __offset;
    bool                     compute_transmitance;
    bool                     compute_inverse;
    int                      fh_micrograph;
    std::vector<std::string> labels;
    double                   stdevFilter = -1;
    Image<unsigned char>       IUChar = {};
    Image<short int>           IShort = {};
    Image<unsigned short int>  IUShort = {};
    Image<int>                 IInt = {};
    Image<unsigned int>        IUInt = {};
    Image<float>               IFloat = {};
public:
    size_t                   Xdim;
    size_t                   Ydim;
    size_t                   Zdim;
    size_t                   Ndim;
    Point                    point1;
    Point                    point2;
    std::vector<Particle_coords> coords;

    void clear();

    int getDatatype() const
    {
        return datatype;
    }
    int getDatatypeDetph() const;

    void open_micrograph(const FileName &fn_micrograph);

    void close_micrograph();

    const FileName& micrograph_name()
    {
        return fn_micrograph;
    }

    void set_micrograph_name(const FileName& fn)
    {
        fn_micrograph = fn;
    }

    // TODO: when needed, change this function to take MDRow&
    void set_ctfparams(const MDRowVec &ctf)
    {
        ctfRow = ctf;
        ctfRow.detach();
    }

    const MDRow& get_ctfparams()
    {
      return ctfRow;
    }

    void setStdevFilter(double d)
    {
        stdevFilter=d;
    }

    void write_coordinates(int label, double minCost, const FileName &fn_coords = "");

    void read_coordinates(int label, const FileName &fn_coords);

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

    void scale_coordinates(const double &c);

    int ParticleNo() const
    {
        return coords.size();
    }

    std::vector<Particle_coords> & Particles()
    {
        return coords;
    }

    void get_Particles(std::vector<Particle_coords> & _coords)
    {
        _coords = coords;
    }

    void set_window_size(int _X_window_size, int _Y_window_size)
    {
        X_window_size = _X_window_size;
        Y_window_size = _Y_window_size;
    }

    void set_transmitance_flag(bool flag_value)
    {
        compute_transmitance =
            flag_value;
    }

    void setDataType(DataType _datatype)
    {
        datatype = _datatype;
    }
    bool read_transmitance_flag(void)
    {
        return compute_transmitance;
    }

    void set_inverse_flag(bool flag_value)
    {
        compute_inverse =
            flag_value;
    }

    bool read_inverse_flag(void)
    {
        return compute_inverse;
    }

    template <typename T>
    int templateScissor(const Image<T> &I,
                        const Particle_coords &P, MultidimArray<double> &result,
                        double Dmin, double Dmax, double scaleX, double scaleY,
                        bool only_check, bool fillBorders)
    {
        result.initZeros(Y_window_size, X_window_size);
        int _i0 = ROUND(scaleY * P.Y) + FIRST_XMIPP_INDEX(Y_window_size);
        int _iF = ROUND(scaleY * P.Y) + LAST_XMIPP_INDEX(Y_window_size);
        int _j0 = ROUND(scaleX * P.X) + FIRST_XMIPP_INDEX(X_window_size);
        int _jF = ROUND(scaleX * P.X) + LAST_XMIPP_INDEX(X_window_size);
        int retval = 1;
        double irange=1.0/(Dmax - Dmin);

        if (!fillBorders && (_i0 < 0 || _iF >= Ydim || _j0 < 0 || _jF >= Xdim))
            retval = 0;
        else
            if (!only_check)
            {

                for (int i = _i0, i_i0=0; i <= _iF; i++, i_i0++)
                {
                    int ifrom=i;
                    if (fillBorders)
                    {
                        if (ifrom<0)
                            ifrom=0;
                        else if (ifrom>=Ydim)
                            ifrom=Ydim-1;
                    }
                    for (int j=_j0, j_j0=0; j<=_jF; j++, j_j0++)
                    {
                        int jfrom=j;
                        if (fillBorders)
                        {
                            if (jfrom<0)
                                jfrom=0;
                            else if (jfrom>=Xdim)
                                jfrom=Xdim-1;
                        }
                        double val=IMGPIXEL(I,ifrom,jfrom);
                        if (compute_transmitance)
                        {
                            double temp;
                            if (val < 1)
                                temp = val;
                            else
                                temp = log10(val);
                            if (compute_inverse)
                                A2D_ELEM(result,i_i0, j_j0) = (Dmax - temp) * irange;
                            else
                                A2D_ELEM(result,i_i0, j_j0) = (temp - Dmin) * irange;
                        }
                        else
                        {
                            if (compute_inverse)
                                A2D_ELEM(result, i_i0, j_j0) = (Dmax - val) * irange;
                            else
                                A2D_ELEM(result, i_i0, j_j0) = val;
                        }
                    }
                }
            }
        return retval;
    }

    int scissor(const Particle_coords &P, MultidimArray<double> &result,
                double Dmin, double Dmax,
                double scaleX = 1, double scaleY = 1, bool only_check = false,
                bool fillBorders = false);

    double operator()(size_t y, size_t x) const
    {
        if (y < 0 || y >= Ydim || x < 0 || x >= Xdim)
            // COSS: REPORT_ERROR(1, "Micrograph::(): index out of range");
            return 0;
        if (datatype == DT_UChar || datatype == DT_UHalfByte)
        {
            return IMGPIXEL(IUChar,y,x);
        }
        else if (datatype == DT_UShort)
        {
            return IMGPIXEL(IUShort,y,x);
        }
        else if (datatype == DT_Short)
        {
            return IMGPIXEL(IShort,y,x);
        }
        else if (datatype == DT_UInt)
        {
            return IMGPIXEL(IUInt,y,x);
        }
        else if (datatype == DT_Int)
        {
            return IMGPIXEL(IInt,y,x);
        }
        else if (datatype == DT_Float)
        {
            return IMGPIXEL(IFloat,y,x);
        }
        else
            REPORT_ERROR(ERR_TYPE_INCORRECT, "Micrograph::(): unknown datatype");
    }

    void computeDoubleMinMax(double &Dmin, double &Dmax) const
    {
        if (datatype == DT_UChar || datatype == DT_UHalfByte)
        {
            return IUChar().computeDoubleMinMax(Dmin,Dmax);
        }
        else if (datatype == DT_UShort)
        {
            return IUShort().computeDoubleMinMax(Dmin,Dmax);
        }
        else if (datatype == DT_Short)
        {
            return IShort().computeDoubleMinMax(Dmin,Dmax);
        }
        else if (datatype == DT_UInt)
        {
            return IUInt().computeDoubleMinMax(Dmin,Dmax);
        }
        else if (datatype == DT_Int)
        {
            return IInt().computeDoubleMinMax(Dmin,Dmax);
        }
        else if (datatype == DT_Float)
        {
            return IFloat().computeDoubleMinMax(Dmin,Dmax);
        }

        else
            REPORT_ERROR(ERR_TYPE_INCORRECT, "Micrograph::computeDoubleMinMax::(): unknown datatype");
    }

    //Dangerous function indeed
    //write in default endiam
    void set_val(int y, int x, double new_val)
    {
        if (datatype == DT_UChar || datatype == DT_UHalfByte)
        {
            IMGPIXEL(IUChar,y,x) = (unsigned char) new_val;
        }
        else if (datatype == DT_UShort)
        {
            IMGPIXEL(IUShort,y,x) = (unsigned short) new_val;
        }
        else if (datatype == DT_Short)
        {
            IMGPIXEL(IShort,y,x) = (short) new_val;
        }
        else if (datatype == DT_UInt)
        {
            IMGPIXEL(IUInt,y,x) = (unsigned int) new_val;
        }
        else if (datatype == DT_Int)
        {
            IMGPIXEL(IInt,y,x) = (int) new_val;
        }
        else if (datatype == DT_Float)
        {
            IMGPIXEL(IFloat,y,x) = (float) new_val;
        }

        else
            REPORT_ERROR(ERR_TYPE_INCORRECT, "Micrograph::set_val::(): unknown datatype");

    }

    void produce_all_images(int label, double minCost, const FileName &fn_root,
                            const FileName &fn_image = "", double ang = 0,
                            //double gamma = 0., double psi = 0.,
                             bool rmStack=false,
                            bool fillBorders=false,
                            bool extractNoise=false,
                            int Nnoise=-1);

    int search_coord_near(int x, int y, int prec = 3) const;

    void invalidate_coord(int n);

    int add_coord(int x, int y, int label, double cost);

    void move_last_coord_to(int x, int y);

    Particle_coords & coord(int n)
    {
        if (n < 0 || n > ParticleNo())
            REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS, "Micrograph::coord(): index out of range");
        return coords[n];
    }

    void get_coord(int n, Particle_coords &_coords)
    {
        _coords = coord(n);
    }

    int add_label(const std::string &label)
    {
        labels.push_back(label);
        return labels.size() - 1;
    }

    int LabelNo()
    {
        return labels.size();
    }

    std::string & get_label(int n)
    {
        if (n < 0 || n > LabelNo())
            REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS, "Micrograph::get_label(): index out of range");
        return labels[n];
    }

    void get_label(int n, std::string &_label)
    {
        _label = get_label(n);
    }

    void size(int &_Xdim, int &_Ydim) const
    {
        _Xdim = Xdim;
        _Ydim = Ydim;
    }

    void resize(int Xdim, int Ydim, const FileName &filename="");

    void write(const FileName &fileName,CastWriteMode castMode=CW_CAST);
};

class TiltPairAligner
{
public:
    std::vector<int> coordU;
    std::vector<int> coordT;
public:
    TiltPairAligner();

    void clear();

    void addCoordinatePair(int _muX, int _muY, int _mtX, int _mtY);

    void adjustPassingMatrix(int _muX, int _muY, int _mtX, int _mtY);

    void calculatePassingMatrix();

    void passToTilted(int _muX, int _muY, int &_mtX, int &_mtY);

    void passToUntilted(int _mtX, int _mtY, int &_muX, int &_muY);

    void computeGamma();

    void computeAngles(double &ualpha, double &talpha, double &ogamma);
public:
    // For tilted-untilted correspondance
    Matrix2D<double>    Au;     // Untilted "positions"
    Matrix2D<double>    Bt;     // Tilted   "positions"
    Matrix2D<double>    Put;    // Passing matrix from untilted to tilted
    Matrix2D<double>    Ptu;    // Passing matrix from tilted to untilted
    int                 Nu;     // Number of points in matrices
    double              gamma;  // Tilting angle in radians
    double              alpha_u;// Angle of axis with X axis in radians
    double              alpha_t;// Anfle of axis with X axis in radians
    Matrix1D<double> m;
    Matrix2D<double> pair_E;
};


#endif