Use shoelace approach to compute normals (#245)

* Changed orthogonal_vector to use shoelace formula

* Added testing for ngons and non-planar faces

* Further attempt to fix use of new orthogonal_vector

* add to_ndim from Makie

* bring back fast path for triangles

* try to derive eltypes

* move coordinates into orthogonal_vector and always return a vector

* add some tests for to_ndim

* add basic tests for orthogonal_vector paths

* fix 1.6

* add tests for missed functions

---------

Co-authored-by: ffreyer <[email protected]>

Author: Kevin-Mattheus-Moerman

src/fixed_arrays.jl

@@ -145,6 +145,13 @@ Base.isnan(p::Union{AbstractPoint,Vec}) = any(isnan, p)
 Base.isinf(p::Union{AbstractPoint,Vec}) = any(isinf, p)
 Base.isfinite(p::Union{AbstractPoint,Vec}) = all(isfinite, p)
 
+function to_ndim(T::Type{<: VecTypes{N, ET}}, vec::VecTypes{N2}, fillval) where {N,ET,N2}
+    T(ntuple(Val(N)) do i
+        i > N2 && return ET(fillval)
+        @inbounds return vec[i]
+    end)
+end
+
 ## Generate aliases
 ## As a text file instead of eval/macro, to not confuse code linter
 

src/geometry_primitives.jl

@@ -140,22 +140,61 @@ function collect_with_eltype!(result::AbstractVector{T}, iter) where {T}
 end
 
 """
-    orthogonal_vector(p1, p2, p3)
+    orthogonal_vector([target_type = Vec3f], points)
 
-Calculates an orthogonal vector `cross(p2 - p1, p3 - p1)` to a plane described
-by 3 points p1, p2, p3. 
+Calculates an orthogonal vector to a polygon defined by a vector of ordered
+`points`. Note that the orthogonal vector to a collection of 2D points needs to
+expand to 3D space.
+
+Note that this vector is not normalized.
 """
-orthogonal_vector(p1, p2, p3) = cross(p2 - p1, p3 - p1)
-orthogonal_vector(::Type{VT}, p1, p2, p3) where {VT} = orthogonal_vector(VT(p1), VT(p2), VT(p3))
+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
+        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}}
+    a, b, c = triangle
+    return cross(to_ndim(VT, b-a, 0), to_ndim(VT, c-a, 0))
+end
+
+# derive target type
+orthogonal_vector(vertices::Ngon{D, T}) where {D, T} = orthogonal_vector(Vec3{T}, vertices)
+function orthogonal_vector(vertices::NTuple{N, VT}) where {N, D, T, VT <: VecTypes{D, T}}
+    return orthogonal_vector(Vec3{T}, vertices)
+end
+function orthogonal_vector(vertices::AbstractArray{VT}) where {D, T, VT <: VecTypes{D, T}}
+    return orthogonal_vector(Vec3{T}, vertices)
+end
+# fallback to Vec3f if vertices is something else
+orthogonal_vector(vertices) = orthogonal_vector(Vec3f, vertices)
 
 """
     normals(positions::Vector{Point3{T}}, faces::Vector{<: NgonFace}[; normaltype = Vec3{T}])
 
 Compute vertex normals from the given `positions` and `faces`.
 
 This runs through all faces, computing a face normal each and adding it to every
-involved vertex. The direction of the face normal is based on winding direction 
-and assumed counter-clockwise faces. At the end the summed face normals are 
+involved vertex. The direction of the face normal is based on winding direction
+and assumed counter-clockwise faces. At the end the summed face normals are
 normalized again to produce a vertex normal.
 """
 function normals(vertices::AbstractVector{Point{3,T}}, faces::AbstractVector{F};
@@ -165,12 +204,12 @@ end
 
 function normals(vertices::AbstractVector{<:Point{3}}, faces::AbstractVector{<: NgonFace},
                  ::Type{NormalType}) where {NormalType}
-                 
+
     normals_result = zeros(NormalType, length(vertices))
     for face in faces
         v = vertices[face]
         # we can get away with two edges since faces are planar.
-        n = orthogonal_vector(NormalType, v[1], v[2], v[3])
+        n = orthogonal_vector(NormalType, v)
         for i in 1:length(face)
             fi = face[i]
             normals_result[fi] = normals_result[fi] .+ n
@@ -188,15 +227,15 @@ Compute face normals from the given `positions` and `faces` and returns an
 appropriate `FaceView`.
 """
 function face_normals(
-        positions::AbstractVector{<:Point3{T}}, fs::AbstractVector{<: AbstractFace}; 
+        positions::AbstractVector{<:Point3{T}}, fs::AbstractVector{<: AbstractFace};
         normaltype = Vec3{T}) where {T}
     return face_normals(positions, fs, normaltype)
 end
-    
+
 @generated function face_normals(positions::AbstractVector{<:Point3}, fs::AbstractVector{F},
         ::Type{NormalType}) where {F<:NgonFace,NormalType}
-    
-    # If the facetype is not concrete it likely varies and we need to query it 
+
+    # If the facetype is not concrete it likely varies and we need to query it
     # doing the iteration
     FT = ifelse(isconcretetype(F), :($F), :(typeof(f)))
 
@@ -206,7 +245,7 @@ end
 
         for (i, f) in enumerate(fs)
             ps = positions[f]
-            n = orthogonal_vector(NormalType, ps[1], ps[2], ps[3])
+            n = orthogonal_vector(NormalType, ps)
             normals[i] = normalize(n)
             faces[i] = $(FT)(i)
         end

src/meshes.jl

@@ -213,9 +213,9 @@ end
 Calculate the signed volume of one tetrahedron. Be sure the orientation of your
 surface is right.
 """
-function volume(triangle::Triangle)
+function volume(triangle::Triangle{3, T}) where {T}
     v1, v2, v3 = triangle
-    sig = sign(orthogonal_vector(v1, v2, v3) ⋅ v1)
+    sig = sign(orthogonal_vector(Vec3{T}, triangle) ⋅ v1)
     return sig * abs(v1 ⋅ (v2 × v3)) / 6
 end
 

src/triangulation.jl

@@ -13,54 +13,42 @@ The above copyright notice and this permission notice shall be included in all c
 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 =#
 """
-    area(vertices::AbstractVector{Point{3}}, face::TriangleFace)
+    area(vertices::AbstractVector{Point{3}}, face::NgonFace)
 
-Calculate the area of one triangle.
+Calculate the area of one face.
 """
-function area(vertices::AbstractVector{<:Point{3,VT}},
-              face::TriangleFace{FT}) where {VT,FT}
-    v1, v2, v3 = vertices[face]
-    return 0.5 * norm(orthogonal_vector(v1, v2, v3))
+function area(vertices::AbstractVector{<:Point}, face::NgonFace)
+    return 0.5 * norm(orthogonal_vector(vertices[face]))
 end
 
 """
-    area(vertices::AbstractVector{Point{3}}, faces::AbstractVector{TriangleFace})
+    area(vertices::AbstractVector{Point}, faces::AbstractVector{NgonFace})
 
-Calculate the area of all triangles.
+Calculate the area of all faces.
 """
-function area(vertices::AbstractVector{Point{3,VT}},
-              faces::AbstractVector{TriangleFace{FT}}) where {VT,FT}
+function area(vertices::AbstractVector{<:Point}, faces::AbstractVector{<:NgonFace})
     return sum(x -> area(vertices, x), faces)
 end
 
 """
-    area(contour::AbstractVector{Point}})
+    area(vertices::AbstractVector{Point}})
 
 Calculate the area of a polygon.
 
 For 2D points, the oriented area is returned (negative when the points are
 oriented clockwise).
 """
-function area(contour::AbstractVector{Point{2,T}}) where {T}
-    length(contour) < 3 && return zero(T)
-    A = zero(T)
-    p = lastindex(contour)
-    for q in eachindex(contour)
-        A += cross(contour[p], contour[q])
-        p = q
-    end
-    return A * T(0.5)
+function area(vertices::AbstractVector{Point{2,T}}) where {T}
+    length(vertices) < 3 && return zero(T)
+    return T(0.5) * orthogonal_vector(vertices)[3]
 end
 
-function area(contour::AbstractVector{Point{3,T}}) where {T}
-    A = zero(eltype(contour))
-    o = first(contour)
-    for i in (firstindex(contour) + 1):(lastindex(contour) - 1)
-        A += cross(contour[i] - o, contour[i + 1] - o)
-    end
-    return norm(A) * T(0.5)
+function area(vertices::AbstractVector{Point{3,T}}) where {T}
+    length(vertices) < 3 && return zero(T)
+    return T(0.5) * norm(orthogonal_vector(vertices))
 end
 
+
 """
     in(point, triangle)
 

test/fixed_arrays.jl

@@ -1,4 +1,5 @@
 using Test
+using GeometryBasics: to_ndim
 
 @testset "Construction and Conversion" begin
     for VT in [Point, Vec]
@@ -19,6 +20,10 @@ using Test
 
     @test convert(Point, (2, 3)) === Point(2, 3)
     @test convert(Point, (2.0, 3)) === Point(2.0, 3.0)
+
+    @test to_ndim(Point3f, Vec2i(1,2), 0) == Point3f(1,2,0)
+    @test to_ndim(Vec4i, (1f0, 2), 0) == Vec4i(1,2,0,0)
+    @test to_ndim(NTuple{2, Float64}, Point3f(1,2,3), 0) == (1.0, 2.0)
 end
 
 @testset "broadcast" begin

test/geometrytypes.jl

@@ -332,6 +332,32 @@ end
     @test length(fv) == 5
 end
 
+@testset "orthogonal_vector" begin
+    tri = Triangle(Point3d(0,0,0), Point3d(1,0,0), Point3d(0,1,0))
+    @test GeometryBasics.orthogonal_vector(tri) === Vec3d(0,0,1)
+    @test GeometryBasics.orthogonal_vector(collect(coordinates(tri))) === Vec3d(0,0,1)
+    @test GeometryBasics.orthogonal_vector(Vec3f, tri) === Vec3f(0,0,1)
+    @test GeometryBasics.orthogonal_vector(Vec3f, collect(coordinates(tri))) === Vec3f(0,0,1)
+
+    quad = GeometryBasics.Quadrilateral(Point2i(0,0), Point2i(1,0), Point2i(1,1), Point2i(0,1))
+    @test GeometryBasics.orthogonal_vector(quad) === Vec3i(0,0,2)
+    @test GeometryBasics.orthogonal_vector(collect(coordinates(quad))) === Vec3i(0,0,2)
+    @test GeometryBasics.orthogonal_vector(Vec3d, quad) === Vec3d(0,0,2)
+    @test GeometryBasics.orthogonal_vector(Vec3d, collect(coordinates(quad))) === Vec3d(0,0,2)
+
+    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)#
+
+    # Maybe the ::Any fallback is too generic...?
+    struct TestType
+        data::Vector{Vec3f}
+    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)
+end
+
 @testset "Normals" begin
     # per face normals
     r = Rect3f(Point3f(0), Vec3f(1))
@@ -351,12 +377,31 @@ end
     ns = normals(c)
     # caps without mantle
     f_ns = face_normals(coordinates(c), filter!(f -> f isa TriangleFace, faces(c)))
-    @test all(n -> n == values(ns)[end-1], values(f_ns)[1:15])
-    @test all(n -> n == values(ns)[end], values(f_ns)[16:end])
+    @test all(n -> n ≈ values(ns)[end-1], values(f_ns)[1:15])
+    @test all(n -> n ≈ values(ns)[end], values(f_ns)[16:end])
     # Mantle without caps
     v_ns = normals(coordinates(c), filter!(f -> f isa QuadFace, faces(c)))[1:end-2]
     @test values(ns)[1:15] ≈ v_ns[1:15]
     @test values(ns)[1:15] ≈ v_ns[16:30] # repeated via FaceView in ns
+
+    # Planar QuadFace with colinear edge
+    v = [Point3d(0,0,0),Point3d(1,0,0),Point3d(2,0,0),Point3d(2,1,0)]
+    f = [QuadFace{Int}(1,2,3,4)]
+    n = normals(v,f)
+    @test all(n_i -> n_i ≈ [0.0,0.0,1.0], n)
+
+    # Planar NgonFace (5-sided) with colinear edge
+    v = [Point3d(0,0,0),Point3d(1,0,0),Point3d(2,0,0),Point3d(2,1,0),Point3d(2,0.5,0)]
+    f = [NgonFace{5,Int}(1,2,3,4,5)]
+    n = normals(v,f)
+    @test all(n_i -> n_i ≈ [0.0,0.0,1.0], n)
+
+    # Non-planar NgonFace (6 sided), features equal up and down variations resulting in z-dir average face_normal
+    t = range(0.0, 2*pi-(2*pi)/6, length = 6)
+    v = [Point{3,Float64}(cos(t[i]),sin(t[i]),iseven(i)) for i in eachindex(t)]
+    f = [NgonFace{6,Int}(1,2,3,4,5,6)]
+    n = normals(v,f)
+    @test all(n_i -> n_i ≈ [0.0,0.0,1.0], n)
 end
 
 @testset "HyperSphere" begin