Merge branch 'master' into sd/simple-mesh

Author: SimonDanisch

.github/workflows/TagBot.yml

@@ -4,6 +4,22 @@ on:
     types:
       - created
   workflow_dispatch:
+    inputs:
+      lookback:
+        default: 3
+permissions:
+  actions: read
+  checks: read
+  contents: write
+  deployments: read
+  issues: read
+  discussions: read
+  packages: read
+  pages: read
+  pull-requests: read
+  repository-projects: read
+  security-events: read
+  statuses: read
 jobs:
   TagBot:
     if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'

Project.toml

@@ -1,17 +1,19 @@
 name = "GeometryBasics"
 uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 authors = ["SimonDanisch <[email protected]>"]
-version = "0.4.5"
+version = "0.4.9"
 
 [deps]
 EarCut_jll = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
+Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 EarCut_jll = "2"
+Extents = "0.1"
 GeoInterface = "1.0.1"
 IterTools = "1.3.0"
 julia = "1.6"

src/GeometryBasics.jl

@@ -2,6 +2,7 @@ module GeometryBasics
 
 using IterTools, LinearAlgebra
 using GeoInterface
+import Extents
 using EarCut_jll
 import Base: *
 

src/geointerface.jl

@@ -2,12 +2,14 @@
 
 GeoInterface.isgeometry(::Type{<:AbstractGeometry}) = true
 GeoInterface.isgeometry(::Type{<:AbstractFace}) = true
+
 GeoInterface.isgeometry(::Type{<:Point}) = true
-GeoInterface.isgeometry(::Type{<:AbstractVector{<:AbstractGeometry}}) = true
-GeoInterface.isgeometry(::Type{<:AbstractVector{<:Point}}) = true
-GeoInterface.isgeometry(::Type{<:AbstractVector{<:LineString}}) = true
-GeoInterface.isgeometry(::Type{<:AbstractVector{<:AbstractPolygon}}) = true
-GeoInterface.isgeometry(::Type{<:AbstractVector{<:AbstractFace}}) = true
+GeoInterface.isgeometry(::Type{<:AbstractMesh}) = true
+GeoInterface.isgeometry(::Type{<:AbstractPolygon}) = true
+GeoInterface.isgeometry(::Type{<:LineString}) = true
+GeoInterface.isgeometry(::Type{<:MultiPoint}) = true
+GeoInterface.isgeometry(::Type{<:MultiLineString}) = true
+GeoInterface.isgeometry(::Type{<:MultiPolygon}) = true
 GeoInterface.isgeometry(::Type{<:Mesh}) = true
 
 GeoInterface.geomtrait(::Point) = PointTrait()
@@ -23,6 +25,7 @@ GeoInterface.geomtrait(::AbstractMesh) = PolyhedralSurfaceTrait()
 # GeoInterface calls this method in `GeoInterface.convert(GeometryBasics, ...)`
 geointerface_geomtype(::GeoInterface.PointTrait) = Point
 geointerface_geomtype(::GeoInterface.MultiPointTrait) = MultiPoint
+geointerface_geomtype(::GeoInterface.LineTrait) = Line
 geointerface_geomtype(::GeoInterface.LineStringTrait) = LineString
 geointerface_geomtype(::GeoInterface.MultiLineStringTrait) = MultiLineString
 geointerface_geomtype(::GeoInterface.PolygonTrait) = Polygon
@@ -59,8 +62,7 @@ GeoInterface.getgeom(::MultiPointTrait, g::MultiPoint, i::Int) = g[i]
 function GeoInterface.ngeom(::MultiLineStringTrait, g::MultiLineString)
     return length(g)
 end
-function GeoInterface.getgeom(::MultiLineStringTrait, g::MultiLineString,
-                              i::Int)
+function GeoInterface.getgeom(::MultiLineStringTrait, g::MultiLineString, i::Int)
     return g[i]
 end
 GeoInterface.ncoord(::MultiLineStringTrait, g::MultiLineString{Dim}) where {Dim} = Dim
@@ -91,13 +93,29 @@ GeoInterface.ngeom(::PolyhedralSurfaceTrait, g::AbstractMesh) = length(g)
 GeoInterface.getgeom(::PolyhedralSurfaceTrait, g::AbstractMesh, i) = g[i]
 
 function GeoInterface.convert(::Type{Point}, type::PointTrait, geom)
