#include <queue>
#include <list>
#include "data/fourier_filter.h"
#include "core/histogram.h"
#include "core/metadata_extension.h"
#include "core/multidim_array.h"
#include "core/transformations.h"
#include "core/xmipp_program.h"
#include "filters.h"
#include "mask.h"
#include "morphology.h"
#include "wavelet.h"

Include dependency graph for filters.cpp:

Classes
struct	Coordinate2D

struct	Coordinate3D

Macros
#define	CHECK_POINT(i, j)

#define	CHECK_POINT_3D(k, i, j)

#define	CHECK_POINT_3D(k, i, j)

#define	CHECK_POINT_DIST(i, j, proposedDistance)

Functions
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, float filling_colour, bool less, int neighbourhood)

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

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

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

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

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

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

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

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

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

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

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

double	OtsuSegmentation (MultidimArray< double > &V)

double	EntropySegmentation (MultidimArray< double > &V)

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

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

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

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

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

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

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

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

template float	bestShift (MultidimArray< float > &, float &, float &, const MultidimArray< int > *, int)

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

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

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

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

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

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

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

double	bestShiftRealSpace (const MultidimArray< double > &I1, MultidimArray< double > &I2, double &shiftX, double &shiftY, const MultidimArray< int > *mask, int maxShift, double shiftStep)

void	computeAlignmentTransforms (const MultidimArray< double > &I, AlignmentTransforms &ITransforms, AlignmentAux &aux, CorrelationAux &aux2)

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

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

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

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

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

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

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

double	fastBestRotation (const MultidimArray< double > &IrefCyl, const MultidimArray< double > &Icyl, CorrelationAux &aux, VolumeAlignmentAux &aux2, double deltaAng)

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

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

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

void	fastBestRotation (const MultidimArray< double > &IrefCylZ, const MultidimArray< double > &IrefCylY, const MultidimArray< double > &IrefCylX, const MultidimArray< double > &I, const MultidimArray< double > &Icurrent, MultidimArray< double > &Ifinal, char axis, Matrix2D< double > &R, CorrelationAux &aux, VolumeAlignmentAux &aux2)

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

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

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

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

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

double	Shah_energy (const MultidimArray< double > &img, const MultidimArray< double > &surface_strength, const MultidimArray< double > &edge_strength, double K, const Matrix1D< double > &W)

double	Update_surface_Shah (MultidimArray< double > &img, MultidimArray< double > &surface_strength, MultidimArray< double > &edge_strength, const Matrix1D< double > &W)

double	Update_edge_Shah (MultidimArray< double > &img, MultidimArray< double > &surface_strength, MultidimArray< double > &edge_strength, double K, const Matrix1D< double > &W)

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)

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)

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, bool limitShift)

void	forcePositive (MultidimArray< double > &V)

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

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

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	matlab_filter (Matrix1D< double > &B, Matrix1D< double > &A, const MultidimArray< double > &X, MultidimArray< double > &Y, MultidimArray< double > &Z)

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

Variables
constexpr double	SHIFT_THRESHOLD = 0.95

constexpr float	ROTATE_THRESHOLD = 1.0

constexpr double	INITIAL_SHIFT_THRESHOLD = SHIFT_THRESHOLD + 1.0

constexpr float	INITIAL_ROTATE_THRESHOLD = ROTATE_THRESHOLD + 1.0

constexpr double	SHAH_CONVERGENCE_THRESHOLD = 0.0001

Macro Definition Documentation

◆ CHECK_POINT

#define CHECK_POINT	(	i,
		j
	)

Value:

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

◆ CHECK_POINT_3D [1/2]

#define CHECK_POINT_3D	(	k,
		i,
		j
	)

Value:

{\
     int I=i; \
        int J=j; \
        int K=k; \
  if (INSIDEXYZ(V_out,J,I,K))  { \
   double voxel=A3D_ELEM(V_out,K,I,J); \
   if (voxel!=filling_colour) \
    if ((less && voxel < stop_colour)|| \
     (!less &&voxel > stop_colour)) { \
     coord.ii=I; \
     coord.jj=J; \
     coord.kk=K; \
     A3D_ELEM (V_out,coord.kk,coord.ii,coord.jj)=filling_colour; \
     iNeighbours.push_front(coord); \
    } \
  }\
 }

◆ CHECK_POINT_3D [2/2]

#define CHECK_POINT_3D	(	k,
		i,
		j
	)

Value:

{\
        int I=i; \
        int J=j; \
        int K=k; \
  if (INSIDEXYZ(V_out,J,I,K))  { \
   if (A3D_ELEM(V_in,K,I,J)==bgColor && A3D_ELEM(V_out,K,I,J)==1) { \
     coord.ii=I; \
     coord.jj=J; \
     coord.kk=K; \
     A3D_ELEM (V_out,coord.kk,coord.ii,coord.jj)=filling_value; \
     iNeighbours.push_front(coord); \
    } \
  }\
 }\

◆ CHECK_POINT_DIST

#define CHECK_POINT_DIST	(	i,
		j,
		proposedDistance
	)

Value:

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

Function Documentation

◆ alignImages()

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

