高斯滤波与双边滤波

2022-10-05

字数统计: 9.2k字 | 阅读时长: 33分

前言

学习一下高斯滤波和双边滤波，及其需要的前置知识，仅记录一下个人理解，先后逻辑比较混乱，参考请谨慎。如有错误请评论区指正，感谢。

文中图片多数来自网络，我尽量擦去了水印，是为了提升观感，参考文献在文末有注明。本文的出发点，大概是在尽量用直观、快速的方式了解这些抽象概念吧。

频域（Frequency domain）、时域（Time domain）

时域（Time domain）是描述数学函数或物理信号对时间的关系。例如一个信号的时域波形可以表达信号随着时间的变化。

频域（Frequency domain）是描述信号在频率方面特性时用到的一种坐标系。频域最重要的性质是：它不是真实的，而是一个数学构造。时域是惟一客观存在的域，而频域是一个遵循特定规则的数学范畴，频域也被一些学者称为上帝视角。

以上是对时域和频域下的初步定义，时域是好理解的，但要真正理解频域，还需要若干补充文字。

补充一：频谱（Spectrum）

频域即频率域，在频域图像（频谱）中，自变量是频率（横轴），纵轴是该频率信号的幅度（振幅）。频谱图描述了信号的频率结构及频率与该频率信号幅度的关系。

所以，以时间作为变量所进行的研究就是时域，以频率作为变量所进行的研究就是频域。下图生动的表示了这一点：

fdandtd

时域中的红色波形为三个正弦波相叠加的结果，频域中的三条蓝色竖线分别对应三个正弦波的频率与振幅关系。

补充二：傅立叶分析（Fourier analysis）

傅立叶分析是分析学中逐渐形成的一个重要分支，它研究并扩展傅立叶级数和傅立叶变换的概念，又称调和分析。

Fourier在1807年发表的传记和1822年出版的《热分析理论》一书中指出：无论函数多么复杂，只要它是周期的，并且满足某些适度的条件，都可以表示为不同频率的正弦和（或）余弦函数之和的形式，每个正弦和（或）余弦函数都乘以不同的系数。

我先简单说明一下傅立叶分析究竟用来干什么，仅对于图像方面。假如我现在有一个波形如下：

sin35

假如告诉你该复杂周期信号是由两个简单周期信号组成的（后面会提到这两个概念），其中一个是$\sin (3x)$，要求你把该信号从图中取出，并绘制出剩下的信号波形。显然，这个任务在时域几乎不可能完成。

但是在频域中，无非是几条竖线而已，取出其中一个再容易不过了。这是频域的意义。但是不要忘了，现实世界中只有时域，于是，傅立叶分析作为时域与频域的桥梁出现了，他允许我们从时域跨越到频域完成任务（取出$\sin (3x)$）。顺便一提，这个任务在工程上称为滤波。

傅立叶分析可以分为两块，傅立叶级数（Fourier series）用于处理连续周期函数，傅立叶变换(Fourier transformation)则可以处理连续非周期函数。这里的函数是数学上的说法，在物理上，称为信号。总而言之，傅立叶提出的两个重要观点如下：

周期信号都可以表示为谐波关系的正弦信号的加权和（傅里叶级数）
非周期信号都可用正弦信号的加权积分表示（傅里叶变换）

傅立叶选这个正余弦函数也是有原因的，具体可以看下面对正弦波的补充。其实信号不只能分解为正余弦信号，也可以是什么幂级数展开、泰勒展开之类的。

傅立叶变换（傅立叶级数是其特例）就是通过数学的方法反向分解复杂的信号，使之成为若干简单的正余弦信号，那么问题又来了，这么多简单的信号，每个正弦信号都在时域图上画出来也没有什么意义，那么人类只需要记住每个正弦信号的相位和角频率就可以了，记住后随时随地我们都可以制造这些简单的信号，把它们在时域上相加，就会再现相同的复杂信号。

