Improve test coverage

docs/make.jl

@@ -13,8 +13,9 @@ makedocs(format=Documenter.HTML(prettyurls=get(ENV, "CI", "false") == "true"),
                 "polygons.md",
                 "meshes.md",
                 "decomposition.md",
+                "boundingboxes.md",
                 "static_array_types.md",
-                "api.md"
+                "api.md",
                ],
          modules=[GeometryBasics])
 

docs/src/boundingboxes.md

@@ -0,0 +1,31 @@
+# Bounding Boxes
+
+You can generate an axis aligned bounding box for any `AbstractGeometry` by calling `Rect(geom)`.
+Depending on the object this will either rely on `coordinates(geom)` or a specialized method.
+You can also create a bounding box of set dimension or type by adding the related parameters.
+
+
+```@repl
+using GeometryBasics
+
+s = Circle(Point2f(0), 1f0)
+Rect(s) # specialized, exact bounding box
+Rect3(s)
+Rect3d(s)
+RectT{Float64}(s)
+Rect(GeometryBasics.mesh(s)) # using generated coordinates in mesh
+```
+
+## Extending
+
+If you want to add a specialized bounding box method you should implement `Rect{N, T}(geom) = ...`.
+All other methods funnel into that one, defaulting to the same `N, T` that the given `AbstractGeometry{N, T}` has.
+GeometryBasics allows the user given dimension `N` to be smaller or equal to that of the geometry.
+This is checked with `GeometryBasics.bbox_dim_check(user_dim, geom_dim)` which you may reuse.
+
+```julia
+function Rect{N, T}(a::HyperSphere{N2}) where {N, N2, T}
+    GeometryBasics.bbox_dim_check(N, N2)
+    return Rect{N, T}(minimum(a), widths(a))
+end
+```

src/GeometryBasics.jl

@@ -41,7 +41,7 @@ export AbstractFace, TriangleFace, QuadFace, GLTriangleFace
 export OffsetInteger, ZeroIndex, OneIndex, GLIndex
 export decompose, coordinates, faces, normals, decompose_uv, decompose_normals,
        texturecoordinates, vertex_attributes
-export expand_faceviews
+export expand_faceviews, split_mesh, remove_duplicates
 export face_normals
 export Tessellation, Normal, UV, UVW
 export AbstractMesh, Mesh, MetaMesh, FaceView
@@ -58,7 +58,7 @@ export uv_mesh, normal_mesh, uv_normal_mesh
 export height, origin, radius, width, widths
 export HyperSphere, Circle, Sphere
 export Cylinder, Pyramid, extremity
-export HyperRectangle, Rect, Rect2, Rect3, Recti, Rect2i, Rect3i, Rectf, Rect2f, Rect3f, Rectd, Rect2d, Rect3d
+export HyperRectangle, Rect, Rect2, Rect3, Recti, Rect2i, Rect3i, Rectf, Rect2f, Rect3f, Rectd, Rect2d, Rect3d, RectT
 export before, during, meets, overlaps, intersects, finishes
 export centered, direction, area, volume, update
 export max_dist_dim, max_euclidean, max_euclideansq, min_dist_dim, min_euclidean

src/basic_types.jl

@@ -321,6 +321,17 @@ function coordinates(polygon::Polygon{N,T}) where {N,T}
     end
 end
 
+function Base.promote_rule(::Type{Polygon{N, T1}}, ::Type{Polygon{N, T2}}) where {N, T1, T2}
+    return Polygon{N, promote_type(T1, T2)}
+end
+
+function Base.convert(::Type{Polygon{N, T}}, poly::Polygon{N}) where {N, T}
+    return Polygon(
+        convert(Vector{Point{N, T}}, poly.exterior),
+        convert(Vector{Vector{Point{N, T}}}, poly.interiors),
+    )
+end
+
 """
     MultiPolygon(polygons::AbstractPolygon)
 
@@ -337,6 +348,7 @@ end
 Base.getindex(mp::MultiPolygon, i) = mp.polygons[i]
 Base.size(mp::MultiPolygon) = size(mp.polygons)
 Base.length(mp::MultiPolygon) = length(mp.polygons)
+Base.:(==)(a::MultiPolygon, b::MultiPolygon) = a.polygons == b.polygons
 
 """
     LineString(points::AbstractVector{<:Point})
@@ -361,6 +373,7 @@ end
 Base.getindex(ms::MultiLineString, i) = ms.linestrings[i]
 Base.size(ms::MultiLineString) = size(ms.linestrings)
 Base.length(mpt::MultiLineString) = length(mpt.linestrings)
