Fix numerical error in voxel latitude intersection

CesiumGS · Nov 15, 2023 · a69ebee · a69ebee
1 parent c744988
commit a69ebee
Show file tree

Hide file tree

Showing 3 changed files with 10 additions and 6 deletions.
diff --git a/packages/engine/Source/Scene/VoxelEllipsoidShape.js b/packages/engine/Source/Scene/VoxelEllipsoidShape.js
@@ -96,6 +96,7 @@ function VoxelEllipsoidShape() {
    */
   this.shaderUniforms = {
     ellipsoidRadiiUv: new Cartesian3(),
+    eccentricitySquared: 0.0,
     ellipsoidInverseRadiiSquaredUv: new Cartesian3(),
     ellipsoidRenderLongitudeMinMax: new Cartesian2(),
     ellipsoidShapeUvLongitudeMinMaxMid: new Cartesian3(),
@@ -427,6 +428,9 @@ VoxelEllipsoidShape.prototype.update = function (
     shapeMaxExtent,
     shaderUniforms.ellipsoidRadiiUv
   );
+  const axisRatio =
+    Cartesian3.minimumComponent(shapeOuterExtent) / shapeMaxExtent;
+  shaderUniforms.eccentricitySquared = 1.0 - axisRatio * axisRatio;
 
   // Used to compute geodetic surface normal.
   shaderUniforms.ellipsoidInverseRadiiSquaredUv = Cartesian3.divideComponents(

diff --git a/packages/engine/Source/Shaders/Voxels/IntersectEllipsoid.glsl b/packages/engine/Source/Shaders/Voxels/IntersectEllipsoid.glsl
@@ -24,6 +24,7 @@
 #if defined(ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE)
     uniform vec2 u_ellipsoidRenderLongitudeMinMax;
 #endif
+uniform float u_eccentricitySquared;
 uniform vec2 u_ellipsoidRenderLatitudeCosHalfMinMax;
 uniform vec2 u_clipMinMaxHeight;
 
@@ -185,13 +186,12 @@ vec3 getConeNormal(in vec3 p, in bool convex) {
 float getLatitudeConeShift(in float sinLatitude) {
     // Find prime vertical radius of curvature: 
     // the distance along the ellipsoid normal to the intersection with the z-axis
-    float axisRatio = u_ellipsoidRadiiUv.z / u_ellipsoidRadiiUv.x;
-    float eccentricitySquared = 1.0 - axisRatio * axisRatio;
-    //float eccentricitySquared = 6.69437999014e-3; // ASSUMES WGS84. Supply as uniform?
-    float primeVerticalRadius = inversesqrt(1.0 - eccentricitySquared * sinLatitude * sinLatitude);
+    float x2 = u_eccentricitySquared * sinLatitude * sinLatitude;
+    //float primeVerticalRadius = inversesqrt(1.0 - x2);
+    float primeVerticalRadius = 1.0 + x2 * (0.5 + x2 * (0.375 + 0.3125 * x2));
 
     // Compute a shift from the origin to the intersection of the cone with the z-axis
-    return primeVerticalRadius * eccentricitySquared * sinLatitude;
+    return primeVerticalRadius * u_eccentricitySquared * sinLatitude;
 }
 
 /**

diff --git a/packages/engine/Source/Shaders/Voxels/convertUvToEllipsoid.glsl b/packages/engine/Source/Shaders/Voxels/convertUvToEllipsoid.glsl
@@ -224,7 +224,7 @@ VoxelCell convertShapeUvToShapeSpace(in VoxelCell voxel) {
         dP.y = dP.y / u_ellipsoidUvToShapeUvLatitude.x;
     #endif
     // Convert from [0, 1] to radians [-pi/2, pi/2]
-    p.y = p.y * czm_pi - czm_piOverTwo - 0.000006;
+    p.y = p.y * czm_pi - czm_piOverTwo;
     dP.y = dP.y * czm_pi;
 
     #if defined(ELLIPSOID_HAS_SHAPE_BOUNDS_HEIGHT_FLAT)