图像处理学习笔记-04-频率域滤波03

下式是二维傅里叶变换$DFT$

$$F(\mu,v) = \sum_{x = 0}^{M - 1}\sum_{y = 0}^{N - 1}f(x,y)e^{-j2\pi(\mu x / M + vy / N)}$$

其中$f(x,y)$是大小为$M \times N$的数字图像，当给出$F(\mu,v)$的呀得到傅里叶反变换$IDFT$：

$$f(x,y) = \frac{1}{MN}\sum_{\mu = 0}^{M - 1}\sum_{v = 0}^{N - 1}F(\mu,v)e^{j2\pi(\mu x / M + vy / N)}$$

二维离散傅里叶变换的一些性质

空间和频率间隔的关系

和一维的类似，假设对连续函数$f(t,z)$取样生成了一幅数字图像$f(x,y)$，分别在$t$和$z$方向所取的$M \times N$个样本点组成，$\Delta T,\Delta Z$表示样本间的间隔：

$$\Delta u = \frac{1}{M \Delta T},\Delta v = \frac{1}{N \Delta Z}$$

平移和旋转

• 平移：用指数项乘以$f(x,y)$将使DFT的原点一道点$(\mu_0,v_0)$，反之，用负指数乘以$F(\mu,v)$将使$f(x,y)$的原点移到$(x_0,y_0)$，平移不影响幅度，原图像平移之后对应的傅里叶变换和原图像的傅里叶变换是相同的

$$f(x,y)e^{j2\pi(\mu_0x/M + v_0y / N)} \Leftrightarrow F(\mu - \mu_0,v - v_0) \\ f(x - x_0,y - y_0) \Leftrightarrow F(\mu,v)e^{-j2\pi(\mu_0x/M + v_0y / N)}$$

• 旋转：
  使用极坐标：

$$x = r\cos \theta,y = r\sin\theta,y = \omega\cos\varphi,v = \omega\sin\varphi$$

可以得到：

$$f(r,\theta + \theta_0) \Leftrightarrow F(\omega,\varphi + \varphi_0)$$

也就是，若$f(x,y)$旋转$\theta_0$角度，那么$F(\mu,v)$也旋转相同的角度，反之相同

周期性

二维傅里叶变换及其反变换在$\mu,v$方向是无限周期的，即：

$$F(u,v) = F(u + k_1M,v) = F(u, v + k_2N) = F(u + k_1M,v + k_2N) \\ f(x,y) = f(x + k_1M,y) = f(x,y+k_2N) = f(x + k_1M,y + k_2N)$$

先考虑一维的情况，在区间$[0,M-1]$中，变换数据由两个在点$M / 2$处碰面的背靠背的半个周期组成，在该区间内有一个完整的周期会更加方便：

$$f(x)e^{j2\pi(ux_0/M)} \Leftrightarrow F(u - u_o)$$

假设$u_0 = M / 2$那么指数项变为$e^{j\pi x}$，又因为$x$是整数，所以指数项为$(-1)^x$：

$$f(x)(-1)^x \Leftrightarrow F(u - M / 2)$$

即用$(-1)^x$乘以$f(x)$将位于原点的数据$F(0)$移动到区间$[0,M - 1]$的中心位置，在二维情况下原理是一样的：

$$f(x,y)(-1)^{x + y} \Leftrightarrow F(u - M / 2,v - N / 2)$$

对称性：

任何实函数或虚函数$w(x,y)$可表示为一个奇数部分和一个偶数部分的和：

$$w(x,y) = w_e(x,y) + w_o(x,y)$$

其中偶数部分和奇数部分定义为：

$$w_e(x,y) = \frac{w(x,y) + w(-x,-y)}{2} \\ w_o(x,y) = \frac{w(x,y) - w(-x,-y)}{2}$$

偶函数是对称的，奇函数是反对称的，因为：

$$w_e(x,y) = w_e(-x,-y) \\ w_o(x,y) = -w_o(-x,-y)$$

因为$DFT,IDFT$中所有指数都是正的，所以是关于序列中点的对称(反对称)：

$$w_e(x,y) = w_e(M - x,N - y) \\ w_o(x,y) = -w_o(M - x,N - y)$$

一个奇函数和一个偶函数的乘积是奇函数，而离散函数为奇函数的唯一方法就是所有的样本的和为0：

