Clean up Rect set operations (union, intersect, diff, overlaps)

.github/workflows/ci.yml

@@ -27,7 +27,7 @@ jobs:
           - os: macOS-latest
             arch: aarch64
             version: '1'
-          - os: macOS-13
+          - os: macOS-15
             arch: x64
             version: '1'
     steps:

src/GeometryBasics.jl

@@ -58,7 +58,7 @@ export height, origin, radius, width, widths
 export HyperSphere, Circle, Sphere, Cone
 export Cylinder, Pyramid, extremity
 export HyperRectangle, Rect, Rect2, Rect3, Recti, Rect2i, Rect3i, Rectf, Rect2f, Rect3f, Rectd, Rect2d, Rect3d, RectT
-export before, during, meets, overlaps, intersects, finishes
+export before, during, meets, overlaps, intersects, finishes, bbox_diff
 export centered, direction, area, volume, update
 export max_dist_dim, max_euclidean, max_euclideansq, min_dist_dim, min_euclidean
 export min_euclideansq, minmax_dist_dim, minmax_euclidean, minmax_euclideansq
@@ -68,7 +68,7 @@ if Base.VERSION >= v"1.8"
     include("precompiles.jl")
 end
 
-# Needed for GeometryBasicsGeoInterfaceExt. 
+# Needed for GeometryBasicsGeoInterfaceExt.
 # In future this can go away as can use Module dispatch.
 function geointerface_geomtype end
 

src/primitives/rectangles.jl

@@ -320,41 +320,74 @@ end
 # set operations
 
 """
-    isempty(h::Rect)
+    isempty(h::Rect[, volumetric = true])
 
-Return `true` if any of the widths of `h` are negative.
+Return `true` if the rectangle is empty.
+
+By default (`volumetric = true`) a `Rect{D}` is considered empty if its
+D-dimensional volume is empty, i.e. if any of the widths are 0. Alternatively
+(`volumetric = false`) all of its widths are zero.
 """
-Base.isempty(h::Rect{N,T}) where {N,T} = any(<(zero(T)), h.widths)
+function Base.isempty(h::Rect{N,T}, volumetric = true) where {N,T}
+    degenerate = any(<=(zero(T)), h.widths)
+    empty = all(<=(zero(T)), h.widths)
+    return ifelse(volumetric, degenerate, empty)
+end
 
 """
     union(r1::Rect{N}, r2::Rect{N})
 
 Returns a new `Rect{N}` which contains both r1 and r2.
 """
 function Base.union(h1::Rect{N}, h2::Rect{N}) where {N}
-    m = min.(minimum(h1), minimum(h2))
-    mm = max.(maximum(h1), maximum(h2))
-    return Rect{N}(m, mm - m)
+    isempty(h1, false) && return h2
+    isempty(h2, false) && return h1
+    mini = min.(minimum(h1), minimum(h2))
+    maxi = max.(maximum(h1), maximum(h2))
+    return Rect{N}(mini, maxi - mini)
 end
 
-# TODO: What should this be? The difference is "h2 - h1", which could leave an
-# L shaped cutout. Should we pad that back out into a full rect?
-# """
-#     diff(h1::Rect, h2::Rect)
+# TODO: Add a diff that returns the slabs created. This could be anywhere between
+# 0 (r2 fully covers r1) to 2*D (?) (r2 is fully inside r1)
+"""
+    bbox_diff(r1::Rect{N}, r2::Rect{N})
+
+Returns the bounding box of the difference "r1 - r2".
+"""
+function bbox_diff(a::Rect{D, T1}, b::Rect{D, T2}) where {D, T1, T2}
+    T = promote_type(T1, T2)
+    cut_left = minimum(a) .>= minimum(b)
+    cut_right = maximum(a) .<= maximum(b)
+    a_fully_inside_b = cut_left .&& cut_right
+    fully_outside = any((minimum(a) .> maximum(b)) .|| (maximum(a) .< minimum(b)))
+    N = sum(a_fully_inside_b)
+    if N == D # intersection is a
+        return Rect{D, T}()
+    end
+
+    mini = ifelse.(a_fully_inside_b, minimum(a), ifelse.(cut_right, minimum(a), maximum(b)))
+    maxi = ifelse.(a_fully_inside_b, maximum(a), ifelse.(cut_left, maximum(a), minimum(b)))
+    widths = maxi - mini
+    if (N == D - 1) && !fully_outside && all(>=(0), widths)
+        # one dimension is not fully cut && shapes instersect &&
+        # b does not bisect a (cut_left and cut_right are both false => mini, maxi become maxi, mini of b)
+        return Rect{D, T}(mini, maxi .- mini)
+    end
+    return Rect{D, T}(a)
+end
 
