三维点云课程—Harris特征点介绍

在讲解三维的 $Harris$之前，先来介绍一下图片二维的 $Harris$的原理

1.图片的Harris特征值介绍

给与一张图片，考虑一个$patch$ $x,y\in \Omega$。此时在$[u,v]$单位下移动这个$patch$，图片的$intensity$变化($intensity$此时为像素)为
$$E(u,v)=\sum _{x,y\in \Omega }w(x,y)[I(x+u,y+v)-I(x,y){]}^2$$
其中$w(x,y)$表示强度变化的权重，通常用均值或者高斯表示，如下

左图为权重按照均值分布，表示当$x,y\in \Omega$，$w(x,y)=1$，否则$w(x,y)=0$。右图为权重按照高斯分布。为了推导方便，本文采用均值分布的方式进行介绍。

对$I(x+u,y+v)$在$x,y$处进行泰勒公式展开
$$I(x+u,y+v)=I(x,y)+u{f}_x(x,y)+v{f}_y(x,y)+\frac{1}{2!}[u^2{f}_{xx}(x,y)+uv{f}_{xy}(x,y)]+v^2{f}_{yy}(x,y)+...$$
取一阶泰勒公式
$$I(x+u,y+v)\approx I(x,y)+u{f}_x(x,y)+v{f}_y(x,y)$$
那么
$$\begin{gathered} E(u,v)=\sum _{x,y\in \Omega }(I(x+u,y+v)-I(x,y){)}^2 \\ \approx \sum _{x,y\in \Omega }(u{f}_x(x,y)+v{f}_y(x,y){)}^2 \\ =\sum _{x,y\in \Omega }{u}^2{I}_x^2+2uv{I}_x{I}_y+v^2{I}_y^2 \\ =\sum _{x,y\in \Omega }\left[\begin{array}{cc}u& v\end{array}\right]\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right] \\ =\left[\begin{array}{cc}u& v\end{array}\right](\sum _{x,y\in \Omega }\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right])\left[\begin{array}{c}u\\ v\end{array}\right] \\ =\left[\begin{array}{cc}u& v\end{array}\right]M\left[\begin{array}{c}u\\ v\end{array}\right] \end{gathered}$$

其中$M=\sum _{x,y\in \Omega }\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]$表示为图片梯度的协方差矩阵。定义$I_i=[I_{xi},I_{yi}{]}^T$,那么$M={\sum }_i{I}_i{I}_i^T$。

注意：由于图像是离散的，因此不能通过求导的方式来求像素的梯度，通过差商来表示像素的梯度，差商分为前向差商$\frac{\partial f(x,y)}{\partial x}=f(x+1,y)-f(x,y)$，后向差商$\frac{\partial f(x,y)}{\partial x}=f(x,y)-f(x-1,y)$，中心差商$\frac{\partial f(x,y)}{\partial x}=\frac{f(x+1,y)-f(x-1,y)}{2}$等方法。

对于上述图片中经常遇到的三种情况，对应的$I_i$如下

Liner Edge：$I_i=[I_{xi},0{]}^T$

Flat：$I_i=[0,0{]}^T$

Corner：$I_i=[I_{xi},I_{yi}{]}^T$

以每一个像素点的$I_x$为$x$轴，$I_y$为$y$轴，得到

分析如下，求解协方差矩阵$M$的两个特征值${\lambda }_1,{\lambda }_2$，如果其中一个特征值远远大于另一个特征值，则表示为Linear Edge，类似于上图椭圆的两个轴；如果两个特征值都比较小，则表示为Flat,类似于上图小圆的半径；如果两个特征值都比较大，则表示为Corner,类似于上图大圆的半径。图形表示如下

那么如何判断一个特征的大小呢，采用响应函数(Response Function)$R$来表达

Harris & Stephens(1988):
$$\begin{gathered} R=det(M)-k(trace(M){)}^2 \\ det(M)={\lambda }_1{\lambda }_2 \\ trace(M)={\lambda }_1+lambd{a}_2 \\ k\in [0.04,0.06] \end{gathered}$$
Kanade & Tomasi(1994):
$$R=min({\lambda }_1,{\lambda }_2)$$
Nobel(1994):
$$R=\frac{det(M)}{trace(M)+\epsilon }$$
在实际应用中采用上述三种的任何一种都可以表示特征值的阈值，但为了方便起见，通常采用第二种方式，记作Tomasi Response。

---