Definition at line 2047 of file filters.cpp.

 {
     I.checkDimension(2);
 
     aux.ARS.initIdentity(3);
     aux.ASR.initIdentity(3);
     aux.IauxSR = I;
     aux.IauxRS = I;
 
     aux.rotationalCorr.resize(2 * IrefTransforms.polarFourierI.getSampleNoOuterRing() - 1);
     aux3.local_transformer.setReal(aux.rotationalCorr);
 
     double shiftXSR=INITIAL_SHIFT_THRESHOLD;
     double shiftYSR=INITIAL_SHIFT_THRESHOLD;
     double bestRotSR=INITIAL_ROTATE_THRESHOLD;
     double shiftXRS=INITIAL_SHIFT_THRESHOLD;
     double shiftYRS=INITIAL_SHIFT_THRESHOLD;
     double bestRotRS=INITIAL_ROTATE_THRESHOLD;
 
     // Align the image with the reference
     for (int i = 0; i < 3; i++)
     {
         // Shift then rotate
         if (((shiftXSR > SHIFT_THRESHOLD) || (shiftXSR < (-SHIFT_THRESHOLD))) ||
             ((shiftYSR > SHIFT_THRESHOLD) || (shiftYSR < (-SHIFT_THRESHOLD))))
         {
             bestNonwrappingShift(Iref, IrefTransforms.FFTI, aux.IauxSR, shiftXSR, shiftYSR, aux2);
             MAT_ELEM(aux.ASR,0,2) += shiftXSR;
             MAT_ELEM(aux.ASR,1,2) += shiftYSR;
             applyGeometry(xmipp_transformation::LINEAR, aux.IauxSR, I, aux.ASR, xmipp_transformation::IS_NOT_INV, wrap);
         }
 
         if (bestRotSR > ROTATE_THRESHOLD)
         {
             polarFourierTransform<true>(aux.IauxSR, aux.polarFourierI, true,
                                             XSIZE(Iref) / 5, XSIZE(Iref) / 2, aux.plans, 1);
 
             bestRotSR = best_rotation(IrefTransforms.polarFourierI, aux.polarFourierI, aux3);
             rotation2DMatrix(bestRotSR, aux.R);
             aux.ASR = aux.R * aux.ASR;
             applyGeometry(xmipp_transformation::LINEAR, aux.IauxSR, I, aux.ASR, xmipp_transformation::IS_NOT_INV, wrap);
         }
 
         // Rotate then shift
         if (bestRotRS > ROTATE_THRESHOLD)
         {
             polarFourierTransform<true>(aux.IauxRS, aux.polarFourierI, true,
                                             XSIZE(Iref) / 5, XSIZE(Iref) / 2, aux.plans, 1);
             bestRotRS = best_rotation(IrefTransforms.polarFourierI, aux.polarFourierI, aux3);
             rotation2DMatrix(bestRotRS, aux.R);
             aux.ARS = aux.R * aux.ARS;
             applyGeometry(xmipp_transformation::LINEAR, aux.IauxRS, I, aux.ARS, xmipp_transformation::IS_NOT_INV, wrap);
         }
 
         if (((shiftXRS > SHIFT_THRESHOLD) || (shiftXRS < (-SHIFT_THRESHOLD))) ||
             ((shiftYRS > SHIFT_THRESHOLD) || (shiftYRS < (-SHIFT_THRESHOLD))))
         {
             bestNonwrappingShift(Iref, IrefTransforms.FFTI, aux.IauxRS, shiftXRS, shiftYRS, aux2);
             MAT_ELEM(aux.ARS,0,2) += shiftXRS;
             MAT_ELEM(aux.ARS,1,2) += shiftYRS;
             applyGeometry(xmipp_transformation::LINEAR, aux.IauxRS, I, aux.ARS, xmipp_transformation::IS_NOT_INV, wrap);
         }
     }
 
     double corrRS = correlationIndex(aux.IauxRS, Iref);
     double corrSR = correlationIndex(aux.IauxSR, Iref);
     double corr;
     if (corrRS > corrSR)
     {
         I = aux.IauxRS;
         M = aux.ARS;
         corr = corrRS;
     }
     else
     {
         I = aux.IauxSR;
         M = aux.ASR;
         corr = corrSR;
     }
     return corr;
 }

◆ alignImagesConsideringMirrors() [1/2]

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

Fast alignment of two images considering mirrors. The transforms of Iref are presumed to be precomputed in IrefTransforms.

Definition at line 2150 of file filters.cpp.

 {
     MultidimArray<double> Imirror;
     Matrix2D<double> Mmirror;
     Imirror = I;
     Imirror.selfReverseX();
     Imirror.setXmippOrigin();
 
     double corr=alignImages(Iref, IrefTransforms, I, M, wrap, aux, aux2, aux3);
     double corrMirror=alignImages(Iref, IrefTransforms, Imirror, Mmirror, wrap, aux, aux2, aux3);
     if (mask!=nullptr)
     {
         corr = correlationIndex(Iref, I, mask);
         corrMirror = correlationIndex(Iref, Imirror, mask);
     }
     double bestCorr = corr;
     if (corrMirror > bestCorr)
     {
         bestCorr = corrMirror;
         I = Imirror;
         M = Mmirror;
         MAT_ELEM(M,0,0) *= -1;
         MAT_ELEM(M,1,0) *= -1;
     }
     return bestCorr;
 }

◆ alignImagesConsideringMirrors() [2/2]

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

Align two images considering mirrors

Definition at line 2180 of file filters.cpp.

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

◆ alignSetOfImages()

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

