From 53293e77a423f49c03cade03a5b576a8af55c264 Mon Sep 17 00:00:00 2001 From: Kevin-Mattheus-Moerman Date: Fri, 7 Mar 2025 10:37:19 +0000 Subject: [PATCH 01/11] Changed orthogonal_vector to use shoelace formula --- src/geometry_primitives.jl | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/geometry_primitives.jl b/src/geometry_primitives.jl index 7f60e481..a63f85dc 100644 --- a/src/geometry_primitives.jl +++ b/src/geometry_primitives.jl @@ -140,13 +140,23 @@ function collect_with_eltype!(result::AbstractVector{T}, iter) where {T} end """ - orthogonal_vector(p1, p2, p3) + orthogonal_vector(vertices) -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 face or polygon defined by a set of +ordered points contained in the point vector `vertices`. """ -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(vertices) + N = length(vertices) + c = zeros(eltype(vertices)) # Inherit vector type from input + @inbounds for i in eachindex(vertices) # Use shoelace approach + c += cross(vertices[i],vertices[mod1(i+1,N)]) # Add each edge contribution + end + return c +end + +function orthogonal_vector(VT,vertices) + return orthogonal_vector(VT.(vertices)) +end """ normals(positions::Vector{Point3{T}}, faces::Vector{<: NgonFace}[; normaltype = Vec3{T}]) @@ -165,12 +175,12 @@ end function normals(vertices::AbstractVector{<:Point{3}}, faces::AbstractVector{<: NgonFace}, ::Type{NormalType}) where {NormalType} - + println("WOOOHOOOO") 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 @@ -206,7 +216,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 From 0de3ee36e8ce3dd451bcdda1c44978c6e99d0903 Mon Sep 17 00:00:00 2001 From: Kevin-Mattheus-Moerman Date: Fri, 7 Mar 2025 15:29:17 +0000 Subject: [PATCH 02/11] Added testing for ngons and non-planar faces --- src/geometry_primitives.jl | 2 +- test/geometrytypes.jl | 23 +++++++++++++++++++++-- 2 files changed, 22 insertions(+), 3 deletions(-) diff --git a/src/geometry_primitives.jl b/src/geometry_primitives.jl index a63f85dc..d395f05c 100644 --- a/src/geometry_primitives.jl +++ b/src/geometry_primitives.jl @@ -175,7 +175,7 @@ end function normals(vertices::AbstractVector{<:Point{3}}, faces::AbstractVector{<: NgonFace}, ::Type{NormalType}) where {NormalType} - println("WOOOHOOOO") + normals_result = zeros(NormalType, length(vertices)) for face in faces v = vertices[face] diff --git a/test/geometrytypes.jl b/test/geometrytypes.jl index 0c06beb2..808d0226 100644 --- a/test/geometrytypes.jl +++ b/test/geometrytypes.jl @@ -351,12 +351,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 = [Point{3,Float64}(0.0,0.0,0.0),Point{3,Float64}(1.0,0.0,0.0),Point{3,Float64}(2.0,0.0,0.0),Point{3,Float64}(2.0,1.0,0.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 = [Point{3,Float64}(0.0,0.0,0.0),Point{3,Float64}(1.0,0.0,0.0),Point{3,Float64}(2.0,0.0,0.0),Point{3,Float64}(2.0,1.0,0.0),Point{3,Float64}(2.0,0.5,0.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,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 From 691f7342547b0bc30f85f67903285451f981e827 Mon Sep 17 00:00:00 2001 From: Kevin-Mattheus-Moerman Date: Fri, 7 Mar 2025 17:46:03 +0000 Subject: [PATCH 03/11] Further attempt to fix use of new orthogonal_vector --- src/geometry_primitives.jl | 21 +++++++++++-------- src/meshes.jl | 5 ++--- src/triangulation.jl | 43 ++++++++++++++------------------------ 3 files changed, 30 insertions(+), 39 deletions(-) diff --git a/src/geometry_primitives.jl b/src/geometry_primitives.jl index d395f05c..1a3e766a 100644 --- a/src/geometry_primitives.jl +++ b/src/geometry_primitives.jl @@ -143,19 +143,22 @@ end orthogonal_vector(vertices) Calculates an orthogonal vector to a face or polygon defined by a set of -ordered points contained in the point vector `vertices`. +ordered points contained in the point vector `vertices`. Note that for 2D points +a scalar is returned (the dot product of the orthogonal vector with the +positive Z-direction). """ -function orthogonal_vector(vertices) - N = length(vertices) - c = zeros(eltype(vertices)) # Inherit vector type from input +function orthogonal_vector(::Type{VT},vertices) where VT + c = zeros(VT) # Inherit vector type from input + j = length(vertices) @inbounds for i in eachindex(vertices) # Use shoelace approach - c += cross(vertices[i],vertices[mod1(i+1,N)]) # Add each edge contribution + c += cross(vertices[j],vertices[i]) # Add each edge contribution + j = i end return c end -function orthogonal_vector(VT,vertices) - return orthogonal_vector(VT.(vertices)) +function orthogonal_vector(vertices) + return orthogonal_vector(eltype(vertices),vertices) end """ @@ -178,7 +181,7 @@ function normals(vertices::AbstractVector{<:Point{3}}, faces::AbstractVector{<: normals_result = zeros(NormalType, length(vertices)) for face in faces - v = vertices[face] + v = coordinates(vertices[face]) # we can get away with two edges since faces are planar. n = orthogonal_vector(NormalType, v) for i in 1:length(face) @@ -215,7 +218,7 @@ end faces = resize!(F[], length(fs)) for (i, f) in enumerate(fs) - ps = positions[f] + ps = coordinates(positions[f]) n = orthogonal_vector(NormalType, ps) normals[i] = normalize(n) faces[i] = $(FT)(i) diff --git a/src/meshes.jl b/src/meshes.jl index 28a8d005..eb9a1358 100644 --- a/src/meshes.jl +++ b/src/meshes.jl @@ -213,9 +213,8 @@ end Calculate the signed volume of one tetrahedron. Be sure the orientation of your surface is right. """ -function volume(triangle::Triangle) - v1, v2, v3 = triangle - sig = sign(orthogonal_vector(v1, v2, v3) ⋅ v1) +function volume(triangle::Triangle) + sig = sign(orthogonal_vector(triangle.points) ⋅ v1) return sig * abs(v1 ⋅ (v2 × v3)) / 6 end diff --git a/src/triangulation.jl b/src/triangulation.jl index b0dbc3cf..f9a280db 100644 --- a/src/triangulation.jl +++ b/src/triangulation.jl @@ -13,54 +13,43 @@ 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) 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} + return T(0.5) * norm(orthogonal_vector(vertices)) end + """ in(point, triangle) From 202d37caecb2610c7b140cd178a61ab171f61db8 Mon Sep 17 00:00:00 2001 From: ffreyer Date: Fri, 7 Mar 2025 21:03:01 +0100 Subject: [PATCH 04/11] add to_ndim from Makie --- src/fixed_arrays.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/fixed_arrays.jl b/src/fixed_arrays.jl index 9e32e77c..1f4a9b6e 100644 --- a/src/fixed_arrays.jl +++ b/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 From 2448713d2e2a9313e1b6f1e8790dc30ef9fc9e47 Mon Sep 17 00:00:00 2001 From: ffreyer Date: Fri, 7 Mar 2025 21:04:10 +0100 Subject: [PATCH 05/11] bring back fast path for triangles --- src/geometry_primitives.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/geometry_primitives.jl b/src/geometry_primitives.jl index 1a3e766a..bf64f8e1 100644 --- a/src/geometry_primitives.jl +++ b/src/geometry_primitives.jl @@ -157,9 +157,13 @@ function orthogonal_vector(::Type{VT},vertices) where VT return c end -function orthogonal_vector(vertices) - return orthogonal_vector(eltype(vertices),vertices) +# 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 +end +orthogonal_vector(vertices) = orthogonal_vector(Vec3f, vertices) """ normals(positions::Vector{Point3{T}}, faces::Vector{<: NgonFace}[; normaltype = Vec3{T}]) From 992314ef27ff36b30002040de876fbc3b3487a49 Mon Sep 17 00:00:00 2001 From: ffreyer Date: Fri, 7 Mar 2025 21:04:46 +0100 Subject: [PATCH 06/11] try to derive eltypes --- src/geometry_primitives.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/geometry_primitives.jl b/src/geometry_primitives.jl index bf64f8e1..de41a1ee 100644 --- a/src/geometry_primitives.jl +++ b/src/geometry_primitives.jl @@ -162,7 +162,13 @@ function orthogonal_vector(::Type{VT}, triangle::Triangle) where {VT <: VecTypes 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::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) """ From e4b25ed65727d601474531a503bfcd79b7098ba8 Mon Sep 17 00:00:00 2001 From: ffreyer Date: Fri, 7 Mar 2025 21:05:43 +0100 Subject: [PATCH 07/11] move coordinates into orthogonal_vector and always return a vector --- src/geometry_primitives.jl | 40 +++++++++++++++++++------------------- src/meshes.jl | 5 +++-- src/triangulation.jl | 11 +++++------ 3 files changed, 28 insertions(+), 28 deletions(-) diff --git a/src/geometry_primitives.jl b/src/geometry_primitives.jl index de41a1ee..8969baa1 100644 --- a/src/geometry_primitives.jl +++ b/src/geometry_primitives.jl @@ -140,19 +140,19 @@ function collect_with_eltype!(result::AbstractVector{T}, iter) where {T} end """ - orthogonal_vector(vertices) + orthogonal_vector([target_type = Vec3f], points) -Calculates an orthogonal vector to a face or polygon defined by a set of -ordered points contained in the point vector `vertices`. Note that for 2D points -a scalar is returned (the dot product of the orthogonal vector with the -positive Z-direction). +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. """ -function orthogonal_vector(::Type{VT},vertices) where VT - c = zeros(VT) # Inherit vector type from input - j = length(vertices) - @inbounds for i in eachindex(vertices) # Use shoelace approach - c += cross(vertices[j],vertices[i]) # Add each edge contribution - j = i +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 @@ -177,8 +177,8 @@ orthogonal_vector(vertices) = orthogonal_vector(Vec3f, vertices) 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}; @@ -188,10 +188,10 @@ 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 = coordinates(vertices[face]) + v = vertices[face] # we can get away with two edges since faces are planar. n = orthogonal_vector(NormalType, v) for i in 1:length(face) @@ -211,15 +211,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))) @@ -228,7 +228,7 @@ end faces = resize!(F[], length(fs)) for (i, f) in enumerate(fs) - ps = coordinates(positions[f]) + ps = positions[f] n = orthogonal_vector(NormalType, ps) normals[i] = normalize(n) faces[i] = $(FT)(i) diff --git a/src/meshes.jl b/src/meshes.jl index eb9a1358..8a106971 100644 --- a/src/meshes.jl +++ b/src/meshes.jl @@ -213,8 +213,9 @@ end Calculate the signed volume of one tetrahedron. Be sure the orientation of your surface is right. """ -function volume(triangle::Triangle) - sig = sign(orthogonal_vector(triangle.points) ⋅ v1) +function volume(triangle::Triangle{3, T}) where {T} + v1, v2, v3 = triangle + sig = sign(orthogonal_vector(Vec3{T}, triangle) ⋅ v1) return sig * abs(v1 ⋅ (v2 × v3)) / 6 end diff --git a/src/triangulation.jl b/src/triangulation.jl index f9a280db..9d085d0e 100644 --- a/src/triangulation.jl +++ b/src/triangulation.jl @@ -17,8 +17,7 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI Calculate the area of one face. """ -function area(vertices::AbstractVector{<:Point}, - face::NgonFace) +function area(vertices::AbstractVector{<:Point}, face::NgonFace) return 0.5 * norm(orthogonal_vector(vertices[face])) end @@ -27,8 +26,7 @@ end Calculate the area of all faces. """ -function area(vertices::AbstractVector{<:Point}, - faces::AbstractVector{<:NgonFace}) +function area(vertices::AbstractVector{<:Point}, faces::AbstractVector{<:NgonFace}) return sum(x -> area(vertices, x), faces) end @@ -41,11 +39,12 @@ For 2D points, the oriented area is returned (negative when the points are oriented clockwise). """ function area(vertices::AbstractVector{Point{2,T}}) where {T} - length(vertices) < 3 && return zero(T) - return T(0.5) * orthogonal_vector(vertices) + length(vertices) < 3 && return zero(T) + return T(0.5) * orthogonal_vector(vertices)[3] end function area(vertices::AbstractVector{Point{3,T}}) where {T} + length(vertices) < 3 && return zero(T) return T(0.5) * norm(orthogonal_vector(vertices)) end From a18eef438bdc0ef64efe2b0f11d9335ababc7095 Mon Sep 17 00:00:00 2001 From: ffreyer Date: Fri, 7 Mar 2025 21:11:40 +0100 Subject: [PATCH 08/11] add some tests for to_ndim --- test/fixed_arrays.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/fixed_arrays.jl b/test/fixed_arrays.jl index 694811b1..cfb00fc0 100644 --- a/test/fixed_arrays.jl +++ b/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 From 34160dcd0e8f3323726b9f7f23cd71dbd773fb88 Mon Sep 17 00:00:00 2001 From: ffreyer Date: Fri, 7 Mar 2025 21:26:51 +0100 Subject: [PATCH 09/11] add basic tests for orthogonal_vector paths --- src/geometry_primitives.jl | 3 +++ test/geometrytypes.jl | 20 +++++++++++++++----- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/src/geometry_primitives.jl b/src/geometry_primitives.jl index 8969baa1..7c5b3035 100644 --- a/src/geometry_primitives.jl +++ b/src/geometry_primitives.jl @@ -165,6 +165,9 @@ 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 diff --git a/test/geometrytypes.jl b/test/geometrytypes.jl index 808d0226..f1f0b6cf 100644 --- a/test/geometrytypes.jl +++ b/test/geometrytypes.jl @@ -332,6 +332,16 @@ 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) + + 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) +end + @testset "Normals" begin # per face normals r = Rect3f(Point3f(0), Vec3f(1)) @@ -358,19 +368,19 @@ end @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 = [Point{3,Float64}(0.0,0.0,0.0),Point{3,Float64}(1.0,0.0,0.0),Point{3,Float64}(2.0,0.0,0.0),Point{3,Float64}(2.0,1.0,0.0)] + # 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 = [Point{3,Float64}(0.0,0.0,0.0),Point{3,Float64}(1.0,0.0,0.0),Point{3,Float64}(2.0,0.0,0.0),Point{3,Float64}(2.0,1.0,0.0),Point{3,Float64}(2.0,0.5,0.0)] + # 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 + # 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,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)] From e073fa8b441d86855de16521279ee206b6fcaa8a Mon Sep 17 00:00:00 2001 From: ffreyer Date: Fri, 7 Mar 2025 21:32:57 +0100 Subject: [PATCH 10/11] fix 1.6 --- test/geometrytypes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/geometrytypes.jl b/test/geometrytypes.jl index f1f0b6cf..44a0f508 100644 --- a/test/geometrytypes.jl +++ b/test/geometrytypes.jl @@ -381,7 +381,7 @@ end @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,6) + 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) From 1ec04f932d67f06171ac0a924e70696b6dd205af Mon Sep 17 00:00:00 2001 From: ffreyer Date: Fri, 21 Mar 2025 13:42:35 +0100 Subject: [PATCH 11/11] add tests for missed functions --- src/geometry_primitives.jl | 13 +++++++++++++ test/geometrytypes.jl | 24 ++++++++++++++++++++---- 2 files changed, 33 insertions(+), 4 deletions(-) diff --git a/src/geometry_primitives.jl b/src/geometry_primitives.jl index 7c5b3035..c70c796e 100644 --- a/src/geometry_primitives.jl +++ b/src/geometry_primitives.jl @@ -145,6 +145,8 @@ end 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. """ function orthogonal_vector(::Type{VT}, vertices) where {VT <: VecTypes{3}} c = zeros(VT) # Inherit vector type from input @@ -157,6 +159,17 @@ function orthogonal_vector(::Type{VT}, vertices) where {VT <: VecTypes{3}} 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 diff --git a/test/geometrytypes.jl b/test/geometrytypes.jl index 44a0f508..f35236ca 100644 --- a/test/geometrytypes.jl +++ b/test/geometrytypes.jl @@ -334,12 +334,28 @@ 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(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(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