基于高斯曲率的mesh扁平化

背景

在Geometry Processing中,我们有时候需要对一个mesh进行参数化。对于一个开放的含边界的mesh,我们一般会采用planar parameterization。可以想象将一层附于物体表面的“皮”摊平在一个2D泡面。但对于一个genus-0的封闭mesh,直接展平不可避免的会出现许多的畸变和扭曲。所以我们考虑将其参数化到同胚的单位球上。
一般来说,复杂的多面体无法通过stereographic projection(球极平面投影)直接投射到单位球上。换个角度理解,即无法从几何体的几何中心无遮挡的看到三角面上的所有顶点。三角面在单位球上的投影会overlap或者flip。这都会导致参数化的失败。一般问题的发生场所都是在mesh的突出部分,所以在本文中,将会利用curvature-driven faltenning技术完成对mesh不改变拓扑的扁平化。从而能更好的进行参数化。

曲率

曲率本身是一个二阶连续的表面的二阶导数。但是对于3D模型来说,它是离散的且只有$C^0$连续。所以我们需要对其进行一个分段的线性近似。

Gaussian曲率的定义

如果我们定义$p_i, p, p_{i+1}$构成的夹角为$\alpha$,那么点p处的angular defect就为
$$\delta =2\pi -\sum _{i} \alpha _{i}$$
然后高斯曲率K为
\[K=\frac{2\pi -\sum\limits _{i} \alpha _{i}}{G}\]
其中G为几何项,上图中的蓝色部分
\[G=\sum\limits _{i} A _{i}( p) /3\]
最后化简为
\[
K=\frac{3\left( 2\pi -\sum\limits _{i} \alpha _{i}\right)}{\sum\limits _{i} A _{i}( p)}
\]
其中$A _{i}(p)$是领接与点p的三角形$Tri_i$的面积

算法思路

首先计算出在mesh上的每个顶点的高斯曲率$K_p$。然后我们挑出曲率最大的一个顶点$p_{max}$,即我们看到的最突出的地方。我们计算能让它的曲率变小的目标位置$p’_{max}$
\[
p’ _ {max}=\frac{\sum _ {p_{i} \in Neigh( p_{max})} p_{i}}{val( p_ {max})}
\]
其中$Neigh( p_{max})$代表$p_{max}$的所有相邻点。$val( p_ {max})$是点$p_{max}$的度。上面的式子其实就代表了相邻顶点的几何中心(或者说算术平均位置)
如果$p_{max}$的目标位置和$p_{max}$的原始位置的距离超过了阈值$dist$,俺么我们就移动$p_{max}$到$p’ _{max}$。否则,我们就跳过$p _{max}$,继续找下一个曲率最大的点,尝试更新位置。该操作一直会持续到找到一个点位置。
当成功更新了一个顶点的位置之后,我们就会更新相应的高斯曲率并开始下一次迭代。

值得一提的是,离散高斯曲率的计算有不同形式,但是这种表示方式是一种local的表示,也即意味着如果一个顶点的位置发生了改版,那我们只用更新它相关的顶点并重新计算即可,具有locality。

不断重复上述的过程,最终的效果是将mesh中突出的部分渐渐打磨成圆角。理想的最终结果是一个接近于球的形状。
不过上面式子中的$K_p$在遇到某个小区域内有大量顶点的时候(通常是模型表现细节的地方),由于每个三角形的面积趋向于0,所以会对迭代算法造成不利的影响,收敛速度会大大减慢。为了解决这个问题,于是Mocanu[3]提出了一个系数$\chi _{s}$
\[
\chi _ {s} =\frac{\sum\limits _ {i} A( Tri_{i})}{n_{Tri}}
\]
这个$\chi _{s}$代表了整个mesh的所有三角面的平均大小。这个纠正因子可以通过加强angular defect的影响而改善迭代的效果。最终修改过的高斯曲率$K_s$为
\[
K_{s} =\frac{2\pi -\sum\limits _ {i} \alpha _ {i}}{\chi _ {s} +\sum\limits _ {i} A_ {i}( p) /3}
\]

这个步骤会不断进行。在结束的时候,会用一些PCA和大小归一化的方法来避免产生萎缩的问题

算法实现细节

Step1

首先先对mesh的全部三角面计算一次平均面积。在此过程中,三角形的面积会cache,这样后面只需更新那些移动顶点位置了的三角的cahcedArea就可以进行进一步迭代了。体现了local操作的性能优势

1
2
3
4
5
6
7
float average_face_area = 0.0;
for (auto f = mesh.facesBegin(); f != mesh.facesEnd(); f++) {
f->cachedArea = f->area();
average_face_area += f->cachedArea;
}
// Get Xs
average_face_area /= mesh.nFaces();

Step2

接着对每个顶点生成一个曲率的记录,主要记录当前曲率和曲率减小的目标点的位置。计算时会需要用到刚才的平均面积

1
2
3
4
5
MutablePriorityQueue<VertexRecord> queue;
for (auto v = mesh.verticesBegin(); v != mesh.verticesEnd(); v++) {
v->record = VertexRecord(v, average_face_area);
queue.insert(v->record);
}

Step3

不断的找出最大曲率的点并尝试更新。
对于那些虽然曲率很大,但是移动距离没有能达到阈值要求的顶点,我们会暂时将他们的曲率设为0,并重新插入优先队列。这样他们就不会在本次迭代中再次被访问和判断。只有在已经成功选出一个更新点并更新位置之后,它们的曲率才会再次得到计算和更新。

1
2
3
4
5
6
7
while (queue.top().offset < threshold) {
VertexRecord maxCurvedV = queue.top();
queue.pop();
// Make vR the least "curved" vert and insert back to queue.
maxCurvedV.curvature = 0.0;
queue.insert(maxCurvedV);
}

Step4

对于满足条件,需要更新的顶点。在完成位置移动之后,我们首先要更新它相邻三角面的面积。同时注意,此时mesh的平均面积也发生了变化,可是我们可以通过increment变量来局部的更新平均面积而不用整个再重新计算一遍。

1
2
3
4
5
6
7
8
9
10
11
// Move p to pMax
maxCurvedV.vert->position = maxCurvedV.optimalPoint;
// Update neighbor faces
auto adjFs = maxCurvedV.vert->AdjFaces();
for (auto adjF : adjFs) {
float newArea = adjF->area();
// Update Xs
average_face_area += (newArea - adjF->cachedArea) / mesh.nFaces();
// Update face area before we use it to calculate our latest curvature
adjF->cachedArea = adjF->area();
}

Step5

最后,别忘了将刚才更新过的顶点remove并re-insert回去,以便下一次能找到曲率最大的顶点。

1
2
3
4
5
6
7
8
9
10
11
12
13
// Update current vert's curvature
maxCurvedV.Update(average_face_area);

// As well as those neighbor verts
auto adjVerts = maxCurvedV.vert->AdjVertices();
for (auto adjV : adjVerts) {
queue.remove(adjV->record);
adjV->record.Update(average_face_area);
queue.insert(adjV->record);
}

// after updated maxCurvedV get re-insert to the queue to get the right order
queue.insert(maxCurvedV);

引用

[1] Discrete Differential-Geometry Operators for Triangulated 2-Manifolds [Meyer et al. 2003]
[2] Compute Harmonic weights on a triangular mesh
[3] Direct Spherical Parameterization of 3D Triangular Meshes Using Local Flattening Operations