Align a set of images. Align a set of images and produce a class average as well as the set of alignment parameters. The output is in Iavg. The metadata is modified by adding the alignment information (transformation and euclidean distance to the average. The process is iterative, first an average is computed. All images are aligned to the average, and this is updated. The process is run for a given number of iterations.

Definition at line 2199 of file filters.cpp.

 {
     Image<double> I;
     MultidimArray<double> InewAvg;
     FileName fnImg;
     AlignmentAux aux;
     CorrelationAux aux2;
     RotationalCorrelationAux aux3;
     Matrix2D<double> M;
     size_t Nimgs;
     size_t Xdim;
     size_t Ydim;
     size_t Zdim;
     getImageSize(MD, Xdim, Ydim, Zdim, Nimgs);
     for (int n = 0; n < Niter; ++n)
     {
         bool lastIteration = (n == (Niter - 1));
         InewAvg.initZeros(Ydim, Xdim);
         InewAvg.setXmippOrigin();
         for (size_t objId : MD.ids())
         {
             MD.getValue(MDL_IMAGE, fnImg, objId);
             I.read(fnImg);
             I().setXmippOrigin();
             double corr;
             if (considerMirror)
                 corr = alignImagesConsideringMirrors(Iavg, I(), M, aux, aux2,
                                                      aux3, xmipp_transformation::WRAP);
             else
                 corr = alignImages(Iavg, I(), M, xmipp_transformation::WRAP, aux, aux2, aux3);
             InewAvg += I();
             if (n == 0)
                 Iavg = InewAvg;
             if (lastIteration)
             {
                 double scale;
                 double shiftx;
                 double shifty;
                 double psi;
                 bool flip;
                 transformationMatrix2Parameters2D(M, flip, scale, shiftx,
                                                   shifty, psi);
                 MD.setValue(MDL_FLIP, flip, objId);
                 MD.setValue(MDL_SHIFT_X, shiftx, objId);
                 MD.setValue(MDL_SHIFT_Y, shifty, objId);
                 MD.setValue(MDL_ANGLE_PSI, psi, objId);
                 MD.setValue(MDL_MAXCC, corr, objId);
             }
         }
         InewAvg /= Nimgs;
         Iavg = InewAvg;
         centerImage(Iavg, aux2, aux3, 4);
     }
 }

◆ bestNonwrappingShift()

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

Translational search (non-wrapping). Assumes that the FFTI1 is already computed.

Definition at line 1903 of file filters.cpp.

 {
     I1.checkDimension(2);
     I2.checkDimension(2);
 
     bestShift(I1, FFTI1, I2, shiftX, shiftY, aux);
     double bestCorr;
     double corr;
     MultidimArray<double> Iaux;
 
     translate(1, Iaux, I1, vectorR2(-shiftX, -shiftY), xmipp_transformation::DONT_WRAP);
     //I1.translate(vectorR2(-shiftX,-shiftY),Iaux, DONT_WRAP);
     bestCorr = corr = fastCorrelation(I2, Iaux);
     double finalX = shiftX;
     double finalY = shiftY;
 #ifdef DEBUG
 
     std::cout << "shiftX=" << shiftX << " shiftY=" << shiftY
     << " corr=" << corr << std::endl;
     ImageXmipp save;
     save()=I1;
     save.write("PPPI1.xmp");
     save()=I2;
     save.write("PPPI2.xmp");
     save()=Iaux;
     save.write("PPPpp.xmp");
 #endif
 
     Iaux.initZeros();
     double testX = (shiftX > 0) ? (shiftX - XSIZE(I1)) : (shiftX + XSIZE(I1));
     double testY = shiftY;
     translate(1, Iaux, I1, vectorR2(-testX, -testY), xmipp_transformation::DONT_WRAP);
     //I1.translate(vectorR2(-testX,-testY),Iaux,DONT_WRAP);
     corr = fastCorrelation(I2, Iaux);
     if (corr > bestCorr)
         finalX = testX;
 #ifdef DEBUG
 
     std::cout << "shiftX=" << testX << " shiftY=" << testY
     << " corr=" << corr << std::endl;
     save()=Iaux;
     save.write("PPPmp.xmp");
 #endif
 
     Iaux.initZeros();
     testX = shiftX;
     testY = (shiftY > 0) ? (shiftY - YSIZE(I1)) : (shiftY + YSIZE(I1));
     translate(1, Iaux, I1, vectorR2(-testX, -testY), xmipp_transformation::DONT_WRAP);
     //I1.translate(vectorR2(-testX,-testY),Iaux,DONT_WRAP);
     corr = fastCorrelation(I2, Iaux);
     if (corr > bestCorr)
         finalY = testY;
 #ifdef DEBUG
 
     std::cout << "shiftX=" << testX << " shiftY=" << testY
     << " corr=" << corr << std::endl;
     save()=Iaux;
     save.write("PPPpm.xmp");
 #endif
 
     Iaux.initZeros();
     testX = (shiftX > 0) ? (shiftX - XSIZE(I1)) : (shiftX + XSIZE(I1));
     testY = (shiftY > 0) ? (shiftY - YSIZE(I1)) : (shiftY + YSIZE(I1));
     translate(1, Iaux, I1, vectorR2(-testX, -testY), xmipp_transformation::DONT_WRAP);
     //I1.translate(vectorR2(-testX,-testY),Iaux,DONT_WRAP);
     corr = fastCorrelation(I2, Iaux);
     if (corr > bestCorr)
     {
         finalX = testX;
         finalY = testY;
     }
 #ifdef DEBUG
     std::cout << "shiftX=" << testX << " shiftY=" << testY
     << " corr=" << corr << std::endl;
     save()=Iaux;
     save.write("PPPmm.xmp");
 #endif
 
     shiftX = finalX;
     shiftY = finalY;
 }

◆ bestShift() [1/4]

template float bestShift	(	MultidimArray< float > &	,
		float &	,
		float &	,
		const MultidimArray< int > *	,
		int
	)

◆ bestShift() [2/4]

template<typename T >

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

Definition at line 1594 of file filters.cpp.

 {
     int imax = INT_MIN;
     int jmax;
     int i_actual;
     int j_actual;
     double xmax;
     double ymax;
     double avecorr;
     double stdcorr;
     double dummy;
     bool neighbourhood = true;
 
     /*
      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 = 0.;
             return -1;
         }
         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((T)0, (T)1);
 
     // Look for maximum shift
     if (maxShift==-1)
         Mcorr.maxIndex(imax, jmax);
     else
     {
         int maxShift2=maxShift*maxShift;
         double bestCorr=std::numeric_limits<T>::lowest();
         for (int i=-maxShift; i<=maxShift; i++)
             for (int j=-maxShift; j<=maxShift; j++)
             {
                 if (i*i+j*j>maxShift2) // continue if the Euclidean distance is too far
                     continue;
                 else if (A2D_ELEM(Mcorr, i, j)>bestCorr)
                 {
                     imax=i;
                     jmax=j;
                     bestCorr=A2D_ELEM(Mcorr, imax, jmax);
                 }
             }
     }
     double max = A2D_ELEM(Mcorr, imax, jmax);
 
     // Estimate n_max around the maximum
     int n_max = -1;
     while (neighbourhood)
     {
         n_max++;
         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 > A2D_ELEM(Mcorr, i_actual, j_actual))
                 {
                     neighbourhood = false;
                     break;
                 }
             }
         }
     }
 
     // We have the neighbourhood => looking for the gravity centre
     xmax = ymax = 0.;
     double sumcorr = 0.;
     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 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 = A2D_ELEM(Mcorr, i_actual, j_actual);
             ymax += i_actual * val;
             xmax += j_actual * val;
             sumcorr += val;
         }
     }
     if (sumcorr != 0)
     {
         shiftX = xmax / sumcorr;
         shiftY = ymax / sumcorr;
     }
     return max;
 }

◆ bestShift() [3/4]

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

Translational search. Assumes that FFTI1 is already computed.

Definition at line 1721 of file filters.cpp.

 {
     I1.checkDimension(2);
     I2.checkDimension(2);
 
     MultidimArray<double> Mcorr;
 
     correlation_matrix(FFTI1, I2, Mcorr, aux);
     return bestShift(Mcorr, shiftX, shiftY, mask, maxShift);
 }

◆ bestShift() [4/4]

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

Translational search. Assumes that FFTI1 and FFTI2 are already computed. Mcorr must already have the right size.

Definition at line 1744 of file filters.cpp.

 {
     correlation_matrix(FFTI1, FFTI2, Mcorr, aux);
     return bestShift(Mcorr, shiftX, shiftY, mask, maxShift);
 }

◆ computeAlignmentTransforms()

void computeAlignmentTransforms	(	const MultidimArray< double > &	I,
		AlignmentTransforms &	ITransforms,
		AlignmentAux &	aux,
		CorrelationAux &	aux2
	)

Definition at line 2035 of file filters.cpp.

 {
     aux2.transformer1.FourierTransform((MultidimArray<double> &)I, ITransforms.FFTI, false);
     polarFourierTransform<true>(I, ITransforms.polarFourierI, false, XSIZE(I) / 5, XSIZE(I) / 2, aux.plans, 1);
 }

◆ computeEdges()

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

Compute edges with Sobel

Definition at line 3517 of file filters.cpp.

 {
     MultidimArray<double> BSpline_coefs;
     BSpline_coefs.initZeros(vol);
     produceSplineCoefficients(3, BSpline_coefs, vol);
     vol_edge.initZeros(vol);
 
     FOR_ALL_ELEMENTS_IN_ARRAY3D(vol)
     {
         double V_dx;
         double V_dy;
         double V_dz;
         V_dx = interpolatedElementBSplineDiffX(BSpline_coefs, j, i, k, 3);
         V_dy = interpolatedElementBSplineDiffY(BSpline_coefs, j, i, k, 3);
         V_dz = interpolatedElementBSplineDiffZ(BSpline_coefs, j, i, k, 3);
         A3D_ELEM(vol_edge,k,i,j) = sqrt(
                                        (V_dx * V_dx) + (V_dy * V_dy) + (V_dz * V_dz));
 
     }
 }

◆ correlationMasked()

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

Masked correlation based on standard deviation values.

Definition at line 1397 of file filters.cpp.

 {
     double mean1;
     double std1;
     double th1;
     I1.computeAvgStdev(mean1,std1);
     th1 = std1;
 
     // Estimate the mean and stddev within the mask
     double N1=0;
     double sumMI1=0;
     double sumMI2=0;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I1)
     {
         double p1=DIRECT_MULTIDIM_ELEM(I1,n);
         double p2=DIRECT_MULTIDIM_ELEM(I2,n);
         if(p1>=th1){
             sumMI1+=p1;
             sumMI2+=p2;
             N1+=1.0;
         }
     }
 
     double sumMI1I2=0.0;
     double sumMI1I1=0.0;
     double sumMI2I2=0.0;
     double iN1=1;
     double avgM1=0;
     double avgM2=0;
     double corrM1M2;
     if (N1>0){
         iN1=1.0/N1;
         avgM1=sumMI1*iN1;
         avgM2=sumMI2*iN1;
     }
     else
         return 0;
 
     double p1a;
     double p2a;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I1)
     {
         double p1=DIRECT_MULTIDIM_ELEM(I1,n);
         double p2=DIRECT_MULTIDIM_ELEM(I2,n);
         if (p1>th1){
             p1a=p1-avgM1;
             p2a=p2-avgM2;
             sumMI1I1 +=p1a*p1a;
             sumMI2I2 +=p2a*p2a;
             sumMI1I2+=p1a*p2a;
         }
     }
     corrM1M2=sumMI1I2/sqrt(sumMI1I1*sumMI2I2);
 
     return corrM1M2;
 }

◆ correlationWeighted()

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

Weighted correlation based on differences between images.

Definition at line 1457 of file filters.cpp.

 {
 
     MultidimArray<double> Idiff;
     MultidimArray<double> I2Aligned;
     Matrix2D<double> M;
 
     I1.setXmippOrigin();
     I2.setXmippOrigin();
     I2Aligned=I2;
     alignImages(I1, I2Aligned, M, false);
     Idiff=I1;
     Idiff-=I2Aligned;
 
     double mean;
     double std;
     Idiff.computeAvgStdev(mean,std);
     Idiff.selfABS();
     double threshold=std;
 
     // Estimate the mean and stddev within the mask
     double N=0;
     double sumWI1=0;
     double sumWI2=0;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Idiff)
     {
         if (DIRECT_MULTIDIM_ELEM(Idiff,n)>threshold)
         {
             double p1=DIRECT_MULTIDIM_ELEM(I1,n);
             double p2Alg=DIRECT_MULTIDIM_ELEM(I2Aligned,n);
             sumWI1+=p1;
             sumWI2+=p2Alg;
             N+=1.0;
         }
     }
 
     double sumWI1I2=0.0;
     double sumWI1I1=0.0;
     double sumWI2I2=0.0;
     double iN;
     double avgW1;
     double avgW2;
     double corrW1W2;
     if(N>0){
         iN=1.0/N;
         avgW1=sumWI1*iN;
         avgW2=sumWI2*iN;
     }
     else
         return -1;
     double p1a;
     double p2a;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I1)
     {
         if (DIRECT_MULTIDIM_ELEM(Idiff,n)>threshold)
         {
             double p1=DIRECT_MULTIDIM_ELEM(I1,n);
             double p2Alg=DIRECT_MULTIDIM_ELEM(I2Aligned,n);
             p1a=p1-avgW1;
             p2a=p2Alg-avgW2;
 
             double w=DIRECT_MULTIDIM_ELEM(Idiff,n);
             double wp1a=w*p1a;
             double wp2a=w*p2a;
 
             sumWI1I1 +=wp1a*p1a;
             sumWI2I2 +=wp2a*p2a;
             sumWI1I2 +=wp1a*p2a;
         }
     }
     corrW1W2=sumWI1I2/sqrt(sumWI1I1*sumWI2I2);
 
     return corrW1W2;
 }

◆ covarianceMatrix()

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

Covariance matrix of an image

Definition at line 1582 of file filters.cpp.

 {
     Matrix2D<double> mI;
     I.copy(mI);
     subtractColumnMeans(mI);
     matrixOperation_AtA(mI,C);
     C*=1.0/(YSIZE(I)-1.0);
 }

◆ fastBestRotation() [1/3]

double fastBestRotation	(	const MultidimArray< double > &	IrefCyl,
		const MultidimArray< double > &	Icyl,
		CorrelationAux &	aux,
		VolumeAlignmentAux &	aux2,
		double	deltaAng
	)

Definition at line 2255 of file filters.cpp.

 {
     correlation_matrix(IrefCyl, Icyl, aux2.corr, aux, false);
     STARTINGZ(aux2.corr) = STARTINGY(aux2.corr) = STARTINGX(aux2.corr) = 0;
     double bestCorr = A3D_ELEM(aux2.corr,0,STARTINGY(aux2.corr),0);
     double bestAngle = 0;
     for (int i = STARTINGY(aux2.corr) + 1; i <= FINISHINGY(aux2.corr); i++)
     {
         double corr = A3D_ELEM(aux2.corr,0,i,0);
         if (corr > bestCorr)
         {
             bestCorr = corr;
             bestAngle = i;
         }
     }
     bestAngle *= deltaAng * 180.0 / PI;
     return -bestAngle;
 }

◆ fastBestRotation() [2/3]

void fastBestRotation	(	const MultidimArray< double > &	IrefCylZ,
		const MultidimArray< double > &	IrefCylY,
		const MultidimArray< double > &	IrefCylX,
		const MultidimArray< double > &	I,
		const MultidimArray< double > &	Icurrent,
		MultidimArray< double > &	Ifinal,
		char	axis,
		Matrix2D< double > &	R,
		CorrelationAux &	aux,
		VolumeAlignmentAux &	aux2
	)

Definition at line 2318 of file filters.cpp.

 {
     double bestAngle;
     if (axis=='Z')
         bestAngle = fastBestRotationAroundZ(IrefCylZ, Icurrent, aux, aux2);
     else if (axis=='Y')
         bestAngle = fastBestRotationAroundY(IrefCylY, Icurrent, aux, aux2);
     else
         bestAngle = fastBestRotationAroundX(IrefCylX, Icurrent, aux, aux2);
     Matrix2D<double> Raux;
     rotation3DMatrix(bestAngle, axis, Raux);
     R=Raux*R;
     applyGeometry(xmipp_transformation::LINEAR, Ifinal, I, R, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);
 }

◆ fastBestRotation() [3/3]

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

Find a ZYZ rotation that transforms I into Iref. The rotation matrix is returned in R. I is modified to be aligned.

Definition at line 2341 of file filters.cpp.

 {
     R.initIdentity(4);
     fastBestRotation(IrefCylZ,IrefCylY,IrefCylX,I,I,aux2.I1,eulerAngles[0],R, aux,aux2);
     fastBestRotation(IrefCylZ,IrefCylY,IrefCylX,I,aux2.I1,aux2.I12,eulerAngles[1],R,aux,aux2);
     fastBestRotation(IrefCylZ,IrefCylY,IrefCylX,I,aux2.I12,aux2.I123,eulerAngles[2],R,aux,aux2);
     I=aux2.I123;
 }

◆ fastBestRotationAroundX()

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

Fast best rotation around X

Definition at line 2304 of file filters.cpp.

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

◆ fastBestRotationAroundY()

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

Fast best rotation around Y

Definition at line 2290 of file filters.cpp.

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

◆ fastBestRotationAroundZ()

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

Fast best rotation around Z

Definition at line 2276 of file filters.cpp.

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

◆ fastCorrentropy()

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

Correntropy with mask.

Definition at line 1244 of file filters.cpp.

 {
     double retvalxy = 0;
     double isigma = 1.0 / sigma;
     int maskSum = 0;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(x)
     {
         if (DIRECT_MULTIDIM_ELEM(mask,n))
         {
             retvalxy += G.getValue(
                             isigma
                             * (DIRECT_MULTIDIM_ELEM(x,n)
                                - DIRECT_MULTIDIM_ELEM(y,n)));
             ++maskSum;
         }
     }
     if (0 == maskSum) {
         REPORT_ERROR(ERR_NUMERICAL, "maskSum is zero (0), which would lead to division by zero");
     }
     return (retvalxy / maskSum);
 }

◆ forceDWTSparsity()

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

Force sparsity. Only eps*100 % of the DWT coefficients are kept. It is assumed that the input volume is cubic.

Definition at line 3539 of file filters.cpp.

 {
     int size0=XSIZE(V);
     int sizeF=(int)NEXT_POWER_OF_2(size0);
     selfScaleToSize(xmipp_transformation::BSPLINE3,V,sizeF,sizeF,sizeF);
     MultidimArray<double> vol_wavelets;
     MultidimArray<double> vol_wavelets_abs;
     set_DWT_type(DAUB12);
     DWT(V,vol_wavelets);
     vol_wavelets_abs=vol_wavelets;
     vol_wavelets_abs.selfABS();
     double *begin=MULTIDIM_ARRAY(vol_wavelets_abs);
     double *end=MULTIDIM_ARRAY(vol_wavelets_abs)+MULTIDIM_SIZE(vol_wavelets_abs);
     std::sort(begin,end);
     double threshold1=DIRECT_MULTIDIM_ELEM(vol_wavelets_abs,
                                            (long int)((1.0-eps)*MULTIDIM_SIZE(vol_wavelets_abs)));
     vol_wavelets.threshold("abs_below", threshold1, 0.0);
     IDWT(vol_wavelets,V);
     selfScaleToSize(xmipp_transformation::BSPLINE3,V,size0, size0, size0);
 }

◆ imedDistance()

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

IMED distance. See Wang, L et al. On the Euclidean distance of images. IEEE PAMI 27: 1334-1339 (2005)

Definition at line 1269 of file filters.cpp.

 {
     // [x,y]=meshgrid([-3:1:3],[-3:1:3])
     // format long  // w=1/sqrt(2*pi)*exp(-0.5*(x.*x+y.*y))
     double *refW;
     double w[49]={
             0.000049233388666,   0.000599785460091,   0.002688051941039,   0.004431848411938,   0.002688051941039,   0.000599785460091,   0.000049233388666,
             0.000599785460091,   0.007306882745281,   0.032747176537767,   0.053990966513188,   0.032747176537767,   0.007306882745281,   0.000599785460091,
             0.002688051941039,   0.032747176537767,   0.146762663173740,   0.241970724519143,   0.146762663173740,   0.032747176537767,   0.002688051941039,
             0.004431848411938,   0.053990966513188,   0.241970724519143,   0.398942280401433,   0.241970724519143,   0.053990966513188,   0.004431848411938,
             0.002688051941039,   0.032747176537767,   0.146762663173740,   0.241970724519143,   0.146762663173740,   0.032747176537767,   0.002688051941039,
             0.000599785460091,   0.007306882745281,   0.032747176537767,   0.053990966513188,   0.032747176537767,   0.007306882745281,   0.000599785460091,
             0.000049233388666,   0.000599785460091,   0.002688051941039,   0.004431848411938,   0.002688051941039,   0.000599785460091,   0.000049233388666};
 
     MultidimArray<double> diffImage = I1 - I2;
     double  imed = 0;
     int imiddle=YSIZE(I1)/2;
     int jmiddle=XSIZE(I1)/2;
     int R2max=imiddle*imiddle;
     int ysize=(int)YSIZE(I1);
     int xsize=(int)XSIZE(I1);
     for (int i=3; i<ysize-3; ++i)
     {
         int i2=(i-imiddle)*(i-imiddle);
         for (int j=3; j<xsize-3; ++j)
         {
             int j2=(j-jmiddle)*(j-jmiddle);
             if (i2+j2>R2max) // Measure only within the maximum circle
                 continue;
             double diffi=DIRECT_A2D_ELEM(diffImage,i,j);
             int index=0;
             for (int ii=-3; ii<=3; ++ii)
             {
                 refW = &w[index];
                 index = index + 7;
                 double *diffj=&DIRECT_A2D_ELEM(diffImage,i+ii,j-3);
                 double  imedAux=(*refW++)*(*diffj++);
                 imedAux+=(*refW++)*(*diffj++);
                 imedAux+=(*refW++)*(*diffj++);
                 imedAux+=(*refW++)*(*diffj++);
                 imedAux+=(*refW++)*(*diffj++);
                 imedAux+=(*refW++)*(*diffj++);
                 imedAux+=(*refW++)*(*diffj++);
 
                 imed+=imedAux*diffi;
             }
         }
     }
     return sqrt(imed);
 }

◆ imedNormalizedDistance()

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

IMED normalized distance.

Definition at line 1321 of file filters.cpp.

 {
     // [x,y]=meshgrid([-3:1:3],[-3:1:3])
     // format long
     // w=1/sqrt(2*pi)*exp(-0.5*(x.*x+y.*y))
     const double w[7][7]={
             {0.000049233388666,   0.000599785460091,   0.002688051941039,   0.004431848411938,   0.002688051941039,   0.000599785460091,   0.000049233388666},
             {0.000599785460091,   0.007306882745281,   0.032747176537767,   0.053990966513188,   0.032747176537767,   0.007306882745281,   0.000599785460091},
             {0.002688051941039,   0.032747176537767,   0.146762663173740,   0.241970724519143,   0.146762663173740,   0.032747176537767,   0.002688051941039},
             {0.004431848411938,   0.053990966513188,   0.241970724519143,   0.398942280401433,   0.241970724519143,   0.053990966513188,   0.004431848411938},
             {0.002688051941039,   0.032747176537767,   0.146762663173740,   0.241970724519143,   0.146762663173740,   0.032747176537767,   0.002688051941039},
             {0.000599785460091,   0.007306882745281,   0.032747176537767,   0.053990966513188,   0.032747176537767,   0.007306882745281,   0.000599785460091},
             {0.000049233388666,   0.000599785460091,   0.002688051941039,   0.004431848411938,   0.002688051941039,   0.000599785460091,   0.000049233388666}
     };
     int imiddle=YSIZE(I1)/2;
     int jmiddle=XSIZE(I1)/2;
     int R2max=imiddle*imiddle;
     // Measure the average of both images
     double avg1=0;
     double avg2=0;
     double N=0;
     int ydim=YSIZE(I1);
     int xdim=XSIZE(I1);
     for (int i=3; i<ydim-3; ++i)
     {
         int i2=(i-imiddle)*(i-imiddle);
         for (int j=3; j<xdim-3; ++j)
         {
             int j2=(j-jmiddle)*(j-jmiddle);
             if (i2+j2>R2max) // Measure only within the maximum circle
                 continue;
             avg1+=DIRECT_A2D_ELEM(I1,i,j);
             avg2+=DIRECT_A2D_ELEM(I2,i,j);
             ++N;
         }
     }
     avg1/=N;
     avg2/=N;
     double diffAvg=avg1-avg2;
 
     double imed = 0;
     double imed1=0;
     double imed2=0;
     for (int i=3; i<ydim-3; ++i)
     {
         int i2=(i-imiddle)*(i-imiddle);
         for (int j=3; j<xdim-3; ++j)
         {
             int j2=(j-jmiddle)*(j-jmiddle);
             if (i2+j2>R2max) // Measure only within the maximum circle
                 continue;
             double subtractedI1ij=DIRECT_A2D_ELEM(I1,i,j)-avg1;
             double subtractedI2ij=DIRECT_A2D_ELEM(I2,i,j)-avg2;
             double diffi=DIRECT_A2D_ELEM(I1,i,j)-DIRECT_A2D_ELEM(I2,i,j)-diffAvg;
             for (int ii=-3; ii<=3; ++ii)
             {
                 int i_ii=i+ii;
                 for (int jj=-3; jj<=3; ++jj)
                 {
                     int j_jj=j+jj;
                     //double diffj=DIRECT_A2D_ELEM(I1,i_ii,j_jj)-DIRECT_A2D_ELEM(I2,i_ii,j_jj)-diffAvg;
                     double diffj=DIRECT_A2D_ELEM(I1,i,j)-DIRECT_A2D_ELEM(I2,i_ii,j_jj)-diffAvg;
                     double wiijj=w[ii+3][jj+3];
                     imed+=wiijj*diffi*diffj;
                     imed1+=wiijj*subtractedI1ij*(DIRECT_A2D_ELEM(I1,i_ii,j_jj)-avg1);
                     imed2+=wiijj*subtractedI2ij*(DIRECT_A2D_ELEM(I2,i_ii,j_jj)-avg2);
                 }
             }
         }
     }
     return imed/sqrt((imed1+imed2)/2);
 }

◆ matlab_filter()

void matlab_filter	(	Matrix1D< double > &	B,
		Matrix1D< double > &	A,
		const MultidimArray< double > &	X,
		MultidimArray< double > &	Y,
		MultidimArray< double > &	Z
	)

Definition at line 4256 of file filters.cpp.

 {
         if (VEC_ELEM(A,0)!=1.0)
         {
             A/=VEC_ELEM(A,0);
             B/=VEC_ELEM(B,0);
         }
 
         size_t input_size = XSIZE(X);
         size_t filter_order = std::max(VEC_XSIZE(A),VEC_XSIZE(B));
         if (VEC_XSIZE(B)!=filter_order)
             B.resize(filter_order);
         if (VEC_XSIZE(A)!=filter_order)
             A.resize(filter_order);
         Y.resize(input_size);
         Z.initZeros(filter_order);
 
         for (size_t i = 0; i < input_size; ++i)
         {
             size_t order = filter_order - 1;
             while (order)
             {
                 if (i >= order)
                 {
                     DIRECT_A1D_ELEM(Z,order - 1) = VEC_ELEM(B,order) * DIRECT_A1D_ELEM(X,i - order) -
                                                    VEC_ELEM(A,order) * DIRECT_A1D_ELEM(Y,i - order) + DIRECT_A1D_ELEM(Z,order);
                 }
                 --order;
             }
             DIRECT_A1D_ELEM(Y,i) = VEC_ELEM(B,0) * DIRECT_A1D_ELEM(X,i) + DIRECT_A1D_ELEM(Z,0);
         }
 }

◆ regionGrowing3D()

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

Region growing for volumes

Given a position inside a volume this function grows a region with (filling_colour) until it finds a border of value (stop_colour). If the point is outside the volume 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.

Definition at line 412 of file filters.cpp.

 {
     V_in.checkDimension(3);
 
     std::list<Coordinate3D> iNeighbours; /* A list for neighbour voxels */
     int iCurrentk;
     int iCurrenti;
     int iCurrentj; /* Coordinates of the current voxel considered */
 
     /* First task is copying the input volume into the output one */
     V_out = V_in;
 
     /**** Then the region growing is done in output volume ****/
     /* Insert at the beginning of the list the seed coordinates */
     Coordinate3D coord;
     coord.ii = i;
     coord.jj = j;
     coord.kk = k;
     iNeighbours.push_front(coord);
 
     /* Fill the seed coordinates */
     A3D_ELEM(V_out, k, 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;
         iCurrentk = coord.kk;
 
         /* a macro for doing exploration of a voxel. If the voxel has a value
          lower than stop_colour, its filled with filling colour and added to the
          list for exploring its neighbours */
 #define CHECK_POINT_3D(k,i,j) \
  {\
      int I=i; \
         int J=j; \
         int K=k; \
   if (INSIDEXYZ(V_out,J,I,K))  { \
    double voxel=A3D_ELEM(V_out,K,I,J); \
    if (voxel!=filling_colour) \
     if ((less && voxel < stop_colour)|| \
      (!less &&voxel > stop_colour)) { \
      coord.ii=I; \
      coord.jj=J; \
      coord.kk=K; \
      A3D_ELEM (V_out,coord.kk,coord.ii,coord.jj)=filling_colour; \
      iNeighbours.push_front(coord); \
     } \
   }\
  }
 
         /* Make the exploration of the pixel neighbours */
         CHECK_POINT_3D(iCurrentk, iCurrenti, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk, iCurrenti, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj);
         CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj);
         CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj + 1);
     }
 }

