Geek Answers Handbook

How does the Richardson-Lucy algorithm work? Code example?

I am trying to understand how deconvolution works. I understand the idea, but I want to understand some of the real algorithms that implement it - algorithms that take a blurred image with its point selection function (blur the core) as input and produce a hidden image as the output.

So far I have found the Richardson-Lucy algorithm , where the math does not seem so difficult, however I cannot understand how the actual algorithm works. Wikipedia says:

This leads to an equation for which it can be solved iteratively according to ...

however, it does not show the actual cycle. Can someone point me to a resource where the actual algorithm is explained. At Google, I manage to find methods that use Richardson-Lucy as one of his steps, but not the real Richardson-Lucy algorithm.

An algorithm in any language or pseudo-code will be nice, however, if it is available in Python, it will be awesome.

Thanx in advance.

Edit

Essentially, what I want to find out is a blurry image (nxm) is given:

x00 x01 x02 x03 .. x0n
x10 x11 x12 x13 .. x1n
...
xm0 xm1 xm2 xm3 .. xmn

and the kernel (ixj), which was used to get the blurry image:

p00 p01 p02 .. p0i
p10 p11 p12 .. p1i
...
pj0 pj1 pj2 .. pji

What are the exact steps of the Richardson-Lucy algorithm to figure out the original image.

+5
source share
4 answers

Here is a very simple Matlab implementation:

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

, , . :

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);

, , - , PSF . - , ! Matlab R-L - PSF ( ) - - , !

, , Matlab image = edgetaper(image, PSF), RL_deconv.

Matlab deconvlucy.m btw - .

+4

t+1 t. :

def iter_step(prev):
  updated_value = <function from Wikipedia>
  return updated_value

def iterate(initial_guess):
  cur = initial_guess
  while True:
    prev, cur = cur, iter_step(cur)
    if difference(prev, cur) <= tolerance:
      break
  return cur

, difference, , . , (, ).

+1

More articles:


All Articles