Как работает алгоритм Ричардсона – Люси?Пример кода? - PullRequest
Новая
       11

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

5 голосов
/ 24 марта 2012

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

До сих пор я нашел алгоритм Ричардсона – Люси , где математика не кажется такой сложной, однако я не могу понять, как работает настоящий алгоритм. В Википедии написано:

Это приводит к уравнению, для которого можно решить итеративно в соответствии ...

однако он не показывает фактический цикл. Может кто-нибудь указать мне на ресурс, где объясняется фактический алгоритм. В Google мне только удается найти методы, которые используют Ричардсона-Люси в качестве одного из своих шагов, но не фактический алгоритм Ричардсона-Люси.

Алгоритм на любом языке или псевдокоде был бы хорош, однако, если бы он был доступен на Python, это было бы удивительно.

Спасибо заранее.

Редактировать

По сути, я хочу получить размытое изображение (nxm): 

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

и ядро ​​(ixj), которое использовалось для получения размытого изображения: 

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

Каковы точные шаги в алгоритме Ричардсона – Люси, чтобы выяснить исходное изображение.

Ответы [ 4 ]

4 голосов
/ 08 февраля 2016

Вот очень простая реализация 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')

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

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

1 голос
/ 19 июля 2012

Уравнение в Википедии дает функцию для итерации 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 голос
/ 29 июля 2013

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

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

Ричардсон Люси - это строительный блок для многих других алгоритмов деконволюции.Например, приведенный выше пример iocbio изменил алгоритм, чтобы лучше справляться с шумом.

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

1 голос
/ 24 марта 2012

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

...