wavelet.cpp File Reference

#include <core/bilib/wavelet.h>
#include <core/args.h>
#include <core/histogram.h>
#include <core/xmipp_image.h>
#include <core/xmipp_fftw.h>
#include "wavelet.h"
#include "numerical_tools.h"
#include "mask.h"

Include dependency graph for wavelet.cpp:

Go to the source code of this file.

Macros
#define	DWT_Scale(i, smax)   ((int)((i==0)?smax-1:(ABS((CEIL(log10((double)(i+1))/log10(2.0))-smax)))))
#define	DWT_Quadrant1D(i, s, smax)   ((s!=smax-1)?'1':((i==0)?'0':'1'))
#define	DWT_QuadrantnD(i, s, sp, smax)   ((s!=sp)?'0':DWT_Quadrant1D(i,s,smax))
#define	DWT_icoef1D(i, s, smax)   ()

Functions
void	Bilib_DWT (const MultidimArray< double > &input, MultidimArray< double > &result, int iterations, int isign)
void	set_DWT_type (int DWT_type)
void	IDWT (const MultidimArray< double > &v, MultidimArray< double > &result)
void	DWT_lowpass2D (const MultidimArray< double > &v, MultidimArray< double > &result)
std::string	Quadrant2D (int q)
std::string	Quadrant3D (int q)
void	Get_Scale_Quadrant (int size_x, int x, int &scale, std::string &quadrant)
void	Get_Scale_Quadrant (int size_x, int size_y, int x, int y, int &scale, std::string &quadrant)
void	Get_Scale_Quadrant (int size_x, int size_y, int size_z, int x, int y, int z, int &scale, std::string &quadrant)
void	clean_quadrant2D (MultidimArray< double > &I, int scale, const std::string &quadrant)
void	clean_quadrant3D (MultidimArray< double > &I, int scale, const std::string &quadrant)
void	soft_thresholding (MultidimArray< double > &I, double th)
void	adaptive_soft_thresholding_block2D (MultidimArray< double > &I, int scale, const std::string &quadrant, double sigma)
double	compute_noise_power (MultidimArray< double > &I)
void	adaptive_soft_thresholding2D (MultidimArray< double > &I, int scale)
void	DWT_keep_central_part (MultidimArray< double > &I, double R)
void	DWT_Bijaoui_denoise_LL2D (MultidimArray< double > &WI, int scale, const std::string &orientation, double mu, double S, double N)
void	DWT_Bijaoui_denoise_LL3D (MultidimArray< double > &WI, int scale, const std::string &orientation, double mu, double S, double N)
void	bayesian_solve_eq_system (const Matrix1D< double > &power, const Matrix1D< double > &Ncoefs, double SNR0, double SNRF, double powerI, double power_rest, bool white_noise, Matrix1D< double > &estimatedS)
Matrix1D< double >	bayesian_wiener_filtering2D (MultidimArray< double > &WI, int allowed_scale, double SNR0, double SNRF, bool white_noise, int tell, bool denoise)
void	bayesian_wiener_filtering2D (MultidimArray< double > &WI, int allowed_scale, Matrix1D< double > &estimatedS)
Matrix1D< double >	bayesian_wiener_filtering3D (MultidimArray< double > &WI, int allowed_scale, double SNR0, double SNRF, bool white_noise, int tell, bool denoise)
void	bayesian_wiener_filtering3D (MultidimArray< double > &WI, int allowed_scale, Matrix1D< double > &estimatedS)
void	phaseCongMono (MultidimArray< double > &I, MultidimArray< double > &Ph, MultidimArray< double > &Or, MultidimArray< double > &Energy, MultidimArray< double > &lowPass, MultidimArray< double > &Radius, MultidimArray< std::complex< double > > &H, int nScale, double minWaveLength, double mult, double sigmaOnf)

Macro Definition Documentation

DWT_icoef1D

#define DWT_icoef1D	(		i,
			s,
			smax
	)		()

Definition at line 246 of file wavelet.cpp.

DWT_Quadrant1D

#define DWT_Quadrant1D	(		i,
			s,
			smax
	)		((s!=smax-1)?'1':((i==0)?'0':'1'))

Definition at line 243 of file wavelet.cpp.

DWT_QuadrantnD

#define DWT_QuadrantnD	(		i,
			s,
			sp,
			smax
	)		((s!=sp)?'0':DWT_Quadrant1D(i,s,smax))

Definition at line 244 of file wavelet.cpp.

DWT_Scale

#define DWT_Scale	(		i,
			smax
	)		((int)((i==0)?smax-1:(ABS((CEIL(log10((double)(i+1))/log10(2.0))-smax)))))

Definition at line 242 of file wavelet.cpp.

Function Documentation

adaptive_soft_thresholding_block2D()

void adaptive_soft_thresholding_block2D	(	MultidimArray< double > &	I,
		int	scale,
		const std::string &	quadrant,
		double	sigma
	)

