Wavelets

Collaboration diagram for Wavelets:

Macros
#define	HAAR   2
#define	DAUB4   4
#define	DAUB12   12
#define	DAUB20   20

Wavelet transform
void	Bilib_DWT (const MultidimArray< double > &input, MultidimArray< double > &result, int iterations, int isign=1)
void	set_DWT_type (int DWT_type)
int	Get_Max_Scale (int size)
template<typename T >
void	DWT (const MultidimArray< T > &v, MultidimArray< double > &result, int isign=1)
void	IDWT (const MultidimArray< double > &v, MultidimArray< double > &result)

Wavelet related functions
void	DWT_lowpass2D (const MultidimArray< double > &v, MultidimArray< double > &result)
template<typename T >
void	SelectDWTBlock (int scale, const MultidimArray< T > &I, const std::string &quadrant, int &x1, int &x2)
template<typename T >
void	SelectDWTBlock (int scale, const MultidimArray< T > &I, const std::string &quadrant, int &x1, int &x2, int &y1, int &y2)
template<typename T >
void	SelectDWTBlock (int scale, const MultidimArray< T > &I, const std::string &quadrant, int &x1, int &x2, int &y1, int &y2, int &z1, int &z2)
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)
#define	DWT_Imin(s, smax, l)
#define	DWT_Imax(s, smax, l)

Wavelet Denoising
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_thresholding2D (MultidimArray< double > &I, int scale)
void	DWT_keep_central_part (MultidimArray< double > &I, double R)
Matrix1D< double >	bayesian_wiener_filtering2D (MultidimArray< double > &WI, int allowed_scale, double SNR0=0.1, double SNRF=0.2, bool white_noise=false, int tell=0, bool denoise=true)
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=0.1, double SNRF=0.2, bool white_noise=false, int tell=0, bool denoise=true)
void	bayesian_wiener_filtering3D (MultidimArray< double > &WI, int allowed_scale, Matrix1D< double > &estimatedS)
void	phaseCongMono (MultidimArray< double > &I, MultidimArray< double > &PC, MultidimArray< double > &FT, MultidimArray< double > &Energy, MultidimArray< double > &lowPass, MultidimArray< double > &Radius, MultidimArray< std::complex< double > > &H, int nScale, double minWaveLength, double mult, double sigmaOnf)

Detailed Description

Macro Definition Documentation

DAUB12

#define DAUB12   12

Definition at line 39 of file wavelet.h.

DAUB20

#define DAUB20   20

Definition at line 40 of file wavelet.h.

DAUB4

#define DAUB4   4

Definition at line 38 of file wavelet.h.

DWT_Imax

#define DWT_Imax	(		s,
			smax,
			l
	)

Value:

static_cast< int>(((l == '0') ? \
        pow(2.0, smax - s - 1) - 1 : pow(2.0, smax - s) - 1))

Definition at line 124 of file wavelet.h.

DWT_Imin

#define DWT_Imin	(		s,
			smax,
			l
	)

Value:

static_cast< int >(((l == '0') ? 0 : \
        pow(2.0, smax - s - 1)))

Definition at line 122 of file wavelet.h.

HAAR

#define HAAR   2

Definition at line 37 of file wavelet.h.

Function Documentation

adaptive_soft_thresholding2D()

void adaptive_soft_thresholding2D	(	MultidimArray< double > &	I,
		int	scale
	)

Adaptive soft thresholding 2D.

Chang, Yu, Betterli. IEEE Int. Conf. Image Processing

Definition at line 384 of file wavelet.cpp.

{
    double sigma = compute_noise_power(I);
    for (int s = 0; s <= scale; s++)
    {
        adaptive_soft_thresholding_block2D(I, s, "01", sigma);
        adaptive_soft_thresholding_block2D(I, s, "10", sigma);
        adaptive_soft_thresholding_block2D(I, s, "11", sigma);
    }
}

bayesian_wiener_filtering2D() [1/2]

Matrix1D< double > bayesian_wiener_filtering2D	(	MultidimArray< double > &	WI,
		int	allowed_scale,
		double	SNR0 = 0.1,
		double	SNRF = 0.2,
		bool	white_noise = false,
		int	tell = 0,
		bool	denoise = true
	)

Bayesian, Wiener filtering.

Bijaoui, Signal Processing 2002, 82: 709-712. The Denoising procedure is applied up to the scale given. SNR0 is the smallest SNR and SNRF is the largest SNR.

This function returns the estimated coefficients for S and N at each scale. If denoise is set to false, then S and N coefficients are estimated but they are not applied to the image.

Definition at line 603 of file wavelet.cpp.

