filters.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'
***************************************************************************/

#include "core/xmipp_image.h"
#include "core/matrix2d.h"
#include "data/numerical_tools.h"
#include "data/polar.h"

#ifndef CORE_FILTERS_H
#define CORE_FILTERS_H
constexpr double LOG2 = 0.693147181;

class XmippProgram;
class MetaData;


void substractBackgroundPlane(MultidimArray<double> &I);

void substractBackgroundRollingBall(MultidimArray<double> &I, int radius);

void detectBackground(const MultidimArray<double> &vol, MultidimArray<double> &mask, double alpha,
                      double &final_mean);

void contrastEnhancement(Image<double>* I);

void regionGrowing2D(const MultidimArray< double >& I_in,
                     MultidimArray< double >& I_out,
                     int i,
                     int j,
                     float stop_colour = 1,
                     float filling_colour = 1,
                     bool less = true,
                     int neighbourhood = 8);

void regionGrowing3D(const MultidimArray< double >& V_in,
                     MultidimArray< double >& V_out,
                     int k,
                     int i,
                     int j,
                     float stop_colour = 1,
                     float filling_colour = 1,
                     bool less = true);


void regionGrowing3DEqualValue(const MultidimArray<double> &V_in,
                                MultidimArray<int> &V_out, int filling_value);

void distanceTransform(const MultidimArray<int> &in,
                       MultidimArray<int> &out, bool wrap=false);

int labelImage2D(const MultidimArray< double >& I,
                 MultidimArray< double >& label,
                 int neighbourhood = 8);

int labelImage3D(const MultidimArray< double >& V, MultidimArray< double >& label);

void removeSmallComponents(MultidimArray< double >& I,
                           int size,
                           int neighbourhood = 8);

void keepBiggestComponent(MultidimArray< double >& I,
                          double percentage = 0,
                          int neighbourhood = 8);

void fillBinaryObject(MultidimArray< double >&I, int neighbourhood = 8);

void varianceFilter(MultidimArray<double> &I, int kernelSize = 10, bool relative=false);

void noisyZonesFilter(MultidimArray<double> &I, int kernelSize = 10);

double giniCoeff(MultidimArray<double> &I, int varKernelSize = 50);

double OtsuSegmentation(MultidimArray<double> &V);

double EntropySegmentation(MultidimArray<double> &V);

double EntropyOtsuSegmentation(MultidimArray<double> &V, double percentil=0.05,
                               bool binarizeVolume=true);

template <typename T>
double correlation(const MultidimArray< T >& x,
                   const MultidimArray< T >& y,
                   const MultidimArray< int >* mask = nullptr,
                   int l = 0,
                   int m = 0,
                   int q = 0)
{
    SPEED_UP_temps;

    double retval = 0; // returned value
    int i;
    int j;
    int k;
    int ip;
    int jp;
    int kp; // indexes
    int Rows;
    int Cols;
    int Slices; // of the volumes

    // do the computation
    Cols = XSIZE(x);
    Rows = YSIZE(x);
    Slices = ZSIZE(x);

    long N = 0;
    for (k = 0; k < Slices; k++)
    {
        kp = k - q;
        if (kp < 0 || kp >= Slices)
            continue;
        for (i = 0; i < Rows; i++)
        {
            ip = i - l;
            if (ip < 0 || ip >= Rows)
                continue;
            for (j = 0; j < Cols; j++)
            {
                jp = j - m;

                if (jp >= 0 && jp < Cols)
                {
                    if (mask != nullptr)
                        if (!DIRECT_A3D_ELEM((*mask), k, i, j))
                            continue;

                    retval += DIRECT_A3D_ELEM(x, k, i, j) *
                              DIRECT_A3D_ELEM(y, kp, ip, jp);
                    ++N;
                }
            }
        }
    }

    return retval / N;
}

template <typename T>
double fastMaskedCorrelation(const MultidimArray< T >& x,
                             const MultidimArray< T >& y,
                             const MultidimArray< int >& mask)
{
    double retval = 0; // returned value
    long N = 0;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(x)
    {
        if (DIRECT_MULTIDIM_ELEM(mask,n))
        {
            retval += DIRECT_MULTIDIM_ELEM(x, n) * DIRECT_MULTIDIM_ELEM(y, n);
            ++N;
        }
    }
    return retval / N;
}

