Fix float stability in orthogonal_vector (#252)

* fix float precision issue?

* make fast path more reachable

* add sanity check for larger polygon face

* test stability

* fix stability

* add back old Syntax for Makie compat

* add comments

Author: ffreyer

src/geometry_primitives.jl

@@ -148,30 +148,32 @@ expand to 3D space.
 
 Note that this vector is not normalized.
 """
-function orthogonal_vector(::Type{VT}, vertices) where {VT <: VecTypes{3}}
-    c = zeros(VT) # Inherit vector type from input
-    prev = to_ndim(VT, last(coordinates(vertices)), 0)
-    @inbounds for p in coordinates(vertices) # Use shoelace approach
-        v = to_ndim(VT, p, 0)
-        c += cross(prev, v) # Add each edge contribution
-        prev = v
-    end
-    return c
-end
-
-function orthogonal_vector(::Type{VT}, vertices::Tuple) where {VT <: VecTypes{3}}
-    c = zeros(VT) # Inherit vector type from input
-    prev = to_ndim(VT, last(vertices), 0)
-    @inbounds for p in vertices # Use shoelace approach
-        v = to_ndim(VT, p, 0)
-        c += cross(prev, v) # Add each edge contribution
+orthogonal_vector(::Type{VT}, vertices) where {VT <: VecTypes{3}} = _orthogonal_vector(VT, coordinates(vertices))
+orthogonal_vector(::Type{VT}, vertices::Tuple) where {VT <: VecTypes{3}} = _orthogonal_vector(VT, vertices)
+
+function _orthogonal_vector(::Type{VT}, vertices) where {VT <: VecTypes{3}}
+    # Using shoelace approach (with N+1 = 1) (see #245)
+    # \sum_i p_i × p_{i+1}
+    # with distance vectors to avoid float precision issues
+    # \sum_i (p_i - p_1) × (p_{i+1} - p_1)
+    # These terms are equal when expanded (extra terms cancel over full sum)
+    c = zero(VT)
+    p0 = first(vertices)
+    prev = zero(VT)
+    for i in eachindex(vertices)
+        # i = 1 and N don't contribute as then produce (q_1 - q_1) terms
+        # we need i = 1 to set up prev though
+        i == lastindex(vertices) && break
+        v = to_ndim(VT, vertices[i+1] - p0, 0)
+        c += cross(prev, v)
         prev = v
     end
     return c
 end
 
 # Not sure how useful this fast path is, but it's simple to keep
-function orthogonal_vector(::Type{VT}, triangle::Triangle) where {VT <: VecTypes{3}}
+orthogonal_vector(p1::VT, p2::VT, p3::VT) where {VT <: VecTypes{3}} = _orthogonal_vector(Vec3f, to_ndim.(Vec3f, (p1, p2, p3), 0))
+function orthogonal_vector(::Type{VT}, triangle::Union{Triangle, NTuple{3, <:VecTypes}, StaticVector{3, <:VecTypes}}) where {VT <: VecTypes{3}}
     a, b, c = triangle
     return cross(to_ndim(VT, b-a, 0), to_ndim(VT, c-a, 0))
 end

test/geometrytypes.jl

@@ -352,7 +352,14 @@ end
 
     t = (Point3f(0), Point3f(1,0,1), Point3f(0,1,0))
     @test GeometryBasics.orthogonal_vector(t) == Vec3f(-1,0,1)
-    @test GeometryBasics.orthogonal_vector(Vec3i, t) == Vec3i(-1,0,1)#
+    @test GeometryBasics.orthogonal_vector(Vec3i, t) == Vec3i(-1,0,1)
+
+    t = Point2f[(0,0), (1,0), (2,1), (1,2), (0,2), (-1, 1)]
+    @test GeometryBasics.orthogonal_vector(Vec3f, t) == Vec3f(0,0,8)
+
+    # Makie uses this syntax
+    @test GeometryBasics.orthogonal_vector(Point3f(0), Point3f(1,0,1), Point3f(0,1,0)) == Vec3f(-1,0,1)
+    @test GeometryBasics.orthogonal_vector((1,2,3), (2,3,1), (3,1,2)) == Vec3f(-3.0)
 
     # Maybe the ::Any fallback is too generic...?
     struct TestType
@@ -361,6 +368,32 @@ end
     GeometryBasics.coordinates(x::TestType) = x.data
     x = TestType([Point3f(1,1,1), Point3f(0,0,0), Point3f(0.5,0,0)])
     @test GeometryBasics.orthogonal_vector(x) == Vec3f(0, -0.5, 0.5)
+
+    @testset "Float Precision" begin
+        # Collection of risky triangles
+        tris = [
+            # with 0s that remove terms from cross(a-b, c-b)
+            (Point3f(-10.180012, 9.928894, 63.42825), Point3f(-10.180012, 9.92804, 63.424137), Point3f(-10.181978, 9.92804, 63.42825)),
+            (Point3f(-11.75692, 6.245064, 86.05223), Point3f(-11.758429, 6.245064, 86.048584), Point3f(-11.758429, 6.2428894, 86.05223)),
+            (Point3f(52.430557, 13.611008, 7.657501), Point3f(52.43055, 13.611008, 7.6574783), Point3f(52.43055, 13.611008, 7.657501)),
+            (Point3f(52.430557, 13.611008, 7.657501), Point3f(52.43055, 13.611008, 7.657501), Point3f(52.45581, 13.61525, 8.18364)),
+            # tiny shift
+            (Point3f(1, 1 + 1e-4, 1), Point3f(1 + 1e-4, 1 + 1e-4, 1 - 1e-4), Point3f(1 - 1e-4, 1, 1 + 1e-4))
+        ]
+
+        for tri in tris
+            a, b, c = tri
+            n64 = Vec3f(cross(Point3d(b) - Point3d(a), Point3d(c) - Point3d(a)))
+            # fast paths
+            n_fast = GeometryBasics.orthogonal_vector(Vec3f, Triangle(a, b, c))
+            @test n64 ≈ n_fast
+            @test n_fast == GeometryBasics.orthogonal_vector(Vec3f, (a, b, c))
+            @test n_fast == GeometryBasics.orthogonal_vector(Vec3f, GeometryBasics.SVector(a, b, c)) # hit by points[face]
+            n_slow = GeometryBasics.orthogonal_vector(Vec3f, TestType([a, b, c]))
+            @test n_slow ≈ n64
+            @test n_slow == GeometryBasics.orthogonal_vector(Vec3f, [a, b, c])
+        end
+    end
 end
 
 @testset "Normals" begin