mask.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:     Carlos Oscar S. Sorzano (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_MASK_H
#define CORE_MASK_H

#include "blobs.h"
#include "core/histogram.h"
#include "core/xmipp_metadata_program.h"

template<typename T>
class MultidimArray;
template<typename T>
class Matrix2D;

void apply_geo_binary_2D_mask(MultidimArray< int > &mask,
                              const Matrix2D< double >& A);
void apply_geo_cont_2D_mask(MultidimArray< double >& mask,
                            const Matrix2D< double >& A);



constexpr int INNER_MASK = 1;
constexpr int OUTSIDE_MASK = 2;
constexpr int NO_ACTIVATE = 0;
constexpr int ACTIVATE = 1;



void RaisedCosineMask(MultidimArray< double >& mask,
                      double r1, double r2, int mode = INNER_MASK, double x0 = 0,
                      double y0 = 0, double z0 = 0);

void RaisedCrownMask(MultidimArray< double >& mask,
                     double r1, double r2, double pix_width,
                     int mode = INNER_MASK,
                     double x0 = 0, double y0 = 0, double z0 = 0);

void KaiserMask(MultidimArray<double> &mask, double delta = 0.01,
                double Deltaw = 1.0 / 12.0);

void SincMask(MultidimArray< double >& mask,
              double omega, int mode = INNER_MASK, double x0 = 0, double y0 = 0, double z0 = 0);

void SincKaiserMask(MultidimArray<double> &mask,
                    double omega, double delta = 0.01, double Deltaw = 1.0 / 12.0);

void BlackmanMask(MultidimArray< double >& mask, int mode = INNER_MASK,
                  double x0 = 0, double y0 = 0, double z0 = 0);

void SincBlackmanMask(MultidimArray< double >& mask,
                      double omega, double power_percentage,
                      double x0 = 0, double y0 = 0, double z0 = 0);

void BinaryCircularMask(MultidimArray< int >& mask,
                        double radius, int mode = INNER_MASK,
                        double x0 = 0, double y0 = 0, double z0 = 0);


void BlobCircularMask(MultidimArray<double> &mask,
                      double r1, blobtype blob, int mode,
                      double x0 = 0, double y0 = 0, double z0 = 0);


void BinaryCrownMask(MultidimArray< int >& mask,
                     double R1, double R2, int mode = INNER_MASK,
                     double x0 = 0, double y0 = 0, double z0 = 0);

void BlobCrownMask(MultidimArray<double> &mask,
                   double r1, double r2, blobtype blob, int mode,
                   double x0 = 0, double y0 = 0, double z0 = 0);

void GaussianMask(MultidimArray< double >& mask,
                  double sigma, int mode = INNER_MASK,
                  double x0 = 0, double y0 = 0, double z0 = 0);

void BinaryDWTCircularMask2D(MultidimArray< int >& mask,
                             double radius, int smin, int smax,
                             const std::string& quadrant);

void SeparableSincKaiserMask2D(MultidimArray<double> &mask,
                               double omega, double delta = 0.01,
                               double Deltaw = 1.0 / 12.0);

void mask2D_4neig(MultidimArray< int >& mask, int value = 1, int center = NO_ACTIVATE);

void mask2D_8neig(MultidimArray< int >& mask, int value1 = 1, int value2 = 1,
                  int center = NO_ACTIVATE);

void BinaryDWTSphericalMask2D(MultidimArray< int >& mask,
                              double radius, int smin, int smax,
                              const std::string& quadrant);

void BinaryCylinderMask(MultidimArray< int >& mask,
                        double R, double H, int mode = INNER_MASK,
                        double x0 = 0, double y0 = 0, int z0 = 0);

void BinaryFrameMask(MultidimArray< int >& mask,
                     int Xrect, int Yrect, int Zrect, int mode = INNER_MASK,
                     double x0 = 0, double y0 = 0, double z0 = 0);

void BinaryConeMask(MultidimArray< int >& mask,
                    double theta, int mode = INNER_MASK, bool centerOrigin=false);

void BinaryWedgeMask(MultidimArray< int >& mask, double theta0, double thetaF,
                     const Matrix2D< double > &A, bool centerOrigin=false);


void mask3D_6neig(MultidimArray< int >& mask, int value = 1, int center = NO_ACTIVATE);

void mask3D_18neig(MultidimArray< int >& mask, int value1 = 1, int value2 = 1,
                   int center = NO_ACTIVATE);

void mask3D_26neig(MultidimArray< int >& mask, int value1 = 1, int value2 = 1,
                   int value3 = 1, int center = NO_ACTIVATE);

class Mask
{
public:

#define NO_MASK                   0
#define BINARY_CIRCULAR_MASK      1
#define BINARY_CROWN_MASK         2
#define BINARY_CYLINDER_MASK      3
#define BINARY_FRAME_MASK         4
#define GAUSSIAN_MASK             5
#define RAISED_COSINE_MASK        6
#define BLACKMAN_MASK             7
#define SINC_MASK                 8
#define SINC_BLACKMAN_MASK        9
#define READ_BINARY_MASK         10
#define READ_REAL_MASK           11
#define RAISED_CROWN_MASK        12
#define BINARY_DWT_CIRCULAR_MASK 13
#define BINARY_DWT_SPHERICAL_MASK 14
#define BINARY_CONE_MASK         15
#define BINARY_WEDGE_MASK        16
#define BLOB_CIRCULAR_MASK       17
#define BLOB_CROWN_MASK          18
#define BINARY_TUBE              19

#define INT_MASK    1
#define DOUBLE_MASK 2
#define ALL_KINDS   INT_MASK | DOUBLE_MASK

    static void defineParams(XmippProgram * program, int allowed_data_types = ALL_KINDS,
                             const char* prefix=nullptr, const char* comment=nullptr, bool moreOptions=false);
    void readParams(XmippProgram * program);


    int type;

    int mode;

    double R1;

    double R2;

    double pix_width;

    int blob_order;
    
    double blob_radius;
    double blob_alpha;

    double H;

    double sigma;

    double omega;

    int Xrect;

    int Yrect;

    int Zrect;

    double z0;

    double y0;

    double x0;

    int smin;

    int smax;

    std::string quadrant;

    FileName fn_mask;

    Matrix2D< double > mask_geo;

    int allowed_data_types;

    MultidimArray< int > imask;

    MultidimArray< double > dmask;

    std::string mask_type;


public:

    Mask(int _allowed_data_type = ALL_KINDS);

    void clear();

    void read(int argc, const char** argv);

    void show() const;

    void usage() const;

    void write_mask(const FileName& fn);

    int datatype()
    {
        if (type == BINARY_CIRCULAR_MASK || type == BINARY_CROWN_MASK ||
            type == BINARY_CYLINDER_MASK || type == BINARY_FRAME_MASK ||
            type == BINARY_TUBE ||
            type == NO_MASK || type == READ_BINARY_MASK ||
            type == BINARY_DWT_CIRCULAR_MASK || type == BINARY_CONE_MASK)
            return INT_MASK;

        else if (type == GAUSSIAN_MASK || type == RAISED_COSINE_MASK ||
                 type == SINC_MASK || type == SINC_BLACKMAN_MASK ||
                 type == BLACKMAN_MASK || type == RAISED_CROWN_MASK ||
                 type == BINARY_WEDGE_MASK || type == BLOB_CIRCULAR_MASK ||
                 type == BLOB_CROWN_MASK|| type == READ_REAL_MASK)
            return DOUBLE_MASK;

        return 0;
    }

    void resize(size_t Xdim);

    void resize(size_t Ydim, size_t Xdim);

    void resize(size_t Zdim, size_t Ydim, size_t Xdim);

    template<typename T>
    void resize(const MultidimArray< T>& m)
    {
        switch (datatype())
        {
        case INT_MASK:
            imask.resizeNoCopy(m);
            break;

        case DOUBLE_MASK:
            dmask.resizeNoCopy(m);
            break;
        }
    }

    void generate_mask(bool apply_geo = false );

    void generate_mask(size_t Zdim, size_t Ydim, size_t Xdim)
    {
        resize(Zdim, Ydim, Xdim);
        generate_mask();
    }

    void generate_mask(int Ydim, int Xdim, bool apply_geo = false)
    {
        resize(Ydim, Xdim);
        generate_mask(apply_geo);
    }

    void generate_mask(int Xdim)
    {
        resize(Xdim);
        generate_mask();
    }

    template<typename T>
    void generate_mask(const MultidimArray< T >& m, bool apply_geo = false)
    {
        resize(m);
        generate_mask(apply_geo);
    }

    template<typename T>
    void apply_mask(const MultidimArray< T >& I, MultidimArray< T >& result,
                    T subs_val = 0, bool apply_geo = false)
    {
        switch (datatype())
        {
        case INT_MASK:
            if (apply_geo)
                apply_geo_binary_2D_mask(imask, mask_geo);

            apply_binary_mask(imask, I, result, subs_val);
            break;

        case DOUBLE_MASK:
            if (apply_geo)
                apply_geo_cont_2D_mask(dmask, mask_geo);

            apply_cont_mask(dmask, I, result);
            break;
        }
    }


    template<typename T>
    void produce_vector(const MultidimArray< T >& I, MultidimArray< T >& result)
    {
        // Resize the output vector
        if (XSIZE(result) == 0)
        {
            int size = 0;
            switch (datatype())
            {
            case INT_MASK:
                FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(imask)
                if (DIRECT_A2D_ELEM(imask, i, j) > 0)
                    size++;
                break;

            case DOUBLE_MASK:
                FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(dmask)
                if (DIRECT_A2D_ELEM(dmask, i, j) > 0)
                    size++;
                break;
            }
            result.initZeros(size);
        }

        int p = 0;
        switch (datatype())
        {
        case INT_MASK:
            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(imask)
            if (DIRECT_A2D_ELEM(imask, i, j) > 0)
                DIRECT_A1D_ELEM(result, p++) = DIRECT_A3D_ELEM(I, k, i, j);
            break;
        case DOUBLE_MASK:
            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(dmask)
            if (DIRECT_A2D_ELEM(dmask, i, j) > 0)
                DIRECT_A1D_ELEM(result, p++) = DIRECT_A3D_ELEM(I, k, i, j);
            break;
        }
    }

    const MultidimArray< int >& get_binary_mask() const
    {
        return imask;
    }

    MultidimArray< int >& get_binary_mask()
    {
        return imask;
    }

    void set_binary_mask(MultidimArray< int >& _imask)
    {
        imask = _imask;
    }

    const MultidimArray< double >& get_cont_mask() const
    {
        return dmask;
    }

    MultidimArray< double >& get_cont_mask()
    {
        return dmask;
    }

    void set_cont_mask(MultidimArray< double >& _dmask)
    {
        dmask = _dmask;
    }

    void force_to_be_continuous()
    {
        if (datatype() == INT_MASK)
        {
            typeCast(imask, dmask);
        }
    }

    void force_to_be_binary()
    {
        if (datatype() == DOUBLE_MASK)
        {
            typeCast(dmask, imask);
        }
    }

};

void apply_geo_binary_2D_mask(MultidimArray< int >& mask,
                              const Matrix2D< double >& A);

void apply_geo_cont_2D_mask(MultidimArray< double >& mask,
                            const Matrix2D< double >& A);

template<typename T1, typename T>
void computeStats_within_binary_mask(const MultidimArray< T1 >& mask,
                                     const MultidimArray< T >& m, double& min_val,
                                     double& max_val,
                                     double& avg, double& stddev)
{
    SPEED_UP_tempsInt;
    double sum1 = 0;
    double sum2 = 0;
    int N = 0;

    max_val = min_val = DIRECT_A3D_ELEM(m, 0, 0, 0);

    FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(mask, m)
    {
        if (A3D_ELEM(mask, k, i, j) > 0)
        {
            N++;

            double aux=A3D_ELEM(m, k, i, j);
            // Minimum and maximum
            if (aux < min_val)
                min_val = aux;
            else if (aux > max_val)
                max_val = aux;

            // cumulative sums for average and standard deviation
            sum1 += aux;
            sum2 += aux*aux;
        }
    }

    // average and standard deviation
    avg  = sum1 / (double) N;
    if (N > 1)
        stddev = sqrt(fabs(sum2 / N - avg * avg) * N / (N - 1));
    else
        stddev = 0;
}

inline void computeStats_within_binary_mask(const MultidimArray< int >& mask,
        const MultidimArrayGeneric &m, double& min_val,
        double& max_val,
        double& avg, double& stddev)
{
#define COMPUTESTATS(type) \
computeStats_within_binary_mask(mask,*((MultidimArray<type>*)m.im),min_val,max_val,avg,stddev);

    SWITCHDATATYPE(m.datatype,COMPUTESTATS);

#undef COMPUTESTATS
}

template<typename T>
void apply_binary_mask(const MultidimArray< int >& mask, const MultidimArray< T >& m_in,
                       MultidimArray< T >& m_out, T subs_val = (T) 0)
{
    m_out.resize(m_in);
    FOR_ALL_ELEMENTS_IN_ARRAY3D(m_out)
    // If in common with the mask
    if (k >= STARTINGZ(mask) && k <= FINISHINGZ(mask) &&
        i >= STARTINGY(mask) && i <= FINISHINGY(mask) &&
        j >= STARTINGX(mask) && j <= FINISHINGX(mask))
        if (A3D_ELEM(mask, k, i, j) == 0)
            A3D_ELEM(m_out, k, i, j) = subs_val;
        else
            A3D_ELEM(m_out, k, i, j) = A3D_ELEM(m_in, k, i, j);
    // It is not in common, leave the original one
    else
        A3D_ELEM(m_out, k, i, j) = A3D_ELEM(m_in, k, i, j);
}

template<typename T>
void apply_cont_mask(const MultidimArray< double >& mask, const MultidimArray< T >& m_in,
                     MultidimArray< T >& m_out)
{
    m_out.resize(m_in);
    FOR_ALL_ELEMENTS_IN_ARRAY3D(m_out)
    // If in common with the mask
    if ((k >= STARTINGZ(mask)) && (k <= FINISHINGZ(mask)) &&
        (i >= STARTINGY(mask)) && (i <= FINISHINGY(mask)) &&
        (j >= STARTINGX(mask) && j <= FINISHINGX(mask)))
    {
        A3D_ELEM(m_out, k, i, j) = (T)(A3D_ELEM(m_in, k, i, j) * A3D_ELEM(mask, k, i, j));
    }
    // It is not in common, leave the original one
    else
        A3D_ELEM(m_out, k, i, j) = A3D_ELEM(m_in, k, i, j);
}

template<typename T>
void compute_hist_within_binary_mask(const MultidimArray< int >& mask,
                                     MultidimArray< T >& v, Histogram1D &hist,
                                     int no_steps)
{
    T min_val;
    T max_val;
    double avg;
    double stddev;

    computeStats_within_binary_mask(mask, v, min_val, max_val, avg, stddev);
    compute_hist_within_binary_mask(mask, v, hist, min_val, max_val, no_steps);
}

template<typename T>
void compute_hist_within_binary_mask(const MultidimArray< int >& mask,
                                     const MultidimArray< T >& v, Histogram1D& hist,
                                     T min, T max, int no_steps)
{
    SPEED_UP_tempsInt;
    hist.init(min, max, no_steps);
    FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(mask, v)
    if (A3D_ELEM(mask, k, i, j) != 0)
    {
        double value=A3D_ELEM(v, k, i, j);
        INSERT_VALUE(hist,value);
    }
}

constexpr int  COUNT_ABOVE = 1;
constexpr int  COUNT_BELOW = 2;
constexpr int  COUNT_BETWEEN = 3;

#define count_with_mask_above(mask, m, th) \
    count_with_mask(mask, m, COUNT_ABOVE, th, 0);

#define count_with_mask_below(mask, m, th) \
    count_with_mask(mask, m, COUNT_BELOW, th, 0);

#define count_with_mask_between(mask, m, th1, th2) \
    count_with_mask(mask, m, COUNT_BETWEEN, th1, th2);

template<typename T>
int count_with_mask(const MultidimArray< int >& mask,
                    const MultidimArray< T >& m, int mode, double th1, double th2)
{
    SPEED_UP_tempsInt;
    int N = 0;
    FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(mask, m)
    if (A3D_ELEM(mask, k, i, j))
        switch (mode)
        {
        case COUNT_ABOVE:
                        if (A3D_ELEM(m, k, i, j) >= th1)
                            N++;
            break;

        case COUNT_BELOW:
                        if (A3D_ELEM(m, k, i, j) <= th1)
                            N++;
            break;

        case COUNT_BETWEEN:
                        if (A3D_ELEM(m, k, i, j) >= th1 && A3D_ELEM(m, k, i, j) <= th2)
                            N++;
            break;
        }
    return N;
}

// Specialization for complex numbers
int count_with_mask(const MultidimArray< int >& mask,
                    const MultidimArray< std::complex< double > > & m, int mode,
                    double th1, double th2);

template<typename T>
void invert_binary_mask(MultidimArray< T >& mask)
{
    T* ptr=NULL;
    unsigned long int n;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(mask,n,ptr)
    *ptr = 1-(*ptr);
}

void rangeAdjust_within_mask(const MultidimArray< double >* mask,
                             const MultidimArray< double >& m1,
                             MultidimArray< double >& m2);

class ProgMask: public XmippMetadataProgram
{
public:

    Mask         mask;
    FileName     fn_mask;
    int          create_mask;
    int          count_above;
    double       th_above;
    int          count_below;
    double       th_below;
    double       subs_val;
    std::string  str_subs_val;
    int          count;
    int          max_length;

    void defineParams();
    void readParams();
    void preProcess();
    void postProcess();

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

#endif