template <typename T>
T fastCorrelation(const MultidimArray< T >& x,
                       const MultidimArray< T >& y)
{
    T * __restrict__ refX = x.data;
    T * __restrict__ refY = y.data;
    size_t nmax=4*(MULTIDIM_SIZE(x)/4);

    T retval = 0;
    // loop unrolling
    for (size_t n=0; n<nmax; n+=4)
    {
        retval += (*refX++)*(*refY++);
        retval += (*refX++)*(*refY++);
        retval += (*refX++)*(*refY++);
        retval += (*refX++)*(*refY++);
    }

    for (size_t n=nmax; n<MULTIDIM_SIZE(x); ++n)
    {
        retval += DIRECT_MULTIDIM_ELEM(x, n)   * DIRECT_MULTIDIM_ELEM(y, n);
    }

    return retval / MULTIDIM_SIZE(x);
}


void covarianceMatrix(const MultidimArray<double> &I, Matrix2D<double> &C);

template <typename T>
double fastCorrentropy(const MultidimArray<T> &x, const MultidimArray<T> &y,
                       double sigma, const GaussianInterpolator &G)
{
    double retval=0;
    double isigma=1.0/sigma;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(x)
    retval+=G.getValue(isigma*(DIRECT_MULTIDIM_ELEM(x,n)-
                               DIRECT_MULTIDIM_ELEM(y,n)));
    retval/=XSIZE(x);
    return retval;
}

double fastCorrentropy(const MultidimArray<double> &x, const MultidimArray<double> &y,
                       double sigma, const GaussianInterpolator &G, const MultidimArray<int> &mask);

template <typename T>
double correntropy(const MultidimArray<T> &x, const MultidimArray<T> &y,
                   double sigma)
{
    double retval=0;
    double K=-0.5/(sigma*sigma);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(x)
    {
        double diff=DIRECT_MULTIDIM_ELEM(x,n)-DIRECT_MULTIDIM_ELEM(y,n);
        retval+=exp(K*diff*diff);
    }
    retval/=XSIZE(x)*YSIZE(x)*ZSIZE(x);
    return retval;
}

double imedDistance(const MultidimArray<double>& I1, const MultidimArray<double>& I2);

double imedNormalizedDistance(const MultidimArray<double>& I1, const MultidimArray<double>& I2);

double correlationMasked(const MultidimArray<double>& I1, const MultidimArray<double>& I2);

double correlationWeighted(MultidimArray<double>& I1, MultidimArray<double>& I2);

double svdCorrelation(const MultidimArray<double>& I1, const MultidimArray<double>& I2, const MultidimArray< int >* mask = nullptr);

double bestShift(const MultidimArray< double >& I1,
               const MultidimArray< double >& I2,
               double& shiftX,
               double& shiftY,
               CorrelationAux &aux,
               const MultidimArray< int >* mask = nullptr,
               int maxShift=-1);

double bestShift(const MultidimArray<double> &I1, const MultidimArray< std::complex<double> > &FFTI1,
               const MultidimArray<double> &I2,
               double &shiftX, double &shiftY, CorrelationAux &aux,
               const MultidimArray<int> *mask=nullptr, int maxShift=-1);

double bestShift(const MultidimArray< std::complex<double> > &FFTI1,
                const MultidimArray< std::complex<double> > &FFTI2,
                MultidimArray<double> &Mcorr,
               double &shiftX, double &shiftY, CorrelationAux &aux,
               const MultidimArray<int> *mask=nullptr, int maxShift=-1);

void bestShift(const MultidimArray<double> &I1, const MultidimArray<double> &I2,
               double &shiftX, double &shiftY, double &shiftZ, CorrelationAux &aux,
               const MultidimArray<int> *mask=nullptr);

template<typename T>
T bestShift(MultidimArray<T> &Mcorr,
               T &shiftX, T &shiftY, const MultidimArray<int> *mask, int maxShift);

void bestNonwrappingShift(const MultidimArray< double >& I1,
                          const MultidimArray< double >& I2,
                          double& shiftX,
                          double& shiftY,
                          CorrelationAux &aux);

void bestNonwrappingShift(const MultidimArray<double> &I1, const MultidimArray< std::complex<double> > &FFTI1,
                          const MultidimArray<double> &I2, double &shiftX, double &shiftY,
                          CorrelationAux &aux);

double bestShiftRealSpace(const MultidimArray<double> &I1, MultidimArray<double> &I2,
               double &shiftX, double &shiftY,
               const MultidimArray<int> *mask=nullptr, int maxShift=5, double shiftStep=1.0);

class AlignmentAux
{
public:
    Matrix2D<double> ARS;
    Matrix2D<double> ASR;
    Matrix2D<double> R;
    MultidimArray<double> IauxSR;
    MultidimArray<double> IauxRS;
    MultidimArray<double> rotationalCorr;
    Polar_fftw_plans *plans;
    Polar< std::complex<double> > polarFourierI;
    AlignmentAux();
    AlignmentAux(const AlignmentAux &)=delete; // Do not use the default copy constructor
    AlignmentAux & operator=(const AlignmentAux&) = delete; // Do not use the default copy assignment
    ~AlignmentAux();
};

