今天要介绍的是Jad Khoury在其硕士论文中介绍的使用CS完成物件、地形网格的渐进式Tessellation方案，这个方案目前被业界很多的引擎所借鉴引用，具有较高知名度。

1.1 Introduction

GPU光栅化在三角面片在屏幕上的投影面积较大（覆盖像素较多）时效率会高一些，当其投影面积小于某个阈值时，就会导致深度buffer的锯齿以及shading rate的明显下降[Riccio 12]，而这会使得复杂场景的渲染存在很多的挑战，因为当物件距离相机越来越远，就会导致三角面片在屏幕上的投影面积变得越来越小，本文给出的算法可以为任意的polygon mesh实现渐进式的精修，从而避免上述问题。（为什么渐进式精修可以避免问题？）

传统的模型渐进式精修是在CPU上借助不断迭代的算法如四叉树或者subdivision surface完成的，但是这种算法的问题在于要想完成渲染需要将数据从CPU传输到GPU，而这个成本是比较高的，因此大家想是不是可以直接在GPU上通过tessellation shader完成这个精修处理，然而tessellation shader存在着如下的两个限制：

1. 只支持最高级的细分
2. 随着细分深度的增加，性能会随之下降（具体原因暂时不明）[AMD 13]

本文给出的算法也是在GPU上完成的，不过可以避免上面tessellation shader的相关问题，而且即使使用再多的细分层级，其内存消耗都会维持在一个稳定水平。整个过程是在Compute Shader中完成的，采用的是一种基于三角面片的implicit subdivision scheme，数据的读取跟写入都是在一个紧凑的双缓冲（double buffered）数组中完成。

1.2 Implicit Triangle Subdivision

1.2.1 Subdivision Rule

面片细分算法的基础是面片拆分规则（subdivision rule），拆分规则描述的是单个面片如何被拆解成多个子面片。本文使用的是一种叫做二分三角形（binary triangle）的拆分规则（简称BTSR），具体可以参考下图（a）中的示意图：

如图所示，BTSR可以将直角等腰三角形拆分成两个直角等腰三角形，我们用序号0和1表示，在图中有对应的位置示意，这两个子三角形的质心空间变换矩阵（barycentric-space transformation matrices）给出如下：

这里需要解释下什么是质心空间变换矩阵，我这边没有搜到比较好的匹配答案，通过摸索推测出，通过这个变换矩阵可以将Cartesian坐标系（笛卡尔坐标系）下的Parent Triangle的三个顶点坐标（分别用(0, 0, 1)， (0, 1, 1)，(1, 0, 1)表示，以直角所在顶点为横纵坐标轴的交点）转换成Child Triangle的三个顶点的笛卡尔坐标（比如1号子面片的三个顶点坐标（与前面Parent Triangle的三个顶点一一对应）分为(1/2, 1/2, 1)，(0, 0, 1)，(1, 0,1)；而0号子面片的三个顶点则为(1/2, 1/2, 1)，(0, 1, 1)，(0, 0,1)）

mat3 bitToXform (in uint bit )
{
  float s = float ( bit ) - 0.5;
  vec3 c1 = vec3 ( s, -0.5, 0);
  vec3 c2 = vec3 ( -0.5 , -s, 0);
  vec3 c3 = vec3 (+0.5 , +0.5 , 1);
  return mat3 (c1 , c2 , c3);
}

上面的伪代码给出了使用GLSL是如何根据输入的序号构建对应的的，从前面的分析可以看到，第N次拆分，就会生成个新的三角面片。

1.2.2 Implicit Representation

根据前面的描述，经过细分后的三角面片都可以通过多次细分时的选择来表达，比如每次细分中选择划分线左边用0划分线右边用1表示，那么如前面图中(b)，(c)小图所示，每个三角形都可以用一个编码来表示，这个编码我们称之为key，其中前者表示的是完全划分（uniform），而后者表示的是部分划分（adaptive），每个三角形三个顶点的坐标则可以通过多次划分是的质心空间坐标转换矩阵连乘来表示，比如key = 0100的三角形的变换矩阵就可以用如下方法得到：