◆ regionGrowing3DEqualValue()

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

Region growing with equal values for volumes

Given a position inside a volume this function grows a region with (filling_value) until it finds a border. If the point is outside the volume then nothing is done.

The growing only considers those voxels with value equal to the selected voxel

Definition at line 499 of file filters.cpp.

 {
     V_in.checkDimension(3);
 
     std::list<Coordinate3D> iNeighbours; /* A list for neighbour voxels */
     int iCurrentk;
     int iCurrenti;
     int iCurrentj; /* Coordinates of the current voxel considered */
 
     /* First task is copying the input volume into the output one */
     V_out.resizeNoCopy(V_in);
     V_out.initConstant(1);
 
     /**** Then the region growing is done in output volume ****/
     /* Insert at the beginning of the list the seed coordinates */
     Coordinate3D coord;
     coord.ii = STARTINGY(V_in);
     coord.jj = STARTINGX(V_in);
     coord.kk = STARTINGZ(V_in);
     iNeighbours.push_front(coord);
 
     /* Fill the seed coordinates */
     double bgColor = A3D_ELEM(V_in, STARTINGZ(V_in), STARTINGY(V_in), STARTINGX(V_in));
     A3D_ELEM(V_out, STARTINGZ(V_out), STARTINGY(V_out), STARTINGX(V_out)) = 0;
 
     while (!iNeighbours.empty())
     {
         /* Take the current pixel to explore */
         coord = iNeighbours.front();
         iNeighbours.pop_front();
         iCurrenti = coord.ii;
         iCurrentj = coord.jj;
         iCurrentk = coord.kk;
 
         /* a macro for doing exploration of a voxel. If the voxel has a value
          lower than stop_colour, its filled with filling colour and added to the
          list for exploring its neighbours */
 #undef CHECK_POINT_3D
 #define CHECK_POINT_3D(k,i,j) \
  {\
         int I=i; \
         int J=j; \
         int K=k; \
   if (INSIDEXYZ(V_out,J,I,K))  { \
    if (A3D_ELEM(V_in,K,I,J)==bgColor && A3D_ELEM(V_out,K,I,J)==1) { \
      coord.ii=I; \
      coord.jj=J; \
      coord.kk=K; \
      A3D_ELEM (V_out,coord.kk,coord.ii,coord.jj)=filling_value; \
      iNeighbours.push_front(coord); \
     } \
   }\
  }\
 
         /* Make the exploration of the pixel neighbours */
         CHECK_POINT_3D(iCurrentk, iCurrenti, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk, iCurrenti, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj);
         CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk, iCurrenti - 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj);
         CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk, iCurrenti + 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti - 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk - 1, iCurrenti + 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti - 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj + 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj - 1);
         CHECK_POINT_3D(iCurrentk + 1, iCurrenti + 1, iCurrentj + 1);
     }
 }

