Merge pull request #80 from JuliaGeometry/sd/polygons

make polygon plotting work

Author: SimonDanisch

Project.toml

@@ -4,6 +4,7 @@ authors = ["SimonDanisch <sdanisch@gmail.com>"]
 version = "0.2.15"
 
 [deps]
+EarCut_jll = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

src/GeometryBasics.jl

@@ -1,6 +1,7 @@
 module GeometryBasics
 
     using StaticArrays, Tables, StructArrays, IterTools, LinearAlgebra
+    using EarCut_jll
 
     using Base: @propagate_inbounds
 
@@ -12,15 +13,15 @@ module GeometryBasics
     include("viewtypes.jl")
     include("geometry_primitives.jl")
     include("rectangles.jl")
-    include("triangulation.jl")
     include("meshes.jl")
+    include("triangulation.jl")
     include("lines.jl")
     include("boundingboxes.jl")
 
     export AbstractGeometry, GeometryPrimitive
     export Mat, Point, Vec
     export LineFace, Polytope, Line, NgonFace, convert_simplex
-    export LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon
+    export LineString, AbstractPolygon, Polygon, MultiPoint, MultiLineString, MultiPolygon
     export Simplex, connect, Triangle, NSimplex, Tetrahedron
     export QuadFace, metafree, coordinates, TetrahedronFace
     export TupleView, SimplexFace, Mesh, meta
@@ -36,7 +37,7 @@ module GeometryBasics
     export AbstractMesh, Mesh, TriangleMesh
     export GLNormalMesh2D, PlainTriangleMesh
     export MetaT, meta_table
-    
+
     # all the different predefined mesh types
     # Note: meshes can contain arbitrary meta information,
     export AbstractMesh, TriangleMesh, PlainMesh, GLPlainMesh, GLPlainMesh2D, GLPlainMesh3D

src/basic_types.jl

@@ -92,7 +92,6 @@ function Polytope(::Type{<: NNgon{N}}, P::Type{<: AbstractPoint{NDim, T}}) where
     return Ngon{NDim, T, N, P}
 end
 
-
 const LineP{Dim, T, P <: AbstractPoint{Dim, T}} = Ngon{Dim, T, 2, P}
 const Line{Dim, T} = LineP{Dim, T, Point{Dim, T}}
 
