Filters

Collaboration diagram for Filters:

Functions
double	denoiseTVenergy (double mu, const MultidimArray< double > &X, const MultidimArray< double > &Y, double s, double q, double lambda, double sigmag, double g)
void	denoiseTVgradient (double mu, const MultidimArray< double > &X, const MultidimArray< double > &Y, double s, double q, double lambda, double sigmag, double g, MultidimArray< double > &gradientI)
void	denoiseTVproj (const MultidimArray< double > &X, const MultidimArray< double > &Y, double theta, MultidimArray< double > &dold)
void	denoiseTVFilter (MultidimArray< double > &xnew, int maxIter)
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	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)
template<typename T >
double	fastMaskedCorrelation (const MultidimArray< T > &x, const MultidimArray< T > &y, const MultidimArray< int > &mask)
template<typename T >
double	fastCorrentropy (const MultidimArray< T > &x, const MultidimArray< T > &y, double sigma, const GaussianInterpolator &G)
template<typename T >
double	correntropy (const MultidimArray< T > &x, const MultidimArray< T > &y, double sigma)
double	bestShift (const MultidimArray< double > &I1, const MultidimArray< double > &I2, 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)
void	bestNonwrappingShift (const MultidimArray< double > &I1, 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)
double	alignImages (const MultidimArray< double > &Iref, MultidimArray< double > &I, Matrix2D< double > &M, bool wrap=xmipp_transformation::WRAP)
double	alignImages (const MultidimArray< double > &Iref, MultidimArray< double > &I, Matrix2D< double > &M, bool wrap, AlignmentAux &aux, CorrelationAux &aux2, RotationalCorrelationAux &aux3)
double	bestRotationAroundZ (const MultidimArray< double > &Iref, 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	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)
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)
void	fourierBesselDecomposition (const MultidimArray< double > &img_in, double r2, int k1, int k2)
template<typename T >
void	sort (T a, T b, T c, MultidimArray< T > &v)
template<typename T >
void	medianFilter3x3 (MultidimArray< T > &m, MultidimArray< T > &out)
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)
template<typename T >
void	pixelDesvFilter (MultidimArray< T > &V, double thresFactor)
template<typename T >
void	logFilter (MultidimArray< T > &V, double a, double b, double c)
void	SmoothResize (byte *srcpic8, byte *destpic8, size_t swide, size_t shigh, size_t dwide, size_t dhigh)

Functions for all multidimensional arrays
template<typename T >
double	correlationIndex (const MultidimArray< T > &x, const MultidimArray< T > &y, const MultidimArray< int > *mask=NULL, MultidimArray< double > *Contributions=NULL)

Detailed Description

Function Documentation

alignImages() [1/2]

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

Align two images

This function modifies I to be aligned with Iref. Translational and rotational alignments are both considered. The matrix transforming I into Iref is returned.

The function returns the correlation between the two aligned images.

Definition at line 2141 of file filters.cpp.

{
    AlignmentAux aux;
    CorrelationAux aux2;
    RotationalCorrelationAux aux3;
    return alignImages(Iref, I, M, wrap, aux, aux2, aux3);
}

alignImages() [2/2]

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

Fast version of align two images

Definition at line 2131 of file filters.cpp.

{
    Iref.checkDimension(2);
    AlignmentTransforms IrefTransforms;
    computeAlignmentTransforms(Iref,IrefTransforms, aux, aux2);
    return alignImages(Iref, IrefTransforms, I, M, wrap, aux, aux2, aux3);
}

alignImagesConsideringMirrors()

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
	)

Align two images considering also the mirrors

This function modifies I to be aligned with Iref. Translational and rotational alignments are both considered. The correlation coefficient between I transformed and Iref is returned. A mask can be supplied for computing this correlation.

Definition at line 2189 of file filters.cpp.

{
    AlignmentTransforms IrefTransforms;
    computeAlignmentTransforms(Iref,IrefTransforms, aux, aux2);
    return alignImagesConsideringMirrors(Iref, IrefTransforms, I, M, aux, aux2, aux3, wrap, mask);
}

bestNonwrappingShift()

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

Translational search (non-wrapping)

This function returns the best interpolated shift for the alignment of two images. You can restrict the shift to a region defined by a mask (the maximum will be sought where the mask is 1).

Definition at line 1895 of file filters.cpp.

{
    aux.transformer1.FourierTransform((MultidimArray<double> &)I1, aux.FFT1, false);
    bestNonwrappingShift(I1,aux.FFT1,I2,shiftX,shiftY,aux);
}

bestRotationAroundZ()

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

Align two volumes by applying a rotation around Z. The rotation must be applied on I to get Iref.

Definition at line 2356 of file filters.cpp.

{
    Iref.checkDimension(3);
    I.checkDimension(3);
    if (!I.sameShape(Iref))
        REPORT_ERROR(ERR_MULTIDIM_SIZE,
                     "Both volumes should be of the same shape");

    double deltaAng = atan(2.0 / XSIZE(I));
    Matrix1D<double> v(3);
    XX(v) = 0;
    YY(v) = 0;
    ZZ(v) = 1;
    volume_convertCartesianToCylindrical(Iref, aux2.IrefCyl, 3, XSIZE(I) / 2, 1,
                                         0, 2 * PI, deltaAng, v);
    return fastBestRotationAroundZ(aux2.IrefCyl, I, aux, aux2);
}

bestShift() [1/2]

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

Translational search

This function returns the best interpolated shift for the alignment of two images. You can restrict the shift to a region defined by a mask (the maximum will be sought where the mask is 1).

To apply these results you must shift I1 by (-shiftX,-shiftY) or I2 by (shiftX, shiftY).

You can limit the maximum achievable shift by using maxShift. If it is set to -1, any shift is valid.

The function returns the maximum correlation found.

Definition at line 1735 of file filters.cpp.

{
    aux.transformer1.FourierTransform((MultidimArray<double> &)I1, aux.FFT1, false);
    return bestShift(I1,aux.FFT1,I2,shiftX,shiftY,aux,mask,maxShift);
}

bestShift() [2/2]

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

Translational search (3D)

This function returns the best interpolated shift for the alignment of two volumes. You can restrict the shift to a region defined by a mask (the maximum will be sought where the mask is 1).

To apply these results you must shift I1 by (-shiftX,-shiftY,-shiftZ) or I2 by (shiftX, shiftY,shiftZ).

Definition at line 1755 of file filters.cpp.

{
    I1.checkDimension(3);
    I2.checkDimension(3);

    int imax;
    int jmax;
    int kmax;
    int i_actual;
    int j_actual;
    int k_actual;
    double max;
    double xmax;
    double ymax;
    double zmax;
    double sumcorr;
    double avecorr;
    double stdcorr;
    double dummy;
    bool neighbourhood = true;
    MultidimArray<double> Mcorr;

    correlation_matrix(I1, I2, Mcorr, aux);

    /*
     Warning: for masks with a small number of non-zero pixels, this routine is NOT reliable...
     Anyway, maybe using a mask is not a good idea at al...
     */

    // Adjust statistics within shiftmask to average 0 and stddev 1
    if (mask != nullptr)
    {
        if ((*mask).sum() < 2)
        {
            shiftX = shiftY = shiftZ = 0.;
            return;
        }
        else
        {
            computeStats_within_binary_mask(*mask, Mcorr, dummy, dummy, avecorr,
                                            stdcorr);
            double istdcorr = 1.0 / stdcorr;
            FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Mcorr)
            if (DIRECT_MULTIDIM_ELEM(*mask, n))
                DIRECT_MULTIDIM_ELEM(Mcorr, n) =
                    (DIRECT_MULTIDIM_ELEM(Mcorr, n) - avecorr)
                    * istdcorr;
            else
                DIRECT_MULTIDIM_ELEM(Mcorr, n) = 0.;
        }
    }
    else
        Mcorr.statisticsAdjust(0., 1.);
    Mcorr.maxIndex(kmax, imax, jmax);
    max = A3D_ELEM(Mcorr, kmax, imax, jmax);

    // Estimate n_max around the maximum
    int n_max = -1;
    while (neighbourhood)
    {
        n_max++;
        for (int k = -n_max; k <= n_max && neighbourhood; k++)
        {
            k_actual = k + kmax;
            if (k_actual < STARTINGZ(Mcorr) || k_actual > FINISHINGZ(Mcorr))
            {
                neighbourhood = false;
                break;
            }
            for (int i = -n_max; i <= n_max && neighbourhood; i++)
            {
                i_actual = i + imax;
                if (i_actual < STARTINGY(Mcorr) || i_actual > FINISHINGY(Mcorr))
                {
                    neighbourhood = false;
                    break;
                }
                for (int j = -n_max; j <= n_max && neighbourhood; j++)
                {
                    j_actual = j + jmax;
                    if (j_actual < STARTINGX(Mcorr)
                        || j_actual > FINISHINGX(Mcorr))
                    {
                        neighbourhood = false;
                        break;
                    }
                    else if (max
                             / 1.414 > A3D_ELEM(Mcorr, k_actual, i_actual, j_actual))
                    {
                        neighbourhood = false;
                        break;
                    }
                }
            }
        }
    }

    // We have the neighbourhood => looking for the gravity centre
    zmax = xmax = ymax = sumcorr = 0.;
    if (kmax-n_max<STARTINGZ(Mcorr))
        n_max=std::min(kmax-STARTINGZ(Mcorr),n_max);
    if (kmax+n_max>FINISHINGZ(Mcorr))
        n_max=std::min(FINISHINGZ(Mcorr)-kmax,n_max);
    if (imax-n_max<STARTINGY(Mcorr))
        n_max=std::min(imax-STARTINGY(Mcorr),n_max);
    if (imax+n_max>FINISHINGY(Mcorr))
        n_max=std::min(FINISHINGY(Mcorr)-imax,n_max);
    if (jmax-n_max<STARTINGY(Mcorr))
        n_max=std::min(jmax-STARTINGX(Mcorr),n_max);
    if (jmax+n_max>FINISHINGY(Mcorr))
        n_max=std::min(FINISHINGX(Mcorr)-jmax,n_max);
    for (int k = -n_max; k <= n_max; k++)
    {
        k_actual = k + kmax;
        for (int i = -n_max; i <= n_max; i++)
        {
            i_actual = i + imax;
            for (int j = -n_max; j <= n_max; j++)
            {
                j_actual = j + jmax;
                double val = A3D_ELEM(Mcorr, k_actual, i_actual, j_actual);
                zmax += k_actual * val;
                ymax += i_actual * val;
                xmax += j_actual * val;
                sumcorr += val;
            }
        }
    }
    if (sumcorr != 0)
    {
        shiftX = xmax / sumcorr;
        shiftY = ymax / sumcorr;
        shiftZ = zmax / sumcorr;
    }
}

bestShiftRealSpace()

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
	)

Translational search (non-wrapping).

Search is performed in real-space

Definition at line 1989 of file filters.cpp.

{
    I1.checkDimension(2);
    I2.checkDimension(2);

    double bestCorr=-1e38;

    MultidimArray<double> alignedI2;
    MultidimArray<double> bestI2;
    int maxShift2=maxShift*maxShift;
    Matrix1D<double> shift(2);
    for (double y=-maxShift; y<=maxShift; y+=shiftStep)
        for (double x=-maxShift; x<=maxShift; x+=shiftStep)
        {
            if (y*y+x*x>maxShift2)
                continue;
            YY(shift)=y;
            XX(shift)=x;
            translate(xmipp_transformation::LINEAR,alignedI2,I2,shift,xmipp_transformation::DONT_WRAP,0.0);

            double corr=correlationIndex(I1,alignedI2,mask);
            if (corr>bestCorr)
            {
                bestI2=alignedI2;
                bestCorr=corr;
                shiftY=y;
                shiftX=x;
            }
        }
    I2=bestI2;

    return bestCorr;
}

boundMedianFilter()

template<typename T >

void boundMedianFilter	(	MultidimArray< T > &	V,
		const MultidimArray< char > &	mask,
		int	n = 0
	)

Remove bad pixels.

A boundaries median filter is applied at those pixels given by the mask. The mask is set to zeros in the process of the tiltering. Keep a copy in case you further use it.

Definition at line 1309 of file filters.h.

{
    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);
}

centerImage()

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

Center an image both translationally and rotationally

Given an image, this function returns the image that has been centered rotationally and translationally. For doing so, it compares this image with its mirrored (X, Y, XY) versions. The image is aligned translationally and then rotationally Niter times.

Definition at line 3277 of file filters.cpp.

