Cesium地形瓦片数据格式——quantized-mesh浅析

2024-01-09 09:19:58

最近因为工作需要，开始接触了一些三维数据格式。项目是基于Cesium.js进行的开发，三维场景免不了使用到地形数据，Cesium官方以及CesiumLab都提供了将DEM数据转换成quantized-mesh地形瓦片的生成工具。但是由这些工具生成的地形瓦片不太好控制（比如说我想更改其中部分三角网的高程），因此参考了官方文档写了一个三维瓦片生成的小工具，其中碰到了一点小坑，这篇文章在官方文档的基础上加入了一些个人理解，留作记录。

1. 数据组织形式

首先需要对数据组织形式有个清晰的认知，这边我先提前处理好了一个DEM，将其处理成地形瓦片的结果如下图所示，由于我在做数据处理的时候设置了最大生成层数为8，因此可以看到最上层是0-8九个文件夹以及layer.json，meta.json，这两个文件比较简单放在文章最后说，点开任意一个文件夹，下层依然是文件夹，再进入某个文件夹终于看到terrain数据，这里文件的组织形式其实就是官方文档的简单多分辨率四叉树网格金字塔（simple multi-resolution quadtree pyramid of meshes）。金子塔遵从TMS（Tile Map Service）切片规则，简单来说就是将全球的经纬度范围（-180~180，-90~90）先划分成两个格网(-180 deg, -90 deg) - (0 deg, 90 deg) ，(0 deg, -90 deg) - (180 deg, 90 deg)，即图中名称为0的文件夹下的文件，然后每个格网内四叉树递归划分......回到terrain数据格式上，理论情况下层级为2的文件夹下应该是0~7，8个文件夹（2^3)，每个文件夹下应该有0~3四个terrain文件，图中则是DEM所在范围的格网。

2. 瓦片格式分析

每个瓦片其实就是dem在这个瓦片范围内的点以及瓦片的顶点构建的三角网。

2.1 瓦片头信息

struct QuantizedMeshHeader
{
    // 瓦片中心点的地心坐标
    double CenterX;
    double CenterY;
    double CenterZ;

    // 这块瓦片的最低和最高的高程
    float MinimumHeight;
    float MaximumHeight;

    // 瓦片的球体包围盒信息，其中xyz与瓦片中心点的地心坐标一致，radius为包围球半径
    double BoundingSphereCenterX;
    double BoundingSphereCenterY;
    double BoundingSphereCenterZ;
    double BoundingSphereRadius;

    // 地平线遮挡点
    // 更多信息：http://cesiumjs.org/2013/04/25/Horizon-culling/ for more information.
    double HorizonOcclusionPointX;
    double HorizonOcclusionPointY;
    double HorizonOcclusionPointZ;
};

在生成头文件的包围盒信息时，我使用的是Ritter的包围球算法，半径结果与官方比较接近。Ritter的包围球算法比较简单，我在文末的参考资料中放了链接，有需要的话可以去学习学习。地平线遮挡点简单来说就是减少瓦片加载和渲染量的标记，每个瓦片都会有一个独特的点，当该点低于视锥的地平线时，该点所代表的瓦片将不会加载和渲染。它的计算方法我在官方文档中没有找到明确的定义，因此我稍微理解了一下视点的使用原理，使用了地心坐标系下坐标除以椭球在对应的xyz方向的半径，计算结果比较接近官方裁切工具裁切出来的实际值，就作近似代替了。

double HorizonOcclusionPoint[3] = { centerX / 6378137, centerY / 6378137, centerZ / 6356752.3142451793 };

2.2 瓦片顶点信息

struct VertexData
{
    unsigned int vertexCount;                // 顶点个数
    unsigned short u[vertexCount];        // 顶点横坐标
    unsigned short v[vertexCount];        // 顶点纵坐标
    unsigned short height[vertexCount]; // 顶点高程值
};

　　紧跟着头信息的是顶点信息，首先以一个无符号整型存储顶点的个数，然后是三个短整型数组存储顶点坐标，为了方便存储和数据传输，这里使用了zig-zag编码方法对数据进行压缩，根据官方文档给出了顶点坐标的解码方法：

