实验五 图像逆滤波处理
数字图像处理实验记录
实验五 图像逆滤波处理
利用逆滤波和其他逆卷积算法对运动模糊或散焦模糊图像进行图像复原，并给出实验结果。
实验步骤：对原图像进行聚焦模糊和运动模糊，之后利用逆滤波复原对图像进行复原
1、python代码见附录：
2、Matlab实验代码见附录；
实验结果：
1、python处理：
2、Matlab处理：
原图：
模糊处理：
逆滤波复原：
振幅谱：
附录：
1、Python代码：
import matplotlib.pyplot as graph
import numpy as np
from numpy import fft
import math
import cv2
def motion_process(image_size, motion_angle):
PSF = np.zeros(image_size)
print(image_size)
center_position = (image_size[0] - 1) / 2
print(center_position)
slope_tan = math.tan(motion_angle * math.pi / 180)
slope_cot = 1 / slope_tan
if slope_tan <= 1:
for i in range(15):
offset = round(i * slope_tan)
PSF[int(center_position + offset), int(center_position - offset)] = 1
return PSF / PSF.sum()
else:
for i in range(15):
offset = round(i * slope_cot)
PSF[int(center_position - offset), int(center_position + offset)] = 1
return PSF / PSF.sum()
def make_blurred(input, PSF, eps):
input_fft = fft.fft2(input)
PSF_fft = fft.fft2(PSF) + eps
blurred = fft.ifft2(input_fft * PSF_fft)
blurred = np.abs(fft.fftshift(blurred))
return blurred
def inverse(input, PSF, eps):
input_fft = fft.fft2(input)
PSF_fft = fft.fft2(PSF) + eps
result = fft.ifft2(input_fft / PSF_fft)
result = np.abs(fft.fftshift(result))
return result
def wiener(input, PSF, eps, K=0.01):
input_fft = fft.fft2(input)
PSF_fft = fft.fft2(PSF) + eps
PSF_fft_1 = np.conj(PSF_fft) / (np.abs(PSF_fft) ** 2 + K)
result = fft.ifft2(input_fft * PSF_fft_1)
result = np.abs(fft.fftshift(result))
return result
image = cv2.imread(‘Jellyfish.jpg’)
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
img_h = image.shape[0]
img_w = image.shape[1]
graph.figure(1)
graph.xlabel(“01Original Image”)
graph.gray()
graph.imshow(image)
graph.figure(2)
graph.gray()
PSF = motion_process((img_h, img_w), 60)
blurred = np.abs(make_blurred(image, PSF, 1e-3))
graph.subplot(231)
graph.xlabel(“Motion blurred”)
graph.imshow(blurred)
result = inverse(blurred, PSF, 1e-3)
graph.subplot(232)
graph.xlabel(“02inverse deblurred”)
graph.imshow(result)
result = wiener(blurred, PSF, 1e-3)
graph.subplot(233)
graph.xlabel(“03wiener deblurred(k=0.01)”)
graph.imshow(result)
blurred_noisy = blurred + 0.1 * blurred.std() *
np.random.standard_normal(blurred.shape)
graph.subplot(234)
graph.xlabel(“04motion & noisy blurred”)
graph.imshow(blurred_noisy)
result = inverse(blurred_noisy, PSF, 0.1 + 1e-3)
graph.subplot(235)
graph.xlabel(“05inverse deblurred”)
graph.imshow(result)
result = wiener(blurred_noisy, PSF, 0.1 + 1e-3)
graph.subplot(236)
graph.xlabel(“06wiener deblurred(k=0.01)”)
graph.imshow(result)
graph.show()
(562, 1000)
280.5
2、matlab代码：
f = imread(‘Jellyfish.jpg’);
f = mat2gray(f,[0 255]);
f_original = f;
[M,N] = size(f);
P = 2M;
Q = 2N;
fc = zeros(M,N);
for x = 1:1:M
for y = 1:1:N
fc(x,y) = f(x,y) * (-1)^(x+y);
end
end
F_I = fft2(fc,P,Q);
figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel(‘a).Original Image’);
subplot(1,2,2);
imshow(log(1 + abs(F_I)),[ ]);
xlabel(‘b).Fourier spectrum of a).’);
%% ------motion blur------------------
H = zeros(P,Q);
a = 0.02;
b = 0.02;
T = 1;
for x = (-P/2):1:(P/2)-1
for y = (-Q/2):1:(Q/2)-1
R = (xa + yb)pi;
if(R == 0)
H(x+(P/2)+1,y+(Q/2)+1) = T;
else H(x+(P/2)+1,y+(Q/2)+1) = (T/R)(sin®)exp(-1iR);
end
end
end
%% ------the atmospheric turbulence modle------------------
H_1 = zeros(P,Q);
k = 0.0025;
for x = (-P/2):1:(P/2)-1
for y = (-Q/2):1:(Q/2)-1
D = (x^2 + y2)(5/6);
D_0 = 60;
H_1(x+(P/2)+1,y+(Q/2)+1) = exp(-kD);
end
end
%% -----------noise------------------
a = 0;
b = 0.2;
n_gaussian = a + b . randn(M,N);
Noise = fft2(n_gaussian,P,Q);
figure();
subplot(1,2,1);
imshow(n_gaussian,[-1 1]);
xlabel(‘a).Gaussian noise’);
subplot(1,2,2);
imshow(log(1 + abs(Noise)),[ ]);
xlabel(‘b).Fourier spectrum of a).’);
%%
G = H .* F_I + Noise;
% G = H_1 .* F_I + Noise;
gc = ifft2(G);
gc = gc(1:1:M+27,1:1:N+27);
for x = 1:1:(M+27)
for y = 1:1:(N+27)
g(x,y) = gc(x,y) .* (-1)^(x+y);
end
end
gc = gc(1:1:M,1:1:N);
for x = 1:1:(M)
for y = 1:1:(N)
g(x,y) = gc(x,y) .* (-1)^(x+y);
end
end
figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel(‘a).Original Image’);
subplot(1,2,2);
imshow(log(1 + abs(F_I)),[ ]);
xlabel(‘b).Fourier spectrum of a).’);
figure();
subplot(1,2,1);
imshow(abs(H),[ ]);
xlabel(‘c).The motion modle H(u,v)(a=0.02,b=0.02,T=1)’);
subplot(1,2,2);
n = 1:1:P;
plot(n,abs(H(400,:)));
axis([0 P 0 1]);grid;
xlabel(‘H(n,400)’);
ylabel(’|H(u,v)|’);
figure();
subplot(1,2,1);
imshow(real(g),[0 1]);
xlabel(‘d).Result image’);
subplot(1,2,2);
imshow(log(1 + abs(G)),[ ]);
xlabel('e).Fourier spectrum of d). ');
%% --------------inverse_filtering---------------------
%F = G ./ H;
%F = G ./ H_1;
for x = (-P/2):1:(P/2)-1
for y = (-Q/2):1:(Q/2)-1
D = (x^2 + y2)(0.5);
if(D < 258)
F(x+(P/2)+1,y+(Q/2)+1) = G(x+(P/2)+1,y+(Q/2)+1) ./ H_1(x+(P/2)+1,y+(Q/2)+1);
% no noise D < 188
% noise D < 56
else F(x+(P/2)+1,y+(Q/2)+1) = G(x+(P/2)+1,y+(Q/2)+1);
end
end
end
% Butterworth_Lowpass_Filters
H_B = zeros(P,Q);
D_0 = 70;
for x = (-P/2):1:(P/2)-1
for y = (-Q/2):1:(Q/2)-1
D = (x^2 + y2)(0.5);
%if(D < 200) H_B(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/D_0)^100);end
H_B(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/D_0)^20);
end
end
F = F .* H_B;
f = real(ifft2(F));
f = f(1:1:M,1:1:N);
for x = 1:1:(M)
for y = 1:1:(N)
f(x,y) = f(x,y) * (-1)^(x+y);
end
end
%% ------show Result------------------
figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel(‘a).Result image’);
subplot(1,2,2);
imshow(log(1 + abs(F)),[ ]);
xlabel(‘b).Fourier spectrum of a).’);
figure();
n = 1:1:P;
plot(n,abs(F(400,:)),‘r-’,n,abs(F(400,:)),‘b-’);
axis([0 P 0 1000]);grid;
xlabel(‘Number of rows(400th column)’);
ylabel(‘Fourier amplitude spectrum’);
legend(‘F_{limit}(u,v)’,‘F(u,v)’);
figure();
n = 1:1:P;
plot(n,abs(H(400,:)),‘g-’);
axis([0 P 0 1]);grid;
xlabel(‘H’’{s}(n,400)’);
ylabel(’|H’’{s}(u,v)|’);
%% ----------Wiener filters-----------
% K = 0.000014;
K = 0.02;
%H_Wiener = ((abs(H_1).2)./((abs(H_1).2)+K)).(1./H_1);
H_Wiener = ((abs(H).2)./((abs(H).2)+K)).(1./H);
F_Wiener = H_Wiener .* G;
f_Wiener = real(ifft2(F_Wiener));
f_Wiener = f_Wiener(1:1:M,1:1:N);
for x = 1:1:(M)
for y = 1:1:(N)
f_Wiener(x,y) = f_Wiener(x,y) * (-1)^(x+y);
end
end
[SSIM_Wiener mssim] = ssim_index(f_Wiener,f_original,[0.01 0.03],ones(8),1);
SSIM_Wiener
%% ------show Result------------------
figure();
subplot(1,2,1);
%imshow(f_Wiener(1:128,1:128),[0 1]);
imshow(f_Wiener,[0 1]);
xlabel(‘d).Result image by Wiener filter’);
subplot(1,2,2);
imshow(log(1+abs(F_Wiener)),[ ]);
xlabel(‘c).Fourier spectrum of c).’);
% subplot(1,2,2);
% %imshow(f(1:128,1:128),[0 1]);
% imshow(f,[0 1]);
% xlabel(‘e).Result image by inverse filter’);
figure();
n = 1:1:P;
plot(n,abs(F(400,:)),‘r-’,n,abs(F_Wiener(400,:)),‘b-’);
axis([0 P 0 500]);grid;
xlabel(‘Number of rows(400th column)’);
ylabel(‘Fourier amplitude spectrum’);
legend(‘F(u,v)’,‘F_{Wiener}(u,v)’);
figure();
subplot(1,2,1);
imshow(log(1 + abs(H_Wiener)),[ ]);
xlabel(‘a).F_{Wiener}(u,v).’);
subplot(1,2,2);
n = 1:1:P;
plot(n,abs(H_Wiener(400,:)));
axis([0 P 0 80]);grid;
xlabel(‘Number of rows(400th column)’);
ylabel(‘Amplitude spectrum’);
%% ------------Constrained_least_squares_filtering---------
p_laplacian = zeros(M,N);
Laplacian = [ 0 -1 0;
-1 4 -1;
0 -1 0];
p_laplacian(1:3,1:3) = Laplacian;
P = 2M;
Q = 2N;
for x = 1:1:M
for y = 1:1:N
p_laplacian(x,y) = p_laplacian(x,y) * (-1)^(x+y);
end
end
P_laplacian = fft2(p_laplacian,P,Q);
F_C = zeros(P,Q);
r = 0.2;
H_clsf = ((H’)./((abs(H).^2)+r.*P_laplacian));
F_C = H_clsf .* G;
f_c = real(ifft2(F_C));
f_c = f_c(1:1:M,1:1:N);
for x = 1:1:(M)
for y = 1:1:(N)
f_c(x,y) = f_c(x,y) * (-1)^(x+y);
end
end
%%
figure();
subplot(1,2,1);
imshow(f_c,[0 1]);
xlabel(‘e).Result image by constrained least squares filter (r = 0.2)’);
subplot(1,2,2);
imshow(log(1 + abs(F_C)),[ ]);
xlabel(‘f).Fourier spectrum of c).’);
[SSIM_CLSF mssim] = ssim_index(f_c,f_original,[0.01 0.03],ones(8),1);
figure();
subplot(1,2,1);
imshow(log(1 + abs(H_clsf)),[ ]);
xlabel(‘a).F_{clsf}(u,v).’);
subplot(1,2,2);
n = 1:1:P;
plot(n,abs(H_clsf(400,:)));
axis([0 P 0 80]);grid;
xlabel(‘Number of rows(400th column)’);
ylabel(‘Amplitude spectrum’);