{
    I.checkDimension(2);

    I.setXmippOrigin();
    double avg = I.computeAvg();
    I -= avg;

    MultidimArray<double> Ix;
    MultidimArray<double> Iy;
    MultidimArray<double> Ixy;
    MultidimArray<double> Iaux;
    Matrix2D<double> A;
    A.initIdentity(3);
    Iaux = I;

    MultidimArray<int> mask;
    mask.initZeros(I);
    BinaryCircularMask(mask, XSIZE(I) / 2);

    MultidimArray<double> lineY;
    MultidimArray<double> lineX;
    lineY.initZeros(YSIZE(I));
    STARTINGX(lineY) = STARTINGY(I);
    lineX.initZeros(XSIZE(I));
    STARTINGX(lineX) = STARTINGX(I);

    Polar_fftw_plans *plans = nullptr;
    Polar<std::complex<double> > polarFourierI;
    Polar<std::complex<double> > polarFourierIx;
    MultidimArray<double> rotationalCorr;
    Matrix2D<double> R;
    for (int i = 0; i < Niter; i++)
    {
        // Mask Iaux
        FOR_ALL_ELEMENTS_IN_ARRAY2D(mask)
        if (!A2D_ELEM(mask,i,j))
            A2D_ELEM(Iaux,i,j) = 0;

        // Center translationally
        Ix = Iaux;
        Ix.selfReverseX();
        Ix.setXmippOrigin();
        Iy = Iaux;
        Iy.selfReverseY();
        Iy.setXmippOrigin();
        Ixy = Ix;
        Ixy.selfReverseY();
        Ixy.setXmippOrigin();

        double meanShiftX = 0;
        double meanShiftY = 0;
        double shiftX;
        double shiftY;
        double Nx = 0;
        double Ny = 0;
        bestNonwrappingShift(Iaux, Ix, shiftX, shiftY, aux);
#ifdef DEBUG

        ImageXmipp save;
        save()=Ix;
        save.write("PPPx.xmp");
        std::cout << "con Ix: " << shiftX << " " << shiftY << std::endl;
#endif

        if (fabs(shiftX) < XSIZE(I) / 3 || !limitShift)
        {
            meanShiftX += shiftX;
            Nx++;
        }
        if (fabs(shiftY) < YSIZE(I) / 3 || !limitShift)
        {
            meanShiftY += shiftY;
            Ny++;
        }
        bestNonwrappingShift(Iaux, Iy, shiftX, shiftY, aux);
#ifdef DEBUG

        save()=Iy;
        save.write("PPPy.xmp");
        std::cout << "con Iy: " << shiftX << " " << shiftY << std::endl;
#endif

        if (fabs(shiftX) < XSIZE(I) / 3 || !limitShift)
        {
            meanShiftX += shiftX;
            Nx++;
        }
        if (fabs(shiftY) < YSIZE(I) / 3 || !limitShift)
        {
            meanShiftY += shiftY;
            Ny++;
        }
        bestNonwrappingShift(Iaux, Ixy, shiftX, shiftY, aux);
#ifdef DEBUG

        save()=Ixy;
        save.write("PPPxy.xmp");
        std::cout << "con Ixy: " << shiftX << " " << shiftY << std::endl;
#endif

        if (fabs(shiftX) < XSIZE(I) / 3 || !limitShift)
        {
            meanShiftX += shiftX;
            Nx++;
        }
        if (fabs(shiftY) < YSIZE(I) / 3 || !limitShift)
        {
            meanShiftY += shiftY;
            Ny++;
        }
        if (Nx > 0)
            meanShiftX /= Nx;
        if (Ny > 0)
            meanShiftY /= Ny;

        MAT_ELEM(A,0,2) += -meanShiftX / 2;
        MAT_ELEM(A,1,2) += -meanShiftY / 2;
        Iaux.initZeros();
        applyGeometry(xmipp_transformation::LINEAR, Iaux, I, A, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);
        FOR_ALL_ELEMENTS_IN_ARRAY2D(mask)
        if (!A2D_ELEM(mask,i,j))
            A2D_ELEM(Iaux,i,j) = 0;

#ifdef DEBUG

        std::cout << "Iter " << i << std::endl;
        std::cout << "shift=" << -meanShiftX << "," << -meanShiftY << std::endl;
        save()=I;
        save.write("PPP.xmp");
        save()=Iaux;
        save.write("PPPshift.xmp");
#endif

        // Center rotationally
        Ix = Iaux;
        Ix.selfReverseX();
        Ix.setXmippOrigin();

        polarFourierTransform<true>(Iaux, polarFourierI, true, XSIZE(I) / 5,
                                        XSIZE(I) / 2, plans);
        rotationalCorr.resizeNoCopy(
            2 * polarFourierI.getSampleNoOuterRing() - 1);
        aux2.local_transformer.setReal(rotationalCorr);

        polarFourierTransform<true>(Ix, polarFourierIx, false,
                                        XSIZE(Ix) / 5, XSIZE(Ix) / 2, plans);
        double bestRot = best_rotation(polarFourierIx, polarFourierI, aux2);
        bestRot = realWRAP(bestRot,0,180);
        if (bestRot > 90)
            bestRot = bestRot - 180;

        rotation2DMatrix(bestRot / 2, R);
        A = R * A;
        Iaux.initZeros();
        applyGeometry(xmipp_transformation::LINEAR, Iaux, I, A, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);
        FOR_ALL_ELEMENTS_IN_ARRAY2D(mask)
        if (!A2D_ELEM(mask,i,j))
            A2D_ELEM(Iaux,i,j) = 0;

#ifdef DEBUG

        std::cout << "rot=" << -bestRot/2 << std::endl;
        save()=Iaux;
        save.write("PPProt.xmp");
#endif

        // Remove horizontal/vertical ambiguity
        lineX.initZeros();
        lineY.initZeros();
        FOR_ALL_ELEMENTS_IN_ARRAY2D(Iaux)
        {
            double val = A2D_ELEM(Iaux,i,j);
            if (j == 0)
                lineY(i) = val;
            else if (lineY(i) < val)
                lineY(i) = val;
            if (i == 0)
                lineX(j) = val;
            else if (lineX(j) < val)
                lineX(j) = val;
        }

        double thX = lineX.computeMin()
                     + 0.75 * (lineX.computeMax() - lineX.computeMin());
        double thY = lineY.computeMin()
                     + 0.75 * (lineY.computeMax() - lineY.computeMin());
        int x0 = STARTINGX(lineX);
        while (lineX(x0) < thX)
            x0++;
        int y0 = STARTINGX(lineY);
        while (lineY(y0) < thY)
            y0++;
        int xF = FINISHINGX(lineX);
        while (lineX(xF) < thX)
            xF--;
        int yF = FINISHINGX(lineY);
        while (lineY(yF) < thY)
            yF--;
        if ((xF - x0) > (yF - y0))
        {
            rotation2DMatrix(-90., R);
            A = R * A;
        }
        applyGeometry(xmipp_transformation::LINEAR, Iaux, I, A, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);
#ifdef DEBUG

        lineX.write("PPPlineX.txt");
        lineY.write("PPPlineY.txt");
        std::cout << "dev X=" << xF-x0 << std::endl;
        std::cout << "dev Y=" << yF-y0 << std::endl;
        save()=Iaux;
        save.write("PPPhorver.xmp");
        std::cout << "Press any key\n";
        char c;
        std::cin >> c;
#endif

    }
    applyGeometry(xmipp_transformation::BSPLINE3, Iaux, I, A, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);
    I = Iaux;
    I += avg;
    delete plans;
    return A;
}

centerImageRotationally()

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

Center an image rotationally

Given an image, this function returns the image that has been centered rotationally. For doing so, it compares this image with its mirrored (X) version.

Definition at line 3248 of file filters.cpp.

{
    I.checkDimension(2);

    MultidimArray<double> Ix = I;
    Ix.selfReverseX();
    Ix.setXmippOrigin();

    Polar_fftw_plans *plans = nullptr;
    Polar<std::complex<double> > polarFourierI;
    Polar<std::complex<double> > polarFourierIx;
    polarFourierTransform<true>(Ix, polarFourierIx, false, XSIZE(Ix) / 5,
                                    XSIZE(Ix) / 2, plans);
    polarFourierTransform<true>(I, polarFourierI, true, XSIZE(I) / 5,
                                    XSIZE(I) / 2, plans);

    MultidimArray<double> rotationalCorr;
    rotationalCorr.resize(2 * polarFourierI.getSampleNoOuterRing() - 1);
    aux.local_transformer.setReal(rotationalCorr);
    double bestRot = best_rotation(polarFourierIx, polarFourierI, aux);

    MultidimArray<double> auxI = I;
    rotate(3, I, auxI, -bestRot / 2, xmipp_transformation::WRAP);
    //I.selfRotateBSpline(3,-bestRot/2,WRAP);
}

centerImageTranslationally()

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

Center an image translationally

Given an image, this function returns the image that has been centered translationally. For doing so, it compares this image with its mirrored (X, Y, XY) versions.

Definition at line 3212 of file filters.cpp.

{
    I.checkDimension(2);

    MultidimArray<double> Ix = I;
    Ix.selfReverseX();
    Ix.setXmippOrigin();
    MultidimArray<double> Iy = I;
    Iy.selfReverseY();
    Iy.setXmippOrigin();
    MultidimArray<double> Ixy = Ix;
    Ixy.selfReverseY();
    Ixy.setXmippOrigin();

    double meanShiftX = 0;
    double meanShiftY = 0;
    double shiftX;
    double shiftY;
    bestNonwrappingShift(I, Ix, meanShiftX, meanShiftY, aux);
    bestNonwrappingShift(I, Iy, shiftX, shiftY, aux);
    meanShiftX += shiftX;
    meanShiftY += shiftY;
    bestNonwrappingShift(I, Ixy, shiftX, shiftY, aux);
    meanShiftX += shiftX;
    meanShiftY += shiftY;
    meanShiftX /= 3;
    meanShiftY /= 3;

    Matrix1D<double> shift(2);
    VECTOR_R2(shift, -meanShiftX, -meanShiftY);
    MultidimArray<double> aux2 = I;
    translate(3, I, aux2, shift);
    //I.selfTranslateBSpline(3,shift);
}

contrastEnhancement()

void contrastEnhancement

(

Image< double > *

I

)

Constrast enhancement

The minimum density value is brought to 0 and the maximum to 255.

Definition at line 327 of file filters.cpp.

{
    (*I)().rangeAdjust(0, 255);
}

correlation()

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
	)

Correlation nD

Definition at line 249 of file filters.h.

{
    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;
}

correlationIndex()

template<typename T >

double correlationIndex	(	const MultidimArray< T > &	x,
		const MultidimArray< T > &	y,
		const MultidimArray< int > *	mask = NULL,
		MultidimArray< double > *	Contributions = NULL
	)

correlationIndex nD

Definition at line 3974 of file multidim_array.h.