int zigzagDecode(unsigned short value) {
    return (value >> 1) ^ (-(value & 1));
}

int u = 0;
int v = 0;
int h = 0;
int divideSum = 32767;
for (unsigned int i = 0; i < vertexCount[0]; i++ ) {
    u += zigzagDecode(uvh[i]);
    v += zigzagDecode(uvh[vertexCount[0] + i]);
    h += zigzagDecode(uvh[vertexCount[0] * 2 + i]);
    double lon = minx + double(u) / divideSum * (maxx - minx);
    double lat = miny + double(v) / divideSum * (maxy - miny);
    double hei = height[0] + double(h) / divideSum * (height[1] - height[0]);
}

　　这个解码是个累计的过程，也即后一个点是前一个点的相对位置（差值），得到的u,v,h并不是最终的经纬度结果，解码完成后还需进一步计算坐标，计算方法是目标= 最小 + (u/32767) * (最大-最小)。这个地方会有点难以理解，可以先从编码的角度看看，编码的代码如下。其中vectorx,vectory,vectorz分别是顶点经度、纬度、高程的数组，在一个瓦片中我们将每个方向均等分成2^15 = 32767段，这样我们在表示任意一个顶点的任意一个方向坐标的时候就可以使用 targetInt = 32767 * (tartgetDouble - minx) / (maxx - minx)公式将其转化为整形，然后Cesium针对每个顶点，使用转化后的顶点坐标与前一个顶点坐标各个方向的差值来尽量压缩这个整型使其绝对值尽可能小，以便使用zigzag编码进行压缩（差值的使用原因可以看看zigzag编码的原理，文末有参考链接）。理解了编码方式以后再回头去看解码方式好懂很多。

double vectorX[vertexCount];
double vectorY[vertexCount];
double vectorZ[vertexCount];

unsigned short zizagEncode(short value) {
	return (value >> 15) ^ (value << 1);
}

short* uList = new short[vertexCount];
short* vList = new short[vertexCount];
short* hList = new short[vertexCount];

for (unsigned int i = 0; i < vertexCount; i++) {
	uList[i] = 32767 * (vectorX[i] - minx) / (maxx - minx);
	vList[i] = 32767 * (vectorY[i] - miny) / (maxy - miny);
	hList[i] = 32767 * (vectorZ[i] - minHeight) / (maxHeight - minHeight);
}

for (unsigned int i = vertexCount -1; i > 0; i--) {
	uList[i] = uList[i] - uList[i - 1];
	vList[i] = vList[i] - vList[i - 1];
	hList[i] = hList[i] - hList[i - 1];
}

unsigned short* uvh = new unsigned short[vertexCount * 3];

for (unsigned int i = 0; i < vertexCount; i++) {
	uvh[i] = zizagEncode(uList[i]);
	uvh[vertexCount + i] = zizagEncode(vList[i]);
	uvh[vertexCount * 2 + i] = zizagEncode(hList[i]);
}

　　2.3 三角面顶点索引

　　顶点之后是用来强制字节对齐的占格符，用来确保IndexData16为2字节对齐和IndexData32为4字节对齐。在二进制文件中，经常会出现强制对齐的现象，本人是这样理解的：前面的变量类型有double,float,int,short，其中double,float,int一定是2或者4的整数倍，short的数组，一定是2的整数倍，但是不一定是4的整数倍，这个对齐的意思就是当定点数大于65536时，为了使数据对齐，就需要填充空格使得前面的字节总数一定是对应数据类型的整数倍。这个地方在解析的时候需要注意，如果漏掉了会出现解析不正确Cesium无法正常加载的情况。

　　然后便是三角面的顶点索引数据，也就是顶点坐标在顶点数组中的位置（下标）。数据三个一组，指定顶点如何链接在一起成三角形。如果tile具有超过65536个顶点，则tile使用IndexData32结构对索引进行编码。否则，它使用IndexData16结构。

struct IndexData16
{
    unsigned int triangleCount;                // 三角形个数
    unsigned short indices[triangleCount * 3]; // 三角形顶点索引
}
 