class AlignmentTransforms
{
public:
    Polar< std::complex<double> > polarFourierI;
    MultidimArray< std::complex< double > > FFTI;
};

double alignImages(const MultidimArray< double >& Iref,
                   MultidimArray< double >& I,
                   Matrix2D< double >&M,
                   bool wrap=xmipp_transformation::WRAP);

double alignImagesConsideringMirrors(const MultidimArray<double>& Iref, MultidimArray<double>& I,
                   Matrix2D<double>&M, bool wrap);

double alignImagesConsideringMirrors(const MultidimArray<double>& Iref, const AlignmentTransforms& IrefTransforms,
                                     MultidimArray<double>& I, Matrix2D<double> &M, AlignmentAux& aux,
                                     CorrelationAux& aux2, RotationalCorrelationAux &aux3, bool wrap,
                                     const MultidimArray<int>* mask=nullptr);

double alignImages(const MultidimArray< double >& Iref,
                   MultidimArray< double >& I,
                   Matrix2D< double >&M,
                   bool wrap,
                   AlignmentAux &aux,
                   CorrelationAux &aux2,
                   RotationalCorrelationAux &aux3);

class VolumeAlignmentAux
{
public:
    MultidimArray<double> IrefCyl;
    MultidimArray<double> Icyl;
    MultidimArray<double> corr;
    MultidimArray<double> I1;
    MultidimArray<double> I12;
    MultidimArray<double> I123;
    Matrix2D<double> R1;
    Matrix2D<double> R2;
    Matrix2D<double> R3;
};

double bestRotationAroundZ(const MultidimArray< double >& Iref,
                           const MultidimArray< double >& I,
                           CorrelationAux &aux,
                           VolumeAlignmentAux &aux2);

void fastBestRotation(const MultidimArray<double>& IrefCylZ,
                      const MultidimArray<double>& IrefCylY,
                      const MultidimArray<double>& IrefCylX,
                      MultidimArray<double>& I,
                      const String &eulerAngles,
                      Matrix2D<double> &R, CorrelationAux &aux, VolumeAlignmentAux &aux2);

double fastBestRotationAroundZ(const MultidimArray<double>& IrefCylZ,
                               const MultidimArray<double>& I, CorrelationAux &aux,
                               VolumeAlignmentAux &aux2);

double fastBestRotationAroundY(const MultidimArray<double>& IrefCylY,
                               const MultidimArray<double>& I, CorrelationAux &aux,
                               VolumeAlignmentAux &aux2);

double fastBestRotationAroundX(const MultidimArray<double>& IrefCylX,
                               const MultidimArray<double>& I, CorrelationAux &aux,
                               VolumeAlignmentAux &aux2);

double alignImagesConsideringMirrors(const MultidimArray< double >& Iref,
                                     MultidimArray< double >& I,
                                     Matrix2D<double> &M,
                                     AlignmentAux &aux,
                                     CorrelationAux &aux2,
                                     RotationalCorrelationAux &aux3,
                                     bool wrap=xmipp_transformation::WRAP,
                                     const MultidimArray< int >* mask = NULL);

void alignSetOfImages(MetaData &MD, MultidimArray< double >& Iavg,
                      int Niter=10, bool considerMirror=true);

double unnormalizedGaussian2D(const Matrix1D<double> &r,
                              const Matrix1D<double> &mu,
                              const Matrix2D<double> &sigmainv);

void estimateGaussian2D(const MultidimArray<double> &I,
                        double &a, double &b, Matrix1D<double> &mu, Matrix2D<double> &sigma,
                        bool estimateMu=true, int iterations=10);

template <typename T>
double euclidianDistance(const MultidimArray< T >& x,
                         const MultidimArray< T >& y,
                         const MultidimArray< int >* mask = nullptr)
{
    SPEED_UP_temps;

    double retval = 0;
    long n = 0;

    FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(x, y)
    {
        if (mask != nullptr)
            if (!(*mask)(k, i, j))
                continue;

        retval += (A3D_ELEM(x, k, i, j) - A3D_ELEM(y, k, i, j)) *
                  (A3D_ELEM(x, k, i, j) - A3D_ELEM(y, k, i, j));
        n++;
    }

    if (n != 0)
        return sqrt(retval);
    else
        return 0;
}

template <typename T>
double mutualInformation(const MultidimArray< T >& x,
                         const MultidimArray< T >& y,
                         int nx = 0,
                         int ny = 0,
                         const MultidimArray< int >* mask = nullptr);