{
    SPEED_UP_tempsInt;

    double retval = 0, aux;
    double mean_x, mean_y;
    double stddev_x, stddev_y;

    long N = 0;

    if (mask == NULL)
    {
        x.computeAvgStdev(mean_x, stddev_x);
        y.computeAvgStdev(mean_y, stddev_y);
    }
    else
    {
        x.computeAvgStdev_within_binary_mask(*mask, mean_x,stddev_x);
        y.computeAvgStdev_within_binary_mask(*mask, mean_y,stddev_y);
    }
    if (ABS(stddev_x)<XMIPP_EQUAL_ACCURACY ||
        ABS(stddev_y)<XMIPP_EQUAL_ACCURACY)
        return 0;

    // If contributions are desired. Please, be careful interpreting individual
    // contributions to the covariance! One pixel value afect others.
    if (Contributions != NULL)
    {
        FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(x, y)
        {
            if (mask != NULL)
                if (!A3D_ELEM(*mask,k, i, j))
                    continue;

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

        FOR_ALL_ELEMENTS_IN_ARRAY3D(*Contributions)
        A3D_ELEM(*Contributions, k, i, j) /= ((stddev_x * stddev_y) * N);
    }
    else
    {
        if (mask==NULL && x.sameShape(y))
        {
            FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(x)
                retval += DIRECT_MULTIDIM_ELEM(x, n)*DIRECT_MULTIDIM_ELEM(y, n);
            N=MULTIDIM_SIZE(x);
            retval-=N*mean_x*mean_y;
        }
        else
        {
            FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(x, y)
            {
                if (mask != NULL)
                    if (!A3D_ELEM(*mask,k, i, j))
                        continue;

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

    if (N != 0)
        return retval / ((stddev_x * stddev_y) * N);
    else
        return 0;
}

correntropy()

template<typename T >

double correntropy	(	const MultidimArray< T > &	x,
		const MultidimArray< T > &	y,
		double	sigma
	)

Correntropy nD

Definition at line 381 of file filters.h.

{
    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;
}

denoiseTVenergy()

double denoiseTVenergy	(	double	mu,
		const MultidimArray< double > &	X,
		const MultidimArray< double > &	Y,
		double	s,
		double	q,
		double	lambda,
		double	sigmag,
		double	g
	)

Compute energy.

• energy function 0.5*|vst(X)-Y|^2 + mu*TV(X) = data_term + mu*TV

X=current image Y=sensed image mu=regularization weight (on TV norm)

Definition at line 3959 of file filters.cpp.

{
    // parameter beta which role is to prevent division with zeros
    const double beta = 0.00001;
    const double beta2 = beta*beta;

    double m = 0;
    double tv = 0;
    const double K1 = ((3.0/8.0)*lambda*lambda + sigmag*sigmag - lambda*g) / (s*s);
    const double K2 = lambda * (q/(s*s));
    const double K3 = 2.0 / lambda;
    double Xij;
    double Yij;
    double dXx;
    double dXy;
    double d;
    double msqrt;
    int i1;
    int j1;

    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(X)
    {
        Xij = DIRECT_A2D_ELEM(X,i,j);
        Yij = DIRECT_A2D_ELEM(Y,i,j);
        i1 = i + 1;
        if (i1==YSIZE(X))
           i1 = 0;
        j1 = j + 1;
        if (j1==XSIZE(X))
           j1 = 0;

        dXx = DIRECT_A2D_ELEM(X,i,j1) - Xij;
        dXy = DIRECT_A2D_ELEM(X,i1,j) - Xij;
        tv += sqrt(dXx*dXx + dXy*dXy + beta2);

        d = K2*Xij + K1;
        msqrt = K3 * sqrt(std::max(0.0, d)) - Yij;
        m += msqrt*msqrt;
    }
    return 0.5*m + mu*tv;
}

denoiseTVFilter()

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

Compute denoising.

Definition at line 4129 of file filters.cpp.

{
    // parameters for denoising
    double lambda = 1.0;    // Poisson scaling factor
    double sigmag = 5.8;  // Gaussian component N(g,sigma^2)
    double g = 0.0;

    // *** Anscombe transformation of input image ***
    double q = 255.0;

    double xnewmin = xnew.computeMin();
    double xnewscale = 255.0 / (xnew.computeMax() - xnewmin);

    // apply generalized Anscombe VST tranformation
    double K1 = ((3.0/8.0)*lambda*lambda + sigmag*sigmag - lambda*g);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(xnew)
    {
        DIRECT_MULTIDIM_ELEM(xnew, n) = (DIRECT_MULTIDIM_ELEM(xnew, n) - xnewmin) * xnewscale;
        DIRECT_MULTIDIM_ELEM(xnew, n) = 2.0 / lambda * sqrt(std::max(0.0, lambda*DIRECT_MULTIDIM_ELEM(xnew, n) + K1));
    }

    double mx = xnew.computeMax();

    // xold = initial degraded image
    MultidimArray<double> xold = xnew;
    MultidimArray<double> grold = xnew;
    MultidimArray<double> grnew = xnew;
    MultidimArray<double> dold = xnew;
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xnew)
    {
        DIRECT_A2D_ELEM(xold,i,j) = DIRECT_A2D_ELEM(xnew,i,j) / mx;
    }
    MultidimArray<double> origInput = xold;

    // *** Spectral Projected Gradient (SPG) optimization ***
    // parameters of optimization
    double tol = 0;
    double theta;
    double thetamin = 0.001;
    double thetamax = 1000;
    double gamma = 0.0001;
    double sigma1 = 0.1;
    double sigma2 = 0.9;
    double mu = 0.03;

    // initial objective function value of xold
    double fold = denoiseTVenergy(mu, xold, origInput, mx, q, lambda, sigmag, g);

    // gradient of objective function
    denoiseTVgradient(mu, xold, origInput, mx, q, lambda, sigmag, g, grold);

    denoiseTVproj(xold, grold, 1.0, dold);

    double delta;
    double ksi;
    double fnew;
    double s_norm;
    double p;
    double s2;
    double xij;
    for (int kk=1; kk <= maxIter; kk++)
    {
        delta = 0;

        // make a step
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xnew)
        {
            DIRECT_A2D_ELEM(xnew,i,j) = DIRECT_A2D_ELEM(xold,i,j) + DIRECT_A2D_ELEM(dold,i,j);
            delta += DIRECT_A2D_ELEM(grold,i,j) * DIRECT_A2D_ELEM(dold,i,j);
        }

        // objective function of xnew
        fnew = denoiseTVenergy(mu, xnew, origInput, mx, q, lambda, sigmag, g);

        ksi = 1.0;

        while (fnew > fold + gamma*ksi*delta)
        {
            double ksitsl = -0.5 * (ksi*ksi) * delta / (fnew - fold - ksi*delta);
            if ((ksitsl >= sigma1)  && (ksitsl <= sigma2*ksi))
                ksi = ksitsl;
            else
                ksi = ksi / 2.0;

            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xnew)
            {
                DIRECT_A2D_ELEM(xnew,i,j) = DIRECT_A2D_ELEM(xold,i,j) + ksi*DIRECT_A2D_ELEM(dold,i,j);
            }

            // because we updated xnew, update also function value
            fnew = denoiseTVenergy(mu, xnew, origInput, mx, q, lambda, sigmag, g);
        }

        denoiseTVgradient(mu, xnew, origInput, mx, q, lambda, sigmag, g, grnew);

        s_norm = -100;
        p = 0;
        s2 = 0;
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xnew)
        {
            xij = (DIRECT_A2D_ELEM(xnew,i,j) - DIRECT_A2D_ELEM(xold,i,j));
            p += xij * (DIRECT_A2D_ELEM(grnew,i,j) - DIRECT_A2D_ELEM(grold,i,j));
            s2 += xij * xij;
            if (xij > s_norm) s_norm = xij;
        }

        if (s_norm < tol)
            break;

        if (p <= 0)
            theta = thetamax;
        else
            theta = std::min(thetamax, std::max(thetamin, s2/p));

        denoiseTVproj(xnew, grnew, theta, dold);

        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(xold)
        {
            DIRECT_A2D_ELEM(xold,i,j) = DIRECT_A2D_ELEM(xnew,i,j);
            DIRECT_A2D_ELEM(grold,i,j) = DIRECT_A2D_ELEM(grnew,i,j);
        }

        // we already computed it, let's use it for next round
        fold = fnew;
    }
}

denoiseTVgradient()

void denoiseTVgradient	(	double	mu,
		const MultidimArray< double > &	X,
		const MultidimArray< double > &	Y,
		double	s,
		double	q,
		double	lambda,
		double	sigmag,
		double	g,
		MultidimArray< double > &	gradientI
	)

Compute gradient.

• Jacobian of Energy_TV 0.5*|vst(X)-Y|^2 + mu*TV(X)

X=current image Y=sensed image mu=regularization weight (on TV norm)

Definition at line 4017 of file filters.cpp.

{
    // parameter beta which role is to prevent division with zeros
    const double beta = 0.00001;
    const double beta2 = beta*beta;

    MultidimArray<double> d = X;
    double K1 = ((3.0/8.0)*lambda*lambda + sigmag*sigmag - lambda*g) / (s*s);
    double K2 = lambda * (q/s*s);
    double K3 = (2.0 / (lambda*lambda)) * (q / (s*s)) * lambda;
    double Xij;
    double Yij;
    double dXx;
    double dXy;
    double dij;
    double d_left;
    double d_up;
    double X_left;
    double X_right;
    double X_up;
    double X_down;
    double dTV;
    double dE;
    int i1;
    int j1;
    int i_1;
    int j_1;

    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(X)
    {
        Xij = DIRECT_A2D_ELEM(X,i,j);
        i1 = i + 1;
        if (i1==YSIZE(X))
           i1 = 0;
        j1 = j + 1;
        if (j1==XSIZE(X))
           j1 = 0;

        dXx = DIRECT_A2D_ELEM(X,i,j1) - Xij;
        dXy = DIRECT_A2D_ELEM(X,i1,j) - Xij;
        DIRECT_A2D_ELEM(d,i,j) = 1.0 / sqrt(dXx*dXx + dXy*dXy + beta2);
    }

    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(d)
    {
        dij = DIRECT_A2D_ELEM(d,i,j);
        Xij = DIRECT_A2D_ELEM(X,i,j);
        Yij = DIRECT_A2D_ELEM(Y,i,j);
        i1 = i + 1;
        if (i1==YSIZE(d))
           i1 = 0;
        j1 = j + 1;
        if (j1==XSIZE(d))
           j1 = 0;
        i_1 = i - 1;
        if (i_1==-1)
           i_1 = YSIZE(d) - 1;
        j_1 = j - 1;
        if (j_1==-1)
           j_1 = XSIZE(d) - 1;

        d_left = DIRECT_A2D_ELEM(d,i,j_1);
        d_up = DIRECT_A2D_ELEM(d,i_1,j);
        X_left = DIRECT_A2D_ELEM(X,i,j_1);
        X_right = DIRECT_A2D_ELEM(X,i,j1);
        X_up = DIRECT_A2D_ELEM(X,i_1,j);
        X_down = DIRECT_A2D_ELEM(X,i1,j);

        dTV = Xij * (2.0*dij + d_left + d_up) - X_left*d_left - X_up*d_up - dij*(X_right + X_down);

        if (K2*Xij + K1 > 0)
            dE = K3 - (q / (s*s)) * Yij / sqrt(Xij * (q / (s*s)) * lambda + K1);
        else
            dE = 0;

        DIRECT_A2D_ELEM(gradientI,i,j) = dE + mu*dTV;
    }
}

denoiseTVproj()

void denoiseTVproj	(	const MultidimArray< double > &	X,
		const MultidimArray< double > &	Y,
		double	theta,
		MultidimArray< double > &	dold
	)

Function which makes projections on feasible set.

Definition at line 4108 of file filters.cpp.

{
    double div;
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(X)
    {
        div = DIRECT_A2D_ELEM(X,i,j) - (DIRECT_A2D_ELEM(Y,i,j)*theta);
        if (div < 0)
           DIRECT_A2D_ELEM(dold,i,j) = 0 - DIRECT_A2D_ELEM(X,i,j);
        else if (div > 1)
           DIRECT_A2D_ELEM(dold,i,j) = 1 - DIRECT_A2D_ELEM(X,i,j);
        else
           DIRECT_A2D_ELEM(dold,i,j) = div - DIRECT_A2D_ELEM(X,i,j);
    }
}

detectBackground()

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

Detect background

This function receives a Matrix3D vol, and try to find the background assuming that all the outside planes contain background, and apply interval confidence, were alpha is the probabity to fail. Mask is of the same size of vol, and is the solution, mask have value 1 if background else value 0

Definition at line 197 of file filters.cpp.

