大家好，又见面了，我是你们的朋友全栈君。

KAZE 特征

本文是在研读论文KAZE Features的基础上以及看了作者源代码后的一个总结，由于能力有限，总结可能存在误差，读者可以在论文作者项目主页上下载源代码，更深入熟悉KAZE特征

一、特征简介

虽然sift,surf等特征在图像特征提取方面已经取得比较好的成绩，但是这类特征都是基于一个通过高斯核进行的线性的尺度空间进行特征检测的，相同尺度下每个点的变换是一样的，由于高斯函数是低通滤波函数，会平滑图像边缘，以至图像损失掉许多细节信息。针对这一问题，作者提出了一种基于非线性尺度空间的特征点检测方法，该非线性尺度空间保证了图像边缘在尺度变化中信息损失量非常少，从而极大地保持了图像细节信息。

二、KAZE 特征检测

1.1 非线性尺度空间构建

1.1.1 非线性尺度滤波
文中非线性尺度空间的构建主要基于非线性扩散滤波原理，非线性扩散滤波的基本公式是:

$$\dfrac{\partial L}{\partial t}=div(c(x,y,t)\cdot \nabla L)$$

其中

div
$div$和

∇
$\nabla$分别表示散度和梯度，而函数

c(x,y,t)
$c(x,y,t)$表示扩散的传导函数，正是该函数的引入使得扩散能够适应图像局部特征，在本文中该传导函数主要由梯度幅值控制，参数

t
$t$是尺度参数。

1.1.2 Perona Malik 扩散方程
非线性滤波在图像方面的应用是由Perona 和 Malik 于1990年文章Scale-space and edge detection using annisotropic diffusion提出来的，为了尽量使平滑在区域内而不是区域间，减少在边缘处的扩散，Perona提出了一个使得传导函数$c$依赖图像梯度幅值的函数：

$$c(x,y,t)=g(\left|{\nabla L_\sigma(x,y,t)}\right|)$$

其中

Lσ
$L_\sigma$是原始图像高斯平滑后的梯度

Perona 提出了两个传导函数：

$$g1=exp\left({-\dfrac{\left|\nabla L_\sigma\right|^2}{k^2}}\right)$$

$$g2=\dfrac{1}{1+\frac{\left|\nabla L_\sigma\right|^2}{k^2}}$$

而后Weircket又提出了另外一个方程:

$$g3=\begin{cases}1& ,\left|\nabla L_\sigma\right|^2=0 \\1-exp\left(-\dfrac{3.315}{\left(\left|\nabla L_\sigma\right| \middle/k \right)^8}\right) &,\left|\nabla L_\sigma\right|^2>0\end{cases}$$

其中参数

k
$k$表示控制扩散的对比因子，值越大，同一点处相对扩散越小。在文中作者主要通过图像梯度百分比完成，一般是k是梯度直方图中70%百分比的值，其求解伪代码如下：

ComputeContrast_K(img, kpercentile=0.7, nbins = 300):
    //calculate gradient
    foreach point in img:
        calculate dx dy in x and y dirctions 
        mog = dx*dx + dy*dy
    hmax = max[mog in img]

    //calculate gradient histogram
    foreach  point in gradient:
        nbin = floor(nbins*(mog/hmax));
        hist[nbin]++

    //find the perc of the histogram percentile
    kp = the bin that make sum[0-kp] > kpercentile * npoints
    k = hmax*((float)(kp)/(float)nbins);
    return k

1.1.3 AOS方法求解非线性扩散方程
对于上述的非线性扩散偏微分方程没有解析解，只能通过一些数学方法近似求解该偏微分方程，本文主要引用了文章Efficient and Reliable Schemesfor Nonlinear Diffusion Filtering中的AOS求解方案。上述方程的离散化表达式如下：

$$\dfrac{L^{i+1} - L^i}{\tau}=\sum_{l=1}^mA_l(L^i)L^{i+1}$$

因此有:

$$L^{i+1}=\left(I-\tau\sum_{l=1}^mA_l(L^i)\right)^{-1} L^i$$



Li+1
$L^{i+1}$可以通过解方程得到,其中

Al
$A_l$是传导函数

c(x,y,t)
$c(x,y,t)$在每一个维度的构造函数，其函数如下:

$$A_{ij}(L^l):=\begin{cases} \\\dfrac{{g_i}^l+{g_j}^l}{2h^2}&\left[j\in N(i)\right], \\- \sum_{n\in N(i)}\dfrac{{g_i}^l+{g_j}^l}{2h^2}&(j=i) \\0 & else\end{cases}$$

该构造函数将使得

Aij
$A_{ij}$变成一个三对角矩阵：

$$B=I-\tau A_l(L^i)=\begin{pmatrix}\alpha_1&\beta_1 &0&&\\\gamma_1&\alpha_2&\beta_2&&\\&\ddots&\ddots&\ddots \\&&\gamma_{N-2}&\alpha_{N-1}&beta_{N-1} \\&&0&\gamma_{N-1}&\alpha_N\end{pmatrix}$$

代入方程式有:

$$BL^{i+1}=L^i$$

该方程式是三角矩阵函数式，其求解方法可以通过
Thomas算法得到求解:

（1）LR分解：将B分解成一个下双对角矩阵L和一个上双对角矩阵R:

$$F=\begin{pmatrix}1&&&& \\l_1&1 \\&\ddots&&\ddots \\&&&l_{N-1}&1\end{pmatrix}$$

$$R=\begin{pmatrix}m_1&\gamma_1&&\\&\ddots&\ddots \\&&m_{N-1}&\gamma_{N-1} \\&&&m_N \end{pmatrix}$$

比较三个矩阵，很容易可以发现

γi=βi
$\gamma_i=\beta_i$,而

mi
$m_i$和

li
$l_i$可以通过以下算法获得：

$$\begin{eqnarray*}&&m_1 :=\alpha_1 \\&&for i=1,2, ...,N-1: \\&&l_i:=\gamma_i /m_i \\&&m_{i+1} :=\alpha_{i+1} - l_i\beta_i \\\end{eqnarray*}$$

原方程可以为：

$$FRL^{i+1}=L^i$$

(2)前向替换:将

RLi+1
$RL^{i+1}$看作一个整体未知数

H
$H$,则可以先解方程



FR=Li

$FR=L^i$,由于矩阵

F
$F$的特殊性，可以很容易解出结果



y

$y$

(3)后向替换:在前向替换的基础上有

RLi+1=y
$RL^{i+1}=y$，而R也是一个特殊矩阵，可以容易计算得

Li+1
$L^{i+1}$

1.1.3 构建尺度空间
文章采取的尺度空间的构建与Sift尺度空间的构建比较相似，都是通过指数步长的系列组合(O个组，S个层)来离散化尺度空间的，但与sift不同的是文章中的各组图像都是对原图像的操作而不是对图像进行降采样，而且在同一尺度下各个像素点的变化是非线性的。各个层之间的尺度关系如下:

$$\sigma_i(o,s)=\sigma_02^{o+s/S} ,o \in [0...O-1], s\in [0...S-1], i \in[0...N]$$

其中

σ0
$\sigma_0$是基本尺度，

o
$o$表示组序号，



s

$s$表示所在组的层序号，

N
$N$表示总的层数，为了进行非线性扩散滤波，需要将尺度空间中的尺度单元转换到时间单元概念，其转换如下:

$$t_i=\dfrac{1}{2}{\sigma_i}^2, i=\left\{0...N\right\}$$

然后根据AOS方法进行非线性扩散滤波得到每一层，其伪代码如下:

gaussianImg = Gaussian_Smooth(img, sigma0)
k = ComputeContrast_K(img)
L[0] = gaussianImg
for i=1 to N-1:
    L[i]=L[i-1]
    calculate L[i].dx L[i].dy
    L[i].flow = g2(L[i].dx L[i].dy, k)      //calculate conductivity
    L[i].etime = simga0*pow(o+s/S, 2) AOS_Step_Scalar(L[i-1], L[i].flow, L[i].etime-L[i-1].etime) //AOS algorithm to alculate
                                                                  //L[i] diffuse from L[i]

其中ComputeContrast_K 函数为1.1.2中所述求对比度函数，AOS_Step_Scalar 为通过AOS算法进行非线性滤波的扩散，1.1.3中已有描述，在实现的时候需要对x,y方向方便进行非信息扩散然后在取平均值

1.2 特征点检测

特征点的检测与Sift特征点检测也非常相似，首先计算各点在本层，上下层3X3的立方体邻域内响应值是否为极值，再去掉重复点，再亚像素化精确位置，所不同的是各个像素点处响应值的计算。文章中响应值是采用多尺度下尺度归一化Hessian行列式计算的如下:

$$L_{hessian}=\sigma^2(L_{xx}L_{yy}-{L_{xy}}^2),$$

其中

Lxx
$L_{xx}$

Lyy
$L_{yy}$

Lxy
$L_{xy}$分别表示x方向2阶导数，y方向2阶导数，xy方向混合导数，其主要过程描述如下：

（1）计算各个像素点处响应值，响应值的计算如上

（2）查找在3×3空间邻域内的极大值

（3）删除重复特征点，两点在相邻层并且距离在尺度范围内认为是重复的。

（4）亚像素化，精确定位，亚像素化的过程与sift中的一样，通过泰勒表达式求导数，倒数为0的地方即为极值所在处，该位置即为偏移量：

$$\hat{X}=-{\dfrac{\partial^2D}{\partial X}}^{-1}\dfrac{\partial D}{\partial X}$$

可以通过解方程:

$${\dfrac{\partial^2D}{\partial X}}\hat{X}=- \dfrac{\partial D}{\partial X}$$

得到

X^(x,y,s)
$\hat{X} (x,y,s)$,而

∂2D∂X
$\dfrac{\partial^2D}{\partial X}$可以即为对应

Hessian(x,y,s)
$Hessian(x,y,s)$矩阵。

三 KAZE特征描述

3.1 计算特征点主方向

文章所用的特征点主方向计算方法与SURF算法如出一辙。主要过程如下：
（1） 以特征点为中心，半径为6s区域内的点求$x,y$方向倒数$dx,dy$,然后以标准差为2.5s的高斯函数进行加权
（2）建立一个$\dfrac{\pi}{3}$的扇形区域，该扇形以0.15弧度为步长进行旋转，然后计算落角度在该扇形区域内点的x,y方向导数和
（3）获取所有扇形区域中使得x,y方向导数和的平方和的最大的扇形，主方向即为该扇形区域x,y方向导数和的反正切结果
伪代码如下：

foreach keypoint in keypoints:
    for all points S around keypoint within a distance of 6*sigma:
        w = gaussian(x,y, 2.5sigma)
        calculate dx*w, dy*w, direction
    for direction_start = 0 to 2pi with a step 0.15:
        direction_end = direction_start  + 3/pi
        for all points S around keypoint within a distance of 6*sigma:
            if direction[p] is in [direction_start, direction_end ]:
                SumX += dx
                SumY += dy
                mog = SumX*SumX + SumY*SumY
                if mog > magMog:
                    angle = arctan(SumY, SumX)
    keypoint.angle = angle              

3.2建立描述子

文章中描述子的建立是采用M-SURF描述子，对于一个尺度为$\sigma_i$的特征，取该特征点邻域$24\sigma_i$x$24\sigma_i$的矩形，再将该矩形分成4×4个$9\sigma_i$x$9\sigma_i$的小矩形区域(每个相邻子区域有$4\sigma_i$个像素是重叠的);计算邻域内所有点的x,y方向导数，并且以以标准差为2.5$\sigma_i$的高斯函数进行加权，然后统计每个子区域内的x,y方向导数sum_dx,sum_dy以及方向导数绝对值的和sum_mdx, sum_mdy，对于每个子区域再通过标准差为1.5$\sigma_i$的高斯函数进行加权，最终每个子区域就能形成一个4维特征,4×4个子区域能够形成64维特征，如果想将64维特征扩展到128维特征，可以对方向导数的正负值分别统计。
伪代码如下：

foreach keypoint in keypoints:
    for all points around keypint in a rectangle of 24s x 24s
    divide the rectange into 4x4 subregions of size 9s x 9s
    foreach subregion in rectangel:
        foreach point in the subregion:
            w1 = gaussian(x,y, 2.5sigma)
            sum_dx = w*dx;
            sum_dy = w*dy;
            sum_mdx = abs(w*dx)
            sum_mdy = abs(w*dy)
        w2 = gaussian(center_x,center_y, 1.5)
        sum_dx  *= w2 
        sum_dy  *= w2 
        sum_mdx  *= w2 
        sum_mdy  *= w2 
collect all sum_dx, sum_dy, sum_mdx, sum_mdy in each subregion as a final descriptor of keypoint