Fix 3D model query intersection inaccuracies

GitOrigin-RevId: 621eaabde48d77614f1f58350eb1d6c40c3f3c08

Author: stephanheigl

3d-style/style/style_layer/model_style_layer.ts

@@ -101,7 +101,7 @@ class ModelStyleLayer extends StyleLayer {
                 const model = modelManager.getModel(modelId, this.scope);
                 if (!model) return false;
 
-                let matrix = mat4.create();
+                let matrix: mat4 = [];
                 const position = new LngLat(0, 0);
                 const id = bucket.canonical;
                 let minDepth = Number.MAX_VALUE;
@@ -111,8 +111,8 @@ class ModelStyleLayer extends StyleLayer {
 
                     const va = instances.instancedDataArray.float32;
                     const translation: [number, number, number] = [va[offset + 4], va[offset + 5], va[offset + 6]];
-                    const pointX = va[offset];
-                    const pointY = va[offset + 1] | 0; // point.y stored in integer part
+                    const pointX = Math.floor(va[offset]); // point.x stored in integer part
+                    const pointY = Math.floor(va[offset + 1]); // point.y stored in integer part
 
                     tileToLngLat(id, position, pointX, pointY);
 

3d-style/util/model_util.ts

@@ -20,6 +20,7 @@ import {Aabb} from '../../src/util/primitives';
 import {polygonIntersectsPolygon} from '../../src/util/intersection_tests';
 import Point from '@mapbox/point-geometry';
 
+import type {vec2} from 'gl-matrix';
 import type Transform from '../../src/geo/transform';
 
 export function rotationScaleYZFlipMatrix(out: mat4, rotation: vec3, scale: vec3) {
@@ -214,6 +215,79 @@ export function convertModelMatrixForGlobe(matrix: mat4, transform: Transform, s
     return modelMatrix;
 }
 
+/**
+ * Computes the convex hull from points using Jarvis's algorithm (Gift Wrapping).
+ * The algorithm finds the leftmost point first and then iterates through
+ * outermost points until the hull is closed.
+ *
+ * @param points Array of 2D points to compute convex hull for
+ * @returns Array of Points forming the convex hull as a closed linear ring,
+ *          or empty array if input has fewer than 3 points
+ * @private
+ */
+// Ported from gl-native: src/mbgl/util/bounding_volumes.cpp - convexPolygonFromPoints
+export function convexPolygonFromPoints(points: Array<vec2>): Array<Point> {
+    assert(points.length >= 3, 'Polygon must have at least 3 points');
+
+    const ring: Array<Point> = [];
+
+    // Find the leftmost and bottommost point first
+    let leftmostIdx = 0;
+    for (let i = 1; i < points.length; i++) {
+        if (points[i][0] < points[leftmostIdx][0] ||
+            (points[i][0] === points[leftmostIdx][0] && points[i][1] < points[leftmostIdx][1])) {
+            leftmostIdx = i;
+        }
+    }
+
+    let current = leftmostIdx;
+    let next: number;
+    const visitedMap = new Uint8Array(points.length);
+
+    do {
+        // Return if the previous point is revisited, otherwise will end up with infinite loop
+        if (visitedMap[current]) {
+            break;
+        }
+
+        ring.push(new Point(points[current][0], points[current][1]));
+        visitedMap[current] = 1;
+
+        // Find the next most outmost point
+        next = (current + 1) % points.length;
+        for (let i = 0; i < points.length; i++) {
+            if ((points[i][0] === points[next][0] && points[i][1] === points[next][1]) ||
+                (points[i][0] === points[current][0] && points[i][1] === points[current][1])) {
+                continue;
+            }
+
+            // Sign of the result of the cross product tells the direction
+            const a: vec2 = [points[i][0] - points[current][0], points[i][1] - points[current][1]];
+            const b: vec2 = [points[next][0] - points[current][0], points[next][1] - points[current][1]];
+            const det = a[0] * b[1] - a[1] * b[0];
+            const dir = a[0] * b[0] + a[1] * b[1];
+
+            // Iterate over the points and select the rightmost point as 'next' relative to 'current'.
+            // If vectors a and b are collinear and point in the same direction, update 'next' to index i if i is
+            // further than next. If vectors a and b are collinear but point in opposite directions, index i is not
+            // considered rightmost but leftmost and is ignored. If next and current are overlapped, update 'next'
+            // to index i.
+            if (det > 0.0 || (det === 0.0 && dir >= 0.0 && (a[0] * a[0] + a[1] * a[1]) > (b[0] * b[0] + b[1] * b[1]))) {
+                next = i;
+            }
+        }
+
+        current = next;
+    } while (current !== leftmostIdx);
+
+    // Close the ring by adding the first point
+    if (ring.length > 0) {
+        ring.push(ring[0]);
+    }
+
+    return ring;
+}
+
 // In case of intersection, returns depth of the closest corner. Otherwise, returns undefined.
 export function queryGeometryIntersectsProjectedAabb(
     queryGeometry: Point[],
@@ -225,38 +299,17 @@ export function queryGeometryIntersectsProjectedAabb(
     const corners = Aabb.projectAabbCorners(aabb, worldViewProjection);
     // convert to screen points
     let minDepth = Number.MAX_VALUE;
-    let closestCornerIndex = -1;
     for (let c = 0; c < corners.length; ++c) {
         const corner = corners[c];
         corner[0] = (0.5 * corner[0] + 0.5) * transform.width;
         corner[1] = (0.5 - 0.5 * corner[1]) * transform.height;
         if (corner[2] < minDepth) {
-            closestCornerIndex = c;
             minDepth = corner[2]; // This is a rough aabb intersection check for now and no need to interpolate over aabb sides.
         }
     }
-    const p = (i: number): Point => new Point(corners[i][0], corners[i][1]);
-
-    let convexPolygon;
-    switch (closestCornerIndex) {
-    case 0:
-    case 6:
-        convexPolygon = [p(1), p(5), p(4), p(7), p(3), p(2), p(1)];
-        break;
-    case 1:
-    case 7:
-        convexPolygon = [p(0), p(4), p(5), p(6), p(2), p(3), p(0)];
-        break;
-    case 3:
-    case 5:
-        convexPolygon = [p(1), p(0), p(4), p(7), p(6), p(2), p(1)];
-        break;
-    default:
-        convexPolygon = [p(1), p(5), p(6), p(7), p(3), p(0), p(1)];
-        break;
-    }
 
-    // eslint-disable-next-line @typescript-eslint/no-unsafe-argument
+    const convexPolygon = convexPolygonFromPoints(corners);
+
     if (polygonIntersectsPolygon(queryGeometry, convexPolygon)) {
         return minDepth;
     }

test/integration/query-tests/model/sf-trees-puck-landmarks-meshopt-terrain/expected.json

@@ -18,7 +18,7 @@
   {
     "geometry": {
       "coordinates": [
-          -122.39465564489365,
+        -122.39465564489365,
         37.79525622982453
       ],
       "type": "Point"
@@ -37,7 +37,24 @@
   {
     "geometry": {
       "coordinates": [
-        -122.39468783140182,
+        -122.39474147558212,
+        37.79522231786474
+      ],
+      "type": "Point"
+    },
+    "id": 5832589281457790,
+    "properties": {
+      "tree_type": "palm"
+    },
+    "source": "trees",
+    "sourceLayer": "trees",
+    "state": {},
+    "type": "Feature"
+  },
+  {
+    "geometry": {
+      "coordinates": [
+          -122.39468783140182,
         37.79514813539846
       ],
       "type": "Point"

test/integration/query-tests/model/sf-trees-puck-landmarks-terrain/expected.json

@@ -18,7 +18,7 @@
   {
     "geometry": {
       "coordinates": [
-          -122.39465564489365,
+        -122.39465564489365,
         37.79525622982453
       ],
       "type": "Point"
@@ -37,7 +37,24 @@
   {
     "geometry": {
       "coordinates": [
-        -122.39468783140182,
+        -122.39474147558212,
+        37.79522231786474
+      ],
+      "type": "Point"
+    },
+    "id": 5832589281457790,
+    "properties": {
+      "tree_type": "palm"
+    },
+    "source": "trees",
+    "sourceLayer": "trees",
+    "state": {},
+    "type": "Feature"
+  },
+  {
+    "geometry": {
+      "coordinates": [
+          -122.39468783140182,
         37.79514813539846
       ],
       "type": "Point"

test/integration/query-tests/model/vector-source-query-geometry-intersection/expected.json

@@ -0,0 +1,17 @@
+[
+  {
+    "geometry": {
+      "coordinates": [
+        -74.49999991804361,
+        39.99999996437472
+      ],
+      "type": "Point"
+    },
+    "properties": {
+      "id": "object-one"
+    },
+    "source": "layer-source",
+    "state": {},
+    "type": "Feature"
+  }
+]

test/integration/query-tests/model/vector-source-query-geometry-intersection/style.json

@@ -0,0 +1,72 @@
+{
+  "version": 8,
+  "metadata": {
+    "test": {
+      "height": 512,
+      "width": 512,
+      "operations": [
+        [
+          "wait"
+        ],
+        [
+          "wait"
+        ],
+        [
+          "wait"
+        ]
+      ],
+      "queryGeometry": [
+        165,
+        279
+      ]
+    }
+  },
+  "sources": {
+    "layer-source": {
+      "type": "geojson",
+      "data": {
+        "type": "Feature",
+        "properties": {
+          "id": "object-one"
+        },
+        "geometry": {
+          "coordinates": [
+            -74.5,
+            40
+          ],
+          "type": "Point"
+        }
+      }
+    }
+  },
+  "transition": {
+    "duration": 0
+  },
+  "models": {
+    "model-glb": "local://models/Box.glb"
+  },
+  "pitch": 50,
+  "bearing": 45,
+  "zoom": 22,
+  "center": [
+    -74.5,
+    40
+  ],
+  "layers": [
+    {
+      "id": "model-layer",
+      "source": "layer-source",
+      "type": "model",
+      "layout": {
+        "model-id": "model-glb"
+      },
+      "paint": {
+        "model-scale": [
+          2,
+          2,
+          2
+        ]
+      }
+    }
+  ]
+}