Tutorial :Blind deconvolution algorithm question


I was studying deconvolution,
and stumbled upon Richardson-Lucy deconvolution,
I was thinking of writing a simple program to do post-processing using this method,
does anybody know where I can find complete implementable algorithms or source code that I can study and play around with?

Preferably in C++ language or matlab.

I have read a few books but they are a little general and too theoretical.

thanks, Charles Mawby.
but i'm still having problems looking for the .m files online,
all i get are reference form a reference not a real file.
really appreciate if you can provide more details.
thanks in advance!


MATLAB has a decent implementation (search for corelucy.m and deconvlucy.m on Google, or download the MATLAB and image processing toolbox demos).

MathWorks has documentation on their website:


The outer loop (setting up the deconvolution point spread function, doing iterations) is here:

The inner loop (the core part of the LR algorithm):

Awefully nice of NASA to host parts of MATLAB!


I would suggest using MATLAB, or the F/OSS alternative GNU Octave. They're much better for this sort of thing, as they have libraries of image processing routines, and convolution is a heavily-optimized built-in function.


If you want to use the Richardson-Lucy deconvolution algorithm that comes in the MATLAB Image Processing toolbox (DECONVLUCY), you need to first get the Image Processing toolbox =). In another answer I gave for an SO question, I mentioned how you can get trials of MATLAB and its various toolboxes, if you don't already have them. Once you get the toolbox, you will probably be able to look at the source code for DECONVLUCY (either .m or .c files) to study the algorithm and figure out how it works.


Here is a very simple Matlab implementation of Richardson-Lucy deconvolution :

function result = RL_deconv(image, PSF, iterations)      % to utilise the conv2 function we must make sure the inputs are double      image = double(image);      PSF = double(PSF);      latent_est = image; % initial estimate, or 0.5*ones(size(image));       PSF_HAT = PSF(end:-1:1,end:-1:1); % spatially reversed psf      % iterate towards ML estimate for the latent image      for i= 1:iterations          est_conv      = conv2(latent_est,PSF,'same');          relative_blur = image./est_conv;          error_est     = conv2(relative_blur,PSF_HAT,'same');           latent_est    = latent_est.* error_est;      end      result = latent_est;    original = im2double(imread('lena256.png'));  figure; imshow(original); title('Original Image')  

enter image description here

hsize=[9 9]; sigma=1;  PSF = fspecial('gaussian', hsize, sigma);  blr = imfilter(original, PSF);  figure; imshow(blr); title('Blurred Image')  

enter image description here

res_RL = RL_deconv(blr, PSF, 1000); toc;  figure; imshow(res_RL2); title('Recovered Image')  

enter image description here

You can also work in the frequency domain instead of in the spatial domain as above. In that case the code would be :

function result = RL_deconv(image, PSF, iterations)  fn = image; % at the first iteration  OTF = psf2otf(PSF,size(image));   for i=1:iterations      ffn = fft2(fn);       Hfn = OTF.*ffn;       iHfn = ifft2(Hfn);       ratio = image./iHfn;       iratio = fft2(ratio);       res = OTF .* iratio;       ires = ifft2(res);       fn = ires.*fn;   end  result = abs(fn);   

To get rid of the artefacts at the edges you could mirror the input image at the edges and then crop away the mirrored bits afterwards or use Matlab's image = edgetaper(image, PSF) before you call RL_deconv.

The native Matlab implementation deconvlucy.m is a bit more complicated btw - the source code of that one can be found here and uses an accelerated version of the basic algorithm.

Note:If u also have question or solution just comment us below or mail us on toontricks1994@gmail.com
Next Post »