2.点云中的Harris特征值介绍

因为$Harris$算法工作在一个连续的函数上面(要进行泰勒公式的展开,对应于2维就是图像的像素)，而对于点云来说，分为有intensity的点云,如激光雷达的反射率,和没有intensity的点云,如点云数据中只有$x,y,z,(r,g,b)$。

2.1 Harris3D With Intensity

假设$intensity$是一个连续的函数
$$I(x,y,z):R^3\to R,[x,y,z]\in \Omega$$
假设存在一个小的"位移" $[u,v,w]$

那么$intensity$的变化如下
$$E(u,v,w)=\sum _{x,y,z\in \Omega }[I(x+u,y+v,z+w)-I(x,y,z){]}^2$$
泰勒公式化简得
$$E(u,v,w)=\left[\begin{array}{ccc}u& v& w\end{array}\right]M\left[\begin{array}{c}u\\ v\\ w\end{array}\right],M=\sum _{x,y,z\in \Omega }\left[\begin{array}{ccc}{I}_x^2& I_x{I}_y& I_x{I}_z\\ I_x{I}_y& I_y^2& I_y{I}_z\\ I_x{I}_z& I_y{I}_z& I_z^2\end{array}\right]$$
**其中$M$是表面Intensity的协方差矩阵。**但是不同于图像的是，三维点云中的点是离散的，那么怎么去求一个点的梯度呢?

• 定义我们要求的一个点$p=[p_x,p_y,p_z{]}^T$
• 在领域$\Omega$内，有一些点 { x i = [ x i , y i , z i ] } \{x_i=[x_i,y_i,z_i]\} {xi​=[xi​,yi​,zi​]}
• $p$周围的$intensity$梯度是一个向量，记作$e=[e_x,e_y,e_z{]}^T$
  • $e$方向为$intensity$增加最大方向
  • $||e||$表示$intensity$变化的速率
• 理论上，我们有

$$\begin{gathered} (x_1-p_x)e_x+(y_1-p_y)e_y+(z_1-p_z)e_z=I(x_1,y_1,z_1)-I(p_x,p_y,p_z) \\ ...=... \\ (x_i-p_x)e_x+(y_i-p_y)e_y+(z_i-p_z)e_z=I(x_i,y_i,z_i)-I(p_x,p_y,p_z) \\ ...=... \end{gathered}$$

标量形式
$$(x_i-p_x)e_x+(y_i-p_y)e_y+(z_i-p_z)e_z=I(x_i,y_i,z_i)-I(p_x,p_y,p_z)$$
向量形式
$$\begin{gathered} x_i^{{}^{\prime }T}e=\Delta I_i \\ {x}_i^{{}^{\prime }}=[x_i^{{}^{\prime }},y_i^{{}^{\prime }},z_i^{{}^{\prime }}{]}^T=[x_i-p_x,y_i-p_y,z_i-p_z{]}^T \end{gathered}$$
矩阵形式
$$Ae=b,A=\left[\begin{array}{ccc}{x}_1^{{}^{\prime }}& y_1^{{}^{\prime }}& z_1^{{}^{\prime }}\\ ...& ...& ...\\ x_i^{{}^{\prime }}& y_i^{{}^{\prime }}& z_i^{{}^{\prime }}\\ ...& ...& ...\end{array}\right],b=\left[\begin{array}{c}\Delta I_1\\ ...\\ \Delta I_i\\ ...\end{array}\right]$$
显然采用最小二乘的解法
$$e=(A^{T}A{)}^{-1}{A}^{T}b=[I_x,I_y,I_z{]}^T$$
其实在实际应用中，需要将上式计算的e投影到点p的表面上

其中点颜色的变化表示该点的$intensity$的变化，红色表示很高，黑色表示很低。蓝线表示上式计算出来的$e$,黑线表示该点的单位法向量$n$,$n=[n_x,n_y,n_z{]}^T$，绿线表示在该表面投影后的向量$e^{’}$
$$e^{’}=e-n(n^{T}e)=e-n(e^{T}n)$$
最后将$e^{’}$代替上式的$e$，这样可以降低噪声的灵感度。

最后一个问题，Tomasi Response怎么确定，是取M矩阵的最小特征值吗？答案并不是

通过上图可知，公路上车道线的终点也是一个特征点，但在平面上$intensity$是按照$x,y$的方向分布的，$z$方向分布为零，即平面上的变化也能产生特征点，因次Tomasi Response $R={\lambda }_2$并不是${\lambda }_3$。

