0
votes

I am trying to implement one of the basic 2D wavelet transform by Haar transformation.

I applied this to the image denoising problem.

My restored result has some black blocks and somw white blocks.

I guess I stuck on the part of soft-thresholding without normalizing.

Here is my Code :

   #include "StdAfx.h"
   #include "WaveletDenoising.h"
   #include <cmath>


   WaveletDenoising::WaveletDenoising(void)
   {
   }

   WaveletDenoising::~WaveletDenoising(void)
   {
   }

   /* Forward Haar wavelet transform: */
   void WaveletDenoising::ForwardHaar1D(double* data, int length)
   {
    const float inv_sqrt2 = 1/sqrt((double)2.0);

    float norm = 1.0f/sqrt((double)length);

    for(int i=0; i < length; i++) {
        data[i] *= norm;
    }

    float *tmp = new float[length];

    while(length > 1) {
        length /= 2;

        for(int i=0; i < length; i++) {
            tmp[i] = (data[2*i] + data[2*i+1]) * inv_sqrt2;
            tmp[length + i] = (data[2*i] - data[2*i+1]) * inv_sqrt2;
        }

        memcpy(data, tmp, length*2*sizeof(float));
    }

    delete [] tmp;
   }

   /* Transpose matrix: */   
   void WaveletDenoising::Transpose(double *data, int width, int height)
   {
    double *B = new double[width*height];

    for(int y=0; y < height; y++) {
        for(int x=0; x < width; x++) {
            B[x*height + y] = data[y*width + x];
        }
    }

    memcpy(data, B, sizeof(double)*width*height);

    delete [] B;
   }

   /* Forward 2d Haar wavelet transform: */
   void WaveletDenoising::ForwardHaar2D(double* data, int width, int height)
   {
    for(int i=0; i < height; i++) 
        ForwardHaar1D(&data[i*width], width);

    Transpose(data, width, height);

    for(int i=0; i < width; i++) 
        ForwardHaar1D(&data[i*height], height);

    Transpose(data, height, width);
   }

   /* Inverse 1d Haar transform */
   void WaveletDenoising::InverseHaar1D(double* data, int length)
   {
    const float inv_sqrt2 = 1/sqrt((double)2.0);
    float inv_norm = sqrt((double)length);

    float *tmp = new float[length];
    int k = 1;

    while(k < length)  {
        for(int i=0; i < k; i++) {
            tmp[2*i] = (data[i] + data[k+i]) * inv_sqrt2;
            tmp[2*i+1] = (data[i] - data[k+i]) * inv_sqrt2;
        }

        memcpy(data, tmp, sizeof(double)*(k*2));

        k *= 2;
    }

    for(int i=0; i < length; i++) {
        data[i] *= inv_norm;
    }

    delete [] tmp;
   }

   /* Inverse 2d Haar wavelet transform */
   void WaveletDenoising::InverseHaar2D(double* data, int width, int height)
   {
    for(int i=0; i < width; i++) {
        InverseHaar1D(&data[i*height], height);
    }

    Transpose(data, height, width);

    for(int i=0; i < height; i++) {
        InverseHaar1D(&data[i*width], width);
    }

    Transpose(data, width, height);
   }

   /* Image denoising by soft-thresholding */
   void WaveletDenoising::WaveletThresholdDenoising(int width, int height, double* src,    double* des, double threshold)
   {
    int i, j, x, y;

    this->ForwardHaar2D(src, width, height);

    double mi = src[0*width+0]; /* find min value */
    double ma = src[0*width+0]; /* find max value */

    for (y=0; y<height; y++)
    {
        for (x=0; x<width; x++)
        {
            if (mi > src[y*width+x])
                mi = src[y*width+x];
            if (ma < src[y*width+x])
                ma = src[y*width+x];
        }
    }

    /* soft-thresholding */
    for (y=0; y<height; y++)
    {
        for (x=0; x<width; x++)
        {
            if (src[y*width+x] < threshold) 
                src[y*width+x] = 0;
            else if (src[y*width+x] > threshold)
                src[y*width+x] = src[y*width+x] - threshold;
            else if (src[y*width+x] < -threshold)
                src[y*width+x] = src[y*width+x] + threshold;


        }
    }

    this->InverseHaar2D(src, width, height);

    for (y=0; y<height; y++)
    {
        for (x=0; x<width; x++)
        {
            // for normalized:
            src[y*width+x] = ((src[y*width+x] - mi) / (ma - mi))*255;

            double temp =  src[y*width+x];
            if (temp < 0) temp = 0;
            else if (temp >255) temp = 255;
            else temp = temp;
            des[y*width+x] = (BYTE) src[y*width+x];
        }
    }
   }

How can i do this ?

1

1 Answers

1
votes

After spending some hours on this code, I finally found the problem of my code. First, I had to change double type instead of float of the temp variable in InverseHaar1D function. Second, adjust the threshold value in the calling function depending on the degree of noise level. Third, get rid of some redundancy lines as the following clear function. You can see the result with this link: https://www.mediafire.com/?v80rslisl7fff6n

   /* Image denoising by soft-thresholding */
   void WaveletDenoising::WaveletThresholdDenoising(int width, int height, double* src, double* des, double threshold)
   {
    int x, y;

    /* Forward 2d Haar transform */
    this->ForwardHaar2D(src, width, height);

    /* soft-thresholding */
    for(y=0; y < height; y++) 
    {
        for(x=0; x < width; x++) 
        {
            if (src[y*width+x] > threshold)
                src[y*width+x] = src[y*width+x] - threshold;
            else if (src[y*width+x] < -threshold)
                src[y*width+x] = src[y*width+x] + threshold;
            else
                src[y*width+x] = 0;
        }
    }

    /* Inverse 2D Haar transform */
    this->InverseHaar2D(src, width, height);

    for (y=0; y<height; y++)
    {
        for (x=0; x<width; x++)
        {
            double temp =  src[y*width+x];
            if (temp < 0) temp = 0;
            else if (temp >255) temp = 255;
            else 
               des[y*width+x] = (BYTE) temp;
        }
    }
   }