图像降噪：Kernel Regression for Image Processing and Reconstruction

Kernel Regression for Image Processing and Reconstruction

本篇论文是一个基于核回归的图像处理框架，其所提出的核回归方法在图像降噪、融合、超分辨等领域均有广泛应用。本篇论文是 Hiroyuki Takeda的博士论文，是对其前后发表的多篇关于核回归的总结：《Image denoising by adaptive kernel regression》、《Kernel regression for image processing and reconstruction》、《Robust kernel regression for restoration and reconstruction of images from sparse noisy data》。

回归

回归是监督学习中的一个重要问题，用于预测输入变量$X$与输出变量$Y$之间的关系，等价于函数拟合：
$$y_i=m(x_i)+{?}_i;m(.)=E(y|x);E(?)=0$$
例如，使用$n?1$阶多项式就可以完美拟合$n$个观测数据$(x_i,y_i)$，当$n$稍大时，观测数据取值范围以外的数据容易出现[龙格现象](http://www.tlu.ee/~tonu/Arvmeet/Runge’s phenomenon.pdf)，无法实现对全局数据的推断和预测，这自然就引出了一套局部化的非参数估计方法。

• 参数回归(Parametric Regression)
  在参数回归中往往假定输入变量$X$与输出变量$Y$之间函数关系的形式$m$是已知的，即通常给出其参数形式的假定，例如最经典的线性回归模型：
  $$y_i=\beta x_i+{?}_i$$
  而在实际应用过程中，函数关系的形式往往是未知的，不正确函数形式的假设往往会造成对全局数据的正确推断和预测。
• 半参数回归(Sime-Parametric Approach)
  半参数回归中包含参数模型，用以把握变量$Y$的大致走势，非参数模型，用以对$Y$做局部调整，使模型更好地拟合输入变量，所以半参数回归模型同时兼顾数据的确定性规律与非确定性规律，具有很好的适应性：
  $$y_i=\beta x_i+m_z(z_i)+{?}_i$$
• 非参数回归
  非参数回归和参数回归是相对的，完全不假设样本的分布函数，由样本本身来决定模型的结构：
  $$y_i=m(x_i)+{?}_i;$$
  在非参数回归中，核回归(Kernel Regression)是用核函数作为权重对输入变量在其局部邻域内的加权平均作为对输出变量的估计，最常用于降噪当中，Nadaraya-Watson Estimator(NWE)为最早被提出的核回归方法：
  $$\begin{array}{rl}y=\hat{m}(x)& =\frac{{\sum }_{i=1}^n{K}_{h_x}\left(x?x_i\right)y_i}{{\sum }_{i=1}^n{K}_{h_x}\left(x?x_i\right)}\\ & =\sum _{i=1}^n{W}_{h_x}\left(x,x_i\right)y_i\end{array}$$

经典核回归

首先针对一维观测数据，将输入变量$X$与输出变量$Y$间的映射关系做如下表示：
$$y_i=z(x_i)+{?}_i;i=1,2,...,P\text{\hspace{1em}(1)}$$
其中$z$的形式未知，$?$为均值为0的噪声，假设$z$是局部平滑的，为了估计点$x$(在$x_i$附近)的输出值$y$，对上式做$N$阶泰勒展开：
$$\begin{array}{rl}z\left(x_i\right)& \approx z(x)+z^{\prime }(x)\left(x_i?x\right)+\frac{1}{2!}{z}^{\prime \prime }(x){\left(x_i?x\right)}^2+?+\frac{1}{N!}{z}^{(N)}(x){\left(x_i?x\right)}^N\\ & ={\beta }_0+{\beta }_1\left(x_i?x\right)+{\beta }_2{\left(x_i?x\right)}^2+?+{\beta }_N{\left(x_i?x\right)}^N\end{array}\text{\hspace{1em}(2)}$$
以上泰勒展开是对函数$z$在点$x_i$处的局部近似，自然离点$x$越远的采样点的权重越小，即求解式(2)的损失函数如下：
min?{βn}n=0N∑i=1P[yi?β0?β1(xi?x)?β2(xi?x)2???βN(xi?x)N]21hK(xi?xh)(3)\min _{\left\{\beta_{n}\right\}_{n=0}^{N}} \sum_{i=1}^{P}\left[y_{i}-\beta_{0}-\beta_{1}\left(x_{i}-x\right)-\beta_{2}\left(x_{i}-x\right)^{2}-\cdots-\beta_{N}\left(x_{i}-x\right)^{N}\right]^{2} \frac{1}{h} K\left(\frac{x_{i}-x}{h}\right)\tag{3} { βn?}n=0N?min?i=1∑P?[yi??β0??β1?(xi??x)?β2?(xi??x)2???βN?(xi??x)N]2h1?K(hxi??x?)(3)
其中$K(.)$为核函数，$h$为的平滑参数，且核函数需要满足以下条件：
$${\int }_{R^1}tK(t)dt=0,{\int }_{R^1}{t}^{2}K(t)dt=c\text{\hspace{1em}(4)}$$
其中$c$为常数，实验表明，核函数的选择对估测结果的准确性不起决定性作用，因此，一般选择计算复杂度较低的高斯核函数。

一些常用的核函数

除了核函数的选取，泰勒展开的阶数也对最终的近似结果起到非常关键的作用，总的来说$N$越小，图像越平滑，$N$越大越多噪声残留。特别的，当$N=0$时，核回归退化为NWE：
$$\hat{z}(x)=\frac{{\sum }_{i=1}^P{K}_h\left(x_i?x\right)y_i}{{\sum }_{i=1}^P{K}_h\left(x_i?x\right)},K_h(t)=\frac{1}{h}K\left(\frac{t}{h}\right)\text{\hspace{1em}(5)}$$
下面是不同展开阶数的核回归的对含噪数据的拟合情况，当噪声水平较低$N=2$阶核回归的拟合性能较好，当噪声水平较高$N=0、1$阶核回归的拟合性能较好。

不同阶数的展开方式对不同信噪比信号的拟合结果

同一维观测数据相似，将二维观测数据做如下表示：
$$y_i=z\left({\mathbf{x}}_i\right)+{\epsilon }_i;{\mathbf{x}}_i={\left[x_{1i},x_{2i}\right]}^T;i=1,2,?,P\text{\hspace{1em}(6)}$$
其中$\mathbf{x}$为坐标，对应的泰勒展开为：
z(xi)≈z(x)+{?z(x)}T(xi?x)+12!(xi?x)T{Hz(x)}(xi?x)+?=z(x)+{?z(x)}T(xi?x)+12vec?T{Hz(x)}vec?{(xi?x)(xi?x)T}+?(7)\begin{aligned} z\left(\mathbf{x}_{i}\right) & \approx z(\mathbf{x})+\{\nabla z(\mathbf{x})\}^{T}\left(\mathbf{x}_{i}-\mathbf{x}\right)+\frac{1}{2 !}\left(\mathbf{x}_{i}-\mathbf{x}\right)^{T}\{\mathcal{H} z(\mathbf{x})\}\left(\mathbf{x}_{i}-\mathbf{x}\right)+\cdots \\ &=z(\mathbf{x})+\{\nabla z(\mathbf{x})\}^{T}\left(\mathbf{x}_{i}-\mathbf{x}\right)+\frac{1}{2} \operatorname{vec}^{T}\{\mathcal{H} z(\mathbf{x})\} \operatorname{vec}\left\{\left(\mathbf{x}_{i}-\mathbf{x}\right)\left(\mathbf{x}_{i}-\mathbf{x}\right)^{T}\right\}+\cdots \end{aligned}\tag{7} z(xi?)?≈z(x)+{ ?z(x)}T(xi??x)+2!1?(xi??x)T{ Hz(x)}(xi??x)+?=z(x)+{ ?z(x)}T(xi??x)+21?vecT{ Hz(x)}vec{ (xi??x)(xi??x)T}+??(7)
其中$?$为梯度向量，$\mathcal{H}$为Hessian矩阵，$vech(.)$为将对称阵的上三角进行矢量化的半矢量算子(Half-Vectorization-Operator)，具体如下：
$$\begin{array}{l}\text{ vech }\left(\left[\begin{array}{ll}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]\right)={\left[\begin{array}{lll}{a}_{11}& a_{21}& a_{22}\end{array}\right]}^T\\ \text{ vech }??\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}??={\left[\begin{array}{llllll}{a}_{11}& a_{21}& a_{31}& a_{23}& a_{32}& a_{33}\end{array}\right]}^T\end{array}$$
考虑到Hessian矩阵的对称性，将式(7)化简为：
$$z\left({\mathbf{x}}_i\right)={\beta }_0+{\mathbf{\beta }}_1^T\left({\mathbf{x}}_i?\mathbf{x}\right)+{\mathbf{\beta }}_2^T\mathrm{vech}\left\{\left({\mathbf{x}}_i?\mathbf{x}\right){\left({\mathbf{x}}_i?\mathbf{x}\right)}^T\right\}+?\text{\hspace{1em}(8)}$$
其中：
$$\begin{array}{l}{\mathbf{\beta }}_0=z(x)\\ {\mathbf{\beta }}_1={\left[{{\frac{?z(\mathbf{x})}{?x_1}|}_{\mathbf{x}={\mathbf{x}}_i}\frac{?z(\mathbf{x})}{?x_2}|}_{\mathbf{x}={\mathbf{x}}_i}\right]}^T,\\ {\mathbf{\beta }}_2=\frac{1}{2}{\left[{{{\frac{{?}^{2}z(\mathbf{x})}{?x_1^2}|}_{\mathbf{x}={\mathbf{x}}_i}2\frac{{?}^{2}z(\mathbf{x})}{?x_1?x_2}|}_{\mathbf{x}={\mathbf{x}}_i}\frac{{?}^{2}z(\mathbf{x})}{?x_2^2}|}_{\mathbf{x}={\mathbf{x}}_i}\right]}^T\end{array}$$
与式(3)一样，可以通过优化求解如下损失函数来得到{βn}\left\{\boldsymbol{\beta}_{n}\right\}{ βn?}：
min?{βn}n=0N∑i=1P[yi?β0?β1T(xi?x)T?β2Tvech?{(xi?x)(xi?x)T}??]2KH(xi?x)KH(t)=1det?(H)K(H?1t)(8)\min _{\left\{\boldsymbol{\beta}_{n}\right\}_{n=0}^{N}} \sum_{i=1}^{P}\left[y_{i}-\beta_{0}-\boldsymbol{\beta}_{1}^{T}\left(\mathbf{x}_{i}-\mathbf{x}\right)^{T}-\boldsymbol{\beta}_{2}^{T} \operatorname{vech}\left\{\left(\mathbf{x}_{i}-\mathbf{x}\right)\left(\mathbf{x}_{i}-\mathbf{x}\right)^{T}\right\}-\cdots\right]^{2} K_{\mathbf{H}}\left(\mathbf{x}_{i}-\mathbf{x}\right)\\ K_{\mathbf{H}}(\mathbf{t})=\frac{1}{\operatorname{det}(\mathbf{H})} K\left(\mathbf{H}^{-1} \mathbf{t}\right)\tag{8} { βn?}n=0N?min?i=1∑P?[yi??β0??β1T?(xi??x)T?β2T?vech{ (xi??x)(xi??x)T}??]2KH?(xi??x)KH?(t)=det(H)1?K(H?1t)(8)
其中$\mathbf{H}$为$2\times 2$的平滑矩阵，将上式写成矩阵形式：
$$\begin{gathered} \arg \underset{\mathbf{b}}{\min }{\Vert \mathbf{y}?{\mathbf{X}}_{\mathbf{x}}\mathbf{b}\Vert }_{{\mathbf{W}}_{\mathbf{x}}}^2=\arg \underset{\mathbf{b}}{\min }{\left(\mathbf{y}?{\mathbf{X}}_{\mathbf{x}}\mathbf{b}\right)}^T{\mathbf{W}}_{\mathbf{x}}\left(\mathbf{y}?{\mathbf{X}}_{\mathbf{x}}\mathbf{b}\right) \\ \begin{array}{l}\mathbf{y}={\left[\begin{array}{llll}{y}_1& y_2& ?& y_P\end{array}\right]}^T,\mathbf{b}={\left[\begin{array}{llll}{\beta }_0& {\mathbf{\beta }}_1^T& {\mathbf{\beta }}_2^T& ?\end{array}\right]}^T\\ {\mathbf{W}}_{\mathbf{x}}=\mathrm{diag}\left[K_{\mathbf{H}}\left({\mathbf{x}}_1?\mathbf{x}\right),K_{\mathbf{H}}\left({\mathbf{x}}_2?\mathbf{x}\right),?,K_{\mathbf{H}}\left({\mathbf{x}}_P?\mathbf{x}\right)\right]\\ {\mathbf{X}}_{\mathbf{x}}=?\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\}& ?\\ 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\}& ?\\ ?& ?& ?& ?\\ 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\}& ?\end{array}?\end{array}\text{\hspace{1em}(9)} \end{gathered}$$
此加权最小二乘解得：
$$\hat{z}(\mathbf{x})={\hat{\beta }}_0={\mathbf{e}}_1^T{\left({\mathbf{X}}_{\mathbf{x}}^T{\mathbf{W}}_{\mathbf{x}}{\mathbf{X}}_{\mathbf{x}}\right)}^{?1}{\mathbf{X}}_{\mathbf{x}}^T{\mathbf{W}}_{\mathbf{x}}\mathbf{y};e_1=[1,0,...,0]\in R^{N+1}\text{\hspace{1em}(10)}$$
可以将上式简单表示成局部线性滤波：
$$\hat{z}(\mathbf{x})=\sum _{i=1}^P{W}_i(\mathbf{x};N,\mathbf{H})y_i\text{\hspace{1em}(11)}$$
将上式(10)中的${\mathbf{X}}_{\mathbf{x}}^T{\mathbf{W}}_{\mathbf{x}}{\mathbf{X}}_{\mathbf{x}}$为$(N+1)\times (N+1)$的分块矩阵：
$${\mathbf{X}}_{\mathbf{x}}^T{\mathbf{W}}_{\mathbf{x}}{\mathbf{X}}_{\mathbf{x}}=?\begin{array}{cccc}{s}_{11}& {\mathbf{s}}_{12}& {\mathbf{s}}_{13}& ?\\ s_{21}& {\mathbf{s}}_{22}& {\mathbf{s}}_{23}& ?\\ s_{31}& {\mathbf{s}}_{32}& {\mathbf{s}}_{33}& ?\\ ?& ?& ?& ?\end{array}?$$
当$N=2$时分块矩阵内的元素如下：
$$\begin{array}{rl}{s}_{11}& =\sum _{i=1}^P{K}_{\mathbf{H}}\left({\mathbf{x}}_i?\mathbf{x}\right)\\ {\mathbf{s}}_{12}& ={\mathbf{s}}_{21}^T=\sum _{i=1}^P{\left({\mathbf{x}}_i?\mathbf{x}\right)}^T{K}_{\mathbf{H}}\left({\mathbf{x}}_i?\mathbf{x}\right)\\ {\mathbf{s}}_{22}& =\sum _{i=1}^P\left({\mathbf{x}}_i?\mathbf{x}\right){\left({\mathbf{x}}_i?\mathbf{x}\right)}^T{K}_{\mathbf{H}}\left({\mathbf{x}}_i?\mathbf{x}\right)\\ {\mathbf{s}}_{13}& ={\mathbf{s}}_{31}^T=\sum _{i=1}^P{\mathrm{vech}}^T\left\{\left({\mathbf{x}}_i?\mathbf{x}\right){\left({\mathbf{x}}_i?\mathbf{x}\right)}^T\right\}{K}_{\mathbf{H}}\left({\mathbf{x}}_i?\mathbf{x}\right)\\ {\mathbf{s}}_{23}& ={\mathbf{s}}_{32}^T=\sum _{i=1}^P\left({\mathbf{x}}_i?\mathbf{x}\right){\mathrm{vech}}^T\left\{\left({\mathbf{x}}_i?\mathbf{x}\right){\left({\mathbf{x}}_i?\mathbf{x}\right)}^T\right\}{K}_{\mathbf{H}}\left({\mathbf{x}}_i?\mathbf{x}\right)\\ {\mathbf{s}}_{33}& =\sum _{i=1}^P\mathrm{vech}\left\{\left({\mathbf{x}}_i?\mathbf{x}\right){\left({\mathbf{x}}_i?\mathbf{x}\right)}^T\right\}{\mathrm{vech}}^T\left\{\left({\mathbf{x}}_i?\mathbf{x}\right){\left({\mathbf{x}}_i?\mathbf{x}\right)}^T\right\}{K}_{\mathbf{H}}\left({\mathbf{x}}_i?\mathbf{x}\right)\end{array}$$
于是式(11)中的权重矩阵可以表示如下：
Wi(x;0,H)=KH(xi?x)s11Wi(x;1,H)={1?s12s22?1(xi?x)}KH(xi?x)s11?s12s22?1s21Wi(x;2,H)=[1?S12S22?1(xi?x)?S13S33?1vech?{(xi?x)(xi?x)T}]KH(xi?x)s11?S12S22?1s21?S13S33?1s31S12=s12?s13s33?1s32,S22=s22?s23s33?1s32S13=s13?s12s22?1s23,S33=s33?s32s22?1s23\begin{array}{l} W_{i}(\mathbf{x} ; 0, \mathbf{H})=\frac{K_{\mathbf{H}}\left(\mathbf{x}_{i}-\mathbf{x}\right)}{s_{11}} \\ W_{i}(\mathbf{x} ; 1, \mathbf{H})=\frac{\left\{1-\mathbf{s}_{12} \mathbf{s}_{22}^{-1}\left(\mathbf{x}_{i}-\mathbf{x}\right)\right\} K_{\mathbf{H}}\left(\mathbf{x}_{i}-\mathbf{x}\right)}{s_{11}-\mathbf{s}_{12} \mathbf{s}_{22}^{-1} \mathbf{s}_{21}}\\ W_{i}(\mathbf{x} ; 2, \mathbf{H})=\frac{\left[1-\mathbf{S}_{12} \mathbf{S}_{22}^{-1}\left(\mathbf{x}_{i}-\mathbf{x}\right)-\mathbf{S}_{13} \mathbf{S}_{33}^{-1} \operatorname{vech}\left\{\left(\mathbf{x}_{i}-\mathbf{x}\right)\left(\mathbf{x}_{i}-\mathbf{x}\right)^{T}\right\}\right] K_{\mathbf{H}}\left(\mathbf{x}_{i}-\mathbf{x}\right)}{s_{11}-\mathbf{S}_{12} \mathbf{S}_{22}^{-1} \mathbf{s}_{21}-\mathbf{S}_{13} \mathbf{S}_{33}^{-1} \mathbf{s}_{31}}\\ \mathbf{S}_{12}=\mathbf{s}_{12}-\mathbf{s}_{13} \mathbf{s}_{33}^{-1} \mathbf{s}_{32}, \quad \mathbf{S}_{22}=\mathbf{s}_{22}-\mathbf{s}_{23} \mathbf{s}_{33}^{-1} \mathbf{s}_{32} \\ \mathbf{S}_{13}=\mathbf{s}_{13}-\mathbf{s}_{12} \mathbf{s}_{22}^{-1} \mathbf{s}_{23}, \quad \mathbf{S}_{33}=\mathbf{s}_{33}-\mathbf{s}_{32} \mathbf{s}_{22}^{-1} \mathbf{s}_{23} \end{array} Wi?(x;0,H)=s11?KH?(xi??x)?Wi?(x;1,H)=s11??s12?s22?1?s21?{ 1?s12?s22?1?(xi??x)}KH?(xi??x)?Wi?(x;2,H)=s11??S12?S22?1?s21??S13?S33?1?s31?[1?S12?S22?1?(xi??x)?S13?S33?1?vech{ (xi??x)(xi??x)T}]KH?(xi??x)?S12?=s12??s13?s33?1?s32?,S22?=s22??s23?s33?1?s32?S13?=s13??s12?s22?1?s23?,S33?=s33??s32?s22?1?s23??
可以看出，传统的核回归都是基于对邻域内数据进行加权来得到最终的预测结果，而高阶的核回归($N>0$)从此物理意义上讲与0阶核回归相比仅仅是具有一个更加复杂的权重函数而已。下图是对均匀采样分别进行0、1、2阶展开的等效权重剖面图，可以观察到，对0阶展开与1阶展开的等效权重是重合的，这一现象可以推广为${\mathbf{X}}_{\mathbf{x}}^T{\mathbf{W}}_{\mathbf{x}}{\mathbf{X}}_{\mathbf{x}}$中的奇数项${\mathbf{s}}_{2j,2k+1}、{\mathbf{s}}_{2k+1,2j}$中所包含的元素均接近于0，对于均匀采样而言，使用偶数阶对局部做展开比使用奇数阶对局部做展开的性价比更高，因为$2N$阶展开的结果与$2N+1$阶相同。在下图中还可以观察到与$N=0、1$不同的是，$N=2$时的等价核包含部分为负数的元素，在频域内其频率传递函数的频率截止速度高于比$N=0、1$，类似巴特沃斯低通滤波(这部分为负数的元素，可能是由Ringing造成的)。

