Cesium原理篇：3最长的一帧之地形(3：STK)

Posted 2020-08-09 法克鸡丝

       有了之前高度图的基础，再介绍STK的地形相对轻松一些。STK的地形是TIN三角网的，基于特征值，坦白说，相比STK而言，高度图属于淘汰技术，但高度图对数据的要求相对简单，而且支持实时构建网格，STK具有诸多好处，但确实有一个不足，计算量比较大，所以必须预先生成。当然，Cesium也提供了一个Online的免费服务，不过因为是国外服务器，所以性能和不稳定因素都不小。好的东西自然得来不易，所以不同的层次，根据具体的情况选择不同的方案，技术并不是唯一决定因素，甚至不是主要因素。

       CesiumTerrainProvider提供了高度图和STK两种地形服务。但目前只提供后者的在线服务，前者其实只是在维护了(HeightmapTessellator在计算Horizon Cull中有一处笔误，错把relativeToCenter写成relativetoCenter了，一直到目前最新的1.25版本还没有修改，说明已经不重视了)。

       下面，我们通过CesiumTerrainProvider类，来详细介绍一下STK地形的细节。

layer.json

       首先，CesiumTerrainProvider在初始化的时候会请求layer.json，这相当于一个配置文件，用来获取地形数据的基本属性信息，截图如下：

 

 

       如上是STK对应的layer.json中的属性，包括版本，format(高度图 or STK)，extension扩展(是否包括法线和水面)，层级数和范围。最主要的是available，里面记录有地形的切片数据，这样，如果该Tile的地形均为0，减少不必要的terrain请求。比如Level 11，在[startX,startY]到[endX,endY]范围下存在地形数据。当然，提供了CesiumTerrainProvider.prototype.getTileDataAvailable来进行这个判断。

STK Structure

       上一篇，我们讲到createHeightmapTerrainData方法，参数对应的是一张高度图，里面就是记录高程值的点串。而createQuantizedMeshTerrainData则不同，它对应的不是原始的，简单抽稀的点串，而是已经预处理的，缓存在服务器中的TIN地形格网。所以，首先我们要清楚该terrain文件的数据结构。

       首先，在Cesium中是以ArrayBuffer形式传输的，也就是二进制形式。它的Header结构如下：

struct QuantizedMeshHeader
{
    // 切片的中心点，对应球面坐标系，米单位
    double CenterX;
    double CenterY;
    double CenterZ;
    
    // 该切片地形高度的最大值和最小值，米单位
    float MinimumHeight;
    float MaximumHeight;

    // BoundingSphere的中心点和半径，米单位
    // 具体计算方式和高度图一文中的算法一致，
    double BoundingSphereCenterX;
    double BoundingSphereCenterY;
    double BoundingSphereCenterZ;
    double BoundingSphereRadius;

    // 水平裁剪的XYZ值，后面会介绍其原理
    double HorizonOcclusionPointX;
    double HorizonOcclusionPointY;
    double HorizonOcclusionPointZ;
};

       这个头介绍了该Tile的大概的位置和范围，之后就是该Tile对应的TIN三角网的顶点数据,该数据里面用到了zigzag编码方式，后面会具体介绍：

struct VertexData
{
    // 顶点数据
    // 分别对应uv和height
    // 里面是相对值而不是真实值
    // 详细内容会在下面介绍
    unsigned int vertexCount;
    unsigned short u[vertexCount];
    unsigned short v[vertexCount];
    unsigned short height[vertexCount];
};

       接着就是顶点索引了，采用了high water mark编码方式，同样后面一并介绍：

// 顶点索引，分为short和int两种情况
// 这里是short的结构，int的同理
struct IndexData16
{
    // 顶点数
    unsigned int triangleCount;
    // 顶点索引数据
    unsigned short indices[triangleCount * 3];
}

       下面则是裙边，上一篇也提到过，主要就是勾勒一下四周，保证无缝合成，后面会给出一个STK的wireframe方式，让大家看清楚高度图和STK的差别，也体验一下STK数据的简单高效：

// 裙边对应的顶点索引
struct EdgeIndices16
{
    //  对应东南西北四个边
    //  顶点数和顶点索引
    unsigned int westVertexCount;
    unsigned short westIndices[westVertexCount];

    unsigned int southVertexCount;
    unsigned short southIndices[southVertexCount];

    unsigned int eastVertexCount;
    unsigned short eastIndices[eastVertexCount];

    unsigned int northVertexCount;
    unsigned short northIndices[northVertexCount];
}

      最后就是一些扩展属性，主要是水面和法线（光照效果），这里就不在赘述。有了这个数据结构，我们就可以按这个结构来解析对应的ArrayBuffer了，这个过程就不在此介绍了，对ArrayBuffer不了解的，可以参考《ArrayBuffer简析》。数据全部解析完后，我们知道了该Tile的范围以及高度的最大最小值，通过OrientedBoundingBox.fromRectangle方法，可以构造出OrientedBoundingBox。思路和高度图中的介绍一样。最终，如上的参数构造了QuantizedMeshTerrainData对象。

