论文阅读：ROBUST KERNEL REGRESSION FOR RESTORATION AND RECONSTRUCTION OF IMAGES FROM SPARSE NOISY DATA

核回归与稀疏、噪声图像数据的重建

ROBUST KERNEL REGRESSION FOR RESTORATION AND RECONSTRUCTION OF IMAGES
FROM SPARSE NOISY DATA

Abstract

解决的问题：

从噪声损坏或稀疏收集的样本中重建图像与信号

解决方案

文章采用了基于非参数估计方法（non-parametric estimation methods）的局部自适应核，它既考虑了可用样本的局部密度，又考虑了这些样本的实际值，以实现对样本“几何形状”和“辐射度”的自适应。方法不依赖于有关噪声或采样分布的特定假设。

效果

• 适用于多种问题，包括图像放大，仅从15％（不规则采样）像素中重建高质量图像，
• 从噪声和不确定较多的数据集中实现超分辨率
• 图像高斯噪声以及其他噪声的去噪，有效去除压缩伪像

key words

Inverse problem, image reconstruction, piecewise
polynomial approximation, nonlinear estimation

1.Introduction

• 经典的参数型图形处理方法依赖于特定的信号处理模型，但非参数方法依赖于数据本身来控制模型的结构，这种模型被称为回归函数。
• 本文使用了名为“核回归”的非参数方法，这种方法归纳了归一化卷积normalized convolution，双边滤波器bilateral
  filter和移动最小二乘法moving least-squares等方法。
• 该方法是通用框架，可在单帧降噪与多帧超分辨率中应用；对错误数据建模有更好的鲁棒性

2.数据自适应核回归

二维估计问题中,在$x_i$点估计$y_i$的值，只需要通过下面的实在就可以完成
$$y_i=z\left({\mathbf{x}}_i\right)+{\epsilon }_i,i=1,\cdots ,P\text{\hspace{1em}(1)}$$
这里$z(x_i)$就是回归函数。${\epsilon }_i$是独立噪声分布量。假设该函数平滑，$x$是样本$x_i$附近的点，那么$z(x_i)$可以在$x$处泰勒展开：
z ( x i ) ≈ z ( x ) + { ∇ z ( x ) } T ( x i − x ) + 1 2 ( x i − x ) T { H z ( x ) } ( x i − x ) + ⋯ = β 0 + β 1 T ( x i − x ) + β 2 T vech ⁡ { ( x i − x ) ( x i − x ) T } + ⋯ (2) \begin{aligned} z\left(\mathrm{x}_{i}\right) \approx & z(\mathrm{x})+\{\nabla z(\mathrm{x})\}^{T}\left(\mathrm{x}_{i}-\mathrm{x}\right) \\ &+\frac{1}{2}\left(\mathrm{x}_{i}-\mathrm{x}\right)^{T}\{\mathcal{H} z(\mathrm{x})\}\left(\mathrm{x}_{i}-\mathrm{x}\right)+\cdots \\=& \beta_{0}+\beta_{1}^{T}\left(\mathrm{x}_{i}-\mathrm{x}\right) \\ &+\beta_{2}^{T} \operatorname{vech}\left\{\left(\mathrm{x}_{i}-\mathrm{x}\right)\left(\mathrm{x}_{i}-\mathrm{x}\right)^{T}\right\}+\cdots \end{aligned} \tag2 z(xi​)≈=​z(x)+{∇z(x)}T(xi​−x)+21​(xi​−x)T{Hz(x)}(xi​−x)+⋯β0​+β1T​(xi​−x)+β2T​vech{(xi​−x)(xi​−x)T}+⋯​(2)
其中$\nabla z(x) $是$z(x)$在$x$处的梯度，$\mathcal{H}$是Hessian算子。vech是半向量化运算符，按字典顺序将矩阵的“下三角”部分排序为列向量。
$$G\left(x_0\right)={\left[\begin{array}{cccc}\frac{{\partial }^{2}f}{\partial x_1^2}& \frac{{\partial }^{2}f}{\partial x_1\partial x_2}& \cdots & \frac{{\partial }^{2}f}{\partial x_1\partial x_n}\\ \frac{{\partial }^{2}f}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f}{\partial x_2^2}& \cdots & \frac{{\partial }^{2}f}{\partial x_2\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}f}{\partial x_n\partial x_1}& \frac{{\partial }^{2}f}{\partial x_n\partial x_2}& \cdots & \frac{{\partial }^{2}f}{\partial x_n^2}\end{array}\right]}_{x_0}\text{\hspace{1em}(3)}$$
(2)式是Hessian算子的一种形式。实际上表示的就是多元函数的各个方向高次偏导数的矩阵。