{

    // 2.1.-Background detection------------------------------------------
    MultidimArray<double> bg; // We create the volumen with
    bg.resizeNoCopy(vol); // -1:not visited 0:mol 1:background
    bg.initConstant(-1); // -2:in the list

    // Ponemos las seis caras de esta variable como visitadas e inicializamos
    // la cola de pixeles por visitar
    std::queue<int> list_for_compute; // Lista del modo [x1,y1,z1,...,xi,yi,zi]
    // que contiene los pixeles por procesar
    std::vector<double> bg_values; // Vector con los valores del background
    size_t xdim = XSIZE(bg);
    size_t ydim = YSIZE(bg);
    size_t zdim = ZSIZE(bg);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(bg)
    {
        if (j == 0 || j == xdim - 1 || i == 0 || i == ydim - 1 || k == 0
            || k == zdim - 1)
        { // Visited coord
            DIRECT_A3D_ELEM(bg,k,i,j) = 1;
            // We introduce the values of the background
            bg_values.push_back(DIRECT_A3D_ELEM(vol,k,i,j));

        }
        if ((j == 1 || j == xdim - 2) && (i != 0) && (i != ydim - 1) && (k != 0)
            && (k != zdim - 1))
        { // Coord for compute for x
            list_for_compute.push(j);
            list_for_compute.push(i);
            list_for_compute.push(k);
        }
        if ((i == 1 || i == ydim - 2) && (j > 1) && (j < xdim - 2) && (k != 0)
            && (k != zdim - 1))
        { // Coord for compute for y
            list_for_compute.push(j);
            list_for_compute.push(i);
            list_for_compute.push(k);
        }
        if ((k == 1 || k == zdim - 2) && (j > 1) && (j < xdim - 2) && (i > 1)
            && (i < ydim - 2))
        { // Coord for compute for y
            list_for_compute.push(j);
            list_for_compute.push(i);
            list_for_compute.push(k);
        }
    } // end of FOR_ALL_ELEMENTS

    // We work until the list_for_compute is empty
    int n = 250; //each 250 pixels renew stats
    int cont = 250; //We start here for compute stat for first time
    double A=0; // A and B are numbers such the interval of confidence is [A,B]
    double B=0; //
    float z = icdf_gauss(1 - alpha / 2);
    while (!list_for_compute.empty())
    {

        //We compute stat when is needed
        if (cont == n)
        {
            // Compute statistics
            double avg=0;
            double stddev=0;
            computeAvgStddev(bg_values, avg, stddev);
            final_mean = avg;
            // Compute confidence interval
            A = avg - (z * stddev);
            B = avg + (z * stddev);
            cont = 0;
        } // end of if
        // Now we start to take coords from the list_for_compute
        int x_coord = list_for_compute.front();
        list_for_compute.pop();
        int y_coord = list_for_compute.front();
        list_for_compute.pop();
        int z_coord = list_for_compute.front();
        list_for_compute.pop();
        // Is visited
        DIRECT_A3D_ELEM(bg,z_coord,y_coord,x_coord) = -2;
        //We take the value
        double value = DIRECT_A3D_ELEM(vol,z_coord,y_coord,x_coord);
        // We see if is background or not
        if (A <= value && value <= B)
        {
            // We now is background
            DIRECT_A3D_ELEM(bg,z_coord,y_coord,x_coord) = 1;
            // We introduce the values of the background
            bg_values.push_back(value);
            // We update the cont variable
            cont++;

            // We add his neighbours in the list_for_compute
            for (int xx = x_coord - 1; xx <= x_coord + 1; xx++)
                for (int yy = y_coord - 1; yy <= y_coord + 1; yy++)
                    for (int zz = z_coord - 1; zz <= z_coord + 1; zz++)
                        //We see if it has been visited
                        if (DIRECT_A3D_ELEM(bg,zz,yy,xx) == -1) // not visited
                        {
                            // So we include it in the list
                            list_for_compute.push(xx);
                            list_for_compute.push(yy);
                            list_for_compute.push(zz);
                            // Is in the list
                            DIRECT_A3D_ELEM(bg,zz,yy,xx) = -2;
                        }
        } // end of if
        else
        {
            // Isn't background
            DIRECT_A3D_ELEM(bg,z_coord,y_coord,x_coord) = 0;
        }
    } // end of while
    // Now we change not visited for mol
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(bg)
    if (DIRECT_A3D_ELEM(bg,k,i,j) == -1)
        DIRECT_A3D_ELEM(bg,k,i,j) = 0;

    // End of 2.1-----------------------------------------------------------
    // 2.2.-Matematical Morphology
    MultidimArray<double> bg_mm; // We create the output volumen
    bg_mm.initZeros(vol);
    closing3D(bg, bg_mm, 26, 0, 1);
    // Output
    //typeCast(bg_mm,mask);
    mask = bg_mm;
}

distanceTransform()

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

L1 distance transform

If wrap is set, the image borders are wrapped around. This is useful if the image coordinates represent angles

Definition at line 584 of file filters.cpp.

{
    std::list<int> toExplore; /* A list of points to explore */

    in.checkDimension(2);

    out.resize(in);
    out.initConstant(XSIZE(in) + YSIZE(in));

#define CHECK_POINT_DIST(i,j,proposedDistance) \
    { \
        int ip=i; \
        int jp=j; \
        if (wrap) \
        { \
            ip=intWRAP(ip,STARTINGY(out),FINISHINGY(out));\
            jp=intWRAP(jp,STARTINGX(out),FINISHINGX(out));\
        } \
        if (ip>=STARTINGY(out) && ip<=FINISHINGY(out) && \
            jp>=STARTINGX(out) && jp<=FINISHINGX(out)) \
            if (out(ip,jp)>proposedDistance) \
            { \
                out(ip,jp)=proposedDistance; \
                toExplore.push_back(ip); \
                toExplore.push_back(jp); \
                toExplore.push_back(proposedDistance); \
            } \
    }

    // Look for all elements in the binary mask and set the corresponding
    // distance to 0
    FOR_ALL_ELEMENTS_IN_ARRAY2D(in)
    if (in(i, j))
    {
        out(i, j) = 0;
        CHECK_POINT_DIST(i-1, j, 1);
        CHECK_POINT_DIST(i+1, j, 1);
        CHECK_POINT_DIST(i, j-1, 1);
        CHECK_POINT_DIST(i, j+1, 1);
    }

    while (!toExplore.empty())
    {
        int i = toExplore.front();
        toExplore.pop_front();
        int j = toExplore.front();
        toExplore.pop_front();
        int proposedDistance = toExplore.front();
        toExplore.pop_front();

        if (proposedDistance == out(i, j))
        {
            // If this is the current distance (i.e., it has not
            // been supersceded by a later value
            CHECK_POINT_DIST(i-1, j, proposedDistance+1);
            CHECK_POINT_DIST(i+1, j, proposedDistance+1);
            CHECK_POINT_DIST(i, j-1, proposedDistance+1);
            CHECK_POINT_DIST(i, j+1, proposedDistance+1);
        }
    }
}

EntropyOtsuSegmentation()

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

Segment an object using a combination of Otsu and Entropy method

The combination aims at minimizing Z(t)=-log10(sigma2B(t))/H(t) Minimizing the intraclass variance in Otsu is the same as maximizing sigma2B. H is the entropy of the two classes in the entropy method.

Then, the lowest percentil of Z is computed. The threshold applied to the volume is the first value in the curve Z(t) falling below this percentil.

The binarization threshold is returned

Definition at line 1145 of file filters.cpp.

{
    //V.checkDimension(3);

    // Compute the probability density function
    Histogram1D hist;
    hist.clear();
    compute_hist(V, hist, 300);
    hist /= hist.sum();

    // Compute the cumulative 0th and 1st order moments
    double x;
    MultidimArray<double> mom0;
    MultidimArray<double> mom1;
    mom0.initZeros(XSIZE(hist));
    mom1.initZeros(XSIZE(hist));
    mom0(0) = hist(0);
    hist.index2val(0, x);
    mom1(0) = hist(0) * x;
    for (size_t i = 1; i < XSIZE(mom0); i++)
    {
        mom0(i) = mom0(i - 1) + hist(i);
        hist.index2val(i, x);
        mom1(i) = mom1(i - 1) + hist(i) * x;
    }

    // Entropy for black and white parts of the histogram
    const double epsilon = 1e-15;
    MultidimArray<double> h1;
    MultidimArray<double> h2;
    h1.initZeros(XSIZE(hist));
    h2.initZeros(XSIZE(hist));
    for (size_t i = 0; i < XSIZE(hist); i++)
    {
        // Entropy h1
        double w1 = mom0(i);
        if (w1 > epsilon)
            for (size_t ii = 0; ii <= i; ii++)
                if (hist(ii) > epsilon)
                {
                    double aux = hist(ii) / w1;
                    h1(i) -= aux * log10(aux);
                }

        // Entropy h2
        double w2 = 1 - mom0(i);
        if (w2 > epsilon)
            for (size_t ii = i + 1; ii < XSIZE(hist); ii++)
                if (hist(ii) > epsilon)
                {
                    double aux = hist(ii) / w2;
                    h2(i) -= aux * log10(aux);
                }
    }

    // Compute sigma2B and H
    MultidimArray<double> sigma2B;
    MultidimArray<double> H;
    MultidimArray<double> HSigma2B;
    sigma2B.initZeros(XSIZE(hist) - 1);
    H.initZeros(XSIZE(hist) - 1);
    HSigma2B.initZeros(XSIZE(hist) - 1);
    for (size_t i = 0; i < XSIZE(hist) - 1; i++)
    {
        double w1 = mom0(i);
        double w2 = 1 - mom0(i);
        double mu1 = mom1(i);
        double mu2 = mom1(XSIZE(mom1) - 1) - mom1(i);
        sigma2B(i) = w1 * w2 * (mu1 - mu2) * (mu1 - mu2);
        H(i) = h1(i) + h2(i);
        HSigma2B(i) = -log10(sigma2B(i)) / H(i);
        // The logic behind this expression is
        // Otsu:    max sigma2B -> max log10(sigma2B) -> min -log10(sigma2B)
        // Entropy: max H       -> max H              -> min 1/H
    }

    // Sort HSigma2B and take a given percentage of it
    MultidimArray<double> HSigma2Bsorted;
    HSigma2B.sort(HSigma2Bsorted);
    int iTh = ROUND(XSIZE(HSigma2B)*percentil);
    double threshold = HSigma2Bsorted(iTh);

    // Find the first value within HSigma2B falling below this threshold
    iTh = 0;
    while (A1D_ELEM(HSigma2B,iTh) >= threshold)
        iTh++;
    iTh--;
    if (iTh <= 0)
        x = threshold;
    else
        hist.index2val(iTh, x);

    if (binarizeVolume)
        V.binarize(x);
    return x;
}

EntropySegmentation()

double EntropySegmentation

(

MultidimArray< double > &

V

)

Segment an object using Entropy method

Entropy method determines a threshold such that the entropy of the two classes is maximized

http://rsbweb.nih.gov/ij/plugins/download/Entropy_Threshold.java

Definition at line 1078 of file filters.cpp.

{
    // V.checkDimension(3);

    // Compute the probability density function
    Histogram1D hist;
    hist.clear();
    compute_hist(V, hist, 200);
    hist /= hist.sum();

    // Compute the cumulative 0th and 1st order moments
    double x;
    MultidimArray<double> mom0;
    mom0.initZeros(XSIZE(hist));
    mom0(0) = hist(0);
    for (size_t i = 1; i < XSIZE(mom0); i++)
        mom0(i) = mom0(i - 1) + hist(i);

    // Entropy for black and white parts of the histogram
    const double epsilon = 1e-15;
    MultidimArray<double> h1;
    MultidimArray<double> h2;
    h1.initZeros(XSIZE(hist));
    h2.initZeros(XSIZE(hist));
    for (size_t i = 0; i < XSIZE(hist); i++)
    {
        // Entropy h1
        double w1 = mom0(i);
        if (w1 > epsilon)
            for (size_t ii = 0; ii <= i; ii++)
                if (hist(ii) > epsilon)
                {
                    double aux = hist(ii) / w1;
                    h1(i) -= aux * log10(aux);
                }

        // Entropy h2
        double w2 = 1 - mom0(i);
        if (w2 > epsilon)
            for (size_t ii = i + 1; ii < XSIZE(hist); ii++)
                if (hist(ii) > epsilon)
                {
                    double aux = hist(ii) / w2;
                    h2(i) -= aux * log10(aux);
                }
    }

    // Find histogram index with maximum entropy
    double Hmax = h1(0) + h2(0);
    size_t iHmax = 0;
    for (size_t i = 1; i < XSIZE(hist) - 1; i++)
    {
        double H = h1(i) + h2(i);
        if (H > Hmax)
        {
            Hmax = H;
            iHmax = i;
        }
    }

    hist.index2val(iHmax, x);
    V.binarize(x);

    return x;
}

estimateGaussian2D()

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

Fit Gaussian spot to an image.

The fitted Gaussian is a*G(r,mu,sigma)+b where G(r,mu,sigma)=exp(-0.5 * (r-mu)^t sigma^-1 (r-mu))

You can choose if the center is estimated or it is assumed to be 0. You can choose the number of iterations for the estiamtion.

Definition at line 2390 of file filters.cpp.

