Check coplanarity using all points and distance tolerance

AUTHORS.md

@@ -12,6 +12,7 @@
 
 ## Contributors
 
+- Anton Soloviev <<anton@praviel.com>> [@antonsoo](https://github.com/antonsoo)
 - Andrew Liew <<a.liew@sheffield.ac.uk>> [@andrewliew](https://github.com/andrewliew)
 - Tomás Méndez Echenagucia <<mendez@arch.ethz.ch>> [@tmsmendez](https://github.com/tmsmendez)
 - Matthias Rippmann <<rippmann@arch.ethz.ch>> [@Rippmann](https://github.com/Rippmann)

CHANGELOG.md

@@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 
+* Fixed `is_coplanar` and `Polygon.is_planar` discarding points, using a triangle-area-dependent tolerance and rejecting noisy coplanar clouds.
 * Changed `Data.__dtype__` to delegate to `Data.__clstype__` (removes duplicated dtype format logic).
 * Changed `pca_numpy` to run the SVD on the centered data matrix instead of the covariance matrix: forming the covariance squares the condition number, so near-degenerate inputs (e.g. almost collinear point clouds) returned wrong principal directions — `bestfit_plane_numpy` normals were off by 11–37 degrees on exactly planar sliver clouds (#1522). Well-conditioned results are unchanged (same eigenvectors; eigenvalues rescaled to keep their variance meaning).
 * Changed `Tolerance` class to no longer use singleton pattern. `Tolerance()` now creates independent instances instead of returning the global `TOL`. 

src/compas/geometry/_core/predicates_3.py

@@ -13,6 +13,7 @@
 from compas.geometry import length_vector
 from compas.geometry import normal_polygon
 from compas.geometry import subtract_vectors
+from compas.geometry.bestfit import bestfit_plane
 from compas.itertools import window
 from compas.tolerance import TOL
 
@@ -101,43 +102,62 @@ def is_coplanar(points, tol=None):
     points : sequence[point]
         A sequence of point locations.
     tol : float, optional
-        A tolerance for planarity validation.
+        Absolute distance tolerance for planarity validation.
+        Defaults to the global absolute tolerance.
 
     Returns
     -------
     bool
         True if the points are coplanar.
         False otherwise.
 
-    See Also
-    --------
-    is_colinear
-    is_colinear_line_line
-
     Notes
     -----
-    Compute the normal vector (cross product) of the vectors formed by the first three points.
-    Taking the first point as base point, include one more vector at a time and check if that vector is perpendicular to the normal vector.
-    If all vectors are perpendicular, the points are coplanar.
+    Construct a stable frame from the farthest baseline endpoint and the point
+    farthest from that line. Points within the tolerance of the baseline are
+    considered collinear. Otherwise, center the plane between the extreme
+    projections of all points. If this plane exceeds the tolerance, fit a plane
+    to all points in the local frame and repeat the centered distance check.
+    The tolerance applies to every point, including the frame's basis points.
 
     """
     if len(points) < 4:
         return True
 
-    temp = points[:]
+    a = points[0]
+    b_index = max(range(1, len(points)), key=lambda i: distance_point_point(a, points[i]))
+    b = points[b_index]
+    if TOL.is_zero(distance_point_point(a, b), tol):
+        return True
 
-    while len(temp) >= 3:
-        a = temp.pop(0)
-        b = temp.pop(0)
-        c = temp.pop(0)
-        if not is_colinear(a, b, c, tol):
-            break
-        if not temp:
-            return True
+    c_index = max(range(1, len(points)), key=lambda i: distance_point_line(points[i], (a, b)))
+    c = points[c_index]
+    if TOL.is_zero(distance_point_line(c, (a, b)), tol):
+        return True
 
     n = cross_vectors(subtract_vectors(b, a), subtract_vectors(c, a))
+    length = length_vector(n)
+
+    distances = [dot_vectors(n, subtract_vectors(point, a)) / length for point in points]
+    if TOL.is_zero(0.5 * (max(distances) - min(distances)), tol):
+        return True
 
-    return all(is_perpendicular_vector_vector(n, subtract_vectors(d, a), tol) for d in temp)
+    # Fit in a stable local frame to avoid cancellation in covariance minors of sliver clouds.
+    n = [value / length for value in n]
+    baseline = subtract_vectors(b, a)
+    baseline_length = length_vector(baseline)
+    u = [value / baseline_length for value in baseline]
+    v = cross_vectors(n, u)
+    local = []
+    for point, distance in zip(points, distances):
+        vector = subtract_vectors(point, a)
+        local.append([dot_vectors(vector, u), dot_vectors(vector, v), distance])
+    _, normal = bestfit_plane(local)
+    world_normal = [normal[0] * u[i] + normal[1] * v[i] + normal[2] * n[i] for i in range(3)]
+    projections = [dot_vectors(point, normal) for point in local]
+
+    # A centered offset permits residuals on the basis points as well as the other samples.
+    return TOL.is_zero(0.5 * (max(projections) - min(projections)) / length_vector(world_normal), tol)
 
 
 # =============================================================================

tests/compas/geometry/test_core_predicates_3.py

@@ -1,4 +1,5 @@
 from compas.geometry import is_coplanar
+from compas.geometry import Polygon
 import pytest
 
 
@@ -10,3 +11,107 @@ def test_is_coplanar():
         is_coplanar(points)
     except IndexError:
         pytest.fail("is_coplanar raised an IndexError")
+
+
+@pytest.mark.parametrize(
+    "points",
+    [
+        [[0, 0, 0], [1, 0, 0], [2, 0, 0], [0, 0, 1], [0, 1, 1], [1, 0, 1]],
+        [[0, 0, 0], [1, 0, 0], [2, 0, 0], [0, 1, 0], [1, 1, 0], [2, 1, 0], [0, 0, 1], [1, 0, 1], [2, 0, 1]],
+        [[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 1], [0, 1, 1], [1, 0, 1]],
+    ],
+)
+@pytest.mark.parametrize("shift", [0, 1, 2])
+def test_is_coplanar_retains_points_from_collinear_triples(points, shift):
+    points = points[shift:] + points[:shift]
+    original = [point[:] for point in points]
+
+    assert not is_coplanar(points)
+    assert not Polygon(points).is_planar
+    assert points == original
+
+
+@pytest.mark.parametrize("tol", [None, 1e-8])
+@pytest.mark.parametrize("reverse", [False, True])
+def test_is_coplanar_uses_well_conditioned_third_point(tol, reverse):
+    epsilon = 0.22 * (1e-9 if tol is None else tol)
+    remaining = [[0, 1, 0.4], [0, 0.5, -0.4]]
+    if reverse:
+        remaining.reverse()
+    points = [[0, 0, 0], [1, 0, 0], [2, 0, 0], [10, 0, 0], [0, epsilon, 0]] + remaining
+    original = [point[:] for point in points]
+
+    assert not is_coplanar(points, tol=tol)
+    assert not Polygon(points).is_planar
+    assert points == original
+
+
+def test_is_coplanar_normalizes_plane_distance():
+    points = [[0, 0, 0], [0.05, 0, 0], [0.1, 0, 0], [0, 1e-4, 0], [0, -1e-4, 0], [0.1, 0, 0], [0, 0, 9e-5]]
+
+    assert not is_coplanar(points)
+    assert not Polygon(points).is_planar
+
+
+@pytest.mark.parametrize("scale", [2.0**-20, 1.0, 2.0**20])
+@pytest.mark.parametrize("tol", [None, 1e-8])
+@pytest.mark.parametrize("height, expected", [(0.5, True), (2.0, False)])
+def test_is_coplanar_tolerance_is_a_distance(scale, tol, height, expected):
+    h = height * (1e-9 if tol is None else tol)
+    points = [
+        [-scale, -scale, 0],
+        [0, 0, 0],
+        [scale, scale, 0],
+        [-scale, scale, 0],
+        [scale, -scale, 0],
+        [0, 0, h],
+        [0, 0, -h],
+    ]
+    original = [point[:] for point in points]
+
+    assert is_coplanar(points, tol=tol) == expected
+    assert points == original
+
+
+@pytest.mark.parametrize("tol", [None, 1e-8])
+@pytest.mark.parametrize("axis", [0, 1, 2])
+@pytest.mark.parametrize("reverse", [False, True])
+def test_is_coplanar_fits_noisy_points(tol, axis, reverse):
+    h = 0.6 * (1e-9 if tol is None else tol)
+    points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [10, 10, h], [10, -10, -h], [-10, 10, -h], [-10, -10, h]]
+    if reverse:
+        points = points[:1] + points[:0:-1]
+    points = [[point[(i + axis) % 3] for i in range(3)] for point in points]
+    original = [point[:] for point in points]
+
+    assert is_coplanar(points, tol)
+    if tol is None:
+        assert Polygon(points).is_planar
+    assert points == original
+
+
+@pytest.mark.parametrize("tol", [None, 1e-8])
+@pytest.mark.parametrize("axis", [0, 1, 2])
+@pytest.mark.parametrize("height, expected", [(1.8, True), (2.2, False)])
+def test_is_coplanar_allows_residuals_for_every_point(tol, axis, height, expected):
+    h = height * (1e-9 if tol is None else tol)
+    points = [[x, y, z] for z in [0, h] for x, y in [(0, 0), (2, 0), (2, 2), (0, 2)]]
+    points = [[point[(i + axis) % 3] for i in range(3)] for point in points]
+
+    assert is_coplanar(points, tol) == expected
+
+
+@pytest.mark.parametrize(
+    "points",
+    [
+        [[0, 0, 0]] * 6,
+        [[x, 0, 0] for x in range(7)],
+        [[0, 0, 0], [1, 0, 0], [2, 0, 0], [2, 1, 0], [1, 1, 0], [0, 1, 0]],
+        [[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 1, 0], [1, 0, 0], [1, 1, 0]],
+    ],
+)
+@pytest.mark.parametrize("offset", [[0, 0, 0], [8, -4, 2]])
+def test_is_coplanar_accepts_planar_degenerate_prefixes(points, offset):
+    translated = [[value + delta for value, delta in zip(point, offset)] for point in points]
+
+    assert is_coplanar(translated)