-    dim = Int(ncoord(geom))
-    return Point{dim, Float64}(GeoInterface.coordinates(geom))
+    x, y = GeoInterface.x(geom), GeoInterface.y(geom)
+    if GeoInterface.is3d(geom)
+        z = GeoInterface.z(geom)
+        T = promote_type(typeof(x), typeof(y), typeof(z))
+        return Point{3,T}(x, y, z)
+    else
+        GeoInterface.x(geom), GeoInterface.y(geom)
+        T = promote_type(typeof(x), typeof(y))
+        return Point{2,T}(x, y)
+    end
 end
 
 function GeoInterface.convert(::Type{LineString}, type::LineStringTrait, geom)
-    dim = Int(ncoord(geom))
-    return LineString([Point{dim, Float64}(GeoInterface.coordinates(p)) for p in getgeom(geom)])
+    g1 = getgeom(geom, 1)
+    x, y = GeoInterface.x(g1), GeoInterface.y(g1)
+    if GeoInterface.is3d(geom)
+        z = GeoInterface.z(g1)
+        T = promote_type(typeof(x), typeof(y), typeof(z))
+        return LineString([Point{3,T}(GeoInterface.x(p), GeoInterface.y(p), GeoInterface.z(p)) for p in getgeom(geom)])
+    else
+        T = promote_type(typeof(x), typeof(y))
+        return LineString([Point{2,T}(GeoInterface.x(p), GeoInterface.y(p)) for p in getgeom(geom)])
+    end
 end
 
 function GeoInterface.convert(::Type{Polygon}, type::PolygonTrait, geom)
@@ -106,22 +124,35 @@ function GeoInterface.convert(::Type{Polygon}, type::PolygonTrait, geom)
     if GeoInterface.nhole(geom) == 0
         return Polygon(exterior)
     else
-        interiors = GeoInterface.convert.(LineString, Ref(t), GeoInterface.gethole(geom))
+        interiors = map(h -> GeoInterface.convert(LineString, t, h), GeoInterface.gethole(geom))
         return Polygon(exterior, interiors)
     end
 end
 
 function GeoInterface.convert(::Type{MultiPoint}, type::MultiPointTrait, geom)
-    dim = Int(ncoord(geom))
-    return MultiPoint([Point{dim, Float64}(GeoInterface.coordinates(p)) for p in getgeom(geom)])
+    g1 = getgeom(geom, 1)
+    x, y = GeoInterface.x(g1), GeoInterface.y(g1)
+    if GeoInterface.is3d(geom)
+        z = GeoInterface.z(g1)
+        T = promote_type(typeof(x), typeof(y), typeof(z))
+        return MultiPoint([Point{3,T}(GeoInterface.x(p), GeoInterface.y(p), GeoInterface.z(p)) for p in getgeom(geom)])
+    else
+        T = promote_type(typeof(x), typeof(y))
+        return MultiPoint([Point{2,T}(GeoInterface.x(p), GeoInterface.y(p)) for p in getgeom(geom)])
+    end
 end
 
 function GeoInterface.convert(::Type{MultiLineString}, type::MultiLineStringTrait, geom)
     t = LineStringTrait()
-    return MultiLineString([GeoInterface.convert(LineString, t, l) for l in getgeom(geom)])
+    return MultiLineString(map(l -> GeoInterface.convert(LineString, t, l), getgeom(geom)))
 end
 
 function GeoInterface.convert(::Type{MultiPolygon}, type::MultiPolygonTrait, geom)
     t = PolygonTrait()
-    return MultiPolygon([GeoInterface.convert(Polygon, t, poly) for poly in getgeom(geom)])
+    return MultiPolygon(map(poly -> GeoInterface.convert(Polygon, t, poly), getgeom(geom)))
 end
+
+function Extents.extent(rect::Rect2)
+    (xmin, ymin), (xmax, ymax) = extrema(rect)
+    return Extents.Extent(X=(xmin, xmax), Y=(ymin, ymax))
+end

src/geometry_primitives.jl

@@ -75,15 +75,33 @@ collect_with_eltype(::Type{T}, vec::Vector{T}) where {T} = vec
 collect_with_eltype(::Type{T}, vec::AbstractVector{T}) where {T} = collect(vec)
 
 function collect_with_eltype(::Type{T}, iter) where {T}