$$\sum_{x = 0}^{M - 1}\sum_{y = 0}^{N - 1}w_e(x,y)w_o(x,y) = 0$$

离散函数例如：

$$f = \{f(0)\ f(1)\ f(2)\ f(3)\}$$

其中$M = 4$，可以看到，如果需要该函数是偶函数那么只需要$f(1) = f(3)$，如果需要该函数为奇函数，那么需要满足条件$f(0) = f(2) = 0,f(1) = -f(3)$
实函数$f(x,y)$的傅里叶变换是共轭对称的：

$$F^*(u,v) = F(-u,-v)$$

虚函数$f(x,y)$是共轭反对称的，$F^*(-u,-v) = -F(u,v)$
下标列出了傅里叶变换的对称性和相关的性质，双箭头表示傅里叶变换对，即对于表中的任意一行，右边的性质可由拥有左边性质的函数的傅里叶变换来满足，反之亦然，例如对于第三条性质，读作：如果$f(x,y)$是实函数，则其$DFT$的实部是偶函数，其虚部是奇函数，类似的，如果一个$DFT$分别具有偶函数的实部和奇函数的虚部，则其$IDFT$是一个实函数

傅里叶谱和相角

通常二维$DFT$是复函数，可以使用极坐标表示：

$$F(u,v) = |F(u,v)|e^{j\phi (u,v)}$$

$R,I$分别是函数$F(u,v)$的实部和虚部，其中，幅度：

$$|F(u,v)| = [R^2(u,v) + I^2(u,v)]^{1/2}$$

称为傅里叶谱或频谱，相角为：

$$\phi(u,v) = \arctan\left[\frac{I(u,v)}{R(u,v)}\right]$$

功率谱定义为：

$$P(u,v) = |F(u,v)|^2 = R^2(u,v) + I^2(u,v)$$

实函数的傅里叶变换是共轭对称的，所以谱是关于原点偶对称的：

$$|F(u,v)| = |F(-u,-v)|$$

相角是关于原点奇对称的：

$$\phi(u,v) = -\phi(-u,-v)$$

有：

$$F(0,0) = \sum_{x = 0}^{M - 1}\sum_{y = 0}^{N - 1}f(x,y)$$

指出零频率项与$f(x,y)$的平均值成正比：

$$F(0,0) = MN\frac{1}{MN}\sum_{x = 0}^{M - 1}\sum_{y = 0}^{N - 1}f(x,y) = MN\overline{f}(x,y)$$

其中$\overline{f}$表示$f$的平均值

• 当二维$DFT$的幅度是一个阵列时，其谱决定了图像中的灰度，相角决定了图像细节，也就是仅用相角信息反傅里叶变换可以重建图像的基本形状，但仅有谱信息不可以

二维卷积定理

$$f(x,y) ★ h(x,y) \Leftrightarrow F(u,v)H(u,v) \\ f(x,y)h(x,y) \Leftrightarrow F(u,v) ★ H(u,v)$$

两个周期函数卷积的时候为了避免缠绕错误必须要对函数进行填充操作，在二维下，令$f(x,y)$和$h(x,y)$分别是大小为$A \times B$和$C \times D$像素的图像阵列，在循环卷积中的缠绕错误可以通过对这两个函数进行零填充来避免：

$$f_p(x,y) = \begin{cases} f(x,y) & 0 \leq x \leq A - 1,0 \leq y \leq B - 1 \\ 0 & A \leq x \leq P,B \leq y \leq Q \end{cases} \\ h_p(x,y) = \begin{cases} h(x,y) & 0 \leq x \leq C - 1,0 \leq y \leq D - 1 \\ 0 & C \leq x \leq P,D \leq y \leq Q \end{cases} \\ P \geq A + C - 1 \\ Q \geq B + D - 1$$

填充后的图像大小为$P \times Q$

频率域滤波基础

假定输入图像为一幅大小为$M \times N$的$f(x,y)$，那么基本滤波公式如下：

$$g(x,y) = \mathfrak{J}^{-1}[H(u,v)F(u,v)]$$

$F(u,v)$是输入图像$f(x,y)$的$DFT$，$H(u,v)$是滤波函数，$g(x,y)$是滤波后的输出图像

频率域滤波步骤小结

