Abstract
讨论了高斯高通滤波器的原理，并实现了高斯高通滤波器：能够指定生成的二维函数的大小$M\times N$和高斯函数中心的二维位置。从原始图像中减去原图像以获得锐化图像，分析其与高斯高通结果的不同。调整高斯低通滤波器的方差，直到通过图像减法获得的结果与高高斯高通的滤波结果相似，并分析原因。
Technical discussion
截止频率处在距频率矩形中心为$D_0$的高斯高通滤波器的传递函数，如下： \(H(u,v)=1-e^{-D^2(u,v)/2D_0^2}\) 以下用GHPF指代高斯高通滤波器。GHPF得到的结果会比一般滤波器更加平滑，即使是对微小物体和细线条使用高斯滤波器滤波，结果也是比较清晰的。
Discussion of results
Highpass Filtering Using a Lowpass Image
(a) Subtract your image in Project 04-03(b) from the original to obtain a sharpened image, as in Eq. (4.4-14). You will note that the resulting image does not resemble the Gaussian highpass results in Fig. 4.26. Explain why this is so.
上图为$D_0=200$的时候，进行的比较。
说实话，两者几乎没有差别，事实上也应该如此。但是可以看到，用原图减去高斯低通滤波处理后的图像会比直接高斯高通滤波更加模糊，而且在取值相同的情况下，原图像减去滤波图像后的输出图像的边缘亮度要高于直接使用高斯滤波器的滤波图像。
但从公式的角度上来说，两者理论上应该是一样的，因为两者是等价的，可能是因为不同的计算方式导致了不同的浮点数误差造成的。
(b) Adjust the variance of your Gaussian lowpass filter until the result obtained by image subtraction looks similar to Fig. 4.26(c). Explain your result.
- $D_0=250$
- $D_0=200$
- $D_0=150$
- $D_0=100$
- $D_0=50$
可以看到在$D_0$取较小的值的时候锐化效果更为明显。
因此当$D_0$的值变小的时候，两幅图像会越来越接近。当$D_0=50$的时候，两幅图像已经接近一样了。
在低通滤波器中，当信号频率低于这个截止频率时，信号得以通过；当信号频率高于这个截止频率时，信号输出将被大幅衰减。这个截止频率即被定义为通带和阻带的界限。也就是说，截止频率越低，原图像减去滤波图像得到的锐化图像，锐化效果更好，边缘更加清晰。
Innovation points
实现了一个高斯高通滤波器，能够更方便地进行比较。
Appendix
主函数：
f = imread("1.tif");
g = GaussianHighPassFilter(f, 50);
h = GaussianLowPassFilter(f, 50);
f = mat2gray(f, [0 255]);
h = double(f) - h;
figure
subplot(1, 3, 1)
imshow(f)
title('原图')
subplot(1, 3, 2)
imshow(g)
title('高斯高通滤波锐化图像')
subplot(1, 3, 3)
imshow(h)
title('（原图-高斯低通滤波）图像')
高斯高通滤波器：
function [g] = GaussianHighPassFilter(f, D0, P, Q, u0, v0)
% 高斯高通滤波器
% f 是输入图像
% D0 是截止频率
% P 和 Q 是滤波器的大小
% u0 和 v0 是滤波器的中心
% 1.给定一幅大小为 M x N 的输入图像f(x,y),得到填充参数 P 和 Q
% 这里的缺省参数是P = 2M, Q = 2N, u0 = P / 2, v0 = Q / 2
[m,n] = size(f);
f = mat2gray(f, [0 255]);
if nargin == 2
P = 2 * m;
Q = 2 * n;
u0 = P / 2;
v0 = Q / 2;
end
% 2.对f(x, y)添加必要数量的0，形成大小为P x Q的填充后的图像fp(x, y)
fp = zeros(P, Q);
fp(1 : m, 1 : n) = f(1 : m, 1 : n);
% 3.用(-1)^[2 * (u0 * x / P + v0 * y / Q)]乘以fp(x, y)移到(u0, v0)
cx = 2 * u0 / P;
cy = 2 * v0 / Q;
for x = 1 : m
for y = 1 : n
fp(x, y) = double(fp(x, y) * (-1) ^ (cx * x + cy * y));
end
end
% 4.计算来自步骤3的图像的DFT，得到F(u,v)
F = fft2(fp, P, Q);
% 5.生成一个实的、对称的滤波函数H(u, v),其大小为P x Q，中心在（u0, v0）处
% 用阵列相乘形成乘积G(u, v) = H(u, v)F(u, v);即G(i, k) = H(i, k)F(i, k)
H = zeros(P, Q);
for u = 1 : P
for v = 1 : Q
D = (u - u0) ^ 2 + (v - v0) ^ 2;
H(u, v) = 1 - exp((-D) / (2 * D0 ^ 2));
end
end
G = H .* F;
% 6.得到处理后的图像
gp = ifft2(G);
gp = real(gp);
for x = 1 : m
for y = 1 : n
gp(x, y) = double(gp(x, y) * (-1) ^ (cx * x + cy * y));
end
end
% 7.通过从gp(x,y)的左上象限提取 M * N 区域，得到最终处理结果g(x,y)
g(1 : m, 1 : n) = gp(1 : m, 1 : n);