struct IndexData32
{
    unsigned int triangleCount;
    unsigned int indices[triangleCount * 3];
}

　　其中索引数据使用了webgl的高水位编码法(high water mark)进行了压缩，解码方式如下：

int highest = 0;
for (var i = 0; i < indices.length; ++i) {
    int code = indices[i];
    indices[i] = highest - code;
    if (code == 0) {
        ++highest;
    }
}

　　2.4 瓦片边缘顶点索引

　　三角面索引后是边缘顶点索引信息，这些信息记录下来主要是为了瓦片避免边缘贴合度不够产生的裂缝问题。边缘顶点由于相邻两个点的差值很大不适合做zigzag压缩编码因此直接存储了原始顶点索引id，根据顶点是否大于65536，定义的两种数据结构如下：

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];
}
 
struct EdgeIndices32
{
    unsigned int westVertexCount;
    unsigned int westIndices[westVertexCount];
 
    unsigned int southVertexCount;
    unsigned int southIndices[southVertexCount];
 
    unsigned int eastVertexCount;
    unsigned int eastIndices[eastVertexCount];
 
    unsigned int northVertexCount;
    unsigned int northIndices[northVertexCount];
}

　2.5 扩展信息

　　紧跟着边缘顶点的是瓦片的扩展信息，虽然说是扩展信息但是如果想要自行做地形切片，这部分基本上必须使用，特别是Meta元信息，因此非常有必要掌握。每个扩展信息按照头信息-主体信息的数据格式进行存储，其中扩展头信息的数据结构如下：

struct ExtensionHeader
{
    unsigned char extensionId;
    unsigned int extensionLength;
}

　　extensionId是当前扩展信息的标识字段，用来告知该扩展信息属于哪一类，目前Cesium定义的扩展信息有三类：①id=1，地形光照。②id=2，水面掩码。③id=4，元信息。其中地形光照和水面掩码本人项目中没有使用到因此也没有单独研究，这里仅说说meta信息。官方文档给出的数据结构为：

struct Metadata
{
    unsigned int jsonLength;
    char json[jsonLength];
}

　　通过名字也能明白meta信息存储的是json转化而来的字符串，刚刚说到如果使用自己做的地形切片，其实不需要对每一层进行切分，因此在小比例尺下，三维场景在远距离下基本上是看不到高低起伏的，因此TMS中一些低层级的瓦片（例如1-7层,对应第一张图中编号1-7的文件夹）其实是不需要的。注意这里说的是1-7，但是第0层是一定需要的，这跟Cesium的解析原理有关，第0层会在meta信息中，存储所有不同层级下的瓦片的id范围信息（即DEM所在TMS切片下的瓦片id范围），Cesium在对瓦片进行解析的时候首先会读取第0层的这个范围从而方便后续的计算（例如，顶点信息压缩时需要使用的瓦片最大最小经纬度）。具体存储格式为：

{ "available": [ [ { "endX": 3, "endY": 1, "startX": 3, "startY": 1 } ], [ { "endX": 6, "endY": 2, "startX": 6, "startY": 2 } ], [ { "endX": 12, "endY": 5, "startX": 12, "startY": 5 } ], [ { "endX": 25, "endY": 10, "startX": 25, "startY": 10 } ] ... ] }

3. json文件分析

　　首先是layer.json,文件内容如下，大部分内容都还比较直观，其中available即瓦片的范围和第0层的meta信息比较类似，这里也就不再展开了，这个validbounds即dem的范围。minzoom和maxzoom为缩放层级，这个怎么改都不会影像文件的解析和显示。meta.json的文件内容就更简单了，这边也贴个截图，其中Bounds和latlonBounds都是DEM的经纬度范围。

　总结：

　　最后简单做个总结，其实Cesium的quantized-mesh数据格式也不复杂，不过中间为了数据压缩加入了一些数据编码算法。另外强制对齐这个也是个坑，之前没有怎么解除过二进制文件，官方在这个地方说明也是一笔带过，第一次做完了切片后数据一直没法被解析，反复检查了几遍才发现问题，需要引起重视。

参考资料：

码农公寓

相关文章