template <typename T>
double rms(const MultidimArray< T >& x,
           const MultidimArray< T >& y,
           const MultidimArray< int >* mask = nullptr,
           MultidimArray< double >* Contributions = nullptr)
{
    SPEED_UP_tempsInt;

    double retval = 0;
    double aux;
    int n = 0;

    // If contributions are desired
    if (Contributions != nullptr)
    {
        FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(x, y)
        {
            if (mask != nullptr)
                if (!(*mask)(k, i, j))
                    continue;

            aux = (A3D_ELEM(x, k, i, j) - A3D_ELEM(y, k, i, j)) *
                  (A3D_ELEM(x, k, i, j) - A3D_ELEM(y, k, i, j));
            A3D_ELEM(*Contributions, k, i, j) = aux;
            retval += aux;
            n++;
        }

        FOR_ALL_ELEMENTS_IN_ARRAY3D(*Contributions)
        A3D_ELEM(*Contributions, k, i, j) = sqrt(A3D_ELEM(*Contributions,
                                            k, i, j) / n);
    }
    else
    {
        FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(x, y)
        {
            if (mask != nullptr)
                if (!(*mask)(k, i, j))
                    continue;

            retval += (A3D_ELEM(x, k, i, j) - A3D_ELEM(y, k, i, j)) *
                      (A3D_ELEM(x, k, i, j) - A3D_ELEM(y, k, i, j));
            n++;
        }
    }

    if (n != 0)
        return sqrt(retval / n);
    else
        return 0;
}

void fourierBesselDecomposition(const MultidimArray< double >& img_in,
                                double r2,
                                int k1,
                                int k2);


//void harmonicDecomposition(const MultidimArray< double >& img_in,
//                           MultidimArray< double >& v_out);

// Function needed by median filtering
template <typename T>
void sort(T a, T b, T c, MultidimArray< T >& v)
{
    if (a < b)
        if (b < c)
        {
            DIRECT_MULTIDIM_ELEM(v,0) = a;
            DIRECT_MULTIDIM_ELEM(v,1) = b;
            DIRECT_MULTIDIM_ELEM(v,2) = c;
        }
        else if (a < c)
        {
            DIRECT_MULTIDIM_ELEM(v,0) = a;
            DIRECT_MULTIDIM_ELEM(v,1) = c;
            DIRECT_MULTIDIM_ELEM(v,2) = b;
        }
        else
        {
            DIRECT_MULTIDIM_ELEM(v,0) = c;
            DIRECT_MULTIDIM_ELEM(v,1) = a;
            DIRECT_MULTIDIM_ELEM(v,2) = b;
        }
    else if (a < c)
    {
        DIRECT_MULTIDIM_ELEM(v,0) = b;
        DIRECT_MULTIDIM_ELEM(v,1) = a;
        DIRECT_MULTIDIM_ELEM(v,2) = c;
    }
    else if (b < c)
    {
        DIRECT_MULTIDIM_ELEM(v,0) = b;
        DIRECT_MULTIDIM_ELEM(v,1) = c;
        DIRECT_MULTIDIM_ELEM(v,2) = a;
    }
    else
    {
        DIRECT_MULTIDIM_ELEM(v,0) = c;
        DIRECT_MULTIDIM_ELEM(v,1) = b;
        DIRECT_MULTIDIM_ELEM(v,2) = a;
    }
}

// Function needed by median filtering
template <typename T>
void mergeSort(MultidimArray< T >& v1, MultidimArray< T >& v2, MultidimArray< T >& v)
{
    int i1 = 0;
    int i2 = 0;
    int i = 0;

    while ((i1 < 3) && (i2 < 3))
    {
        if (v1(i1) < v2(i2))
            v(i++) = v1(i1++);
        else
            v(i++) = v2(i2++);
    }

    while (i1 < 3)
        v(i++) = v1(i1++);

    while (i2 < 3)
        v(i++) = v2(i2++);
}