a为均匀采样的样本；b为分别对其进行0、1、2解展开的等效权重

但是上式结论对于非均匀采样的样本数据并不成立，非均匀采样的等价权重如下图所示，由于奇数阶分块子矩阵不接近0，$N=0、1$的等价权重不再相同。

a为非均匀采样的样本；b~d为分别对其进行0、1、2解展开的等效权重

对于经典的核回归算法，平滑矩阵$\mathbf{H}$也在很大程度上影响最终的估计结果，$\mathbf{H}$负责将式(3)中的损失函数的权重更多的分配给有效样本，对于二维图像数据，将$\mathbf{H}$定义为一个$2\times 2$的对称阵。

数据自适应的核回归

在式(3)中定义的全图统一的核函数$K_{\mathbf{H}}$的基础上引入由灰度距离控制的新的权重值，为局部自适应核回归，如下图所示，为了在对边缘或者强纹理的区域能够做到保护边缘信息的同时很好的去除噪声，比较直观的思路是根据图像的边缘信息来各向异性地分配核函数的权重，即根据样本的分布来自适应核函数中的平滑矩阵$H$

a：传统的基于样本距离的核函数；b：数据自适应的核函数

为达成以上数据自适应的分配核函数权重的目的，首先引入双边滤波(BLF)的权重函数作为核函数，此时式(3)中的损失函数如下：
$$\mathrm{K}\left({\mathrm{x}}_i?\mathrm{x},y_i?y\right)\equiv K_{{\mathbf{H}}_s}\left({\mathbf{x}}_i?\mathbf{x}\right)K_{h_r}\left(y_i?y\right)$$
其中${\mathbf{H}}_s=\left[\begin{array}{cc}{h}_s& 0\\ 0& h_s\end{array}\right]$，当$N=0$时，BLF可以看作是自适应的NWE：
$$\hat{z}(\mathbf{x})=\frac{{\sum }_{i=1}^P{K}_{{\mathbf{H}}_s}\left({\mathbf{x}}_i?\mathbf{x}\right)K_{h_r}\left(y_i?y\right)y_i}{{\sum }_{i=1}^P{K}_{{\mathbf{H}}_s}\left({\mathbf{x}}_i?\mathbf{x}\right)K_{h_r}\left(y_i?y\right)}$$
此时的损失函数如下：
min?{βn}n=0N∑i=1P[yi?β0?β1T(xi?x)T?β2Tvech?{(xi?x)(xi?x)T}??]2K(xi?x,yi?y)(12)\min _{\left\{\boldsymbol{\beta}_{n}\right\}_{n=0}^{N}} \sum_{i=1}^{P}\left[y_{i}-\beta_{0}-\boldsymbol{\beta}_{1}^{T}\left(\mathbf{x}_{i}-\mathbf{x}\right)^{T}-\boldsymbol{\beta}_{2}^{T} \operatorname{vech}\left\{\left(\mathbf{x}_{i}-\mathbf{x}\right)\left(\mathbf{x}_{i}-\mathbf{x}\right)^{T}\right\}-\cdots\right]^{2} \mathrm{~K}\left(\mathbf{x}_{i}-\mathbf{x}, y_{i}-y\right)\tag{12} { βn?}n=0N?min?i=1∑P?[yi??β0??β1T?(xi??x)T?β2T?vech{ (xi??x)(xi??x)T}??]2 K(xi??x,yi??y)(12)
但是这种直接将核函数$K$分解为两个核函数的做法并不一定能带来核回归估计的性能提升，因为当噪声较大时，灰度差都比较大，灰度核的权重趋于0，最终导致核函数内大部分权重趋于0，失去了核函数对回归过程的约束作用。基于此，作者引用Steering Kernel来对经典核回归方法中各向同性的核函数进行优化，其中最核心的思想就是根据像素点在其邻域内的梯度分布来度量图像在局部邻域内的Dominant Orientation，即梯度变化的主要方向，利用此主方向来进一步控制局部核函数的权重分配，如下图所示：