{
    I.checkDimension(2);

    MultidimArray<double> z(I);

    // Estimate b
    Histogram1D hist;
    compute_hist(z, hist, 100);
    b = hist.percentil(5);

    // Iteratively estimate all parameters
    for (int n = 0; n < iterations; n++)
    {
        // Reestimate z
        FOR_ALL_ELEMENTS_IN_ARRAY2D(z)
        z(i, j) = XMIPP_MAX(I(i,j)-b,0);

        // Sum of z
        double T = z.sum();

        // Estimate center
        mu.initZeros(2);
        if (estimateMu)
        {
            FOR_ALL_ELEMENTS_IN_ARRAY2D(z)
            {
                double val = z(i, j);
                XX(mu) += val * j;
                YY(mu) += val * i;
            }
            mu /= T;
        }

        // Estimate sigma
        sigma.initZeros(2, 2);
        FOR_ALL_ELEMENTS_IN_ARRAY2D(z)
        {
            double val = z(i, j);
            double j_mu = j - XX(mu);
            double i_mu = i - YY(mu);
            sigma(0, 0) += val * j_mu * j_mu;
            sigma(0, 1) += val * i_mu * j_mu;
            sigma(1, 1) += val * i_mu * i_mu;
        }
        sigma(1, 0) = sigma(0, 1);
        sigma /= T;

        // Estimate amplitude
        Matrix2D<double> sigmainv = sigma.inv();
        Matrix1D<double> r(2);
        double G2 = 0;
        a = 0;
        FOR_ALL_ELEMENTS_IN_ARRAY2D(z)
        {
            XX(r) = j;
            YY(r) = i;
            double G = unnormalizedGaussian2D(r, mu, sigmainv);
            a += z(i, j) * G;
            G2 += G * G;
        }
        a /= G2;

        // Reestimate b
        FOR_ALL_ELEMENTS_IN_ARRAY2D(z)
        {
            XX(r) = j;
            YY(r) = i;
            double G = unnormalizedGaussian2D(r, mu, sigmainv);
            z(i, j) = I(i, j) - a * G;
        }
        compute_hist(z, hist, 100);
        b = hist.percentil(5);
    }
}

euclidianDistance()

template<typename T >

double euclidianDistance	(	const MultidimArray< T > &	x,
		const MultidimArray< T > &	y,
		const MultidimArray< int > *	mask = nullptr
	)

euclidian distance nD

Definition at line 672 of file filters.h.

{
    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;
}

fastCorrentropy()

template<typename T >

double fastCorrentropy	(	const MultidimArray< T > &	x,
		const MultidimArray< T > &	y,
		double	sigma,
		const GaussianInterpolator &	G
	)

Fast Correntropy 1D

Definition at line 361 of file filters.h.

{
    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;
}

fastMaskedCorrelation()

template<typename T >

double fastMaskedCorrelation	(	const MultidimArray< T > &	x,
		const MultidimArray< T > &	y,
		const MultidimArray< int > &	mask
	)

Correlation nD

Definition at line 310 of file filters.h.

{
    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;
}

fillBinaryObject()

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

Fill object

Everything that is not background is assumed to be object.

Definition at line 758 of file filters.cpp.

{
    I.checkDimension(2);

    MultidimArray<double> label;
    FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
    I(i, j) = 1 - I(i, j);
    labelImage2D(I, label, neighbourhood);
    double l0 = label(STARTINGY(I), STARTINGX(I));
    FOR_ALL_ELEMENTS_IN_ARRAY2D(label)
    if (label(i, j) == l0)
        I(i, j) = 0;
    else
        I(i, j) = 1;
}

fillTriangle()

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

Fill a triangle defined by three points

The points are supplied as a pointer to three integer positions. They can be negative

Definition at line 3012 of file filters.cpp.

{
    img.checkDimension(2);

    /*
     * Order in y values
     */
    int y1 = 0;
    while (!y1)
    {
        y1 = 1;
        for (int y2 = 0; y2 < 2; y2++)
        {
            if (ty[y2] > ty[y2 + 1]
                || (ty[y2] == ty[y2 + 1] && tx[y2] < tx[y2 + 1]))
            {
                int x1 = ty[y2];
                ty[y2] = ty[y2 + 1];
                ty[y2 + 1] = x1;
                x1 = tx[y2];
                tx[y2] = tx[y2 + 1];
                tx[y2 + 1] = x1;
                y1 = 0;
            }
        }
    }

    int dx1 = tx[1] - tx[0];
    int dx2 = tx[2] - tx[0];
    int dy1 = ty[1] - ty[0];
    int dy2 = ty[2] - ty[0];

    int sx1 = SGN0(dx1);
    int sx2 = SGN0(dx2);
    int sy1 = SGN0(dy1);

    int ix1 = ABS(dx1);
    int ix2 = ABS(dx2);
    int iy1 = ABS(dy1);
    int iy2 = ABS(dy2);

    int inc1 = XMIPP_MAX(ix1, iy1);
    int inc2 = XMIPP_MAX(ix2, iy2);

    int x1;
    int x2;
    int y2;
    int xl;
    int xr;
    x1 = x2 = y1 = y2 = 0;
    xl = xr = tx[0];
    int y = ty[0];

    while (y != ty[1])
    {
        int doit1 = 0;
        int doit2 = 0;

        while (!doit1)
        { /* Wait until y changes */
            x1 += ix1;
            y1 += iy1;
            if (x1 > inc1)
            {
                x1 -= inc1;
                xl += sx1;
            }
            if (y1 > inc1)
            {
                y1 -= inc1;
                y += sy1;
                doit1 = 1;
            }
        }

        while (!doit2)
        { /* Wait until y changes */
            x2 += ix2;
            y2 += iy2;
            if (x2 > inc2)
            {
                x2 -= inc2;
                xr += sx2;
            }
            if (y2 > inc2)
            {
                y2 -= inc2;
                doit2 = 1;
            }
        }

        for (int myx = xl; myx <= xr; myx++)
            img(y, myx) = color;
    }

    dx1 = tx[2] - tx[1];
    dy1 = ty[2] - ty[1];

    sx1 = SGN0(dx1);
    sy1 = SGN0(dy1);

    ix1 = ABS(dx1);
    iy1 = ABS(dy1);

    inc1 = XMIPP_MAX(ix1, iy1);
    xl = tx[1];
    x1 = 0;

    while (y != ty[2])
    {
        int doit1 = 0;
        int doit2 = 0;

        while (!doit1)
        { /* Wait until y changes */
            x1 += ix1;
            y1 += iy1;
            if (x1 > inc1)
            {
                x1 -= inc1;
                xl += sx1;
            }
            if (y1 > inc1)
            {
                y1 -= inc1;
                y += sy1;
                doit1 = 1;
            }
        }

        while (!doit2)
        { /* Wait until y changes */
            x2 += ix2;
            y2 += iy2;
            if (x2 > inc2)
            {
                x2 -= inc2;
                xr += sx2;
            }
            if (y2 > inc2)
            {
                y2 -= inc2;
                doit2 = 1;
            }
        }

        for (int myx = xl; myx <= xr; myx++)
            img(y, myx) = color;
    }
}

forcePositive()

void forcePositive

(

MultidimArray< double > &

V

)

Force positive.

A median filter is applied at those negative values. Positive values are untouched.

Force positive -----------------------------------------------------—

Definition at line 3506 of file filters.cpp.

{
    MultidimArray<char> mask(ZSIZE(V), YSIZE(V), XSIZE(V));
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(V)
    {
        double x = DIRECT_MULTIDIM_ELEM(V, n);
        DIRECT_MULTIDIM_ELEM(mask, n) = (x <= 0);
    }
    boundMedianFilter(V, mask);
}

fourierBesselDecomposition()

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

Fourier-Bessel decomposition

The Fourier-Bessel decomposition of those pixels in img_in whose radius is between r1 and r2 is computed. r1 and r2 are supposed to fit in the image shape, and the image logical origin is used for the decomposition. k1 and k2 determines the harmonic coefficients to be computed.

Definition at line 2469 of file filters.cpp.

{
    img_in.checkDimension(2);

    for (int k = k1; k <= k2; k++)
    {
        int k_1 = k - 1;

        // Compute h and a,b coefficients
        double h = 0;
        double my5 = 0;
        if (k_1 != 0)
        {
            double my = 1 + PI * r2 / 2 / k_1;
            double my4 = my * k_1;
            double ntot = 4 * my4;
            h = 2 * PI / ntot;
        }
        else
        {
            double my = 1 + PI * r2 / 2;
            double my4 = my;
            double ntot = 4 * my4;
            h = 2 * PI / ntot;
        }

        MultidimArray<double> sine(CEIL(my5));
        FOR_ALL_ELEMENTS_IN_ARRAY1D(sine)
        sine(i) = sin((i + 1) * h);

    }
}

giniCoeff()

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

Returns the Gini coefficient of an image. This is related to the Entropy. It also applies a variance filter to the input Image

Definition at line 972 of file filters.cpp.

{
    MultidimArray<double> im = I;

    // std::cout << " - Starting fft filtering " << std::endl;

    FourierFilter filter;
    filter.FilterShape = REALGAUSSIAN;
    filter.FilterBand = LOWPASS;
    filter.w1 = 4;
    filter.applyMaskSpace(im);
   
    // Image<double> imG(im);
    // imG.write("I_Gauss.mrc");

    // std::cout << " - Calling varianceFilter() " << std::endl;
    varianceFilter(I, varKernelSize, true);
    im = I;

    // std::cout << " - Starting 2nd fft filtering " << std::endl;
    filter.w1 = varKernelSize/8;
    filter.applyMaskSpace(I);

    // Image<double> imGV(im);
    // imGV.write("I_Gauss_Var.mrc");
   
    im -= im.computeMin();
    im /= im.computeMax();

    // Image<double> imGVN(im);
    // imGVN.write("I_Gauss_Var_Norm.mrc");
   
    // std::cout << " - Starting histogram analysis " << std::endl;
    Histogram1D hist;
    hist.clear();
    compute_hist(im, hist, 256);

    // std::cout << " - Computing Gini coeff " << std::endl;
    MultidimArray<double> sortedList=hist;
    hist.sort(sortedList);
    double height=0;
    double area=0;
    for (int i=0; i<XSIZE(sortedList); i++)
    {
        height += DIRECT_MULTIDIM_ELEM(sortedList,i);
        area += height - DIRECT_MULTIDIM_ELEM(sortedList,i)/2.0;
    }
        
    double fair_area = height*XSIZE(hist)/2.0;

    double giniValue = (fair_area-area)/fair_area;

    return giniValue;
}

inertiaMoments()

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

Inertia moments

They are measured with respect to the center of the image, and not with respect to the center of mass. For an image there are only two inertia moments. v_out contains the inertia moments while the columns of u contain the directions of the principal axes.

Definition at line 2970 of file filters.cpp.

{
    img.checkDimension(2);

    // Prepare some variables
    double m_11 = 0;
    double m_02 = 0;
    double m_20 = 0;
    double normalize_x = 2.0 / XSIZE(img);
    double normalize_y = 2.0 / YSIZE(img);

    // Compute low-level moments
    FOR_ALL_ELEMENTS_IN_ARRAY2D(img)
    {
        if (mask != nullptr)
            if (!(*mask)(i, j))
                continue;
        double val = img(i, j);
        double dx = j * normalize_x;
        double dy = i * normalize_y;
        double dx2 = dx * dx;
        double dy2 = dy * dy;
        m_11 += val * dx * dy;
        m_02 += val * dy2;
        m_20 += val * dx2;
    }

    // Compute the eigen values of the inertia matrix
    // [m_02 m_11
    //  m_11 m_20]
    Matrix2D<double> A(2, 2);
    A(0, 0) = m_02;
    A(0, 1) = A(1, 0) = m_11;
    A(1, 1) = m_20;
    Matrix2D<double> v;
    svdcmp(A, u, v_out, v);
    v_out = v_out.sort();
}

keepBiggestComponent()

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

Keep the biggest connected component

If the biggest component does not cover the percentage required (by default, 0), more big components are taken until this is accomplished.

Definition at line 713 of file filters.cpp.

{
    MultidimArray<double> label;
    int imax;
    if (ZSIZE(I)==1)
        imax=labelImage2D(I, label, neighbourhood);
    else
        imax=labelImage3D(I, label);
    MultidimArray<int> nlabel(imax + 1);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(label)
    {
        int l=(int)DIRECT_MULTIDIM_ELEM(label,n);
        if (l>0)
            A1D_ELEM(nlabel,l)++;
    }

    MultidimArray <int> best;
    nlabel.indexSort(best);
    best -= 1;
    int nbest = XSIZE(best) - 1;
    double total = nlabel.sum();
    double explained = nlabel(best(nbest));
    while (explained < percentage * total)
    {
        nbest--;
        explained += nlabel(best(nbest));
    }

    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(label)
    {
        int l=(int)DIRECT_MULTIDIM_ELEM(label,n);
        bool among_the_best = false;
        for (int k = nbest; k < imax + 1; k++)
            if (l == A1D_ELEM(best,k))
            {
                among_the_best = true;
                break;
            }
        if (!among_the_best)
            DIRECT_MULTIDIM_ELEM(I,n)=0;
    }
}

labelImage2D()

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

Label a binary image

This function receives a binary volume and labels all its connected components. The background is labeled as 0, and the components as 1, 2, 3 ...

Definition at line 648 of file filters.cpp.

