几何建模与处理-极小曲面局部方法

欢迎关注公众号:sumsmile /专注图形渲染、Android的开发老兵~~

内容主要参考中科大刘利刚老师的课程《几何建模与处理基础》, 代码参考谢聪同学的作业

运行效果

说几句题外话

几何处理涉及的数学知识比较多，微分几何、离散微分几何、拓扑学、数值处理、凸优化等等，入门课很多学校放在数学专业来讲，深入的课程比如微分流行是放在研究生阶段讲，要掌握好，需要一定的高数基础。

为了更好的理解《几何处理》这门课，我中途停下来去学习《微分几何》、《离散微分几何》，所以中断四五个月没有发文。

课程链接：

《微分几何》-西北工业大学 陈航老师

www.bilibili.com/video/BV1g7…
这门课使用参考书微分几何-彭家贵

《Discrete Differential Geometry》-卡耐基梅隆大学的图形大佬 Keenan Crane

www.youtube.com/playlist?li… (B站上也有)

这两门课都讲的很好，而且都有完整的ppt提供。有兴趣的朋友，建议先从《微分几何》开始，Keenan的课是英文没有翻译，我花了三倍的时间逐句翻译才勉强看完。

进入正题~~

曲率

曲率是几何处理最重要的概念之一。微分几何、离散微分几何课程至少有三分之一的内容都在讲曲率的计算与引用。

极小曲面的实现中，核心的逻辑也要依赖曲率的计算。

平均曲率是个向量值，很像梯度的概念，你站在山上沿着下山最快的方向跑，那个方向就是平均曲率的方向，平均曲率的值就是下山的速度。

曲率的定义：

曲率公式的推导和变换要点时间，初次接触，建议只要了解公式的使用即可，上面提到的《微分几何》课程里有详细的曲率推导过程。

什么是极小曲面

严格的定义：平均曲率处处为0的曲面

看上面这张图，极小曲面也是非常光滑的曲面

这篇文章实现将一个立体形状逐步变形成极小曲面。实现的效果参考文章开头的三段视频。

极小曲面原理

先从一般形式入手

迭代的基本原理非常简单：

1. 多个点${Q_k}$和 ${P_{old}}$相连，求点集${Q_k}$的平均值${Q_{mean}}$，向量$\overrightarrow{PQ_{mean}}$方向就是P点移动的方向
2. 方向上乘以系数${\lambda}$，以保证每次是微小的迭代，不至于出现很突兀的、奇异的变形
3. 每个点都按照P点的逻辑移动，一次迭代，整个图形就会有微小的平滑
4. 重复1 2 3步骤，进行多次迭代，直至收敛，整个曲面外观达到理想的平滑状态

和机器学习里算法训练很像，都是沿着最快收敛方向迭代，每次一小步，逐步收敛。如曾国藩说：日拱一卒，功不唐捐。

换成平均曲率

按照周围点的平均值来计算是不严谨的，因为周围点的几何属性没有体现，换言之，每个周围点的权重是不一样的。

合理的方式是，按照平均曲率的方向来计算移动的方向，平均曲率方向更贴近点P的凸起方向。

和普通形式类似，只是把L换成Hn，Hn即平均曲率，是个向量值。

Hn如何算呢？用余切Laplace-Beltrami和平均曲率公式可以求出Hn和几何网格的关系

1. 根据余切 Laplace-Beltrami 算子的定义:

2. 离散平均曲率的计算(离散平均曲率可以由连续形式的平局曲率推导出来)

划重点！！ 这个公式是极小曲面计算的灵魂！！前面的都可以不记得，记住这个就可以开干了！

再次友情提示！这里涉及的公式，能理解多少就算多少，初次接触学会使用即可，不必强求能理解背后的数学推导。

• $\bigtriangledown_p$ : Lalpace-Beltrami算子,详细的推导有点复杂，暂不展开
• A：点P附近小块区域的面积，可以理解为和点P相连的三角形的面积，A的计算比较重要，下文展开讲

代码实现

• 开发语言、环境：C++ CLion
• 依赖库：ImGUI pmp-library openGL
• 代码(-from 谢聪)，上传到github方便大家下载：

github.com/summer-go/G…

使用开源库

使用pmp-library库来处理网格计算，我们只需要关注处理流程

常见的几何网格处理库,选一个合适的即可

• CGAL:www.cgal.org
• Libigl:github.com/libigl/libi…
• MeshLab:www.meshlab.net
• PCL(Point Cloud Library):www.pointclouds.org
• TriMesh:graphics.stanford.edu/software/tr…
• DGtal:dgtal.org

操作说明

pmp-library集成了ImGui，实现了基础的界面交互，按空格键可切换渲染模式

最小曲面迭代计算分隐式计算(解析方程)、显示计算(逐步迭代)，最下面计算高斯曲率和平均曲率的开关

voronoi_area计算

平均曲率涉及邻域面积"A"的计算 有很多种方法，这里介绍两种

1. Mixed Voronoi

即在Voronoi的基础上，针对钝角做些调整，将外心点替换为对边的中心，详细可看下面这种图中的对比说明

2. Barycentric + Voronoi

即根据三角形的特点，分别采用Barycentric 或 Voronoi，区分判断条件如下：

1. 顶点和边夹角位钝角，使用Barycentric cell，即取各边的中点围成的面积，即该小网格的$\frac{1}{2}$
2. 对应的三角网格各个夹角均为锐角的，采用Voronoi Cell计算 $A_{M_{jk}}=\frac{1}{8}(cot\beta*|\mathbf{x}_i-\mathbf{x}_j|+cot\alpha*|\mathbf{x}_i-\mathbf{x}_k|)$
3. 注意,如果临边的夹角为钝角，取该小网格面积的$\frac{1}{4}$，看起来像是减小该网格的权重

本文采用第二种，对应的代码在Common.cpp

double voronoi_area(const SurfaceMesh& mesh, Vertex v)

隐式迭代

和显式迭代的逻辑相反，假设${P_{new}} - \lambda * \bigtriangledown_X$ 可以反向得到${P_{old}}$

建立方程组，一次求出所有${P_{new}}$

$I * \overrightarrow{P_{new}} - \lambda * \bigtriangledown_X = \overrightarrow{P_{old}}$

${I}$是单位矩阵，$\bigtriangledown_X$中也有${P_{new}}$因子 最后可以把P抽出来，得到方程租

$（I - \lambda * \bigtriangledown_{cot-weight}）* \overrightarrow{P_{old}} = \overrightarrow{P_{}}$

参考资料

[1]微分几何-彭家贵: 

www.wenqujingdian.com/Public/edit…

[2]pmp-library: 

github.com/pmp-library…

欢迎关注公众号:sumsmile /专注图形渲染、Android的开发老兵~~