Decode

       Cesium在顶点数据中采用了zigzag的编码方式，在顶点索引中采用了high water mark方式，两个编码的解码方式都很简单，算法如下：

// zigZagDecode
decoded = (encodedValue >> 1) ^ (-(encodedValue & 1))

// High water mark decoding
var highest = 0;
for (var i = 0; i < indices.length; ++i) {
    var code = indices[i];
    indices[i] = highest - code;
    if (code === 0) {
        ++highest;
    }
}

       解码算法比较简单，但肯定会有一个疑问，为什么需要编码，还需要解码，这不是多此一举吗？当然，这样做是有原因的，或者说，通过这种方式可以更高效的压缩数据（整数）。当然要理解这些，还得先从varint这个类型说起。

Varint编码

       Varint 是一种紧凑的表示数字的方法。它用一个或多个字节来表示一个数字，值越小的数字使用越少的字节数。这能减少用来表示数字的字节数。不过相应的，对于大数就要使用更多的字节去存储。在统计学上，一般消息中的数字以小数为主，所以用它可以省空间。如下图所示：

 

 

       在一些通讯协议中，比如Google的ProtoBuffer，采用varint的方式来减少传输。如果整数大小在256以内，则只需要一个字节就可以存储该整数，则可以省去3个字节。因此，为了更好的发挥varint编码的优势，STK中则需要有一种策略，保证数值越小越好。

       首先就是ZigZag编码，该编码会将有符号整型映射为无符号整型，以便绝对值较小的负数仍然可以有较小的varint编码值。因为对于负数，最高位是1，那么就相当于一个很大的整数，如果用varint，那么就很浪费空间了。ZigZag编码的原理就是按照绝对值大小来重新解析二进制。如下是一些具体数值经过zigzag编码后的值，转为无符号整型。

 

 

      另外，对于顶点索引这样的整数，都为正数，而且遵循一定的顺序，是否也能让他们变小一些，这样对处理后的数进行varint编码，也会提高顶点索引的压缩比。于是就有了high water mark算法，经过优化后的顶点索引，要么是你之前见过的，要么是你见过的最大值+1。这样就不需要对真实值（较大）编码，而是对相对值（较小）进行编码就可以。参考如下：

 

 

        如上就是对zigzag和high water mark算法思路和作用的一个意会。详细的编码解码在网上也能找到，这里就不多介绍了。另外，可以看到他们解码的代价是很低的，因此在复杂度和压缩比之间，还是会有很大的受益，特别适合Web下的传输，当然一切的前提是你得使用varint这种编码方式，或者你在传输的过程中采用了体现这种价值的压缩方式。

       另外，Cesium的顶点数据保存的是delta（当前值-前一个值），这个过程可能会产生负值，这是采用zigzag的主因。解码代码对应的也是相对值，注意是+=而不是=：

for (i = 0; i < vertexCount; ++i) {
            u += zigZagDecode(uBuffer[i]);
            v += zigZagDecode(vBuffer[i]);
            height += zigZagDecode(heightBuffer[i]);

            uBuffer[i] = u;
            vBuffer[i] = v;
 }

createMesh

       数据解析完毕后，同高度图一样，把上面的参数，通过Workers技术，在线程中开始构建格网。这个过程是在QuantizedMeshTerrainData.createMesh中完成了，格网的构建是在createVerticesFromQuantizedTerrainMesh方法中实现。这个过程和高度图是完全一样的，所以在此略去。

       我们先看看STK的格网效果，如下图，相比高度图，这个看上去很简单，比较稀疏，所以数据量药小很多。

 

              我们换一个视角再看看，就会发现，虽然稀疏，但把好身材（前凸后翘）都体现出来了，并且裙边也处理的很严谨：

 

 

       举一个不恰当的例子，有一个性感的模特，一个裁缝手艺不到家，就做了一件紧身衣，来突出模特的身材；另一位手艺了得，根据模特的身材，量身制作了一件轻易的霓裳羽衣，穿在身上，模特的身材不仅表现的淋漓尽致，配上衣服更有一番若隐若现的韵味。这就是高度图和TIN在三角网处理上的差距。

       先对比了两者之间的差别后，我们正式进入构建网格的环节。因为STK的数据是预处理的，而且数据比较稀疏。如果很好的理解了之前的数据结构，这一块理解起来也不是问题，主要是把顶点数据解析成球面坐标下对应的真实值。简单来说，如上图所示，此时，格网中的每一个节点，目前是按照0~32767(short的最大值)的范围，保存的一个正整数（比例系数），现在就要根据该Tile的地理范围，实际的高度，通过这个比例系数，还原成真实的，具有地理意义的真实值。不知道大家是否理解这个意思，一言不合就上代码：