{
    I.checkDimension(2);

    label = I;
    int colour = 32000;
    FOR_ALL_ELEMENTS_IN_ARRAY2D(label)
    {
        if (label(i, j) != 1)
            continue;
        regionGrowing2D(label, label, i, j, 0, colour, false, neighbourhood);
        colour++;
    }
    FOR_ALL_ELEMENTS_IN_ARRAY2D(label)
    if (label(i, j) != 0)
        label(i, j) = label(i, j) - 31999;
    return colour - 32000;
}

labelImage3D()

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

Label a binary volume

This function receives a binary image and labels all its connected components. The background is labeled as 0, and the components as 1, 2, 3 ...

Definition at line 669 of file filters.cpp.

{
    V.checkDimension(3);

    label = V;
    int colour = 32000;
    FOR_ALL_ELEMENTS_IN_ARRAY3D(label)
    {
        if (label(k, i, j) != 1)
            continue;
        regionGrowing3D(label, label, k, i, j, 0, colour, false);
        colour++;
    }
    FOR_ALL_ELEMENTS_IN_ARRAY3D(label)
    if (A3D_ELEM(label,k, i, j) != 0)
        A3D_ELEM(label,k, i, j) = A3D_ELEM(label,k, i, j) - 31999;
    return colour - 32000;
}

localThresholding()

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

Local thresholding

This function implements the local thresholding as described in http://homepages.inf.ed.ac.uk/rbf/HIPR2/adpthrsh.htm A mask can be supplied to limit the image area.

The procedure employed is the following:

• Convolving the image with the statistical operator, i.e. the mean or median (the size of the convolution kernel is dimLocal)
• Subtracting the original from the convolved image.
• Thresholding the difference image with C.
• Inverting the thresholded image.

Definition at line 3164 of file filters.cpp.

{

    // Convolve the input image with the kernel
    MultidimArray<double> convolved;
    convolved.initZeros(img);
    FOR_ALL_ELEMENTS_IN_ARRAY2D(convolved)
    {
        if (mask != nullptr)
            if (!(*mask)(i, j))
                continue;
        int ii0 = XMIPP_MAX(STARTINGY(convolved), FLOOR(i - dimLocal));
        int jj0 = XMIPP_MAX(STARTINGX(convolved), FLOOR(j - dimLocal));
        int iiF = XMIPP_MIN(FINISHINGY(convolved), CEIL(i + dimLocal));
        int jjF = XMIPP_MIN(FINISHINGX(convolved), CEIL(j + dimLocal));
        double N = 0;
        for (int ii = ii0; ii <= iiF; ii++)
            for (int jj = jj0; jj <= jjF; jj++)
            {
                if (mask == nullptr)
                {
                    convolved(i, j) += img(ii, jj);
                    ++N;
                }
                else if ((*mask)(i, j))
                {
                    convolved(i, j) += img(ii, jj);
                    ++N;
                }
            }
        if (N != 0)
            convolved(i, j) /= N;
    }

    // Subtract the original from the convolved image and threshold
    result.initZeros(img);
    FOR_ALL_ELEMENTS_IN_ARRAY2D(convolved)
    {
        if (mask != nullptr)
            if (!(*mask)(i, j))
                continue;
        if (img(i, j) - convolved(i, j) > C)
            result(i, j) = 1;
    }
}

logFilter()

template<typename T >

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

Compute logarithm.

apply transformation a+b*ln(x+c). i.e: 4.431-0.4018*LN(P1+336.6)

Definition at line 1411 of file filters.h.

{

    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));
    }
}

medianFilter3x3()

template<typename T >

void medianFilter3x3	(	MultidimArray< T > &	m,
		MultidimArray< T > &	out
	)

Median_filter with a 3x3 selfWindow

Definition at line 1088 of file filters.h.

{
    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;
}

mutualInformation()

template<typename T >

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

mutual information nD

Return the mutual information: MI = sum [ P(x,y)*log2{P(x,y)/(P(x)*P(y))} ] in the common positions. P(x), P(y) are 1D-histograms of the values of matrix x and y. P(x,y) is the 2D-histogram, i.e. the count of times that a certain combination of values in matrices x and y has occurred. The sum runs over all histogram bins.

The histograms are calculated using the number of bins nx and ny. If no values (or zeros) are given, a Gaussian distribution of the values in the matrices is assumed, and the number of bins is calculated as: log2(n)+1. (according to: Tourassi et al. (2001) Med. Phys. 28 pp. 2394-2402.)

Definition at line 4291 of file filters.cpp.

{
    SPEED_UP_temps;

    long n = 0;
    Histogram1D histx;
    Histogram1D histy;
    Histogram2D histxy;
    MultidimArray< T > aux_x;
    MultidimArray< T > aux_y;
    MultidimArray< double > mx;
    MultidimArray< double > my;
    MultidimArray< double > mxy;
    int xdim;
    int ydim;
    int zdim;
    double retval = 0.0;

    xdim=XSIZE(x);
    ydim=YSIZE(x);
    zdim=ZSIZE(x);
    aux_x.resize(xdim * ydim * zdim);
    xdim=XSIZE(y);
    ydim=YSIZE(y);
    zdim=ZSIZE(y);
    aux_y.resize(xdim * ydim * zdim);

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

        aux_x(n) = A3D_ELEM(x, k, i, j);
        aux_y(n) = A3D_ELEM(y, k, i, j);
        n++;
    }

    aux_x.resize(n);
    aux_y.resize(n);

    if (n != 0)
    {
        if (nx == 0)
            //Assume Gaussian distribution
            nx = (int)((log((double) n) / LOG2) + 1);

        if (ny == 0)
            //Assume Gaussian distribution
            ny = (int)((log((double) n) / LOG2) + 1);

        compute_hist(aux_x, histx, nx);
        compute_hist(aux_y, histy, ny);
        compute_hist(aux_x, aux_y, histxy, nx, ny);

        mx = histx;
        my = histy;
        mxy = histxy;
        for (int i = 0; i < nx; i++)
        {
            double histxi = (histx(i)) / n;
            for (int j = 0; j < ny; j++)
            {
                double histyj = (histy(j)) / n;
                double histxyij = (histxy(i, j)) / n;
                if (histxyij > 0)
                    retval += histxyij * log(histxyij / (histxi * histyj)) /
                              LOG2;
            }
        }

        return retval;
    }
    else
        return 0;
}

noisyZonesFilter()

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

Transforms I to a binary mask with 0 where both variance and mean are high.

Definition at line 870 of file filters.cpp.

{
    int kernelSize_2 = kernelSize/2;
    MultidimArray<double> kernel;
    kernel.resize(kernelSize,kernelSize);
    kernel.setXmippOrigin();

    MultidimArray<double> mAvg=I;
    MultidimArray<double> mVar=I;
    double stdKernel;
    double varKernel;
    double avgKernel;
    double min_val;
    double max_val;
    int x0;
    int y0;
    int xF;
    int yF;
    
    for (int i=kernelSize_2; i<(int)YSIZE(I); i+=kernelSize_2)
        for (int j=kernelSize_2; j<(int)XSIZE(I); j+=kernelSize_2)
            {
                x0 = j-kernelSize_2;
                y0 = i-kernelSize_2;
                xF = j+kernelSize_2-1;
                yF = i+kernelSize_2-1;

                if (x0 < 0)
                    x0 = 0;
                if (xF > XSIZE(I))
                    xF = XSIZE(I);
                if (y0 < 0)
                    y0 = 0;
                if (yF > YSIZE(I))
                    yF = YSIZE(I);

                I.window(kernel, y0, x0, yF, xF);
                kernel.computeStats(avgKernel, stdKernel, min_val, max_val);
                varKernel = stdKernel*stdKernel;

                DIRECT_A2D_ELEM(mAvg, i, j) = avgKernel*avgKernel;
                DIRECT_A2D_ELEM(mVar, i, j) = varKernel;
            }

    // filtering to fill the matrices (convolving with a Gaussian)
    FourierFilter filter;
    filter.FilterShape = REALGAUSSIAN;
    filter.FilterBand = LOWPASS;
    filter.w1 = kernelSize_2;
    filter.applyMaskSpace(mAvg);
    filter.applyMaskSpace(mVar);

    // // Draw to debug
    // Image<double> imAvg(mAvg), imVar(mVar);
    // imAvg.write("AvgFilter.mrc");
    // imVar.write("VarFilter.mrc");

    // Working in a auxilary windows to avoid borders bad defined
    auto ysize = static_cast<int>YSIZE(I);
    auto xsize = static_cast<int>XSIZE(I);
    MultidimArray<double> mAvgAux(ysize-kernelSize,xsize-kernelSize);
    MultidimArray<double> mVarAux(ysize-kernelSize,xsize-kernelSize);
    mAvgAux.setXmippOrigin();
    mVarAux.setXmippOrigin();
    mAvg.window(mAvgAux,STARTINGY(mAvg)+kernelSize_2, STARTINGX(mAvg)+kernelSize_2,
                       FINISHINGY(mAvg)-kernelSize_2, FINISHINGX(mAvg)-kernelSize_2);
    mVar.window(mVarAux,STARTINGY(mVar)+kernelSize_2, STARTINGX(mVar)+kernelSize_2,
                       FINISHINGY(mVar)-kernelSize_2, FINISHINGX(mVar)-kernelSize_2);

    // Refiltering to get a smoother distribution
    // filter.w1 = XSIZE(I)/40;
    // filter.applyMaskSpace(mAvgAux);
    // filter.applyMaskSpace(mVarAux);

    // Binarization
    MultidimArray<double> mAvgAuxBin = mAvgAux;
    MultidimArray<double> mVarAuxBin = mVarAux;
    EntropySegmentation(mVarAuxBin);
    // EntropySegmentation(mAvgAuxBin);
    float th = EntropySegmentation(mAvgAuxBin);
    mAvgAuxBin.binarize(th*0.92);
    mAvgAuxBin = 1-mAvgAuxBin;
    // std::cout << "binarize threshold = " << th << std::endl;

    // Returning to the previous windows size
    MultidimArray<double> mAvgBin = mAvg;
    MultidimArray<double> mVarBin = mVar;
    mAvgAuxBin.window(mAvgBin, STARTINGY(mVar), STARTINGX(mVar),
                              FINISHINGY(mVar), FINISHINGX(mVar));
    mVarAuxBin.window(mVarBin, STARTINGY(mVar), STARTINGX(mVar),
                              FINISHINGY(mVar), FINISHINGX(mVar));

    // // Draw to debug
    // Image<double> imAvgBin(mAvgBin), imVarBin(mVarBin);
    // imAvgBin.write("noisyZoneFilter_AVGmask.mrc");
    // imVarBin.write("noisyZoneFilter_VARmask.mrc");

    // Combining both masks
    I = 1-(mVarBin*mAvgBin);
}

OtsuSegmentation()

double OtsuSegmentation

(

MultidimArray< double > &

V

)

Segment an object using Otsu's method

Otsu's method determines a threshold such that the variance of the two classes is minimized

http://www.biomecardio.com/matlab/otsu.html

Definition at line 1028 of file filters.cpp.

{
    // V.checkDimension(3);

    // Compute the probability density function
    Histogram1D hist;
    hist.clear();
    compute_hist(V, hist, 200);
    hist /= hist.sum();

    // Compute the cumulative 0th and 1st order moments
    double x;
    MultidimArray<double> mom0;
    MultidimArray<double> mom1;
    mom0.initZeros(XSIZE(hist));
    mom1.initZeros(XSIZE(hist));
    mom0(0) = hist(0);
    hist.index2val(0, x);
    mom1(0) = hist(0) * x;
    for (size_t i = 1; i < XSIZE(mom0); i++)
    {
        mom0(i) = mom0(i - 1) + hist(i);
        hist.index2val(i, x);
        mom1(i) = mom1(i - 1) + hist(i) * x;
    }

    // Maximize sigma2B
    double bestSigma2B = -1;
    int ibestSigma2B = -1;
    for (size_t i = 0; i < XSIZE(hist) - 1; i++)
    {
        double w1 = mom0(i);
        double w2 = 1 - mom0(i);
        double mu1 = mom1(i);
        double mu2 = mom1(XSIZE(mom1) - 1) - mom1(i);
        double sigma2B = w1 * w2 * (mu1 - mu2) * (mu1 - mu2);
        if (sigma2B > bestSigma2B)
        {
            bestSigma2B = sigma2B;
            ibestSigma2B = i;
        }
    }

    hist.index2val(ibestSigma2B, x);
    V.binarize(x);

    return x;
}

pixelDesvFilter()

template<typename T >

void pixelDesvFilter	(	MultidimArray< T > &	V,
		double	thresFactor
	)

Remove bad pixels.

A boundaries median filter is applied at those pixels whose value is out of range given by thresFactor * std.

