最近因為工作需要，開始接觸了一些三維數據格式。項目是基於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層的這個范圍從而方便后續的計算（例如，頂點信息壓縮時需要使用的瓦片最大最小經緯度）。具體存儲格式為：
3. json文件分析
首先是layer.json,文件內容如下，大部分內容都還比較直觀，其中available即瓦片的范圍和第0層的meta信息比較類似，這里也就不再展開了，這個validbounds即dem的范圍。minzoom和maxzoom為縮放層級，這個怎么改都不會影像文件的解析和顯示。meta.json的文件內容就更簡單了，這邊也貼個截圖， 其中Bounds和latlonBounds都是DEM的經緯度范圍。
總結：
最后簡單做個總結，其實Cesium的quantized-mesh數據格式也不復雜，不過中間為了數據壓縮加入了一些數據編碼算法。另外強制對齊這個也是個坑，之前沒有怎么解除過二進制文件，官方在這個地方說明也是一筆帶過，第一次做完了切片后數據一直沒法被解析，反復檢查了幾遍才發現問題，需要引起重視。
參考資料：