Cesium原理篇：3最长的一帧之地形(2：高度图)

Posted 2020-08-09 法克鸡丝

       这一篇，接着上一篇，内容集中在高度图方式构建地球网格的细节方面。

       此时，Globe对每一个切片（GlobeSurfaceTile）创建对应的TileTerrain类，用来维护地形切片的相关逻辑；接着，在requestTileGeometry中，TileTerrain会请求对应该切片的地形数据。如果读者对这部分有疑问的话，可以阅读《Cesium原理篇：1最长的一帧之渲染调度》；最后，如果你是采用的高度图的地形服务，地形数据对应的是HeightmapTerrainData类，最终，该TerrainData形成了一个TerrainMesh网格。下面，我们就详细的介绍一下最后一步的相关内容。

高度图

       首先，怎么理解高度图？通常一个Tile都会对应一个256*256的影像切片，代表该Tile对应的XYZ范围下对应的影像内容，高度图也是一样的思路，只是此时，前者每一个像素代表当前位置对应的颜色，而后者代表当前位置对应的高度

       一般情况下，高度图是一个缩略图，比如在Cesium中，在没有真实地形数据下，高度图的宽高是16*16大小，每个点对应的值都是0，在有真实地形数据下，高度图是65*65的大小。可见，高度图是有抽稀的，而不是一一对应，一来没必要，不然就是点云了，二来构网的计算量很大，也是效果和效率的一个折中。

Workers线程

       Cesium的调度是基于状态的变化，看似简单，但个人认为非常精髓的。相比基于事件驱动的策略，基于状态可以更好的实时的处理大数据，逻辑上也简单清晰，当然这是题外话，我们继续回到地形本身。

       有了数据，TileTerrain的状态由RECEIVING变为RECEIVED，自然也就进入了下一环节transform：将原始的地形数据（HeightmapTerrainData）转换为格网（TerrainMesh）的过程。

这个过程涉及到不小的计算量，因此，Cesium采用Promise + Workers技术，把计算量放到线程中，这样保证界面操作的流畅。对Workers感兴趣的可以参考《Cesium原理篇：4Web Workers剖析》。

       HeightmapTerrainData.prototype.createMesh方法提供了构建格网的方法，内部正是采用Workers线程的方式，下面我们进入主题，详细介绍高度图构网的细节。

HeightmapTessellator

       如上是一个流程示意图，横线以上的是主线程，调用createMesh，创建线程，把buffer（高度值数组），宽高（width&height），tile的范围（rectangle）和中心点（center）等作为createVerticesFromHeightmap函数的参数，这样，每一个Tile都会创建一个Worker线程，并在线程中实现网格的构建。

Paremeters

       网格构建的算法则封装在HeightmapTessellator.computeVertices函数中，我们先详细了解里面的参数：

• Heightmap
  用于构建格网的高度图点串
• Width&height
  高度图的像素宽高
• skirtHeight
  俗称裙边，每一个Tile四周会围成一个栅栏，指定该栅栏的高度，保证和相邻的Tile拼接时没有间隙
• nativeRectangle
  该Tile的范围，如果是WGS坐标系，单位是度，如果是墨卡托，单位是米
• exaggeration
  地形高度的缩放系数，通常为1，现实真实的地形高度
• rectangle
  该Tile对应的地理范围，单位是弧度，rectangle和nativeRectangle至少要有一个，如果两个参数都有，则互相是匹配的
• isGeographic
  true则为WGS坐标，false为墨卡托
• relativetoCenter
  该Tile对应的中心点，单位是基于球心的笛卡尔坐标，单位为米
• ellipsoid
  椭球体类，提供一些计算和换算方法
• structure
  高度图数据结构，后续再说，感觉有点鸡肋

       以上就是必须传入的参数（也有一些不准确，如果没有设置这些参数，则会采用默认值），当然还有一些可选参数。然后就正式开始构网了。构网的过程主要分为四个部分：

• 构建网格
• 计算BoundingSphere
• 计算HorizonCulling
• Encode

