Cesiumで標高計算処理を他言語に移植するために関連コードを調べた時のメモ書きその２です。
その１はこちら：
本記事ではheightmap形式の地形データを使ったときの標高計算について調べました。
heightmapデータを使う場合HeightmapTerrainData#interpolateHeightで標高の計算が行われます。
コードを読むとencodingがLERCの場合とそうではない場合で大きく処理が分岐していることがわかります。
encodingがLERCの場合:
encodingがLERCの場合に最終的に呼ばれるのはinterpolateMeshHeight関数です。
ただしinterpolateMeshHeight関数を使うには、まずmeshデータ（TerrainMesh型）を作成する必要があります。
meshデータ未作成の状態でinterpolateMeshHeightが呼ばれると
- interpolateMeshHeightはundefinedの結果を返す
- 呼び出し側コード(sampleTerrain.js)は結果がundefinedの場合、TerrainData#createMeshを呼び出し、meshデータ作成後に再びinterpolateHeightを呼び出す
- interpolateMeshHeightはmeshデータをもとにinterpolateMeshHeightで標高を計算する
といった流れで標高を計算しているようです。
encodingがLERCではない場合：
それ以外の場合に最終的に呼ばれるのはinterpolateHeight関数です（※）
※HeightmapTerrainDataのprototypeメソッドのinterpolateHeightと関数名が同じなので紛らわしいのですが、こちらは上記とは別に定義されているファイル内スコープのプライベート関数です。
紛らわしいので、本記事では以後こちらの関数をinterpolateHeight（mesh不要版）と呼ぶことにします。
どちらが呼ばれるのか？
CesiumTerrainProviderのcreateHeightmapTerrainDataのコードを読む限り、CesiumTerrainProviderから生成されたHeightmapTerrainDataのencodingがLERCになることは無さそうです。
以後の調査ではinterpolateHeight（mesh不要版）に焦点を当てることにします。
interpolateHeight（mesh不要版）
引数
まずは引数を見ていきます。
sourceHeights～heightまではheightmapのデータやデータの構造に関する情報でlongitude, latitudeが高さを求めたい地点の座標になります。
function interpolateHeight( sourceHeights, // 高さデータのバッファ elementsPerHeight, // データあたりの要素数 elementMultiplier, // 高さデータの乗数 stride, // 隣接する高さのデータへのステップ数 isBigEndian, // ビッグエンディアンか否か sourceRectangle, // TerrainDataのタイル範囲（地理座標） width, // TerrainDataの２次元配列の幅 height, // TerrainDataの２次元配列の高さ longitude, // 高さを求めたい地点の経度 latitude // 高さを求めたい地点の緯度 )
「simple array of 16-bit int, little-endian」という記述があるので
- elementsPerHeight ... 2
- elementMultiplier ... 256
- stride ... 1
- isBigEndian ... false
また「The total size of the post data is bytes.65 * 65 * 2 = 8450」という記述があるので
- width, height ... 65, 65
になるのかなと思います。
CesiumTerrainProviderのデフォルトのheightmapStructureの定義もそのようになっていました。
cesium/packages/engine/Source/Core/CesiumTerrainProvider.js at 1.115 · CesiumGS/cesium (github.com)
if (data.format === "heightmap-1.0") {
isHeightmap = true;
if (!defined(terrainProviderBuilder.heightmapStructure)) {
terrainProviderBuilder.heightmapStructure = {
heightScale: 1.0 / 5.0,
heightOffset: -1000.0,
elementsPerHeight: 1,
stride: 1,
elementMultiplier: 256.0,
isBigEndian: false,
lowestEncodedHeight: 0,
highestEncodedHeight: 256 * 256 - 1,
};
}
hasWaterMask = true;
terrainProviderBuilder.requestWaterMask = true;
cesium/packages/engine/Source/Core/CesiumTerrainProvider.js at 1.115 · CesiumGS/cesium (github.com)
function TerrainProviderBuilder(options) {
...
this.heightmapWidth = 65;
またsourceRectangleにはterrainProviderのtilingSchemeのtileXYToRectangleで計算した矩形範囲が渡ってきます。
該当箇所：
https://github.com/CesiumGS/cesium/blob/1.115/packages/engine/Source/Core/sampleTerrain.js#L229
処理の流れ
相対位置の計算
longitude とlatitude のsourceRectangle 内での相対位置を計算する。
（fromWest, fromSouth）
1の相対位置から最も近い整数の格子座標（南西、南東、北西、北東の四角形の格子座標）を計算する。
（westInteger, eastInteger, southInteger, northInteger）
また、(fromWest, fromSource)の四角形内での相対座標も計算する。
(dx, dy)
getHeight関数を使用して、2の四角形の南西、南東、北西、北東の標高を取得する。
(southwestHeight, southeastHeight, northwestHeight, northeastHeight)
※実装上の留意点としては、heightmapは左上（北西）から右下（南東）の順番に高さデータが格納されているのでgetHeightを呼ぶ前にsourceIntegerとnorthIntegerの上下をひっくり返しています。
triangleInterpolateHeight関数を使用して、四角形内の高さを補間する。
triangleInterpolateHeight関数処理
cesium/packages/engine/Source/Core/HeightmapTerrainData.js at 1.115 · CesiumGS/cesium (github.com)
function triangleInterpolateHeight( dX, dY, southwestHeight, southeastHeight, northwestHeight, northeastHeight )
(dx, dy)の位置が四角形の左上側の三角形に含まれるか右下側の三角形に含まれるかを判別する。
（dY < dXなら左上、それ以外の場合は右下）
線形補間で(dx, dy)の位置の高さを計算する。
左上側の三角形の場合、southwestHeight, northwestHeight, northeastHeightの高さを使う。
右下側の三角形の場合、southwestHeight, southeastHeight, northeastHeightの高さを使う。