+Base.:(==)(a::MultiLineString, b::MultiLineString) = a.linestrings == b.linestrings
 
 """
     MultiPoint(points::AbstractVector{AbstractPoint})
@@ -639,6 +652,65 @@ function Base.:(==)(a::Mesh, b::Mesh)
            (faces(a) == faces(b)) && (a.views == b.views)
 end
 
+"""
+    strictly_equal_face_vertices(a::Mesh, b::Mesh)
+
+Checks whether mesh a and b are equal in terms of vertices used in their faces.
+This allows for vertex data and indices to be synchronously permuted. For
+example, this will recognize
+```
+a = Mesh([a, b, c], [GLTriangleFace(1,2,3)])
+b = Mesh([a, c, b], [GLTriangleFace(1,3,2)])
+```
+as equal, because while the positions and faces have different orders the vertices
+in the face are the same:
+```
+[a, c, b][[1, 3, 2]] == [a, b, c] == [a, b, c][[1,2,3]]
+```
+
+This still returns false if the order of faces is permuted, e.g.
+`Mesh(ps, [f1, f2]) != Mesh(ps, [f2, f1])`. It also returns false if vertices are
+cyclically permuted within a face, i.e. `ps[[1,2,3]] != ps[[2,3,1]]`.
+"""
+function strictly_equal_face_vertices(a::Mesh, b::Mesh)
+    # Quick checks
+    if propertynames(a) != propertynames(b) || length(faces(a)) != length(faces(b))
+        return false
+    end
+
+    N = length(faces(a))
+    # for views we want to ignore empty ranges (they don't represent geometry)
+    # and treat 1:N as no range (as that is used interchangeably)
+    views1 = filter(view -> length(view) > 0 && (minimum(view) > 1 || maximum(view) < N), a.views)
+    views2 = filter(view -> length(view) > 0 && (minimum(view) > 1 || maximum(view) < N), b.views)
+    views1 != views2 && return false
+
+    # TODO: Allow different face orders & cyclic permutation within faces.
+    # E.g. use hash.(data[face]), cyclically permute min to front, hash result
+    # and add them to heaps (or sets?) so we can compare them at the end
+    # That should probably be another function as it's probably a significant
+    # step up in overhead?
+    for (attrib1, attrib2) in zip(vertex_attributes(a), vertex_attributes(b))
+        if attrib1 isa FaceView
+            if !(attrib2 isa FaceView) || length(faces(attrib1)) != length(faces(attrib2))
+                return false
+            end
+            for (f1, f2) in zip(faces(attrib1), faces(attrib2))
+                values(attrib1)[f1] == values(attrib2)[f2] || return false
+            end
+        else
+            if attrib2 isa FaceView
+                return false
+            end
+            for (f1, f2) in zip(faces(a), faces(b))
+                attrib1[f1] == attrib2[f2] || return false
+            end
+        end
+    end
+
+    return true
+end
+
 function Base.iterate(mesh::Mesh, i=1)
     return i - 1 < length(mesh) ? (mesh[i], i + 1) : nothing
 end
@@ -686,6 +758,12 @@ function Mesh(points::AbstractVector{<:Point}, faces::AbstractVector{<:Integer},
     return Mesh(points, connect(faces, facetype, skip))
 end
 
+# the method above allows Mesh(..., Face(...), ...) to work, but produces bad results
+# explicitly error here
+function Mesh(points::AbstractVector{<:Point}, faces::AbstractFace, args...; kwargs...)
+    throw(MethodError(Mesh, (points, faces, args...)))
+end
+
 function Mesh(; kwargs...)
     fs = faces(kwargs[:position]::FaceView)
     va = NamedTuple{keys(kwargs)}(map(keys(kwargs)) do k

src/boundingboxes.jl

@@ -1,18 +1,34 @@
-function Rect(geometry::AbstractArray{<:Point{N,T}}) where {N,T}
-    return Rect{N,T}(geometry)
+# Boundingbox-like Rect constructors
+
+Rect(p::AbstractGeometry{N, T}) where {N, T} = Rect{N, T}(p)
+RectT{T}(p::AbstractGeometry{N}) where {N, T} = Rect{N, T}(p)
+Rect{N}(p::AbstractGeometry{_N, T}) where {N, _N, T} = Rect{N, T}(p)
+
+Rect(p::AbstractArray{<: VecTypes{N, T}}) where {N, T} = Rect{N, T}(p)
+RectT{T}(p::AbstractArray{<: VecTypes{N}}) where {N, T} = Rect{N, T}(p)
+Rect{N}(p::AbstractArray{<: VecTypes{_N, T}}) where {N, _N, T} = Rect{N, T}(p)
+
+# Implementations
+# Specialize fully typed Rect constructors
+Rect{N, T}(p::AbstractGeometry) where {N, T} = Rect{N, T}(coordinates(p))
+
+function bbox_dim_check(trg, src::Integer)
+    @assert trg isa Integer "Rect{$trg, $T1} is invalid. This may have happened due to calling Rect{$N1}(obj) to get a bounding box."
+    if trg < src
+        throw(ArgumentError("Cannot construct a $trg dimensional bounding box from $src dimensional Points. ($trg must be ≥ $src)"))
+    end
 end
 
 """