网格节点

       构建网格的代码很长，但仔细读一下其实不难理解，结合下面一张图，先跟大家解释一下思路。

       这里先把这个网格想象成正方形的（就像把地球仪平铺成一幅地图），这就是一个Tile对应的格网。通过之前传入的参数，我们已知Tile的长宽（rectangle），行列的格子数（width&height），不难计算出每一个节点的位置（经纬度），说白了，就是两个for循环嘛，伪代码如下：

       当然这是一段伪代码，如果这样写确实也很短，但Cesium认为构建网格的计算量大，也很频繁，所以在此处进行了优化，是一个简化版的cartographicToCartesian函数，这一块说复杂也复杂，需要你对椭球体能有所了解，说简单也可以，因为即使不了解，你也可以直接套公式。大致的图示和公式如下：

       其中，B是纬度，L是经度，N是长半轴，这里是地球半径6378137米，而N（1-e^2）是椭圆的短半轴，这里取值6356752.3142451793米。另外，网格中节点数和高度图的宽高是一致的，这样每一个节点都对应高度图中的一个高度heightSample，这样，套用上面的公式，对应实现代码和注释如下：

 

      这样，对应格网中每一个节点，我们可以计算出positions和heights这两个数组，同理，再次把这个Tile网格想象成平面的，每一个Tile也对应一张影像切片，假设把这两张半透明的纸叠在一起，下面那张的是网格，上面那张是影像图，就是如下的这个效果：

 

      通常影像切片是256*256像素大小，我们把[0,256]的像素范围映射到[0,1]的比例中，这样，也能够计算出每一个节点对应[0,1]的比值，也就是通常说的uv（OpenGL里面的纹理坐标，渲染纹理时需要用到该参数）。Cesium中实现uv的代码如下：

 

      此外，在计算网格节点时，还计算了每个节点距离该Tile中心点relativetoCenter的距离，这个在下面计算boundingsphere时会需要。 这需要掌握图形学和矩阵方面的一些数学基础：

       如上，通过两个遍历，我们得到了和网格节点一一对应的三个数组：positions，heights，uvs以及距离中心点的最大值maximum和最小值minimum。

Cull裁剪

       如果只是单纯的网格构建，工作已经完成，但实际中还远远不够，最直接的一个问题是你不知道是否需要显示这个Tile。

       如果做过渲染优化的人有会有这样一个共鸣吧，渲染一个物体最快的方式就是不去渲染它。看上去很正确，但做起来其实是一个很严肃的问题。把这句话转译成程序员的语言就是，判断这个物体的范围是否在当前可视范围内。在结合我们正在讨论的Cesium地形，考虑到大规模，频繁的渲染环境下，在相机的视锥体下，如何快速，简单的判断当前地形格网是否可见，这是一个严肃的问题。而Cesium在这个问题上，做到了极致，让我深为叹服。

       首先，Cesium主要采用了两种裁剪方式：

• Frustum Cull
• Horizon Cull

       因为里面涉及到很多算法，坦白说，每一个单独的细节，展开讲都很有学问，所以下面主要是思路和个人的理解，我也尽量把条理说的清楚一些，让大家能够一个完整的认识。由于篇幅过长，本篇主要介绍锥体裁剪部分，水平面裁剪在TIN地形的时候在涉及。

Frustum Cull

       首先，当一个物体不在视锥体范围内，自然就不需要显示了。视锥体的大小是清楚的，所以，剩下的就是如何计算该物体的Bounds。同时，从世纪角度来看，这个判断的过程一定不能超过渲染该物体的时间，否则也是没有意义的。因此，构建这个Bounds的关键就在于快速和有效之间的平衡。即能够较快的构建出一个近似准确的Bounds，同时这个Bounds也能高效的较为准确的判断是否可见。

BoundingSphere

       Cesium最先提供的是BoundingSphere，如下图，就是一个模型的BoundingSphere，也就是一个物体的外接圆。

       现在，我们理解BoundingSphere的概念，那么我们有一堆点串，如何实现BoundingSphere.fromPoints这个函数呢？在阅读下面的内容前，希望大家也琢磨一下这个问题。

       BoundingSphere就是一个球，所以这个问题就是获取球的球心和半径。之前我遇到过类似的问题，有一堆点串，已知中心点的情况下，计算其半径。我的思路是遍历所有点，计算每个点和中心点距离，取最大值作为半径。

       现在，这个中心点是未知的，所以我们需要先遍历所有点，找到XYZ三个方向的Max和Min，即X（Min，Max），Y（Min，Max），Z（Min，Max），然后计算Min和Max的均值，作为中心点，即：P = （Min+Max） / 2。这样有了中心点，前面也给出了计算半径的思路，我们就实现了BoundingSphere.fromPoints。

这个算法不难理解，也是最简单最快速的方式，在Cesium中，称这个算法为Naïve Method，看到Naïve，不知道有几个人会会脱口而出“图样图森破”？但这个算法也有一个缺点，这个球通常都不是最优的，就像你穿了一件大一号的衣服，略微不太优雅。

       接着，Cesium对比了Jack Ritter算法。这个算法和Naïve算法相似，也是需要遍历两遍，第一次遍历后，估算出一个初始的球，然后再次遍历，如果点在这个球内则什么也不做，如果点在球外，则调整中心点和半径，确保该点在球内。调整算法如下：

       Naïve和Jack Ritter相比，第一次遍历过程基本一致，但第二次遍历时，Naïve只修改半径，而后者对中心点和半径都会调整。Jack Ritter自己测试，两者在计算量上相当，但后者要多5%的准确性。但Cesium自己测试发现，19%的情况下，效果会比前者差，而11%的情况下，效果会比前者好，说明第一次估算的球和添加点的顺序也会影响Jack Ritter算法的结果。

       如上是一个测试数据对比（参考），最后在Cesium里面会同时执行Naïve和Ritter算法，以半径最大的值作为最后的结果，这个思路是可取的，两次遍历的计算内容99%都一样，就像拼车一样，举手之劳，受益良多。有了BoundingSphere，如何判断是否在视锥体范围内呢，在《Cesium原理篇：2最长的一帧之网格划分》里面有详细解释，这里就不赘述。

