定义问题
给定一组3D点{pi},我们想找到这组3D点满足的平面参数,即平面的法向量n和中心q :
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)对q求偏导并且令导数为零:
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)式中得如下结果:
cost(n;q∗)≜n⊺A(q∗)A(q∗)⊺n=n⊺B(q∗)n(6)
其中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 是特征变换矩阵。要求解的是最优的正交变换矩阵A,它使得新特征ξi 的方差达到极致。正交变换保证了新特征之间不相关,而新特征的方差越大,则样本在该维特征上的差异就越大,因而这一特征就越重要。
考虑第一个新特征ξ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 的协方差矩阵, 可以用样本来估计, E[.] 是数学期望。要在(10)式约束下最大化ξ1 的方差,这等价于求下面拉格朗日函数的极值:
f(a1)=aT1∑a1−v(aT1a1−1)(14)
其中,v 是拉格朗日乘子, 将(14)式对a1 求导并令其为零,得到最优解a1 满足
∑a1=va1(15)
这是协方差矩阵∑ 的特征值, 即a1 是矩阵∑ 的特征向量, v 是对应的特征值。 把(15)式带入(13)式中,可得:
var(ξ1)=aT1∑a1=vaT1a1=v(16)
因此,最后的a1 是矩阵∑ 最大特征值对应的特征向量。 ξ1 称为第一主成分, 它在原始特征的所有线性组合中方差是最大的。
协方差矩阵∑ 共有 p 个特征值λi,i=1,...,p (包括可能相等的特征值和可能为零的特征值),把他们降序排列。按照上面相同的方式,可以得出有对应这些特征向量构造的p 个主成分ξi,i=i,...,p ,全部主成分的方差之和为:
∑i=1pvar(ξi)=∑i=1pλi(17)
变换矩阵A 的各个列向量是由∑ 的正交特征向量组成的, 因此AT=A−1 , 即A 是正交矩阵。 从ξ 到x 的逆变换为:
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