基于排序的网格截面算法应用

  本文讨论的范围不包含几何重建、纹理重建，仅包含基于已有的人体三维三角网格模型(下称人体三角网格)的三维测量算法。围度测量实际上是对三维模型与特定平面求交，并对交集进一步分析其交集的特定闭合部分(若存在多个闭合部分)的轮廓周长和面积等数据，即求三维网格的切片。国标GB/T 16160-2008《服装用人体测量的部位及方法》中给出了人体各尺寸特征的与人身高的比例数据(下称特征比例)，求出特定比例处与身高方向构成的点法式平面与三维网格的切片即可进一步得到特定特征的具体数值。

1.1 人体三角网格的坐标系矫正

  假定人体三角网格的坐标系中Z轴方向为人角的中点与头顶中心构成的方向向量，而手持扫描仪获取到的人体三角网格的坐标系取决于首次扫描的朝向，因此网格的Z轴方向与正确的坐标系有误差，需要通过额外的工作进行坐标轴的矫正。

  本文采用的核心方法为PCA(Principal Component Analysis)。对于手臂与身体躯干的面元粘连问题，可通过扫描目标张开手臂进行规避，而此时扫描目标人站立并张开手臂时的高低肩会影响PCA坐标矫正的效果。影响PCA算法的主要部分为手臂部分三维点，去除手臂部分三维点进行二次PCA算法即可。简单地可以通过特定的比例关系选取，亦可通过三维点到Z轴的距离选取。

1.2 三角网格切片算法

1.2.1 基于排序的三角网格切片算法

  基于三角网格的切片算法可分为两类。一是基于半边数据结构的三角网格切片算法，假定半边结构半边对象为：

struct HE_edge{
	HE_vert* from;
	HE_vert* to;
	HE_edge* pair;   // oppositely oriented adjacent half-edge 
	HE_edge* next;   // next half-edge around the face
	...   			  // other members
};

并且此时获取到了平面plane与指定片面相交的起始半边s, 可知算法的结果为从s出发的一个半边序列,那么算法可以简单描述为HE_Slice(s->pair,*)，通过循环去递归可进一步优化:

void HE_Slice(HE_edge* curr, std::vector<glm::vec3>& res){
	if(curr->next == s || curr->next->next == s) return ;
	for(auto e: {e->next, e->next->next}) {
		if(is_intersect(e,plane)) {
			... // calculate and store the intersection
			HE_Slice(e->pair, res);
		}
	}
}

1.2.2 基于排序的三角网格切片算法

  基于半边数据结构的切片算法是高效的，但半边数据结构的维持需要较多的内存支持且构建耗时_O(n*log(n))_、只提供了局部拓扑信息不能支持按身高比例高效获取特定半边。上述两个缺点均可通过下面介绍的面元排序解决。考虑到人体三角网格与平面plane的交集可以为多个闭合部分(手与躯干、双脚),半边结构能且只能获取到给定的初始半边所在闭合部分，不能获得指定平面与三角网格的交集的闭合部分的数量(下称闭合数)。算法框架如下:

1. 对共享顶点结构的三角网格进行排序
2. 通过二分查找确定需要与平面求交的有序面元的最小子序列
3. 遍历子序列，生成平面与三角形面元的边的交的集合,集合元素为{A,B,k},其中A,B为边的顶点索引，k为相交系数。
4. 通过选择排序对3.生成的集合完成有序化，并得到闭合部分的划分。

  参与切片算法的最小几何元素为边，因此对三角网格排序应当等价于对三角形面元排序。排序算法简单描述如下:

// 1. 对共享顶点结构的三角网格进行排序
void sortByVector(TriMesh *mesh, glm::vec3 x) {
    std::sort(mesh->f.begin(), mesh->f.end(), [=](glm::ivec3& e1, glm::ivec3& e2){
        return glm::dot(mesh->v[e1[0]] - mesh->v[e2[0]], x) > 0;
    });
}

  易知实际中参与切片算法的面元是排序后的面元序列的一个子序列，当给出切片平面(点法式平面，法向量为Z轴)和面元在Z轴上的最大投影距离_gap_时，子序列提取算法如下：