如果是在边缘区域，此时的局部自适应核函数会将权重更多的沿着边缘方向分配，在去噪的过程中有效地保护图像边缘信息，使得去噪后的图像拥有更加丰富的细节信息，这一方法类似于常用的方向滤波算法，但此方法能够表示的方向更多，更精细，关于Steering Kernel的更多细节可以参考论文《Multiscale principal components analysis for image local orientation estimation》，将此时的自适应核函数做如下表示：
$$\mathrm{K}\left({\mathbf{x}}_i?\mathbf{x},y_i?y\right)\equiv K_{{\mathbf{H}}_i^{\mathrm{s}}}\left({\mathbf{x}}_i?\mathbf{x}\right)$$
不同于灰度平滑矩阵${\mathbf{H}}_s$，${\mathbf{H}}_i^{\mathrm{s}}=h{\mu }_i{\mathbf{C}}_i^{?\frac{1}{2}}$是一个$2\times 2$的稠密矩阵，其中${\mathbf{C}}_i$可以简单的用局部灰度值的协方差矩阵来近似，即局部边缘的主方向与梯度协方差是强相关或者说近似等价的：
$${\hat{\mathbf{C}}}_i\approx \left[\begin{array}{cc}{\sum }_{{\mathbf{x}}_j\in w_i}{z}_{x_1}\left({\mathbf{x}}_j\right)z_{x_1}\left({\mathbf{x}}_j\right)& {\sum }_{{\mathbf{x}}_j\in w_i}{z}_{x_1}\left({\mathbf{x}}_j\right)z_{x_2}\left({\mathbf{x}}_j\right)\\ {\sum }_{{\mathbf{x}}_j\in w_i}{z}_{x_1}\left({\mathbf{x}}_j\right)z_{x_2}\left({\mathbf{x}}_j\right)& {\sum }_{{\mathbf{x}}_j\in w_i}{z}_{x_2}\left({\mathbf{x}}_j\right)z_{x_2}\left({\mathbf{x}}_j\right)\end{array}\right]$$
其中$w_i$表示局部邻域窗口，局部梯度的主方向与其对应的协方差矩阵的特征向量相关，这一相关性可以通过对局部邻域内多个梯度向量进行主成分分析(Principal components analysis，PCA)来得到。首先定义梯度向量矩阵$G$，$G$是一个$N\times 2$的矩阵，对$G$进行奇异值分解，分解的结果如下：
$${\mathbf{G}}_i=?\begin{array}{cc}?& ?\\ z_{x_1}\left({\mathbf{x}}_j\right)& z_{x_2}\left({\mathbf{x}}_j\right)\\ ?& ?\end{array}?={\mathbf{U}}_i{\mathbf{S}}_i{\mathbf{V}}_i^T,\mathbf{x}\in w_i$$
分解的结果中$S$的对角线元素为特征值$s_1、s_2$，$V$的第一列表示梯度的主方向，第二列表示灰度的主方向(这两个方向正交)，根据此奇异值分解的结果，计算决定平滑矩阵$H$的3个参数：