// Function needed by median filtering
// This UGLY function performs a fast merge sort for the case of vectors of 3
// elements. This way is guaranteed a minimum number of comparisons (maximum
// number of comparisons to perform the sort, 5)
template <typename T>
void fastMergeSort(MultidimArray< T >& x, MultidimArray< T >& y, MultidimArray< T >& v)
{
    if (DIRECT_MULTIDIM_ELEM(x,0) < DIRECT_MULTIDIM_ELEM(y,0))
    {
        DIRECT_MULTIDIM_ELEM(v,0) = DIRECT_MULTIDIM_ELEM(x,0);
        if (DIRECT_MULTIDIM_ELEM(x,1) < DIRECT_MULTIDIM_ELEM(y,0))
        {
            DIRECT_MULTIDIM_ELEM(v,1) = DIRECT_MULTIDIM_ELEM(x,1);
            if (DIRECT_MULTIDIM_ELEM(x,2) < DIRECT_MULTIDIM_ELEM(y,0))
            {
                DIRECT_MULTIDIM_ELEM(v,2) = DIRECT_MULTIDIM_ELEM(x,2);
                DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(y,0);
                DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(y,1);
                DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(y,2);
            }
            else
            {
                DIRECT_MULTIDIM_ELEM(v,2) = DIRECT_MULTIDIM_ELEM(y,0);
                if (DIRECT_MULTIDIM_ELEM(x,2) < DIRECT_MULTIDIM_ELEM(y,1))
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(x,2);
                    DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(y,1);
                    DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(y,2);
                }
                else
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(y,1);
                    if (DIRECT_MULTIDIM_ELEM(x,2) < DIRECT_MULTIDIM_ELEM(y,2))
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(x,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(y,2);
                    }
                    else
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(y,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(x,2);
                    }
                }
            }
        }
        else
        {
            DIRECT_MULTIDIM_ELEM(v,1) = DIRECT_MULTIDIM_ELEM(y,0);
            if (DIRECT_MULTIDIM_ELEM(x,1) < DIRECT_MULTIDIM_ELEM(y,1))
            {
                DIRECT_MULTIDIM_ELEM(v,2) = DIRECT_MULTIDIM_ELEM(x,1);
                if (DIRECT_MULTIDIM_ELEM(x,2) < DIRECT_MULTIDIM_ELEM(y,1))
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(x,2);
                    DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(y,1);
                    DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(y,2);
                }
                else
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(y,1);
                    if (DIRECT_MULTIDIM_ELEM(x,2) < DIRECT_MULTIDIM_ELEM(y,2))
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(x,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(y,2);
                    }
                    else
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(y,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(x,2);
                    }
                }
            }
            else
            {
                DIRECT_MULTIDIM_ELEM(v,2) = DIRECT_MULTIDIM_ELEM(y,1);
                if (DIRECT_MULTIDIM_ELEM(x,1) < DIRECT_MULTIDIM_ELEM(y,2))
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(x,1);
                    if (DIRECT_MULTIDIM_ELEM(x,2) < DIRECT_MULTIDIM_ELEM(y,2))
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(x,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(y,2);
                    }
                    else
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(y,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(x,2);
                    }
                }
                else
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(y,2);
                    DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(x,1);
                    DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(x,2);
                }
            }
        }
    }
    else
    {
        DIRECT_MULTIDIM_ELEM(v,0) = DIRECT_MULTIDIM_ELEM(y,0);
        if (DIRECT_MULTIDIM_ELEM(x,0) < DIRECT_MULTIDIM_ELEM(y,1))
        {
            DIRECT_MULTIDIM_ELEM(v,1) = DIRECT_MULTIDIM_ELEM(x,0);
            if (DIRECT_MULTIDIM_ELEM(x,1) < DIRECT_MULTIDIM_ELEM(y,1))
            {
                DIRECT_MULTIDIM_ELEM(v,2) = DIRECT_MULTIDIM_ELEM(x,1);
                if (DIRECT_MULTIDIM_ELEM(x,2) < DIRECT_MULTIDIM_ELEM(y,1))
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(x,2);
                    DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(y,1);
                    DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(y,2);
                }
                else
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(y,1);
                    if (DIRECT_MULTIDIM_ELEM(x,2) < DIRECT_MULTIDIM_ELEM(y,2))
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(x,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(y,2);
                    }
                    else
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(y,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(x,2);
                    }
                }
            }
            else
            {
                DIRECT_MULTIDIM_ELEM(v,2) = DIRECT_MULTIDIM_ELEM(y,1);
                if (DIRECT_MULTIDIM_ELEM(x,1) < DIRECT_MULTIDIM_ELEM(y,2))
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(x,1);
                    if (DIRECT_MULTIDIM_ELEM(x,2) < DIRECT_MULTIDIM_ELEM(y,2))
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(x,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(y,2);
                    }
                    else
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(y,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(x,2);
                    }
                }
                else
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(y,2);
                    DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(x,1);
                    DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(x,2);
                }
            }
        }
        else
        {
            DIRECT_MULTIDIM_ELEM(v,1) = DIRECT_MULTIDIM_ELEM(y,1);
            if (DIRECT_MULTIDIM_ELEM(x,0) < DIRECT_MULTIDIM_ELEM(y,2))
            {
                DIRECT_MULTIDIM_ELEM(v,2) = DIRECT_MULTIDIM_ELEM(x,0);
                if (DIRECT_MULTIDIM_ELEM(x,1) < DIRECT_MULTIDIM_ELEM(y,2))
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(x,1);
                    if (DIRECT_MULTIDIM_ELEM(x,2) < DIRECT_MULTIDIM_ELEM(y,2))
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(x,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(y,2);
                    }
                    else
                    {
                        DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(y,2);
                        DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(x,2);
                    }
                }
                else
                {
                    DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(y,2);
                    DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(x,1);
                    DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(x,2);
                }
            }
            else
            {
                DIRECT_MULTIDIM_ELEM(v,2) = DIRECT_MULTIDIM_ELEM(y,2);
                DIRECT_MULTIDIM_ELEM(v,3) = DIRECT_MULTIDIM_ELEM(x,0);
                DIRECT_MULTIDIM_ELEM(v,4) = DIRECT_MULTIDIM_ELEM(x,1);
                DIRECT_MULTIDIM_ELEM(v,5) = DIRECT_MULTIDIM_ELEM(x,2);
            }
        }
    }
}