• 给定一幅大小为$M \times N$的输入图像$f(x,y)$，选择填充参数$P = 2M,Q = 2N$
• 对$f(x,y)$添加必要数量的0，形成大小为$P \times Q$的填充后的图像$f_p(x,y)$
• 用$(-1)^{x + y}$乘以$f_p(x,y)$移到其变换的中心
• 计算变换后的图像的傅里叶变换得到$F(u,v)$
• 生成一个实的、对称的滤波函数$H(u,v)$，其大小为$P \times Q$，中心在$(P/2,Q/2)$处
• 用阵列相乘形成乘积$G(u,v) = H(u,v)F(u,v)$，即$G(i,k) = H(i,k)F(i,k)$
• 得到处理后的图像：

$$g_p(x,y) = \left\{real\left[\mathfrak{J}^{-1}[G(u,v)]\right]\right\}(-1)^{x+y}$$

其中为了忽略由于计算不准确导致的寄生复分量，选择了实部，下标$p$指出了我们处理的是填充后的阵列

• 通过从$g_p(x,y)$的左上象限提取$M \times N$的区域，得到最终处理结果$g(x,y)$

空间和频率域滤波间的对应

• 傅里叶变换对：空间滤波器和频率域滤波器

$$h(x,y) \Leftrightarrow H(u,v)$$

• 在实践中使用空间滤波器，因为实现时速度快而且容易，但是滤波的概念在频率域更加直观，有些在空间域异常困难的或不可能直接用公式表达的任务在频率域却变得很容易，综合两个域的特性的优点的一个方法是在频率域规定一个滤波器，计算它的$IDFT$得到全尺寸的空间滤波器
• 一个高斯函数的正、反傅里叶变换都是实高斯函数，考虑一维高斯滤波器，用$H(u)$表示一维频率域高斯滤波器：

$$H(u) = Ae^{\frac{-u^2}{2\sigma^2}}$$

空间域中相应滤波器为上式的傅里叶反变换：

$$h(x) = \sqrt{2\pi}\sigma Ae^{-2\pi^2\sigma^2x^2}$$

使用频率滤波器平滑图像

一幅图像的边缘和其他尖锐的灰度转变如噪声对其傅里叶变换的高频内容有贡献，滤波函数$H(u,v)$可以理解为大小为$P \times Q$的离散函数，离散频率变量范围为$u = 0,1,\cdots,P - 1,v = 0,1,2,\cdots,Q - 1$

理想低通滤波器ILPF

$$H(u,v) = \begin{cases} 1 & D(u,v) \leq D_0 \\ 0 & D(u,v) \gt D_0 \end{cases}$$

其中：

$$D(u,v) = \left[(u - P / 2)^2 + (v - Q / 2)^2\right]^{1/2}$$

$D_0$是一个常数，称为截止频率，这一滤波器在半径为$D_0$的圆内所有频率无衰减的通过，而在此圆之外的所有频率则完全的被衰减，这种滤波器不能使用电子元件来实现，但是可以在计算机上面模拟

截止频率和规定的总图像功率值$P_T$的关系，$u = 0,1,\cdots,P - 1,v = 0,1,\cdots,Q - 1$:

$$P_T = \sum_{u = 0}^{P - 1}\sum_{v = 0}^{Q - 1}P(u,v) \\ P(u,v) = |F(u,v)|^2 = R^2(u,v) + I^2(u,v)$$

如果DFT已经被中心化了，那么原点位于频率矩形的中心处，半径为$D_0$的圆将包含$\alpha\%$的功率：

$$\alpha = 100\left[\sum_u\sum_vP(u,v)/P_T\right]$$

振铃现象是理想滤波器的一种特性

布特沃斯低通滤波器BLPF

截止频率位于距原点$D_0$处的$n$阶布特沃斯低通滤波器$BLPF$的传递函数定义为：

$$H(u,v) = \frac{1}{1 + \left[D(u,v)/D_0\right]^{2n}}$$

具有平滑传递函数的滤波器的截止频率的定义：使得$H(u,v)$下降为其最大的某个百分比的点，从上式可以得出截止频率点是$D(u,v) = D_0$，也就是$H(u,v) = 0.5$的点，这种滤波器大大减小了振铃现象，在低阶的时候振铃现象很难察觉，但是随着$n$的增加，振铃现象会逐渐明显，二阶的BLPF是有效的低通滤波和可接受的振铃特性之间的好的折中：

高斯低通滤波器GLPF