// 2. 通过二分查找确定需要与平面求交的面元的最小有序子序列
std::array<std::vector<glm::ivec3>::iterator,2> getSliceInterval(TriMesh* mesh,glm::vec3 n, float d, float gap) {
    std::array<std::vector<glm::ivec3>::iterator,2> result = {
        std::lower_bound(mesh->f.begin(),mesh->f.end(), d-gap,
            [=](glm::ivec3& e1,float v){
        return v < glm::dot(mesh->v[e1[0]], n);
    }), std::lower_bound(mesh->f.begin(),mesh->f.end(), d+gap,
            [=](glm::ivec3& e1,float v){
        return v < glm::dot(mesh->v[e1[0]], n);
    })};

    if(result[0] > result[1]) std::swap(result[0],result[1]);
    return result;
}

  确定子序列的边界，切片算法便只需要极少的面元与特定平面求交。面元与平面的相交判定与交线等几何计算如下:

bool isFaceInersected(TriMesh* mesh,glm::ivec3 f, glm::vec3 n, float d){
    int flags = ((glm::dot( mesh->v[f[0]], n ) > d)+(glm::dot( mesh->v[f[1]], n ) > d)+(glm::dot( mesh->v[f[2]], n ) > d));
    return flags ==  1 || flags == 2;
}

std::array<glm::vec3,2> getFaceIntersection(TriMesh* mesh,glm::ivec3 f, glm::vec3 n, float d){
    std::array<glm::vec3,2> result;
	int size = 0;
    // dot(N, P1) + dot(N, t*(P1-P2) ) = d  // 0 < t < 1
    for(int i:{0,1,2})  {
        double numerator   =  d - glm::dot( mesh->v[f[i]], n);
        double denominator = glm::dot( mesh->v[f[i==2?0:i+1]] - mesh->v[f[i]], n );
        if(abs(denominator) > 1e-8) {
            numerator /= denominator;
            if(numerator >= 0 && numerator <= 1.0) {
                result[size++] = glm::vec3(f[i],f[i==2?0:i+1],numerator);
            }
        }
    }
    return result;
}

  3.得到的线段是无序的，不能够满足后续操作，如闭合数、特定闭合部分的平滑、特定部分的周长面积等几何特征的计算，故需要对3.的结果有序化。考虑到实际网格规模下的交线的规模往往低于10000,算法设计采用了选择排序，算法复杂度为_O(n*log(n))_,简单描述如下：

//4. 通过选择排序对3.生成的集合完成有序化，并得到闭合部分的划分。
std::vector<std::array< std::vector<std::array<glm::vec3,2>>::iterator,2>> sortContours(std::vector<std::array<glm::vec3,2>>&& intersections) {
    std::vector<std::array< std::vector<std::array<glm::vec3,2>>::iterator,2>> intervals = {{intersections.begin(),intersections.begin()}};
    for(auto it = intersections.begin(); it != intersections.end(); it++) {
        auto find_it = std::find_if(std::next(it), intersections.end(), [=](std::array<glm::vec3,2>& e){
            return  (e[0][0] == (*it)[1][1] && e[0][1] == (*it)[1][0]) || (e[0][0] == (*it)[0][1] && e[0][1] == (*it)[0][0])
                    ||(e[1][0] == (*it)[0][1] && e[1][1] == (*it)[0][0]) || (e[1][0] == (*it)[1][1] && e[1][1] == (*it)[1][0]);
        });
        if(find_it != intersections.end()) {
            std::swap(*std::next(it),*find_it);
        }else{
            (*std::prev(intervals.end()))[1] = std::next(it);
            if(it != std::prev(intersections.end()))
                intervals.push_back({std::next(it),std::next(it)});
        }
    }
    return intervals;
}	

  到此，基于排序的切片算法基本原理和实现已介绍完毕。具体的代码实现可见本人git仓库yubaMesh