在Matlab中创建高通滤波器

Question

我正在尝试在Matlab中创建一个高通滤波器。我使用以下代码生成高斯核：

function kernel = compute_kernel(sigma,size)
[x,y] = meshgrid(-size/2:size/2,-size/2:size/2);
constant = 1/(2*pi*sigma*sigma);
kernel = constant*exp( -(y.^2 + x.^2 )/(2 * sigma * sigma));
kernel = (kernel - min(kernel(:)))./(max(kernel(:)) - min(kernel(:)));
end

接下来，我使用所创建的内核来为图像（变量im2）创建一个低通滤波器：

g = compute_kernel(9,101);
im2_low = conv2(im2,g,'same');

我理解的是，我可以将经过滤波的图像从原始图像(在频域中)中减去，以提取高频率，使之成为高通滤波器的等效。

F = fft2(im2_low);
IM2 = fft2(im2);
IM2_high = IM2 - F;

figure; fftshow(IM2_high);
im2_high = ifft2(IM2_high);
figure; imshow(im2_high,[]);

然而，这似乎有些问题。当我查看高通滤波图像时，它似乎是一个颜色反转模糊的图像，而不是在线上看到的边缘定义的图像。我不确定我的过程是否有误或者我是否只是使用了错误的高斯核数值。

- user3927312

Answer 1

任何你想用来维护图像特征的内核（即你不想得出某些东西的大小，而是让图像看起来像人类可识别的图像），你需要确保对内核进行一个操作：标准化。

你似乎已经尝试过这个操作，但你错误地解释了内核中的标准化含义。它们不需要在 [0-1] 范围内，它们的和应该等于 1。

因此，考虑到你的代码:

im2=imread('https://upload.wikimedia.org/wikipedia/en/2/24/Lenna.png');
im2=double(rgb2gray(im2));

sigma=9;
sizei=101;
[x,y] = meshgrid(-sizei/2:sizei/2,-sizei/2:sizei/2);
constant = 1/(2*pi*sigma*sigma);
kernel = constant*exp( -(y.^2 + x.^2 )/(2 * sigma * sigma));
%%%%%% NORMALIZATION
kernel=kernel/sum(kernel(:));
%%%%%%

im2_low = conv2(im2,kernel,'same');

F = fft2(im2_low);
IM2 = fft2(im2);
IM2_high = IM2 - F;


im2_high = ifft2(IM2_high);
figure; imshow(im2_high,[]);

但是，正如CrisLuengo所提到的，减法是一个在傅立叶域中不改变的操作，因此答案是

im2_high=im2-im2_low

- Ander Biguri

---

现在正在尝试。谢谢。 - user3927312

@user3927312 为了让自己信服，尝试使用colorbar绘制您的“im2_low”，并检查极值的值。如果您想要有意义地添加/减去它们，两个图像应该在相同的范围内。 - Ander Biguri

你的回答解决了我的问题。不过，我对你上面的评论有点难以理解。如果你能再详细解释一下，那会很有帮助。正如你所看到的，我对图像处理还很新。 - user3927312

在计算卷积核后，执行figure();imshow(im2_low,[]);colorbar;figure();imshow(im2,[]);colorbar;f。注意颜色条的值。在您的情况下，低通滤波器的值为数千，而原始图像仅到255。如果要保持值的范围，则需要将卷积核归一化为总和为1。 - Ander Biguri

我现在明白了。谢谢。 - user3927312

傅里叶域中的减法与空间域中的减法相同。调用fft和ifft是多余的。 - Cris Luengo

Answer 2

这是一个对简短问题的详细回答。如果想要学到更多，可以阅读下面的内容。

低通滤波器和高通滤波器都是线性滤波器。 线性滤波器可以在空间域通过卷积应用，也可以在频率域（即傅里叶域）中作为乘法应用。

事实上，在傅里叶域中，低通滤波器核与全通滤波器之间的差异就是高通滤波器：

high_pass_filter = identity_filter - low_pass_filter

身份过滤器是每个元素都为1的内核。通过乘法应用滤波器，因此