{ β n } n = 1 N \left\{\boldsymbol{\beta}_{n}\right\}_{n=1}^{N} {βn​}n=1N​将提供n阶导数的局部信息，(1)式将问题转换成了根据数据估计系数 { β n } n = 0 N \left\{\boldsymbol{\beta}_{n}\right\}_{n=0}^{N} {βn​}n=0N​的值，这也对数据进行几何加权（近的样本权重大）

但是辐射加权（基于局部边界的相对位置）也同样合适。提出这个问题是为了解决以下优化问题：
min ⁡ { β n } n = 0 N ∑ i = 1 P ∣ y i − β 0 − β 1 T ( x i − x ) − β 2 T vech ⁡ { ( x i − x ) ( x i − x ) T } − ⋯   ∣ m K ( x i − x , y i − y ) (4) \begin{array}{c}{\min _{\left\{\boldsymbol{\beta}_{n}\right\}_{n=0}^{N}} \sum_{i=1}^{P} | y_{i}-\beta_{0}-\boldsymbol{\beta}_{1}^{T}\left(\mathbf{x}_{i}-\mathbf{x}\right)-} \\ {\boldsymbol{\beta}_{2}^{T} \operatorname{vech}\left\{\left(\mathbf{x}_{i}-\mathbf{x}\right)\left(\mathbf{x}_{i}-\mathbf{x}\right)^{T}\right\}-\left.\cdots\right|^{m} \mathrm{K}\left(\mathbf{x}_{i}-\mathbf{x}, y_{i}-y\right)}\end{array} \tag4 min{βn​}n=0N​​∑i=1P​∣yi​−β0​−β1T​(xi​−x)−β2T​vech{(xi​−x)(xi​−x)T}−⋯∣mK(xi​−x,yi​−y)​(4)
$K(x，y)$就是核函数，作为几何距离和辐射距离的惩罚函数。$m$是惩罚参数，常取2，这样，(3)式就成了加权最小二乘问题。m 取1时可以显着提高针对异常值的鲁棒性。该值实际上在核回归框架中合并了鲁棒的$l1$范数估计量。

在矩阵形式下(4)中的问题转化为加权范数下的求解：
$$\hat{\mathbf{b}}=\arg \underset{\mathbf{b}}{\min }{\Vert \mathbf{y}-{\mathbf{X}}_{\mathbf{x}}\mathbf{b}\Vert }_{{\mathbf{W}}_{\mathbf{x}}}^m$$
其中
$$\mathbf{y}={\left[y_1,y_2,\cdots ,y_p\right]}^T,\mathbf{b}={\left[{\beta }_0,{\mathbf{\beta }}_1^T,\cdots ,{\mathbf{\beta }}_N^T\right]}^T$$

$${\mathbf{W}}_{\mathbf{x}}=\mathrm{diag}\left[\mathrm{K}\left({\mathrm{x}}_1-\mathrm{x},y_1-y\right),\cdots ,\mathrm{K}\left({\mathrm{x}}_p-\mathrm{x},y_p-y\right)\right]$$

$${\mathbf{X}}_{\mathbf{x}}=\left[\begin{array}{cccc}1& {\left({\mathbf{x}}_1-\mathbf{x}\right)}^T{\mathrm{vech}}^T\left\{\left({\mathbf{x}}_1-\mathbf{x}\right){\left({\mathbf{x}}_1-\mathbf{x}\right)}^T\right\}& \cdots & \\ 1& {\left({\mathbf{x}}_2-\mathbf{x}\right)}^T{\mathrm{vech}}^T\left\{\left({\mathbf{x}}_2-\mathbf{x}\right){\left({\mathbf{x}}_2-\mathbf{x}\right)}^T\right\}& \cdots & \\ ⋮& ⋮& ⋮& ⋮\\ 1& {\left({\mathbf{x}}_p-\mathbf{x}\right)}^T{\mathrm{vech}}^T\left\{\left({\mathbf{x}}_p-\mathbf{x}\right){\left({\mathbf{x}}_p-\mathbf{x}\right)}^T\right\}& \cdots & \end{array}\right]$$