这时，我们就需要考虑如何记录信号的相位和角频率。如我们所知，频谱只能反映角频率和振幅间的关系，但并未包含$\mathrm{A} \sin (\omega x + \phi)$中的全部信息。因此，还需要相位谱的参与。

补充三：正弦波（Sine wave）

首先，正弦波就是一个圆周运动在一条直线上的投影。

s-w

正弦波是频域中唯一存在的波形，这是频域中最重要的规则，即正弦波是对频域的描述，因为频域中的任何波形都可用正弦波合成。这是正弦波的一个非常重要的性质。然而，它并不是正弦波的独有特性，还有许多其他的波形也有这样的性质。正弦波有四个性质使它可以有效地描述其他任一波形：

（1）频域中的任何波形都可以由正弦波的组合完全且惟一地描述。

（2）任何两个频率不同的正弦波都是正交的。如果将两个正弦波相乘并在整个时间轴上求积分，则积分值为零。这说明可以将不同的频率分量相互分离开。

（3）正弦波有精确的数学定义。

（4）正弦波及其微分值处处存在，没有上下边界。

使用正弦波作为频域中的函数形式有它特别的地方。如果变换到频域并使用正弦波描述，有时会比仅仅在时域中能更快地得到答案。

这里一个重要的意识是：正弦和余弦可以通过相位的转变相互转化。暂时把正余弦信号统称为正弦波。

正弦曲线可表示为：

$\mathrm{y}=\mathrm{A} \sin (\omega x + \phi)+\mathrm{k}$

$k$、$\omega$ 和 $\phi$ 是常数 $(k, \omega, \phi \in R$ 且 $\omega \neq 0)$

A —— 振幅，当物体作轨迹符合正弦曲线的直线往复运动时，其值为行程的 $\frac{1}{2}$ 。
$(\omega x + \phi)$ —— 相位（Phase），反映变量y所处的状态。
$\phi$ —— 初相, $x=0$ 时的相位，反映在坐标系上则为图像的左右移动。
$\mathrm{k}$ —— 偏距, 反映在坐标系上则为图像的上移或下移。
$\omega$ —— 角速度，或角频率，控制正弦周期(单位弧度内震动的次数)。

由于正余弦信号只含有一个频率成分 $\omega$ ，当横轴上的$x$表示时间时，上式不仅是其时域描述也是其频域描述，无需进行变换，可见时域与频域是合二为一的，因此适合将其做为合成其他任意信号的基本信号。

补充四：傅立叶级数与其频谱、相位谱

首先，周期信号可以分为 简单周期信号 和 复杂周期信号，两者的最大区别就是频率结构上分别是单频和多频。对下面的例子来说，正余弦信号是单频信号，而方波信号是多频信号。

将若干个正余弦信号叠加，可得到方波信号，而通过傅立叶级数即可逆向将方波信号分解为多个正余弦信号之和。

f-series

但若要形成标准的90度角矩形波，需要无穷多个正弦波叠加。我们需要接受一个设定：不仅仅是矩形，任何波形都是可以如此方法用正弦波叠加起来获得。

下面再换一个角度看正弦波叠加：

rec-wave

在这几幅图中，最前面黑色的线就是所有正弦波叠加而成的总和，也就是越来越接近矩形波的那个图形。

后面依不同颜色排列而成的正弦波就是组合为矩形波的各个分量，称为频率分量。这些正弦波按照频率从低到高从前向后排列开来，且每一个波的振幅都是不同的。频率越高，波形越密集。

每两个正弦波之间都还有一条直线，那是振幅为0的正弦波。也就是说，为了组成特殊的曲线，有些正弦波成分是不需要的。

如果我们把第一个频率最低的频率分量看作基本单元，该正弦波的角频率为$w_0$，那么频域的基本单元就是$\omega_0$。假设$\omega_0 = 10$，回到频谱中，可以看到该频谱的基本单元正是10。

