第6章: SSEによるLOD制御

/** タイルの Geometric Error（メートル単位の地表幅） */
export function tileGeometricError(z: number): number {
  // 赤道での1タイル幅 [m] = 2*PI * R / 2^z
  return (2 * Math.PI * WGS84.a) / (1 << z);
}

/** Screen Space Error を計算 */
export function computeSSE(
  geometricError: number,
  distance: number,
  screenHeight: number,
  fovRad: number
): number {
  if (distance <= 0) return Infinity;
  return (geometricError / distance) *
    (screenHeight / (2 * Math.tan(fovRad / 2)));
}

           (0,0,0)          z=0
          /   |       (0,0,1) (1,0,1) ...     z=1
    / | (0,0,2) (1,0,2) ...         z=2

export function selectVisibleTiles(
  cameraPosition: { x: number; y: number; z: number },
  fovRad: number,
  screenHeight: number,
  frustumPlanes: FrustumPlane[],
  sseThreshold = 400,
  maxZoom = 6,
  maxTiles = 128
): VisibleTile[] {
  const result: Array<VisibleTile & { dist: number }> = [];

  // ホライズンカリング用の事前計算（ループ外で1回だけ）
  const camDist = Math.sqrt(
    cameraPosition.x ** 2 +
    cameraPosition.y ** 2 +
    cameraPosition.z ** 2
  );

  function traverse(x, y, z): void {
    const center = tileCenterScene(x, y, z);
    const radius = tileBoundingSphereRadius(x, y, z);

    // 1. ホライズンカリング（内積1回で済む）
    if (!isAboveHorizon(center, radius, cameraPosition, camDist))
      return;

    // 2. Frustumカリング（6平面テスト）
    if (!sphereInFrustum(center, radius, frustumPlanes))
      return;

    // 3. SSE判定
    const distance = ...; // カメラまでの距離
    const sse = computeSSE(
      tileGeometricError(z), distance, screenHeight, fovRad
    );

    if (sse <= sseThreshold || z >= maxZoom) {
      result.push({ x, y, z, dist: distance }); // 採用
      return;
    }

    // 4. 4つの子タイルに分割して再帰
    traverse(x*2,   y*2,   z+1);
    traverse(x*2+1, y*2,   z+1);
    traverse(x*2,   y*2+1, z+1);
    traverse(x*2+1, y*2+1, z+1);
  }

  traverse(0, 0, 0);
  // maxTiles超過時は距離ソートして近いもののみ
  return result;
}

traverse(0,0,0)
  |
  +-- ホライズンカリング --> 地球の裏側? --> return
  |
  +-- Frustum外? --> return（カリング）
  |
  +-- SSE <= 閾値 or z >= maxZoom? --> result.push（採用）
  |
  +-- SSE > 閾値 --> 4分割して再帰
       +-- traverse(0,0,1)
       +-- traverse(1,0,1)
       +-- traverse(0,1,1)
       +-- traverse(1,1,1)

if (result.length > maxTiles) {
  result.sort((a, b) => a.dist - b.dist);
  result.length = maxTiles;
}

          カメラ *
         /       |
        /   見える|
       /         |
------*----------|---- ホライズン（水平線）
     接点    見えない
              *
             タイル

/** カメラ方向の楕円体半径を算出する。
 * 楕円体方程式: (x²+z²)/a² + y²/b² = 1 */
function ellipsoidRadiusInCameraDir(
  dirX: number, dirY: number, dirZ: number
): number {
  const a = WGS84.a;
  const b = WGS84.b;
  return 1 / Math.sqrt(
    (dirX * dirX + dirZ * dirZ) / (a * a)
    + (dirY * dirY) / (b * b)
  );
}

export function isAboveHorizon(
  tileCenter: { x: number; y: number; z: number },
  tileRadius: number,
  cameraPosition: { x: number; y: number; z: number },
  cameraDistance: number
): boolean {
  // カメラ方向単位ベクトル
  const invD = 1 / cameraDistance;
  const camDirX = cameraPosition.x * invD;
  const camDirY = cameraPosition.y * invD;
  const camDirZ = cameraPosition.z * invD;

  // カメラ方向の楕円体半径（球体近似ではなく動的計算）
  const R = ellipsoidRadiusInCameraDir(
    camDirX, camDirY, camDirZ
  );

  // カメラが楕円体内部にある場合はカリングしない
  if (cameraDistance <= R) return true;

  // タイル中心のカメラ方向への射影（内積）
  const proj = tileCenter.x * camDirX
             + tileCenter.y * camDirY
             + tileCenter.z * camDirZ;

  // ホライズンカットオフ距離: R^2 / d
  const horizonCut = R * R * invD;

  // バウンディングスフィア全体が裏側なら不可視
  return proj + tileRadius >= horizonCut;
}

        d
*-----------------* カメラ
        | R^2/d
        *--- 接点
       /
      / R（カメラ方向の楕円体半径）
     /
    * 原点

export function tileBoundingSphereRadius(
  x: number, y: number, z: number
): number {
  const bounds = tileBounds(x, y, z);
  const latCenter = (bounds.north + bounds.south) / 2;
  const cosLat = Math.cos(degToRad(latCenter));
  return ((tileGeometricError(z) * cosLat) / 2) * Math.SQRT2;
}

/** 球がFrustum内にあるかチェック */
function sphereInFrustum(
  center: { x: number; y: number; z: number },
  radius: number,
  planes: FrustumPlane[]
): boolean {
  for (const plane of planes) {
    const dist =
      plane.normal.x * center.x +
      plane.normal.y * center.y +
      plane.normal.z * center.z +
      plane.constant;
    if (dist < -radius) return false;
  }
  return true;
}

export function extractFrustumPlanes(
  camera: THREE.PerspectiveCamera
): FrustumPlane[] {
  const frustum = new THREE.Frustum();
  const matrix = new THREE.Matrix4().multiplyMatrices(
    camera.projectionMatrix,
    camera.matrixWorldInverse
  );
  frustum.setFromProjectionMatrix(matrix);

  return frustum.planes.map((plane) => ({
    normal: { x: plane.normal.x, y: plane.normal.y, z: plane.normal.z },
    constant: plane.constant
  }));
}

z=0 --- SSE高い --> 分割
  z=1 -+-- 日本含む象限: SSE高い --> 分割
        |    z=2 --- 日本周辺: SSE高い --> さらに分割...
        |           z=3 --- 東京: SSEまだ高い --> 分割
        |                  z=4 --- SSE <= 閾値 --> 採用!
        +-- 大西洋: Frustum外 --> カリング
        +-- 南極: SSE低い --> z=1で採用
        +-- 太平洋: SSE低い --> z=1で採用

const tiles = selectVisibleTiles(
  cameraPos,       // カメラ位置
  fovRad,          // FOV（ラジアン）
  containerHeight, // 描画領域の高さ(px)
  frustumPlanes,   // 6平面
  400,             // SSE閾値(px)
  6,               // 最大ズーム
  128              // 最大タイル数
);

// 毎フレームの処理
camera.updateMatrixWorld();

if (cameraMatrixChanged()) {
  const needed = selectVisibleTiles(...);

  // 不要になったタイルを削除+dispose
  for (const [key, mesh] of activeTiles) {
    if (!needed.has(key)) {
      scene.remove(mesh);
      mesh.geometry.dispose();
      mesh.material.dispose();
    }
  }

  // 新規タイルを生成+追加
  for (const tile of needed) {
    if (!activeTiles.has(key)) {
      // computeTileGeometry + シーン追加
    }
  }
}