使用最速下降法求解这个问题：
$${\hat{\mathbf{b}}}^{(n+1)}={\hat{\mathbf{b}}}^{(n)}+\alpha {\mathbf{X}}_{\mathbf{x}}^T{\mathbf{W}}_{\mathbf{x}}\mathrm{sign}\left(\mathbf{y}-{\mathbf{X}}_{\mathbf{x}}{\hat{\mathbf{b}}}^{(n)}\right)\odot {\left|\mathbf{y}-{\mathbf{X}}_{\mathbf{x}}{\hat{\mathbf{b}}}^{(n)}\right|}^{m-1}$$
其中$\alpha$是定义梯度方向上步长的标量，而$\odot$是逐元素乘法运算符。

回归的阶数（N）影响局部逼近的复杂度。线性和二次近似（分别对应于N = 0、1、2）比较常用。特别地，选择m = 2的局部常数估计，可获得局部线性自适应滤波器，这就是Nadaraya-Watson估计器（NWE）。
通常，较低的近似值（例如NWE）因为自由度较小会导致图像更平滑（较大的偏差和较小的方差）。高阶会导致过拟合和较小的偏差、较大的估计方差。在第4节的实验中用了二阶（N = 2）近似值。

3.核函数的选择

思路：

非自适应的核函数–>推广导出两种自适应核函数

1.经典核函数

方法原理

$$\mathrm{K}\left({\mathrm{x}}_i-\mathrm{x},y_i-y\right)\equiv K_{{\mathrm{H}}_i}\left({\mathrm{x}}_i-\mathrm{x}\right)$$

这种核函数只关注采样点之间的距离，距离越远值越小其中$K_{{\mathrm{H}}_i}$的定义如下：

$$K_{{\mathbf{H}}_i}(\mathbf{t})=\frac{1}{\det \left({\mathbf{H}}_i\right)}K\left({\mathbf{H}}_i^{-1}\mathbf{t}\right)$$

$H_i$也被称作“距离平滑矩阵”，尺寸为$2\times 2$。一般取
$${\mathbf{H}}_i=h{\mu }_i{\mathbf{I}}_2$$
其中${\mu }_i$是捕获数据样本局部密度的标量，$h$是全局平滑参数，拓展核函数包含足够的样本。

一般情况下，小尺寸核可用于有较多有效信息的样本，大尺寸核适用于稀疏样本。

方法局限性

• 经典的回归算法在图像平滑区域可以获得近似最佳的滤波效果，但对图像边缘附近的像素估计时，容易造成边缘模糊。这是由于，在图像边缘附近存在很多跳变点，因此核窗口中包含不连续的样本点，使得核估计的结果存在较大偏差。

2.双边核函数

考虑经典核回归估计, 是数据的局部线性组合,尽管具有很好的性质，且相对容易理解，但因为数据的局部线性条件而存在固有限制，为了更有效的适用于各种数据，使模型具有非线性，可以考虑在经典方法以像素位置决定权值的基础上，引入像素灰度值。即核函数在计算权值时考虑两个因素：空间距离和灰度距离，称为自适应控制核回归。

方法原理

$$\mathrm{K}\left({\mathrm{x}}_i-\mathrm{x},y_i-y\right)\equiv K_{{\mathrm{H}}_i}\left({\mathrm{x}}_i-\mathrm{x}\right)K_{h_r}\left(y_i-y\right)\text{\hspace{1em}(11)}$$

其中自适应核函数由两部分组成：

• $K_{{\mathrm{H}}_i}$$(x_i-x)$采样点的空间距离，与上一节中所述一致
• $K_{h_r}\left(y_i-y\right)$灰度距离（辐射距离）

其中$h_r$被称作“辐射平滑矩阵”，

自适应核回归方法不仅依赖于采样位置和密度，也依赖于这些采样点的灰度。因此，回归核的有效尺寸和形状是与图像局部特征相关的。下面两张图比较了经典核和自适应核回归边缘处的不同分布情况。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-8v4pLR7f-1582978320334)(C:\Users\sgb6372\AppData\Roaming\Typora\typora-user-images\image-20191128112214485.png)]

方法局限性

• 双边核的应用仅限于降噪问题，因为在仲裁位置$x$上的像素值$y$可能无法从数据中获得。但是，可以通过适当的插值技术使用$y$的初始估计来克服此限制。
• 将$K()$在双边情况下的应用分为空间核及灰度核会削弱估计器的性能，因为它会限制自由度，并且会忽略像素位置及其值之间的相关性。
• 在噪声比较严重的数据中，将K分解为空间核及灰度核将影响估计性能。因为灰度差别都很大，导致灰度权都趋于0从而达不到引入灰度核的效果