Definition at line 1378 of file filters.h.

{
    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);
    }
}

regionGrowing2D()

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
	)

Region growing for images

Given a position inside an image, this function grows a region with (filling_colour) until it finds a border of value (stop_colour). If the point is outside the image then nothing is done.

If less is true the region is grown in sucha a way that all voxels in its border are greater than the region voxels. If less is false the region is grown so that all voxels on its border are smaller than the region voxels.

Valid neighbourhoods are 4 or 8.

Definition at line 340 of file filters.cpp.

{
    I_in.checkDimension(2);

    std::list<Coordinate2D> iNeighbours; /* A list for neighbour pixels */
    int iCurrenti;
    int iCurrentj; /* Coordinates of the current pixel considered */

    /* First task is copying the input image into the output one */
    I_out = I_in;

    /**** Then the region growing is done ****/
    /* Insert at the beginning of the list the seed coordinates */
    Coordinate2D coord;
    coord.ii = i;
    coord.jj = j;
    iNeighbours.push_front(coord);

    /* Fill the seed coordinates */
    A2D_ELEM(I_out, i, j) = filling_colour;

    while (!iNeighbours.empty())
    {
        /* Take the current pixel to explore */
        coord = iNeighbours.front();
        iNeighbours.pop_front();
        iCurrenti = coord.ii;
        iCurrentj = coord.jj;

#define CHECK_POINT(i,j) \
    { \
  int I=i; \
        int J=j; \
  if (INSIDEXY(I_out,J,I))  { \
   double pixel=A2D_ELEM(I_out,I,J);\
   if (pixel!=filling_colour) \
    if ((less && pixel < stop_colour) || \
     (!less && pixel > stop_colour)) { \
     coord.ii=I; \
     coord.jj=J; \
     A2D_ELEM (I_out,coord.ii,coord.jj)=filling_colour; \
     iNeighbours.push_front(coord); \
    } \
  } \
    }

        /* Make the exploration of the pixel's neighbours */
        CHECK_POINT(iCurrenti, iCurrentj - 1);
        CHECK_POINT(iCurrenti, iCurrentj + 1);
        CHECK_POINT(iCurrenti - 1, iCurrentj);
        CHECK_POINT(iCurrenti + 1, iCurrentj);
        if (neighbourhood == 8)
        {
            CHECK_POINT(iCurrenti - 1, iCurrentj - 1);
            CHECK_POINT(iCurrenti - 1, iCurrentj + 1);
            CHECK_POINT(iCurrenti + 1, iCurrentj - 1);
            CHECK_POINT(iCurrenti + 1, iCurrentj + 1);
        }
    }
}

removeSmallComponents()

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

Remove connected components

Remove connected components smaller than a given size. They are set to 0.

Definition at line 689 of file filters.cpp.

{
    MultidimArray<double> label;
    int imax;
    if (ZSIZE(I)==1)
        imax=labelImage2D(I, label, neighbourhood);
    else
        imax=labelImage3D(I, label);
    MultidimArray<int> nlabel(imax + 1);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(label)
    {
        int l=(int)DIRECT_MULTIDIM_ELEM(label,n);
        DIRECT_A1D_ELEM(nlabel,l)++;
    }
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(label)
    {
        int l=(int)DIRECT_MULTIDIM_ELEM(label,n);
        if (DIRECT_A1D_ELEM(nlabel,l)<size)
            DIRECT_MULTIDIM_ELEM(I,n)=0;
    }
}

rms()

template<typename T >

double rms	(	const MultidimArray< T > &	x,
		const MultidimArray< T > &	y,
		const MultidimArray< int > *	mask = nullptr,
		MultidimArray< double > *	Contributions = nullptr
	)

RMS nD

Definition at line 725 of file filters.h.

{
    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;
}

rotationalInvariantMoments()

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

Rotational invariant moments

The mask and the image are supposed to be of the same shape. If no mask is provided, the moments are computed on the whole image. The moments are measured with respect to the origin of the image.

These moments have been taken from http://www.cs.cf.ac.uk/Dave/Vision_lecture/node36.html (moments 1 to 5).

Definition at line 2900 of file filters.cpp.

{
    img.checkDimension(2);

    // Prepare some variables
    double m_11 = 0;
    double m_02 = 0;
    double m_20 = 0;
    double m_12 = 0;
    double m_21 = 0;
    double m_03 = 0;
    double m_30 = 0;
    double normalize_x = 2.0 / XSIZE(img);
    double normalize_y = 2.0 / YSIZE(img);

    // Compute low-level moments
    FOR_ALL_ELEMENTS_IN_ARRAY2D(img)
    {
        if (mask != nullptr)
            if (!(*mask)(i, j))
                continue;
        double val = img(i, j);
        double dx = j * normalize_x;
        double dy = i * normalize_y;
        double dx2 = dx * dx;
        double dx3 = dx2 * dx;
        double dy2 = dy * dy;
        double dy3 = dy2 * dy;
        // m_00+=val;
        m_11 += val * dx * dy;
        m_02 += val * dy2;
        m_20 += val * dx2;
        m_12 += val * dx * dy2;
        m_21 += val * dx2 * dy;
        m_03 += val * dy3;
        m_30 += val * dx3;
    }
    //m_11/=m_00; m_02/=m_00; m_20/=m_00;
    //m_12/=m_00; m_21/=m_00; m_03/=m_00; m_30/=m_00;

    // Compute high-level rotational invariant moments
    v_out.resize(5);
    v_out(0) = m_20 + m_02;
    v_out(1) = (m_20 - m_02) * (m_20 - m_02) + 4 * m_11 * m_11;
    v_out(2) = (m_30 - 3 * m_12) * (m_30 - 3 * m_12)
               + (3 * m_21 - m_03) * (3 * m_21 - m_03);
    v_out(3) = (m_30 + m_12) * (m_30 + m_12) + (m_21 + m_03) * (m_21 + m_03);
    v_out(4) =
        (m_30 - 3 * m_12) * (m_30 + m_12)
        * ((m_30 + m_12) * (m_30 + m_12)
           - 3 * (m_21 + m_03) * (m_21 + m_03))
        + (3 * m_21 - m_03) * (m_21 + m_03)
        * (3 * (m_30 + m_12) * (m_30 + m_12)
           - (m_21 + m_03) * (m_21 + m_03));
    /*
     v_out( 5)=(m_20+m_02)*(
     (m_30+m_12)*(m_30+m_12)
     -3*(m_21+m_03)*(m_21+m_03))
     +4*m_11*(m_30+m_12)*(m_03+m_21);
     v_out( 6)=(3*m_21-m_03)*(m_12+m_30)*(
     (m_30+m_12)*(m_30+m_12)
     -3*(m_21+m_03)*(m_21+m_03))
     -(m_30-3*m_12)*(m_12+m_03)*(
     3*(m_30+m_12)*(m_30+m_12)
     -(m_21+m_03)*(m_21+m_03));
     */
}

smoothingShah()

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
	)

Mumford-Shah smoothing

This function simultaneously smooths and segments an image using non-linear diffusion. Mumford-&-Shah's functional minimization algorithm is used to detect region boundaries and relax image smoothness constraints near these discontinuities. The functional minimized is:

E = W0*(f-d)*(f-d) (data matching)

• W1*(fx*fx + fy*fy)*(1-s)*(1-s) (1st deriv smooth)
• W2*(s*s) (edge strengths)
• W3*(sx*sx + sy*sy) (edge smoothness)

The program diffusion from KUIM (developed by J. Gaush, U. Kansas) was used as the "seed".

Paper: Teboul, et al. IEEE-Trans. on Image Proc. Vol. 7, 387-397.

Definition at line 2720 of file filters.cpp.

{

    img.checkDimension(2);

    typeCast(img, surface_strength);
    if (adjust_range)
        surface_strength.rangeAdjust(0, 1);
    edge_strength.resizeNoCopy(img);

    for (int k = 1; k <= RefinementLoops; k++)
    {
        // Initialize Edge Image.
        edge_strength.initZeros();

        double diffsurface = MAXFLOAT; // Reset surface difference
        for (int i = 0;
             ((i < OuterLoops) && OuterLoops)
             || ((diffsurface > SHAH_CONVERGENCE_THRESHOLD)
                 && !OuterLoops); i++)
        {

            /* std::cout << "Iteration ..." << i+1;*/
            /* Iteratively update surface estimate */
            for (int j = 0; j < InnerLoops; j++)
                diffsurface = Update_surface_Shah(img, surface_strength,
                                                  edge_strength, W);

            /* Iteratively update edge estimate */
            for (int j = 0; j < InnerLoops; j++)
                Update_edge_Shah(img, surface_strength, edge_strength, k, W);

            /* Calculate new functional energy */
            Shah_energy(img, surface_strength, edge_strength, k, W);
        }
    }
}

SmoothResize()

void SmoothResize	(	byte *	srcpic8,
		byte *	destpic8,
		size_t	swide,
		size_t	shigh,
		size_t	dwide,
		size_t	dhigh
	)

Resize of a 8-bit image

Definition at line 93 of file xvsmooth.cpp.

                                                                                                         {
    /* generic interface to Smooth and ColorDither code.
     given an 8-bit-per, swide * shigh image with colormap rmap,gmap,bmap,
     will generate a new 8-bit-per, dwide * dhigh image, which is dithered
     using colors found in rdmap, gdmap, bdmap arrays */

    /* returns ptr to a dwide*dhigh array of bytes, or NULL on failure */
    byte * picSmooth = Smooth(picSrc8, swide, shigh, dwide, dhigh);
    if (picSmooth) {
        DoColorDither(picSmooth, destpic8, dwide, dhigh);
        free(picSmooth);
    }
}

sort()

template<typename T >

void sort	(	T	a,
		T	b,
		T	c,
		MultidimArray< T > &	v
	)

Harmonic decomposition

Definition at line 798 of file filters.h.

{
    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;
    }
}

substractBackgroundPlane()

void substractBackgroundPlane

(

MultidimArray< double > &

I

)

Subtract background

The background is computed as the plane which best fits all density values, then this plane is subtracted from the image.

Definition at line 40 of file filters.cpp.

{

    I.checkDimension(2);

    Matrix2D<double> A(3, 3);
    Matrix1D<double> x(3);
    Matrix1D<double> b(3);

    // Solve the plane 'x'
    A.initZeros();
    b.initZeros();
    FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
    {
        A(0, 0) += j * j;
        A(0, 1) += j * i;
        A(0, 2) += j;
        A(1, 1) += i * i;
        A(1, 2) += i;
        A(2, 2) += 1;
        b(0) += j * A2D_ELEM(I, i, j);
        b(1) += i * A2D_ELEM(I, i, j);
        b(2) += A2D_ELEM(I, i, j);
    }
    A(1, 0) = A(0, 1);
    A(2, 0) = A(0, 2);
    A(2, 1) = A(1, 2);
    solve(A, b, x);

    // Now subtract the plane
    FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
    A2D_ELEM(I, i, j) -= x(0) * i + x(1) * j + x(2);
}

substractBackgroundRollingBall()

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

Subtract background

The background is computed as a rolling ball operation with a ball with this radius. The radius is typically Xdim/10.

This code has been implemented after the one of "Subtract background" in ImageJ.

Definition at line 75 of file filters.cpp.