Definition at line 321 of file wavelet.cpp.

{
    // Compute block variance
    Matrix1D<int> corner1(2), corner2(2);
    Matrix1D<double> r(2);
    SelectDWTBlock(scale, I, quadrant,
                   XX(corner1), XX(corner2), YY(corner1), YY(corner2));
    double dummy, avg, stddev;
    I.computeStats(avg, stddev, dummy, dummy, corner1, corner2);

    // Now denoise
    double th = sigma * sigma / stddev;
    FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(corner1, corner2)
    {
        if (ABS(I(r)) > th)
            if (I(r) > 0)
                I(r) -= th;
            else
                I(r) += th;
        else
            I(r) = 0;
    }
}

bayesian_solve_eq_system()

void bayesian_solve_eq_system	(	const Matrix1D< double > &	power,
		const Matrix1D< double > &	Ncoefs,
		double	SNR0,
		double	SNRF,
		double	powerI,
		double	power_rest,
		bool	white_noise,
		Matrix1D< double > &	estimatedS
	)

Definition at line 476 of file wavelet.cpp.

{

    int scale_dim = power.size();

    Matrix2D<double> A;
    int extra_constraints = 0;
    if (white_noise)
        extra_constraints += (scale_dim - 1);
    A.initZeros(2*(scale_dim - 1) + 2*scale_dim + 2 + extra_constraints, 2*scale_dim);
    for (int i = 1;i < scale_dim;i++)
    {
        A(i - 1, i - 1) = 1;
        A(i - 1, i) = -1;
        A(i - 1 + scale_dim - 1, i - 1 + scale_dim) = 1;
        A(i - 1 + scale_dim - 1, i + scale_dim) = -1;
    }
    for (int i = 0;i < 2*scale_dim;i++)
        A(i + 2*(scale_dim - 1), i) = -1;

    // Constraints on the SNR
    Matrix1D<double> aux0coefs(scale_dim);
    for (int j = 0;j < scale_dim;j++)
        aux0coefs(j) = Ncoefs(j) * SNR0;
    Matrix1D<double> auxFcoefs(scale_dim);
    for (int j = 0;j < scale_dim;j++)
        auxFcoefs(j) = Ncoefs(j) * SNRF;

    //initializing the second last row of A
    for (int j = 0;j < scale_dim;j++)
        A(2*(scale_dim - 1) + 2*scale_dim, j) = (-1) * auxFcoefs(j);
    for (int j = scale_dim;j < 2*scale_dim;j++)
        A(2*(scale_dim - 1) + 2*scale_dim, j) = Ncoefs(j - scale_dim);

    //initializing the last row of A
    for (int j = 0;j < scale_dim;j++)
        A(2*(scale_dim - 1) + 2*scale_dim + 1, j) = aux0coefs(j);
    for (int j = scale_dim;j < 2*scale_dim;j++)
        A(2*(scale_dim - 1) + 2*scale_dim + 1, j) = (-1) * Ncoefs(j - scale_dim);

    // White noise constraints
    if (white_noise)
        for (int i = 0; i < scale_dim - 1; i++)
        {
            A(A.mdimy - (scale_dim - 1) + i, i)  = -1.01;
            A(A.mdimy - (scale_dim - 1) + i, i + 1) = 1;
        }

    //initialize the matrix b
    Matrix1D<double> b(MAT_YSIZE(A));

    // Initialize Aeq matrix
    Matrix2D<double> Aeq;
    Aeq.initZeros(1, 2*scale_dim);
    for (int j = 0;j < scale_dim;j++)
    {
        Aeq(0, j) = Ncoefs(j);
        Aeq(0, j + scale_dim) = Ncoefs(j);
    }

    //initialize beq matrix
    Matrix1D<double> beq;
    beq.initZeros(1);
    beq(0) = powerI - power_rest;

    //initialization of Matrix C (cost matrix)
    Matrix2D<double> C;
    C.initZeros(scale_dim, 2*scale_dim);
    for (int j = 0;j < scale_dim;j++)
    {
        C(j, j) = 1;
        C(j, j + scale_dim) = 1;
    }

    // initialise the estimatedS which will contain the solution vector
    estimatedS.initZeros(2*scale_dim);
#ifdef DEBUG
    //Writing the matrices to ASCII files for comparing with matlab
    C.write("./matrices/C.txt");
    power.write("./matrices/power.txt");
    A.write("./matrices/A.txt");
    b.write("./matrices/b.txt");
    Aeq.write("./matrices/Aeq.txt");
    beq.write("./matrices/beq.txt");

    std::cout << "Equation system Cx=d\n"
    << "C=\n" << C << std::endl
    << "d=" << (power / Ncoefs).transpose() << std::endl
    << "Constraints\n"
    << "Ax<=b\n"
    << "A=\n" << A << std::endl
    << "b=" << b.transpose() << std::endl
    << "Aeq x=beq\n"
    << "Aeq=\n" << Aeq << std::endl
    << "beq=" << beq.transpose() << std::endl;
#endif

    // Solve the system
    Matrix1D<double> bl, bu;
    leastSquare(C, power / Ncoefs, A, b, Aeq, beq, bl, bu, estimatedS);
    // COSS
    estimatedS /= 2;

#ifdef DEBUG

    std::cout << "scale_dim :: " << scale_dim << std::endl;
    std::cout << "--------estimatedS -------- \n";
    std::cout << estimatedS;
    std::cout << "--------------------------- \n";
    std::cout << "Inequality constraints agreement" << std::endl
    << (A*estimatedS).transpose() << std::endl;
    std::cout << "Equality constraints agreement" << std::endl
    << (Aeq*estimatedS).transpose() << std::endl;
    std::cout << "Goal function value: " << (C*estimatedS).transpose() << std::endl;
#endif
}

