The goal of this session is to implement several 2D deconvolution algorithms. The implementation will be done entirely in the Fourier domain. To this end, we provide the class, FourierSpace which contains useful methods.
|Public Members||Description of the FourierSpace class||Example|
|double real||double array containing the real coefficients of the Fourier transformation||G.real[fx][fy] = 10.0;|
|double imag||double array containing the imaginary coefficients of the Fourier transformation||G.imag[fx][fy] = 10.0;|
|int nx, ny||Size in X and Y of the Fourier space||int sizex = G.nx;|
|FourierSpace(nx, ny)||Creates an empty (zeros in real and imaginary) Fourier space of size [nx,ny] corresponding to [0..2π,0..2π]||FourierSpace G = new FourierSpace(nx, ny);|
|FourierSpace(image)||Creates a Fourier space containing the Fourier transformation of a ImageAcces object||FourierSpace G = new FourierSpace(image);|
|FourierSpace(double function)||Creates a Fourier space and fills the whole [0..2π][0..2π] Fourier space based on a function defined in [0..π][0..π]|
|void transform(image)||Fills the Fourier space with the transform of an ImageAccess object||G.transform(image);|
|ImageAccess inverse()||Returns the inverse Fourier transform.||ImageAccess g = G.inverse();|
|void showModulus(String title)||Shows the modulus||G.showModule("G");|
1. Coding the Fourier operations
1.1 Implementing some filters in the frequency domain
Write the following filters in the frequency domain using the class FourierSpace.
|Java routine name in the class FourierFilter||Freq. response (ω is the radial frequency: ω2 = ωx2 + ωy2 )|
||These 2 filters are only used as operator in the Regularized Inverse Algorithm. The corresponding codes are provided.|
|gaussian(nx, ny, sigma)||
σ is the standard deviation of the Gaussian
|airy3(nx, ny, mu)
Simulate a Airy pattern by a sum of three sigmoid functions U.
μ acts as a cutoff frequency in radian. μ1 = 0.6*μ and μ2 = 0.4*μ. K=nx/2.
|blurMotion(nx, ny, L)
Model that simulates a horizontal motion blur.
|This filter is easier to implement in the space domain by drawing a line of intensity 1 of the center of an image.
Then, one can compute the Fourier transform using the following command
return new FourierSpace(Display_Filter.shift(f));
The length of the line is given by the paramenter L.
Important note: it is sufficient to compute the symmetric filter part of the function F(ω) = F(ωx,ωy) between [0..Π][0..Π]. The Java construction FourierSpace(function) will fill the Fourier space based on the hermitian symmetry constraint. Thus, these routines should create a 2D array function[nx/2+1][ny/2+1]: function = F(0,0), ... function[nx/2][ny/2] = F(Π,Π)
Observe the impulse response and the frequency response of your filters using the plugin Display Filter. Fill the report.
1.2. Complex arithmetic operations
Get inspiration from the routine add(FourierSpace f1, FourierSpace f2) to write the code of the following complex operations:
Understand the code of the convolve() method that multiplies the Fourier transform of an image and a filter H defined in the Fourier domain. The convolve() method can be called using the plugin Convolve.
Check your code by convolving images (impulse, retina) for the five filters with various settings. Convolve the retina image with Identity, Laplacian, Gaussian (σ=9), Airy (μ=0.5), and BlurMotion (L=30). Fill the report. Keep the five convolved images for the next question.
1.4. Deconvolution by inverse filtering
Write the method deconvolve() that divides the Fourier transform of an image by a specified filter. The deconvolve() method can be called using the plugin Deconvolve, that uses the "Inverse Filtering" algorithm.
Check your code by deconvolving the five images generated in question 1.3. Measure the SNR with respect to the original image (retina) using the plugin SNR. Explain why the SNR of the Laplacian reconstruction is lower than the SNR of the Gaussian reconstruction. Fill the report.
2. Deconvolution algorithms
For this question, we assume that an image is blurred with a Gaussian filter Hσ(ω), σ=4.0. The goal is to evaluate the performance of three deconvolution algorithms on noiseless and noisy images when the model is exactly known "Known Model" H4 (Gaussian, σ=4) and in more realistic conditions, when the model is not perfectly known "Perturbed Model" H4.5 (Gaussian σ=4.5). The input image is an ImageAccess object g(x) and G(ω) is the Fourier transform of g.
From the image moon, simulate the physical convolution process by blurring the image with the plugin Convolve, Gaussian (σ=4.0) and create 3 degraded test-images:
2.1 Inverse filtering restoration
Using the plugin Deconvolve with "Inverse Filtering", apply the inverse filter of the question 1.3.The routine deconvolve() returns f(x) which is the inverse Fourier transform of F(ω):
Deconvolve the three test images with a "Known Model" and the "Pertubated Model" and measure the SNR using the plugin SNR. Compare the results and fill the report.
2.2 Regularized least-square restoration
Implement the inverse filter deconvolveRegularized().
This routine has to return f(x) which is the inverse Fourier transform of F(ω):
L(ω) is the Laplacian filter.
Deconvolve the three test images with a "Known Model" and the "Pertubated Model" (λ=0.5) and measure the SNR. Compare the results and fill the report.
2.3 Landweber algorithm
Implement the Landweber iterative filter deconvolveLandweber(g, H, γ, n) that performs an iterative least-square restoration (n is the total number of iterations).
This routine has to return f(x), the inverse Fourier transform of F(k)(ω).with
Note: If you insert a line of code and check "Display iterations" the interface will display all intermediate images in a stack of images. Be careful, it can slow down the processing and requiring a lot of memory.
Deconvolve_.stepDisplay(F.inverse()); // step is a boolean that indicates if the checkbox is true.
Deconvolve the three test images with a "Known Model" and the "Perturbated Model" (25 iterations and γ=0.5) and measure the SNR. Compare the results and fill in the report.
3.1 Settings of the Landweber algorithm
Iterative algorithms (as Landweber) tend to give better results when they are stopped after a suitable number of iterations. The main difficulty is to find the parameters of the loop: the relaxation parameter γ and the number of iteration n.
Measure the SNR of each iterations with respect to the original image using the plugin SNR. Try to find experimentally n giving the highest SNR for the deconvolved image2 by observing the peak of the SNR both for the "Known Model" H4 and for the "Perturbed Model" H4.5. Here, We choose γ=1. Fill the report.
3.2 Application: deconvolution of microscope images
A biologist has acquired a blurred image. We assume that the microscope is approximated by a Pinhole model where μ is is unknown. Try to restore the blurred image golgi to improve the visual quality of the image. Choose the best way to deconvolve this image with the tools that you developped in this lab. Try to be the closest to the reference image (ground truth). The results are evaluated visually. Report your settings and insert the resulting image in report.
|Reference image||Blurred image|
3.3 Application: restore a motion blur
Try to restore the blurred image image-motion-blur and to read the text of the image. Choose the best way to deconvolve this image with the tools that you developped in this lab. The results are evaluated visually. Report your settings and insert the resulting image in report.