-# Perform a difference between two Rects.
-# """
-# diff(h1::Rect, h2::Rect) = h1
 
 """
     intersect(h1::Rect, h2::Rect)
 
 Perform a intersection between two Rects.
 """
-function Base.intersect(h1::Rect{N}, h2::Rect{N}) where {N}
+function Base.intersect(h1::Rect{N, T1}, h2::Rect{N, T2}) where {N, T1, T2}
+    T = promote_type(T1, T2)
+    overlaps(h1, h2) || return Rect{N, T}()
     m = max.(minimum(h1), minimum(h2))
     mm = min.(maximum(h1), maximum(h2))
-    return Rect{N}(m, mm - m)
+    return Rect{N, T}(m, mm - m)
 end
 
 function update(b::Rect{N,T}, v::VecTypes{N,T2}) where {N,T,T2}
@@ -364,12 +397,12 @@ end
 function update(b::Rect{N,T}, v::VecTypes{N,T}) where {N,T}
     m = min.(minimum(b), v)
     maxi = maximum(b)
-    mm = if any(isnan, maxi)
+    ws = if any(isnan, maxi)
         v - m
     else
         max.(v, maxi) - m
     end
-    return Rect{N,T}(m, mm)
+    return Rect{N,T}(m, ws)
 end
 
 # Min maximum distance functions between hrectangle and point for a given dimension
@@ -439,6 +472,8 @@ function minmax_euclidean(rect::Rect{N,T}, p::Union{VecTypes{N,T},Rect{N,T}}) wh
 end
 
 # http://en.wikipedia.org/wiki/Allen%27s_interval_algebra
+# These don't really make sense for rectangles though?
+
 function before(b1::Rect{N}, b2::Rect{N}) where {N}
     for i in 1:N
         maximum(b1)[i] < minimum(b2)[i] || return false
@@ -448,16 +483,23 @@ end
 
 meets(b1::Rect{N}, b2::Rect{N}) where {N} = maximum(b1) == minimum(b2)
 
-function overlaps(b1::Rect{N}, b2::Rect{N}) where {N}
-    for i in 1:N
-        maximum(b2)[i] > maximum(b1)[i] > minimum(b2)[i] &&
-            minimum(b1)[i] < minimum(b2)[i] || return false
-    end
-    return true
+"""
+    overlaps(a::Rect, b::Rect)
+
+Returns true if the given rectangles overlap, i.e. if their intersection is not
+empty.
+"""
+function overlaps(a::Rect{N}, b::Rect{N}) where {N}
+    mini1 = minimum(a)
+    maxi1 = maximum(a)
+    mini2 = minimum(b)
+    maxi2 = maximum(b)
+
+    return all(mini1 .< maxi2) && all(maxi1 .> mini2)
 end
 
 function starts(b1::Rect{N}, b2::Rect{N}) where {N}
-    return if minimum(b1) == minimum(b2)
+    if minimum(b1) == minimum(b2)
         for i in 1:N
             maximum(b1)[i] < maximum(b2)[i] || return false
         end

test/geometrytypes.jl

@@ -493,183 +493,6 @@ end
     @test centered(Circle{Float64}) == Circle(Point2(0.0), 0.5)
 end
 