以下部分提供了克服此缺点的解决方案。

3.转向核函数

方法原理

观察到式(11)中$K(y_i-y)$的效果实际上是一个像素在某个邻域内的局部梯度估计，并使用这个估计结果来对数据加权。基于这个直观的想法，提出一个两步方法：

• 第一步进行图像局部结构的初始估计。

• 第二步，利用这个初始估计来自适应的"控制”局部核，可以形成一个形状类似沿局部边缘结构方向扩散的椭圆轮廓的核函数。

有了这些局部自适应核，就可以有效的恢复图像的边缘，使得去噪后的细节更加丰富。自适应核的形式如下

$$\mathrm{K}\left({\mathrm{x}}_i-\mathrm{x},y_i-y\right)\equiv K_{{\mathrm{H}}_i^{\mathrm{s}}}\left({\mathrm{x}}_i-\mathrm{x}\right)$$

平滑矩阵${\mathbf{H}}_i^{\mathrm{s}}$是自适应的全矩阵，定义如下：

$${\mathbf{H}}_i^{\mathrm{s}}=h{\mu }_i{\mathbf{C}}_i^{-\frac{1}{2}}$$

$C_i$是基于不同局部灰度值的协方差矩阵，由数据集合决定。$ \mu_i $是一个表征数据局部密度的量（通常设为1），$h$取值较大时会取得较好的降噪效果，但也会让图像变模糊，这时候可以选择合适的$C_1$来调和这种情况。

当$K$取高斯核时，转向核函数就表示为下面的式子：
$$K_{{\mathbf{H}}_i^{\mathfrak{s}}}\left({\mathbf{x}}_i-\mathbf{x}\right)=\frac{\sqrt{\det \left({\mathbf{C}}_i\right)}}{2\pi h^2}\exp \left\{-\frac{{\left({\mathbf{x}}_i-\mathbf{x}\right)}^T{\mathbf{C}}_i\left({\mathbf{x}}_i-\mathbf{x}\right)}{2h^2}\right\}$$

局部边缘结构与梯度协方差相关，原始的协方差矩阵估计在Kernel Regression for Image Processing and Reconstruction[1]这篇文章中

简化后的式子：

$${\mathbf{C}}_i={\gamma }_i{\mathbf{U}}_{{\theta }_i}{\Lambda }_i{\mathbf{U}}_{{\theta }_i}^T$$

$${\mathbf{U}}_{{\theta }_i}=\left[\begin{array}{rr}\cos {\theta }_i& \sin {\theta }_i\\ -\sin {\theta }_i& \cos {\theta }_i\end{array}\right],{\mathbf{\Lambda }}_i=\left[\begin{array}{cc}{\sigma }_i& 0\\ 0& {\sigma }_i^{-1}\end{array}\right]$$

其中$U_{{\theta }_i}$是旋转矩阵，$ \Lambda_i $是伸缩矩阵。协方差矩阵由$\gamma_{i} \ \ \theta_i \ \sigma_i $三个参数决定，分别是尺度，旋转和伸缩参数。迭代求取这些值得方法也在那篇[1]文章中。

这个公式与Robust Fusion of Irregularly Sampled Data Using Adaptive Normalized Convolution这篇文章中提到的归一化卷积公式很接近

• 初始圆形核由$ \Lambda_i $伸缩，长半轴和短半轴由$\sigma_i$给出。
• 之后，由矩阵$U_{{\theta }_i}$对核进行旋转，最后，由尺度参数${\gamma }_i$对内核进行缩放。

如图所示：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-JcZWajfb-1582978320337)(C:\Users\sgb6372\AppData\Roaming\Typora\typora-user-images\image-20191128151454305.png)]

实验

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-mTRoCJXH-1582978320338)(C:\Users\sgb6372\AppData\Roaming\Typora\typora-user-images\image-20191127110121583.png)]

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-bsPtAqDH-1582978320340)(C:\Users\sgb6372\AppData\Roaming\Typora\typora-user-images\image-20191127110144030.png)]

结论

• 核回归可以作为图像去噪和插值算法的有效通用框架
• $l1$范数惩罚纳入内核回归框架中，可以实现针对数据和噪声模型中异常值的鲁棒性