{
    /*Calculate the power of the wavelet transformed image */
    double powerI = WI.sum2();

    /*Number of pixels and some constraints on SNR*/
    int Xdim = XSIZE(WI);
    int max_scale = ROUND(log(double(Xdim)) / log(2.0));

#ifdef DEBUG

    std::cout << "powerI= " << powerI << std::endl;
    double powerWI = WI.sum2();
    std::cout << "powerWI= " << powerWI << std::endl;
    std::cout << "Ydim = " << Ydim << "  Xdim = " << Xdim << "\n";
    std::cout << "Ncoef_total= " << Ncoef_total << std::endl;
    std::cout << "max_scale = " << max_scale << "\n";
#endif

    /*Calculate the power at each band*/
    //init the scale vector
    Matrix1D<int> scale(XMIPP_MIN(allowed_scale + 1, max_scale - 1));
    FOR_ALL_ELEMENTS_IN_MATRIX1D(scale) scale(i) = i;
    int scale_dim = scale.size();

    //define some vectors
    Matrix1D<double> power(scale_dim), average(scale_dim), Ncoefs(scale_dim);
    Matrix1D<int> x0(2), xF(2), r(2);
    std::vector<std::string> orientation;
    orientation.push_back("01");
    orientation.push_back("10");
    orientation.push_back("11");
    int orientationSize=orientation.size();
    int jmax=scale.size();
    for (int j = 0;j < jmax;j++)
    {
        for (size_t k = 0; k < orientation.size(); k++)
        {
            SelectDWTBlock(scale(j), WI, orientation[k],
                           XX(x0), XX(xF), YY(x0), YY(xF));
            FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(x0, xF)
            {
                double aux=DIRECT_A2D_ELEM(WI,YY(r),XX(r));
                VEC_ELEM(power,j) += aux*aux;
                VEC_ELEM(average,j) += aux;
            }
        }
        VEC_ELEM(Ncoefs,j) = (int)pow(2.0, 2 * (max_scale - VEC_ELEM(scale,j) - 1)) * orientationSize;
        VEC_ELEM(average,j) = VEC_ELEM(average,j) / VEC_ELEM(Ncoefs,j);
    }

    /*Evaluate the power of the unconsidered part of the image */
    double power_rest = 0.0;
    int Ncoefs_rest = 0;
    SelectDWTBlock(scale(scale_dim - 1), WI, "00", XX(x0), XX(xF), YY(x0), YY(xF));
    FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(x0, xF)
    {
        double aux=DIRECT_A2D_ELEM(WI,YY(r),XX(r));
        power_rest += aux*aux;
    }
    Ncoefs_rest = (int)pow(2.0, 2 * (max_scale - 1 - scale(scale_dim - 1)));

    if (tell)
    {
        std::cout << "scale = " << std::endl << scale << std::endl;
        std::cout << "power= " << std::endl << power << "\n";
        std::cout << "average= " << std::endl << average << "\n";
        std::cout << "Ncoefs= " << std::endl << Ncoefs << "\n";
        std::cout << "power_rest= " << power_rest << "\n";
        std::cout << "Ncoefs_rest= " << Ncoefs_rest << "\n";
        std::cout << "powerI= " << powerI << std::endl;
        std::cout << "Total sum of powers = " << power.sum() + power_rest << std::endl;
        std::cout << "Difference powers = " << powerI - power.sum() - power_rest << std::endl;
    }

    /*Solve the Equation System*/
    Matrix1D<double> estimatedS;
    bayesian_solve_eq_system(power, Ncoefs,
                             SNR0, SNRF, powerI, power_rest, white_noise, estimatedS);

    if (tell)
    {
        std::cout << "estimatedS =\n" << estimatedS << std::endl;
        double S = 0, N = 0;
        for (int i = 0; i < scale_dim; i++)
        {
            N += Ncoefs(i) * estimatedS(i);
            S += Ncoefs(i) * estimatedS(scale_dim + i);
        }
        std::cout << "SNR value=" << S / N << std::endl << std::endl;
    }

    /* Apply the Bijaoui denoising to all scales >= allowed_scale */
    if (denoise)
        bayesian_wiener_filtering2D(WI, allowed_scale, estimatedS);

    return estimatedS;
}

bayesian_wiener_filtering2D() [2/2]

void bayesian_wiener_filtering2D	(	MultidimArray< double > &	WI,
		int	allowed_scale,
		Matrix1D< double > &	estimatedS
	)

Bayesian, Wiener filtering.

This is the function that really denoise.

Definition at line 704 of file wavelet.cpp.

{
    std::vector<std::string> orientation;
    orientation.push_back("01");
    orientation.push_back("10");
    orientation.push_back("11");

    int max_scale = ROUND(log(double(XSIZE(WI))) / log(2.0));
    Matrix1D<int> scale(XMIPP_MIN(allowed_scale + 1, max_scale - 1));
    FOR_ALL_ELEMENTS_IN_MATRIX1D(scale) scale(i) = i;

    for (size_t i = 0;i < scale.size();i++)
    {
        double N = estimatedS(i);
        double S = estimatedS(i + scale.size());
        for (size_t k = 0; k < orientation.size(); k++)
            DWT_Bijaoui_denoise_LL2D(WI, scale(i), orientation[k], 0, S, N);
    }
}

bayesian_wiener_filtering3D() [1/2]

Matrix1D< double > bayesian_wiener_filtering3D	(	MultidimArray< double > &	WI,
		int	allowed_scale,
		double	SNR0 = 0.1,
		double	SNRF = 0.2,
		bool	white_noise = false,
		int	tell = 0,
		bool	denoise = true
	)

Bayesian, Wiener filtering.

Bijaoui, Signal Processing 2002, 82: 709-712. The denoising procedure is applied up to the scale given. SNR0 is the smallest SNR and SNRF is the largest SNR.

This function returns the estimated coefficients for S and N at each scale. If denoise is set to false, then S and N coefficients are estimated but they are not applied to the image.

