 Biomedical Imaging Group
IP-LAB: Image Processing Laboratories
BIG > Teaching > IP-LAB > Wavelets and Applications 20

# Wavelets and Applications We propose to implement Daubechies wavelet order 0, 1, and 2 using a tree-structured scheme given, here in 1D. The order 0 correspond to the identity, the order 1 to the Haar wavelet, and the order 2 is usually called the DB2 wavelet. The proposed code is a textbook implementation where each module is written in a separated method and tested individually. Here, we assume that the image is square with a size in power that is a power of 2.

1. Convolution

The core of the wavelet transform is the 2D convolution by the 2D filters Ha, Ga, Hs, and Gs that we want to implement in a separable manner. Here, the components of the 2D filters are 1D filters defined by a short mask, an odd size, centered on the middle of the mask, and stored as a 1D (double) array.

Write the method convolve(ImageAccess image, double wx[], double wy[]) that returns the convolved image with the mask wx for the rows and with the mask wy for the columns. We consider periodic boundary conditions.

The method convolve is called in testConvolution with predefined masks for the filter. Apply it to the images impulse and lena using the plugin Test Convolution. Fill the report.

2. Implementation of the Filters

The scaling Ha and the wavelet Ga filters for the analysis part are defined by their odd-size masks as follows:

01 /L /L L=1Provided as example
13 [1  1  0]/L[-1  1  0]/L L=2½To be coded
25 [1-t  3-t  3+t  1+t  0]/L[-1-t  3+t  t-3  1-t  0]/L L=4 * 2½To be coded (t=3½)

a. Write the methods getScalingFunction(int order) and getWaveletFunction(int order) that return a 1D array containing the coefficients of the filters.

b. The synthesis filters Hs and Gs are obtained by reflecting the analysis filters Ha and Ga. Write the method reflect() that returns a reflected version of an 1D array.

c. Check the effect of the 2D filter Ha by applying the scaling 1D mask to rows, then to columns. Write the method testScalingFilter that shows the filtered image and apply this routine to the images impulse and lena. What type of filter is it? Fill the report.

d. Check the effect of the 2D filter Ga by applying the wavelet 1D mask to rows, then to columns. Write the method testWaveletFilter that shows the filtered image and apply this routine to the images impulse and lena. What type of filter is it? Fill the report.

3. Downsampling and Upsampling

1. Write the method downsample that returns an image downsampled by a factor of 2. Apply this routine to lena using the plugin Test Downsampling.
2. Write the method upsample that returns an image upsampled by a factor of 2. Apply this routine to lena using the plugin Test Upsampling.
3. Apply two times an upsampling, then two times a downsampling by 2 on the image lena. Fill the report.
4. Apply two times a downsampling, then two times an upsampling by 2 on the image lena. Fill the report.

4. Analysis The provided method analyze() performs the n-scale wavelet transform of an image; it calls n times the method analysis1().

Write the method analyze1() that performs a single iteration of the wavelet transform of an image for a given order of the Daubechies wavelet. Take advantage of your methods convolve() and downsample() to perform this task. Subimages can be handled using the ImageAccess methods getSubImage() and putSubImage(). The wavelet coefficients are stored in an ImageAccess output (of the same size as the input) in four quadrants:

• LL is obtained by convolving rows by Ha and columns by Ha, followed by downsampling
• HL is obtained by convolving rows by Ga and columns by Ha, followed by downsampling
• LH is obtained by convolving rows by Ha and columns by Ga, followed by downsampling
• HH is obtained by convolving rows by Ga and columns by Ga, followed by downsampling

Apply the analysis part of the wavelet transform to lena using the plugin Wavelet Daubechies, order 2, by choosing 1 scale and clicking on the button "Analysis". Fill the report.

5. Synthesis The synthesis part reconstructs an image that was decomposed by the analysis part. Your assignment is to code the two methods that perform the inverse wavelet transform.

• The method synthesize() should perform n iterations of the inverse wavelet transform of an image.
• The method synthesize1() should perform a single iteration of the inverse wavelet transform of an image.

Apply the synthesis part of the wavelet transform to mit-coef using the plugin Wavelet Daubechies, order 2, by choosing 1 scale and clicking on the button "Synthesis". Fill the report.