$$H(u,v) = e^{-D^2(u,v)/2\sigma^2}$$

$\sigma$是关于中心的扩展度的度量，通过令$\sigma = D_0$：

$$H(u,v) = e^{-D^2(u,v)/2D_0^2}$$

$D_0$是截止频率，此时$D(u,v) = D_0,H(u,v) = 0.607$，GLPF的傅里叶反变换也是高斯的，也就是通过IDFT得到的空间滤波器没有振铃现象：

使用频率域滤波器锐化图像

高通滤波会衰减傅里叶变换中的低频分量而不会扰乱高频信息,滤波函数$H(u,v)$可以理解为大小为$P \times Q$的离散函数，离散频率变量范围为$u = 0,1,\cdots,P - 1,v = 0,1,2,\cdots,Q - 1$，一个高通滤波器可以从给定的低通滤波器从下式得到：

$$H_{HP}(u,v) = 1 - H_{LP}(u,v)$$

理想高通滤波器IHPF

$$H(u,v) = \begin{cases} 0 & D(u,v) \leq D_0 \\ 1 & D(u,v) \gt D_0 \end{cases}$$

IHPF和ILPF是相对的，ILPF将半径$D_0$的圆内的所有频率置零，毫无衰减的通过圆外的所有频率，在物理上也是没有办法实现的，也会产生失真的情况

布特沃斯高通滤波器BHPF

$$H(u,v) = \frac{1}{1 + [D_0/D(u,v)]^{2n}}$$

高斯高通滤波器GHPF

$$H(u,v) = 1 - e^{-D^2(u,v)/2D_0^2}$$

频率域的拉普拉斯算子

$$H(u,v) = -4\pi^2(u^2 + v^2)$$

或者，关于频率矩形的中心，使用如下滤波器：

$$H(u,v) = -4\pi^2\left[(u - P/2)^2 + (v - Q/2)^2\right] = -4\pi^2D^2(u,v)$$

其中$D(u,v)$是距离函数，拉普拉斯图像由下式得到：

$$\nabla^2f(x,y) = \mathfrak{J}^{-1}\left\{H(u,v)F(u,v)\right\}$$

增强可以用下式实现：

$$g(x,y) = f(x,y) + c\nabla^2f(x,y)$$

这里因为$H(u,v)$是负的，所以$c = -1$，在这里$f(x,y)$和$\nabla^2f(x,y)$可能不是在一个量级，处理该问题的一个方法是在计算$f(x,y)$的DFT之前将$f(x,y)$的值归一化到区间$[0,1]$内，并用它的最大值除以$\nabla^2f(x,y)$

钝化模板、高提升滤波和高频强调滤波

频率域中的模板：

$$g_{mask}(x,y) = f(x,y) - f_{LP}(x,y)$$

其中：

$$f_{LP}(x,y) = \mathfrak{J}^{-1}[H_{LP}(u,v)F(u,v)]$$

$F_{LP}(x,y)$是平滑之后的图像，类似于$\overline{f}(x,y)$，有：

$$g(x,y) = f(x,y) + k \times g_{mask}(x,y)$$

$k = 1$时的钝化模板，$k \gt 1$时高提升滤波器，高频强调滤波器：

$$\begin{aligned} g(x,y) &= \mathfrak{J}^{-1}\left\{\left[1 + k \times \left[1 - H_{LP}(u,v)\right]\right]F(u,v)\right\} \\ &= \mathfrak{J}^{-1}\left\{\left[1 + k \times H_{HP}(u,v)\right]F(u,v)\right\} \end{aligned}$$

高通滤波器将直流项设置为0，这样，就把滤波后的图像的平均灰度减小为0，高频强调滤波器不存在这一个问题，因为在高通滤波器上加了1，常数$k$给出了影响最终结果的高频的比例，高频强调滤波的更一般的公式为：

$$g(x,y) = \mathfrak{J}^{-1}\left\{\left[k_1 + k_2 \times H_{HP}(u,v)\right]F(u,v)\right\}$$

同态滤波

一幅图像$f(x,y)$可以表示成照射分量$i(x,y)$和反射分量$r(x,y)$的乘积：

$$f(x,y) = i(x,y)r(x,y)$$

因为：

$$\mathfrak{J}[f(x,y)] \neq \mathfrak{J}[i(x,y)]\mathfrak{J}[r(x,y)]$$

所以假定：