Definition at line 726 of file wavelet.cpp.

{
    /*Calculate the power of the wavelet transformed image */
    double powerI = WI.sum2();

    /*Number of pixels and some constraints on SNR*/
    size_t Xdim=ZSIZE(WI);
    int max_scale = ROUND(log(double(Xdim)) / log(2.0));

#ifdef DEBUG

    std::cout << "powerI= " << powerI << std::endl;
    double powerWI = WI.sum2();
    std::cout << "powerWI= " << powerWI << std::endl;
    std::cout << "Zdim= " << Zdim << " Ydim = " << Ydim << "  Xdim = " << Xdim << "\n";
    std::cout << "Ncoef_total= " << Ncoef_total << std::endl;
    std::cout << "max_scale = " << max_scale << "\n";
#endif

    /*Calculate the power at each band*/
    //init the scale vector
    Matrix1D<int> scale(allowed_scale + 1);
    FOR_ALL_ELEMENTS_IN_MATRIX1D(scale) scale(i) = i;
    int scale_dim = scale.size();

    //define some vectors
    Matrix1D<double> power(scale_dim), average(scale_dim), Ncoefs(scale_dim);
    Matrix1D<int> x0(3), xF(3), r(3);
    std::vector<std::string> orientation;
    orientation.push_back("001");
    orientation.push_back("010");
    orientation.push_back("011");
    orientation.push_back("100");
    orientation.push_back("101");
    orientation.push_back("110");
    orientation.push_back("111");
    for (size_t j = 0;j < scale.size();j++)
    {
        for (size_t k = 0; k < orientation.size(); k++)
        {
            SelectDWTBlock(scale(j), WI, orientation[k],
                           XX(x0), XX(xF), YY(x0), YY(xF), ZZ(x0), ZZ(xF));
            FOR_ALL_ELEMENTS_IN_ARRAY3D_BETWEEN(x0, xF)
            {
                power(j) += WI(r) * WI(r);
                average(j) += WI(r);
            }
        }
        Ncoefs(j) = (int)pow(2.0, 3 * (max_scale - scale(j) - 1)) * orientation.size();
        average(j) = average(j) / Ncoefs(j);
    }

    /*Evaluate the power of the unconsidered part of the image */
    double power_rest = 0.0;
    int Ncoefs_rest = 0;
    SelectDWTBlock(scale(scale_dim - 1), WI, "000", XX(x0), XX(xF), YY(x0), YY(xF),
                   ZZ(x0), ZZ(xF));
    FOR_ALL_ELEMENTS_IN_ARRAY3D_BETWEEN(x0, xF)
    power_rest += WI(r) * WI(r);
    Ncoefs_rest = (int)pow(2.0, 3 * (max_scale - 1 - scale(scale_dim - 1)));

    if (tell)
    {
        std::cout << "scale = " << std::endl << scale << std::endl;
        std::cout << "power= " << std::endl << power << "\n";
        std::cout << "average= " << std::endl << average << "\n";
        std::cout << "Ncoefs= " << std::endl << Ncoefs << "\n";
        std::cout << "power_rest= " << power_rest << "\n";
        std::cout << "Ncoefs_rest= " << Ncoefs_rest << "\n";
        std::cout << "powerI= " << powerI << std::endl;
        std::cout << "Total sum of powers = " << power.sum() + power_rest << std::endl;
        std::cout << "Difference powers = " << powerI - power.sum() - power_rest << std::endl;
    }

    /*Solve the Equation System*/
    Matrix1D<double> estimatedS;
    bayesian_solve_eq_system(power, Ncoefs,
                             SNR0, SNRF, powerI, power_rest, white_noise, estimatedS);
    if (tell)
    {
        std::cout << "estimatedS =\n" << estimatedS << std::endl;
        double S = 0, N = 0;
        for (int i = 0; i < scale_dim; i++)
        {
            N += Ncoefs(i) * estimatedS(i);
            S += Ncoefs(i) * estimatedS(scale_dim + i);
        }
        std::cout << "SNR value=" << S / N << std::endl << std::endl;
    }

    /* Apply the Bijaoui denoising to all scales >= allowed_scale */
    if (denoise)
        bayesian_wiener_filtering3D(WI, allowed_scale, estimatedS);
    return estimatedS;
}

bayesian_wiener_filtering3D() [2/2]

void bayesian_wiener_filtering3D	(	MultidimArray< double > &	WI,
		int	allowed_scale,
		Matrix1D< double > &	estimatedS
	)

Bayesian, Wiener filtering.

This is the function that really denoise.

Definition at line 825 of file wavelet.cpp.

{
    std::vector<std::string> orientation;
    orientation.push_back("001");
    orientation.push_back("010");
    orientation.push_back("011");
    orientation.push_back("100");
    orientation.push_back("101");
    orientation.push_back("110");
    orientation.push_back("111");

    int max_scale = ROUND(log(double(XSIZE(WI))) / log(2.0));
    MultidimArray<int> scale(XMIPP_MIN(allowed_scale + 1, max_scale - 1));
    FOR_ALL_ELEMENTS_IN_ARRAY1D(scale) scale(i) = i;

    for (size_t i = 0;i < XSIZE(scale);i++)
    {
        double N = estimatedS(i);
        double S = estimatedS(i + XSIZE(scale));
        for (size_t k = 0; k < orientation.size(); k++)
            DWT_Bijaoui_denoise_LL3D(WI, scale(i), orientation[k], 0, S, N);
    }
}

