Jak działa algorytm Richardsona-Lucy? Przykład kodu?

Question

Próbuję dowiedzieć się, jak działa dekonwolucja. Rozumiem ideę, która się za tym kryje, ale chcę zrozumieć niektóre z rzeczywistych algorytmów, które ją implementują - algorytmy, które biorą jako dane wejściowe zamazany obraz z jego funkcją próbki punktowej (rozmycie jądra) i wytwarzają jako wynik ukryty obraz.

Do tej pory znalazłem algorytm Richardson–Lucy, w którym matematyka nie wydaje się być tak trudna, ale nie wiem, jak działa dany algorytm. W Wikipedia mówi:

Prowadzi to do równania, dla których może być rozwiązany iteracyjnie według ...

jednak nie pokazuje faktyczną pętlę. Czy każdy może wskazać mi zasób, w którym wyjaśniono rzeczywisty algorytm. W Google I udało mi się znaleźć tylko metody, które wykorzystują Richardson-Lucy jako jeden z jego kroków, ale nie rzeczywisty algorytm Richardsona-Lucy.

Algorytm w dowolnym języku lub pseudo-kodzie byłby fajny, jednak jeśli jest dostępny w Pythonie, byłoby to niesamowite.

Thanx z góry.

Edit

Zasadniczo co chcę dowiedzieć podano zamazany obraz (NXM):

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

i jądra (ixj), który został użyty w celu uzyskania zamazany obraz:

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

Jakie są dokładne kroki algorytmu Richardsona-Lucy w celu ustalenia oryginalnego obrazu.

Answer 1

Oto open source implementacja Pythona: http://code.google.com/p/iocbio/wiki/IOCBioMicroscope

Answer 2

Równanie Wikipedia daje funkcję iteracji t+1 pod względem iteracji t. można wdrożyć ten typ algorytmu iteracyjnego w następujący sposób:

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 

Oczywiście, trzeba będzie zaimplementować własną difference funkcję, która jest prawidłowa dla niezależnie od typu danych, które pracują z. Musisz również zająć się przypadkiem, w którym nigdy nie osiągnięto zbieżności (na przykład ograniczyć liczbę iteracji).

Answer 3

Jeśli to pomaga tutaj jest realizacja pisałem, że zawiera jakąś dokumentację ....

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

Richardson Lucy jest budulcem dla wielu innych algorytmów dekonwolucji. Na przykład powyższy przykład iocbio zmodyfikował algorytm, aby lepiej radzić sobie z hałasem.

Jest to względnie prosty algorytm (jak to się dzieje) i jest punktem wyjścia dla bardziej skomplikowanych algorytmów, dzięki czemu można znaleźć wiele różnych implementacji.

Answer 4

Tutaj jest bardzo prosta implementacja 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); toc; 
figure; imshow(res_RL2); title('Recovered Image') 

Można także pracować w dziedzinie częstotliwości zamiast w dziedzinie przestrzeni, jak powyżej. W takim przypadku kod byłby następujący:

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

Jedyne, czego nie bardzo rozumiem, to jak działa to przestrzenne odwrócenie PSF i do czego służy. Jeśli ktokolwiek mógłby mi to wyjaśnić, byłoby fajnie! Poszukuję również prostej implementacji Matlab R-L dla przestrzennie wariantowych PSF (tj. Przestrzennie niehomogenicznych funkcji rozprzestrzeniania się punktów) - jeśli ktokolwiek miałby taką, proszę dać mi znać!

Aby pozbyć się artefaktów na krawędziach, można odbić obraz wejściowy na krawędziach, a następnie odciąć odbitki lustrzane lub użyć Matlaba image = edgetaper(image, PSF), zanim zadzwonisz pod numer RL_deconv.

Natywna implementacja programu Matlab deconvlucy.m jest nieco bardziej skomplikowana - kod źródłowy tego kodu można znaleźć pod numerem here i używa on accelerated version of the basic algorithm.