-    # TODO we could be super smart about allocating the right length
-    # but its kinda annoying, since e.g. T == Triangle and first(iter) isa Quad
-    # will need double the length etc - but could all be figured out ;)
-    result = T[]
+    isempty(iter) && return T[]
+    # We need to get `eltype` information from `iter`, it seems to be `Any`
+    # most of the time so the eltype checks here don't actually work
+    l = if Base.IteratorSize(iter) isa Union{Base.HasShape,Base.HasLength}
+        if Base.IteratorEltype(iter) isa Base.HasEltype && isconcretetype(eltype(iter))
+            # Work out the exact length
+            length(convert_simplex(T, first(iter))) * length(iter)
+        else
+            # We know it is at least the length of iter,
+            # after that we will `push!` if we have to
+            length(iter)
+        end
+    else
+        0
+    end
+    n = 0
+    result = Vector{T}(undef, l)
     for element in iter
         # convert_simplex always returns a tuple,
         # so that e.g. convert(Triangle, quad) can return 2 elements
         for telement in convert_simplex(T, element)
-            push!(result, telement)
+            n += 1
+            if n > l
+                push!(result, telement)
+            else
+                result[n] = telement
+            end
         end
     end
     return result

src/meshes.jl

@@ -45,8 +45,25 @@ function triangle_mesh(primitive::AbstractGeometry{N})::TriangleMesh{N} where {N
     return mesh(primitive; pointtype=Point{N, Float32})
 end
 
-function triangle_mesh(primitive::AbstractVector{<: Point2})::TriangleMesh{2}
-    return mesh(Polygon(primitive))
+
+pointtype(x::Mesh) = eltype(decompose(Point, x))
+facetype(x::Mesh) = eltype(faces(x))
+
+function triangle_mesh(primitive::Mesh{N}) where {N}
+    # already target type:
+    if pointtype(primitive) === Point{N,Float32} && GLTriangleFace === facetype(primitive)
+        return primitive
+    else
+        return mesh(primitive; pointtype=Point{N,Float32}, facetype=GLTriangleFace)
+    end
+end
+
+function triangle_mesh(primitive::Meshable{N}; nvertices=nothing) where {N}
+    if nvertices !== nothing
+        @warn("nvertices argument deprecated. Wrap primitive in `Tesselation(primitive, nvertices)`")
+        primitive = Tesselation(primitive, nvertices)
+    end
+    return mesh(primitive; pointtype=Point{N,Float32}, facetype=GLTriangleFace)
 end
 
 function uv_mesh(primitive::AbstractGeometry{N,T}) where {N,T}
@@ -99,12 +116,17 @@ function Base.merge(meshes::AbstractVector{<:Mesh})
     elseif length(meshes) == 1
         return meshes[1]
     else
-        m1 = meshes[1]
-        ps = copy(coordinates(m1))
-        fs = copy(faces(m1))
+        ps = reduce(vcat, coordinates.(meshes))
+        fs = reduce(vcat, faces.(meshes))
+        idx = length(faces(meshes[1]))
+        offset = length(coordinates(meshes[1]))
         for mesh in Iterators.drop(meshes, 1)
-            append!(fs, map(f -> f .+ length(ps), faces(mesh)))
-            append!(ps, coordinates(mesh))
+            N = length(faces(mesh))
+            for i = idx .+ (1:N)
+                fs[i] = fs[i] .+ offset
+            end
+            idx += N
+            offset += length(coordinates(mesh))
         end
         return Mesh(ps, fs)
     end

test/geointerface.jl

@@ -1,17 +1,27 @@
 @testset "Basic types" begin
     point = Point(2, 3)
+    @test geomtrait(point) isa PointTrait
     @test testgeometry(point)
     @test ncoord(point) == 2
     @test getcoord(point, 2) == 3
     @test GeoInterface.coordinates(point) == [2, 3]
 
+    line = Line(Point(2, 3), Point(4, 5))
+    @test geomtrait(line) isa LineTrait
+    @test testgeometry(line)
+    @test ngeom(line) == 2
+    @test getgeom(line, 2) == Point(4, 5)
+    @test GeoInterface.coordinates(line) == [[2, 3], [4, 5]]
+
     mp = MultiPoint([point, point])
+    @test geomtrait(mp) isa MultiPointTrait
     @test testgeometry(mp)
     @test ngeom(mp) == 2
     @test getgeom(mp, 2) == point
     @test GeoInterface.coordinates(mp) == [[2, 3], [2, 3]]
 
     linestring = LineString(Point{2,Int}[(10, 10), (20, 20), (10, 40)])
+    @test geomtrait(linestring) isa LineStringTrait
     @test testgeometry(linestring)
     @test ngeom(linestring) == 3
     @test ncoord(linestring) == 2
@@ -21,22 +31,26 @@
     @test GeoInterface.coordinates(linestring) == [[10, 10], [20, 20], [10, 40]]
 
     multilinestring = MultiLineString([linestring, linestring])
+    @test geomtrait(multilinestring) isa MultiLineStringTrait
     @test testgeometry(multilinestring)
     @test GeoInterface.coordinates(multilinestring) ==
           [[[10, 10], [20, 20], [10, 40]], [[10, 10], [20, 20], [10, 40]]]
     @test ncoord(multilinestring) == 2
 
     poly = Polygon(rand(Point{2,Float32}, 5), [rand(Point{2,Float32}, 5)])
+    @test geomtrait(poly) isa PolygonTrait
     @test testgeometry(poly)
     @test length(GeoInterface.coordinates(poly)) == 2
     @test length(GeoInterface.coordinates(poly)[1]) == 5
 
     triangle = Triangle(point, point, point)
+    @test geomtrait(triangle) isa PolygonTrait # ?? should it be a Triangle trait
     @test testgeometry(triangle)
     @test length(GeoInterface.coordinates(triangle)) == 1
     @test length(GeoInterface.coordinates(triangle)[1]) == 3
 
     polys = MultiPolygon([poly, poly])
+    @test geomtrait(polys) isa MultiPolygonTrait
     @test testgeometry(polys)
     @test length(GeoInterface.coordinates(polys)) == 2
     @test length(GeoInterface.coordinates(polys)[1]) == 2
@@ -79,20 +93,21 @@ end
     multilinestring_gb = GeoInterface.convert(GeometryBasics, multilinestring_json)
     multipolygon_gb = GeoInterface.convert(GeometryBasics, multipolygon_json)
     multipolygon_hole_gb = GeoInterface.convert(GeometryBasics, multipolygon_hole_json)
-
-    @test point_gb === Point{2, Float64}(30.1, 10.1)
-    @test point_3d_gb === Point{3, Float64}(30.1, 10.1, 5.1)
+    
+    @test point_gb === Point{2,Float32}(30.1, 10.1)
+    @test point_3d_gb === Point{3,Float32}(30.1, 10.1, 5.1)
     @test linestring_gb isa LineString
     # TODO, what should we do exactly with linestrings?
     # @test length(linestring_gb) == 2
     # @test eltype(linestring_gb) == Line{2, Float64}
+
     @test polygon_gb isa Polygon
     @test isempty(polygon_gb.interiors)
     @test polygon_hole_gb isa Polygon
     @test length(polygon_hole_gb.interiors) == 1
     @test multipoint_gb isa MultiPoint
     @test length(multipoint_gb) == 4
-    @test multipoint_gb[4] === Point{2, Float64}(30.1, 10.1)
+    @test multipoint_gb[4] === Point{2,Float32}(30.1, 10.1)
     @test multilinestring_gb isa MultiLineString
     @test length(multilinestring_gb) == 2
     @test multipolygon_gb isa MultiPolygon
@@ -102,3 +117,10 @@ end
     @test length(multipolygon_hole_gb[1].interiors) == 0
     @test length(multipolygon_hole_gb[2].interiors) == 1
 end
+
+@testset "Extent" begin
+    rect = Rect2f(Vec2f(0), Vec2f(1.0))
+    ext = extent(rect)
+    @test ext.X == (0.0f0, 1.0f0)
+    @test ext.Y == (0.0f0, 1.0f0)
+end

test/geometrytypes.jl

@@ -73,7 +73,7 @@ using Test, GeometryBasics
         @test faces == decompose(TriangleFace{Int}, Tesselation(s, 8))
 
         m = triangle_mesh(Tesselation(s, 8))
-
+        @test m === triangle_mesh(m)
         @test GeometryBasics.faces(m) == faces
         @test GeometryBasics.coordinates(m) ≈ positions
         m = normal_mesh(s)# just test that it works without explicit resolution parameter

test/runtests.jl

@@ -4,6 +4,7 @@ using LinearAlgebra
 using GeometryBasics: MetaMesh, add_meta, pop_meta
 using GeoInterface
 using GeoJSON
+using Extents
 
 @testset "GeometryBasics" begin
 @testset "algorithms" begin
@@ -147,6 +148,8 @@ end
 end
 
 @testset "decompose/triangulation" begin
+    @test isempty(decompose(Vec3f, []))
+    @test decompose(Vec3f, []) isa Vector{Vec3f}
     primitive = Sphere(Point3f(0), 1)
     @test ndims(primitive) === 3
     mesh = triangle_mesh(primitive)