如果说频域的基本单元是“1”，那么频域的“0”就是一个周期无限长的正弦波（角频率为0），也就是一条直线，0频率也被称为直流分量。

上面已经把频域的事研究的差不多了，但如前文所述，频谱中的信息只是时域中信息的一部分，我们还需要反映相位的相位谱。那么相位谱在哪呢？请看下图，鉴于正弦波是周期的，需要设定一个用来标记正弦波位置的东西，即红色点，红点是距离频率轴最近的波峰。粉色点是红点在下平面的投影，粉点用于标注波峰距离频率轴的距离。

phase-detail

如我们所见，如果说频谱为侧面投影的话，相位谱就是波形在下方的“投影”。但是这个投影不是直接投影，而是间接的。请看下图：

phase-map

说间接投影，就是因为相位不能直接用时间差表示，我们可以看到，时间差是粉点到时域纵轴的距离，我们记为$\varDelta t$，又可以得到当前频率分量的周期$T$，则有相位$p = \frac{\varDelta t}{T} \times 2\pi$，最终相位谱上的纵轴值（相位）是依据时间差和周期算出的。

最终的图如下所示，其实，频域图像也可以叫做幅度频谱，对应的，相位谱叫做相位频谱，因为他们的定义域都是频域，只是值域不同。

tfp

补充五：傅立叶变换与其频谱

前面我们说到，傅立叶级数是傅立叶变换的一个特例，也就是说，傅立叶变换推广了傅立叶级数，这个“推广“体现在傅立叶变换可以处理时域上的非周期的连续信号，将其转换为一个在频域上的非周期连续信号。如下图：

f-t

图片右上角的Frequency resolution是频率分辨率，即数字信号处理系统将相距最近的两个频率分量区分开的能力，它与周期成反比，当一个连续信号是非周期信号时，可以将其周期看作无限大，则该信号近似看为周期无限大的周期连续信号。这代表非周期信号可以通过傅立叶变换分解为无限多个频率无限接近的正余弦信号的和。

试想一下，上面离散的频率分量逐渐靠近…直至连在一起，无穷无尽的分量铺满了空间。

变成了波涛汹涌的大海，这就是连续的频谱下发生的事情。那么，计算时离散频谱下分量的累加，自然也就变成了连续频谱上的积分。

至此，对于研究图像处理来说，似乎已经足够，数学上的东西，需要时再深究也不迟。

补充六：虚数（Imaginary number）和欧拉公式（Euler’s formula）

我们知道，在实数轴上，乘$-1$其实就是使线段逆时针旋转了$180$度，那么乘一个$i$呢？显然，旋转$90$度。

i-r

同时，我们获得了一个垂直的虚数轴。实数轴与虚数轴共同构成了一个复数的平面，也称复平面。这样我们就了解到，乘虚数i的一个功能——旋转。

欧拉公式如下（关于欧拉公式为何被称为“上帝创造的公式”此处不再赘述）：

$e^{i x}=\cos x+i \sin x$

当$x = \pi$时，就有$e^{i x} + 1 = 0$。

这个公式的关键作用之一，是将正弦波统一成了简单的指数形式。

euler

欧拉公式所描绘的，是一个随着时间变化，在复平面上做圆周运动的点，随着时间的改变，在时间轴上就成了一条螺旋线。如果只看它的实数部分，也就是螺旋线在左侧的投影，就是一个最基础的余弦函数。而右侧的投影则是一个正弦函数。

欧拉公式的另一种形式是：

$\sin x=\frac{e^{i x}-e^{-i x}}{2 i}, \cos x=\frac{e^{i x}+e^{-i x}}{2}$

补充七：数学上的傅立叶级数和傅立叶变换

傅立叶级数的三角函数形式如下：

$f(t)=a_0+\sum_{n=1}^{\infty}\left[a_n \cos \left(n \omega t\right)+b_n \sin \left(n \omega t\right)\right], \\ 其中 \omega = \frac{2 \pi}{T}为基波频率, a_n和b_n为傅立叶系数，分别代表余弦分量振幅和正弦分量振幅.$

