The day will be what you make it,
so rise, like the sun, and burn
—— William C. Hannan
最近接触到的基于体素的地形渲染方案频次较高，其中频繁提到的一个概念或者问题是，如何将体素转化为网格以得到平滑的渲染效果，基于这个问题，这里做了一些调研与学习，并输出了这篇文章。
1. 体素网格化方案列表
基于对材料的学习，大致整理出如下的一些网格体素化方案：
- Marching Cube
方案优点：简单、高效
方案不足：效果的二义性与voxel内部特征模拟的精度损失 - Dual Contouring
- Naive Surface Net
- Marching Tetrahedra
下面对上面的几个方案做详细的介绍，希望能够从实现原理、优劣分析等方向给出一些个人见解。
2. Marching Cube
Marching Cube可以算得上是体素网格化方案中最普遍使用，最广为人知的了，大概的实现原理可以理解为，对于每个voxel，其上有8个顶点，基于每个顶点上的数值，我们可以将顶点划分为两个阵营，之后使用若干个面（最多4个，最少0个）来将这两个阵营区分开来，从而实现voxel到面片的转换。
上面介绍的是大致的原理，其中还有许多细节是被掩盖的，需要逐一澄清：
- 每个顶点的数值怎么得到，通常我们只有voxel本身的数据（如voxel的填充率），怎么得到voxel上8个顶点的数据？
- 如何完成对顶点阵营的划分，哪些顶点适合划分到同一个阵营，哪些不适合？
- 对于8个顶点共2^8=256种划分情况来说，我们如何计算出对应的划分面？另外，对于同样的阵营划分结果来说，如果顶点上的数值存在一些差别，最终输出的划分面数据是否是相同的？
针对上面的问题，下面我们来逐一分析澄清。
2.1 顶点数据获取
如何从voxel数据中获取到顶点的数据，这个实际上会取决于对应的应用场景，有的场景就能够直接输出对应顶点的一个数值，有的则是只给了voxel的数值，对于前者，我们就已经有了数据，对于后者来说，我们可以通过对voxel数据进行处理来得到顶点上的数据。
处理方法可以采用最简单直观的，比如对于某个顶点，会有8个相邻的voxel数据，那么可以直接对这8个voxel的数值进行平均，就得到了顶点数据
2.2 顶点阵营划分
在得到了各个顶点的数值之后，下一个要面临的问题在于，我们如何知道哪些顶点需要分配在同一阵营呢？
通常的做法是选择一个目标值，小于这个数值的顶点归于第一阵营，大于等于这个数值的顶点归于第二阵营。
voxel是离散数据，顶点也是，在顶点之间的广袤空间中还分布着无数的点，我们通常会将具有相同数值的点连接组成的曲面称之为等值面，而我们接下来要做的就是基于刚才的目标值来计算此值对应的等值面，即前面说过的划分面，当然这里我们是简化计算的，因为我们只考虑平面情况，不考虑曲面情况，且平面情况中，我们考虑的也是肯定会与voxel的各个面存在相交线的平面情况（这里的平面是有边界的平面，而非数学中的无限平面）。
2.3 划分平面求解
这一块是Marching Cube的核心，在其他关于Marching Cube的文章中，也主要是围绕这个话题展开。
我们这里需要解决的问题在于，对于8个顶点（已知数值），在给定了目标值的情况下，我们要如何求取出voxel上的等值面。
前面说了，针对8个顶点的排列组合，我们共有256种阵营划分方法，Marching Cube这边采取的策略是，将划分结果直接存入到一个尺寸为256的数组中，之后将8个顶点的阵营关系对应到一个uint8的对应bit上，从而得到一个索引值，根据索引值可以拿出对应的划分结果。
上图展示了去重的15种等值面的排布关系，其他的241种划分结果都可以通过对这15种排布关系进行旋转来得到。
可以看到，这里的等值面，实际上是经过如下操作得到的：
- 在voxel的12条边上，根据其两端顶点的阵营的情况，我们进行处理，选定阵营不同的边，在边上放置一个点（图中绘制的是在中点，但实际情况并不一定在中点），这个点肯定是等值面与此边的交点
- 通过对交点进行连线，就得到了等值面
这里就产生了新的问题：
- 交点的位置是如何确定的？
- 哪些交点应该通过线段连在一起？
对于第一个问题，我们可以通过插值公式来计算得到，举个例子，对于某条边AB来说，A顶点的数值为0，B顶点的数值为1，等值面的数值为0.5，那么交点肯定就在AB的中点，如果等值面的数值为0.3，那么交点就在0.3AB的位置，也就是说，通过调整等值面的数值，我们是可以得到不同的交点乃至等值面结果的，而等值面的值可以在项目里根据需要进行设置，比如，如果没有特殊需要，emtpy voxel的数值为-1，full voxel的数值为1，那么此数值可以设定为中间值0，插值公式根据比例关系计算，给出如下：
对于第二个问题，这里为了保证结果的一致性，以避免voxel等值面输出后的连接断裂问题，这里是直接从已经设定好的面片数组中获取的（查表法），比如对一个如下的顶点阵营分布情况，我们得到的就是如下图对应的面片分布：
但实际上，这里还存在其他的连线方式，比如在上述连线的基础上，添加左边面片上两个交点的连线与右侧面片上两个交点的连线，从而得到不同的等值面结果（且这种异构等值面还不止一种），而实际上我们并不能确定实际情况中哪一种等值面是最能准确反映出真实的voxel的数值的：
换成2D视角，可能更容易理解：
这就引出了Marching Cube的一个问题，即等值面分布的二义性，即对于相同的一套顶点阵营分布，输出的等值面结果并不唯一，且我们也无法确定哪种分布是最为合理的。而将这个效果应用到3D地形上，则会导致一个问题就是孔洞，即在某个voxel上经过处理后在右边面片上有一条边，但是在右边的voxel上处理结果在左边的面片上却缺失了这条边（不知道如果通过对查表中的等值面分布进行特殊设计，是否可以确保两个voxel的等值面能正确衔接从而消除这种问题？）
另外，由于marching cube采用查表法实现的等值面都是横跨或者说完整切割掉voxel的，不能模拟出等值面终止于voxel内部的情况，因此会导致一些voxel内部的棱角信息的丢失，如下图所示：
下图是真实希望达到的效果，而上图则是Marching Cube算法的结果，可以看到，右侧的一个角被直接切掉，用一条边来代替了，效果上存在一定的损失。
3. Dual Contouring
针对前面Marching Cube算法所面临的一些问题，有人提出了一种叫做Dual Contouring的算法。
下面介绍一下这种算法的大致实现思路。
给定一个用于表达原始模型的函数f(x)（每个函数都可以理解为一个等值面）以及其对应的梯度f'(x)（梯度在3D中就表示各个点上的法线），我们也有一个与marching cube算法中类似的voxel grid。
对于与f(x)相交的所有voxel，都可以在其中放置一个顶点，那么我们就得到了若干个顶点。
对于voxel上与f(x)相交的边（边相邻的四个voxel肯定与f(x)相交），我们可以将相邻voxel中的四个顶点按照某种顺序连接起来，从而得到两个三角形（这两个三角形可能共面）。
经过上面处理之后得到的三角形list就是对于voxel网格化的结果，用2D来表示，示意如下：
红色连线就相当于voxel的网格化结果，可以看到，在Dual Contouring方法中，如果用voxel中的顶点来代表voxel的话，那么可以认为voxel就成为了网格上的顶点，而顶点与voxel本身是互偶关系，这就是dual contouring方法名字的由来。
上面介绍了大致试下思路，这里会有一些问题需要解决：
- 如何从f(x)与f'(x)中得到各个voxel中的顶点
- 实际计算中的一些异常情况要如何处理
3.1 顶点计算
对于给定等值面（或者说函数）而言，我们可以直接知道等值面与voxel是否相交（判定过程跟marching cube类似），以2D为例，假设我们得到了某个quad（对应于3D的voxel）与等值面的相交情况（如下图的绿色点）：
我们希望的结果是，我们拟合出来的模型表面在相交点上除了函数值要契合，其导数（法线）也要契合，从一个直观的想象中，我们在quad中放置的顶点就需要满足到各条边上交点的连线能够与交点上的法线相垂直，即需要满足，对于任何相交点p_i（法线n_i)有：
对于一些只知道函数f(x)不知道法线f'(x)的情况而言，可以通过微分来计算出这个数值：
对于2D而言，如果只与quad的两条边相交，那么我们总能找到这样一个顶点，但是如果同时与quad的四条边相交，那么可能就不一定能找到这样一个点了，这个情况在下面问题处理一节中有介绍。
在求得顶点之后，将相邻quad中的顶点直接相连就得到了对应的contouring线，而对于3D而言，我们就得到了对应的2个三角形：
3.2 问题处理
- 基于垂直法则进行顶点求取，如果因为存在多条约束无法得到最优解要怎么办？
这个对应的就是前面介绍的，函数与quad的相交边数多于2，可能导致任意两个交点的最优顶点不重合，如下图所示：
在这种情况下，采取的解决方法是找不到最优解，那就寻找极优解，前面说过，我们希望顶点到各个交点处的连线与交点的法线是垂直的，那么问题就转换为，对于每个交点，我们可以为quad中的每个点施加一个惩罚函数，最后只需要将每个交点的惩罚函数加在一起求取在quad中的最小值即可，如下图所示：
最终，我们将问题转化为求解二次误差函数（QEF，quadratic error function）的极小值求解：
这是一个较为常见的极小值求解问题，可以通过QR矩阵分解来求取，具体可以参考[3]，在一些语言如python中提供了直接求解的函数库，可以省去自己实现的时间消耗[4]。
- 求取出来的顶点落在quad之外怎么办？
基于前面的极值求取方法，在某些情况下，如两个交点的法线接近平行（即等值面是一个平面）时，会导致极值点落在quad之外，如下图所示：
针对这种问题，有两种解决方案：
- 第一种是对QEF的解添加约束，将之限定在quad内部
- 第二种则是对QEF函数做修正，比如添加一个统一的偏移，使得这个函数的解永远落在quad内部，通常这个偏移函数需要设计成对于原始解在quad内部的情况施加极低或者不施加影响，而只对原始解落在quad外的情况添加一个往quad中心拖拽的偏移
具体可以参考[4]中的一些实现。
模型或者函数是自相交的怎么办？
模型存在面片穿插的情况，会导致Dual Contouring的结果存在异常，这个问题的解决方案可以参考 “Intersection-free Contouring on An Octree Grid” Ju and Udeshi 2006的实现。voxel尺寸过大导致模型的一些细节无法被展示怎么办？
这里描述的是我们在voxel中通过一个顶点无法处理的情况，比如某个模型在voxel中存在两个不想交的切割面，这时候如果只在voxel中放置一个顶点显然就不能还原这种效果，缩小voxel尺寸可以解决，但是成本会比较高，这个问题称之为manifold问题：
解决这个问题需要对Dual Contouring方法做一些扩展，具体可以参考这个github实现 或者Manifold Dual Contouring
4. Naive Surface Net
Naive Surface Net在生成模型的面数上要远远低于Marching Cube，这里有个Demo可以做一下参考。
这个方法的原理就是，在前面marching cube中，我们拿到了等值面与每个voxel上对应边的一个交点，那么对于每个voxel，我们可以计算出这些边的交点的重心（交点坐标算数平均），之后将这个作为Dual Contouring算法中的顶点，完成后续的网格化操作，下面看下效果对比：
这个方法的优点在于：
- 比Dual Contouring方法简单
- 输出的面数要低于Marching Cube方法，计算消耗也要更低
- 实现起来比较简单
5. Marching Tetrahedra
Marching Tetrahedra方案的实现思路跟Marching Cube可以说是如出一辙，不同的是Marching Cube是针对cube与等值面的256种相交情况输出对应的mesh结果，而Marching Tetrahedra则是针对正四面体与等值面相交的16（2^4)种情况进行查表得到结果，下面看下对应的表现：
相对于Marching Cube，这个方案有如下的一些好处：
- 这个方法没有二义性
- 每个四面体输出的面数要小于cube上输出的面数，使得这个方法在GS中更为合适。
- 相交情况去除重复的话只有三种，考虑重复也只有16种（相对于marching cube的15 & 256种），内存占用更少
它的不足在于，得到的总面数大约是marching cube的4x之多，在实现性能上也略有不如，这也是为什么marching cube更为通用的重要原因。
贴图
前面介绍的都是如何得到平滑的拓扑网格，在实际开发中，我们除了网格之外，还需要为之添加对应的材质（或者说贴图效果），假设我们为给定voxel的六个面（其实是三个方向）赋予了不同的材质效果，那么在进行网格化后的效果要如何计算呢？
这里[6]采用的方法是Triplanar Mapping方案，即对于每个像素，计算其在三个方向上对于材质的采样结果，之后基于这一点的世界法线在三个方向上的投影进行加权混合，比如世界法线与Z轴（向上）平行，那么就只取上表面的采样结果即可。
上面给的是单个voxel中面片的着色，如果考虑到整个网格，可能还需要在不同材质之间进行混合，这里通常就跟地形的材质混合一样，为每个voxel添加多个材质，并在材质之前添加一个渐变权重即可。
参考文献
[1]. 水泡动画模拟(Marching Cubes)
[2]. 网格生成之marching cube算法学习笔记
[3]. 三维等值面提取算法（Dual Contouring）
[4]. Dual Contouring Tutorial
[5]. Smooth Voxel Terrain (Part 2)
[6]. Smooth Voxel Mapping: a Technical Deep Dive on Real-time Surface Nets and Texturing