Обесцвечивание изображений с помощью 2d вейвлет-преобразования Хаара с помощью мягкого порога (C ++)

Я пытаюсь реализовать одно из базовых 2D вейвлет-преобразований с помощью преобразования Хаара.

Я применил это к проблеме шумоподавления изображения.

Мой восстановленный результат имеет несколько черных блоков и некоторые белые блоки.

Я думаю, что застрял со стороны мягкого порога без нормализации.

Вот мой код:

   #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];
}
}
}

Как я могу это сделать ?