key是通过一个uint32来表示的，比如0100的key的二进制表示方法为：

可以看到前面有一个多出来的1，实际上为了知道整个key从哪里开始，需要添加一个起始位，这个是固定的，再往前的下划线（通常应该填充0）对应的是无关的位数。

mat3 keyToXform (in uint key )
{
  mat3 xf = mat3 (1) ;
  while ( key > 1u) 
  {
    xf = bitToXform ( key & 1u) * xf;
    key = key >> 1u;
  }
  return xf;
}

上面的GLSL代码片段给出了对应key的三角面片的变换矩阵计算方法。

单独一个uint32算上根节点level最多只能表达32-1个划分level，要想得到更多的划分层级，需要使用一个vector，比如使用一个uvec2可以提供63个level的划分。

1.2.3 Iterative Construction

面片拆分过程需要通过递归来完成，但是GPU天生是不支持递归的，因此需要一种策略来解决这个问题，原文中给出的算法是使用一个叫做subdivision buffer的double-buffers，通过ping-pong方法从BufferA读取写入BufferB，再从BufferB读取写入BufferA，循环往复直到所有的key bit都被处理完了，这个过程是在Compute Shader中完成的（CS线程是如何分配的？到后面sub triangle数目变多之后，一轮处理是不是不足以覆盖所有的triangle？list中包含多个二进制编码的三角形，之后每个线程直接从中取用即可？）。

每次CS调用可以允许实现三种处理：

1. 计算出下一级的相关数据
2. 回退到上一级的数据
3. 维持原样

因为用key作为面片表达方法的原因，这些操作都可以轻松完成，下面这张图给出了一个key以及其对应的child triangle和parent triangle的key的数据表示：

Parent Key跟Child Key的计算逻辑可以用如下的GLSL代码表示，因为是些基本的移位操作以及逻辑操作，所以成本很低：

uint parentKey (in uint key )
{
  return ( key >> 1u);
}
void childrenKeys (in uint key , out uint children [2])
{
  children [0] = ( key << 1u) | 0u;
  children [1] = ( key << 1u) | 1u;
}

下面给出的是在CS中对subdivision buffer进行更新的GLSL代码，简单来说，如果某个key需要被拆分，那就会生成两个Child Key并塞入到subdivision buffer中（ping-pong），同时抛弃当前key；而如果两个Child Key需要合并成一个Parent Key，这时候就会丢掉这两个Key，并将Parent Key塞入Buffer，为了避免某个key被重复塞入，通常会选择在0号Child上进行相应处理。

buffer keyBufferOut { uvec2 u_SubdBufferOut []; };
uniform atomic_uint u_SubdBufferCounter ;
// write a key to the subdivision buffer
void writeKey ( uint key )
{
  uint idx = atomicCounterIncrement ( u_SubdBufferCounter );
  u_SubdBufferOut [ idx ] = key ;
}

// general routine to update the subdivision buffer
void updateSubdBuffer ( uint key , int targetLod )
{
  // extract subdivision level associated to the key
  int keyLod = findMSB ( key );
  // update the key accordingly
  if (/* subdivide ? */ keyLod < targetLod && ! isLeafKey ( key )) 
  {
    uint children [2]; childrenKeys (key , children );
    writeKey ( children [0]) ;
    writeKey ( children [1]) ;
  } 
  else if (/* keep ? */ keyLod == targetLod ) 
  {
    writeKey ( key);
  }
   else /* merge ? */ 
  {
    if (/* is root ? */ isRootKey ( key )) 
    {
      writeKey ( key );
    }
    else if (/* is zero child ? */ isChildZeroKey ( key)) 
    {
      writeKey ( parentKey (key));
    }
  }
}

下面这段代码用于检测当前节点是否属于其Parent的0号子节点

bool isChildZeroKey (in uint key ) { return (key & 1u == 0u); }

下面这段代码用于判断对应节点是叶子节点还是根节点

bool isRootKey (in uint key ) { return ( key == 1u); }
bool isLeafKey (in uint key ) { return findMSB (key) == 31; }