Bilib_DWT()

void Bilib_DWT	(	const MultidimArray< double > &	input,
		MultidimArray< double > &	result,
		int	iterations,
		int	isign = 1
	)

B-spline Wavelet transform of a vector.

The B-spline wavelet transform of the input array is computed. The size of the array must be so as to allow the downsampling by 2 as many times as the number of iterations. For instance, if iterations is 1, then it must be a multiple of 2. If iterations is 2, then it must be a multiple of 4. If iterations if 3, then it must be a multiple of 8. And so on.

If the isign=-1 then the inverse wavelet transform is performed.

Definition at line 43 of file wavelet.cpp.

{
    if (iterations < 1)
        REPORT_ERROR(ERR_INDEX_OUTOFBOUNDS, "Bilib_DWT: iterations must be >=1");
    auto size_multiple = (int)pow(2.0, (double) iterations);
    if (XSIZE(input) % size_multiple != 0)
        REPORT_ERROR(ERR_MULTIDIM_SIZE,
                     (std::string)"Bilib_DWT: Xsize must be a multiple of " +
                     integerToString(size_multiple));
    if (YSIZE(input) > 1 && YSIZE(input) % size_multiple != 0)
        REPORT_ERROR(ERR_MULTIDIM_SIZE,
                     (std::string)"Bilib_DWT 2D: Ysize must be a multiple of " +
                     integerToString(size_multiple));
    if (ZSIZE(input) > 1 && ZSIZE(input) % size_multiple != 0)
        REPORT_ERROR(ERR_MULTIDIM_SIZE,
                     (std::string)"Bilib_DWT 3D: Zsize must be a multiple of " +
                     integerToString(size_multiple));

    result.initZeros(input);
    TWaveletStruct TW;
    if (isign == 1)
        TW.Operation = "Analysis";
    else if (isign == -1)
        TW.Operation = "Synthesis";
    else
        REPORT_ERROR(ERR_VALUE_INCORRECT, (std::string)"waveletTransform: unrecognized isign");
    TW.Filter = "Orthonormal Spline";
    TW.BoundaryConditions = "Mirror";
    TW.Order = "1";
    TW.Alpha = 0;

    if (isign == 1)
    {
        // First iteration
        TW.Input = MULTIDIM_ARRAY(input);
        TW.NxOutput = TW.NxInput = XSIZE(input);
        TW.NyOutput = TW.NyInput = YSIZE(input);
        TW.NzOutput = TW.NzInput = ZSIZE(input);
        TW.Output = MULTIDIM_ARRAY(result);
        Wavelet(&TW);

        // Subsequent iterations
        for (int i = 1; i < iterations; i++)
        {
            // Size of the Lowest subband
            int xsize = XMIPP_MAX(1, XSIZE(input) / (int)pow(2.0, (double)i));
            int ysize = XMIPP_MAX(1, YSIZE(input) / (int)pow(2.0, (double)i));
            int zsize = XMIPP_MAX(1, ZSIZE(input) / (int)pow(2.0, (double)i));

            // Pick the Lowest subband
            MultidimArray<double> input_aux, result_aux;
            input_aux.resize(zsize, ysize, xsize);
            result_aux.resize(zsize, ysize, xsize);
            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(input_aux)
                DIRECT_A3D_ELEM(input_aux, k, i, j) = DIRECT_A3D_ELEM(result, k, i, j);

            // DWT
            TW.Input = MULTIDIM_ARRAY(input_aux);
            TW.NxOutput = TW.NxInput = XSIZE(input_aux);
            TW.NyOutput = TW.NyInput = YSIZE(input_aux);
            TW.NzOutput = TW.NzInput = ZSIZE(input_aux);
            TW.Output = MULTIDIM_ARRAY(result_aux);
            Wavelet(&TW);

            // Return the subband to the output
            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(input_aux)
                DIRECT_A3D_ELEM(result, k, i, j) = DIRECT_A3D_ELEM(result_aux, k, i, j);
        }
    }
    else if (isign == -1)
    {
        // Subsequent iterations
        for (int i = iterations - 1; i >= 1; i--)
        {
            // Size of the Lowest subband
            int xsize = XMIPP_MAX(1, XSIZE(input) / (int)pow(2.0, (double)i));
            int ysize = XMIPP_MAX(1, YSIZE(input) / (int)pow(2.0, (double)i));
            int zsize = XMIPP_MAX(1, ZSIZE(input) / (int)pow(2.0, (double)i));

            // Pick the Lowest subband
            MultidimArray<double> input_aux, result_aux;
            input_aux.resize(zsize, ysize, xsize);
            result_aux.resize(zsize, ysize, xsize);
            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(input_aux)
                DIRECT_A3D_ELEM(input_aux, k, i, j) = DIRECT_A3D_ELEM(input, k, i, j);

            // DWT
            TW.Input = MULTIDIM_ARRAY(input_aux);
            TW.NxOutput = TW.NxInput = XSIZE(input_aux);
            TW.NyOutput = TW.NyInput = YSIZE(input_aux);
            TW.NzOutput = TW.NzInput = ZSIZE(input_aux);
            TW.Output = MULTIDIM_ARRAY(result_aux);
            Wavelet(&TW);

            // Return the subband to the output
            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(input_aux)
                DIRECT_A3D_ELEM(input, k, i, j) = DIRECT_A3D_ELEM(result_aux, k, i, j);
        }

        // First iteration
        TW.Input = MULTIDIM_ARRAY(input);
        TW.NxOutput = TW.NxInput = XSIZE(input);
        TW.NyOutput = TW.NyInput = YSIZE(input);
        TW.NzOutput = TW.NzInput = ZSIZE(input);
        TW.Output = MULTIDIM_ARRAY(result);
        Wavelet(&TW);
    }
}