// Function needed by median filtering
template <typename T>
void median(MultidimArray< T >& x, MultidimArray< T >& y, T& m)
{
    if (DIRECT_MULTIDIM_ELEM(x,0) < DIRECT_MULTIDIM_ELEM(y,1))
        if (DIRECT_MULTIDIM_ELEM(x,1) < DIRECT_MULTIDIM_ELEM(y,1))
            if (DIRECT_MULTIDIM_ELEM(x,2) < DIRECT_MULTIDIM_ELEM(y,1))
                m = DIRECT_MULTIDIM_ELEM(y,1);
            else
                m = DIRECT_MULTIDIM_ELEM(x,2);
        else
            if (DIRECT_MULTIDIM_ELEM(x,1) < DIRECT_MULTIDIM_ELEM(y,2))
                m = DIRECT_MULTIDIM_ELEM(y,2);
            else
                m = DIRECT_MULTIDIM_ELEM(x,1);
    else
        if (DIRECT_MULTIDIM_ELEM(x,0) < DIRECT_MULTIDIM_ELEM(y,2))
            if (DIRECT_MULTIDIM_ELEM(x,1) < DIRECT_MULTIDIM_ELEM(y,2))
                m = DIRECT_MULTIDIM_ELEM(y,2);
            else
                m = DIRECT_MULTIDIM_ELEM(x,1);
        else
            if (DIRECT_MULTIDIM_ELEM(x,0) < DIRECT_MULTIDIM_ELEM(y,3))
                m = DIRECT_MULTIDIM_ELEM(y,3);
            else
                if (DIRECT_MULTIDIM_ELEM(x,0) < DIRECT_MULTIDIM_ELEM(y,4))
                    m = DIRECT_MULTIDIM_ELEM(x,0);
                else
                    m = DIRECT_MULTIDIM_ELEM(y,4);
}

template <typename T>
void medianFilter3x3(MultidimArray< T >&m, MultidimArray< T >& out)
{
    int backup_startingx = STARTINGX(m);
    int backup_startingy = STARTINGY(m);

    STARTINGX(m) = STARTINGY(m) = 0;
    MultidimArray< T > v1(3);
    MultidimArray< T > v2(3);
    MultidimArray< T > v3(3);
    MultidimArray< T > v4(3);
    MultidimArray< T > v(6);

    // Set the output matrix size
    out.initZeros(m);

    // Set the initial and final matrix indices to explore
    int initialY = 1;
    int initialX = 1;
    int finalY = YSIZE(m) - 2;
    int finalX = XSIZE(m) - 2;

    // For every row
    for (int i = initialY; i <= finalY; i++)
    {
        // For every pair of pixels (mean is computed obtaining
        // two means at the same time using an efficient method)
        for (int j = initialX; j <= finalX; j += 2)
        {
            // If we are in the border case
            if (j == 1)
            {
                // Order the first and second vectors of 3 elements
                sort(DIRECT_A2D_ELEM(m, i - 1, j - 1), DIRECT_A2D_ELEM(m, i, j - 1),
                     DIRECT_A2D_ELEM(m, i + 1, j - 1), v1);
                sort(DIRECT_A2D_ELEM(m, i - 1, j), DIRECT_A2D_ELEM(m, i, j),
                     DIRECT_A2D_ELEM(m, i + 1, j), v2);
            }
            else
            {
                // Simply take ordered vectors from previous
                v1 = v3;
                v2 = v4;
            }

            // As we are computing 2 medians at the same time, if the matrix has
            // an odd number of columns, the last column isn't calculated. It is
            // done here
            if (j == finalX)
            {
                v1 = v3;
                v2 = v4;
                sort(DIRECT_A2D_ELEM(m, i - 1, j + 1), DIRECT_A2D_ELEM(m, i, j + 1),
                     DIRECT_A2D_ELEM(m, i + 1, j + 1), v3);
                fastMergeSort(v2, v3, v);
                median(v1, v, DIRECT_A2D_ELEM(out,i, j));
            }
            else
            {
                // Order the third and fourth vectors of 3 elements
                sort(DIRECT_A2D_ELEM(m, i - 1, j + 1), DIRECT_A2D_ELEM(m, i, j + 1),
                     DIRECT_A2D_ELEM(m, i + 1, j + 1), v3);
                sort(DIRECT_A2D_ELEM(m, i - 1, j + 2), DIRECT_A2D_ELEM(m, i, j + 2),
                     DIRECT_A2D_ELEM(m, i + 1, j + 2), v4);

                // Merge sort the second and third vectors
                fastMergeSort(v2, v3, v);

                // Find the first median and assign it to the output
                median(v1, v, DIRECT_A2D_ELEM(out, i, j));

                // Find the second median and assign it to the output
                median(v4, v, DIRECT_A2D_ELEM(out, i, j + 1));
            }
        }
    }

    STARTINGX(m) = STARTINGX(out) = backup_startingx;
    STARTINGY(m) = STARTINGY(out) = backup_startingy;
}