6. Processing of the wavelet coefficients  The processing of the wavelet coefficient can either be checked as a separate module by clicking on the button "Processing", or be incorporated in the whole chain (analysis + processing + synthesis) by clicking on the button "Run all".

1. Write the method identity() that returns an exact copy of the coefficients.
2. Write the method keepLowpass() which only keeps the low-pass coefficients (LLn) and set to zero the other coefficents.
3. Write the method doSoftThreshold(coef, T) which returns coefficients t(x) = sgn(w(x)) max(0, |w(x)|-T). Note that the w(x) are the input wavelet coefficients and T is a parameter given by the user through the dialog box.
4. For each processing, test the whole chain on the doisneau image for n=5 with the Daubechies wavelet, order 2. Fill the report.

7. Application: Automatic Denoising

The denoising is achieved by applying a soft threshold to the coefficients of the wavelet transform. A rule of thumbs is to choose T = σHH1 as threshold for the denoising. σHH1 is the standard deviation of the noise which can be estimated in the highpass band (HH1).

Write the method denoiseAuto() that performs, the analysis, the soft threshold, and the synthesis using the db2 wavelet with 3 scales. The threshold T is automatically determined, based on the value σHH1 that you have computed in this method (ImageAccess has a method to compute the standard deviation). Print T using IJ.log();.

Apply this routine to doisneau-noise, doisneau, cells, and lowlight using the plugin Automatic Denoise. Fill the report.

Code.java

```import ij.*;

public class Code {

private ImageAccess convolve(ImageAccess in, double wx[], double wy[]) {
// TODO add you code here
return in;
}

public void testConvolution(ImageAccess in) {
double[] wx = {-2, -1, 0, 1, 2};
double[] wy = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
ImageAccess s1 = convolve(in, wx, wy);
s1.show("Test Convolution of " + in.getTitle());
}

public double[] getScalingFunction(int order) {
if (order == 0) {
return new double[] {1};
}
// TODO complete for order 1 and 2
else
throw new ArithmeticException("Unvalid order");
}

public double[] getWaveletFunction(int order) {
if (order == 0) {
return new double[] {1};
}
// TODO complete for order 1 and 2
else
throw new ArithmeticException("Unvalid order");
}

private double[] reflect(double arr[]) {
int n = arr.length;
double[] out = new double[n];
// TODO add you code here
return out;
}

public void testScalingFilter(ImageAccess in, int order) {
// TODO add you code here
in.show("Scaling db" + order);
}

public void testWaveletFilter(ImageAccess in, int order) {
// TODO add you code here
in.show("Scaling db" + order);
}

public ImageAccess analyze(ImageAccess image, int n, int order) {
int nx = image.getWidth();
int ny = image.getHeight();
ImageAccess output = image.duplicate();
for (int i = 0; i < n; i++) {
ImageAccess sub = output.getSubImage(0, 0, nx, ny);
sub = analyze1(sub, order);
output.putSubImage(0, 0, sub);
nx = nx / 2;
ny = ny / 2;
}
return output;
}

private ImageAccess analyze1(ImageAccess in, int order) {
// TODO add you code here
return in;
}

public ImageAccess synthesize(ImageAccess in, int n, int order) {
// TODO add you code here
return in;
}

private ImageAccess synthesize1(ImageAccess coef, int order) {
// TODO add you code here
return coef;
}

public ImageAccess downsample(ImageAccess in) {
// TODO add you code here
return in;
}

public ImageAccess upsample(ImageAccess in) {
// TODO add you code here
return in;
}
public ImageAccess identity(ImageAccess in, int n) {
// TODO add you code here
return in;
}

public ImageAccess keepLowpass(ImageAccess in, int n) {
// TODO add you code here
return in;
}

public ImageAccess doSoftThreshold(ImageAccess in, double threshold) {
// TODO add you code here
return in;
}

public ImageAccess denoiseAuto(ImageAccess in) {
// TODO add you code here
// IJ.log("Threshold " + threshold);
return in;
}
}

```

webmaster.big@epfl.ch • 09.04.2020