定义问题
给定一组3D点 {pi} ,我们想找到这组3D点满足的平面参数,即平面的法向量 n 和中心
n⊺(pi−q)=0,∀i(1)
当然在实际数据中(1)式一般是不可能严格满足的,因此可以定义如下函数:
dist(pi;n,q)≜n⊺(pi−q)(2)
上式表示点到面的符号距离,即就是这个距离可能为正也可能为负
注意: n 表示平面的法向量,并且为单位向量:
nTn=1
求解问题
现在这个问题可以转化为一个最小二乘问题:
cost(n,q)≜∑idist2(pi;n,q)=∑i(n⊺(pi−q))2=n⊺[⋯,pi−q,⋯][⋯,p⊺i−q⊺,⋯]⊺n=n⊺A(q)A(q)⊺n(3)
whereA(q)≜[⋯,pi−q,⋯]
第一步先固定平面法向量求 q ,则(3)对
0=∂cost(n,q)∂q≡∑i(2nn⊺q−2nn⊺pi).(4)
求解(4)式可以得到最优的平面中心 q∗ :
q∗=1|{pi}|∑ipi.(5)
这个结果和我们直观理解一致,并且这个结果与平面法向量的取值无关。
或者可以从另外一个角度去理解,直观使用(5)式就可以求出一组3D点的中心,(4)式从另外一个角度验证了(5)式的正确性
第二步求解平面的法向量 n , 我们可以将(5)式带入(3)式中得如下结果:
其中 B(q∗)≜A(q∗)A(q∗)⊺ ,由于 A(q) 是一个 3×|{pi}| 维的矩阵,所以 B(q) 是一个 3×3 维的正定矩阵,现在问题可以表示为如下形式:
n∗=argminns.t.n⊺B(q∗)nn⊺n=1(7)
上式的形式和PCA的目标函数非常相似:
maxa1s.t.a⊺1∑a1a⊺1a1=1(8)
唯一的区别就在于PCA求的是目标函数的最大值,而我们的问题求的是最小值。所以只需要对 B(q∗) 进行分解,最小特征值对应的特征向量就是平面的法向量。当然我们也可以从直观角度去理解,对于一组3D的点平面的法向量一定是最不重要的那个投影向量(理论情况下所有平面上的点在法向量上的投影为零)。
代码
使用PCA拟合平面是我在看冯晨的《Fast Plane Extraction in Organized Point Clouds Using Agglomerative Hierarchical Clustering》论文时看到的,大家有兴趣的可以看看,代码地址,运行结果如下图所示
PCA
下面简单说说主成分分析(PCA),主要是是为了备忘。
记 x1,...,xp 为p个原始特征,设新的特征 ξi,i=1,...,p 是原始特征的线性组合
ξi=∑j=1p=aijxj=aTix(9)
为了统一 ξi 的尺度,不妨要求线性组合系数的模为1,即:
aTiai=1(10)
将(9)式写成矩阵的形式:
ξ=ATx(11)
其中, ξ 是由新特征 ξi 组成的向量, A 是特征变换矩阵。要求解的是最优的正交变换矩阵
考虑第一个新特征 ξi :
ξ1=∑j=1p=a1jxj=aT1x(12)
它的方差为
var(ξ1)=E[ξ21]−E2[ξ1]=E[aT1xxTa1]−E[aT1x]E[aT1x]=aT1E[xxT]a1−aT1E[x]E[xT]a1=aT1∑a1(13)
其中, ∑ 是 x 的协方差矩阵, 可以用样本来估计,
f(a1)=aT1∑a1−v(aT1a1−1)(14)
其中, v 是拉格朗日乘子, 将(14)式对
∑a1=va1(15)
这是协方差矩阵 ∑ 的特征值, 即 a1 是矩阵 ∑ 的特征向量, v 是对应的特征值。 把(15)式带入(13)式中,可得:
因此,最后的 a1 是矩阵 ∑ 最大特征值对应的特征向量。 ξ1 称为第一主成分, 它在原始特征的所有线性组合中方差是最大的。
协方差矩阵 ∑ 共有 p 个特征值
∑i=1pvar(ξi)=∑i=1pλi(17)
变换矩阵 A 的各个列向量是由
x=Aξ(18)
通常把主成分零均值化,即
ξ=AT(x−μ)x=Aξ+μ(20)
这种平移并不影响主成分的方向。
参考资料
[1] Chen Feng, Fast Plane Extraction in Organized Point Clouds Using Agglomerative Hierarchical Clustering
[2] 主成分分析(PCA)原理详解
[3] 深入理解拉格朗日乘子法(Lagrange Multiplier) 和KKT条件
[4] 张学工, 模式识别(第三版), 主成分分析, 163-164