deconvlucy (Image Processing Toolbox)

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,READOUT)
J = deconvlucy(I,PSF,NUMIT,DAMPAR,READOUT,WEIGHT)
deconvlucy(I,PSF,NUMIT,DAMPAR,READOUT,WEIGHT,SUBSMPL)

Description

J = deconvlucy(I,PSF) restores image I, 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 being an instance of the original image I under Poisson statistics. The input array, I, can be a numeric array or cell array. (Use a cell array when you want to be able to perform additional deconvolutions that start where your initial deconvolution finished. See Resuming Deconvolution for more information.)

To improve the restoration, deconvlucy supports several optional parameters. Use [] as a place holder 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, one could 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 read-out 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) prior calling deconvlucy.

Resuming Deconvolution

If input I is a cell array, the output J becomes a cell array of size 1-by-4, 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

The input cell array could contain one numerical array (the blurred image), or four numerical arrays if it was the output from the previous run of deconvlucy.

Class Support

I can be of class uint8, uint16, or double. The DAMPAR and READOUT arguments have to be of the same class as the input image. Other inputs have to be of class double. Output image (or the first array of the output cell) is of the same class as the input image.

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