Graphic STI
logo EPFL
Text EPFL
english only
Biomedical Imaging Group
IP-LAB: Image Processing Laboratories
BIG > Teaching > IP-LAB > Deconvolution 19

Deconvolution

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 MembersDescription of the FourierSpace classExample
double real[][]double array containing the real coefficients of the Fourier transformationG.real[fx][fy] = 10.0;
double imag[][]double array containing the imaginary coefficients of the Fourier transformationG.imag[fx][fy] = 10.0;
int nx, nySize in X and Y of the Fourier spaceint sizex = G.nx;
Constructors
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 objectFourierSpace 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..π]
Public Methods
void transform(image)Fills the Fourier space with the transform of an ImageAccess objectG.transform(image);
ImageAccess inverse()Returns the inverse Fourier transform.ImageAccess g = G.inverse();
void showModulus(String title)Shows the modulusG.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 )
identity(nx, ny)
identity These 2 filters are only used as operator in the Regularized Inverse Algorithm. The corresponding codes are provided.
laplacian(nx, ny) laplacian
gaussian(nx, ny, sigma) gaussian
σ is the standard deviation of the Gaussian
airy3(nx, ny, mu)

Simulate a Airy pattern by a sum of three sigmoid functions U.
pinhole
μ 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(ωxy) 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[0][0] = 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:

1.3. Convolution

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(ω):
inverse

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(ω):
regularized

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)(ω).

regularized   with   regularized

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. Application

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 imageBlurred 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.



webmaster.big@epfl.ch • 16.05.2019