从现有的信息来看，当前给出的算法可以很好的适应GPU的特性，因此使得我们可以实现此前示意图中所展示的部分划分（Adaptive）方案。 不过需要注意的是，由于每轮迭代，我们只能对每个key执行一次细化或者粗化处理，如果我们需要进行更多次处理的话，就需要依赖多buffer的迭代方案，在原方案的渲染实现中，会在每帧开始之前就进行一次单buffer的迭代（这是因为至少需要调用一次迭代吗？）

1.2.4 Conversion to Explicit Geometry

经过前面的处理，我们拿到了代表不同尺寸与细分结果的keys，下一步要做的就是将key转换成渲染所使用的三角面片数据。这里是通过GPU Instancing来完成的。

对于subdivision buffer中的每个key，这里都会创建出一个对应的三角面片，每个面片的顶点坐标可以通过下面给出的GLSL代码来计算，顶点上的其他属性比如法线等也可以通过类似机制计算出来，出于阐述简单考虑，这里就没有展示了：

// Get 3D Coordinates According To Vector Weights
vec3 berp (in vec3 v[3] , in vec2 u)
{
  return v [0] + u.x * (v [1] - v [0]) + u.y * (v [2] - v [0]) ;
}
// subdivision routine ( vertex position only )
void subd (in uint key , in vec3 v_in [3] , out vec3 v_out [3])
{
  mat3 xf = keyToXform ( key);
  vec2 u1 = (xf * vec3 (0, 0, 1)).xy;
  vec2 u2 = (xf * vec3 (1, 0, 1)).xy;
  vec2 u3 = (xf * vec3 (0, 1, 1)).xy;
  v_out [0] = berp (v_in , u1);
  v_out [1] = berp (v_in , u2);
  v_out [2] = berp (v_in , u3);
}

1.3 Adaptive Subdivision on the GPU

1.3.1 Overview

前面介绍了implicit subdivision scheme，下面我们将介绍如何使用这个方案来对场景中的物件进行曲面细分。总的来说，就是为每个面片创建一个对应的adaptive subdivision，从而控制其在屏幕空间的规模，避免文章开头说到的sub-pixel projection问题（指的是通过将多个sub-pixel面片合并成一个大的面片，从而避免这个问题吧？）

上图给出了本文算法在OpenGL中的实现流程，可以看到这是一个递归调用的过程，每一帧执行一次调用，完成后再将subdBufferIn跟subdBufferOut互换。其中绿色表示的是GPU中的buffer数据，红色表示GPU中的stage，而灰色表示的是CPU中的stage。可以看到，GPU中的执行stage主要有三个（使用OpenGL 4.5版本），这三个都是在Compute Shader中完成：

1. LoadKernel，这个stage会使用之前介绍过的implicit subdivision scheme对CS中的subdivision buffer进行更新，完成曲面细分或者合并；此外，为了减轻后面的工作量，还会触发一个frustum culling调用，使用一个原子计数器对subdivision buffer中的key进行裁剪，只将可见的keys填入到CulledSubdBuffer中。
2. IndirectBatcherkernel，这个stage的目的是为下一次LoadKernel的调用做准备，比如完成相应数据（DispatchIndirectBuffer）的准备工作，同时也完成下一Stage的数据（DrawIndirectBuffer）准备工作。
3. RenderKernel，调用IndirectDraw指令完成geometry数据到framebuffer的渲染过程。这个过程会创建一系列的triangle instance（InstancedGeometryBuffers），前面CulledSubdBuffer中的每个key分别对应一个instance。

1.3.2 LOD Function

前面的implicit subdivision scheme介绍了如何细分与合并，这里来介绍在什么情况下面片需要细分，什么情况下需要合并。为了使得生成的面片对于后续的光栅化流程是友好的，因此使用到相机的距离作为判定准则。在透视投影下，image plane（垂直于相机视角的2D平面）的尺寸s（用水平方向的宽度表征）与相机距离z之间的关系可以用如下公式来表征：

