频域滤波中默认的边界条件——补零与不补零（答作者问）

这个问题源自于Rafael Gonzalez的《数字图像处理》中的这幅图，为什么他频域滤波要将图像零延拓到二倍尺寸？

完全没有没要，既浪费计算，又浪费空间。

廖老师的问题是图像滤波涉及到源图像和滤波器相卷，卷积结果尺寸要大于源图像尺寸，因此要考虑边界处图像卷积处理，而通过延拓扩大的源图像对应的频域变换也就发生了相应变化。不同的延拓，对应的频域变换也是不一样的。

答：如果 $M$ 点信号与 $N$ 点冲激响应卷积，对于线性卷积而言，那结果的长度是 $M + N - 1$ ，如果要完全线性卷积，那就是需要补零延拓。然而对图像来说，我们通常需要滤波后结果与滤波前图像尺寸相同。所以，不用补零也可以，只是默认周期延拓，我更偏向于这样的结果。零延拓不会影响频谱（DTFT），而延拓也不是可选择的，是默认的。
在这里插入图片描述
注：图(d)错误，见这里。
这个他给的解释，并没有从理论高度做解释。

三个基础知识

卷积边界延拓的方式

在这里插入图片描述

离散傅里叶变换的周期性

在这里插入图片描述

线性卷积和循环卷积等效的条件

离散傅里叶变换的卷积定理

循环卷积对应于离散傅里叶变换的乘积，即离散傅里叶变换的卷积定理。

如果 $X(k) = DFT\{x(n)\}$ 和 $Y(k) = DFT\{y(n)\}$ ，那么

$DFT\{\{x(n)\} \circledast \{y(n)\}\}$

这里， $\circledast$ 表示由下式定义的循环卷积：

$\{x(n)\} \circledast \{y(n)\} = \sum_{m=0}^{N-1} x(m)y((n-m) \mod N).$

$\mod N)$ 表示循环移位。

Conventional Shift

在这里插入图片描述

Circular Shift

在这里插入图片描述

线性卷积是这样的， $h (- n)$ 移位，计算与 $f (n)$ 对位相乘再求和，一直下去就是线性卷积的结果 $g (n)$ 。
在这里插入图片描述
然而，不幸的是，离散傅里叶变换的周期性造成时域是循环卷积。某人优点很突出，缺点也很突出，我很欣赏他的优点，怎么办呢？那只能解决掉他的缺点。

零延拓，补零的边界条件

循环卷积没办法了，那就想办法让循环卷积与线性卷积等价，于是就有了零延拓。对信号补零到线性卷积结果的长度。图中虽然循环移位，但是由于信号 $f (n)$ 补零了，对应的那两个点就是零，相乘也是零。换句话说就是把位置给留出来。这样仍然是线性卷积的结果。
在这里插入图片描述

由于是线性卷积，那默认的边界延拓方式就是零延拓。看图说话，相当于 $f (n)$ 补零后周期延拓， $h (- n)$ 移位，那就是零延拓的线性卷积。
在这里插入图片描述

看图像的结果，这里看不出大小。

在这里插入图片描述
因为fft2本身可以做任意点的傅里叶变换，所以无需提前延拓，也就没显示。

空域滤波器，频域实现。

blurtype = 'gaussian';
sigma = 2;
hsize = 2 * ceil(3 * sigma) + 1;
psf_gauss = fspecial('gaussian', hsize, sigma);
offset = floor(size(psf_gauss)/2);
S = size(Img) + size(psf_gauss) - 1;
F = fft2(Img, S(1), S(2));
H = fft2(psf_gauss, S(1), S(2));G=F.* H;
g=im2uint8(real(ifft2(G)));

那Gonzalez为什么要补零为二倍呢？这是由于他在频域直接对模拟滤波器采样为 $M\times N$ 点，对应的空域就是 $M\times N$ 点（频域滤波器不这样设计，很离谱，暂且不论）。所以为了不产生混叠，延拓的图像尺寸至少是 $(2M-1)\times (2N-1)$ 。他的大错不少，小错不断，他认为是 $2M \times 2N$ 也就见怪不怪了。

频域滤波器

在这里插入图片描述

filename = 'cameraman';
suffix = '.tif';
Img = im2double(im2gray(imread([filename, suffix])));
[M,N]=size(Img);
D0 = 0.15;
[U,V] = freqspace([2*M 2*N],'meshgrid');
d = sqrt(U.^2+V.^2);
H = exp(-(d/(2*D0.^2)));F = fftshift(fft2(Img,2*M,2*N));
G = F.* H;
g=im2uint8(real(ifft2(ifftshift(G))));

不补零，周期延拓的边界条件

回头再看不零延拓的话，会怎么样？头部数据混叠，尾部数据丢失。
在这里插入图片描述
但是图像处理的滤波通常是考虑same size的结果，可以先做循环移位。

在这里插入图片描述
与线性卷积的结果比较，中间的五个值是相同的，这是因为边界受周期延拓影响了。

>> gcgc =1.4641    1.5738    0.9545    0.5790    0.3512    0.2130    0.8149
>> gg =0.7071    1.4289    1.5738    0.9545    0.5790    0.3512    0.2130    0.1078    0.0352

所以MATLAB提出了psf2otf函数，先做循环移位，再计算离散傅里叶变换。如果有空域的卷积核，通过这个函数实现频域滤波。相当于周期延拓的边界条件。

空域滤波器，频域实现。

在这里插入图片描述

[M,N]=size(Img);
F = fft2(Img,M,N);blurtype = 'gaussian';
sigma = 2;
hsize = 2 * ceil(3 * sigma) + 1;
offset = (hsize-1)/2;
psf_gauss = fspecial('gaussian', hsize, sigma);
H = psf2otf(psf_gauss, [M,N]); %这里使用fft2就会产生错误，很多刊物这里有误。
G=F.*H;
g=im2uint8(real(ifft2(G)));

注意psf2otf的地方。

频域滤波器

在这里插入图片描述

[M,N]=size(Img);
D0 = 0.3;
[U,V] = freqspace([M N],'meshgrid');
d = sqrt(U.^2+V.^2);
H = exp(-(d/(2*D0.^2)));F = fftshift(fft2(Img));
G = F.* H;g=im2uint8(real(ifft2(ifftshift(G))));