Model polygon clipping shader impl

packages/engine/Source/Renderer/AutomaticUniforms.js

@@ -1008,6 +1008,27 @@ const AutomaticUniforms = {
     },
   }),
 
+  /**
+   * An automatic GLSL uniform containing the geodetic longitude (<code>x</code>),
+   * latitude (<code>y</code>) in radians and height (<code>z</code>) in meters of the
+   * eye (camera). The <code>z</code> component matches {@link czm_eyeHeight}.
+   * This uniform is only valid when the {@link SceneMode} is <code>SCENE3D</code>.
+   *
+   * @example
+   * // GLSL declaration
+   * uniform vec3 czm_eyeCartographic;
+   *
+   * // Example
+   * float cameraLatitude = czm_eyeCartographic.y;
+   */
+  czm_eyeCartographic: new AutomaticUniform({
+    size: 1,
+    datatype: WebGLConstants.FLOAT_VEC3,
+    getValue: function (uniformState) {
+      return uniformState.eyeCartographic;
+    },
+  }),
+
   /**
    * An automatic GLSL uniform containing height (<code>x</code>) and height squared (<code>y</code>)
    * in meters of the eye (camera) above the 2D world plane. This uniform is only valid
@@ -1050,6 +1071,27 @@ const AutomaticUniforms = {
     },
   }),
 
+  /**
+   * An automatic GLSL uniform containing a 3x3 rotation from eye coordinates to an
+   * east-north-up coordinate system centered at the position on the ellipsoid below
+   * the camera.
+   * This uniform is only valid when the {@link SceneMode} is <code>SCENE3D</code>.
+   *
+   * @example
+   * // GLSL declaration
+   * uniform mat3 czm_eyeToEnu;
+   *
+   * // Example
+   * vec3 enu = czm_eyeToEnu * positionEC;
+   */
+  czm_eyeToEnu: new AutomaticUniform({
+    size: 1,
+    datatype: WebGLConstants.FLOAT_MAT3,
+    getValue: function (uniformState) {
+      return uniformState.eyeToEnu;
+    },
+  }),
+
   /**
    * An automatic GLSL uniform containing the transform from model coordinates
    * to an east-north-up coordinate system centered at the position on the

packages/engine/Source/Renderer/UniformState.js

@@ -162,8 +162,10 @@ function UniformState() {
   this._frustum2DWidth = 0.0;
   this._eyeHeight = 0.0;
   this._eyeHeight2D = new Cartesian2();
+  this._eyeCartographic = new Cartesian3();
   this._eyeEllipsoidNormalEC = new Cartesian3();
   this._eyeEllipsoidCurvature = new Cartesian2();
+  this._eyeToEnu = new Matrix3();
   this._modelToEnu = new Matrix4();
   this._enuToModel = new Matrix4();
   this._pixelRatio = 1.0;
@@ -715,6 +717,19 @@ Object.defineProperties(UniformState.prototype, {
     },
   },
 
+  /**
+   * The geodetic longitude (<code>x</code>), latitude (<code>y</code>) in radians and
+   * height (<code>z</code>) in meters of the eye (camera).
+   * This is only valid when the {@link SceneMode} is <code>SCENE3D</code>.
+   * @memberof UniformState.prototype
+   * @type {Cartesian3}
+   */
+  eyeCartographic: {
+    get: function () {
+      return this._eyeCartographic;
+    },
+  },
+
   /**
    * The height (<code>x</code>) and the height squared (<code>y</code>)
    * in meters of the eye (camera) above the 2D world plane. This uniform is only valid
@@ -751,6 +766,19 @@ Object.defineProperties(UniformState.prototype, {
     },
   },
 
+  /**
+   * A 3x3 rotation from eye coordinates to an east-north-up coordinate system
+   * centered at the position on the ellipsoid below the camera.
+   * This uniform is only valid when the {@link SceneMode} is <code>SCENE3D</code>.
+   * @memberof UniformState.prototype
+   * @type {Matrix3}
+   */
+  eyeToEnu: {
+    get: function () {
+      return this._eyeToEnu;
+    },
+  },
+
   /**
    * A transform from model coordinates to an east-north-up coordinate system
    * centered at the position on the ellipsoid below the camera
@@ -1241,6 +1269,7 @@ function setInfiniteProjection(uniformState, matrix) {
 
 const surfacePositionScratch = new Cartesian3();
 const enuTransformScratch = new Matrix4();
+const enuRotationScratch = new Matrix3();
 
 function setCamera(uniformState, camera) {
   Cartesian3.clone(camera.positionWC, uniformState._cameraPosition);
@@ -1266,6 +1295,12 @@ function setCamera(uniformState, camera) {
     );
   } else {
     uniformState._eyeHeight = positionCartographic.height;
+    uniformState._eyeCartographic = Cartesian3.fromElements(
+      positionCartographic.longitude,
+      positionCartographic.latitude,
+      positionCartographic.height,
+      uniformState._eyeCartographic,
+    );
     uniformState._eyeEllipsoidNormalEC =
       ellipsoid.geodeticSurfaceNormalCartographic(
         positionCartographic,
@@ -1307,6 +1342,17 @@ function setCamera(uniformState, camera) {
     uniformState._modelToEnu,
   );
 
+  const enuToWorldRotation = Matrix4.getRotation(
+    enuToWorld,
+    enuRotationScratch,
+  );
+  const enuToView = Matrix3.multiply(
+    uniformState._viewRotation,
+    enuToWorldRotation,
+    enuRotationScratch,
+  );
+  uniformState._eyeToEnu = Matrix3.transpose(enuToView, uniformState._eyeToEnu);
+
   if (
     !CesiumMath.equalsEpsilon(
       ellipsoid._radii.x,

packages/engine/Source/Scene/ClippingPolygonCollection.js

@@ -330,7 +330,7 @@ Object.defineProperties(ClippingPolygonCollection.prototype, {
   /**
    * Returns a number encapsulating the state for this ClippingPolygonCollection.
    *
-   * Clipping mode is encoded in the sign of the number, which is just the total position count.
+   * Clipping mode is encoded in the sign of the number, and the magnitude is the polygon count.
    * If this value changes, then shader regeneration is necessary.
    *
    * @memberof ClippingPolygonCollection.prototype
@@ -340,7 +340,7 @@ Object.defineProperties(ClippingPolygonCollection.prototype, {
    */
   clippingPolygonsState: {
     get: function () {
-      return this.inverse ? -this.extentsCount : this.extentsCount;
+      return this.inverse ? -this.length : this.length;
     },
   },
 });

