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

时间:2021-10-13 00:10:48

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

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

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

layer.json

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

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

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

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

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

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

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

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

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

如上就是对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的格网效果,如下图,相比高度图,这个看上去很简单,比较稀疏,所以数据量药小很多。

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

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

Cesium原理篇:3最长的一帧之地形(3: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。如下这个图很好的说明了水平裁剪的作用。

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

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

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

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

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