目前网上能找到的所有关于“使用积分图计算点云法相”的资料都是来自PCL的,只有code,没有原理
这篇笔记来自论文 Adaptive neighborhood selection for real-time surface normal estimation from organized point cloud data using integral images
积分图定义
将一张图片中的每一个像素点的取值替换为其(在原始图片中)左上角区域的所有像素值之和,就是积分图。公式如下:
I
O
(
m
,
n
)
=
∑
i
=
0
m
∑
j
=
0
n
O
(
i
,
j
)
\mathcal{I}_{\mathcal{O}}(m, n)=\sum_{i=0}^{m} \sum_{j=0}^{n} \mathcal{O}(i, j)
IO(m,n)=i=0∑mj=0∑nO(i,j)
也可通过迭代的方式计算,这样只用从左上角开始遍历一次图片就可以获得它的积分图:
I
O
(
m
,
n
)
=
O
(
m
,
n
)
+
I
O
(
m
−
1
,
n
)
+
I
O
(
m
,
n
−
1
)
−
I
O
(
m
−
1
,
n
−
1
)
\begin{aligned} \mathcal{I}_{\mathcal{O}}(m, n)=& \mathcal{O}(m, n)+\mathcal{I}_{\mathcal{O}}(m-1, n) +\mathcal{I}_{\mathcal{O}}(m, n-1) -\mathcal{I}_{\mathcal{O}}(m-1, n-1) \end{aligned}
IO(m,n)=O(m,n)+IO(m−1,n)+IO(m,n−1)−IO(m−1,n−1)
计算一个方形区域的平均值可以用如下公式:
S
(
I
O
,
m
,
n
,
r
)
=
1
4
r
2
⋅
(
I
O
(
m
+
r
,
n
+
r
)
−
I
O
(
m
−
r
,
n
+
r
)
−
I
O
(
m
+
r
,
n
−
r
)
+
I
O
(
m
−
r
,
n
−
r
)
)
\begin{aligned} S\left(\mathcal{I}_{\mathcal{O}}, m, n, r\right)=\frac{1}{4 r^{2}} \cdot (& \mathcal{I}_{\mathcal{O}}(m+r, n+r) \\ &-\mathcal{I}_{\mathcal{O}}(m-r, n+r) \\ &-\mathcal{I}_{\mathcal{O}}(m+r, n-r) \\ &+\mathcal{I}_{\mathcal{O}}(m-r, n-r)) \end{aligned}
S(IO,m,n,r)=4r21⋅(IO(m+r,n+r)−IO(m−r,n+r)−IO(m+r,n−r)+IO(m−r,n−r))
其中
(
m
,
n
)
⊤
(m, n)^{\top}
(m,n)⊤表示方形区域的中心坐标,
r
r
r为边长的一半。这里计算方形区域是为了后面用的,当然也可以计算任意矩形区域,将
r
r
r换为
a
,
b
a,b
a,b做相应的改变就好啦
参考资料:积分图
平滑处理(预处理)
由于RGB-D相机采集的数据(即深度图)会有很多噪声,常用平滑处理来降噪。
平滑处理的一个很重要的参数是平滑窗口的大小。
在深度图中,不能使用固定大小的平滑窗口,因为1)不同深度的地方,信噪比不一样;2)边缘区域不能用大的平滑窗口,否则边缘信息都被抹掉了。
所以,需要考虑上述两个因素,分别计算每一个地方的平滑窗口大小,计算细节如下:
- 首先,用公式表示深度值
d
d
d与最小深度变化(minimum depth change)的关系:
f
D
C
(
d
)
=
α
⋅
d
2
f_{D C}(d)=\alpha \cdot d^{2}
fDC(d)=α⋅d2也就是下图中的红线
实验结果表明, α = 0.0028 \alpha=0.0028 α=0.0028 - 计算Depth-Dependent Smoothing Area Map: B ( m , n ) = β ⋅ f D C ( D ( m , n ) ) \mathcal{B}(m, n)=\beta \cdot f_{D C}(\mathcal{D}(m, n)) B(m,n)=β⋅fDC(D(m,n))其中, β \beta β由用户设定,用来控制平滑窗口的大小, D ( m , n ) \mathcal{D}(m, n) D(m,n)为深度图上点 ( m , n ) ⊤ (m, n)^{\top} (m,n)⊤处的值(即深度)
- 计算Depth Change Indication Map: C ( m , n ) = { 1 { if ∥ δ D x ( m , n ) ∥ ≥ t D C ( D ( m , n ) ) or ∥ δ D y ( m , n ) ∥ ≥ t D C ( D ( m , n ) ) 0 otherwise. \mathcal{C}(m, n)=\left\{\begin{array}{l} 1 \quad\left\{\begin{array}{l} \text { if }\left\|\delta \mathcal{D}_{x}(m, n)\right\| \geq t_{D C}(\mathcal{D}(m, n)) \\ \text { or }\left\|\delta \mathcal{D}_{y}(m, n)\right\| \geq t_{D C}(\mathcal{D}(m, n)) \end{array}\right. \\ 0 \text { otherwise. } \end{array}\right. C(m,n)=⎩⎨⎧1{ if ∥δDx(m,n)∥≥tDC(D(m,n)) or ∥δDy(m,n)∥≥tDC(D(m,n))0 otherwise. 其中 δ D x ( m , n ) = D ( m + 1 , n ) − D ( m , n ) , δ D y ( m , n ) = D ( m , n + 1 ) − D ( m , n ) \delta \mathcal{D}_{x}(m, n)=\mathcal{D}(m+1, n)-\mathcal{D}(m, n) , \delta \mathcal{D}_{y}(m, n)=\mathcal{D}(m, n+1)-\mathcal{D}(m, n) δDx(m,n)=D(m+1,n)−D(m,n),δDy(m,n)=D(m,n+1)−D(m,n)分别表示在x、y正方向的深度变化, t D C ( d ) = γ ⋅ f D C ( d ) t_{D C}(d)=\gamma \cdot f_{D C}(d) tDC(d)=γ⋅fDC(d)表示深度变化检测的阈值,通过 γ \gamma γ来控制用户想要的检测灵敏度
- 计算Final Smoothing Area Map:也就是每个位置的滑动窗口大小 R ( m , n ) = min ( B ( m , n ) , T ( m , n ) 2 ) \mathcal{R}(m, n)=\min \left(\mathcal{B}(m, n), \frac{\mathcal{T}(m, n)}{\sqrt{2}}\right) R(m,n)=min(B(m,n),2T(m,n))其中, T ( m , n ) \mathcal{T}(m, n) T(m,n)为距离变换算法(Distance Transformations)计算获得的值,表示点 ( m , n ) ⊤ (m, n)^{\top} (m,n)⊤的到下一个depth change的2D距离; 2 \sqrt{2} 2是因为这里用的是方形区域,而非圆。
方法1:只使用一张积分图
最简单的计算法向量的方法是将上下左右的点(红色)连线叉乘,Kinect中用的就是这种方法。公式为:
n
⃗
p
=
v
⃗
p
,
h
×
v
⃗
p
,
v
(1)
\vec{n}_{p}=\vec{v}_{p, h} \times \vec{v}_{p, v} \tag{1}
np=vp,h×vp,v(1)
但是这样计算太粗糙了,很容易受到噪声影响。可以使用单张积分图来计算公式(1)中的两个向量
v
⃗
p
,
h
,
x
=
P
x
(
m
+
r
,
n
)
−
P
x
(
m
−
r
,
n
)
2
v
⃗
p
,
h
,
y
=
P
y
(
m
+
r
,
n
)
−
P
y
(
m
−
r
,
n
)
2
v
⃗
p
,
h
,
z
=
S
(
I
P
z
,
m
+
1
,
n
,
r
−
1
)
2
−
S
(
I
P
z
,
m
−
1
,
n
,
r
−
1
)
2
v
⃗
p
,
v
,
x
=
P
x
(
m
,
n
+
r
)
−
P
x
(
m
,
n
−
r
)
2
v
⃗
p
,
v
,
y
=
P
y
(
m
,
n
+
r
)
−
P
y
(
m
,
n
−
r
)
2
v
⃗
p
,
v
,
z
=
S
(
I
P
z
,
m
,
n
+
1
,
r
−
1
)
2
−
S
(
I
P
z
,
m
,
n
−
1
,
r
−
1
)
2
\begin{aligned} \vec{v}_{p, h, x} &=\frac{\mathcal{P}_{x}(m+r, n)-\mathcal{P}_{x}(m-r, n)}{2} \\ \vec{v}_{p, h, y} &=\frac{\mathcal{P}_{y}(m+r, n)-\mathcal{P}_{y}(m-r, n)}{2} \\ \vec{v}_{p, h, z} &=\frac{S\left(\mathcal{I}_{\mathcal{P}_{z}}, m+1, n, r-1\right)}{2} -\frac{S\left(\mathcal{I}_{\mathcal{P}_{z}}, m-1, n, r-1\right)}{2} \\ \vec{v}_{p, v, x} &=\frac{\mathcal{P}_{x}(m, n+r)-\mathcal{P}_{x}(m, n-r)}{2} \\ \vec{v}_{p, v, y} &=\frac{\mathcal{P}_{y}(m, n+r)-\mathcal{P}_{y}(m, n-r)}{2} \\ \vec{v}_{p, v, z} &= \frac{S\left(\mathcal{I}_{\mathcal{P}_{z}}, m, n+1, r-1\right)}{2} -\frac{S\left(\mathcal{I}_{\mathcal{P}_{z}}, m, n-1, r-1\right)}{2} \end{aligned}
vp,h,xvp,h,yvp,h,zvp,v,xvp,v,yvp,v,z=2Px(m+r,n)−Px(m−r,n)=2Py(m+r,n)−Py(m−r,n)=2S(IPz,m+1,n,r−1)−2S(IPz,m−1,n,r−1)=2Px(m,n+r)−Px(m,n−r)=2Py(m,n+r)−Py(m,n−r)=2S(IPz,m,n+1,r−1)−2S(IPz,m,n−1,r−1)其中
P
x
,
P
y
,
P
z
\mathcal{P}_{x},\mathcal{P}_{y},\mathcal{P}_{z}
Px,Py,Pz分别表示图像上的二维点经内参矩阵变换到三维点的每一个通道(x,y,z)组成的“图片”;
I
P
z
\mathcal{I}_{\mathcal{P}_{z}}
IPz表示z对应的“图片”的积分图;
r
=
R
(
m
,
n
)
r=\mathcal{R}(m, n)
r=R(m,n)。
然后就可以用公式(1)计算法向量了。
方法2:使用多张积分图
另一种计算法向量的思路是计算领域的协方差,然后求得特征值最小的特征向量作为法向量。领域的选择可以为最近领域搜索or一定半径内的范围,但是太费计算量了。可以使用九张积分图来计算协方差然后获得法向量:
- 先计算三张积分图: I P x , I P y , I P z \mathcal{I}_{\mathcal{P}_{x}},\mathcal{I}_{\mathcal{P}_{y}},\mathcal{I}_{\mathcal{P}_{z}} IPx,IPy,IPz,与方法1中倒数第二段中的一样,分别为图像上的二维点经内参矩阵变换到三维点的每一个通道(x,y,z)组成的“图片”的积分图
- 然后计算六张积分图: I P x x , I P x y , I P x z , I P y y , I P y z , I P z z \mathcal{I}_{\mathcal{P}_{x x}}, \mathcal{I}_{\mathcal{P}_{x y}},\mathcal{I}_{\mathcal{P}_{x z}}, \mathcal{I}_{\mathcal{P}_{y y}}, \mathcal{I}_{\mathcal{P}_{y z}}, \mathcal{I}_{\mathcal{P}_{z z}} IPxx,IPxy,IPxz,IPyy,IPyz,IPzz。其中 I P a b \mathcal{I}_{\mathcal{P}_{a b}} IPab表示 I P a , I P b \mathcal{I}_{\mathcal{P}_{a}},\mathcal{I}_{\mathcal{P}_{b}} IPa,IPb的点积,即对应位置相乘(类似向量的点乘)
- 点 ( m , n ) ⊤ (m, n)^{\top} (m,n)⊤处的协方差公式为 C p = ( c x x c x y c x z c y x c y y c y z c z x c z y c z z ) − ( c x c y c z ) ( c x c y c z ) ⊤ \mathcal{C}_{p}=\left(\begin{array}{ccc} c_{x x} & c_{x y} & c_{x z} \\ c_{y x} & c_{y y} & c_{y z} \\ c_{z x} & c_{z y} & c_{z z} \end{array}\right)-\left(\begin{array}{c} c_{x} \\ c_{y} \\ c_{z} \end{array}\right)\left(\begin{array}{c} c_{x} \\ c_{y} \\ c_{z} \end{array}\right)^{\top} Cp=⎝⎛cxxcyxczxcxycyyczycxzcyzczz⎠⎞−⎝⎛cxcycz⎠⎞⎝⎛cxcycz⎠⎞⊤其中 c x x = S ( I P x x , m , n , R ( m , n ) ) c x y = c y x = S ( I P x y , m , n , R ( m , n ) ) c x z = c z x = S ( I P x z , m , n , R ( m , n ) ) c y y = S ( I P y y , m , n , R ( m , n ) ) c y z = c z y = S ( I P y z , m , n , R ( m , n ) ) c z z = S ( I P z z , m , n , R ( m , n ) ) c x = S ( I P x , m , n , R ( m , n ) ) c y = S ( I P y , m , n , R ( m , n ) ) c z = S ( I P z , m , n , R ( m , n ) ) \begin{aligned} c_{x x} &=\mathcal{S}\left(\mathcal{I}_{\mathcal{P}_{x x}}, m, n, \mathcal{R}(m, n)\right) \\ c_{x y}=c_{y x} &=\mathcal{S}\left(\mathcal{I}_{\mathcal{P}_{x y}}, m, n, \mathcal{R}(m, n)\right) \\ c_{x z}=c_{z x} &=\mathcal{S}\left(\mathcal{I}_{\mathcal{P}_{x z}}, m, n, \mathcal{R}(m, n)\right) \\ c_{y y} &=\mathcal{S}\left(\mathcal{I}_{\mathcal{P}_{y y}}, m, n, \mathcal{R}(m, n)\right) \\ c_{y z}=c_{z y} &=\mathcal{S}\left(\mathcal{I}_{\mathcal{P}_{y z}}, m, n, \mathcal{R}(m, n)\right) \\ c_{z z} &=\mathcal{S}\left(\mathcal{I}_{\mathcal{P}_{z z}}, m, n, \mathcal{R}(m, n)\right)\\ c_{x} &=\mathcal{S}\left(\mathcal{I}_{\mathcal{P}_{x}}, m, n, \mathcal{R}(m, n)\right) \\ c_{y} &=\mathcal{S}\left(\mathcal{I}_{\mathcal{P}_{y}}, m, n, \mathcal{R}(m, n)\right) \\ c_{z} &=\mathcal{S}\left(\mathcal{I}_{\mathcal{P}_{z}}, m, n, \mathcal{R}(m, n)\right) \end{aligned} cxxcxy=cyxcxz=czxcyycyz=czyczzcxcycz=S(IPxx,m,n,R(m,n))=S(IPxy,m,n,R(m,n))=S(IPxz,m,n,R(m,n))=S(IPyy,m,n,R(m,n))=S(IPyz,m,n,R(m,n))=S(IPzz,m,n,R(m,n))=S(IPx,m,n,R(m,n))=S(IPy,m,n,R(m,n))=S(IPz,m,n,R(m,n))
- 最后,计算 C p \mathcal{C}_{p} Cp的最小特征值对应的特征向量就可以作为法向量
虽然多个积分图的计算量大于单张积分图的,但是它可以直接获得领域的平面性 (planarity),可用于平面拟合、边缘/角点检测
实验
- 实验平台:2.26 GHz Intel® Core™2 Quad CPU, 4 GB of RAM
- 速度:方法1(SDC) VS 方法2(CM) VS KNN(也就是将方法2中领域用KNN找,再计算协方差矩阵,最后求得特征向量)
可以看到SDC的时间最少,而且CM和SDC的复杂度(几乎)为O(1) - 误差:
可以看到SDC的误差最少
PCL中的积分图计算法向量方法
- COVARIANCE_MATRIX:从邻域的协方差矩阵中创建9个积分图去计算一个点的法线
- AVERAGE_3D_GRADIENT:创建了6个积分图去计算3D梯度里面垂直和水平方向的光滑部分,同时利用两个梯度的叉积来计算法线。
- AVERAGE_DEPTH_CHANGE:只创建了一个单一的积分图,从平均深度的变化中计算法线。
- SIMPLE_3D_GRADIENT:
无任何介绍
参考资料:
PCL官方文档中IntegralImageNormalEstimation模板类里的NormalEstimationMethod枚举