◆ Shah_energy()

double Shah_energy	(	const MultidimArray< double > &	img,
		const MultidimArray< double > &	surface_strength,
		const MultidimArray< double > &	edge_strength,
		double	K,
		const Matrix1D< double > &	W
	)

Definition at line 2511 of file filters.cpp.

 {
     img.checkDimension(2);
 
     int Ydim1 = YSIZE(img) - 1;
     int Xdim1 = XSIZE(img) - 1;
 
     double Kinv = 1.0 / K;
 
     /* Calculate surface energy */
     double E1 = 0.0;
     double E2 = 0.0;
     double E3 = 0.0;
     double E4 = 0.0;
     double w0 = VEC_ELEM(W,0);
     double w1 = VEC_ELEM(W,1);
     double w2 = VEC_ELEM(W,2);
     double w3 = VEC_ELEM(W,3);
     for (int i = 1; i < Ydim1; i++)
     {
         int ip1 = i + 1;
         int im1 = i - 1;
         for (int j = 1; j < Xdim1; j++)
         {
             int jp1 = j + 1;
             int jm1 = j - 1;
             /* Calculate data matching terms */
             double D = dAij(img, i, j);
             double F = dAij(surface_strength, i, j);
             double S = dAij(edge_strength, i, j);
             double diff = D - F;
             E1 += w0 * diff * diff;
             E3 += w2 * K * S * S;
 
             /* Calculate first derivative terms */
             double Fx = 0.5
                         * (dAij(surface_strength, i, jp1)
                            - dAij(surface_strength, i, jm1));
             double Fy = 0.5
                         * (dAij(surface_strength, ip1, j)
                            - dAij(surface_strength, im1, j));
             double Sx =
                 0.5
                 * (dAij(edge_strength, i, jp1)
                    - dAij(edge_strength, i, jm1));
             double Sy =
                 0.5
                 * (dAij(edge_strength, ip1, j)
                    - dAij(edge_strength, im1, j));
             double S_1 = 1 - S;
             E2 += w1 * S_1 * S_1 * (Fx * Fx + Fy * Fy);
             E4 += w3 * Kinv * (Sx * Sx + Sy * Sy);
         }
     }
 
     return E1 + E2 + E3 + E4; // Total energy
 }

