说明

最近需要使用c++实现L1-median骨架提取算法^[1]，所以干脆直接开一个系列记录一下算法具体实现中的一些要点的理论分析和代码实现。

本篇分析的是算法中针对采样点(sample point)设定的一个概念——有向度(directionality degree)。在论文提到，带权重的PCA方法被用于提取sample point的有向度，那么有向度的具体计算和PCA的关系为何？具体是如何使用PCA的呢？

PCA方法简述

例子来自于郑申海的知乎文章^[2]：
有m个2维数据如下：

待处理数据

需要使用PCA方法投影成m个1维数据，理论上需要遵循最大方差原则，使得投影的一维值之间区分度最高。接下来补充几个前置知识：

• 协方差矩阵
  
  协方差矩阵
  
  协方差矩阵具有奇妙的特性：对角线上为每一维数据的方差，其它元素则是对应的协方差，当只有一维时，显然此矩阵只有对角线上有值，这种状态被称为对角化。

• 特征值分解矩阵
  
  特征值分解矩阵
  
  根据特征值定义可以写出如图公式，为特征向量矩阵，C为被求特诊向量的矩阵，最后结果为对角化的特征值矩阵。

PCA的过程可以概括为：对于m个n维数据矩阵，我们只需要找一个使得之后映射的有一维，也就是让的协方差矩阵对角化。而和存在如下关系：

结合特征值分解矩阵的形式我们可以得知，使等同于的特诊向量矩阵的转置，则可以保证对角化。根据此思路，PCA过程可以描述为如下过程:

对于m个n维数据矩阵：
1) 的每一行减去均值（零均值化）
2) 求的协方差矩阵及特诊向量
3) 将特征向量根据特征值的大小降序排列成矩阵
4) 根据计算映射后的

引申出来对三维数据的协方差矩阵：

可进行应用来获取三维数据的分离特征

为什么选择PCA

在PCA的执行过程中，我们可以明显看到协方差矩阵具有对我们十分有意义的特性——单维数据的协方差矩阵对角化。

• directionality degree定义
  
  L1-median原文中定义

在有向度的定义中使用的是的类协方差矩阵格式，具有如下性质：

代表到的矢量，每个组成的矩阵只有一个非0特征值
因而，如果neighborhood size内只有一个其它点的话，
，
所以，当neighborhood size内的点方向越整齐，每个（实际上是，此处为简写）加和组成的矩阵的特征值就越单一，最终的就越接近1

代码思路

对于特征向量和特征根的计算使用了Eigen库，具体内容如下：

double GlobalFun::computeDirectionalityDegree(vector<PointXYZ> diff) {
    if (diff.size() < 1) return 0;
    Matrix3d mm;
    for (PointXYZ d : diff) {
        Vector3d v(d.x, d.y, d.z);
        mm += v * v.adjoint();
    }
    Vector3cd eigens = mm.eigenvalues();
    Vector3d eigenVec(eigens(0).real(), eigens(1).real(), eigens(2).real());
    return eigenVec.maxCoeff() / eigens.sum().real();
}

输入的diff按照定义是由集合中的邻居计算而来，每次邻居的更新使用了pcl库提供的kdtree进行快速更新：

// 根据radius计算X中每个点的neighbors
vector<SamplePoint> updateSelfNeighbors(PcPtr xc, double radius, vector<SamplePoint> &info) {
    KdTreeFLANN<PointXYZ> kdtree;
    kdtree.setInputCloud(xc);
    int index = 0;
    for (PointXYZ xpt : xc->points) {
        vector<int> neighIndex;
        vector<float> neighDistance;
        int num = kdtree.radiusSearch(xpt, radius, neighIndex, neighDistance);
        info[index].setSelfIndics(neighIndex);
        info[index].setSelfDistances(neighDistance);
        index++;
    }
    return info;
}
// 根据neighbors计算diff，然后调用computeDirectionalityDegree计算有向度
double SamplePoint::computeSigma(pi::PcPtr cloud) {
    // 根据s_indic计算
    vector<PointXYZ> delta;
    for (int i : this->s_indics) {
        delta.push_back(PointXYZ(
            this->pos.x - cloud->points[i].x,
            this->pos.y - cloud->points[i].y,
            this->pos.z - cloud->points[i].z
        ));
    }
    this->sigma = GlobalFun::computeDirectionalityDegree(delta);
    
    return this->sigma;
}

总结

其实虽然论文中说的是带权PCA，实际使用的更多的是协方差矩阵的数学性质，到了代码层面重点更是变成了特征值求取和带权协方差矩阵构建。可能稍微有些 “挂羊头，卖狗肉” 的意味吧。

1. Huanghui, Wushihao, Cohenordaniel, et al. L1-medial skeleton of point cloud[J]. ACM Transactions on Graphics, 2013. ↩

2. https://zhuanlan.zhihu.com/p/21580949 ↩