compute_noise_power()

double compute_noise_power

(

MultidimArray< double > &

I

)

Definition at line 346 of file wavelet.cpp.

{
    // Compute histogram of the absolute values of the DWT coefficients
    // at scale=0
    Histogram1D hist;
    double avg, stddev, min_val=0., max_val=0.;
    I.computeStats(avg, stddev, min_val, max_val);
    hist.init(0, XMIPP_MAX(ABS(min_val), ABS(max_val)), 100);

    Matrix1D<int> corner1(2), corner2(2);
    Matrix1D<double> r(2);
    SelectDWTBlock(0, I, "01",
                   XX(corner1), XX(corner2), YY(corner1), YY(corner2));
    FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(corner1, corner2)
    {
        double value=fabs(I(r));
        INSERT_VALUE(hist,value);
    }

    SelectDWTBlock(0, I, "10",
                   XX(corner1), XX(corner2), YY(corner1), YY(corner2));
    FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(corner1, corner2)
    {
        double value=fabs(I(r));
        INSERT_VALUE(hist,value);
    }

    SelectDWTBlock(0, I, "11",
                   XX(corner1), XX(corner2), YY(corner1), YY(corner2));
    FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(corner1, corner2)
    {
        double value=fabs(I(r));
        INSERT_VALUE(hist,value);
    }

    return hist.percentil(50) / 0.6745;
}

DWT_Bijaoui_denoise_LL2D()

void DWT_Bijaoui_denoise_LL2D	(	MultidimArray< double > &	WI,
		int	scale,
		const std::string &	orientation,
		double	mu,
		double	S,
		double	N
	)

Definition at line 414 of file wavelet.cpp.

{
    Matrix1D<int> x0(2), xF(2), r(2);
    SelectDWTBlock(scale, WI, orientation, XX(x0), XX(xF), YY(x0), YY(xF));

    double SN = S + N;
    double S_N = S / SN;
    if (S < 1e-6 && N < 1e-6)
    {
        FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(x0, xF)
        A2D_ELEM(WI,YY(r),XX(r)) = 0;
    }
    else
    {
        double iSN=1.0/SN;
        double iN=1.0/N;
        FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(x0, xF)
        {
            double &WI_r=A2D_ELEM(WI,YY(r),XX(r));
            double y = WI_r;
            double ymu = y - mu;
            double ymu2 = -0.5 * ymu * ymu;
            double expymu2SN = exp(ymu2 * iSN);
            double den = exp(ymu2 * iN) + expymu2SN;
            if (den > 1e-10)
                WI_r *= S_N * expymu2SN / den;
        }
    }
}

DWT_Bijaoui_denoise_LL3D()

void DWT_Bijaoui_denoise_LL3D	(	MultidimArray< double > &	WI,
		int	scale,
		const std::string &	orientation,
		double	mu,
		double	S,
		double	N
	)

Definition at line 446 of file wavelet.cpp.

{
    Matrix1D<int> x0(3), xF(3), r(3);
    SelectDWTBlock(scale, WI, orientation, XX(x0), XX(xF), YY(x0), YY(xF),
                   ZZ(x0), ZZ(xF));

    double SN = S + N;
    double S_N = S / SN;
    if (S < 1e-6 && N < 1e-6)
    {
        FOR_ALL_ELEMENTS_IN_ARRAY3D_BETWEEN(x0, xF) WI(r) = 0;
    }
    else
    {
        FOR_ALL_ELEMENTS_IN_ARRAY3D_BETWEEN(x0, xF)
        {
            double y = WI(r);
            double ymu = y - mu;
            double ymu2 = -0.5 * ymu * ymu;
            double expymu2SN = exp(ymu2 / SN);
            double den = exp(ymu2 / N) + expymu2SN;
            if (den > 1e-10)
                WI(r) = S_N * expymu2SN / den * y;
        }
    }
}