OrientedBoundingBox

       就这样相安无事了不久，人们对性能的追求始终没有停止。之前的BoundingSphere发现还是不够精确，你也看到了，这个球里面有不少的空白区域，造成了过度的渲染，那有没有更精准的Bound呢，这个就是OrientedBoundingBox了。

       BoundingBox是指包围盒，再加上Oriented，顾名思义就是有朝向的包围盒了。上一个对比效果图，可以看到这个其实就是在本地坐标系下的一个包围盒，如下是一个OrientedBoundingBox和对比图，可见，确实范围要小很多。

       同样，我们要考虑两个问题，获取这个Bounds的成本，以及判断Bounds是否可见的成本。

       首先，对应一个地形Tile，总会有一个中心点，也就是参数relativetoCenter，该点对应球面的切面+法线，构成了这个local coordination（NEU：north east up）的XYZ轴，这样一个相对NEU坐标的正交geometry，相对于球心的笛卡尔坐标系就是一个斜geometry了，这就不像BoundingSphere那样具有更好的对称性，可以很直白的用参数化的方式来构造了。OrientedBoundingBox默认是一个2*2*2的正方体，center是该包围盒的中心点，而还有一个矩阵halfAxes用来记录包围盒按照中心点旋转和缩放信息。

       下面简单介绍一下如何获取一个地形Tile对应的OrientedBoundingBox，也就是Cesium.OrientedBoundingBox.fromRectangle函数。

       首先，我们知道该Tile对应的relativetoCenter，然后构造出EllipsoidTangentPlane对象，也就是该点对应的椭球切面，这个过程其实就是从笛卡尔坐标转为NEU坐标的过程，进而获取该点对应椭球的法线方向，点+法线 = 切面。如下图，红线是一个二维椭圆切线示意图，对应三维椭球就是一个切面：

       此时，我们只是获取了OrientedBoundingBox的XYZ的三个朝向，并在XYZ三个方向上无限延伸，但无法确定具体的范围。如何计算这个范围呢？

       先举个简单的例子，上面在一个黑暗的屋子里，你手拿这个足球正对着一面墙，在球心放一盏灯，。假设黑色的部分是透光的，而白色的不透光，这样墙上会有黑色切面的影子。我们把球不断的靠近这面墙，直到刚刚好贴在墙面上，这时每一个切面对应的影子是最小的。

       很显然，这个切面是一个弧面，而墙是一个平面，刚才的过程其实就是把XYZ的三维体投影到XY平面的过程。Cesium也是用同样的思路，通过EllipsoidTangentPlane.prototype.projectPointsToNearestOnPlane方法计算每一个Tile地形在XY上的范围（此时并不是计算所有点，而是类似九宫格，计算Tile对应左中右上中下这九个点来配准），而Z的范围可以简单的理解为当前地形数据对应的最高点和最低点，从而得到该Tile对应的minX, maxX, minY, maxY, minZ, maxZ（该值是以relativetoCenter为原点，米单位）。

       最后通过如上获取该包围盒准确的中心点以及相对偏移量和缩放比，结合NEU矩阵最终构造出halfAxes矩阵。

       这样，我们就找到了地形Tile对应的OrientedBoundingBox，坦白说，理解上要比BoundingSphere复杂很多，但在计算量上，因为只针对特征点来计算，其实性能要好，当然，个人没有测试，只是推测。

       如何判断包围盒和视锥体的位置关系，这个和BoundingSphere的算法非常相似。那这种方式效果是否有改进，还得看疗效，如下是Cesium自己提供的对比效果（参考）：

       可见，改善效果还是很不错的，而判断是否相交的性能上差别不大：

       实际上，两种Bounds方式在Cesium中都在使用，而且，计算格网positions和两个bounds（BoundingSphere，OrientedBoundingBox）中，有一些重复计算的部分，所以还是有一定的优化空间。但这主要是在编程技巧上的对比，从逻辑和算法上，Cesium已经非常专业了，我个人觉得在研究源码时，在方方面面，都受益匪浅。

       如上，我们计算了positions，heights，uvs以及bounds后，我们基本完成了HeightmapTerrainData.CreateMesh的过程，也是地形中最关键的环节，下一步，就是开始加载到显卡中，通过shader渲染了，我们在后续会介绍。同时，由于篇幅的问题，临时决定把水平裁剪和Encode的部分取消，后面找个合适的机会在介绍。