packages/engine/Source/Scene/Model/Model.js

@@ -2633,6 +2633,22 @@ function updateShowCreditsOnScreen(model) {
   }
 }
 
+/**
+ * Determines whether a model is wholly clipped away by inverse clipping. In
+ * inverse mode a model with no polygon geometry lies entirely outside every
+ * polygon, so all of it is clipped and it should not be drawn.
+ * @param {Model} model
+ * @returns {boolean}
+ * @private
+ */
+function isModelClippedAwayByInversePolygons(model) {
+  return (
+    model.isClippingPolygonsEnabled() &&
+    model._clippingPolygons.inverse &&
+    (model._clippingPolygonData?.polygonRings.length ?? 0) === 0
+  );
+}
+
 function submitDrawCommands(model, frameState) {
   // Check that show is true after draw commands are built;
   // we want the user to be able to instantly see the model
@@ -2653,7 +2669,8 @@ function submitDrawCommands(model, frameState) {
     model._show &&
     model._computedScale !== 0 &&
     displayConditionPassed &&
-    (!invisible || silhouette);
+    (!invisible || silhouette) &&
+    !isModelClippedAwayByInversePolygons(model);
 
   const passes = frameState.passes;
   const submitCommandsForPass =

packages/engine/Source/Scene/Model/ModelClippingPolygonsPipelineStage.js

@@ -1,7 +1,11 @@
 import combine from "../../Core/combine.js";
-import ModelClippingPolygonsStageVS from "../../Shaders/Model/ModelClippingPolygonsStageVS.js";
 import ModelClippingPolygonsStageFS from "../../Shaders/Model/ModelClippingPolygonsStageFS.js";
+import ModelClippingPolygonsStageVS from "../../Shaders/Model/ModelClippingPolygonsStageVS.js";
 import ShaderDestination from "../../Renderer/ShaderDestination.js";
+import VectorCommon from "../../Shaders/VectorCommon.js";
+import Cartesian2 from "../../Core/Cartesian2.js";
+import CesiumMath from "../../Core/Math.js";
+import Rectangle from "../../Core/Rectangle.js";
 
 /**
  * The model clipping planes stage is responsible for applying clipping planes to the model.
@@ -14,6 +18,10 @@ const ModelClippingPolygonsPipelineStage = {
   name: "ModelClippingPolygonsPipelineStage", // Helps with debugging
 };
 
+const scratchCameraUv = new Cartesian2();
+const scratchRectangleInverseSize = new Cartesian2();
+const defaultRectangle = Rectangle.MAX_VALUE;
+
 /**
  * Process a model for polygon clipping. This modifies the following parts of the render resources:
  *
@@ -54,40 +62,67 @@ ModelClippingPolygonsPipelineStage.process = function (
     );
   }
 
-  shaderBuilder.addDefine(
-    "CLIPPING_POLYGON_REGIONS_LENGTH",
-    clippingPolygons.extentsCount,
-    ShaderDestination.BOTH,
-  );
+  shaderBuilder.addVarying("vec2", "v_clippingUv");
 
   shaderBuilder.addUniform(
-    "sampler2D",
-    "model_clippingDistance",
-    ShaderDestination.FRAGMENT,
+    "vec2",
+    "u_clippingCameraUv",
+    ShaderDestination.VERTEX,
   );
 
   shaderBuilder.addUniform(
-    "sampler2D",
-    "model_clippingExtents",
+    "vec2",
+    "u_clippingRectangleInverseSize",
     ShaderDestination.VERTEX,
   );
 
-  shaderBuilder.addVarying("vec2", "v_clippingPosition");
-  shaderBuilder.addVarying("int", "v_regionIndex", "flat");
   shaderBuilder.addVertexLines(ModelClippingPolygonsStageVS);
+  shaderBuilder.addFragmentLines(VectorCommon);
   shaderBuilder.addFragmentLines(ModelClippingPolygonsStageFS);
 
   const uniformMap = {
-    model_clippingDistance: function () {
+    // The UV coordinates of the camera within the model's clipping rectangle.
+    u_clippingCameraUv: function () {
+      const rectangle =
+        model._clippingPolygonData?.rectangle ?? defaultRectangle;
+      const halfWidth = rectangle.width * 0.5;
+      const centerLongitude = rectangle.west + halfWidth;
+      const carto = frameState.camera.positionCartographic;
+
+      const longitudeOffset =
+        CesiumMath.negativePiToPi(carto.longitude - centerLongitude) +
+        halfWidth;
+      return Cartesian2.fromElements(
+        longitudeOffset / rectangle.width,
+        (carto.latitude - rectangle.south) / rectangle.height,
+        scratchCameraUv,
+      );
+    },
+    u_clippingRectangleInverseSize: function () {
+      const rectangle =
+        model._clippingPolygonData?.rectangle ?? defaultRectangle;
+      return Cartesian2.fromElements(
+        1.0 / rectangle.width,
+        1.0 / rectangle.height,
+        scratchRectangleInverseSize,
+      );
+    },
+    u_clippingEdgeTexture: function () {
+      return (
+        model._clippingPolygonData?.polygonEdgeTexture ??
+        frameState.context.defaultTexture
+      );
+    },
+    u_clippingEdgePrimitiveIndicesTexture: function () {
       return (
-        // The later should never happen during a render pass, see https://github.com/CesiumGS/cesium/issues/12725
-        clippingPolygons.clippingTexture ?? frameState.context.defaultTexture
+        model._clippingPolygonData?.polygonEdgePrimitiveIndicesTexture ??
+        frameState.context.defaultTexture
       );
     },
-    model_clippingExtents: function () {
+    u_clippingGridCellIndicesTexture: function () {
       return (
-        // The later should never happen during a render pass, see https://github.com/CesiumGS/cesium/issues/12725
-        clippingPolygons.extentsTexture ?? frameState.context.defaultTexture
+        model._clippingPolygonData?.polygonGridCellIndicesTexture ??
+        frameState.context.defaultTexture
       );
     },
   };

packages/engine/Source/Shaders/Builtin/Functions/eyeToCartographicDelta.glsl

@@ -0,0 +1,74 @@
+/**
+ * Computes the geodetic offset (delta longitude, latitude, and height) from a reference cartographic position
+ * to a point given in eye coordinates.
+ *
+ * This is designed to preserve precision. Rather than converting the point's absolute world position to
+ * cartographic -- which is too large to process precisely at 32 bits -- it works entirely with the small
+ * delta between the point and the camera. By projecting the eye-space offset onto the ellipsoid's equatorial
+ * and meridional planes, it derives the change in (longitude, latitude, height) as small, precisely-representable
+ * quantities. The delta gets smaller and more precise as one zooms in.
+ * <br /><br />
+ * This assumes an ellipsoid of revolution (equatorial radii equal, as with WGS84), so that longitude is exact.
+ * The latitude calculation is only first-order accurate, since the meridian is an ellipse rather than a circle.
+ *
+ * @name czm_eyeToCartographicDelta
+ * @glslFunction
+ *
+ * @param {vec3} positionEC The position, in eye coordinates, to measure to.
+ *
+ * @returns {vec3} The geodetic offset from the camera to <code>positionEC</code>, as (delta longitude, delta latitude in radians, delta height in meters).
+ */
+vec3 czm_eyeToCartographicDelta(vec3 positionEC)
+{
+    // A vector representing the camera-to-vertex offset, in an ENU oriented reference frame (centered at the camera)
+    vec3 cameraToVertex = czm_eyeToEnu * positionEC;
+
+    float cosLatitude = cos(czm_eyeCartographic.y);
+    float sinLatitude = sin(czm_eyeCartographic.y);
+
+    // To derive longitude, project the camera and vertex onto the equatorial plane, in a frame such that the camera lies along the +x axis. In this frame,
+    // the vertex's (delta) longitude is simply the atan of its x and y components.
+    float primeVerticalRadius = 1.0 / czm_eyeEllipsoidCurvature.x;
+    vec2 cameraEquatorialPos = vec2((primeVerticalRadius + czm_eyeCartographic.z) * cosLatitude, 0.0);
+    vec2 vertexEquatorialPos = cameraEquatorialPos + vec2(-cameraToVertex.y * sinLatitude + cameraToVertex.z * cosLatitude, cameraToVertex.x);
+    float deltaLongitude = atan(vertexEquatorialPos.y, vertexEquatorialPos.x);
+
+    // Deriving latitude is a bit harder: we can't directly project the vertex onto the camera's meridian — the latitude projection is dependent on the longitude.
+    // Instead we can rotate the vertex (by -deltaLongitude) onto the camera's meridional plane.  (Note: (unlike the exact longitude case) this is only first-order accurate because the meridian is an ellipse rather than a circle)
+    // Using a 2D rotation formula introduces precision issues (subtraction of large-magnitude quantities), so instead we can calculate the vector difference
+    // between the vertex and its rotated version, and apply that offset to the cameraToVertex vector. Then, the cameraToVertex vector accurately
+    // reflects the difference between the camera and the _rotated_ vertex, so we can then project the camera onto the meridional plane and apply this offset - just as we did for deltaLongitude, above.
+    // Best of all, we can do this all with small delta quantities which preserve precision.
+    //
+    // (I suggest drawing this out -- with the vertex and camera vectors projected onto the equatorial plane, with the camera on the +x axis)
+    // Mathematically: if you compare (subtract) vertexEquatorialPos and the same vector rotated onto the camera's meridional plane, you get
+    // |dx| = |vertexEquatorialPos| - vertexEquatorialPos.x = (r - x) = r * (1 - cos(deltaLongitude))
+    // |dy| = cameraToVertex.x (the east component)
+    // (To avoid precision issues, we'll use the identity (1 - cos(x) = 2 * sin^2(x/2)))
+    //
+    // Since these offsets were produced in the equatorial plane, and cameraToVertex is in the camera's ENU frame, we need to deconstruct along the camera's north and up axes. And we only care about
+    // dx, since dy is in the camera's east direction, and that component gets zeroed out when projecting onto the camera's meridional plane.
+    float sinHalfLongitude = sin(deltaLongitude * 0.5);
+    float dx = length(vertexEquatorialPos) * 2.0 * sinHalfLongitude * sinHalfLongitude;
+    vec3 meridionalOffset = vec3(
+        0.0,                                 // east
+        cameraToVertex.y - dx * sinLatitude, // north
+        cameraToVertex.z + dx * cosLatitude  // up
+    );
+
+    // Reframe the camera in a meridional plane, where it lies along the +z axis, and apply the meridionalOffset to get the vertex's position in that plane.
+    // Then, deltaLatitude is simply the atan of its x and y components.
+    float meridionalRadius = 1.0 / czm_eyeEllipsoidCurvature.y;
+    vec2 cameraMeridionalPos = vec2(meridionalRadius + czm_eyeCartographic.z, 0.0);
+    vec2 vertMeridionalPos = cameraMeridionalPos + vec2(meridionalOffset.z, meridionalOffset.y);
+    float deltaLatitude = atan(vertMeridionalPos.y, vertMeridionalPos.x);
+
+    // Finally, derive the change in height above the ellipsoid. This is the meridional-plane analogue of the dx step above:
+    // there we rotated the vertex (in the equatorial plane) to the camera's longitude; here we rotate it (in the meridional plane, by -deltaLatitude)
+    // to the camera's latitude, aligning it with the camera's radial (up) direction.
+    float sinHalfLatitude = sin(deltaLatitude * 0.5);
+    float dz = length(vertMeridionalPos) * 2.0 * sinHalfLatitude * sinHalfLatitude;
+    float deltaHeight = meridionalOffset.z + dz;
+
+    return vec3(deltaLongitude, deltaLatitude, deltaHeight);
+}