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 H_{a}, G_{a}, H_{s}, and G_{s} 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 w_{x} for the rows and with the mask w_{y} 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 H_{a} and the wavelet G_{a} filters for the analysis part are defined by their odd-size masks as follows:
Order | Taps | Scaling mask (H_{a}) | Wavelet mask (G_{a}) | Normalisation factor | Java code |
---|---|---|---|---|---|
0 | 1 | [1]/L | [1]/L | L=1 | Provided as example |
1 | 3 | [1 1 0]/L | [-1 1 0]/L | L=2^{½} | To be coded |
2 | 5 | [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 H_{s} and G_{s} 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 H_{a} 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 G_{a} 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
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:
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.
Take advantage of your methods convolve() and upsample() to perform this task.
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".
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; } }