◆ svdCorrelation()

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

SVD correlation.

Definition at line 1535 of file filters.cpp.

 {
     // Copy input images into matrix2Ds
     Matrix2D<double> mI1;
     Matrix2D<double> mI2;
     I1.copy(mI1);
     I2.copy(mI2);
 
     // SVD Decomposition of I1
     Matrix2D<double> U1;
     Matrix2D<double> V1;
     Matrix1D<double> W1;
     svdcmp(mI1,U1,W1,V1);
 
     // Express I2 in the basis of I1
     Matrix2D<double> aux;
     Matrix2D<double> W2;
     matrixOperation_AB(mI2,V1,aux);
     matrixOperation_AtB(U1,aux,W2);
 
     // Make W2 to be diagonal
     FOR_ALL_ELEMENTS_IN_MATRIX2D(W2)
     if (i!=j)
         MAT_ELEM(W2,i,j)=0.0;
 
     // Recover I2p
     matrixOperation_ABt(W2,V1,aux);
     matrixOperation_AB(U1,aux,mI2);
 
     // Measure correlation
     MultidimArray<double> I2p;
     I2p=mI2;
 //  Image<double> save;
 //  save()=I1;
 //  save.write("PPPI1.xmp");
 //  save()=I2;
 //  save.write("PPPI2.xmp");
 //  save()=I2p;
 //  save.write("PPPI2p.xmp");
 //  std::cout << "Correlation= " << correlationIndex(I1,I2p,mask) << std::endl;
 //  std::cout << "Press any key\n";
 //  char c; std::cin >> c;
     // double ratio=I2p.computeStddev()/I2.computeStddev();
     return correlationIndex(I1,I2p,mask); //*ratio;
 }

