Вот очень простая реализация 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')
Вы также можете работать в частотной области, а не в пространственной области, как указано выше. В этом случае код будет:
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);
Чтобы избавиться от артефактов на краях, вы можете отразить входное изображение по краям, а затем обрезать зеркальные биты или использовать Matlab image = edgetaper(image, PSF) перед вызовом RL_deconv.
Собственная реализация Matlab deconvlucy.m немного сложнее, кстати, ее исходный код можно найти здесь и использует ускоренную версию базового алгоритма .