Wavelets are
mathematical functions that cut up data into different frequency components,
and then study each component with a resolution matched to its scale. They have
advantages over traditional Fourier methods in analyzing physical situations
where the image contains discontinuities and sharp spikes.
Wavelets works in the following way. When you decompose a data set using wavelets,
you use filters that act as averaging filters and others that
produce details. Some of the resulting wavelet coefficients
correspond to details in the data set. If the details are small, they might be
omitted without substantially affecting the main features of the data set. The
idea of thresholding, then, is to set to zero all coefficients
that are less than a particular threshold. These coefficients are used in an
inverse wavelet transformation to reconstruct the data set.
Sample code:
Denoise :
void wavelet_denoise(float **noisy,float
**denoised,int Ni,int
Nj,float threshold)
{
int i,j,k;
int oc_max,oc_r,oc_c;
float **buffer_t,**buffer;
int **normalization;
int
shift_arr_r[MAX_ARR_SIZE],shift_arr_c[MAX_ARR_SIZE];
// 4 level wavelet transform.
const int levs=4;
int Nl,Nh;
float *lp,*hp;
// Choose wavelet filter. See wav_filters.h and wav_trf.c.
choose_filter('D',9);
// Main buffer for operations.
buffer_t=allocate_2d_float(Ni,Nj,0);
// Buffer on which results are accumulated.
buffer=allocate_2d_float(Ni,Nj,1);
// To keep track of the number of overcomplete results
// on each pixel. Legacy stuff, this is really not
necessary
// when doing things image-wide.
normalization=allocate_2d_int(Ni,Nj,1);
// Variable that controls the number of overcomplete
computations.
// oc_max=1<
// oc_max=2 overcomplete at the first level, then regular
(4 x redundant)...
// see comments below.
oc_max=2;
for(oc_r=0;oc_r
// At every level of a 1-D wavelet transform there are two transforms possible
// obtained by the original bank or its one pixel shifted
version.
// In L levels there are 2^L possibilities with each of
these possibilities
// forming a complete transform (i.e., invertible, etc.).
If we evaluate
// more than one then what we end up with is an
overcomplete transform.
// In 2-D there are 4 possibilities at each level and hence
4^L total
// possibilities.
//
// At level l, this code selects either the
"regular" (shift_arr_r[l]=0)
// or the one shifted bank (shift_arr_r[l]=1) over rows and
likewise
// over columns using shift_arr_c[l].
//
// In order to generate the fully overcomplete transform we
effectively
// put the binary representation of numbers from 0 to 2^L-1
(represented
// by oc_max) in shift_arr_r and shift_arr_c to traverse
all 4^L possibilities.
//
// In fully overcomplete denoising example, we evaluate
each of the 4^L transforms,
// hard threshold coefficients, inverse transform, and
finally average the results.
// Convert oc_r into binary and put each digit (0-1) into
shift_arr_r.
// Each digit causes the filters at the resulting wavelet
level
// to be shifted by one (1) or remain unshifted (0).
// For filter shifts over rows.
k=oc_r;
for(i=levs-1;i>=0;i--) {
shift_arr_r[i]=k>>i;
k-=shift_arr_r[i]<
}
for(oc_c=0;oc_c
// For filter shifts over columns.
k=oc_c;
for(i=levs-1;i>=0;i--) {
shift_arr_c[i]=k>>i;
k-=shift_arr_c[i]<
}
// Copy the "noisy" region into the buffer.
for(i=0;i
for(j=0;j
buffer_t[i][j]=noisy[i][j];
}
}
// Select the forward bank of filters.
lp=MFLP;Nl=Nflp;
hp=MFHP;Nh=Nfhp;
if(PACKET) {
// Packet transform.
wavpack2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c);
// Threshold.
hard_threshold(buffer_t,Ni,Nj,threshold);
}
else {
// Regular wavelets.
wav2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,1,shift_arr_r,shift_arr_c);
// Threshold.
hard_threshold(buffer_t,Ni,Nj,threshold);
}
// Select the inverse bank of filters.
lp=MILP;Nl=Nilp;
hp=MIHP;Nh=Nihp;
// Inverse transform.
if(PACKET)
wavpack2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,0,shift_arr_r,shift_arr_c);
else
wav2d_inpl(buffer_t,Ni,Nj,levs,lp,Nl,hp,Nh,0,shift_arr_r,shift_arr_c);
// Accumulate results.
for(i=0;i
for(j=0;j
buffer[i][j]+=buffer_t[i][j];
normalization[i][j]++;
}
}
}
for(i=0;i
for(j=0;j
// Copy final results into place.
if(normalization[i][j])
denoised[i][j]=buffer[i][j]/normalization[i][j];
else
denoised[i][j]=noisy[i][j];
}
// Free buffers.
free_2d_float(buffer_t,Ni);
free_2d_float(buffer,Ni);
free_2d_int(normalization,Ni);
}
Add Noise
/***********************************************************/
void add_uniform_noise(float **orig,int Ni,int Nj,float **noisy,float mean,float
stdev)
{
int i,j;
static int kilroy=0;
// The way uniform noise is generated has variance 1/12.
// Multiply with sqrt(12)*stdev to generate noise with
variance=stdev^2.
float scaler=(float)(sqrt(12)*stdev);
if(!kilroy) {
kilroy=1;
srand(0);
}
for(i=0;i
for(j=0;j
noisy[i][j]=orig[i][j]+(float)((rand()/(RAND_MAX+1.0)-.5)*scaler+mean);
}