其中，对应的是水平方向的FOV。根据这个关系，我们可以按照如下的算法来确定每个key的细分层级k：

float distanceToLod(float z)
{
  // smaller the tmp, higher the subdivision level
  // targetPixelSize refers to the pixels number that the final subdivided triangle covered
  float tmp = s(z) * targetPixelSize / screenResolution;
  return -log2(clamp(tmp, 0.0, 1.0));
}

函数中的参数z指的是从相机到key所对应的subtriangle的距离，从代码看来，其实现逻辑跟这里给出的不太一样：

首先需要计算对应三角形到相机的距离d，之后将d乘上一个lod常量（比如CPU给定的LOD偏移之类）得到三角形的lod距离，之后通过一个log函数就可以算出当前三角形对应的细分层级了。

下面给出LoadKernel的完整实现：

buffer VertexBuffer { vec3 u_VertexBuffer []; };
buffer IndexBuffer { uint u_IndexBuffer []; };
buffer SubdBufferIn { uvec2 u_SubdBufferIn []; };
void main ()
{
  // get threadID ( each key is associated to a thread )
  int threadID = gl_GlobalInvocationID .x;
  // get coarse triangle associated to the key
  uint primID = u_SubdBufferIn [ threadID ].y;
  vec3 v_in [3] = vec3 [3](
  u_VertexBuffer [ u_IndexBuffer [ primID * 3 ]],
  u_VertexBuffer [ u_IndexBuffer [ primID * 3 + 1]] ,
  u_VertexBuffer [ u_IndexBuffer [ primID * 3 + 2]] ,
  );
  // compute distance - based LOD
  uint key = u_SubdBufferIn [ threadID ].x;
  vec3 v [3]; 
  // what does this process do? Get The subdivision triangle's coordinates
  subd (key , v_in , v);
  float z = distance ((v [1] + v [2]) / 2.0 , camPos );
  int targetLod = int( distanceToLod (z));
  // write to u_SubdBufferOut
  updateSubdBuffer (key , targetLod );
}

1.3.3 T-Junction Removal

跟其他的自适应网格细分算法一样，本文所给出的算法也是会出现T-junctions问题的，如下图红色三角形跟绿色三角形交界位置所示，T-junctions在顶点高度调整后会出现裂缝问题。不过如果能够保证相邻两个Child Triangle的细分级别相差不大于1，那么就有希望消除T-junctions问题，如图中蓝色三角形与绿色三角形的交界处所示（但是即使将红色三角形继续细分一级，其与绿色三角形就能达到相差1级的标准，但这种情况下还是会存在T-junctions啊？下面会说，当细分标准以斜边中心点到相机的距离为判断标准，那么红色三角形就还会继续拆分，从而保证可以消除T-Junctions）

为了实现T-junctions的消除，首先需要保证相邻两个Child Triangle的key level相差不大于1，这个是如何做到的呢，首先需要先检测出会出现T-Junctions的三角形，这个是通过对当前三角形的爷爷的兄弟的LOD与当前三角形的LOD做比对完成的，如下图中0101的爷爷就是01，其兄弟就是00，比对两者的LOD就可以知道是否存在T-junctions。

另外还需要以三角形的斜边（hypotenuse）中心点到相机的距离作为三角形是否细分的判定依据，在这两个标准的处理下，下图红色三角形就会被进行两次细分（或者对绿色三角形进行合并），从而避免T-junctions。

buffer VertexBuffer { vec3 u_VertexBuffer []; };
buffer IndexBuffer { uint u_IndexBuffer []; };
in vec2 i_InstancedVertex ;
in uvec2 i_PerInstanceKey ;
void main () 
{
  // get coarse triangle associated to the key
  uint primID = i_PerInstanceKey .y;
  vec3 v_in [3] = vec3 [3](
  u_VertexBuffer [ u_IndexBuffer [ primID * 3 ]],
  u_VertexBuffer [ u_IndexBuffer [ primID * 3 + 1]] ,
  u_VertexBuffer [ u_IndexBuffer [ primID * 3 + 2]] ,
  );
  // compute vertex location
  uint key = i_PerInstanceKey .x;
  vec3 v [3]; 
  subd (key , v_in , v);
  //berp, Get The World Coorindates of Vertex
  vec3 finalVertex = berp (v, i_InstancedVertex );
  // displace , deform , project , etc .
}

