wavelet.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:     Carlos Oscar S. Sorzano (coss@cnb.csic.es)
 *              Antonio Jose Rodriguez Sanchez (ajr@cnb.csic.es)
 *              Arun Kulshreshth        (arun_2000_iitd@yahoo.com)
 *
 * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 * 02111-1307  USA
 *
 *  All comments concerning this program package may be sent to the
 *  e-mail address 'xmipp@cnb.csic.es'
 ***************************************************************************/

#ifndef CORE_WAVELET_H
#define CORE_WAVELET_H

#include <core/multidim_array.h>
#include <core/numerical_recipes.h>


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



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

void set_DWT_type(int DWT_type);

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

template<typename T>
void DWT(const MultidimArray< T >& v, MultidimArray< double >& result, int isign = 1)
{
    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);
}

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



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

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

template<typename T>
void SelectDWTBlock(int scale,
                    const MultidimArray< T >& I,
                    const std::string& quadrant,
                    int& x1,
                    int& x2)
{
    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);
}

template<typename T>
void SelectDWTBlock(int scale,
                    const MultidimArray< T >& I,
                    const std::string& quadrant,
                    int& x1,
                    int& x2,
                    int& y1,
                    int& y2)
{
    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);
}


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


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


#endif