Reconstruction of a blurred and noisy image

71 visualizaciones (últimos 30 días)
Benedetta
Benedetta el 2 de Abr. de 2024 a las 10:40
Comentada: Benedetta el 9 de Abr. de 2024 a las 14:04
I have a blurred and noisy image, which is smeared by an isotropic Gaussian point-spreadfunction (rms-variance pixels) and corrupted by different type of noise (data are stored in the ASCII file blurredimageRplusGaussianplusSPnoise.dat contained in the .zip file, consisting in a rectangular matrix or array I = [ ] of size 387 x 516; the values account for the grey levels of pixels.
I have already applied a appropriate filter to eliminate the salt-and-pepper noise.
But then I have to apply one or more filters, in a single step or by implementing an iterative scheme, to reconstruct the uncorrupted image (for iterative schemes, I can monitor the improvement at each step i by measuring the total power P or total variance of the reconstructed image) and what i get is a huge explosion of the noise! How can I fix this problem? I tried to remove the noise, but the the image is too blurred and I can't recover it...
This is the code:
%% Step 1: Read the data from the ASCII file into MATLAB
load("blurredimageRplusGaussianplusSPnoise.dat")
image_data = blurredimageRplusGaussianplusSPnoise;
clear blurredimageRplusGaussianplusSPnoise
%% Remove salt&pepper
J = medfilt2(image_data);
% imshow(J)
figure(1)
imshowpair(image_data,J,'montage')
%% Remove noise and deblurring
% filter
epsilon = 1e-5;
max_iterations = 500;
denoised_image = wiener2(J, [5 5]);
for iter = 1:max_iterations
denoised_image = wiener2(denoised_image, [5 5]);
% Monitor improvement using total variance
if iter > 1
variance_change = abs(variance_prev - var(denoised_image(:)));
if variance_change < epsilon
break; % Stop if improvement is below threshold
end
end
variance_prev = var(denoised_image(:));
if iter == max_iterations
disp('max iterations reached')
end
end
J = denoised_image;
figure(2)
imshow(J)
PSF = fspecial('gaussian',[5 5],4);
% PSF = fspecial('disk',8);
Y = fft2(J);
% zero pad the psf to match the size of the blurred image
newh = zeros(size(J));
psfsize = size(PSF);
newh(1:psfsize(1),1:psfsize(2)) = PSF;
H = fft2(newh); % Fourier transform of the noise (we have 0 noise)
% use simple X = Y/H to get back original image
% show how much noise affects it
deblurred = ifft2(Y./H); % inverse Fourier transform
% plot deblurred image
figure(3)
imshow(deblurred)
  2 comentarios
DGM
DGM el 2 de Abr. de 2024 a las 19:06
The rest I don't know about, but to start with, maybe don't use a simple local median filter. Use something that actually does noise discrimination. That way you're only filling in the damaged pixels, without changing everything in between.
This is a list of threads containing examples to both fixed and adaptive noise removal filters based on local median (and other order statistics).
Benedetta
Benedetta el 9 de Abr. de 2024 a las 14:04
Thank you!

Iniciar sesión para comentar.

Respuesta aceptada

Image Analyst
Image Analyst el 2 de Abr. de 2024 a las 19:47
When I was in grad school, decades ago, my professor and I published a paper on "median window enhancement". The theory is that you can sharpen the image (steepen edges) but then that edge enhancement produces ringing (oscillatory noise). But then you can clip off the ringing, without affecting the sharpness of the edge by running the image through a median filter. Then just keep going for as many iterations as you want.
If you know the PSF, or just want to try trial and error to see which PSF width looks best, you can try built-in restoration routines that undo the blur and reduce noise. See Image Processing Toolbox functions like deconvreg, deconvwnr, deconvlucy, and deconvblind.
  1 comentario
DGM
DGM el 2 de Abr. de 2024 a las 23:42
I'm just going to throw this comment here, since I'm deferring to IA on the topic, but normally when I see something like
PSF = fspecial('gaussian',[5 5],4);
A flag goes off in my brain, telling me that this kernel doesn't have an appropriate amount of support. I don't know what that means for the task of deconvolution.

Iniciar sesión para comentar.

Más respuestas (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by