• 伸缩矩阵
  $${\mathbf{\Lambda }}_i=\left[\begin{array}{cc}{\rho }_i& 0\\ 0& {\rho }_i^{?1}\end{array}\right];{\rho }_i=\frac{s_1+{\lambda }^{\prime }}{s_2+{\lambda }^{\prime }},{\lambda }^{\prime }\ge 0$$

• 旋转矩阵
  $${\mathbf{U}}_{{\theta }_i}=\left[\begin{array}{cc}\cos {\theta }_i& \sin {\theta }_i\\ ?\sin {\theta }_i& \cos {\theta }_i\end{array}\right]$$

• 尺度参数

$${\gamma }_i=\left(\frac{s_1{s}_2+{\lambda }^{\prime \prime }}{N}\right)$$

其对原核函数的作用过程如图所示：

方向自适应的核函数的变换过程

最终的平滑矩阵为：
$${\mathbf{C}}_i={\gamma }_i{\mathbf{U}}_{{\theta }_i}{\mathbf{\Lambda }}_i{\mathbf{U}}_{{\theta }_i}^T$$
此时的自适应核函数为：
$$K_{{\mathbf{H}}_i^{\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{e}\mathrm{r}}}\left({\mathbf{x}}_i?\mathbf{x}\right)=\frac{\sqrt{\det \left({\mathbf{C}}_i\right)}}{2\pi h^2{\mu }_i^2}\exp \left\{?\frac{{\left({\mathbf{x}}_i?\mathbf{x}\right)}^T{\mathbf{C}}_i\left({\mathbf{x}}_i?\mathbf{x}\right)}{2h^2{\mu }_i^2}\right\}\text{\hspace{1em}(13)}$$
本文的核心精神即是用式(13)中的自适应核函数替代经典核回归算法中的高斯核函数，以到达对未知数据更好的预测效果。