Как работает алгоритм Ричардсона-Люси? Пример кода?

Вот очень простая реализация Matlab:

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

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

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

Вы также можете работать в частотной области, а не в пространственной области, как указано выше. В этом случае код будет:

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 RL для пространственно-вариантных PSF (т.е. пространственно-неоднородных функций разброса точек) - если у кого-нибудь будет такая, пожалуйста, дайте мне знать!

Чтобы избавиться от артефактов по краям, вы можете отразить входное изображение по краям, а затем обрезать зеркальные биты или использовать image = edgetaper(image, PSF) Matlab перед вызовом RL_deconv.

Нативная реализация deconvlucy.m для Matlab немного сложнее, кстати — исходный код этой реализации можно найти здесь и использует ускоренная версия базового алгоритма.

Уравнение в Википедии дает функцию для итерации 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, корректную для любого типа данных, с которыми вы работаете. Вам также необходимо обработать случай, когда сходимость никогда не достигается (например, ограничить количество итераций).

Вот реализация Python с открытым исходным кодом: http://code.google.com/p/iocbio/wiki/IOCBioMicroscope

Если это поможет, вот реализация, которую я написал, которая включает в себя некоторую документацию....

https://github.com/bnorthan/projects/blob/master/truenorthJ/ImageJ2Plugins/functions/src/main/java/com/truenorth/functions/fft/filters/RichardsonLucyFilter.java

Richardson Lucy является строительным блоком для многих других алгоритмов деконволюции. Например, приведенный выше пример iocbio модифицировал алгоритм, чтобы лучше справляться с шумом.

Это относительно простой алгоритм (как это бывает) и является отправной точкой для более сложных алгоритмов, поэтому вы можете найти множество различных реализаций.