---

2.2 Harris3D Without Intensity

点云中的某点$p$的的曲面方程为
$$f(x,y,z)=0$$
同样地，我们构建了一个代价函数
$$E(u,v,w)=\sum _{x,y,z\in \Omega }[f(x+u,y+v,z+w)-f(x,y,z){]}^2$$
点$p$的表面会变成一个平面

此时平面方程为
$$ax+by+cz+d=0$$
此时点$p$的法向量$n=[n_x,n_y,n_z{]}^T$,$a=n_x,b=n_y,c=n_z$。

那么
$$f(x+u,y+v,z+w)=a(x+u)+b(y+v)+c(z+w)+d=a{x}^{’}+b{y}^{’}+c{z}^{’}+d$$
其实
$$dist(point,plane)=\frac{a{x}^{\prime }+b{y}^{\prime }+c{z}^{\prime }+d}{\sqrt{a^2+b^2+c^2}}=\frac{a{x}^{\prime }+b{y}^{\prime }+c{z}^{\prime }+d}{1}$$
$f(x+u,y+v,z+w)$可以直观地理解为，点$p=[x,y,z]$在移动$[u,v,w]$后，距离以点$p$构建的平面的距离。

那么，结合代价函数和$f(x,y,z)$表示，可以将代价函数$E(u,v,w)$变换为
$$\begin{gathered} E(u,v,w)=\sum _{x,y,z\in \Omega }[f(x+u,y+v,z+w)-f(x,y,z){]}^2 \\ f(x,y,z)=ax+by+cz+d \\ f(x+u,y+v,z+w)=a(x+u)+b(y+v)+c(z+w)+d \end{gathered}$$
推导出
$$E(u,v,w)=\left[\begin{array}{ccc}u& v& w\end{array}\right]M\left[\begin{array}{c}u\\ v\\ w\end{array}\right]=\left[\begin{array}{ccc}u& v& w\end{array}\right]\sum _{x,y,z\in \Omega }\left[\begin{array}{ccc}{n}_x^2& n_x{n}_y& n_x{n}_z\\ n_x{n}_y& n_y^2& n_y{n}_z\\ n_x{n}_z& n_y{n}_z& n_z^2\end{array}\right]\left[\begin{array}{c}u\\ v\\ w\end{array}\right]$$
此时$M$是曲面法线在曲面上的协方差矩阵，同样的此时Tomasi Response $R={\lambda }_3$

---

2.3 扩展 Harris 6D

结合上述有$intensity$和无$intensity$的$M$矩阵，其中$I=[I_x,I_y,I_z,n_x,n_y,n_z]$得
$$M=\sum _{x,y,z\in \Omega }\left[\begin{array}{cccccc}{I}_x^2& I_x{I}_y& I_x{I}_z& I_x{n}_x& I_x{n}_y& I_x{n}_z\\ I_x{I}_y& I_y^2& I_y{I}_z& I_y{n}_x& I_y{n}_y& I_y{n}_z\\ I_x{I}_z& I_y{I}_z& I_z^2& I_z{n}_x& I_z{n}_y& I_z{n}_z\\ n_x{I}_x& n_x{I}_y& n_x{I}_z& n_x^2& n_x{n}_y& n_x{n}_z\\ n_y{I}_x& n_y{I}_y& n_y{I}_z& n_x{n}_y& n_y^2& n_y{n}_z\\ n_z{I}_x& n_z{I}_y& n_z{I}_z& n_x{n}_z& n_y{n}_z& n_z^2\end{array}\right]$$
此时Tomasi Response $R={\lambda }_4$

因为$R={\lambda }_3$表示$Harris6D$是{$Harris3D$无$intensity$}的子集，此时并没有{$Harris3D$有$intensity$}的信息，条件太松懈；$R={\lambda }_5$表示一个特征点必须同时是一个几何点($Harris3D$ 无$intensity$)和强度点($Harris3D$ 有$intensity$),条件太苛刻。故$R={\lambda }_4$。

3.总结

	输入	协方差矩阵	响应函数R
$Harris$ $Image$	Image	像素梯度	$R={\lambda }_2$
$Harris3D$	PC with Intensity	Intensity 梯度	$R={\lambda }_2$
$Harris3D$	PC without Intensity	法向量	$R={\lambda }_3$
$Harris6D$	PC with Intensity	Intensity 梯度+法向量	$R={\lambda }_4$