diff --git a/Project.toml b/Project.toml index 14b2ec39..68bece68 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeometryBasics" uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326" authors = ["SimonDanisch "] -version = "0.5.5" +version = "0.5.6" [deps] EarCut_jll = "5ae413db-bbd1-5e63-b57d-d24a61df00f5" diff --git a/src/basic_types.jl b/src/basic_types.jl index 996b0800..48e776d8 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -31,14 +31,6 @@ abstract type AbstractNgonFace{N,T} <: AbstractFace{N,T} end abstract type AbstractSimplex{Dim,T} <: Polytope{Dim,T} end -@propagate_inbounds function Base.getindex(points::AbstractVector{P}, face::F) where {P<: Point, F <: AbstractFace} - return Polytope(P, F)(map(i-> points[i], face.data)) -end - -@propagate_inbounds function Base.getindex(elements::AbstractVector, face::F) where {F <: AbstractFace} - return map(i-> elements[i], face.data) -end - @fixed_vector SimplexFace = AbstractSimplexFace const TetrahedronFace{T} = SimplexFace{4,T} @@ -597,7 +589,12 @@ Mutating these will change the mesh. """ vertex_attributes(mesh::Mesh) = getfield(mesh, :vertex_attributes) -Base.getindex(mesh::Mesh, i::Integer) = mesh.position[mesh.faces[i]] +function Base.getindex(mesh::Mesh, i::Integer) + f = mesh.faces[i] + P = Polytope(eltype(mesh.position), typeof(f)) + return P(map(j -> mesh.position[j], f.data)) +end + Base.length(mesh::Mesh) = length(mesh.faces) function Base.:(==)(a::Mesh, b::Mesh) 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 diff --git a/src/geometry_primitives.jl b/src/geometry_primitives.jl index 7f60e481..c70c796e 100644 --- a/src/geometry_primitives.jl +++ b/src/geometry_primitives.jl @@ -140,13 +140,52 @@ 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}]) @@ -154,8 +193,8 @@ orthogonal_vector(::Type{VT}, p1, p2, p3) where {VT} = orthogonal_vector(VT(p1), 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 diff --git a/src/meshes.jl b/src/meshes.jl index 28a8d005..cb38a721 100644 --- a/src/meshes.jl +++ b/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 @@ -390,6 +390,8 @@ function expand_faceviews(mesh::Mesh) end end +expand_faceviews(m::MetaMesh) = MetaMesh(expand_faceviews(Mesh(m)), meta(m)) + function merge_vertex_indices( faces::NTuple{N_Attrib, <: AbstractVector{FT}}, vertex_index_counter::Integer = T(1) diff --git a/src/triangulation.jl b/src/triangulation.jl index b0dbc3cf..9d085d0e 100644 --- a/src/triangulation.jl +++ b/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) 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 diff --git a/test/geometrytypes.jl b/test/geometrytypes.jl index 0c06beb2..51746602 100644 --- a/test/geometrytypes.jl +++ b/test/geometrytypes.jl @@ -298,10 +298,15 @@ end ps = rand(Point2f, 10) f = GLTriangleFace(1, 2, 3) - @test ps[f] == Triangle(ps[[1,2,3]]...) + @test ps[f] == GeometryBasics.@SVector [ps[1], ps[2], ps[3]] + data = [string(i) for i in 1:10] f = QuadFace(3, 4, 7, 8) - @test data[f] == ("3", "4", "7", "8") + @test data[f] == ["3", "4", "7", "8"] + + @test (f .== 4) == QuadFace(false, true, false, false) + @test f[f .== 4] isa Vector{Int} + @test f[f .== 4] == [4] @test GeometryBasics.cyclic_hash(f) != GeometryBasics.cyclic_hash(QuadFace(1,2,3,4)) @test GeometryBasics.cyclic_hash(f) == GeometryBasics.cyclic_hash(QuadFace(3,4,7,8)) @@ -332,6 +337,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 +382,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 diff --git a/test/meshes.jl b/test/meshes.jl index e6c0b0e8..92d396ad 100644 --- a/test/meshes.jl +++ b/test/meshes.jl @@ -54,6 +54,11 @@ end @test normals(m2) != normals(m) @test normals(m2) == [only(values(normals(m))) for _ in 1:4] @test isempty(m2.views) + + mm = MetaMesh(m, Dict(:test => 1, :a => "a")) + mm2 = GeometryBasics.expand_faceviews(mm) + @test mm2.meta == mm.meta + @test Mesh(mm2) == m2 end @testset "Duplicate face removal" begin