傅立叶变换相关的，先放放吧。可以参考这三篇。

从信号的角度看图像

图像是信号，这是直觉上极难理解的。

首先，我们需要摆脱信号和时间“绑定”的固有直觉，信号只是由自变量和因变量构成的，自变量不一定是时间。

实际上，图像是一个沿空间分布的信号，他的定义域是图像的$x$轴或$y$轴。图像既然是沿着空间分布的信号，说明它是一个二维信号。如果一维数字信号可以看作是一组数字序列，那么二维数字信号则可以看作是一组数字阵列。图像中沿着水平方向的任何一行可以看作是像素点随着$x$轴变化的信号。而沿着竖直方向的任何一列则可以看作是像素点随着$y$轴变化的信号。单独提取一行或者是一列数据出来，按照像素值的大小随轴的变化画出来，就可以得到我们熟悉的信号波形，如下。

pic-signal

可以明显的看出，图片右方的波形代表黑色横线，因为除了塔尖以外像素变化不大。而下方波形则对应竖线，由上到下。

图像既然是信号，那么自然有频率、幅度和相位的特征。

图像的频率对应到的是图像细节的多少，比如城市的部分频率高于天空的部分。再比如在灰度图中，图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。

图像的幅度对应到的是图像中像素的大小，比如天空中太阳照亮的区域幅度高于其他区域。

假设有图片A和图片B，分别获取A、B的幅度频谱和相位频谱，然后将图片A的幅度谱与图片B的相位谱结合做逆傅立叶变换，最终结果的视觉效果会是图片B加上一些噪声，通过对比实验可以得到一个结论：图像的相位谱中，保留了图像的边缘以及整体结构的信息。 而错误的幅度谱看起来则像是噪声覆盖在原图上，但并没有影响图像的内容本身。

同时，和我们在一维信号中看到的一样，如果想要通过修改频域中相位或振幅的方式来间接修改原图像，则需要同时保留幅度频谱和相位频谱，来通过逆傅立叶变换得到修改后的图像。

空域（Spatial domain）

空间域也叫空域，即所说的像素域，在空域的处理就是在像素级的处理，即直接对图像上的像素值进行增加或减少。对应的另一种操作是在频域上的，空域通过傅立叶变换后，得到的是图像的频谱。表示图像的能量梯度。

图像变换技术

为了有效和快速地对图像进行处理和分析，需要将原定义在图像空间的图像以某种形式转换到另外的空间，利用空间的特有性质方便地进行一定的加工，最后再转换回图像空间以得到所需的效果。如：空域与频域间的相互转化。

信号分析在图像处理中的应用

很多图像处理事实上就是信号处理。比如图像去噪 (Denoise)，一般就是应用数字低通滤波器对图像进行滤波。图像增强 (Detail enhancement)，一般就是应用数字高通滤波器得到图像的高频信号，并对高频信号进行增强。对比度增强 (Contrast enhancement)，一般就是参考画面的亮暗程度 (图像的幅度)，并人为修改亮暗的一种处理。相位的概念一般会在图像的缩放 (Scaling) 中使用到。

图像的灰度、亮度、强度

图像灰度(Image grayscale)：把白色与黑色之间按对数关系分为若干等级，称为灰度，灰度分为256阶，从0到255。用灰度表示的图像称为灰度图，实际上，灰度表征的是单色的亮暗程度。

图象亮度(Image brightness)：指画面的明亮程度，单位是堪德拉每平米(cd/m2) 或称 nits。图象亮度是从白色表面到黑色表面的感觉连续体，由反射系数决定，亮度侧重物体，重在“反射”。在灰度图像中，亮度等于灰度，图像运算处理方式相同。但是在彩色图像中，亮度和对比度相关，即通过对RGB颜色分量的增加（增加亮度）或减少（减少亮度）相同的增量来显示，亮度的调整就是给每个分量乘以一个百分比值。