-    Rect(points::AbstractArray{<: Point})
+    Rect(points::AbstractArray{<: VecTypes})
 
 Construct a bounding box containing all the given points.
 """
-function Rect{N1,T1}(geometry::AbstractArray{PT}) where {N1,T1,PT<:Point}
-    N2, T2 = length(PT), eltype(PT)
-    @assert N1 >= N2
+function Rect{N1, T1}(points::AbstractArray{<: VecTypes{N2, T2}}) where {N1, T1, N2, T2}
+    bbox_dim_check(N1, N2)
     vmin = Point{N2,T2}(typemax(T2))
     vmax = Point{N2,T2}(typemin(T2))
-    for p in geometry
+    for p in points
         vmin, vmax = _minmax(p, vmin, vmax)
     end
     o = vmin
@@ -25,29 +41,25 @@ function Rect{N1,T1}(geometry::AbstractArray{PT}) where {N1,T1,PT<:Point}
     end
 end
 
+
 """
     Rect(primitive::GeometryPrimitive)
 
 Construct a bounding box for the given primitive.
 """
-function Rect(primitive::GeometryPrimitive{N,T}) where {N,T}
-    return Rect{N,T}(primitive)
-end
-
-function Rect{T}(primitive::GeometryPrimitive{N,T}) where {N,T}
-    return Rect{N,T}(primitive)
-end
-
-function Rect{T}(a::Pyramid) where {T}
-    w, h = a.width / T(2), a.length
+function Rect{N, T}(a::Pyramid) where {N, T}
+    bbox_dim_check(N, 3)
+    w, h = a.width, a.length
     m = Vec{3,T}(a.middle)
-    return Rect{T}(m .- Vec{3,T}(w, w, 0), m .+ Vec{3,T}(w, w, h))
+    return Rect{N, T}(m .- Vec{3,T}(w / T(2), w / T(2), 0), Vec{3,T}(w, w, h))
 end
 
-function Rect{T}(a::Sphere) where {T}
-    mini, maxi = extrema(a)
-    return Rect{T}(mini, maxi .- mini)
+function Rect{N, T}(a::HyperSphere{N2}) where {N, N2, T}
+    bbox_dim_check(N, N2)
+    return Rect{N, T}(minimum(a), widths(a))
 end
 
-Rect{T}(a) where {T} = Rect{T}(coordinates(a))
-Rect{N,T}(a) where {N,T} = Rect{N,T}(coordinates(a))
+# TODO: exact implementation that doesn't rely on coordinates
+# function Rect{N, T}(a::Cylinder) where {N, T}
+#     return Rect{N, T}(...)
+# end

src/fixed_arrays.jl

@@ -57,7 +57,8 @@ macro fixed_vector(name_parent)
         function $(name){S}(x::T) where {S,T <: Tuple}
             return $(name){S,StaticArrays.promote_tuple_eltype(T)}(x)
         end
-        $(name){S,T}(x::StaticVector) where {S,T} = $(name){S,T}(Tuple(x))
+        $(name){S,T}(x::StaticVector{S}) where {S,T} = $(name){S,T}(Tuple(x))
+        $(name){S,T}(x::StaticVector) where {S,T} = $(name){S,T}(ntuple(i -> x[i], S))
 
         @generated function (::Type{$(name){S,T}})(x::$(name)) where {S,T}
             idx = [:(x[$i]) for i in 1:S]
@@ -139,7 +140,7 @@ const VecTypes{N,T} = Union{StaticVector{N,T},NTuple{N,T}}
 const Vecf{N} = Vec{N,Float32}
 const PointT{T} = Point{N,T} where N
 const Pointf{N} = Point{N,Float32}
-    
+
 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)
@@ -177,9 +178,9 @@ export Vecf, Pointf
     Vec{N, T}(args...)
     Vec{N, T}(args::Union{AbstractVector, Tuple, NTuple, StaticVector})
 
-Constructs a Vec of length `N` from the given arguments. 
+Constructs a Vec of length `N` from the given arguments.
 
-Note that Point and Vec don't follow strict mathematical definitions. Instead 
+Note that Point and Vec don't follow strict mathematical definitions. Instead
 we allow them to be used interchangeably.
 
 ## Aliases
@@ -197,9 +198,9 @@ Vec
     Point{N, T}(args...)
     Point{N, T}(args::Union{AbstractVector, Tuple, NTuple, StaticVector})
 
-Constructs a Point of length `N` from the given arguments. 
+Constructs a Point of length `N` from the given arguments.
 
-Note that Point and Vec don't follow strict mathematical definitions. Instead 
+Note that Point and Vec don't follow strict mathematical definitions. Instead
 we allow them to be used interchangeably.
 
 ## Aliases

src/lines.jl

@@ -7,7 +7,7 @@ Returns `(intersection_found::Bool, intersection_point::Point)`
 """
 # 2D Line-segment intersection algorithm by Paul Bourke and many others.
 # http://paulbourke.net/geometry/pointlineplane/