-@testset "Rectangles" begin
-    rect = Rect2f(0, 7, 20, 3)
-    @test (rect + 4) == Rect2f(4, 11, 20, 3)
-    @test (rect + Vec(2, -2)) == Rect2f(2, 5, 20, 3)
-
-    @test (rect - 4) == Rect2f(-4, 3, 20, 3)
-    @test (rect - Vec(2, -2)) == Rect2f(-2, 9, 20, 3)
-
-    base = Vec3f(1, 2, 3)
-    wxyz = Vec3f(-2, 4, 2)
-    rect = Rect3f(base, wxyz)
-    @test (rect + 4) == Rect3f(base .+ 4, wxyz)
-    @test (rect + Vec(2, -2, 3)) == Rect3f(base .+ Vec(2, -2, 3), wxyz)
-
-    @test (rect - 4) == Rect3f(base .- 4, wxyz)
-    @test (rect - Vec(2, -2, 7)) == Rect3f(base .- Vec(2, -2, 7), wxyz)
-
-    rect = Rect2f(0, 7, 20, 3)
-    @test (rect * 4) == Rect2f(0, 7 * 4, 20 * 4, 3 * 4)
-    @test (rect * Vec(2, -2)) == Rect2f(0, -7 * 2, 20 * 2, -3 * 2)
-
-    base = Vec3f(1, 2, 3)
-    wxyz = Vec3f(-2, 4, 2)
-    rect = Rect3f(base, wxyz)
-    @test (rect * 4) == Rect3f(base .* 4, wxyz .* 4)
-    @test (rect * Vec(2, -2, 3)) == Rect3f(base .* Vec(2, -2, 3), wxyz .* Vec(2, -2, 3))
-
-    rect1 = Rect(Vec(0.0, 0.0), Vec(1.0, 2.0))
-    rect2 = Rect(0.0, 0.0, 1.0, 2.0)
-    @test rect1 isa GeometryBasics.HyperRectangle{2,Float64}
-    @test rect1 == rect2
-
-    split1, split2 = GeometryBasics.split(rect1, 2, 1)
-    @test widths(split1) == widths(split2)
-    @test origin(split1) == Vec(0, 0)
-    @test origin(split2) == Vec(0, 1)
-    @test in(split1, rect1) && in(split2, rect1)
-    @test !(in(rect1, split1) || in(rect1, split2))
-
-    rect1 = Rect(Vec(0.0, 0.0, -1.0), Vec(1.0, 2.0, 1.0))
-    split1, split2 = GeometryBasics.split(rect1, 1, 0.75)
-    @test widths(split1) == Vec(0.75, 2, 1)
-    @test widths(split2) == Vec(0.25, 2, 1)
-    @test origin(split1) == Vec(0, 0, -1)
-    @test origin(split2) == Vec(0.75, 0, -1)
-    @test in(split1, rect1) && in(split2, rect1)
-    @test !(in(rect1, split1) || in(rect1, split2))
-
-    prim = Rect(0.0, 0.0, 1.0, 1.0)
-    @test length(prim) == 2
-
-    @test width(prim) == 1.0
-    @test height(prim) == 1.0
-
-    b1 = Rect2(0.0, 0.0, 2.0, 2.0)
-    b2 = Rect2(0, 0, 2, 2)
-    @test isequal(b1, b2)
-
-    pt = Point(1.0, 1.0)
-    b1 = Rect(0.0, 0.0, 1.0, 1.0)
-    @test in(pt, b1)
-
-    rect = Rect(0.0, 0.0, 1.0, 1.0)
-    @test GeometryBasics.positive_widths(rect) isa GeometryBasics.HyperRectangle{2,Float64}
-
-    h1 = Rect(0.0, 0.0, 1.0, 1.0)
-    h2 = Rect(1.0, 1.0, 2.0, 2.0)
-    @test union(h1, h2) isa GeometryBasics.HyperRectangle{2,Float64}
-    # @test GeometryBasics.diff(h1, h2) == h1
-    @test GeometryBasics.intersect(h1, h2) isa GeometryBasics.HyperRectangle{2,Float64}
-
-    b = Rect(0.0, 0.0, 1.0, 1.0)
-    v = Vec(1, 2)
-    @test update(b, v) isa GeometryBasics.HyperRectangle{2,Float64}
-    v = Vec(1.0, 2.0)
-    @test update(b, v) isa GeometryBasics.HyperRectangle{2,Float64}
-
-    @testset "euclidean distances" begin
-        p = Vec(5.0, 4.0)
-        rect = Rect(0.0, 0.0, 1.0, 1.0)
-        @test min_dist_dim(rect, p, 1) == 4.0
-        @test min_dist_dim(rect, p, 2) == 3.0
-        @test max_dist_dim(rect, p, 1) == 5.0
-        @test max_dist_dim(rect, p, 2) == 4.0
-        @test minmax_dist_dim(rect, p, 1) == (4.0, 5.0)
-
-        rect1 = Rect(0.0, 0.0, 1.0, 1.0)
-        rect2 = Rect(3.0, 1.0, 4.0, 2.0)
-        @test min_dist_dim(rect1, rect2, 1) == 2.0
-        @test min_dist_dim(rect1, rect2, 2) == 0.0
-        @test max_dist_dim(rect1, rect2, 1) == 7.0
-        @test max_dist_dim(rect1, rect2, 2) == 3.0
-        @test minmax_dist_dim(rect1, rect2, 1) == (2.0, 7.0)
-
-        r = Rect2f(-1, -1, 2, 3)
-        p = Point2f(1, 2) + Point2f(3, 4)
-        @test min_euclidean(r, p) == 5f0
-        @test max_euclidean(r, p) ≈ sqrt(5*5 + 7*7)
-
-        r2 = Rect2f(0, 0, 2, 3)
-        @test min_euclidean(r, r2) == 0f0
-        @test max_euclidean(r, r2) == 5f0
-        @test minmax_euclidean(r, r2) == (0f0, 5f0)
-    end
-
-    @test !before(rect1, rect2)
-    rect1 = Rect(0.0, 0.0, 1.0, 1.0)
-    rect2 = Rect(3.0, 2.0, 4.0, 2.0)
-    @test before(rect1, rect2)
-
-    @test !meets(rect1, rect2)
-    rect2 = Rect(1.0, 1.0, 4.0, 2.0)
-    @test meets(rect1, rect2)
-
-    rect1 = Rect(1.0, 1.0, 2.0, 2.0)
-    rect2 = Rect(0.0, 0.0, 2.0, 1.0)
-    @test !overlaps(rect1, rect2)
-    rect1 = Rect(1.0, 1.0, 2.0, 2.0)
-    rect2 = Rect(1.5, 1.5, 2.0, 2.0)
-    @test overlaps(rect1, rect2)
-
-    rect1 = Rect(1.0, 1.0, 2.0, 2.0)
-    rect2 = Rect(0.0, 0.0, 2.0, 1.0)
-    @test !GeometryBasics.starts(rect1, rect2)
-    rect2 = Rect(1.0, 1.0, 1.5, 1.5)
-    @test !GeometryBasics.starts(rect1, rect2)
-    rect2 = Rect(1.0, 1.0, 3.0, 3.0)
-    @test GeometryBasics.starts(rect1, rect2)
-
-    rect1 = Rect(1.0, 1.0, 2.0, 2.0)
-    rect2 = Rect(0.0, 0.0, 4.0, 4.0)
-    @test during(rect1, rect2)
-    rect1 = Rect(0.0, 0.0, 2.0, 3.0)
-    rect2 = Rect(1.0, 1.0, 4.0, 2.0)
-    @test !during(rect1, rect2)
-
-    rect1 = Rect(1.0, 1.0, 2.0, 2.0)
-    rect2 = Rect(0.0, 0.0, 4.0, 4.0)
-    @test !finishes(rect1, rect2)
-    rect1 = Rect(1.0, 0.0, 1.0, 1.0)
-    rect2 = Rect(0.0, 0.0, 2.0, 1.0)
-    @test !finishes(rect1, rect2)
-    rect1 = Rect(1.0, 1.0, 1.0, 2.0)
-    rect2 = Rect(0.0, 0.0, 2.0, 3.0)
-    @test finishes(rect1, rect2)
-
-    rect1 = @inferred Rect(1, 2, 3, 4, 5, 6, 7, 8)
-    rect2 = Rect(Vec(1, 2, 3, 4), Vec(5, 6, 7, 8))
-    @test rect1 == rect2
-
-    @testset "Matrix Multiplications" begin
-        r = Rect2f(-1, -2, 4, 3)
-
-        # TODO: this seems quite dangerous: We pad points with ones which makes
-        #       sense for translations if we go to D+1, but is nonsense if we
-        #       go higher dimensions than that.
-        M = rand(Mat4f)
-        ps = Point2f[M * Point(p..., 1, 1) for p in coordinates(r)]
-        @test Rect2f(ps) ≈ M * r
-
-        M = Mat2f(0.5, -0.3, 0.7, 1.5)
-        ps = Point2f[M * p for p in coordinates(r)]
-        @test Rect2f(ps) ≈ M * r
-
-        r = Rect3f(-1, -2, -3, 2, 4, 1)
-        M = rand(Mat4f)
-        ps = Point3f[M * Point(p..., 1) for p in coordinates(r)]
-        @test Rect3f(ps) ≈ M * r
-    end
-
-    # TODO: this is effectively 0-indexed... should it be?
-    M = reshape(collect(11:100), 10, 9)[1:9, :]
-    r = Rect2i(2, 4, 2, 4)
-    @test M[r] == [53 63 73 83; 54 64 74 84]
-
-end
-
 @testset "LineStrings" begin
     ps1 = rand(Point2f, 10)
     ls1 = LineString(ps1)