void smoothingShah(MultidimArray< double >& img,
                   MultidimArray< double >& surface_strength,
                   MultidimArray< double >& edge_strength,
                   const Matrix1D< double >& W,
                   int OuterLoops,
                   int InnerLoops,
                   int RefinementLoops,
                   bool adjust_range = true);

double tomographicDiffusion(MultidimArray< double >& V,
                            const Matrix1D< double >& alpha, double lambda);

void rotationalInvariantMoments(const MultidimArray< double >& img,
                                const MultidimArray< int >* mask,
                                MultidimArray< double >& v_out);

void inertiaMoments(const MultidimArray< double >& img,
                    const MultidimArray< int >* mask,
                    Matrix1D< double >& v_out,
                    Matrix2D< double >& u);

void fillTriangle(MultidimArray< double >&img, int* tx, int* ty, double color);

void localThresholding(MultidimArray< double >& img,
                       double C,
                       double dimLocal,
                       MultidimArray< int >& result,
                       MultidimArray< int >* mask = nullptr);

void centerImageTranslationally(MultidimArray<double> &I,
                                CorrelationAux &aux);

void centerImageRotationally(MultidimArray<double> &I, RotationalCorrelationAux &aux);

Matrix2D<double> centerImage(MultidimArray<double> &I, CorrelationAux &aux,
                 RotationalCorrelationAux &aux2,
                 int Niter=10, bool limitShift=true);


void forcePositive(MultidimArray<double> &V);


template <typename T>
void boundMedianFilter(MultidimArray< T > &V, const MultidimArray<char> &mask, int n=0)
{
    bool badRemaining;
    T neighbours[125];
    T aux;
    int N = 0;
    int index;

    do
    {
        badRemaining=false;

        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(V)
        if (DIRECT_A3D_ELEM(mask, k, i, j) != 0)
        {
            N = 0;
            for (int kk=-2; kk<=2; kk++)
            {
                size_t kkk=k+kk;
                if (kkk<0 || kkk>=ZSIZE(V))
                    continue;
                for (int ii=-2; ii<=2; ii++)
                {
                    size_t iii=i+ii;
                    if (iii<0 || iii>=YSIZE(V))
                        continue;
                    for (int jj=-2; jj<=2; jj++)
                    {
                        size_t jjj=j+jj;
                        if (jjj<0 || jjj>=XSIZE(V))
                            continue;

                        if (DIRECT_A3D_ELEM(mask, kkk, iii, jjj) == 0)
                        {
                            index = N++;
                            neighbours[index] = DIRECT_A3D_ELEM(V, kkk,iii,jjj);
                            //insertion sort
                            while (index > 0 && neighbours[index-1] > neighbours[index])
                            {
                                SWAP(neighbours[index-1], neighbours[index], aux);
                                --index;
                            }
                        }
                    }
                }
                if (N == 0)
                    badRemaining = true;
                else
                {
                    //std::sort(neighbours.begin(),neighbours.end());
                    if (N % 2 == 0)
                        DIRECT_A3D_ELEM(V, k, i, j) = (T)(0.5*(neighbours[N/2-1]+ neighbours[N/2]));
                    else
                        DIRECT_A3D_ELEM(V, k, i, j) = neighbours[N/2];
                    DIRECT_A3D_ELEM(mask, k, i, j) = false;
                }
            }
        }
    }
    while (badRemaining);
}