◆ unnormalizedGaussian2D()

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

Unnormalized 2D gaussian value using covariance

This function returns the value of a multivariate (2D) gaussian function at the point r (column vector of dimension 2).

G(r,mu,sigma)=exp(-0.5 * (r-mu)^t sigma^-1 (r-mu))

Definition at line 2379 of file filters.cpp.

 {
     double x = XX(r) - XX(mu);
     double y = YY(r) - YY(mu);
     return exp(
                -0.5
                * (sigmainv(0, 0) * x * x + 2 * sigmainv(0, 1) * x * y
                   + sigmainv(1, 1) * y * y));
 }

◆ Update_edge_Shah()

double Update_edge_Shah	(	MultidimArray< double > &	img,
		MultidimArray< double > &	surface_strength,
		MultidimArray< double > &	edge_strength,
		double	K,
		const Matrix1D< double > &	W
	)

Definition at line 2657 of file filters.cpp.

 {
     img.checkDimension(2);
 
     double Diff = 0.0;
     double Norm = 0.0;
     size_t Ydim1 = YSIZE(img) - 1;
     size_t Xdim1 = XSIZE(img) - 1;
     double Kinv = 1.0 / K;
 
     /* Update edge estimate */
     double w1 = VEC_ELEM(W,1);
     double w2 = VEC_ELEM(W,2);
     double w3 = VEC_ELEM(W,3);
     for (size_t i = 1; i < Ydim1; i++)
     {
         int ip1 = i + 1;
         int im1 = i - 1;
         for (size_t j = 1; j < Xdim1; j++)
         {
             int jp1 = j + 1;
             int jm1 = j - 1;
 
             /* Calculate first and second derivative terms */
             double Fx = 0.5
                         * (dAij(surface_strength, i, jp1)
                            - dAij(surface_strength, i, jm1));
             double Fy = 0.5
                         * (dAij(surface_strength, ip1, j)
                            - dAij(surface_strength, im1, j));
             double Constant = w1 * (Fx * Fx + Fy * Fy);
 
             /* Calculate weights for central pixel and neighbors */
             double Central = w2 * K + w3 * Kinv * 4;
             double Neighbors = w3 * Kinv
                                * (dAij(edge_strength, im1, j) + dAij(edge_strength, ip1, j)
                                   + dAij(edge_strength, i, jm1)
                                   + dAij(edge_strength, i, jp1));
 
             /* Calculate new S value */
             double Old_edge_strength = dAij(edge_strength, i, j);
             double S = (Constant + Neighbors) / (Constant + Central);
             if (S < 0)
                 dAij(edge_strength, i, j) *= 0.5;
             else if (S > 1)
                 dAij(edge_strength, i, j) = 0.5
                                             * (dAij(edge_strength, i, j) + 1);
             else
                 dAij(edge_strength, i, j) = S;
 
             // Compute the difference.
             Diff += fabs(dAij(edge_strength, i, j) - Old_edge_strength);
             Norm += fabs(Old_edge_strength);
         }
     }
     return Diff / Norm; // Return the relative difference.
 }

