Come funziona l'algoritmo di Richardson-Lucy? Esempio di codice?

Question

Sto cercando di capire come funziona la deconvoluzione. Capisco l'idea alla base, ma voglio capire alcuni degli algoritmi che lo implementano - algoritmi che prendono come input un'immagine sfocata con la sua funzione di campionamento punto (sfocatura del kernel) e producono come output l'immagine latente.

Finora ho trovato l'algoritmo Richardson–Lucy in cui la matematica non sembra essere così difficile, ma non riesco a capire come funziona l'algoritmo attuale. Wikipedia dice:

Questo porta ad un'equazione per il quale può essere risolto iterativamente secondo ...

tuttavia non mostra il ciclo attuale. Qualcuno può indicarmi una risorsa in cui è spiegato l'algoritmo attuale. Su Google riesco solo a trovare metodi che utilizzano Richardson-Lucy come uno dei suoi passaggi, ma non l'algoritmo di Richardson-Lucy.

Algoritmo in qualsiasi linguaggio o pseudo-codice sarebbe bello, tuttavia se uno è disponibile in Python, sarebbe fantastico.

Grazie in anticipo.

Modifica

In sostanza quello che voglio capire è data un'immagine sfocata (nxm):

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

e il kernel (IXJ) che è stato utilizzato al fine di ottenere l'immagine sfocata:

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

Quali sono i passaggi esatti nell'algoritmo di Richardson-Lucy per poter calcolare l'immagine originale.

Answer 1

Ecco un'implementazione open source Python: http://code.google.com/p/iocbio/wiki/IOCBioMicroscope

Answer 2

L'equazione Wikipedia fornisce una funzione per l'iterazione t+1 in termini di iterazione t. È possibile implementare questo tipo di algoritmo iterativo nel seguente modo:

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 

Naturalmente, si dovrà implementare la propria funzione difference che sia corretta per qualsiasi tipo di dati che si sta lavorando. È inoltre necessario gestire il caso in cui la convergenza non viene mai raggiunta (ad esempio, limitare il numero di iterazioni).

Answer 3

Se aiuta: ecco un implementazione ho scritto che include una certa documentazione ....

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

Richardson Lucy è un blocco di costruzione per molti altri algoritmi di deconvoluzione. Ad esempio, l'esempio di iocbio di cui sopra ha modificato l'algoritmo per gestire meglio il rumore.

È un algoritmo relativamente semplice (come queste cose vanno) ed è un punto di partenza per algoritmi più complicati in modo da poter trovare molte diverse implementazioni.

Answer 4

Ecco una semplice implementazione 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') 

Puoi anche lavorare nel dominio della frequenza invece che nel dominio spaziale come sopra. In tal caso, il codice sarebbe:

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

L'unica cosa che non capisco è come questa inversione spaziale delle opere PSF e che cosa è per. Se qualcuno potesse spiegarmi che per me sarebbe bello! Sto anche cercando una semplice implementazione Matlab R-L per PSF spazialmente varianti (cioè funzioni di diffusione punti spaziale non omogenee) - se qualcuno ne avesse una, fatemelo sapere!

Per eliminare gli artefatti ai bordi, è possibile riflettere l'immagine di input ai bordi e ritagliare successivamente i bit con mirroring oppure utilizzare image = edgetaper(image, PSF) di Matlab prima di chiamare RL_deconv.

L'implementazione nativa di Matlab deconvlucy.m è un po 'più complicata btw - il codice sorgente di quello può essere trovato here e utilizza un accelerated version of the basic algorithm.