图像强度(Image intensity)：表示单通道图像像素的强度（值的大小）。在灰度图像中，它是图像的灰度。在RGB颜色空间中，可以理解把它为是R通道的像素灰度值，G通道的像素灰度值，或是B通道的像素灰度值，也就是RGB中含三个图像强度。在其他颜色空间类似，也就是每个通道的图像的像素灰度值。

这里多说一下，为什么RGB也有灰度？实际上，RGB各单通道的图像也都是灰度图，只是最终成像时将他们分别放置在红、绿、蓝的单色通道中再叠加罢了，也就是说，一个RGB彩色图像，是由三个表示对应位置亮与暗程度的图像，通过放置在三个单色通道中再叠加形成的（这里属于理糙话不糙了）。

认识到所有图片都是灰度图或其叠加后，很多图像处理方法的应用都可以举一反三了。

图像滤波

图像滤波，就是对图像的频率进行过滤。本文关注两类，即高通滤波和低通滤波。这两种滤波方式不难理解，从名字就可以看出，高通滤波就是减弱或阻隔低频信号，保留高频信号。而低通滤波则是减弱或阻隔高频信号，保留低频信号。

现在关键的问题是，图像中的高频信号和低频信号分别代表什么？结合前面所提到的，高（低）频信号也叫做高（低）频分量，图像的频率是灰度值变化剧烈程度的标准，可以推出高频信号就是灰度值变化剧烈的地方，同理，低频信号是图像灰度值变化平缓的地方。

具体些说，高频信号就是相邻区域之间灰度相差很大的地方，例如一个影像与背景的边缘部位，通常会有明显的差别，灰度变化很快，也就是变化频率高的部位。同样的，图像的细节处也是属于灰度值急剧变化的区域。另外，噪声（即噪点）也是这样，噪点之所以是噪点，就是因为它与附近的像素点灰度不一样了。

相对应的，低频信号就是那些连续渐变的区域。人眼对图像中的高频信号更为敏感，举例来说，在一张白纸上有一行字，那么人眼会直接聚焦在文字上，而不会太在意白纸本身，这里文字就是高频信号，而白纸就是低频信号。

频谱中心化

研究数字图像有时需要变换到频域做处理，比如滤波等。但直接对数字图像进行二维DFT（离散傅立叶变换）得到的频谱图是高频在中间，低频在四角。频谱图比较亮的地方就是低频，因为图像的能量一般都是集中在低频部分。为了把能量集中起来便于使用滤波器，可以利用二维DFT的平移性质对频谱进行中心化。如下：

s-c

经二维傅立叶变换的平移性质推导，结论是频谱中心化只需对数字图像的每个像素点的取值直接乘以$(-1)^{x+y}$就可以了，其中$x，y$为像素坐标。

频域滤波

频域滤波的基本过程如下：
1、对原图像做频谱中心化处理，
2、对第一步的结果进行DFT，
3、使用某滤波器乘第二步的结果，
4、对第三步结果做反DFT，再做一次频谱中心化，得到滤波后图像。

空域滤波

空域滤波是指利用像素及像素邻域组成的空间进行图像增强的方法。之所以用“滤波”这个词，是因为借助了频域里的概念。并且，频率域的卷积与空间域的乘积存在对应关系，由卷积定理可知所有频域的滤波理论上都可以转化为空域的卷积操作。

空域滤波是在图像空间通过邻域操作完成的，邻域操作常借助模板运算来实现。由此也可以看出，空域滤波的算法比较简单，处理速度快。而频域滤波算法复杂，计算慢。（这里只是相比于做DFT来说，现在好像还有FFT可以加速运算，具体如何先不深究。）

图像的空域滤波还可以按照运算的方式分为线性滤波和非线性滤波，如果运算只是加权求和之类的操作，则为线性滤波。而若是在模版内进行求最值、绝对值等操作，则属于非线性滤波。