这里给出的是RenderKernel Stage中的VS代码。

1.3.4 Results

为了验证实现效果，原文作者给出了一个使用displacement-mapping算法实现的地形渲染demo，链接中给出了对应的git地址，感兴趣的同学可以自行下载，下面是渲染的表现：

下面给出的是在Intel i7-8700k CPU, 3.70 GHz,NVidia GTX1080 GPU，8GiB显存配置下使用1080p分辨率绘制时的时间消耗，相机上使用的是固定相机朝向，但是允许缩放的配置，这里的时间消耗数据的每一行分别对应三个处理步骤的时间消耗，另外带有stdev后缀的指的是随着相机不断缩进（放大）过程中5000帧的一个标准差（standard deviation）。

从数据上来看，性能消耗很低，而且十分稳定（不会出现尖刺之类的问题），因为地形渲染消耗对于shading会有比较大的依赖，因此为了重点验证细分算法，这里使用的是constant color来进行着色处理，因此上面的消耗基本上就可以理解为地形渲染中顶点处理与subdivision算法的时间消耗。

1.4 Discussion

本文给出了一个通过Compute Shader完成的，细分层级上限不受限制的地形网格细分方案，在后续的工作中，可能会尝试一些更为复杂的细分scheme，比如Catmull-Clark Scheme，此外这里还给出了一些本文实现中的一些细节考虑：

需要为包含subdivision keys的buffer分配多大的内存空间？
这个取决于屏幕空间的三角面片密度，一个三角形作为一个节点，且最高的划分层级为max_level的话，那么最少需要分配3 * max_level + 1个节点，这个数值对应的是之前介绍过的adaptive subdivision中相邻两个Child Triangle层级相差不超过1的情况，最大需要分配个节点，这个数值对应的是完全细分的情况（uniform split）。

本文算法是否会受到浮点数精度的限制？
implicit subdivision scheme是使用整数表示的，本身不存在问题，但是前面通过多次累乘输出变换矩阵上会存在一定的风险，这边给出的结论是使用31级或者低于31级的情况不会有问题，但是如果最大细分级别高于31的话，是可能会因为浮点精度不足导致表现问题的。这里给出的一个解决方案是使用双精度浮点数（OpenGL4.0+之后的硬件支持），虽然不能完全消除这个问题，但是在出问题之前的最大细分级别应该是足以满足所有需要了。

是否能够将这个算法用在tessellation shader中来消除tessellation shader 的限制呢？
原文作者在github的demo-isubd-terrain demo中给出了实现源码，开发者可以根据软硬件条件决定使用哪种方案更好。

要实现对地形网格的细分有两种方案，分别对应本文的subdivision scheme以及传统的LOD（如CDLOD）方案，这两种方案优劣对比是怎样的？
原文作者认为表现取决于具体的平台，之前的github源码给出了instanced triangle grid进行tessellation的工具，因此可以很好的比对两种方案，下图给出了在NVidia GTX1080上不同sublevel下的per-instance性能消耗，可以看到，随着sublevel的增加，每个instance绘制的时间消耗是逐渐降低的（不论是LOD Stage还是Render Stage，不过这种对比貌似也不能说明什么问题，因为虽然per-instance消耗降低了，但是instance数目增加了）

本文方法是否可以用于除地形网格细化之外的其他区域，比如模型细分？
implicit subdivision scheme提供了跟tessellation shader相同的功能，因此凡是tessellation shader使用的地方都可以使用。原文作者使用此方法分别应用[Vlachos et al. 01] 中的PN-triangles算法以及[Boubekeur and Alexa 08] 的Phong Tessellation算法完成了对模型的面片优化，结果见下图：

参考文献

[1] Adaptive GPU Tessellation with Compute Shaders
[2] GPU Tessellation with Compute Shaders