-function intersects(a::Line{2,T1}, b::Line{2,T2}) where {T1,T2}
+function intersects(a::Line{2,T1}, b::Line{2,T2}; eps = 0) where {T1,T2}
     T = promote_type(T1, T2)
     p0 = zero(Point2{T})
 
@@ -34,7 +34,7 @@ function intersects(a::Line{2,T1}, b::Line{2,T2}) where {T1,T2}
 
     # Values between [0, 1] mean the intersection point of the lines rests on
     # both of the line segments.
-    if 0 <= unknown_a <= 1 && 0 <= unknown_b <= 1
+    if eps <= unknown_a <= 1-eps && eps <= unknown_b <= 1-eps
         # Substituting an unknown back lets us find the intersection point.
         x = x1 + (unknown_a * (x2 - x1))
         y = y1 + (unknown_a * (y2 - y1))
@@ -62,27 +62,39 @@ end
 """
     self_intersections(points::AbstractVector{<:Point})
 
-Finds all self intersections of polygon `points`
+Finds all self intersections of in a continuous line described by `points`.
+Returns a Vector of indices where each pair `v[2i], v[2i+1]` refers two
+intersecting line segments by their first point, and a Vector of intersection
+points.
+
+Note that if two points are the same, they will generate a self intersection
+unless they are consecutive segments. (The first and last point are assumed to
+be shared between the first and last segment.)
 """
-function self_intersections(points::AbstractVector{<:Point})
+function self_intersections(points::AbstractVector{<:VecTypes{D, T}}) where {D, T}
+    ti, sections = _self_intersections(points)
+    # convert array of tuples to flat array
+    return [x for t in ti for x in t], sections
+end
+
+function _self_intersections(points::AbstractVector{<:VecTypes{D, T}}) where {D, T}
     sections = similar(points, 0)
-    intersections = Int[]
-
-    wraparound(i) = mod1(i, length(points) - 1)
-
-    for (i, (a, b)) in enumerate(consecutive_pairs(points))
-        for (j, (a2, b2)) in enumerate(consecutive_pairs(points))
-            is1, is2 = wraparound(i + 1), wraparound(i - 1)
-            if i != j &&
-               is1 != j &&
-               is2 != j &&
-               !(i in intersections) &&
-               !(j in intersections)
-                intersected, p = intersects(Line(a, b), Line(a2, b2))
-                if intersected
-                    push!(intersections, i, j)
-                    push!(sections, p)
-                end
+    intersections = Tuple{Int, Int}[]
+
+    N = length(points)
+
+    for i in 1:N-3
+        a = points[i]; b = points[i+1]
+        # i+1 == j describes consecutive segments which are always "intersecting"
+        # at point i+1/j. Skip those (start at i+2)
+        # Special case: We assume points[1] == points[end] so 1 -> 2 and N-1 -> N
+        # always "intersect" at 1/N. Skip this too (end at N-2 in this case)
+        for j in i+2 : N-1 - (i == 1)
+            a2 = points[j]; b2 = points[j+1]
+            intersected, p = intersects(Line(a, b), Line(a2, b2))
+            if intersected
+                push!(intersections, (i, j))
+                push!(sections, p)
             end
         end
     end
@@ -95,15 +107,14 @@ end
 Splits polygon `points` into it's self intersecting parts. Only 1 intersection
 is handled right now.
 """
-function split_intersections(points::AbstractVector{<:Point})
-    intersections, sections = self_intersections(points)
+function split_intersections(points::AbstractVector{<:VecTypes{N, T}}) where {N, T}
+    intersections, sections = _self_intersections(points)
     return if isempty(intersections)
         return [points]
-    elseif length(intersections) == 2 && length(sections) == 1
-        a, b = intersections
+    elseif length(intersections) == 1 && length(sections) == 1
+        a, b = intersections[1]
         p = sections[1]
-        a, b = min(a, b), max(a, b)
-        poly1 = simple_concat(points, (a + 1):(b - 1), p)
+        poly1 = simple_concat(points, (a + 1):b, p)
         poly2 = simple_concat(points, (b + 1):(length(points) + a), p)
         return [poly1, poly2]
     else