image restoration matlab code

270 次查看(过去 30 天)
hello i am trying to implement this code to get the result in figure 5.26(b) but still something wrong could you help me in that
thanks,
refernce: digital image processing gonzalez third edition
here is my code
c = im2double(imread('Fig0526(a)(original_DIP).tif'));
figure,imshow(c)
a = 0.1;
b = 0.1;
T=1;
[M, N] = size(c);
h = zeros(M,N);
for u = 1:M
for v = 1:N
h(u,v) = (T/(pi*(u*a + v*b)))*...
sin(pi*(u*a + v*b))*...
exp(-1i*pi*(u*a + v*b));
end
end
cf=(fft2(c));
c1 = cf.*h;
c1a=(ifft2(c1));
figure,imshow(abs(log(c1a)+1), []);

采纳的回答

Image Analyst
Image Analyst 2021-12-2
The origin was not at the right place. Try this:
clc; % Clear the command window.
fprintf('Beginning to run %s.m ...\n', mfilename);
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;
grayImage = im2double(imread('logo.tif'));
% Make sure it's 688 like the book says.
grayImage = imresize(grayImage, [688, 688]);
subplot(2, 3, 1);
imshow(grayImage, [])
impixelinfo
axis('on', 'image')
title('Original Image', 'FontSize', fontSize)
% Make a filter for the spectral domain.
a = 0.1;
b = a;
T=1;
[rows, columns] = size(grayImage);
h = zeros(rows, columns);
for u = 1 : rows
for v = 1 : columns
rx = v - columns/2;
ry = u - rows/2;
h(u,v) = (T/(pi*(ry*a + rx*b))) *...
sin(pi*(ry*a + rx*b))*...
exp(-1i*pi*(ry*a + rx*b));
end
end
% Print the max value of h.
absh = real(h);
fprintf('Max h = %f.\n', max(absh(:)))
subplot(2, 3, 2);
imshow(log(absh), []);
impixelinfo;
axis('on', 'image')
title('h Frequency Domain Filter', 'FontSize', fontSize)
% Compute fft of the image.
fftOrig = fftshift(fft2(grayImage));
subplot(2, 3, 3);
imshow(log(abs(fftOrig)), []);
impixelinfo;
axis('on', 'image')
title('FFT of original image', 'FontSize', fontSize)
% Multiply it by the filter. fftOrig and h have DC position at center of image.
c1 = fftOrig .* h;
% There are some nan's which mess up the inverse transform. Set nan's to zero.
rc = real(c1);
ic = imag(c1);
rc(isnan(rc)) = 0;
ic(isnan(ic)) = 0;
c1 = rc + 1i * ic;
subplot(2, 3, 4);
imshow(abs(log(c1)+1), []);
title('Filtered Spectrum', 'FontSize', fontSize)
% Inverse transform from frequency domain back to space domain.
% After multiplication, shift back.
c1Spatial = ifft2(ifftshift(c1));
subplot(2, 3, 5);
imshow(real(c1Spatial), []);
title('Filtered Image', 'FontSize', fontSize)
g = gcf;
g.WindowState = 'maximized'
  2 个评论
aziz alfares
aziz alfares 2021-12-2
thanks a lot it worked just fine
i have another question if you dont mind here i am trying to implement inverse filtering and weiner filtering
in this figure both of gaussian noise and motion blure were added after that both of uinverse filtering and weiner filtering were applied
note : the value of mean and vaiance are 0 and 650 respectively but when i enter the 650 to the variance it doesnt give me any result so i changed it to 0.01
this is where i came so far
I = im2double(imread('Fig0526(a)(original_DIP).tif'));
imshow(I);
title('Original Image (courtesy of MIT)');
LEN = 90;
THETA = 45;
PSF = fspecial('motion', LEN, THETA);
blurred = imfilter(I, PSF, 'conv', 'circular');
figure, imshow(blurred)
noise_mean = 0;
noise_var =0.01 ;
blurred_noisy = imnoise(blurred, 'gaussian', noise_mean, noise_var);
figure, imshow(blurred_noisy)
title('Simulate Blur and Noise')
estimated_nsr = noise_var / var(I(:));
wnr3 = deconvwnr(blurred_noisy, PSF, estimated_nsr);
figure, imshow(wnr3)
title('Restoration of Blurred, Noisy Image Using Estimated NSR');
vaishnavi padhar
vaishnavi padhar 2022-4-13
can you please provide the code for inverse filtering

请先登录,再进行评论。

更多回答(2 个)

Image Analyst
Image Analyst 2021-12-2
Why are you doing
exp(-1i*pi*(u*a + v*b));
??? Just use fft2().
  1 个评论
aziz alfares
aziz alfares 2021-12-2
the degradtion function as shown in (5.6-11) must be applied so i am trying to use the same equation to get the result

请先登录,再进行评论。


vaishnavi padhar
vaishnavi padhar 2022-4-13
Can you please provide the block diagram for the same.

产品


版本

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by