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

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

       此时,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

9

       如上是一个流程示意图,横线以上的是主线程,调用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

网格节点

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

2

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

10

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

11

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

12

 

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

13

 

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

14

 

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

image

       如上,通过两个遍历,我们得到了和网格节点一一对应的三个数组: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,也就是一个物体的外接圆。

image

       现在,我们理解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算法相似,也是需要遍历两遍,第一次遍历后,估算出一个初始的球,然后再次遍历,如果点在这个球内则什么也不做,如果点在球外,则调整中心点和半径,确保该点在球内。调整算法如下:

image

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

image

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

OrientedBoundingBox

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

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

craterlake_comp

       同样,我们要考虑两个问题,获取这个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坐标的过程,进而获取该点对应椭球的法线方向,点+法线 = 切面。如下图,红线是一个二维椭圆切线示意图,对应三维椭球就是一个切面:

image

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

image

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

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

image

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

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

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

image

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

image

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

image

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

posted @ 2016-09-16 13:38  fu*k  阅读(20314)  评论(11编辑  收藏  举报