精读Bi-Exponential Edge-Preserving Smoother

注：本文全称为 Bi-Exponential Edge-Preserving Smoother ，原文链接为：http://bigwww.epfl.ch/publications/thevenaz1202.pdf。

文章脉络：

首先是摘要，接下来第一节综述了以前较为流行的几种边缘保留平滑滤波器，第二节则直截了当地用数学公式描述了Bi-Exponential Edge-Preserving Smoother（下面简称为BEEPS）原理，第三节讨论了 BEEPS 与双边滤波器的理论关系，第四节通过实验结果定性及定量比较了 BEEPS 和双边滤波器，第五节则展示了在何种情况下 BEEPS 能完成双边滤波器完成不了的任务，第六节则是惯常的总结。

摘要部分

文章摘要部分主要概述了 BEEPS 特点，重点在于最后一句：> **The cost per pixel is constant and dependsneither on the data nor on the filter parameters, not even on the degree of smoothing.**图像的计算花费（主要体现在运算速度上）只与图像本身的大小有关，而和图像数据、滤波器参数、最终的平滑度无关，这的确是一个重大的改进，具体怎么实现的我们要往后面看。

Introduction

这一部分的主线内容十分清晰明了；最开始提到图像平滑，人们最常用的是 space-invariant linear smoothers（线性空间不变滤波器），出现的问题是 indiscriminate blurring（不加选择地平滑），于是大家开始研究非线性滤波器，早期这类滤波器的问题是 do not adpat the data（不能有效适应数据本身），直到出现了双指数边缘滤波器（bilateral filter），然后大家的重点就在于如何加速这种滤波器。之前的办法集中在优化 range filter ，这样 spatial filter 则不受约束。不过，本文有一个相反的思路：> **we constrain s to take the form of a bi-exponential filter. This allows for the design of a substitute to the 1-D bilateral filter that runs as a pair of one-tap recursivefilters.**选择约束 spatial filter ，这使得双边滤波器可以优化成一个头尾相连递归的 1-D bilateral filter 。
—附注：双边滤波器的精髓在于它由空域滤波器和值域滤波器两部分组成，而一般的模糊滤波器（如高斯模糊滤波器）只考虑空域，原因是我们通常认为位置相邻的两个像素点灰度值近似相等，但图像边缘等高频信息处却并非如此，所以高斯模糊等一般方法会使得边缘一样变得模糊，起不到保边的效果。—

Proposed algorithm

数学原理

这一部分主要在于从数学公式入手描述 BEEPS的原理。其与双边滤波器类似的结构是两个辅助序列（就像双边滤波器也是由空域、值域两个滤波器组成），其中一个辅助序列如下图所示，是一个向前递进的序列：

而另一个辅助序列则是一个向后递归的序列：

然后将递进序列和回归序列按下式合并成输出序列 y:

其中关于 r ，原文在此的解释是一种 range filter ，而在随后的下文edge preservation 中，对此有具体的描述：

We consider now the general use of the BEEPS, with some nonvanishing degree of smoothing 0 <λ<1 and a nontrivial range filter r. Although r can be chosen freely, it is customary to assume that it takes the shape of a centered bump function. In particular, a prototypical instance is the denormalized Gaussian function

即 r 在原文处选择采用下面这种非正规化的高斯函数：

（σ 是图像平滑系数）
因为是一个递归式，所以必然需要一个初始值，而原文作者做了如下的初始化操作：

并且据原文所说这样初始化有利于提升整体速度和性能。（存疑）

至此，BEEPS 的原理已经讲完了，原文在之后有一个原理推导，但是笔者在看这部分时发现推导过程较短，而且有点强行靠拢结论的感觉（主要思想是用 Z 变换来做），所以在此暂且放下不提，待以后搞明白再写。当然这只是个人感觉，也有可能是鄙人愚钝，想要了解这个算法怎么来的读者也不用担心，后面原文有详细介绍 BEEPS 和双边滤波器联系的部分，可以看作是从双边滤波器出发推导的过程。

Edge Preservation

原文本部分分三段递次向下叙述：1. 从滤波器核函数角度描述 BEEPS 2. 大致讲述了双边滤波器3. 详细讲述 BEEPS 与双边滤波器的联系因为本文是精读，所以我们将主要讲述第一和第三部分，至于双边滤波器的部分，读者可自行查阅补充。### BEEPS这一段先讲述了选择前文讲到的非常规高斯核的考量，其次用文字推导的方式着墨讲述了　BEEPS　的保边原理：当处于非边缘区域时，我们有下面的先验：

即认为相邻像素点的灰度值近似相同。将该式带回到数学原理部分的公式，可以得到下式：

并最终推导出：

（在原文中其实是假设了最后一个式子，倒推回第一个式子，与先验想印证。而最后一个式子其实也体现了在非边缘区域图像平滑的思想。）不过这样的结论明显不适合在边缘处进行迭代，那怎么办呢？

In the presence of an edge between indices (k−1) and k, however, the value of ｅ[k] drops and ϕ[k] is properly reinitialized before the recursion is resumed．

