deconvlucy() - Image Processing
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.I can be a N-dimensional array.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.
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)
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)');
Output / Return Value
Limitations
Alternatives / See Also
Reference