{

    I.checkDimension(2);

    // Build the ball
    int arcTrimPer;
    int shrinkFactor;
    if (radius <= 10)
    {
        shrinkFactor = 1;
        arcTrimPer = 24; // trim 24% in x and y
    }
    else if (radius <= 30)
    {
        shrinkFactor = 2;
        arcTrimPer = 24; // trim 24% in x and y
    }
    else if (radius <= 100)
    {
        shrinkFactor = 4;
        arcTrimPer = 32; // trim 32% in x and y
    }
    else
    {
        shrinkFactor = 8;
        arcTrimPer = 40; // trim 40% in x and y
    }

    double smallballradius = radius / shrinkFactor;
    if (smallballradius < 1)
        smallballradius = 1;
    double r2 = smallballradius * smallballradius;
    int xtrim = (int) (arcTrimPer * smallballradius) / 100; // only use a patch of the rolling ball
    int halfWidth = ROUND(smallballradius - xtrim);
    int ballWidth = 2 * halfWidth + 1;
    MultidimArray<double> ball(ballWidth, ballWidth);
    ball.setXmippOrigin();
    FOR_ALL_ELEMENTS_IN_ARRAY2D(ball)
    {
        double temp = r2 - i * i - j * j;
        ball(i, j) = temp > 0. ? sqrt(temp) : 0.;
    }

    // Shrink the image: each point in the shrinked image is the
    // minimum of its neighbourhood
    int sXdim = (XSIZE(I) + shrinkFactor - 1) / shrinkFactor;
    int sYdim = (YSIZE(I) + shrinkFactor - 1) / shrinkFactor;
    MultidimArray<double> shrinkI(sYdim, sXdim);
    shrinkI.setXmippOrigin();
    for (int ySmall = 0; ySmall < sYdim; ySmall++)
    {
        for (int xSmall = 0; xSmall < sXdim; xSmall++)
        {
            double minVal = 1e38;
            for (int j = 0, y = shrinkFactor * ySmall;
                 j < shrinkFactor && y < (int)YSIZE(I); j++, y++)
                for (int k = 0, x = shrinkFactor * xSmall;
                     k < shrinkFactor && x < (int)XSIZE(I); k++, x++)
                {
                    double thispixel = DIRECT_A2D_ELEM(I,y,x);
                    if (thispixel < minVal)
                        minVal = thispixel;
                }
            DIRECT_A2D_ELEM(shrinkI,ySmall,xSmall) = minVal;
        }
    }

    // Now roll the ball
    radius = ballWidth / 2;
    MultidimArray<double> Irolled;
    Irolled.resizeNoCopy(shrinkI);
    Irolled.initConstant(-500);
    for (int yb = -radius; yb < (int)YSIZE(shrinkI) + radius; yb++)
    {
        // Limits of the ball
        int y0 = yb - radius;
        if (y0 < 0)
            y0 = 0;
        int y0b = y0 - yb + radius; //y coordinate in the ball corresponding to y0
        int yF = yb + radius;
        if (yF >= (int)YSIZE(shrinkI))
            yF = (int)YSIZE(shrinkI) - 1;

        for (int xb = -radius; xb < (int)XSIZE(shrinkI) + radius; xb++)
        {
            // Limits of the ball
            int x0 = xb - radius;
            if (x0 < 0)
                x0 = 0;
            int x0b = x0 - xb + radius;
            int xF = xb + radius;
            if (xF >= (int)XSIZE(shrinkI))
                xF = (int)XSIZE(shrinkI) - 1;

            double z = 1e38;
            for (int yp = y0, ybp = y0b; yp <= yF; yp++, ybp++)
                for (int xp = x0, xbp = x0b; xp <= xF; xp++, xbp++)
                {
                    double zReduced = DIRECT_A2D_ELEM(shrinkI,yp,xp)
                                      - DIRECT_A2D_ELEM(ball,ybp,xbp);
                    if (z > zReduced)
                        z = zReduced;
                }
            for (int yp = y0, ybp = y0b; yp <= yF; yp++, ybp++)
                for (int xp = x0, xbp = x0b; xp <= xF; xp++, xbp++)
                {
                    double zMin = z + DIRECT_A2D_ELEM(ball,ybp,xbp);
                    if (DIRECT_A2D_ELEM(Irolled,yp,xp) < zMin)
                        DIRECT_A2D_ELEM(Irolled,yp,xp) = zMin;
                }
        }
    }

    // Now rescale the background
    MultidimArray<double> bgEnlarged;
    scaleToSize(xmipp_transformation::LINEAR, bgEnlarged, Irolled, XSIZE(I), YSIZE(I));
    bgEnlarged.copyShape(I);
    I -= bgEnlarged;
}

tomographicDiffusion()

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

Tomographic diffusion

The direction of the tilt axis must be taken into account in the definition of the diffusion constants alpha.

The function returns the value of the regularization term.

Definition at line 2764 of file filters.cpp.

{
    V.checkDimension(3);

    double alphax = XX(alpha);
    double alphay = YY(alpha);
    double alphaz = ZZ(alpha);
    double diffx;
    double diffy;
    double diffz;

    // Compute regularization error
    double regError = 0;
    for (size_t z = 1; z < ZSIZE(V) - 1; z++)
        for (size_t y = 1; y < YSIZE(V) - 1; y++)
            for (size_t x = 1; x < XSIZE(V) - 1; x++)
            {
                diffx = DIRECT_A3D_ELEM(V,z,y,x+1) - DIRECT_A3D_ELEM(V,z,y,x-1);
                diffy = DIRECT_A3D_ELEM(V,z,y+1,x) - DIRECT_A3D_ELEM(V,z,y-1,x);
                diffz = DIRECT_A3D_ELEM(V,z+1,y,x) - DIRECT_A3D_ELEM(V,z-1,y,x);
                regError += sqrt(
                                alphax * diffx * diffx + alphay * diffy * diffy
                                + alphaz * diffz * diffz);
            }
    regError *= 0.5;

    // Compute the gradient of the regularization error
    MultidimArray<double> gradient;
    gradient.initZeros(V);
    for (size_t z = 2; z < ZSIZE(V) - 2; z++)
        for (size_t y = 2; y < YSIZE(V) - 2; y++)
            for (size_t x = 2; x < XSIZE(V) - 2; x++)
            {
                // First term
                double V000 = DIRECT_A3D_ELEM(V,z,y,x);
                double V_200 = DIRECT_A3D_ELEM(V,z,y,x-2);
                double V_110 = DIRECT_A3D_ELEM(V,z,y+1,x-1);
                double V_1_10 = DIRECT_A3D_ELEM(V,z,y-1,x-1);
                double V_101 = DIRECT_A3D_ELEM(V,z+1,y,x-1);
                double V_10_1 = DIRECT_A3D_ELEM(V,z-1,y,x-1);
                diffx = V000 - V_200;
                diffy = V_110 - V_1_10;
                diffz = V_101 - V_10_1;
                double t1 = diffx
                            / sqrt(
                                alphax * diffx * diffx + alphay * diffy * diffy
                                + alphaz * diffz * diffz);

                // Second term
                double V200 = DIRECT_A3D_ELEM(V,z,y,x+2);
                double V110 = DIRECT_A3D_ELEM(V,z,y+1,x+1);
                double V1_10 = DIRECT_A3D_ELEM(V,z,y-1,x+1);
                double V101 = DIRECT_A3D_ELEM(V,z+1,y,x+1);
                double V10_1 = DIRECT_A3D_ELEM(V,z-1,y,x+1);
                diffx = V200 - V000;
                diffy = V110 - V1_10;
                diffz = V101 - V10_1;
                double t2 = diffx
                            / sqrt(
                                alphax * diffx * diffx + alphay * diffy * diffy
                                + alphaz * diffz * diffz);

                // Third term
                double V0_20 = DIRECT_A3D_ELEM(V,z,y-2,x);
                double V0_11 = DIRECT_A3D_ELEM(V,z+1,y-1,x);
                double V0_1_1 = DIRECT_A3D_ELEM(V,z-1,y-1,x);
                diffx = V1_10 - V_1_10;
                diffy = V000 - V0_20;
                diffz = V0_11 - V0_1_1;
                double t3 = diffy
                            / sqrt(
                                alphax * diffx * diffx + alphay * diffy * diffy
                                + alphaz * diffz * diffz);

                // Fourth term
                double V020 = DIRECT_A3D_ELEM(V,z,y+2,x);
                double V011 = DIRECT_A3D_ELEM(V,z+1,y+1,x);
                double V01_1 = DIRECT_A3D_ELEM(V,z-1,y+1,x);
                diffx = V110 - V_110;
                diffy = V020 - V000;
                diffz = V011 - V01_1;
                double t4 = diffy
                            / sqrt(
                                alphax * diffx * diffx + alphay * diffy * diffy
                                + alphaz * diffz * diffz);

                // Fifth term
                double V00_2 = DIRECT_A3D_ELEM(V,z-2,y,x);
                diffx = V10_1 - V_10_1;
                diffy = V01_1 - V0_1_1;
                diffz = V000 - V00_2;
                double t5 = diffz
                            / sqrt(
                                alphax * diffx * diffx + alphay * diffy * diffy
                                + alphaz * diffz * diffz);

                // Sixth term
                double V002 = DIRECT_A3D_ELEM(V,z+2,y,x);
                diffx = V101 - V_101;
                diffy = V011 - V0_11;
                diffz = V002 - V000;
                double t6 = diffz
                            / sqrt(
                                alphax * diffx * diffx + alphay * diffy * diffy
                                + alphaz * diffz * diffz);

                // Compute gradient
                DIRECT_A3D_ELEM(gradient,z,y,x) = 0.5
                                                  * (alphax * (t1 - t2) + alphay * (t3 - t4)
                                                     + alphaz * (t5 - t6));
            }
#ifdef DEBUG
    VolumeXmipp save;
    save()=V;
    save.write("PPPvolume.vol");
    save()=gradient;
    save.write("PPPgradient.vol");
    std::cout << "Press any key\n";
    char c;
    std::cin >> c;
#endif

    // Update volume
    for (size_t z = 2; z < ZSIZE(V) - 2; z++)
        for (size_t y = 2; y < YSIZE(V) - 2; y++)
            for (size_t x = 2; x < XSIZE(V) - 2; x++)
                DIRECT_A3D_ELEM(V,z,y,x) -= lambda
                                            * DIRECT_A3D_ELEM(gradient,z,y,x);

    // Finish
    return regError;
}

varianceFilter()

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

Applays a variance filter to an image

If relative=True, the filter is normalized to the mean (coeficient of variation)

Definition at line 775 of file filters.cpp.

{
    int kernelSize_2 = kernelSize/2;
    MultidimArray<double> kernel;
    kernel.resize(kernelSize,kernelSize);
    kernel.setXmippOrigin();

    // std::cout << " Creating the variance matrix " << std::endl;
    MultidimArray<double> mVar(YSIZE(I),XSIZE(I));
    mVar.setXmippOrigin();
    double stdKernel;
    double avgKernel;
    double min_val;
    double max_val;
    int x0;
    int y0;
    int xF;
    int yF;

    // I.computeStats(avgImg, stdImg, min_im, max_im);    
    
    for (int i=kernelSize_2; i<=(int)YSIZE(I)-kernelSize_2; i+=kernelSize_2)
        for (int j=kernelSize_2; j<=(int)XSIZE(I)-kernelSize_2; j+=kernelSize_2)
            {
                x0 = j-kernelSize_2;
                y0 = i-kernelSize_2;
                xF = j+kernelSize_2-1;
                yF = i+kernelSize_2-1;

                if (x0 < 0)
                    x0 = 0;
                if (xF > XSIZE(I))
                    xF = XSIZE(I);
                if (y0 < 0)
                    y0 = 0;
                if (yF > YSIZE(I))
                    yF = YSIZE(I);

                I.window(kernel, y0, x0, yF, xF);
                kernel.computeStats(avgKernel, stdKernel, min_val, max_val);
                
                DIRECT_A2D_ELEM(mVar, i, j) = stdKernel;
            }

    // filtering to fill the matrices (convolving with a Gaussian)
    FourierFilter filter;
    filter.FilterShape = REALGAUSSIAN;
    filter.FilterBand = LOWPASS;
    filter.w1 = kernelSize_2;
    filter.applyMaskSpace(mVar);

    if (relative) // normalize to the mean of the variance
    {    
        double avgVar;
        double stdVar;
        double minVar;
        double maxVar;
        mVar.computeStats(avgVar, stdVar, minVar, maxVar);
        mVar = mVar/avgVar;
    }

    I = mVar;

    // filling the borders with the nearest variance value
    // the corners only are filled once (with Ycoord)
    for (int i=0; i<YSIZE(I); ++i)
        for (int j=0; j<kernelSize; j++)
        {
            if (i<kernelSize)
                DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, kernelSize, kernelSize);
            else if (i>YSIZE(I)-kernelSize)
                DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, YSIZE(I)-kernelSize, kernelSize);
            else
                DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, i, kernelSize);
        }
    for (int i=0; i<YSIZE(I); ++i)
        for (int j=XSIZE(I)-kernelSize; j<XSIZE(I); j++)
        {
            if (i<kernelSize)
                DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, kernelSize, XSIZE(I)-kernelSize);
            else if (i>YSIZE(I)-kernelSize)
                DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, YSIZE(I)-kernelSize, XSIZE(I)-kernelSize);
            else
                DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, i, XSIZE(I)-kernelSize);
        }
    for (int j=kernelSize; j<XSIZE(I)-kernelSize; ++j)
        for (int i=0; i<kernelSize; i++)
            DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, kernelSize, j);
    for (int j=kernelSize; j<XSIZE(I)-kernelSize; ++j)
        for (int i=YSIZE(I)-kernelSize; i<YSIZE(I); i++)
            DIRECT_A2D_ELEM(I, i, j) = DIRECT_A2D_ELEM(mVar, YSIZE(I)-kernelSize, j);
    
}