模版运算

模板也称为核（kernel）。模板运算的基本思路是将某个像素的值替换为它本身灰度值和其相邻像素灰度值的函数值。

模板一般是$n \times n$的方阵（$n$通常是3、5、7、9等很小的奇数）。当$n$为奇数时，可以定义模板的半径$r = \frac{(n-1)}{2}$。模板中有一个锚点（anchor point），通常是矩阵中心点，和原图像中待计算点对应。整个模板对应的区域，就是原图像中像素点的相邻区域。一个$n \times n$的模板最多可有$n \times n$个系数，该模板的功能由这些系数的取值所决定。

以模版运算中最常用的模版卷积举例，其在空域实现的主要步骤如下：
1、将模板在图中漫游，并将模板锚点与图中某个像素位置重合（待计算点），
2、将模板上的各个系数与模板下各对应像素的灰度值相乘，
3、将所有乘积相加（为保持灰度范围，常将结果再除以模板的系数个数），
4、将上述运算结果（模板的输出响应）赋给待计算点的像素。

高斯滤波（Gaussian filter）

高斯滤波是一种根据高斯函数的形状来选择权值的线性平滑滤波，适用于消除高斯噪声，广泛应用于图像处理的减噪过程。通俗的讲，高斯滤波就是对整幅图像进行加权平均的过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到。

高斯噪声就是概率密度函数服从高斯分布的一类噪声。如果一个噪声，它的幅度分布服从高斯分布，而它的功率谱密度（频谱的绝对值平方？）又是均匀分布的，则称它为高斯白噪声。高斯白噪声的二阶矩不相关，一阶矩为常数，是指先后信号在时间上的相关性。

高斯滤波的空域实现步骤很简单，先获得高斯核，然后对图像做模版卷积即可（可以做padding）。关键还是在高斯核上，或者说，在于高斯核上的权值如何计算上。

当然了，高斯滤波是指用高斯函数作为滤波函数，具体上实现的是模糊还是锐化效果，要根据是高通还是低通来判别。高斯锐化的高斯核就是$1 - 高斯模糊卷积核$。

高斯核

一维和二维的高斯函数如下，标准差$\sigma$的大小决定了高斯函数的宽度：

gaussian

根据二维高斯函数图像来理解高斯滤波就非常的直观了，假设我们要实现高斯模糊（低通），在确定高斯核时，首先需要将锚点（核的中心点）作为高斯图像中的原点（山峰的中心点），然后确定一个标准差$\sigma$和核的半径$r$，根据$r$我们可以计算出高斯核中各点的坐标，通过坐标$(x,y)$和$\sigma$依次计算二维高斯函数值，并填入对应的位置，然后，还需要进行归一化处理，即将核中各值除去权重总和，使核的权重综合为$1$，这样就得到了最终的高斯核。

现在的问题是，如何选择$\sigma$和$r$呢？

对于标准差，可以从二维高斯函数图像来分析，当标准差越小时，图像中的“山峰”越“瘦高”，权重分布越向中心集中，平滑效果越差；当标准差越大时，“山峰”越“矮胖”，权重分布越均匀，则平滑效果越明显。（有理论说空域上的标准差和频域的标准差成负相关，以此可以进一步解释高斯核的高低通效果，但是笔者目前还没看明白，留待日后补充。）

半径的选择同样可以从图像来分析。理论上，高斯分布在定义域的所有地方都有非负值，这就需要一个无限大的核。然而钟形曲线在区间$(\mu - \sigma, \mu + \sigma)$范围内的面积占曲线下总面积的$68\%$，在区间$(\mu - 2\sigma, \mu + 2\sigma)$范围内的面积占$95\%$，在区间$(\mu - 3\sigma, \mu + 3\sigma)$范围内的面积占$99.7\%$，通常情况下$3\sigma$以外的区域所占面积已经很小，可以被忽略（认为该段分布近似包含所有情况）。所以，当半径取$3\sigma$时，就基本上可以近似使用了，此时核的大小为$(6\sigma + 1) \times (6\sigma + 1)$。