template <typename T>
void pixelDesvFilter(MultidimArray< T > &V, double thresFactor)
{
    if (thresFactor > 0 )
    {
        double avg;
        double stddev;
        double high;
        double low;
        T dummy;
        MultidimArray<char> mask(ZSIZE(V), YSIZE(V), XSIZE(V));
        avg = stddev = low = high = 0;
        V.computeStats(avg, stddev, dummy, dummy);//min and max not used
        low  = (avg - thresFactor * stddev);
        high = (avg + thresFactor * stddev);

        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(V)
        {
            double x = DIRECT_MULTIDIM_ELEM(V, n);
            DIRECT_MULTIDIM_ELEM(mask, n) = (x < low || x > high) ? 1 : 0;
        }
        boundMedianFilter(V, mask);
    }
}

//#include <math.h>

template <typename T>
void logFilter(MultidimArray< T > &V, double a, double b, double c)
{

    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(V)
    {
        double x = DIRECT_MULTIDIM_ELEM(V, n);
        //this casting is kind of risky
        DIRECT_MULTIDIM_ELEM(V, n) = (T)(a-b*log(x+c));
    }
}

void denoiseTVFilter(MultidimArray<double> &V, int maxIter);


void computeEdges (const MultidimArray <double>& vol, MultidimArray<double> &vol_edge);

void forceDWTSparsity(MultidimArray<double> &V, double eps);

class XmippFilter
{
public:
    virtual void readParams(XmippProgram *program)
    {}//do nothing by default
    virtual void apply(MultidimArray<double> &img) = 0;

    virtual ~XmippFilter()
    {}
    ;

    virtual void show()
    {}
    ;
};

class BadPixelFilter: public XmippFilter
{
public:
    typedef enum { NEGATIVE, MASK, OUTLIER } BadPixelFilterType;

    BadPixelFilterType type; //type of filter
    double factor;    //for the case of outliers bad pixels
    Image<char> *mask; //for the case of mask bad pixels

    static void defineParams(XmippProgram * program);
    void readParams(XmippProgram * program);
    void apply(MultidimArray<double> &img);
};

class LogFilter: public XmippFilter
{
public:
    //4.431-0.4018*LN(ABS(P1+336.6))
    //a-b*ln(x+c)
    double a;
    double b;
    double c;
    static void defineParams(XmippProgram * program);
    void readParams(XmippProgram * program);
    void apply(MultidimArray<double> &img);
};

class DenoiseTVFilter: public XmippFilter
{
public:
    int maxIter;
    static void defineParams(XmippProgram * program);
    void readParams(XmippProgram * program);
    void apply(MultidimArray<double> &img);
};

class BackgroundFilter: public XmippFilter
{
    typedef enum { PLANE, ROLLINGBALL } BackgroundType;
    BackgroundType type;
    int radius;

public:
    static void defineParams(XmippProgram * program);
    void readParams(XmippProgram * program);
    void apply(MultidimArray<double> &img);
};

class MedianFilter: public XmippFilter
{
public:
    static void defineParams(XmippProgram * program);
    void readParams(XmippProgram * program);
    void apply(MultidimArray<double> &img);
};

class DiffusionFilter: public XmippFilter
{
    int Shah_outer;

    int Shah_inner;

    int Shah_refinement;

    Matrix1D< double > Shah_weight;

    bool Shah_edge;
public:
    static void defineParams(XmippProgram * program);
    void readParams(XmippProgram * program);
    void show();
    void apply(MultidimArray<double> &img);
};

class BasisFilter: public XmippFilter
{
    FileName fnBasis;

    int Nbasis;

    // Stack with the basis
    Image<double> basis;
public:
    static void defineParams(XmippProgram * program);
    void readParams(XmippProgram * program);
    void show();
    void apply(MultidimArray<double> &img);
};

class RetinexFilter: public XmippFilter
{
public:
    double percentile; // Percentile of noise derivative energies filtered out
    double eps; 
    Image<int> *mask; // for the case of mask bad pixels

    static void defineParams(XmippProgram * program);
    void readParams(XmippProgram * program);
    void apply(MultidimArray<double> &img);
    void laplacian(const MultidimArray<double> &img, MultidimArray< std::complex<double> > &fimg, bool direct);
};

/* Matlab 1D filter function, y=filter(b,a,x);
 * a and b may be modified either by dividing by a0, or by resizing it.
 * Z is an auxiliary vector that will help to calculate the filter
 * */
void matlab_filter(Matrix1D<double> &B, Matrix1D<double> &A,
            const MultidimArray<double> &X, MultidimArray<double> &Y, MultidimArray<double> &Z);

/* Matlab 1D filtfilt function, y=filtfilt(b,a,x);
 * a(0) is supposed to be 1
 * */

#endif