IM2 * high_pass_filter = IM2 * ( identity_filter - low_pass_filter )

这与

IM2 * high_pass_filter = IM2 - IM2 * low_pass_filter

在这里，与问题中一样，IM2是图像im2的傅里叶域表示；所有黄色边框中的内容都是指方程，但使用伪代码编写，并使用*符号表示乘法。

因此，OP希望在傅里叶域中应用低通滤波器并减去输入图像，以获得高通滤波后的图像。

然而，傅里叶变换的一个基本属性是它是线性变换。这意味着：

F(a*f + b*g) == a * F(f) + b * F(g)

（其中F（。）为傅里叶变换，a和b是常数，f和g是函数）。设置a = 1和b = -1，并将g设置为低通滤波后的图像，将f设置为输入图像，我们得到

F(im2 - im2_low) == F(im2) - F(im2_low)

也就是说，在空间域和傅里叶域中的减法是等效的。因此，如果在空间域中计算im2_low，则无需转到傅里叶域进行减法。这两个代码段产生相同的结果（精度上略有差异）：

F = fft2(im2_low);
IM2 = fft2(im2);
IM2_high = IM2 - F;
im2_high = ifft2(IM2_high);

im2_high = im2 - im2_low;

此外，卷积也是线性的。这意味着，如果您将上述方程中的F(.)视为卷积，则这些方程仍然成立。您可以进行以下操作：

conv(f, h) - f == conv(f, h) - conv(f, 1) == conv(f, h-1)

这直接导致了空间域中高通滤波器的定义：

g = - compute_kernel(9,101);
g(51,51) = g(51,51) + 1;
im2_high2 = conv2(im2,g,'same');

您会发现max(max(abs(im2_high-im2_high2)))得出的值非常接近0。

关于计算高斯滤波的注意事项：

本问题中发布的compute_kernel函数通过直接计算2D高斯来计算2D滤波器内核。生成的滤波器内核为101x101像素，这意味着计算卷积需要101*101*N个乘加操作（MADs），其中N是过滤图像中的像素数。然而，高斯滤波是可分离的，这意味着可以仅用101*2*N个MAD（少50倍！）获得相同的结果。此外，对于sigma=9，也可以使用较小的内核。

1. Gaussian kernel size:

   The Gaussian function never reaches zero, but it reaches very close to zero quite quickly. When cutting it off at 3*sigma, very little of it is lost. I find 3 sigma to be a good balance. In the case of sigma = 9, the 3 sigma cutoff leads to a kernel with 55 pixels (3*sigma * 2 + 1).

2. Gaussian separability:

   The multi-dimensional Gaussian can be obtained by multiplying 1D Gaussians together:

   exp(-(y.^2+x.^2)/(2*sigma*sigma)) == exp(-(x.^2)/(2*sigma*sigma)) * exp(-(y.^2)/(2*sigma*sigma))
   

   This leads to a much more efficient implementation of the convolution:

   conv(f,h1*h2) == conv( conv(f,h1), h2 )
   

   That is, convolving an image with a column filter h1 and then convolving the result with a row filter h2 is the same as convolving the image with a 2D filter h1*h2. In code:

   sigma = 9;
   sizei = ceil(3*sigma); % 3 sigma cutoff
   g = exp(-(-sizei:sizei).^2/(2*sigma.^2)); % 1D Gaussian kernel
   g = g/sum(g(:)); % normalize kernel
   im2_low = conv2(g,g,im2,'same');
   
   g2d = g' * g;
   im2_low2 = conv2(im2,g2d,'same');
   

   The difference is numerical imprecision:

   max(max(abs(im2_low-im2_low2)))
   

   ans =
      1.3927e-12
   

你可以在我的博客上找到更详细的高斯滤波描述，还有使用MATLAB图像处理工具箱时可能遇到的一些问题。

点击这里查看有关高斯滤波的更详细说明，点击这里了解使用MATLAB图像处理工具箱时可能会遇到的一些问题。

- Cris Luengo

---

太棒了的回答。 - Ander Biguri