clean_quadrant2D()

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

Remove all information within a quadrant and scale.

Definition at line 288 of file wavelet.cpp.

{
    Matrix1D<int> corner1(2), corner2(2);
    Matrix1D<double> r(2);
    SelectDWTBlock(scale, I, quadrant, XX(corner1), XX(corner2), YY(corner1), YY(corner2));
    FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(corner1, corner2)
        I(r) = 0;
}

clean_quadrant3D()

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

Remove all information within a quadrant and scale.

Definition at line 297 of file wavelet.cpp.

{
    int x1, y1, z1, x2, y2, z2;
    SelectDWTBlock(scale, I, quadrant, x1, x2, y1, y2, z1, z2);
    Matrix1D<int> corner1(3), corner2(3);
    Matrix1D<double> r(3);
    SelectDWTBlock(scale, I, quadrant, XX(corner1), XX(corner2),
                   YY(corner1), YY(corner2), ZZ(corner1), ZZ(corner2));
    FOR_ALL_ELEMENTS_IN_ARRAY3D_BETWEEN(corner1, corner2) I(r) = 0;
}

DWT()

template<typename T >

void DWT	(	const MultidimArray< T > &	v,
		MultidimArray< double > &	result,
		int	isign = 1
	)

DWT of a MultidimArray

The output vector can be the same as the input one. Previously the type of DWT must be set with set_DWT_type. If isign=1 the direct DWT is performed, if isign=-1 the inverse DWT is done.

Definition at line 83 of file wavelet.h.

{
    unsigned long int nn[3];
    unsigned long int *ptr_nn = nn - 1;
    int dim;
    nn[2] = ZSIZE(v);
    nn[1] = YSIZE(v);
    nn[0] = XSIZE(v);

    if (YSIZE(v) == 1 && ZSIZE(v) == 1)
        dim = 1;
    else if (ZSIZE(v) == 1)
        dim = 2;
    else
        dim = 3;

    typeCast(v, result);
    double* ptr_result = MULTIDIM_ARRAY(result) - 1;
    wtn(ptr_result, ptr_nn, dim, isign, pwt);
}

DWT_keep_central_part()

void DWT_keep_central_part	(	MultidimArray< double > &	I,
		double	R
	)

Keep central part 2D.

Keep those coefficients in a certain radius.

Definition at line 396 of file wavelet.cpp.

{
    Mask mask(INT_MASK);
    mask.type = BINARY_DWT_CIRCULAR_MASK;
    if (R == -1)
        mask.R1 = (double)XSIZE(I) / 2 + 1;
    else
        mask.R1 = R;
    mask.smin = 0;
    mask.smax = Get_Max_Scale(XSIZE(I));
    mask.quadrant = "xx";
    mask.resize(I);
    mask.generate_mask();
    mask.apply_mask(I, I);
}

DWT_lowpass2D()

void DWT_lowpass2D	(	const MultidimArray< double > &	v,
		MultidimArray< double > &	result
	)

DWT Low pass versions.

This function returns the low pass versions at different scales. The low pass version of the image at scale s is stored in the 01 quadrant of that scale.

Definition at line 165 of file wavelet.cpp.

{
    MultidimArray<double> dwt, aux;
    result.initZeros(YSIZE(v), XSIZE(v) / 2);
    DWT(v, dwt);
    int Nx = Get_Max_Scale(XSIZE(v));
    for (int s = 0; s < Nx; s++)
    {
        // Perform the inverse DWT transform of the low pass
        dwt.resize(XSIZE(dwt) / 2, YSIZE(dwt) / 2);
        IDWT(dwt, aux);
        // Copy the result to the 01 quadrant of the result
        int x1, y1, x2, y2, x, y, i, j;
        SelectDWTBlock(s, v, "01", x1, x2, y1, y2);
        for (y = y1, i = 0; y <= y2; y++, i++)
            for (x = x1, j = 0; x <= x2; x++, j++)
                A2D_ELEM(result, y, x) = A2D_ELEM(aux, i, j);
    }
}

Get_Max_Scale()

int Get_Max_Scale ( int  size )

inline

Get maximum scale.

This function returns the maximum scale achievable by the DWT transform of a given size.

Definition at line 71 of file wavelet.h.

{
    return ROUND(log10(static_cast< double >(size)) / log10(2.0));
}

Get_Scale_Quadrant() [1/3]

void Get_Scale_Quadrant	(	int	size_x,
		int	x,
		int &	scale,
		std::string &	quadrant
	)

Get scale and quadrant 1D.

Given a point and the maximum size of the image, this routine returns the scale and quadrant it belongs.

Definition at line 248 of file wavelet.cpp.

