Normalization of images and volumes

Collaboration diagram for Normalization of images and volumes:

Classes
class	ProgNormalize

Image normalization procedures
This functions implement the normalization of a single image. They should be called with all images in the corresponding SelFile. In the following documentation m(x) is the mean of x, v(x) is the variance, bg(x) is its background. The original image is X and is supposed to be related to each single projection by I=a*(X+n)+b where a and b are different for every projection. Noise is assumed to follow a gaussian distribution N(0,sqrt(v(n))) In general the background mask is used only to compute statistics, while the mask is the one which is really applied to the image. Supply NULL if you don't want any mask to be applied When b is used it is measured as the mean in the background, while a*sqrt(v(n)) is the standard deviation in the same area.
void	normalize_OldXmipp (MultidimArray< double > &I)
void	normalize_Near_OldXmipp (MultidimArray< double > &I, const MultidimArray< int > &bg_mask)
void	normalize_OldXmipp_decomposition (MultidimArray< double > &I, const MultidimArray< int > &bg_mask, const MultidimArray< double > *mask=nullptr)
void	normalize_tomography (MultidimArray< double > &I, double tilt, double &mui, double &sigmai, bool tiltMask, bool tomography0=false, double mu0=0, double sigma0=1)
void	normalize_Michael (MultidimArray< double > &I, const MultidimArray< int > &bg_mask)
void	normalize_NewXmipp (MultidimArray< double > &I, const MultidimArray< int > &bg_mask)
void	normalize_NewXmipp2 (MultidimArray< double > &I, const MultidimArray< int > &bg_mask)
void	normalize_Robust (MultidimArray< double > &I, const MultidimArray< int > &bg_mask, bool clip)
void	normalize_ramp (MultidimArray< double > &I, MultidimArray< int > *bg_mask=nullptr)
void	normalize_remove_neighbours (MultidimArray< double > &I, const MultidimArray< int > &bg_mask, const double &threshold)

Detailed Description

Function Documentation

normalize_Michael()

void normalize_Michael	(	MultidimArray< double > &	I,
		const MultidimArray< int > &	bg_mask
	)

Michael's normalization.

Formula:

I'=(I-b)/b

Properties:

m(I')=0                             m(bg(I'))=-a*m(x)/b
v(I')=a^2*(v(X)+v(n))/b^2           v(bg(I'))=a^2*v(n)/b^2

Comments: it's not bad but positivity constraints cannot be imposed and the statistical properties are not so good.

Definition at line 230 of file normalize.cpp.

{
    double avg;
    double stddev;
    double min=0.;
    double max;
    double avgbg;
    double stddevbg;
    double minbg;
    double maxbg;
    I.computeStats(avg, stddev, min, max);
    computeStats_within_binary_mask(bg_mask, I, minbg, maxbg, avgbg,
                                    stddevbg);
    if (avgbg > 0)
    {
        I -= avgbg;
        I /= avgbg;
    }
    else
    { // To avoid the contrast inversion
        I -= (avgbg - min);
        I /= (avgbg - min);
    }
}

normalize_Near_OldXmipp()

void normalize_Near_OldXmipp	(	MultidimArray< double > &	I,
		const MultidimArray< int > &	bg_mask
	)

Near_OldXmipp normalization.

Formula:

I'=(I-m(I))/sqrt(v(bg))

Properties:

m(I')=0                             m(bg(I'))=-m(x)/sqrt(v(n))
v(I')=(v(X)+v(n))/v(n)              v(bg(I'))=1

Comments: it's not bad but positivity constraints cannot be imposed

Definition at line 43 of file normalize.cpp.

{
    double avg=0.;
    double stddev;
    double min;
    double max;
    double avgbg;
    double stddevbg;
    double minbg;
    double maxbg;
    I.computeStats(avg, stddev, min, max);
    computeStats_within_binary_mask(bg_mask, I, minbg, maxbg, avgbg,
                                    stddevbg);
    I -= avg;
    I /= stddevbg;
}

normalize_NewXmipp()

void normalize_NewXmipp	(	MultidimArray< double > &	I,
		const MultidimArray< int > &	bg_mask
	)

NewXmipp's normalization.

Formula:

I'=(I-b)/a*sqrt(v(n))
// I''=(I'>0)? I':0
// I'''=I''-a*sqrt(v(n)/2*PI)
I''''=I'''*mask