利用 r 核函数对所有参数进行 reinitialize。下面一小节讲述双边滤波器的部分有句话我觉得对 BEEPS 和双边滤波器都适用，顺便也贴上来：

The governing principle is to locally adapt the coefficients of a linear smoothing filter to perform smoothing only over neighboring data of similar photometry (i.e., intensity).

只对相邻数据做平滑处理，理解这一点，双边滤波器和 BEEPS 都不会很难懂。### Link Between the BEEPS and the Bilateral Filter既然要谈 BEEPS 与双边滤波器的联系，那自然要提一下双边滤波器。双边滤波器的离散型输出如下：

其中 r 函数表征图像相邻像素灰度值的相似度，空域滤波器的 s 函数则决定了双边滤波器的空间范围。通常我们假设 ，当 r 较小时会导致几乎不进行平滑，反之则在允许的空间范围内进行平滑。

By design, the denominator of (17) ensures a proper normalization— observe that, through r, this normalization depends on data in a nonlinear fashion and must therefore be computed anew at every coordinate k.

当我们选用合适的 s 函数时，上式的分母会进行适当的归一化，该归一化以非线性的方式由数据映射而来。而下面为了方便理解，我们将在一维讨论，所以先将上式改写：

其中 $\Lambda$ 是 $\lambda$ 的归一化因子：

B₀ 是双边滤波器的归一化因子：

下面就有一个重要的推理过程——当图像进行弱光滑处理时，即 $\lambda$大于等于0，远远小于1，此时 n 大于1的项的值很小（极小项），可视作余项处理：

相应的输出 y₀ 为：

同样的推理方式，BEEPS　的　$\psi$ 和 $\varphi$ 也可改写为：

相应的输出 y 为：

注意上式中的 $\alpha$₀、$\alpha$_0P、$\alpha$_0R以及对应的$\alpha$、$\alpha$_P、$\alpha$_R，我们将它们对应求比例：

当我们设定 r _k,k = 1 时，我们得到结论：

根据这个结论，当 $\lambda$ = 0 时，上面三个极限均等于 1 ，也就是说此条件下　BEEPS　就是双边滤波器。

In doing so, it is crucial to remember that the various coefficients {$\alpha$_P,$\alpha$,$\alpha$_R,$\alpha$_0P,$\alpha$₀,$\alpha$_0R} are not constant-valued but depend on the data.

本节最后一点也提到了此时的参数　$\alpha$₀、$\alpha$_0P、$\alpha$_0R、$\alpha$、$\alpha$_P、$\alpha$_R　均依旧依赖于数据本身。### Adjustment of the Contra-Decay因为 BEEPS 不限制使用固定类型的 range filter ，所以当我们希望 BEEPS 和双边滤波器效果相似时，需要是 BEEPS 和双边滤波器的 range filter 成比例，即类型相似。然而根据 Link Between the BEEPS and the Bilateral Filter 一节的第一个式子，双边滤波器的 normalization 是明确的（该节的第二个式子和第三个式子），但是 BEEPS　的是隐含的（该节中 BEEPS 的输出式），所以我们需要用这个式子来 normalization。

为了使 BEEPS 效果和双边滤波器效果相似，且保持前者的 spatial filter 使用bi-exponential，后者使用 Gaussian，所以 BEEPS 反衰减 $\lambda$ 必须按照双边滤波器中的高斯滤波的标准差调整。

在此多提一句，衰减的来源来自：

Application to Several Dimensions

这一节的重点是提出了一种新的二维运算策略：

先把图像按行展开，和先按列展开计算，再把两方面结果取平均值。优点有以下几点：

our range filter can be chosen and tuned freely, while the weight-adaptive method of [16] restricts itself to the use of curvilinear lengths in a space-intensity manifold

1. 可自由选择调整 range filter。

the recursive feedback offered by (2) and (4) to make the adaptation of the range filter less local than that proposed in [16]

2. 使 range filter 的自适应性比原有的的局部递归方案区域更小。（？？？这里笔者确实不太懂，不过原文意思如果是BEEPS在计算每一个像素点的值只依赖于上一个像素点的值，比正常双边滤波器的卷积运算计算量少的话，倒确实是BEEPS高速的原因）。

we apply the two recursions (1) and (3) in parallel and merge them in (5), which ensures that our filter is anisotropic in one dimension3. 确保滤波器在两个维度上保持各向异性。

3. 最重要的一点是当开发语言允许的情况下，这样的两个分支可以并行计算，互不干扰。

缺点则比较明显：

This strategy comes at the price of doubling our computational cost because we consider two branches in the sketch of Fig. 3, but improves over the early solution for the separability of the bilateral filter proposed in [21] where a single branch was followed.

双倍的计算代价。

EXPERIMENTS

一句话就是　BEEPS　比传统双边滤波器快，效果和双边滤波器一样（其实笔者实际测试，在一些光线明暗变化较大的地方，取较好的参数后，BEEPS效果更好）。------最后，根据传统贴出我的复现代码，此代码用MATLAB编写，且需要用此运行，受限于MATLAB本身的特性，笔者并没有进行并行计算处理，敬请谅解，并希望大家在使用之余，多多指教本文与代码。GITHUB：https://github.com/EatonL/BEEPS