{
    double Nx = Get_Max_Scale(size_x);
    quadrant = "x";
    scale = DWT_Scale(x, Nx);
    quadrant[0] = DWT_Quadrant1D(x, scale, Nx);
}

Get_Scale_Quadrant() [2/3]

void Get_Scale_Quadrant	(	int	size_x,
		int	size_y,
		int	x,
		int	y,
		int &	scale,
		std::string &	quadrant
	)

Get scale and quadrant 2D.

Given a point and the maximum size of the image, this routine returns the scale and quadrant it belongs.

Definition at line 257 of file wavelet.cpp.

{
    double Nx = Get_Max_Scale(size_x);
    double Ny = Get_Max_Scale(size_y);
    quadrant = "xy";
    double scalex = DWT_Scale(x, Nx);
    double scaley = DWT_Scale(y, Ny);
    scale = (int)(XMIPP_MIN(scalex, scaley));
    quadrant[1] = DWT_QuadrantnD(y, scaley, scale, Ny);
    quadrant[0] = DWT_QuadrantnD(x, scalex, scale, Nx);
}

Get_Scale_Quadrant() [3/3]

void Get_Scale_Quadrant	(	int	size_x,
		int	size_y,
		int	size_z,
		int	x,
		int	y,
		int	z,
		int &	scale,
		std::string &	quadrant
	)

Get scale and quadrant 3D.

Given a point and the maximum size of the image, this routine returns the scale and quadrant it belongs.

Definition at line 270 of file wavelet.cpp.

{
    double Nx = Get_Max_Scale(size_x);
    double Ny = Get_Max_Scale(size_y);
    double Nz = Get_Max_Scale(size_z);
    quadrant = "xyz";
    double scalex = DWT_Scale(x, Nx);
    double scaley = DWT_Scale(y, Ny);
    double scalez = DWT_Scale(z, Nz);
    scale = (int)(XMIPP_MIN(scalez, XMIPP_MIN(scalex, scaley)));
    quadrant[2] = DWT_QuadrantnD(z, scalez, scale, Nz);
    quadrant[1] = DWT_QuadrantnD(y, scaley, scale, Ny);
    quadrant[0] = DWT_QuadrantnD(x, scalex, scale, Nx);
}

IDWT()

void IDWT	(	const MultidimArray< double > &	v,
		MultidimArray< double > &	result
	)

IDWT of a MultidimArray.

The output volume can be the same as the input one. Previously the type of DWT must be set with set_DWT_type.

Definition at line 159 of file wavelet.cpp.

{
    DWT(v, result, -1);
}

phaseCongMono()

void phaseCongMono	(	MultidimArray< double > &	I,
		MultidimArray< double > &	PC,
		MultidimArray< double > &	FT,
		MultidimArray< double > &	Energy,
		MultidimArray< double > &	lowPass,
		MultidimArray< double > &	Radius,
		MultidimArray< std::complex< double > > &	H,
		int	nScale,
		double	minWaveLength,
		double	mult,
		double	sigmaOnf
	)

Phase congruency filtering.

Phase congruency of an image using monogenic filters.

   I                     - Image to be processed
   PC                    - Phase congruency indicating edge significance
   FT                    - Local weighted mean phase angle at every point in the
                           image.  A value of pi/2 corresponds to a bright line, 0
                         corresponds to a step and -pi/2 is a dark line.
nscale            5    - Number of wavelet scales, try values 3-6
 minWaveLength    3    - Wavelength of smallest scale filter.
 mult             2.1  - Scaling factor between successive filters.
 sigmaOnf         0.55 - Ratio of the standard deviation of the Gaussian
                         describing the log Gabor filter's transfer function
                         in the frequency domain to the filter center frequency.

References:

Peter Kovesi, "Image Features From Phase Congruency". Videre: A Journal of Computer Vision Research. MIT Press. Volume 1, Number 3, Summer 1999 http://mitpress.mit.edu/e-journals/Videre/001/v13.html

Michael Felsberg and Gerald Sommer, "A New Extension of Linear Signal Processing for Estimating Local Properties and Detecting Features". DAGM Symposium 2000, Kiel

Michael Felsberg and Gerald Sommer. "The Monogenic Signal" IEEE Transactions on Signal Processing, 49(12):3136-3144, December 2001

Peter Kovesi, "Phase Congruency Detects Corners and Edges". Proceedings DICTA 2003, Sydney Dec 10-12

Definition at line 850 of file wavelet.cpp.