◆ Update_surface_Shah()

double Update_surface_Shah	(	MultidimArray< double > &	img,
		MultidimArray< double > &	surface_strength,
		MultidimArray< double > &	edge_strength,
		const Matrix1D< double > &	W
	)

Definition at line 2577 of file filters.cpp.

 {
     img.checkDimension(2);
 
     double Diff = 0.0;
     double Norm = 0.0;
     size_t Ydim1 = YSIZE(img) - 1;
     size_t Xdim1 = XSIZE(img) - 1;
 
     /* Update surface estimate */
     double w0 = VEC_ELEM(W,0);
     double w1 = VEC_ELEM(W,1);
     for (size_t i = 1; i < Ydim1; i++)
     {
         int ip1 = i + 1;
         int im1 = i - 1;
         for (size_t j = 1; j < Xdim1; j++)
         {
             int jp1 = j + 1;
             int jm1 = j - 1;
             /* Calculate edge partial derivative terms */
             double S = dAij(edge_strength, i, j);
             double Sx =
                 0.5
                 * (dAij(edge_strength, i, jp1)
                    - dAij(edge_strength, i, jm1));
             double Sy =
                 0.5
                 * (dAij(edge_strength, ip1, j)
                    - dAij(edge_strength, im1, j));
 
             double nS = 1 - S;
             double nS2 = nS * nS;
 
             /* Calculate surface partial derivative terms (excluding central pixel) */
             double F;
             double D;
             F = D = dAij(img, i, j);
             double SS_i_jp1 = dAij(surface_strength, i, jp1);
             double SS_i_jm1 = dAij(surface_strength, i, jm1);
             double SS_ip1_j = dAij(surface_strength, ip1, j);
             double SS_im1_j = dAij(surface_strength, im1, j);
             double Fx = 0.5 * (SS_i_jp1 - SS_i_jm1);
             double Fy = 0.5 * (SS_ip1_j - SS_im1_j);
             double Fxx = SS_i_jp1 + SS_i_jm1;
             double Fyy = SS_ip1_j + SS_im1_j;
 
             /* Calculate surface partial derivative weights */
             double wFx = 4 * w1 * nS * Sx;
             double wFy = 4 * w1 * nS * Sy;
             double wFxx = -2 * w1 * nS2;
             double wFyy = -2 * w1 * nS2;
 
             /* Calculate new surface value */
             double Constant = -2 * w0 * D;
             double Central = -2 * w0 + 2 * wFxx + 2 * wFyy;
             double Neighbors = wFx * Fx + wFy * Fy + wFxx * Fxx + wFyy * Fyy;
 
             if (fabs(Central) > XMIPP_EQUAL_ACCURACY)
                 F = (Constant + Neighbors) / Central;
             F = CLIP(F, 0, 1);
 
             // Compute the difference.
             double SS_i_j = dAij(surface_strength, i, j);
             Diff += fabs(SS_i_j - F);
             Norm += fabs(SS_i_j);
 
             // Update the new value.
             dAij(surface_strength, i, j) = F;
         }
     }
     return Diff / Norm; // Return the relative difference.
 }

Variable Documentation

◆ INITIAL_ROTATE_THRESHOLD

constexpr float INITIAL_ROTATE_THRESHOLD = ROTATE_THRESHOLD + 1.0

Definition at line 2045 of file filters.cpp.

◆ INITIAL_SHIFT_THRESHOLD

constexpr double INITIAL_SHIFT_THRESHOLD = SHIFT_THRESHOLD + 1.0

Definition at line 2044 of file filters.cpp.

◆ ROTATE_THRESHOLD

constexpr float ROTATE_THRESHOLD = 1.0

Definition at line 2043 of file filters.cpp.

◆ SHAH_CONVERGENCE_THRESHOLD

constexpr double SHAH_CONVERGENCE_THRESHOLD = 0.0001

Definition at line 2719 of file filters.cpp.

◆ SHIFT_THRESHOLD

constexpr double SHIFT_THRESHOLD = 0.95

Definition at line 2042 of file filters.cpp.

Classes

Macros

Functions

Variables

Macro Definition Documentation

◆ CHECK_POINT

◆ CHECK_POINT_3D [1/2]

◆ CHECK_POINT_3D [2/2]

◆ CHECK_POINT_DIST

Function Documentation

◆ alignImages()

◆ alignImagesConsideringMirrors() [1/2]

◆ alignImagesConsideringMirrors() [2/2]

◆ alignSetOfImages()

◆ bestNonwrappingShift()

◆ bestShift() [1/4]

◆ bestShift() [2/4]

◆ bestShift() [3/4]

◆ bestShift() [4/4]

◆ computeAlignmentTransforms()

◆ computeEdges()

◆ correlationMasked()

◆ correlationWeighted()

◆ covarianceMatrix()

◆ fastBestRotation() [1/3]

◆ fastBestRotation() [2/3]

◆ fastBestRotation() [3/3]

◆ fastBestRotationAroundX()

◆ fastBestRotationAroundY()

◆ fastBestRotationAroundZ()

◆ fastCorrentropy()

◆ forceDWTSparsity()

◆ imedDistance()

◆ imedNormalizedDistance()

◆ matlab_filter()

◆ regionGrowing3D()

◆ regionGrowing3DEqualValue()

◆ Shah_energy()

◆ svdCorrelation()

◆ unnormalizedGaussian2D()

◆ Update_edge_Shah()

◆ Update_surface_Shah()

Variable Documentation

◆ INITIAL_ROTATE_THRESHOLD

◆ INITIAL_SHIFT_THRESHOLD

◆ ROTATE_THRESHOLD

◆ SHAH_CONVERGENCE_THRESHOLD

◆ SHIFT_THRESHOLD