deconvlucy (Image Processing Toolbox User's Guide)

Image Processing Toolbox User's Guide

Restore image using the Lucy-Richardson algorithm

Syntax

J = deconvlucy(I,PSF)
J = deconvlucy(I,PSF,NUMIT)
J = deconvlucy(I,PSF,NUMIT,DAMPAR)
J = deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT)
J = deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT)
J = deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT,SUBSMPL)

Description

J = deconvlucy(I,PSF) restores image I that was degraded by convolution with a point-spread function PSF and possibly by additive noise. The algorithm is based on maximizing the likelihood of the resulting image J's being an instance of the original image I under Poisson statistics. The input array I can be a numeric array (the blurred image) or a cell array.

If I is a cell array, it can contain a single numerical array (the blurred image) or it can be the output from a previous run of deconvlucy.

When you pass a cell array to deconvlucy as input, it returns a 1-by-4 cell array J, where

J{1} contains I, the original image.

J{2} contains the result of the last iteration.

J{3} contains the result of the next-to-last iteration.

J{4} is an array generated by the iterative algorithm.

To improve the restoration, deconvlucy supports several optional parameters. Use [] as a placeholder if you do not specify an intermediate parameter.

J = deconvlucy(I,PSF,NUMIT) specifies the number of iterations the deconvlucy function performs. If this value is not specified, the default is 10.

J = deconvlucy(I,PSF,NUMIT,DAMPAR) specifies the threshold deviation of the resulting image from the image I (in terms of the standard deviation of Poisson noise) below which damping occurs. Iterations are suppressed for pixels that deviate beyond the DAMPAR value from their original value. This suppresses the noise generation in such pixels, preserving necessary image details elsewhere. The default value is 0 (no damping).

J = deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT) specifies the weight to be assigned to each pixel to reflect its recording quality in the camera. A bad pixel is excluded from the solution by assigning it zero weight value. Instead of giving a weight of unity for good pixels, you can adjust their weight according to the amount of flat-field correction. The default is a unit array of the same size as input image I.

J = deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT) specifies a value corresponding to the additive noise (e.g., background, foreground noise) and the variance of the readout camera noise. READOUT has to be in the units of the image. The default value is 0.

J = deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT,SUBSMPL), where SUBSMPL denotes subsampling and is used when the PSF is given on a grid that is SUBSMPL times finer than the image. The default value is 1.

Note The output image J could exhibit ringing introduced by the discrete Fourier transform used in the algorithm. To reduce the ringing, use I = edgetaper(I,PSF) before calling deconvlucy.

Class Support

I and PSF can be uint8, uint16, int16, double, or single. DAMPAR and READOUT must have the same class as the input image. Other inputs have to be double. The output image J (or the first array of the output cell) has the same class as the input image I.

Example

I = checkerboard(8);
PSF = fspecial('gaussian',7,10);
V = .0001;
BlurredNoisy = imnoise(imfilter(I,PSF),'gaussian',0,V);
WT = zeros(size(I));
WT(5:end-4,5:end-4) = 1;
J1 = deconvlucy(BlurredNoisy,PSF);
J2 = deconvlucy(BlurredNoisy,PSF,20,sqrt(V));
J3 = deconvlucy(BlurredNoisy,PSF,20,sqrt(V),WT);

subplot(221);imshow(BlurredNoisy);
title('A = Blurred and Noisy');
subplot(222);imshow(J1);
title('deconvlucy(A,PSF)');
subplot(223);imshow(J2);
title('deconvlucy(A,PSF,NI,DP)');
subplot(224);imshow(J3);
title('deconvlucy(A,PSF,NI,DP,WT)');

See Also

deconvblind, deconvreg, deconvwnr, otf2psf, padarray, psf2otf

deconvblind deconvreg