Cesium3DTileset地面阴影面积计算：逻辑与代码疑问

我想咨询如何计算Cesium3DTileset在地面的阴影面积，当前我遵循以下逻辑：

1. 获取Cesium3DTileset的BoundingSphere。
2. 使用Simon1994PlanetaryPositions计算太阳位置。
3. 计算从太阳指向Cesium3DTileset中心的光线向量。
4. 计算光线方向向量与物体的交点并投影到地面。
5. 计算这些点构成的多边形面积。

const calculate3DTilesShadowArea = (
  tileset: Cesium.Cesium3DTileset, julianDate: Cesium.JulianDate, viewer: Cesium.Viewer,
) => {
  // Get BoundingSphere
  const boundingSphere = tileset.boundingSphere;

  // Get Sun Position
  const sunPosition = Cesium.Simon1994PlanetaryPositions.computeSunPositionInEarthInertialFrame(julianDate);

  // Vector light direction to bounding sphere center
  const lightDirection = Cesium.Cartesian3.normalize(
    Cesium.Cartesian3.subtract(
      sunPosition,
      boundingSphere.center,
      new Cesium.Cartesian3(),
    ),
    new Cesium.Cartesian3(),
  );

  // Create a matrix to transform points from the local space of the bounding sphere to world space
  const transformMatrix = Cesium.Transforms.eastNorthUpToFixedFrame(boundingSphere.center);

  // Generate points on the surface of the bounding sphere to simulate the shadow outline
  const numPoints = 64; // Adjust for accuracy
  let boundingVolumePoints = [];
  for (let i = 0; i < numPoints; i++) {
    const angle = (i / numPoints) * Cesium.Math.TWO_PI;
    const x = boundingSphere.radius * Math.cos(angle);
    const y = boundingSphere.radius * Math.sin(angle);
    const z = 0;

    // Transform the point to world coordinates
    const point = Cesium.Matrix4.multiplyByPoint(transformMatrix, new Cesium.Cartesian3(x, y, z), new Cesium.Cartesian3());

    // Project the point in the direction of the light onto the ground plane (z = 0)
    const t = -point.z / lightDirection.z;
    const shadowPoint = new Cesium.Cartesian3(
      point.x + t * lightDirection.x,
      point.y + t * lightDirection.y,
      0, // z = 0 for ground plane
    );

    boundingVolumePoints.push(new Cesium.Cartesian2(shadowPoint.x, shadowPoint.y));
  }
  if (boundingVolumePoints.length < 3) {
    console.error('Not enough points for shadow calculation');
    return 0;
  }

  const area = computePolygonArea(boundingVolumePoints);

  return area;
};

请问我的逻辑或代码存在哪些问题？

• 坐标系不匹配导致光线方向错误
  Simon1994PlanetaryPositions.computeSunPositionInEarthInertialFrame返回的太阳位置属于地球惯性坐标系（ECI），而tileset.boundingSphere.center是**地球固定坐标系（ECEF）**下的坐标。直接对不同坐标系的点做向量运算，得到的光线方向完全不符合实际，必须先通过Cesium.Transforms.transformPoint将太阳位置转换到ECEF坐标系后再计算。

• BoundingSphere采样范围严重不足
  代码仅在BoundingSphere的局部z=0平面采样点，完全没有覆盖球体的上下半球。当太阳从顶部、底部或斜向照射时，这些采样点无法反映物体的真实轮廓，会导致阴影面积计算出现极大偏差。正确做法是在球体表面均匀采样（比如覆盖多个纬度圈），或者使用BoundingSphere的外接正多面体顶点生成轮廓点。

• 地面平面假设不符合Cesium的地球模型
  代码将地面视为z=0的平面，但Cesium使用的是地球椭球面模型。在高纬度或远离球心的区域，平面投影会产生严重失真，应该计算光线与地球椭球面的交点来得到阴影的地面位置。

• 光线方向定义错误
  太阳距离地球极远，阴影投影应使用平行光模型，光线方向应为指向太阳的反方向。当前代码计算的是从物体中心指向太阳的向量，正确的平行光方向应该取该向量的反方向（或直接使用太阳的单位向量）。

• 面积计算的坐标系未转换
  直接用ECEF坐标系下的Cartesian2点计算面积，结果没有实际地理意义。需要将投影后的点转换为大地坐标系（经纬度），再转成平面投影坐标系（如UTM）后计算面积，否则会因地球曲率和坐标系缩放导致结果错误。

内容的提问来源于stack exchange，提问作者user25091635