Move curvature calculations to Ellipsoid

CesiumGS · Dec 1, 2023 · 2448e52 · 2448e52
1 parent 3e63f3a
commit 2448e52
Show file tree

Hide file tree

Showing 3 changed files with 92 additions and 23 deletions.
diff --git a/packages/engine/Source/Core/Ellipsoid.js b/packages/engine/Source/Core/Ellipsoid.js
@@ -1,3 +1,4 @@
+import Cartesian2 from "./Cartesian2.js";
 import Cartesian3 from "./Cartesian3.js";
 import Cartographic from "./Cartographic.js";
 import Check from "./Check.js";
@@ -694,6 +695,50 @@ Ellipsoid.prototype.getSurfaceNormalIntersectionWithZAxis = function (
   return result;
 };
 
+const scratchEndpoint = new Cartesian3();
+
+/**
+ * Computes the ellipsoid curvatures at a given position on the surface.
+ *
+ * @param {Cartesian3} surfacePosition The position on the ellipsoid surface where curvatures will be calculated.
+ * @param {Cartesian2} [result] The cartesian to which to copy the result, or undefined to create and return a new instance.
+ * @returns {Cartesian2} The local curvature of the ellipsoid surface at the provided position, in east and north directions.
+ *
+ * @exception {DeveloperError} position is required.
+ */
+Ellipsoid.prototype.getLocalCurvature = function (surfacePosition, result) {
+  //>>includeStart('debug', pragmas.debug);
+  Check.typeOf.object("surfacePosition", surfacePosition);
+  //>>includeEnd('debug');
+
+  if (!defined(result)) {
+    result = new Cartesian2();
+  }
+
+  const primeVerticalEndpoint = this.getSurfaceNormalIntersectionWithZAxis(
+    surfacePosition,
+    0.0,
+    scratchEndpoint
+  );
+  const primeVerticalRadius = Cartesian3.distance(
+    surfacePosition,
+    primeVerticalEndpoint
+  );
+  // meridional radius = (1 - e^2) * primeVerticalRadius^3 / a^2
+  // where 1 - e^2 = b^2 / a^2,
+  // so meridional = b^2 * primeVerticalRadius^3 / a^4
+  //   = (b * primeVerticalRadius / a^2)^2 * primeVertical
+  const radiusRatio =
+    (this.minimumRadius * primeVerticalRadius) / this.maximumRadius ** 2;
+  const meridionalRadius = primeVerticalRadius * radiusRatio ** 2;
+
+  return Cartesian2.fromElements(
+    1.0 / primeVerticalRadius,
+    1.0 / meridionalRadius,
+    result
+  );
+};
+
 const abscissas = [
   0.14887433898163,
   0.43339539412925,

diff --git a/packages/engine/Source/Renderer/UniformState.js b/packages/engine/Source/Renderer/UniformState.js
@@ -1119,7 +1119,6 @@ function setInfiniteProjection(uniformState, matrix) {
 }
 
 const surfacePositionScratch = new Cartesian3();
-const primeVerticalScratch = new Cartesian3();
 const enuTransformScratch = new Matrix4();
 
 function setCamera(uniformState, camera) {
@@ -1197,26 +1196,8 @@ function setCamera(uniformState, camera) {
     return;
   }
 
-  const primeVerticalEndpoint = ellipsoid.getSurfaceNormalIntersectionWithZAxis(
+  uniformState._eyeEllipsoidCurvature = ellipsoid.getLocalCurvature(
     surfacePosition,
-    0.0,
-    primeVerticalScratch
-  );
-  const primeVerticalRadius = Cartesian3.distance(
-    surfacePosition,
-    primeVerticalEndpoint
-  );
-  // meridional radius = (1 - e^2) * primeVerticalRadius^3 / a^2
-  // where 1 - e^2 = b^2 / a^2,
-  // so meridional = b^2 * primeVerticalRadius^3 / a^4
-  //   = (b * primeVerticalRadius / a^2)^2 * primeVertical
-  const radiusRatio =
-    (ellipsoid.minimumRadius * primeVerticalRadius) /
-    ellipsoid.maximumRadius ** 2;
-  const meridionalRadius = primeVerticalRadius * radiusRatio ** 2;
-  uniformState._eyeEllipsoidCurvature = Cartesian2.fromElements(
-    1.0 / primeVerticalRadius,
-    1.0 / meridionalRadius,
     uniformState._eyeEllipsoidCurvature
   );
 }

diff --git a/packages/engine/Specs/Core/EllipsoidSpec.js b/packages/engine/Specs/Core/EllipsoidSpec.js
@@ -1,6 +1,11 @@
-import { Cartesian3, Cartographic, Ellipsoid, Rectangle } from "../../index.js";
-
-import { Math as CesiumMath } from "../../index.js";
+import {
+  Cartesian2,
+  Cartesian3,
+  Cartographic,
+  Ellipsoid,
+  Rectangle,
+  Math as CesiumMath,
+} from "../../index.js";
 
 import createPackableSpecs from "../../../../Specs/createPackableSpecs.js";
 
@@ -721,6 +726,44 @@ describe("Core/Ellipsoid", function () {
     );
   });
 
+  it("getLocalCurvature throws with no position", function () {
+    expect(function () {
+      Ellipsoid.WGS84.getLocalCurvature(undefined);
+    }).toThrowDeveloperError();
+  });
+
+  it("getLocalCurvature returns expected values at the equator", function () {
+    const ellipsoid = Ellipsoid.WGS84;
+    const cartographic = Cartographic.fromDegrees(0.0, 0.0);
+    const cartesianOnTheSurface = ellipsoid.cartographicToCartesian(
+      cartographic
+    );
+    const returnedResult = ellipsoid.getLocalCurvature(cartesianOnTheSurface);
+    const expectedResult = new Cartesian2(
+      1.0 / ellipsoid.maximumRadius,
+      ellipsoid.maximumRadius /
+        (ellipsoid.minimumRadius * ellipsoid.minimumRadius)
+    );
+    expect(returnedResult).toEqualEpsilon(expectedResult, CesiumMath.EPSILON8);
+  });
+
+  it("getLocalCurvature returns expected values at the north pole", function () {
+    const ellipsoid = Ellipsoid.WGS84;
+    const cartographic = Cartographic.fromDegrees(0.0, 90.0);
+    const cartesianOnTheSurface = ellipsoid.cartographicToCartesian(
+      cartographic
+    );
+    const returnedResult = ellipsoid.getLocalCurvature(cartesianOnTheSurface);
+    const semiLatusRectum =
+      (ellipsoid.maximumRadius * ellipsoid.maximumRadius) /
+      ellipsoid.minimumRadius;
+    const expectedResult = new Cartesian2(
+      1.0 / semiLatusRectum,
+      1.0 / semiLatusRectum
+    );
+    expect(returnedResult).toEqualEpsilon(expectedResult, CesiumMath.EPSILON8);
+  });
+
   it("ellipsoid is initialized with _squaredXOverSquaredZ property", function () {
     const ellipsoid = new Ellipsoid(4, 4, 3);