¿Cómo funciona el algoritmo Richardson-Lucy? ¿Ejemplo de código?

Question

Estoy tratando de descubrir cómo funciona la deconvolución. Entiendo la idea detrás de esto, pero quiero entender algunos de los algoritmos reales que lo implementan: algoritmos que toman como entrada una imagen borrosa con su función de muestra de punto (desenfocar kernel) y producen como salida la imagen latente.

Hasta ahora encontré el algoritmo Richardson–Lucy donde las matemáticas no parecen ser tan difíciles, pero no puedo entender cómo funciona el algoritmo real. En Wikipedia dice:

Esto lleva a una ecuación para que puede resolverse de forma iterativa de acuerdo ...

sin embargo, no muestra el bucle real. ¿Alguien puede dirigirme a un recurso donde se explica el algoritmo real? En Google, solo logro encontrar métodos que usan Richardson-Lucy como uno de sus pasos, pero no el algoritmo real de Richardson-Lucy.

Algoritmo en cualquier idioma o pseudocódigo sería bueno, sin embargo, si uno está disponible en Python, sería increíble.

Gracias por adelantado.

Editar

Esencialmente lo que quiero averiguar es dada imagen borrosa (nxm):

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

y el núcleo (ixj) que se utilizó con el fin de obtener la imagen borrosa:

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

¿Cuáles son los pasos exactos en el algoritmo de Richardson-Lucy con el fin de averiguar la imagen original.

Answer 1

Aquí es una implementación de código abierto Python: http://code.google.com/p/iocbio/wiki/IOCBioMicroscope

Answer 2

La ecuación en Wikipedia ofrece una función para la iteración t+1 en términos de iteración t. Puede implementar este tipo de algoritmo iterativo de la siguiente manera:

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 

Por supuesto, usted tendrá que poner en práctica su propia función difference que es correcto para cualquier tipo de datos que se está trabajando. También debe manejar el caso donde nunca se alcanza la convergencia (por ejemplo, limitar el número de iteraciones).

Answer 3

Si se ayuda a que aquí hay una aplicación que escribí que incluye algún tipo de documentación ....

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

Richardson Lucy es un bloque de construcción para muchos otros algoritmos de deconvolución. Por ejemplo, el ejemplo iocbio anterior modificó el algoritmo para tratar mejor el ruido.

Es un algoritmo relativamente simple (como funcionan estas cosas) y es un punto de partida para algoritmos más complicados, por lo que puede encontrar muchas implementaciones diferentes.

Answer 4

Aquí es una aplicación muy simple 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') 

También puede trabajar en el dominio de la frecuencia en lugar de en el dominio espacial como se indicó anteriormente. En ese caso, el código sería:

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

Lo único que no entiendo muy bien es cómo esta inversión espacial de las obras de PSF y para qué sirve. ¡Si alguien pudiera explicar eso para mí eso sería genial! También estoy buscando una implementación simple de Matlab R-L para PSF espacialmente variantes (es decir, funciones de dispersión de puntos espacialmente no homogéneas). Si alguien tuviera una, por favor, ¡hágamelo saber!

Para deshacerse de los artefactos en los bordes, puede reflejar la imagen de entrada en los bordes y luego recortar los bits reflejados después o utilizar image = edgetaper(image, PSF) de Matlab antes de llamar al RL_deconv.

La implementación deconlablucy.m nativa de Matlab es un poco más complicada, por cierto, el código fuente de ese se puede encontrar here y usa un accelerated version of the basic algorithm.