最后，对高斯核进行归一化处理的一个重要目的是确保卷积运算后，像素值处于$[0, 255]$的范围之内。

实验-高斯模糊

光说不练，不是好文章，下面就实验一下实现高斯滤波。

import cv2
import numpy as np


def gaussian_filter(img, sigma=0.5):
    if len(img.shape) == 3:
        H, W, C = img.shape
    else:
        img = np.expand_dims(img, axis=-1)
        H, W, C = img.shape

    # 计算半径，取4倍
    radius = int(sigma * 4)
    # 计算高斯核大小
    n = 2 * radius + 1
    # padding大小
    pad = n // 2
    # padding的地方置为0
    out = np.zeros((H + pad * 2, W + pad * 2, C), dtype=np.float64)
    out[pad: pad + H, pad: pad + W] = img.copy().astype(np.float64)
    # 构造高斯核坐标
    x_index, y_index = np.mgrid[-radius:radius + 1, -radius:radius + 1]
    # 计算核值
    k = np.exp(-(x_index ** 2 + y_index ** 2) / (2 * sigma ** 2))
    # 归一化
    k = k / k.sum()
    # 原图像深拷贝
    origin = out.copy()
    # 高斯核卷积，计算有点慢
    for y in range(H):
        for x in range(W):
            for c in range(C):
                out[pad + y, pad + x, c] = np.sum(k * origin[y: y + n, x: x + n, c])
    # 切除填充部分
    out = out[pad: pad + H, pad: pad + W].astype(np.uint8)
    return out


if __name__ == '__main__':
    # 读取图像
    img = cv2.imread("./data/bridge.jpeg")
    # 滤波
    result = gaussian_filter(img, sigma=2.5)
    # 保存结果
    cv2.imwrite("./data/bridge-blur.jpeg", result)
    # 显示
    cv2.imshow("bridge", img)
    cv2.imshow("bridge-blur", result)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

输入（笔者在重庆随手拍的一个桥，叫什么忘记了）：
bridge
输出：
blur-bridge

双边滤波（Bilateral filter）

开局先甩一个定义：

双边滤波是一种非线性的保边滤波器，通常在计算机视觉中用作工作流中的简单降噪阶段。它将每个输出像素的强度计算为输入图像中附近像素的强度值的加权平均值。至关重要的是，权重不仅取决于当前像素和相邻像素之间的欧几里得距离，还取决于它们之间的辐射差异（例如颜色强度差异）。结果是保留了边缘，同时平滑了具有相似强度的区域。

只看定义，好像懂了，又好像没懂，那么下面就来详细的解释一下。众所周知，图像的噪声和边缘都属于高频分量，高斯滤波在降低噪声的同时没有很好的保留边缘，而是一起“降”了，这对于人眼观看（或者特征提取）就不怎么友好了，形成的图片模糊成了一团。

这时，有人想了，能不能在降低噪声的同时，保持图像的边缘呢？这就是双边滤波干的事。双边滤波同样是基于模版运算的，他的特点在于模版的权值同时考虑了空域上距离的差别和灰度值上数值的差别。这里后者也可以叫做灰度距离，与前者的欧式距离相呼应（没错前者实际上就是欧式距离）。

回顾高斯滤波中的高斯核，它的权值分配是按照二维高斯图像的形状分配的，这实际上是一种考虑空域上距离差别的方式，远离中心点的像素点会被分配到更少的权值。

那么为什么要考虑灰度值的差别呢？这是由于边缘的表现其实是灰度值的突变，如果想要保留边缘，则不能将高灰度值的区域与低灰度值的区域相“混合”。传统的高斯核无法考虑当前覆盖区域的灰度值信息。