{

//    #define DEBUG
    double epsilon= .0001; // Used to prevent division by zero.
    //First we set the image origin in the image center
    I.setXmippOrigin();
    //Image size
    size_t NR, NC,NZ, NDim;
    I.getDimensions(NC,NR,NZ,NDim);

    if ( (NZ!=1) || (NDim!=1) )
        REPORT_ERROR(ERR_MULTIDIM_DIM,(std::string)"ZDim and NDim has to be equals to one");

    MultidimArray< std::complex <double> > fftIm;
    // Fourier Transformer
    FourierTransformer ftrans(FFTW_BACKWARD);//, ftransh(FFTW_BACKWARD), ftransf(FFTW_BACKWARD);
    ftrans.FourierTransform(I, fftIm, false);

    if ( (lowPass.xinit == 0) & (lowPass.yinit == 0) & (Radius.xinit == 0) & (Radius.yinit == 0)
          & (H.xinit == 0) & (H.yinit == 0) )
    {

        H.resizeNoCopy(I);
        lowPass.resizeNoCopy(I);
        Radius.resizeNoCopy(I);

        double cutoff = .4;
        double n = 10;
        double wx,wy;

        for (int i=STARTINGY(I); i<=FINISHINGY(I); ++i)
        {
            FFT_IDX2DIGFREQ(i,YSIZE(I),wy);
            double wy2=wy*wy;
            for (int j=STARTINGX(I); j<=FINISHINGX(I); ++j)
            {
                FFT_IDX2DIGFREQ(j,XSIZE(I),wx);

                double wx2=wx*wx;
                double radius=sqrt(wy2+wx2);
                if (radius < 1e-10) radius=1;
                A2D_ELEM(Radius,i,j)=radius;
                auto *ptr=(double*)&A2D_ELEM(H,i,j);
                *ptr=wy/radius;
                *(ptr+1)=wx/radius;
                A2D_ELEM(lowPass,i,j)= 1.0/(1.0+std::pow(radius/cutoff,n));
            }

        }

    }

    MultidimArray<double> logGabor;
    logGabor.resizeNoCopy(I);
    //Bandpassed image in the frequency domain and in the spatial domain.
    //h is a Bandpassed monogenic filtering, real part of h contains
    //convolution result with h1, imaginary part
    //contains convolution result with h2.
    MultidimArray< std::complex <double> > fftImF, fullFftIm, f,Hc,h;

    fftImF.resizeNoCopy(I);
    Hc.resizeNoCopy(I);
    f.resizeNoCopy(I);
    h.resizeNoCopy(I);

    MultidimArray<double> An,AnTemp;
    MultidimArray<double> temp;
    MultidimArray<double> F;
    MultidimArray<double> h1;
    MultidimArray<double> h2;
    MultidimArray<double> tempMat;

    temp.resizeNoCopy(I);
    An.resizeNoCopy(I);
    F.resizeNoCopy(I);
    h1.resizeNoCopy(I);
    h2.resizeNoCopy(I);
    AnTemp.resizeNoCopy(I);
    ftrans.getCompleteFourier(fullFftIm);
    CenterFFT(fullFftIm,false);

    for (int num = 0; num < nScale; ++num)
    {
        double waveLength = minWaveLength;

        for(int numMult=0; numMult < num;numMult++)
            waveLength*=mult;

        double fo = 1.0/waveLength;
        FOR_ALL_ELEMENTS_IN_ARRAY2D(fullFftIm)
        {
            double temp1 =(std::log(A2D_ELEM(Radius,i,j)/fo));
            double temp2 = std::log(sigmaOnf);
            A2D_ELEM(logGabor,i,j) = std::exp(-(temp1*temp1) /(2 * temp2*temp2))*A2D_ELEM(lowPass,i,j);
            if (A2D_ELEM(Radius,i,j)<=1e-10) (A2D_ELEM(logGabor,i,j)=0);

            auto *ptr= (double*)&A2D_ELEM(fullFftIm,i,j);
            temp1 = *ptr;
            temp2 = *(ptr+1);
            ptr= (double*)&A2D_ELEM(fftImF,i,j);
            *ptr=temp1*A2D_ELEM(logGabor,i,j);
            *(ptr+1)=temp2*A2D_ELEM(logGabor,i,j);
        }

        CenterFFT(fftImF,true);
        Hc = fftImF*H;
        ftrans.inverseFourierTransform(fftImF,f);
        ftrans.inverseFourierTransform(Hc,h);

        //ftransf.setFourier(fftImF);
        //ftransf.inverseFourierTransform();
        //ftransh.setFourier(Hc);
        //ftransh.inverseFourierTransform();

        f.getReal(temp);
        F += temp;
        AnTemp = temp*temp;

        h.getReal(temp);
        h1 += temp;
        AnTemp += temp*temp;

        h.getImag(temp);
        h2 += temp;
        AnTemp += temp*temp;

        AnTemp.selfSQRT();
        An += AnTemp;
    }

    Or.resizeNoCopy(I);
    Or.setXmippOrigin();

    Ph.resizeNoCopy(I);
    Ph.setXmippOrigin();

    Energy.resizeNoCopy(I);
    Energy.setXmippOrigin();

    FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
    {
        double temph1 = A2D_ELEM(h1,i,j);
        double temph2 = A2D_ELEM(h2,i,j);
        double tempF =  A2D_ELEM(F,i,j);

        A2D_ELEM(Or,i,j) = std::atan2(temph1,temph2);
        A2D_ELEM(Ph,i,j)=  std::atan2(tempF, std::sqrt(temph1*temph1+
                           temph2*temph2));
        A2D_ELEM(Energy,i,j) = std::sqrt(tempF*tempF+temph1*temph1
                             + temph2*temph2)+epsilon;
    }

#ifdef DEBUG

    FileName fpName1    = "test1.txt";
    FileName fpName2    = "test2.txt";
    FileName fpName3    = "test3.txt";
    FileName fpName4    = "test4.txt";

    Or.write(fpName1);
    Ph.write(fpName2);
    Energy.write(fpName3);

    #endif
}

Quadrant2D()

std::string Quadrant2D

(

int

q

)