@@ -111,6 +110,17 @@ const Quadrilateral{Dim, T} = Ngon{Dim, T, 4, P} where P <: AbstractPoint{Dim, T
 Base.show(io::IO, x::Quadrilateral) = print(io, "Quad(", join(x, ", "), ")")
 Base.summary(io::IO, x::Type{<: Quadrilateral}) = print(io, "Quad")
 
+function coordinates(lines::AbstractArray{LineP{Dim, T, PointType}}) where {Dim, T, PointType}
+    if lines isa Base.ReinterpretArray
+        return coordinates(lines.parent)
+    else
+        result = PointType[]
+        for line in lines
+            append!(result, coordinates(line))
+        end
+        return result
+    end
+end
 
 """
 A `Simplex` is a generalization of an N-dimensional tetrahedra and can be thought
@@ -181,11 +191,11 @@ struct LineString{
     points::V
 end
 
-coordinates(x::LineString) = x.points
+coordinates(x::LineString) = coordinates(x.points)
 
 Base.copy(x::LineString) = LineString(copy(x.points))
-Base.size(x::LineString) = size(coordinates(x))
-Base.getindex(x::LineString, i) = getindex(coordinates(x), i)
+Base.size(x::LineString) = size(getfield(x, :points))
+Base.getindex(x::LineString, i) = getindex(getfield(x, :points), i)
 
 function LineString(points::AbstractVector{LineP{Dim, T, P}}) where {Dim, T, P}
     return LineString{Dim, T, P, typeof(points)}(points)
@@ -276,6 +286,27 @@ function Polygon(exterior::AbstractVector{P}, faces::AbstractVector{<: LineFace}
     return Polygon(LineString(exterior, faces))
 end
 
+function Polygon(exterior::AbstractVector{P}, interior::AbstractVector{<:AbstractVector{P}}) where P <: AbstractPoint{Dim, T} where {Dim, T}
+    ext = LineString(exterior)
+    # We need to take extra care for empty interiors, since
+    # if we just map over it, it won't infer the element type correctly!
+    int = typeof(ext)[]
+    foreach(x-> push!(int, LineString(x)), interior)
+    return Polygon(ext, int)
+end
+
+function coordinates(polygon::Polygon{N, T, PointType}) where {N, T, PointType}
+    exterior = coordinates(polygon.exterior)
+    if isempty(polygon.interiors)
+        return exterior
+    else
+        result = PointType[]
+        append!(result, exterior)
+        foreach(x-> append!(result, coordinates(x)), polygon.interiors)
+        return result
+    end
+end
+
 """
     MultiPolygon(polygons::AbstractPolygon)
 """

src/interfaces.jl

@@ -75,13 +75,12 @@ faces(tesselation::Tesselation) = faces(tesselation.primitive, nvertices(tessela
 normals(tesselation::Tesselation) = normals(tesselation.primitive, nvertices(tesselation))
 texturecoordinates(tesselation::Tesselation) = texturecoordinates(tesselation.primitive, nvertices(tesselation))
 
-
 ## Decompose methods
 # Dispatch type to make `decompose(UV{Vec2f0}, primitive)` work
 # and to pass through tesselation information
 
 # Types that can be converted to a mesh via the functions below
-const Meshable{Dim, T} = Union{Tesselation{Dim, T}, Mesh{Dim, T},
+const Meshable{Dim, T} = Union{Tesselation{Dim, T}, Mesh{Dim, T}, AbstractPolygon{Dim, T},
                                GeometryPrimitive{Dim, T}, AbstractVector{<: AbstractPoint{Dim, T}}}
 
 struct UV{T} end
@@ -108,23 +107,14 @@ function decompose(::Type{Point}, primitive::Meshable{Dim, T}) where {Dim, T}
     return collect_with_eltype(Point{Dim, T}, metafree(coordinates(primitive)))
 end
 
-function decompose(::Type{T}, primitive) where {T}
-    return collect_with_eltype(T, primitive)
+function decompose(::Type{Point}, primitive::LineString{Dim, T}) where {Dim, T}
+    return collect_with_eltype(Point{Dim, T}, metafree(coordinates(primitive)))
 end
 
-function decompose(::Type{P}, pol::Polygon) where {P<:AbstractPoint}
-    if isempty(pol.interiors)
-        return decompose(P, pol.exterior)
-    else
-        arr = copy(decompose(P, pol.exterior))
-        for i in pol.interiors
-            append!(arr, decompose(P, i))
-        end
-        return arr
-    end
+function decompose(::Type{T}, primitive) where {T}
+    return collect_with_eltype(T, primitive)
 end
 
-decompose(::Type{P}, ls::LineString) where {P<:AbstractPoint} = ls.points.parent.data
 decompose_uv(primitive) = decompose(UV(), primitive)
 decompose_uvw(primitive) = decompose(UVW(), primitive)
 decompose_normals(primitive) = decompose(Normal(), primitive)

src/meshes.jl

@@ -14,7 +14,6 @@ function normals(mesh::AbstractMesh)
     return nothing
 end
 
-
 const GLTriangleElement = Triangle{3, Float32}
 const GLTriangleFace = TriangleFace{GLIndex}
 const PointWithUV{Dim, T} = PointMeta{Dim, T, Point{Dim, T}, (:uv,), Tuple{Vec{2, T}}}
@@ -144,14 +143,19 @@ Polygon triangluation!
 function mesh(polygon::AbstractVector{P}; pointtype=P, facetype=GLTriangleFace,
               normaltype=nothing) where {P<:AbstractPoint{2}}
 
+    return mesh(Polygon(polygon); pointtype, facetype, normaltype)
+end
+
+function mesh(polygon::AbstractPolygon{Dim, T}; pointtype=Point{Dim, T}, facetype=GLTriangleFace,
+              normaltype=nothing) where {Dim, T}
+
     faces = decompose(facetype, polygon)
     positions = decompose(pointtype, polygon)
 
     if normaltype !== nothing
         n = normals(positions, faces; normaltype=normaltype)
         positions = meta(positions; normals=n)
     end
-
     return Mesh(positions, faces)
 end
 

src/triangulation.jl

@@ -169,3 +169,66 @@ function decompose(::Type{FaceType}, points::AbstractArray{P}) where {P<:Abstrac
 
     return result
 end
+
+
+function earcut_triangulate(polygon::Vector{Vector{Point{2, Float64}}})
+    lengths = map(x-> UInt32(length(x)), polygon)
+    len = UInt32(length(lengths))
+    array = ccall(
+        (:u32_triangulate_f64, libearcut),
+        Tuple{Ptr{GLTriangleFace}, Cint},
+        (Ptr{Ptr{Float64}}, Ptr{UInt32}, UInt32),
+        polygon, lengths, len
+    )
+    return unsafe_wrap(Vector{GLTriangleFace}, array[1], array[2])
+end
+
+function earcut_triangulate(polygon::Vector{Vector{Point{2, Float32}}})
+    lengths = map(x-> UInt32(length(x)), polygon)
+    len = UInt32(length(lengths))
+    array = ccall(
+        (:u32_triangulate_f32, libearcut),
+        Tuple{Ptr{GLTriangleFace}, Cint},
+        (Ptr{Ptr{Float32}}, Ptr{UInt32}, UInt32),
+        polygon, lengths, len
+    )
+    return unsafe_wrap(Vector{GLTriangleFace}, array[1], array[2])
+end
+
+function earcut_triangulate(polygon::Vector{Vector{Point{2, Int64}}})
+    lengths = map(x-> UInt32(length(x)), polygon)
+    len = UInt32(length(lengths))
+    array = ccall(
+        (:u32_triangulate_i64, libearcut),
+        Tuple{Ptr{GLTriangleFace}, Cint},
+        (Ptr{Ptr{Int64}}, Ptr{UInt32}, UInt32),
+        polygon, lengths, len
+    )
+    return unsafe_wrap(Vector{GLTriangleFace}, array[1], array[2])
+end
+
+function earcut_triangulate(polygon::Vector{Vector{Point{2, Int32}}})
+    lengths = map(x-> UInt32(length(x)), polygon)
+    len = UInt32(length(lengths))
+    array = ccall(
+        (:u32_triangulate_i32, libearcut),
+        Tuple{Ptr{GLTriangleFace}, Cint},
+        (Ptr{Ptr{Int32}}, Ptr{UInt32}, UInt32),
+        polygon, lengths, len
+    )
+    return unsafe_wrap(Vector{GLTriangleFace}, array[1], array[2])
+end
+
+best_earcut_eltype(x) = Float64
+best_earcut_eltype(::Type{Float64}) = Float64
+best_earcut_eltype(::Type{<:AbstractFloat}) = Float32
+best_earcut_eltype(::Type{Int64}) = Int64
+best_earcut_eltype(::Type{Int32}) = Int32
+best_earcut_eltype(::Type{<:Integer}) = Int64
+
+function faces(polygon::Polygon{Dim, T}) where {Dim, T}
+    PT = Point{Dim, best_earcut_eltype(T)}
+    points = [decompose(PT, polygon.exterior)]
+    foreach(x-> push!(points, decompose(PT, x)), polygon.interiors)
+    return earcut_triangulate(points)
+end

src/viewtypes.jl

@@ -36,6 +36,8 @@ struct TupleView{T, N, Skip, A} <: AbstractVector{T}
     connect::Bool
 end
 
+coordinates(tw::TupleView) = coordinates(tw.data)
+
 function Base.size(x::TupleView{T, N, M}) where {T, N, M}
     nitems = length(x.data) ÷ (N - (N - M))
     nitems = nitems - max(N - M, 0)

test/runtests.jl

@@ -56,7 +56,7 @@ using GeometryBasics: attributes
         end
 
     end
-   
+
     @testset "polygon with metadata" begin
         polys = [Polygon(rand(Point{2, Float32}, 20)) for i in 1:10]
         pnames = [randstring(4) for i in 1:10]
@@ -68,7 +68,7 @@ using GeometryBasics: attributes
         multipoly = MultiPolygonMeta(polys, name = pnames, value = numbers, category = bin)
         @test multipoly isa AbstractVector
         @test poly isa GeometryBasics.AbstractPolygon
-        
+
         @test GeometryBasics.getcolumn(poly, :name) == pnames[1]
         @test GeometryBasics.MetaFree(PolygonMeta) == Polygon
 
@@ -92,7 +92,7 @@ using GeometryBasics: attributes
         @test metafree(pm) === p
         @test propertynames(pm) == (:position, :a, :b)
     end
-    
+
     @testset "MultiPoint with metadata" begin
         p = collect(Point{2, Float64}(x, x+1) for x in 1:5)
         @test p isa AbstractVector
@@ -144,18 +144,18 @@ end
         multipoly = MetaT(multipol, name = pnames, value = numbers, category = bin)
         @test multipoly isa MetaT
         @test poly isa MetaT
-        
+
         @test GeometryBasics.getcolumn(poly, :name) == pnames[1]
         @test GeometryBasics.getcolumn(multipoly, :name) == pnames
-        
+
         meta_p = MetaT(polys[1], boundingbox=Rect(0, 0, 2, 2))
         @test meta_p.boundingbox === Rect(0, 0, 2, 2)
         @test GeometryBasics.metafree(meta_p) == polys[1]
         @test GeometryBasics.metafree(poly) == polys[1]
         @test GeometryBasics.metafree(multipoly) == multipol
         @test GeometryBasics.meta(meta_p) == (boundingbox = GeometryBasics.HyperRectangle{2,Int64}([0, 0], [2, 2]),)
         @test GeometryBasics.meta(poly) == (name = pnames[1], value = 0.0, category = bin[1])
-        @test GeometryBasics.meta(multipoly) == (name = pnames, value = numbers, category = bin)    
+        @test GeometryBasics.meta(multipoly) == (name = pnames, value = numbers, category = bin)
     end
 
     @testset "MetaT{Point}" begin
@@ -171,7 +171,7 @@ end
         @test GeometryBasics.metafree(pm) == p
         @test GeometryBasics.meta(pm) == (a = 1, b = 2)
     end
-    
+
     @testset "MetaT{MultiPoint}" begin
         p = collect(Point{2, Float64}(x, x+1) for x in 1:5)
         @test p isa AbstractVector
@@ -623,7 +623,7 @@ end
 
 @testset "MetaT and heterogeneous data" begin
     ls = [LineString([Point(i, (i+1)^2/6), Point(i*0.86,i+5), Point(i/3, i/7)]) for i in 1:10]
-    mls = MultiLineString([LineString([Point(i+1, (i)^2/6), Point(i*0.75,i+8), Point(i/2.5, i/6.79)]) for i in 5:10]) 
+    mls = MultiLineString([LineString([Point(i+1, (i)^2/6), Point(i*0.75,i+8), Point(i/2.5, i/6.79)]) for i in 5:10])
     poly = Polygon(Point{2, Int}[(40, 40), (20, 45), (45, 30), (40, 40)])
     geom = [ls..., mls, poly]
     prop = Any[(country_states = "India$(i)", rainfall = (i*9)/2) for i in 1:11]
@@ -637,17 +637,14 @@ end
     @test propertynames(sa) === (:main, :country_states, :rainfall)
     @test getproperty(sa, :country_states) isa Array{Any}
     @test getproperty(sa, :main) == geom
-    
-    @test GeometryBasics.getnamestypes(typeof(feat[1])) == 
+
+    @test GeometryBasics.getnamestypes(typeof(feat[1])) ==
     (LineString{2,Float64,Point{2,Float64},Base.ReinterpretArray{GeometryBasics.Ngon{2,Float64,2,Point{2,Float64}},1,Tuple{Point{2,Float64},Point{2,Float64}},TupleView{Tuple{Point{2,Float64},Point{2,Float64}},2,1,Array{Point{2,Float64},1}}}},
     (:country_states, :rainfall), Tuple{String,Float64})
-    
-    @test StructArrays.staticschema(typeof(feat[1])) == 
-    NamedTuple{(:main, :country_states, :rainfall),Tuple{LineString{2,Float64,Point{2,Float64},Base.ReinterpretArray{GeometryBasics.Ngon{2,Float64,2,Point{2,Float64}},1,Tuple{Point{2,Float64},Point{2,Float64}},TupleView{Tuple{Point{2,Float64},Point{2,Float64}},2,1,Array{Point{2,Float64},1}}}},
-    String,Float64}}
+
 
     @test StructArrays.createinstance(typeof(feat[1]), LineString([Point(1, (2)^2/6), Point(1*0.86,6), Point(1/3, 1/7)]), "Mumbai", 100) isa typeof(feat[1])
-    
+
     @test Base.getindex(feat[1], 1) isa Line
     @test Base.size(feat[1]) == (2,)
 end