Deblurring Images Using the Lucy-Richardson Algorithm

Step 1: Read Image

The example reads in an RGB image and crops it to be 256-by-256-by-3. The deconvlucy function can handle arrays of any dimension.

I = imread("board.tif");
I = I(50+(1:256),2+(1:256),:);
imshow(I)
title("Original Image")

Step 2: Simulate a Blur and Noise

Simulate a real-life image that could be blurred due to camera motion or lack of focus. The image could also be noisy due to random disturbances. The example simulates the blur by convolving a Gaussian filter with the true image (using imfilter). The Gaussian filter then represents a point-spread function, PSF.

PSF = fspecial("gaussian",5,5);
Blurred = imfilter(I,PSF,"symmetric","conv");
imshow(Blurred)
title("Blurred")

The example simulates the noise by adding a Gaussian noise of variance V to the blurred image (using imnoise). The noise variance V is used later to define a damping parameter of the algorithm.

V = .002;
BlurredNoisy = imnoise(Blurred,"gaussian",0,V);
imshow(BlurredNoisy)
title("Blurred & Noisy")

Step 3: Restore the Blurred and Noisy Image

Restore the blurred and noisy image providing the PSF and using only 5 iterations (default is 10). The output is an array of the same type as the input image.

luc1 = deconvlucy(BlurredNoisy,PSF,5);
figure;
imshow(luc1);
title("Restored Image, NUMIT = 5");

Step 4: Iterate to Explore the Restoration

The resulting image changes with each iteration. To investigate the evolution of the image restoration, you can do the deconvolution in steps: do a set of iterations, see the result, and then resume the iterations from where they were stopped. To do so, the input image has to be passed as a part of a cell array. For example, start the first set of iterations by passing in {BlurredNoisy} instead of BlurredNoisy as the input image parameter.

luc1_cell = deconvlucy({BlurredNoisy},PSF,5);

In that case the output, luc1_cell, becomes a cell array. The cell output consists of four numeric arrays, where the first is the BlurredNoisy image, the second is the restored image of class double, the third array is the result of the one-before-last iteration, and the fourth array is an internal parameter of the iterated set. The second numeric array of the output cell-array, image luc1_cell{2}, is identical to the output array of the Step 3 image, luc1, with a possible exception of their class (the cell output always gives the restored image of class double).

To resume the iterations, take the output from the previous function call, the cell-array luc1_cell, and pass it into the deconvlucy function. Use the default number of iterations (NUMIT = 10). The restored image is the result of a total of 15 iterations.

luc2_cell = deconvlucy(luc1_cell,PSF);
luc2 = im2uint8(luc2_cell{2});
imshow(luc2)
title("Restored Image, NUMIT = 15")

Step 5: Control Noise Amplification by Damping

The latest image, luc2, is the result of 15 iterations. Although it is sharper than the earlier result from 5 iterations, the image develops a "speckled" appearance. The speckles do not correspond to any real structures (compare it to the true image), but instead are the result of fitting the noise in the data too closely.

To control the noise amplification, use the damping option by specifying the DAMPAR parameter. DAMPAR has to be of the same class as the input image. The algorithm dampens changes in the model in regions where the differences are small compared with the noise. The DAMPAR used here equals 3 standard deviations of the noise. Notice that the image is smoother.

DAMPAR = im2uint8(3*sqrt(V));
luc3 = deconvlucy(BlurredNoisy,PSF,15,DAMPAR);
imshow(luc3)
title("Restored Image with Damping, NUMIT = 15")

The next part of this example explores the weight and subsample input parameters of the deconvlucy function, using a simulated star image (for simplicity & speed).

Step 6: Create Sample Image

The example creates a black/white image of four stars.

I = zeros(32);
I(5,5) = 1;
I(10,3) = 1;
I(27,26) = 1;
I(29,25) = 1;

imshow(1-I,[],InitialMagnification="fit");
ax = gca;
ax.Visible = "on";
ax.XTickLabel = [];
ax.YTickLabel = [];
ax.XTick = [7 24];
ax.XGrid = "on";
ax.YTick = [5 28];
ax.YGrid = "on";
title("Data")

Step 7: Simulate a Blur

The example simulates a blur of the image of the stars by creating a Gaussian filter, PSF, and convolving it with the true image.

PSF = fspecial("gaussian",15,3);
Blurred = imfilter(I,PSF,"conv","sym");

Now simulate a camera that can only observe part of the stars' images (only the blur is seen). Create a weighting function array, WT, that consists of ones in the central part of the Blurred image ("good" pixels, located within the dashed lines) and zeros at the edges ("bad" pixels - those that do not receive the signal).

WT = zeros(32);
WT(6:27,8:23) = 1;
CutImage = Blurred.*WT;

To reduce the ringing associated with borders, apply the edgetaper function with the given PSF.

CutEdged = edgetaper(CutImage,PSF);

imshow(1-CutEdged,[],"InitialMagnification","fit");
ax = gca;
ax.Visible = "on";
ax.XTickLabel = [];
ax.YTickLabel = [];
ax.XTick = [7 24];
ax.XGrid = "on";
ax.YTick = [5 28];
ax.YGrid = "on";
title("Observed");

Step 8: Provide the WEIGHT Array

The algorithm weights each pixel value according to the weight array while restoring the image. In our example, only the values of the central pixels are used (where WT = 1), while the "bad" pixel values are excluded from the optimization. However, the algorithm can place the signal power into the location of these "bad" pixels, beyond the edge of the camera's view. Notice the accuracy of the resolved star positions.

luc4 = deconvlucy(CutEdged,PSF,300,0,WT);

imshow(1-luc4,[],InitialMagnification="fit");
ax = gca;
ax.Visible = "on";
ax.XTickLabel = [];
ax.YTickLabel = [];
ax.XTick = [7 24];
ax.XGrid = "on";
ax.YTick = [5 28];
ax.YGrid = "on";
title("Restored")

Step 9: Provide a finer-sampled PSF

deconvlucy can restore undersampled image given a finer sampled PSF (finer by subsample times). To simulate the poorly resolved image and PSF, the example bins the Blurred image and the original PSF, two pixels in one, in each dimension.

Binned = squeeze(sum(reshape(Blurred,[2 16 2 16])));
BinnedImage = squeeze(sum(Binned,2));
Binned = squeeze(sum(reshape(PSF(1:14,1:14),[2 7 2 7])));
BinnedPSF = squeeze(sum(Binned,2));

imshow(1-BinnedImage,[],InitialMagnification="fit")
ax = gca;
ax.Visible = "on";
ax.XTick = [];
ax.YTick = [];
title("Binned Observed")

Restore the undersampled image, BinnedImage, using the undersampled PSF, BinnedPSF. Notice that the luc5 image distinguishes only 3 stars.

luc5 = deconvlucy(BinnedImage,BinnedPSF,100);

imshow(1-luc5,[],InitialMagnification="fit")
ax = gca;
ax.Visible = "on";
ax.XTick = [];
ax.YTick = [];
title("Poor PSF")

The next example restores the undersampled image (BinnedImage), this time using the finer PSF (defined on a subsample-times finer grid). The reconstructed image (luc6) resolves the position of the stars more accurately. Note how it distributes power between the two stars in the lower right corner of the image. This hints at the existence of two bright objects, instead of one, as in the previous restoration.

luc6 = deconvlucy(BinnedImage,PSF,100,[],[],[],2);

imshow(1-luc6,[],InitialMagnification="fit")
ax = gca;
ax.Visible = "on";
ax.XTick = [];
ax.YTick = [];
title("Fine PSF")