低通巴特沃斯Lbtfilter(image_in, D0, n),高通巴特沃斯Hbtfilter(image_in, D0, n),两者只有传递函数不同,其余部分相同。 下列程序仅展示低通巴特沃斯。
function [image_out,r] = Lbtfilter(image_in, D0, n)
% 输入为需要进行滤波的灰度图像,Butterworth滤波器的截止频率D0,阶数N
% r返回功率谱比
[x, y] = size(image_in);
P = 2 * x;Q = 2 * y;
f1 = zeros(P, Q);
image_in=double(image_in);
%对图像填充0,并且移到变换中心
for i = 1 : x
for j = 1 : y
f1(i, j) =image_in(i, j) * (-1)^(i+j);
end
end
% 对填充后的图像进行傅里叶变换
F1 = fft2(f1);
% 生成Butterworth滤波函数
H = zeros(P, Q);G=zeros(P, Q);
a = D0^(2 * n);
for u = 1 : P
for v = 1 : Q
D1 = (u-(x+1.0))^2 + (v-(y+1.0))^2;
H(u, v) = 1 / (1 + (D1^n) / a);
%高通巴特沃斯
%H(u, v) = 1 / (1 + a/(D1^n));
G(u,v)=F1(u,v)*H(u,v);
end
end
% 反傅里叶变换
Gn=ifft2(G);
image_out = zeros(x, y, 'uint8');
Gn = real(Gn);
g = zeros(x, y);
%中心频率移到左上角
for i = 1 : x
for j = 1 : y
g(i, j) = Gn(i, j) * (-1)^(i+j);
end
end
%归一化
mmax = max(g(:));mmin = min(g(:));
for i = 1 : x
for j = 1 : y
image_out(i,j) = uint8(255 * (g(i, j)-mmin) / (mmax-mmin));
end
end
%%计算功率谱比
s=0;s1=0
for u=1:P
for v=1:Q
L1=(abs(G(u,v)))^2;
s1=s1+L1;
L=(abs(F1(u,v)))^2;
s=s+L;
end
end
r=s1/s;
end
低通高斯 Lgaussfilter(image_in, D0),高通高斯 Hgaussfilter(image_in, D0),两者只有传递函数不同,其余部分相同。 下列程序仅展示低通高斯。
function [image_out,r] = Lgaussfilter(image_in, D0)
% 输入为需要进行滤波的灰度图像,Butterworth滤波器的截止频率D0
% r为功率谱比
[x, y] = size(image_in);
P = 2 * x;
Q = 2 * y;
f1 = zeros(P, Q);
image_in=double(image_in);
%对图像填充0,并且移到变换中心
for i = 1 : x
for j = 1 : y
f1(i, j) =image_in(i, j) * (-1)^(i+j);
end
end
% 对填充后的图像进行傅里叶变换
F1 = fft2(f1);
% 生成Butterworth滤波函数
H = zeros(P, Q);G=zeros(P, Q);
a = D0^2 ;
for u = 1 : P
for v = 1 : Q
D1 = (u-(x+1.0))^2 + (v-(y+1.0))^2;
H(u, v) = exp(-D1/2/a);
%高通高斯
%H(u, v) = 1-exp(-D1/2/a);
G(u,v)=F1(u,v)*H(u,v);
end
end
% 反傅里叶变换
Gn=ifft2(G);
image_out = zeros(x, y, 'uint8');
Gn = real(Gn);
%中心频率移到左上角
g = zeros(x, y);
for i = 1 : x
for j = 1 : y
g(i, j) = Gn(i, j) * (-1)^(i+j);
end
end
mmax = max(g(:));mmin = min(g(:));
for i = 1 : x
for j = 1 : y
image_out(i,j) = uint8(255 * (g(i, j)-mmin) / (mmax-mmin));
end
end
%%计算功率谱比
s=0;s1=0
for u=1:P
for v=1:Q
L1=(abs(G(u,v)))^2;
s1=s1+L1;
L=(abs(F1(u,v)))^2;
s=s+L;
end
end
r=s1/s;
end
function [image_out] = Laplace(image_in, c)
% 输入为需要进行滤波的灰度图像,Butterworth滤波器的截止频率D0,阶数N
[x, y] = size(image_in);
P = 2 * x;
Q = 2 * y;
f1 = zeros(P, Q);
image_in=double(image_in);
%对图像填充0,并且移到变换中心
for i = 1 : x
for j = 1 : y
f1(i, j) =image_in(i, j) * (-1)^(i+j);
end
end
% 对填充后的图像进行傅里叶变换
F1 = fft2(f1);
% 生成Butterworth滤波函数
H = zeros(P, Q);G=zeros(P, Q);
for u = 1 : P
for v = 1 : Q
D1 = (u-(x+1.0))^2 + (v-(y+1.0))^2;
H(u,v)=1+c*4*pi^2*D1;
G(u,v)=F1(u,v)*H(u,v);
end
end
% 反傅里叶变换
Gn=ifft2(G);
image_out = zeros(x, y, 'uint8');
Gn = real(Gn);
g = zeros(x, y);
for i = 1 : x
for j = 1 : y
g(i, j) = Gn(i, j) * (-1)^(i+j);
end
end
mmax = max(g(:));mmin = min(g(:));
for i = 1 : x
for j = 1 : y
image_out(i,j) = uint8(255 * (g(i, j)-mmin) / (mmax-mmin));
end
end
end
function [image_out] = Unmask(image_in, D0, k1,k2)
% 输入为需要进行滤波的灰度图像,Butterworth滤波器的截止频率D0,阶数N
[x, y] = size(image_in);
P = 2 * x;
Q = 2 * y;
f1 = zeros(P, Q);
image_in=double(image_in);
%对图像填充0,并且移到变换中心
for i = 1 : x
for j = 1 : y
f1(i, j) =image_in(i, j) * (-1)^(i+j);
end
end
% 对填充后的图像进行傅里叶变换
F1 = fft2(f1);
% 生成Butterworth滤波函数
H = zeros(P, Q);G=zeros(P, Q);
a = D0^2 ;
for u = 1 : P
for v = 1 : Q
D1 = (u-(x+1.0))^2 + (v-(y+1.0))^2;
H(u,v)=1-exp(-D1/(2*a));
G(u,v)=(k1+k2*H(u,v))*F1(u,v);
end
end
% 反傅里叶变换
Gn=ifft2(G);
image_out = zeros(x, y, 'uint8');
Gn = real(Gn);
g = zeros(x, y);
for i = 1 : x
for j = 1 : y
g(i, j) = Gn(i, j) * (-1)^(i+j);
end
end
%归一化
mmax = max(g(:));mmin = min(g(:));
for i = 1 : x
for j = 1 : y
image_out(i,j) = uint8(255 * (g(i, j)-mmin) / (mmax-mmin));
end
end
end
clc;clear;
test1=imread('C:\Users\杨德宇\Desktop\数字图像处理\第5次作业\test1.pgm');
test2=imread('C:\Users\杨德宇\Desktop\数字图像处理\第5次作业\test2.tif');
%%低通butterworth
[btest1_1,b11]=Lbtfilter(test1,25,2);[btest2_1,b21]=Lbtfilter(test2,25,2);
[btest1_2,b12]=Lbtfilter(test1,50,2);[btest2_2,b22]=Lbtfilter(test2,50,2);
[btest1_3,b13]=Lbtfilter(test1,75,2);[btest2_3,b23]=Lbtfilter(test2,75,2);
%%低通gauss
[gtest1_1,g11]=Lgaussfilter(test1,25);[gtest2_1,g21]=Lgaussfilter(test2,25);
[gtest1_2,g12]=Lgaussfilter(test1,50);[gtest2_2,g22]=Lgaussfilter(test2,50);
[gtest1_3,g13]=Lgaussfilter(test1,75);[gtest2_3,g23]=Lgaussfilter(test2,75);
%%绘图1
figure(1);imshow(test1);title('原图');
figure(2);imshow(btest1_1);title('D0=25 低通butterworth');
figure(3);imshow(btest1_2);title('D0=50 低通butterworth');
figure(4);imshow(btest1_3);title('D0=75 低通butterworth');
figure(5);imshow(btest2_1);title('D0=25 低通butterworth');
figure(6);imshow(btest2_2);title('D0=50 低通butterworth');
figure(7);imshow(btest2_3);title('D0=75 低通butterworth');
figure(8);imshow(test2);title('原图');
%%绘图高斯
figure(9);imshow(gtest1_1);title('D0=25 低通gauss');
figure(10);imshow(gtest1_2);title('D0=50 低通gauss');
figure(11);imshow(gtest1_3);title('D0=75 低通gauss');
figure(12);imshow(gtest2_1);title('D0=25 低通gauss');
figure(13);imshow(gtest2_2);title('D0=50 低通gauss');
figure(14);imshow(gtest2_3);title('D0=75 低通gauss');
clc;clear;
test3=imread('C:\Users\杨德宇\Desktop\数字图像处理\第5次作业\test3_corrupt.pgm');
test4=imread('C:\Users\杨德宇\Desktop\数字图像处理\第5次作业\test4 copy.bmp');
%%高通butterworth
[btest1_1,b11]=Hbtfilter(test3,25,2);[btest2_1,b21]=Hbtfilter(test4,25,2);
[btest1_2,b12]=Hbtfilter(test3,50,2);[btest2_2,b22]=Hbtfilter(test4,50,2);
[btest1_3,b13]=Hbtfilter(test3,75,2);[btest2_3,b23]=Hbtfilter(test4,75,2);
%%高通gauss
[gtest1_1,g11]=Hgaussfilter(test3,25);[gtest2_1,g21]=Hgaussfilter(test4,25);
[gtest1_2,g12]=Hgaussfilter(test3,50);[gtest2_2,g22]=Hgaussfilter(test4,50);
[gtest1_3,g13]=Hgaussfilter(test3,75);[gtest2_3,g23]=Hgaussfilter(test4,75);
%%绘图1
figure(2);imshow(btest1_1);title('D0=25 高通butterworth');
figure(3);imshow(btest1_2);title('D0=50 高通butterworth');
figure(4);imshow(btest1_3);title('D0=75 高通butterworth');
figure(5);imshow(btest2_1);title('D0=25 高通butterworth');
figure(6);imshow(btest2_2);title('D0=50 高通butterworth');
figure(7);imshow(btest2_3);title('D0=75 高通butterworth');
%%绘图高斯
figure(9);imshow(gtest1_1);title('D0=25 高通gauss');
figure(10);imshow(gtest1_2);title('D0=50 高通gauss');
figure(11);imshow(gtest1_3);title('D0=75 高通gauss');
figure(12);imshow(gtest2_1);title('D0=25 高通gauss');
figure(13);imshow(gtest2_2);title('D0=50 高通gauss');
figure(14);imshow(gtest2_3);title('D0=75 高通gauss');
clear;clc;
test3=imread('C:\Users\杨德宇\Desktop\数字图像处理\第5次作业\test3_corrupt.pgm');
test4=imread('C:\Users\杨德宇\Desktop\数字图像处理\第5次作业\test4 copy.bmp');
%拉普拉斯
test31=Laplace(test3,1);test32=Laplace(test3,2);test33=Laplace(test3,0.5);
test41=Laplace(test4,1);test42=Laplace(test4,2);test43=Laplace(test4,0.5);
%Unmask
test3u1=Unmask(test3,75,1,1);test3u2=Unmask(test3,75,1,2);test3u3=Unmask(test3,75,1,0.5);
test4u1=Unmask(test4,75,1,1);test4u2=Unmask(test4,75,1,2);test4u3=Unmask(test4,75,1,0.5);
%绘图
figure(1)
subplot(2,2,1);imshow(test3);title('原图');
subplot(2,2,2);imshow(test33);title('拉普拉斯滤波 c=0.5');
subplot(2,2,3);imshow(test31);title('拉普拉斯滤波 c=1');
subplot(2,2,4);imshow(test32);title('拉普拉斯滤波 c=2');
figure(2)
subplot(2,2,1);imshow(test4);title('原图');
subplot(2,2,2);imshow(test43);title('拉普拉斯滤波 c=0.5');
subplot(2,2,3);imshow(test41);title('拉普拉斯滤波 c=1');
subplot(2,2,4);imshow(test42);title('拉普拉斯滤波 c=2');
figure(3)
subplot(2,2,1);imshow(test3);title('原图');
subplot(2,2,2);imshow(test3u1);title('Unmask k1=1 k2=1');
subplot(2,2,3);imshow(test3u2);title('Unmask k1=1 k2=2');
subplot(2,2,4);imshow(test3u3);title('Unmask k1=1 k2=0.5');
figure(4)
subplot(2,2,1);imshow(test4);title('原图');
subplot(2,2,2);imshow(test4u1);title('Unmask k1=1 k2=1');
subplot(2,2,3);imshow(test4u2);title('Unmask k1=1 k2=2');
subplot(2,2,4);imshow(test4u3);title('Unmask k1=1 k2=0.5');