Properties:

m(I')=m(X)/sqrt(v(n))               m(bg(I'))=0
v(I')=(v(X)+v(n))/v(n)              v(bg(I'))=1

Comments: In general, we cannot assure that mass projects into positive numbers, so the "denoising" capability directly on the images is disabled. However, a positivity constraint can be applied on the 3D volume.

Definition at line 255 of file normalize.cpp.

{
    double avgbg;
    double stddevbg;
    I.computeAvgStdev_within_binary_mask(bg_mask, avgbg, stddevbg);
    double istddevbg=1.0/stddevbg;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I)
    DIRECT_MULTIDIM_ELEM(I,n)=(DIRECT_MULTIDIM_ELEM(I,n)-avgbg)*istddevbg;
}

normalize_NewXmipp2()

void normalize_NewXmipp2	(	MultidimArray< double > &	I,
		const MultidimArray< int > &	bg_mask
	)

NewXmipp 2's normalization.

Formula:

I'=(I-m(bg))/(m(I)-m(bg))

Properties:

m(I')=m(X)/sqrt(v(n))               m(bg(I'))=0
v(I')=(v(X)+v(n))/v(n)              v(bg(I'))=1

Comments: In general, we cannot assure that mass projects into positive numbers, so the "denoising" capability directly on the images is disabled. However, a positivity constraint can be applied on the 3D volume.

Definition at line 315 of file normalize.cpp.

{
    double avg=0;
    double stddev;
    double min;
    double max;
    double avgbg=0;
    double stddevbg;
    double minbg;
    double maxbg;
    I.computeStats(avg, stddev, min, max);
    computeStats_within_binary_mask(bg_mask, I, minbg, maxbg, avgbg,
                                    stddevbg);
    double K=1.0/fabs(avg - avgbg);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I)
    DIRECT_MULTIDIM_ELEM(I,n)=(DIRECT_MULTIDIM_ELEM(I,n)-avgbg)*K;
}

normalize_OldXmipp()

void normalize_OldXmipp

(

MultidimArray< double > &

I

)

OldXmipp normalization.

Formula:

I'=(I-m(I))/sqrt(v(I))

Properties:

m(I')=0                             m(bg(I'))=-m(x)/sqrt((v(X)+v(n)))
v(I')=1                             v(bg(I'))=v(n)/(v(X)+v(n))

Comments: it's not bad but positivity constraints cannot be imposed

Definition at line 33 of file normalize.cpp.

{
    double mean;
    double std;
    I.computeAvgStdev(mean,std);
    double istd=1.0/std;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I)
    DIRECT_MULTIDIM_ELEM(I,n)=(DIRECT_MULTIDIM_ELEM(I,n)-mean)*istd;
}

normalize_OldXmipp_decomposition()

void normalize_OldXmipp_decomposition	(	MultidimArray< double > &	I,
		const MultidimArray< int > &	bg_mask,
		const MultidimArray< double > *	mask = nullptr
	)

OldXmipp decomposition.

Formula:

I'=(I-b)/a*sqrt(v(n))
I''=I'*mask
I'''=(I''-m(I''))/sqrt(v(I''))

Properties:

m(I')=m(X)/sqrt(v(n))               m(bg(I'))=0
v(I')=(v(X)+v(n))/v(n)              v(bg(I'))=1

Comments: it's not bad but positivity constraints cannot be imposed. If no mask is applied, then this formula is a beautiful decomposition of the OldXmipp method in two steps.

Definition at line 60 of file normalize.cpp.

{
    double avgbg;
    double stddevbg;
    double minbg;
    double maxbg;
    computeStats_within_binary_mask(bg_mask, I, minbg, maxbg, avgbg,
                                    stddevbg);
    I -= avgbg;
    I /= stddevbg;
    if (mask != nullptr)
        I *= *mask;
    normalize_OldXmipp(I);
}

normalize_ramp()

void normalize_ramp	(	MultidimArray< double > &	I,
		MultidimArray< int > *	bg_mask = nullptr
	)

Removal of inclined background densities (ramps).

fitting of a least squares plane through the pixels in the bg_mask, then subtraction of the plane, and division by the standard deviation of the pixels in the bg_mask

Definition at line 333 of file normalize.cpp.

{
    int Npoints=0;              // # points in mask.
    double pA;
    double pB;
    double pC;          // Least squares coefficients.

    // Only 2D ramps implemented
    I.checkDimension(2);

    // Check if mask is NULL.
    if (bg_mask == nullptr)
    {
        Npoints = I.xdim*I.ydim;
    }
    // Check if mask is not full of 0's.
    else
    {
        Npoints=(int)bg_mask->sum();
    }

    // Fit a least squares plane through the background pixels
    I.setXmippOrigin();
    if (bg_mask == nullptr)
    {
        least_squares_plane_fit_All_Points(I, pA, pB, pC);
    }
    else
    {
        int idx=0;
        bg_mask->setXmippOrigin();
        auto *allpoints=new FitPoint[Npoints];

        FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
        {
            if (A2D_ELEM( *bg_mask, i, j))
            {
                FitPoint &p=allpoints[idx++];
                p.x = j;
                p.y = i;
                p.z = A2D_ELEM(I, i, j);
                p.w = 1.;
            }
        }

        least_squares_plane_fit(allpoints, Npoints, pA, pB, pC);
        delete [] allpoints;
    }

    // Subtract the plane from the image and compute stddev within mask
    double sum1 = 0;
    double sum2 = 0;
    if (bg_mask == nullptr)
    {
        double *ref;
        for (int i=STARTINGY(I); i<=FINISHINGY(I); i++)
        {
            double aux=pB * i + pC;
            ref = &A2D_ELEM(I, i, STARTINGX(I));
            for (int j=STARTINGX(I); j<=FINISHINGX(I); j++)
            {
                (*ref) -= pA * j + aux;
                sum1 += (*ref);
                sum2 += (*ref)*(*ref);
                ref++;
            }
        }
    }
    else
    {
        for (int i=STARTINGY(I); i<=FINISHINGY(I); i++)
        {
            double aux=pB * i + pC;
            for (int j=STARTINGX(I); j<=FINISHINGX(I); j++)
            {
                A2D_ELEM(I, i, j) -= pA * j + aux;
                if (A2D_ELEM( *bg_mask, i, j))
                {
                    double aux=A2D_ELEM(I,i,j);
                    sum1 += aux;
                    sum2 += aux*aux;
                }
            }
        }
    }

    double avgbg  = sum1 / Npoints;
    double stddevbg = sqrt(fabs(sum2 / Npoints - avgbg * avgbg) * Npoints / (Npoints - 1));
    if (stddevbg>1e-6)
    {
        I *= 1.0/stddevbg;
    }
}

normalize_remove_neighbours()

void normalize_remove_neighbours	(	MultidimArray< double > &	I,
		const MultidimArray< int > &	bg_mask,
		const double &	threshold
	)

Removal of neighbouring particles.

....

Definition at line 427 of file normalize.cpp.

{
    double             pA;
    double             pB;
    double             pC;
    double             avgbg;
    double             stddevbg;
    double             minbg;
    double             maxbg;
    double             aux;
    double             newstddev;
    double             sum1 = 0.;
    double             sum2 = 0;
    int                N = 0;

    // Only implemented for 2D arrays
    I.checkDimension(2);

    // Fit a least squares plane through the background pixels
    auto Npoints=(int)bg_mask.sum();
    auto *allpoints=new FitPoint[Npoints];
    I.setXmippOrigin();

    // Get initial statistics
    computeStats_within_binary_mask(bg_mask, I, minbg, maxbg, avgbg,stddevbg);

    // Fit plane through those pixels within +/- threshold*sigma
    int idx=0;
    FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
    {
        if (A2D_ELEM(bg_mask, i, j))
        {
            if ( fabs(avgbg - A2D_ELEM(I, i, j)) < threshold * stddevbg)
            {
                FitPoint &p=allpoints[idx++];
                p.x = j;
                p.y = i;
                p.z = A2D_ELEM(I, i, j);
                p.w = 1.;
            }
        }
    }
    least_squares_plane_fit(allpoints, Npoints, pA, pB, pC);
    delete []allpoints;

    // Subtract the plane from the image
    FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
    {
        A2D_ELEM(I, i, j) -= pA * j + pB * i + pC;
    }

    // Get std.dev. of the background pixels within +/- threshold*sigma
    FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
    {
        if (A2D_ELEM(bg_mask, i, j))
        {
            if ( ABS(A2D_ELEM(I, i, j)) < threshold * stddevbg)
            {
                N++;
                sum1 +=  (double) A2D_ELEM(I, i, j);
                sum2 += ((double) A2D_ELEM(I, i, j)) *
                        ((double) A2D_ELEM(I, i, j));
            }
        }
    }
    // average and standard deviation
    if (0 == N) {
        REPORT_ERROR(ERR_NUMERICAL, "N is zero (0), which would lead to division by zero");
    }
    aux = sum1 / (double) N;
    newstddev = sqrt(fabs(sum2 / N - aux*aux) * N / (N - 1));

    // Replace pixels outside +/- threshold*sigma by samples from
    // a gaussian with avg-plane and newstddev
    FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
    {
        if (A2D_ELEM(bg_mask, i, j))
        {
            if ( fabs(A2D_ELEM(I, i, j)) > threshold * stddevbg)
            {
                // get local average
                aux = pA * j + pB * i + pC;
                A2D_ELEM(I, i, j)=rnd_gaus(aux, newstddev );
            }
        }
    }

    // Divide the entire image by the new background
    I /= newstddev;
}

normalize_Robust()

void normalize_Robust	(	MultidimArray< double > &	I,
		const MultidimArray< int > &	bg_mask,
		bool	clip
	)

Definition at line 265 of file normalize.cpp.

{
    std::vector<double> voxel_vector;
    SPEED_UP_temps;

    if (bg_mask.computeMax() == 0)
    {
        Image<double> mask;
        mask() = I;
        static_cast<void>(EntropySegmentation(mask()));
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mask())
        {
            if (DIRECT_MULTIDIM_ELEM(mask(), n) == 0)
                DIRECT_MULTIDIM_ELEM(bg_mask,n) = 1;
        }
    }
    
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(bg_mask)
    {
        if (DIRECT_MULTIDIM_ELEM(bg_mask, n) == 0)
            voxel_vector.push_back(DIRECT_MULTIDIM_ELEM(I,n));
    }

    std::sort(voxel_vector.begin(), voxel_vector.end());

    double medianBg;
    double p99;
    double ip99;
    int idx;
    I.computeMedian_within_binary_mask(bg_mask, medianBg);
    idx = (int)(voxel_vector.size() * 0.99);
    p99 = voxel_vector[idx];
    ip99 = 1 / p99;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I)
        DIRECT_MULTIDIM_ELEM(I,n)=(DIRECT_MULTIDIM_ELEM(I,n) - medianBg) * ip99;

    if (clip)
    {
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I)
            if (DIRECT_MULTIDIM_ELEM(I,n) > 1.3284)
            {
                DIRECT_MULTIDIM_ELEM(I,n) = 1.3284;
            }
            else if (DIRECT_MULTIDIM_ELEM(I,n) < -1.3284)
            {
                DIRECT_MULTIDIM_ELEM(I,n) = -1.3284;
            }           
    }
}

normalize_tomography()

void normalize_tomography	(	MultidimArray< double > &	I,
		double	tilt,
		double &	mui,
		double &	sigmai,
		bool	tiltMask,
		bool	tomography0 = false,
		double	mu0 = 0,
		double	sigma0 = 1
	)

Tomography normalization.

This is similar to the OldXmipp normalization, but the mean and standard deviation of the images are computed only within a region determined by the tilt angle. Formula for tomography:

I'=(I-m(I))/(sqrt(v(I))*cos^2(tilt))

Formula for tomography0:

I'=(I/cos(tilt)-m0)/(sqrt(v(I))*cos(tilt))

The estimated mean of the image and the local variance are returned in sigmai and mui.

Definition at line 77 of file normalize.cpp.

{
    // Tilt series are always 2D
    I.checkDimension(2);

    const int L=2;

    // Build a mask using the tilt angle
    I.setXmippOrigin();
    MultidimArray<int> mask;
    mask.initZeros(I);
    int Xdimtilt=(int)XMIPP_MIN(FLOOR(0.5*(XSIZE(I)*cos(DEG2RAD(tilt)))),
                           0.5*(XSIZE(I)-(2*L+1)));
    int N=0;
    for (int i=STARTINGY(I); i<=FINISHINGY(I); i++)
        for (int j=-Xdimtilt; j<=Xdimtilt;j++)
        {
            A2D_ELEM(mask,i,j)=1;
            N++;
        }

    // Estimate the local variance
    MultidimArray<double> localVariance;
    localVariance.initZeros(I);
    double meanVariance=0;
    FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
    {
        // Center a mask of size 5x5 and estimate the variance within the mask
        double meanPiece=0;
        double variancePiece=0;
        double Npiece=0;
        for (int ii=i-L; ii<=i+L; ii++)
        {
            if (INSIDEY(I,ii))
                for (int jj=j-L; jj<=j+L; jj++)
                {
                    if (INSIDEX(I,jj))
                    {
                        double pixval=A2D_ELEM(I,ii,jj);
                        meanPiece+=pixval;
                        variancePiece+=pixval*pixval;
                        ++Npiece;
                    }
                }
        }
        meanPiece/=Npiece;
        variancePiece=variancePiece/(Npiece-1)-
                      Npiece/(Npiece-1)*meanPiece*meanPiece;
        A2D_ELEM(localVariance,i,j)=variancePiece;
        meanVariance+=variancePiece;
    }
    if (0 == N) {
        REPORT_ERROR(ERR_NUMERICAL, "N is zero (0), which would lead to division by zero");
    }
    meanVariance*=1.0/N;

    // Test the hypothesis that the variance in this piece is
    // the same as the variance in the whole image
    double iFu=1/icdf_FSnedecor(4*L*L+4*L,N-1,0.975);
    double iFl=1/icdf_FSnedecor(4*L*L+4*L,N-1,0.025);
    double iMeanVariance=1.0/meanVariance;
    FOR_ALL_ELEMENTS_IN_ARRAY2D(localVariance)
    {
        if (A2D_ELEM(localVariance,i,j)==0)
            A2D_ELEM(mask,i,j)=-2;
        if (A2D_ELEM(mask,i,j)==1)
        {
            double ratio=A2D_ELEM(localVariance,i,j)*iMeanVariance;
            double thl=ratio*iFu;
            double thu=ratio*iFl;
            if (thl>1 || thu<1)
                A2D_ELEM(mask,i,j)=-1;
        }
    }
#ifdef DEBUG
    Image<double> save;
    save()=I;
    save.write("PPP.xmp");
    save()=localVariance;
    save.write("PPPLocalVariance.xmp");
    Image<int> savemask;
    savemask()=mask;
    savemask.write("PPPmask.xmp");
#endif

    // Compute the statistics again in the reduced mask
    double avg;
    double stddev;
    double min;
    double max;
    computeStats_within_binary_mask(mask, I, min, max, avg, stddev);
    double cosTilt=cos(DEG2RAD(tilt));
    double iCosTilt=1.0/cosTilt;
    if (tomography0)
    {
        double iadjustedStddev=1.0/(sigma0*cosTilt);
        FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
        switch (A2D_ELEM(mask,i,j))
        {
        case -2:
            A2D_ELEM(I,i,j)=0;
            break;
        case 0:
            if (tiltMask)
                A2D_ELEM(I,i,j)=0;
            else
                A2D_ELEM(I,i,j)=(A2D_ELEM(I,i,j)*iCosTilt-mu0)*iadjustedStddev;
            break;
        case -1:
        case 1:
            A2D_ELEM(I,i,j)=(A2D_ELEM(I,i,j)*iCosTilt-mu0)*iadjustedStddev;
            break;
        }
    }
    else
    {
        double adjustedStddev=stddev*cosTilt;
        double iAdjustedStddev=1.0/adjustedStddev;
        FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
        switch (A2D_ELEM(mask,i,j))
        {
        case -2:
            A2D_ELEM(I,i,j)=0;
            break;
        case 0:
            if (tiltMask)
                A2D_ELEM(I,i,j)=0;
            else
                A2D_ELEM(I,i,j)=(A2D_ELEM(I,i,j)-avg)*iAdjustedStddev;
            break;
        case -1:
        case 1:
            A2D_ELEM(I,i,j)=(A2D_ELEM(I,i,j)-avg)*iAdjustedStddev;
            break;
        }
    }

    // Prepare values for returning
    mui=avg;
    sigmai=sqrt(meanVariance);
#ifdef DEBUG

    save()=I;
    save.write("PPPafter.xmp");
    std::cout << "Press any key\n";
    char c;
    std::cin >> c;
#endif
}