$$z(x,y) = \ln f(x,y) = \ln i(x,y) + \ln r(x,y)$$

有：

$$\mathfrak{J}\{z(x,y)\} = \mathfrak{J}\{\ln f(x,y) \} = \mathfrak{J}\{\ln i(x,y) \} + \mathfrak{J}\{\ln r(x,y) \} \\ Z(u,v) = F_i(u,v) + F_r(u,v)$$

$F_i(u,v),F_r(u,v)$分别是$\ln i(x,y),\ln r(x,y)$的傅里叶变换，可以利用$H(u,v)$对$Z(u,v)$进行滤波：

$$S(u,v) = H(u,v)Z(u,v) = H(u,v)F_i(u,v) + H(u,v)F_r(u,v)$$

在空间域，滤波之后的图像为：

$$s(x,y) = \mathfrak{J}^{-1}\{S(u,v)\} = \mathfrak{J}^{-1}\{H(u,v)F_i(u,v)\} + \mathfrak{J}^{-1}\{H(u,v)F_r(u,v)\}$$

定义：

$$i'(x,y) = \mathfrak{J}^{-1}\{H(u,v)F_i(u,v)\} \\ r'(x,y) = \mathfrak{J}^{-1}\{H(u,v)F_r(u,v)\}$$

所以有：

$$s(x,y) = i'(x,y) + r'(x,y)$$

因为$s(x,y)$是原图像取自然对数产生的，所以输出之后的图像为：

$$g(x,y) = e^{s(x,y)} = e^{i'(x,y)}e^{r'(x,y)} = i_0(x,y)r_0(x,y)$$

其中：

$$i_0(x,y) = e^{i'(x,y)} \\ r_0(x,y) = e^{r'(x,y)}$$

是处理后的图像的照射分量和反射分量，图像的照射分量通常由慢的空间变化来表征，而反射变量引起突变，也就是低频成分与照射相联系，高频成分与反射相联系，使用同态滤波器可以更好的控制照射分量和反射分量，例如下面就是一个同态滤波器：

可以看到上面的滤波器趋向于衰减低频照射的贡献，增强了高频反射的贡献，最终可以达到动态范围压缩和对比度增强的效果，类似于高频强调滤波器

选择性滤波

目的：之前讨论的滤波器在整个频率矩形上面操作，但是在很多应用中，需要处理指定的频段或频率矩形的小区域，第一类滤波器称为带阻滤波器或带通滤波器，第二类称为陷波滤波器

带阻滤波器和带通滤波器

• 带阻滤波器：
  
• 带通滤波器：从带阻滤波器中得到

$$H_{BP}(u,v) = 1 - H_{BR}(u,v)$$

陷波滤波器

陷波滤波器拒绝或通过实现定义的关于频率矩形中心的一个领域的频率，零相移滤波器必须是关于原点对称的，因此一个中心位于$(u_0,v_0)$的陷波在位置$(-u_0,-v_0)$必须有一个相应的陷波，陷波滤波器可以用中心已被平移到陷波滤波器中心的高通滤波器的乘积来构造，一般形式为：

$$H_{NR}(u,v) = \prod_{k = 1}^QH_k(u,v)H_{-k}(u,v)$$

其中$H_k(u,v),H_{-k}(u,v)$是高通滤波器，他们的中心分别位于$(u_k,v_k),(-u_k,-v_k)$处，这些中心是根据频率矩形的中心$(M/2,N/2)$来确定的，对于每个滤波器，距离的计算由下式执行：

$$D_k(u,v) = [(u - M / 2 - u_k)^2 + (v - N / 2 - v_k)^2]^{1/2} \\ D_{-k}(u,v) = [(u - M / 2 + u_k)^2 + (v - N / 2 + v_k)^2]^{1/2}$$

例如一个包含三个陷波对的$n$阶布特沃斯陷波带阻滤波器：

$$H_{NR}(u,v) = \prod_{k = 1}^3\left[\frac{1}{1 + [D_{0k}/D_k(u,v)]^{2n}}\right]\left[\frac{1}{1 + [D_{0k}/D_{-k}(u,v)]^{2n}}\right]$$

陷波带通滤波器可由陷波带通滤波器得到：

$$H_{NP}(u,v) = 1 - H_{NR}(u,v)$$

陷波滤波的一个主要应用就是选择性的修改DFT的局部区域