Given a quadrant number it returns the string associated to it.

That is nothing more than its corresponding binary representation.

Definition at line 187 of file wavelet.cpp.

{
    switch (q)
    {
    case 0:
        return "00";
        break;
    case 1:
        return "01";
        break;
    case 2:
        return "10";
        break;
    case 3:
        return "11";
        break;
    default:
        REPORT_ERROR(ERR_ARG_INCORRECT,"Unknown quadrant");
    }
}

Quadrant3D()

std::string Quadrant3D

(

int

q

)

Given a quadrant number it returns the string associated to it.

That is nothing more than its corresponding binary representation.

Definition at line 208 of file wavelet.cpp.

{
    switch (q)
    {
    case 0:
        return "000";
        break;
    case 1:
        return "001";
        break;
    case 2:
        return "010";
        break;
    case 3:
        return "011";
        break;
    case 4:
        return "100";
        break;
    case 5:
        return "101";
        break;
    case 6:
        return "110";
        break;
    case 7:
        return "111";
        break;
    default:
        REPORT_ERROR(ERR_ARG_INCORRECT,"Unknown quadrant");
    }
}

SelectDWTBlock() [1/3]

template<typename T >

void SelectDWTBlock	(	int	scale,
		const MultidimArray< T > &	I,
		const std::string &	quadrant,
		int &	x1,
		int &	x2
	)

Select Block 1D.

Given the scale (s=0 is the finest) and the quadrant "0" (Lower frequencies) or "1"(Higher frequencies) this routine returns the indices that should be explored for this block (x1 and x2 should be included in the for).

Definition at line 134 of file wavelet.h.

{
    double Nx = Get_Max_Scale(XSIZE(I));
    I.toLogical(DWT_Imin(scale, Nx, quadrant[0]), x1);
    I.toLogical(DWT_Imax(scale, Nx, quadrant[0]), x2);
}

SelectDWTBlock() [2/3]

template<typename T >

void SelectDWTBlock	(	int	scale,
		const MultidimArray< T > &	I,
		const std::string &	quadrant,
		int &	x1,
		int &	x2,
		int &	y1,
		int &	y2
	)

Select Block 2D.

Given the scale (s=0 is the finest) and the quadrant "xy"="00" (Upper left), "01" (Upper right), "10" (Lower left), "11" (Lower right). This routine returns the indices that should be explored for this block (the extremes should be included in the for).

Definition at line 153 of file wavelet.h.

{
    double Nx = Get_Max_Scale(XSIZE(I));
    double Ny = Get_Max_Scale(YSIZE(I));

    x1 = DWT_Imin(scale, Nx, quadrant[0]);
    y1 = DWT_Imin(scale, Ny, quadrant[1]);
    x2 = DWT_Imax(scale, Nx, quadrant[0]);
    y2 = DWT_Imax(scale, Ny, quadrant[1]);

    I.toLogical(y1, x1, y1, x1);
    I.toLogical(y2, x2, y2, x2);
}

SelectDWTBlock() [3/3]

template<typename T >

void SelectDWTBlock	(	int	scale,
		const MultidimArray< T > &	I,
		const std::string &	quadrant,
		int &	x1,
		int &	x2,
		int &	y1,
		int &	y2,
		int &	z1,
		int &	z2
	)

Select Block 3D.

Given the scale (s=0 is the finest) and the quadrant "xyz"="000", "001", "010", "011", "100", "101", "110", "111". This routine returns the indices that should be explored for this block (the extremes should be included in the for).

Definition at line 182 of file wavelet.h.

{
    double Nx = Get_Max_Scale(XSIZE(I));
    double Ny = Get_Max_Scale(YSIZE(I));
    double Nz = Get_Max_Scale(ZSIZE(I));

    x1 = DWT_Imin(scale, Nx, quadrant[0]);
    y1 = DWT_Imin(scale, Ny, quadrant[1]);
    z1 = DWT_Imin(scale, Nz, quadrant[2]);
    x2 = DWT_Imax(scale, Nx, quadrant[0]);
    y2 = DWT_Imax(scale, Ny, quadrant[1]);
    z2 = DWT_Imax(scale, Nz, quadrant[2]);

    I.toLogical(z1, y1, x1, z1, y1, x1);
    I.toLogical(z2, y2, x2, z2, y2, x2);
}

set_DWT_type()

void set_DWT_type

(

int

DWT_type

)

Set DWT type.

The DWT type should be set before starting making transforms. Valid types are: HAAR, DAUB4, DAUB12, DAUB20

Definition at line 154 of file wavelet.cpp.

{
    pwtset(DWT_type);
}

soft_thresholding()

void soft_thresholding	(	MultidimArray< double > &	I,
		double	th
	)

Soft thresholding .

Subtract a value from all coefficients, if the the value is greater than the absolute value of the coefficient, that coefficient is set to 0.

Definition at line 308 of file wavelet.cpp.

{
    FOR_ALL_ELEMENTS_IN_ARRAY3D(I)
        if (ABS(A3D_ELEM(I, k, i, j)) > th)
        if (A3D_ELEM(I, k, i, j) > 0)
            A3D_ELEM(I, k, i, j) -= th;
        else
            A3D_ELEM(I, k, i, j) += th;
    else
        A3D_ELEM(I, k, i, j) = 0;
}