由此，双边滤波的核函数是空域核与值域（灰度值域）核的综合：在图像的平坦区域，也就是灰度值变化小的区域，值域核权重趋于$1$，此时空域核权重起主要作用，与高斯模糊类似；在图像的含边缘区域，灰度值变化很大，则值域核权重会增大，从而起到保留边缘的效果。如下图：

b-f

图片左侧输入表示输入图像的灰度值形状，可以明显看到图像左侧整体灰度值低而右侧整体灰度值高，值变化很大，图像中部存在一个断层，这表示该区域属于含边缘区域。而断层两侧的部分，也有起起伏伏的小土包，这表示灰度值变化很小的区域，即平坦区域。当核移动到图中红色箭头所指位置的正上方时，空域核权重与值域核权重的情况如方框中所示，显然，该位置属于含边缘区域，值域核在形状上也表征为断层。它与空域核结合后的结果就是最终的核形状，即一个被切了一半的高斯图像，这个核的权重也很容易想象，整体上，左侧权重低而右侧权重高，起到了保留边缘的作用。

以上，应该能够理解双边滤波是如何做的了，那具体怎么做呢？接着往下看。

双边滤波器中，最终输出像素的值$g$同样依赖于领域像素值的加权组合：

$g(i, j)=\frac{\sum_{k, l} f(k, l) w(i, j, k, l)}{\sum_{k, l} w(i, j, k, l)}$

这里，$(i,j)$是待计算位置的坐标，$(k,l)$是领域像素位置坐标，$函数f$表示计算灰度。核权重$w$就是重点了，它取决于空域核权重$s$与值域核权重$r$：

$s(i, j, k, l)=\exp \left(-\frac{(i-k)^2+(j-l)^2}{2 \sigma_s^2}\right), \\ r(i, j, k, l)=\exp \left(- \frac{(f(i,j)-f(k,l))^2}{2 \sigma_r^2}\right)$

$\sigma_s$ 和 $\sigma_r$ 分别是空域高斯（Spatial Gaussian）和值域高斯（Range Gaussian）的标准差，代表了对图片的滤波程度。$w$取决于$s$和$r$的乘积：

$w(i, j, k, l) = \exp \left(-\frac{(i-k)^2+(j-l)^2}{2 \sigma_s^2} - \frac{(f(i,j)-f(k,l))^2}{2 \sigma_r^2}\right)$

根据该权重函数，就可以获得双边滤波的核了，然后再进行模版卷积，归一化处理即可。

正式的一个公式应该是：

$B F[I]_{\mathrm{p}}=\frac{1}{W_{\mathrm{p}}} \sum_{\mathbf{q} \in \mathcal{S}} G_{\sigma_s}(\|\mathbf{p}-\mathbf{q}\|) G_{\sigma_r}\left(\lvert I_{\mathbf{p}}-I_{\mathbf{q}} \rvert \right) I_{\mathbf{q}}$

其中$W_{\mathrm{p}}$是一个归一化因子（normalization factor）：

$W_{\mathrm{p}} = \sum_{\mathbf{q} \in \mathcal{S}} G_{\sigma_s}(\|\mathbf{p}-\mathbf{q}\|) G_{\sigma_r}\left(\lvert I_{\mathbf{p}}-I_{\mathbf{q}} \rvert \right)$

这里，$B F[I]_{\mathrm{p}}$为像素$p$滤波后的灰度值；$I_{\mathbf{p}}$为像素$p$原本的灰度值；$\mathrm{p}$则表示像素$p$的空间坐标；$G_{\sigma_s}$和$G_{\sigma_r}$分别为空域高斯和值域高斯；$|\mathbf{p}-\mathbf{q}|$表示取模，即欧式距离；$\lvert I_{\mathbf{p}}-I_{\mathbf{q}} \rvert$为灰度值差的绝对值；像素$q$为领域点；$\mathcal{S}$表示整个核区域内的点。

这样解读完以后，其实和一开始的公式是一样的。

参考文献

前言