var maxShort = 32767;
// quantizedVertexCount为顶点数据的总数
for (var i = 0; i < quantizedVertexCount; ++i) {
    //  uv为该点对应该Tile下[0,1]范围内的位置
    var u = uBuffer[i] / maxShort;
    var v = vBuffer[i] / maxShort;
    // 高度值，米单位
    var height = CesiumMath.lerp(minimumHeight, maximumHeight, heightBuffer[i] / maxShort);

    // 通过插值算法，计算对应的经纬度值，弧度单位
    cartographicScratch.longitude = CesiumMath.lerp(west, east, u);
    cartographicScratch.latitude = CesiumMath.lerp(south, north, v);
    cartographicScratch.height = height;

    // 经纬度转换为以地球球心为原点的笛卡尔坐标系的值，米单位
    var position = ellipsoid.cartographicToCartesian(cartographicScratch);

    uvs[i] = new Cartesian2(u, v);
    heights[i] = height;
    positions[i] = position;
    
    // ……
}

TerrainEncoding

       格网对应的节点构建完成后，Cesium还做了TerrainEncoding编码，主要是看数据是否可以压缩，比如把两个float值压缩为一个float，这样来降低显存的占用。首先，格那句XYZ三个轴的最大最小值，构造出aaBox和TerrainEncoding对象：

var aaBox = new AxisAlignedBoundingBox(minimum, maximum, center);
var encoding = new TerrainEncoding(aaBox, hMin, maximumHeight, fromENU, hasVertexNormals);

       前者是一个包围盒，后者是后面需要用的压缩工具，如果包围盒三个维度中最大的差距在2^12 = 4096 范围内，则把两个float值压缩为一个float，如果差距太大，超过这个范围，认为这种压缩对精度的损失比较大，泽不进行编码，还是有原始数据。这里我有一个疑问，为什么范围设定在4096，一个float是32位，如果是二合一，理论上最大范围可以设定在2^16。我不确定我的推断是否正确：因为float的精度原因，尾数占23bit，精度只有7位，所以4096倍的放大取整就足够了。当然，一拆二的过程（把这一个float还原成两个float）是在shader中计算的，性能没有影响。通过这种方式，对小范围的Tile而言，可以减少50%顶点数据的显存占用，也是一种优化吧。下面给出二合一的压缩算法，不难理解就不多解释了：

AttributeCompression.compressTextureCoordinates = function(textureCoordinates) {
    var x = textureCoordinates.x === 1.0 ? 4095.0 : (textureCoordinates.x * 4096.0) | 0;
    var y = textureCoordinates.y === 1.0 ? 4095.0 : (textureCoordinates.y * 4096.0) | 0;
    return 4096.0 * x + y;
};

       至此，我们的Mesh就构建完成了，细细的和高度图createMesh的过程对比 ，真的可以用疏归同途来形容，方法不一，但最终都能构建出TerrainMesh这个对象。

Horizon Cull

       STK的TIN三角网方式就介绍到这了，上一篇因为篇幅问题，没有将水平裁剪，在这里补充一下，毕竟无论是高度图还是TIN，都用到了这个技术。下面介绍一下它的主要思路和关键实现点，详细内容可以访问Cesium。如下这个图很好的说明了水平裁剪的作用。

 

 

       在高度图中我们提到了Frustum Cull（视锥体裁剪），无论是BoundingSphere还是OrientedBoundingBox，都是用来判断图片中红点区域的，而蓝点虽然也在椎体内，却在背面，其实可以不用渲染，而Frustum Cull不能判断出这种情况。而Horizon Cull就是用来解决这种情况的。

       首先的问题是，怎样判断这个点是否在球的背面，这时一个挺有意思的话题，之前研究OpenWebGlobe的源码，他也有水平裁剪的判断，不过相比而言计算量就有点大，因为他计算的是该点和圆心的法线，以及相机和圆心的向量之间的角度，如果大于90度，则是背面。我觉得Cesium的数据功底真的很强，理论上打造的很扎实，不是随便凑凑，灵机一动就能想到的解决技术。如下图：

 

       一个复杂的向量计算，转化成了一个相似三角形的问题，即判断VA和VP之间的大小。具体的推导过程就不再这里额外展开了，Cesium的官网博文里面有很详细的推导，而且源码里面也有对应的判断。下面就是第二个问题，如何把每一个Tile面抽象成一个点，不然需要判断该TIle对应的所有点，这个计算量也是不实际的。同样，Cesium中也提供了推导过程和相关代码，请参考。额外说一下，在推导中，你需要注意你的球面参数是用的椭球和圆球，参数的不同，也会导致结果，有时候这种差异不是误差，而是错误，而Cesium中采用的是椭球，长轴、短轴和曲率相关的参数。