From 6d6a0f097b4e0472198402ca5753f45d365f3f0c Mon Sep 17 00:00:00 2001
From: ffreyer <frederic481994@hotmail.de>
Date: Sat, 31 Jan 2026 15:56:36 +0100
Subject: [PATCH 01/10] add rect diff

---
 src/GeometryBasics.jl        |  4 ++--
 src/primitives/rectangles.jl | 37 ++++++++++++++++++++++++-------
 test/geometrytypes.jl        | 43 +++++++++++++++++++++++++++++++-----
 3 files changed, 69 insertions(+), 15 deletions(-)

diff --git a/src/GeometryBasics.jl b/src/GeometryBasics.jl
index d01865ae..57276380 100644
--- a/src/GeometryBasics.jl
+++ b/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
 
diff --git a/src/primitives/rectangles.jl b/src/primitives/rectangles.jl
index 912161a4..e9ee6184 100644
--- a/src/primitives/rectangles.jl
+++ b/src/primitives/rectangles.jl
@@ -337,14 +337,35 @@ function Base.union(h1::Rect{N}, h2::Rect{N}) where {N}
     return Rect{N}(m, mm - m)
 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)
-
-# Perform a difference between two Rects.
-# """
-# diff(h1::Rect, h2::Rect) = h1
+# 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
+
 
 """
     intersect(h1::Rect, h2::Rect)
diff --git a/test/geometrytypes.jl b/test/geometrytypes.jl
index bab22981..ff2b4a77 100644
--- a/test/geometrytypes.jl
+++ b/test/geometrytypes.jl
@@ -558,11 +558,44 @@ end
     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}
+    @testset "union, intersect, diff" begin
+        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}
+
+        # Cases:
+        # 1. always got 1+ negative?
+        # 2. mixed, always 0 trues
+        # 3. always positive?
+
+        a = Rectf(-1, -1, -1, 2, 2, 2)
+        # no intersection -> return a
+        @test bbox_diff(a, Rectf(2, -1, -1, 2, 2, 2)) == a
+        @test bbox_diff(a, Rectf(2, 2, 2, 3, 3, 3)) == a
+        @test bbox_diff(a, Rectf(-3, -3, -3, 1, 1, 1)) == a
+        @test bbox_diff(a, Rectf(-8, -4, -4, 4, 10, 8)) == a
+
+        # intersection is a chunk < a in every dimension -> return a
+        @test bbox_diff(a, Rectf(0, 0, 0, 2, 2, 2)) == a
+        @test bbox_diff(a, Rectf(-2, -4, -1, 2, 3.8, 1)) == a
+        @test bbox_diff(a, Rectf(-0.5, -0.5, -1, 1, 1, 2)) == a
+
+        # intersection is matching the size of a in just one dimension -> return a
+        @test bbox_diff(a, Rectf(-2, 0, 0, 4, 2, 2)) == a
+        @test bbox_diff(a, Rectf(-0.5, -0.5, -2, 1, 1, 10)) == a
+
+        # intersection is a slab matching the size of a in two dimension -> reduced
+        @test bbox_diff(a, Rectf(-2, -2, 0, 4, 4, 2)) == Rectf(-1, -1, -1, 2, 2, 1)
+        @test bbox_diff(a, Rectf(-5, -2, -2, 5.2, 4, 4)) ≈ Rectf(0.2, -1, -1, 0.8, 2, 2)
+        # edge case: bisection
+        @test bbox_diff(a, Rectf(-0.5, -2, -2, 1, 4, 4)) == a
+
+        # intersection is a -> empty Rect
+        @test bbox_diff(a, Rectf(-2, -2, -2, 4, 4, 4)) == Rect3f()
+        @test bbox_diff(a, Rectf(-1, -1, -1, 2, 2, 2)) == Rect3f()
+    end
 
     b = Rect(0.0, 0.0, 1.0, 1.0)
     v = Vec(1, 2)

From 400ea6895f821dbb9f765be8fa4f49b9f9610bd9 Mon Sep 17 00:00:00 2001
From: ffreyer <frederic481994@hotmail.de>
Date: Sat, 31 Jan 2026 23:31:09 +0100
Subject: [PATCH 02/10] move and extend Rect tests

---
 test/geometrytypes.jl | 210 --------------------------------
 test/rectangles.jl    | 271 ++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl      |   4 +
 3 files changed, 275 insertions(+), 210 deletions(-)
 create mode 100644 test/rectangles.jl

diff --git a/test/geometrytypes.jl b/test/geometrytypes.jl
index ff2b4a77..c7c59519 100644
--- a/test/geometrytypes.jl
+++ b/test/geometrytypes.jl
@@ -493,216 +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}
-
-    @testset "union, intersect, diff" begin
-        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}
-
-        # Cases:
-        # 1. always got 1+ negative?
-        # 2. mixed, always 0 trues
-        # 3. always positive?
-
-        a = Rectf(-1, -1, -1, 2, 2, 2)
-        # no intersection -> return a
-        @test bbox_diff(a, Rectf(2, -1, -1, 2, 2, 2)) == a
-        @test bbox_diff(a, Rectf(2, 2, 2, 3, 3, 3)) == a
-        @test bbox_diff(a, Rectf(-3, -3, -3, 1, 1, 1)) == a
-        @test bbox_diff(a, Rectf(-8, -4, -4, 4, 10, 8)) == a
-
-        # intersection is a chunk < a in every dimension -> return a
-        @test bbox_diff(a, Rectf(0, 0, 0, 2, 2, 2)) == a
-        @test bbox_diff(a, Rectf(-2, -4, -1, 2, 3.8, 1)) == a
-        @test bbox_diff(a, Rectf(-0.5, -0.5, -1, 1, 1, 2)) == a
-
-        # intersection is matching the size of a in just one dimension -> return a
-        @test bbox_diff(a, Rectf(-2, 0, 0, 4, 2, 2)) == a
-        @test bbox_diff(a, Rectf(-0.5, -0.5, -2, 1, 1, 10)) == a
-
-        # intersection is a slab matching the size of a in two dimension -> reduced
-        @test bbox_diff(a, Rectf(-2, -2, 0, 4, 4, 2)) == Rectf(-1, -1, -1, 2, 2, 1)
-        @test bbox_diff(a, Rectf(-5, -2, -2, 5.2, 4, 4)) ≈ Rectf(0.2, -1, -1, 0.8, 2, 2)
-        # edge case: bisection
-        @test bbox_diff(a, Rectf(-0.5, -2, -2, 1, 4, 4)) == a
-
-        # intersection is a -> empty Rect
-        @test bbox_diff(a, Rectf(-2, -2, -2, 4, 4, 4)) == Rect3f()
-        @test bbox_diff(a, Rectf(-1, -1, -1, 2, 2, 2)) == Rect3f()
-    end
-
-    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)
diff --git a/test/rectangles.jl b/test/rectangles.jl
new file mode 100644
index 00000000..01fe9fab
--- /dev/null
+++ b/test/rectangles.jl
@@ -0,0 +1,271 @@
+@testset "Arithmetics" 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))
+end
+
+@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
+
+@testset "Constructors & Getters" begin
+    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
+
+    prim = Rect(0.0, 0.0, 1.0, 1.0)
+    @test length(prim) == 2
+
+    @test width(prim) == 1.0
+    @test height(prim) == 1.0
+
+    rect = Rect(0.0, 0.0, 1.0, 1.0)
+    @test GeometryBasics.positive_widths(rect) == rect
+    @test GeometryBasics.positive_widths(Rect2f(0, 0, -1, -2)) == Rect2f(-1, -2, 1, 2)
+
+    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
+end
+
+@testset "checks" begin
+    @testset "isempty" begin
+        @test isempty(Rect2f())
+        @test isempty(Rect3f())
+        @test isempty(Rect2i())
+        @test isempty(Rect(0, 0, 0, 0))
+        @test isempty(Rect(0, 0, 0, 0, 0, 0))
+        @test isempty(Rect(0, 0, 0, 0, 0, 0))
+    end
+
+    @testset "equality" begin
+        b1 = Rect2(0.0, 0.0, 2.0, 2.0)
+        b2 = Rect2(0, 0, 2, 2)
+        @test isequal(b1, b2)
+        @test !isequal(b1, Rect2f(0.1, 0, 2, 2))
+    end
+
+    @testset "in" begin
+        pt = Point(1.0, 1.0)
+        b1 = Rect(0.0, 0.0, 1.0, 1.0)
+        @test in(pt, b1)
+        @test in(Rect2f(0, 0, 1, 1), Rect2f(-1, -1, 3, 3))
+        @test in(Rect2f(0, 0, 1, 1), Rect2f(0, 0, 3, 3))
+        @test !in(Rect2f(0, 0, 1, 1), Rect2f(0.1, 0, 3, 3))
+    end
+
+    @testset "overlaps" begin
+        @test !overlaps(Rect2f(-1, -1, 2, 2), Rect2f(-2, -2, 1, 1))
+        @test !overlaps(Rect2f(-1, -1, 2, 2), Rect2f(-2, -2, 1, 1.2))
+        @test overlaps(Rect2f(-1, -1, 2, 2), Rect2f(-2, -2, 1.2, 1.2))
+        @test overlaps(Rect2f(-1, -1, 2, 2), Rect2f(-2, -2, 4, 4))
+        @test overlaps(Rect2f(-1, -1, 2, 2), Rect2f(-0.5, -0.5, 1, 1))
+        @test overlaps(Rect2f(-1, -1, 2, 2), Rect2f(0.5, 0.5, 1, 1))
+        @test !overlaps(Rect2f(-1, -1, 2, 2), Rect2f(1, 1, 1, 1))
+        @test !overlaps(Rect2f(-1, -1, 2, 2), Rect2f(2, 2, 1, 1))
+
+        @test !overlaps(Rect2f(-2, -0.1, 0.5, 0.2), Rect2f(-1, -1, 2, 2))
+        @test !overlaps(Rect2f(-2, -0.1, 1, 0.2), Rect2f(-1, -1, 2, 2))
+        @test overlaps(Rect2f(-2, -0.1, 1.1, 0.2), Rect2f(-1, -1, 2, 2))
+        @test overlaps(Rect2f(-2, -0.1, 3, 0.2), Rect2f(-1, -1, 2, 2))
+        @test overlaps(Rect2f(-1, -0.1, 3, 0.2), Rect2f(-1, -1, 2, 2))
+        @test overlaps(Rect2f(0, -0.1, 3, 0.2), Rect2f(-1, -1, 2, 2))
+        @test !overlaps(Rect2f(1, -0.1, 3, 0.2), Rect2f(-1, -1, 2, 2))
+        @test !overlaps(Rect2f(2, -0.1, 3, 0.2), Rect2f(-1, -1, 2, 2))
+    end
+end
+
+@testset "Set Operations" begin
+    @testset "split" begin
+        rect1 = Rect(Vec(0.0, 0.0), Vec(1.0, 2.0))
+        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))
+    end
+
+    @testset "union" begin
+        h1 = Rect(0.0, 0.0, 1.0, 1.0)
+        h2 = Rect(1.0, 1.0, 2.0, 2.0)
+        @test union(h1, h2) == Rect2d(0, 0, 3, 3)
+        @test union(Rect2f(), Rect2f(0, 0, 1, 1)) == Rect2f(0, 0, 1, 1)
+        @test union(Rect2f(0, 0, 1, 1), Rect2f()) == Rect2f(0, 0, 1, 1)
+        @test union(Rect3f(0,0,0,1,1,1), Rect3f(-2,-2,-2,1,1,1)) == Rect3f(-2,-2,-2,3,3,3)
+    end
+
+    @testset "intersection" begin
+        h1 = Rect(0.0, 0.0, 1.0, 1.0)
+        h2 = Rect(1.0, 1.0, 2.0, 2.0)
+        @test intersect(h1, h2) == Rect()
+        @test intersect(Rect2f(0,0,2,1), Rect2f(1,0,1,2)) == Rect2f(1, 0, 1, 1)
+        @test intersect(Rect2f(0,0,2,2), Rect2f(1,1,2,2)) == Rect2f(1, 1, 1, 1)
+        @test intersect(Rect2f(0,0,2,3), Rect2f(1,1,2,1)) == Rect2f(1, 1, 1, 1)
+        @test intersect(Rect2f(0,0,3,3), Rect2f(1,1,1,1)) == Rect2f(1, 1, 1, 1)
+    end
+
+    @testset "diff" begin
+        a = Rectf(-1, -1, -1, 2, 2, 2)
+        # no intersection -> return a
+        @test bbox_diff(a, Rectf(2, -1, -1, 2, 2, 2)) == a
+        @test bbox_diff(a, Rectf(2, 2, 2, 3, 3, 3)) == a
+        @test bbox_diff(a, Rectf(-3, -3, -3, 1, 1, 1)) == a
+        @test bbox_diff(a, Rectf(-8, -4, -4, 4, 10, 8)) == a
+
+        # intersection is a chunk < a in every dimension -> return a
+        @test bbox_diff(a, Rectf(0, 0, 0, 2, 2, 2)) == a
+        @test bbox_diff(a, Rectf(-2, -4, -1, 2, 3.8, 1)) == a
+        @test bbox_diff(a, Rectf(-0.5, -0.5, -1, 1, 1, 2)) == a
+
+        # intersection is matching the size of a in just one dimension -> return a
+        @test bbox_diff(a, Rectf(-2, 0, 0, 4, 2, 2)) == a
+        @test bbox_diff(a, Rectf(-0.5, -0.5, -2, 1, 1, 10)) == a
+
+        # intersection is a slab matching the size of a in two dimension -> reduced
+        @test bbox_diff(a, Rectf(-2, -2, 0, 4, 4, 2)) == Rectf(-1, -1, -1, 2, 2, 1)
+        @test bbox_diff(a, Rectf(-5, -2, -2, 5.2, 4, 4)) ≈ Rectf(0.2, -1, -1, 0.8, 2, 2)
+        # edge case: bisection
+        @test bbox_diff(a, Rectf(-0.5, -2, -2, 1, 4, 4)) == a
+
+        # intersection is a -> empty Rect
+        @test bbox_diff(a, Rectf(-2, -2, -2, 4, 4, 4)) == Rect3f()
+        @test bbox_diff(a, Rectf(-1, -1, -1, 2, 2, 2)) == Rect3f()
+    end
+end
+
+@testset "update" begin
+    b = Rect(0.0, 0.0, 1.0, 1.0)
+    @test update(b, Vec(1, 2)) == Rect2d(0, 0, 1, 2)
+    @test update(b, Vec(2.0, 2.0)) == Rect2d(0, 0, 2, 2)
+    @test update(b, Vec(-1, 2)) == Rect2d(-1, 0, 2, 2)
+end
+
+@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
+
+# Does this make sense to extend to Rects?
+@testset "Allen's interval algebra" begin
+    rect1 = Rect(0.0, 0.0, 1.0, 1.0)
+    rect2 = Rect(3.0, 1.0, 4.0, 2.0)
+    @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)
+end
+
+@testset "Utilities" begin
+    # 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
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 5fb98ded..917fba4e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -344,6 +344,10 @@ end
     include("geometrytypes.jl")
 end
 
+@testset "Rectangles and Bounding Boxes" begin
+    include("rectangles.jl")
+end
+
 @testset "Point & Vec type" begin
     include("fixed_arrays.jl")
 end

From 3fc4445878f9b8b5fbde763e0f356ec0579fb768 Mon Sep 17 00:00:00 2001
From: ffreyer <frederic481994@hotmail.de>
Date: Sat, 31 Jan 2026 23:31:56 +0100
Subject: [PATCH 03/10] consider 0 width empty too

---
 src/primitives/rectangles.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/primitives/rectangles.jl b/src/primitives/rectangles.jl
index e9ee6184..9ff4f3df 100644
--- a/src/primitives/rectangles.jl
+++ b/src/primitives/rectangles.jl
@@ -324,7 +324,7 @@ end
 
 Return `true` if any of the widths of `h` are negative.
 """
-Base.isempty(h::Rect{N,T}) where {N,T} = any(<(zero(T)), h.widths)
+Base.isempty(h::Rect{N,T}) where {N,T} = any(<=(zero(T)), h.widths)
 
 """
     union(r1::Rect{N}, r2::Rect{N})

From 0336572fccf13415a5bf3a6ac56719d55b70f228 Mon Sep 17 00:00:00 2001
From: ffreyer <frederic481994@hotmail.de>
Date: Sat, 31 Jan 2026 23:33:38 +0100
Subject: [PATCH 04/10] make overlaps() test overlap instead of duplicating
 `in`

---
 src/primitives/rectangles.jl | 22 ++++++++++++++++------
 1 file changed, 16 insertions(+), 6 deletions(-)

diff --git a/src/primitives/rectangles.jl b/src/primitives/rectangles.jl
index 9ff4f3df..bb12b975 100644
--- a/src/primitives/rectangles.jl
+++ b/src/primitives/rectangles.jl
@@ -469,12 +469,22 @@ 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(mini2 .<= mini1 .< maxi2) || # minimum(a) in b
+        all(mini2 .< maxi1 .<= maxi2) || # maximum(a) in b
+        all(mini1 .<= mini2 .< maxi1) || # minimum(b) in a
+        all(mini1 .< maxi2 .<= maxi1) # maximum(b) in a
 end
 
 function starts(b1::Rect{N}, b2::Rect{N}) where {N}

From 08c1dbc21b66521956506527110c0a6e7d350b25 Mon Sep 17 00:00:00 2001
From: ffreyer <frederic481994@hotmail.de>
Date: Sat, 31 Jan 2026 23:35:15 +0100
Subject: [PATCH 05/10] check for empty Rects in union

---
 src/primitives/rectangles.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/primitives/rectangles.jl b/src/primitives/rectangles.jl
index bb12b975..c0138409 100644
--- a/src/primitives/rectangles.jl
+++ b/src/primitives/rectangles.jl
@@ -332,9 +332,11 @@ Base.isempty(h::Rect{N,T}) where {N,T} = any(<=(zero(T)), h.widths)
 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) && return h2
+    isempty(h2) && return h1
+    mini = min.(minimum(h1), minimum(h2))
+    maxi = max.(maximum(h1), maximum(h2))
+    return Rect{N}(mini, maxi - mini)
 end
 
 # TODO: Add a diff that returns the slabs created. This could be anywhere between

From 2e7b91017955a9cca2faf16d88cb388de17f21f5 Mon Sep 17 00:00:00 2001
From: ffreyer <frederic481994@hotmail.de>
Date: Sat, 31 Jan 2026 23:36:31 +0100
Subject: [PATCH 06/10] return empty Rect in `intersect` when inputs don't
 intersect

---
 src/primitives/rectangles.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/primitives/rectangles.jl b/src/primitives/rectangles.jl
index c0138409..7b364e58 100644
--- a/src/primitives/rectangles.jl
+++ b/src/primitives/rectangles.jl
@@ -374,10 +374,12 @@ end
 
 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}

From 5b7d2379fd22f76bf829a39090b80dc0098cb1a1 Mon Sep 17 00:00:00 2001
From: ffreyer <frederic481994@hotmail.de>
Date: Sat, 31 Jan 2026 23:36:45 +0100
Subject: [PATCH 07/10] some cleanup

---
 src/primitives/rectangles.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/primitives/rectangles.jl b/src/primitives/rectangles.jl
index 7b364e58..40e59b40 100644
--- a/src/primitives/rectangles.jl
+++ b/src/primitives/rectangles.jl
@@ -389,12 +389,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
@@ -464,6 +464,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
@@ -492,7 +494,7 @@ function overlaps(a::Rect{N}, b::Rect{N}) where {N}
 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

From 63cd4e0c14d04913f27ab73983f26c67a2da25ea Mon Sep 17 00:00:00 2001
From: ffreyer <frederic481994@hotmail.de>
Date: Sat, 31 Jan 2026 23:38:42 +0100
Subject: [PATCH 08/10] switch to macOS 15

---
 .github/workflows/ci.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 16ade19d..57165685 100644
--- a/.github/workflows/ci.yml
+++ b/.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:

From 74be2c2b88672fb78fdf227f3849422151484e94 Mon Sep 17 00:00:00 2001
From: ffreyer <frederic481994@hotmail.de>
Date: Sun, 1 Feb 2026 16:27:41 +0100
Subject: [PATCH 09/10] fix overlap

---
 src/primitives/rectangles.jl | 5 +----
 test/rectangles.jl           | 2 ++
 2 files changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/primitives/rectangles.jl b/src/primitives/rectangles.jl
index 40e59b40..da5905de 100644
--- a/src/primitives/rectangles.jl
+++ b/src/primitives/rectangles.jl
@@ -487,10 +487,7 @@ function overlaps(a::Rect{N}, b::Rect{N}) where {N}
     mini2 = minimum(b)
     maxi2 = maximum(b)
 
-    return all(mini2 .<= mini1 .< maxi2) || # minimum(a) in b
-        all(mini2 .< maxi1 .<= maxi2) || # maximum(a) in b
-        all(mini1 .<= mini2 .< maxi1) || # minimum(b) in a
-        all(mini1 .< maxi2 .<= maxi1) # maximum(b) in a
+    return all(mini1 .< maxi2) && all(maxi1 .> mini2)
 end
 
 function starts(b1::Rect{N}, b2::Rect{N}) where {N}
diff --git a/test/rectangles.jl b/test/rectangles.jl
index 01fe9fab..1d4f5e83 100644
--- a/test/rectangles.jl
+++ b/test/rectangles.jl
@@ -111,6 +111,8 @@ end
         @test overlaps(Rect2f(0, -0.1, 3, 0.2), Rect2f(-1, -1, 2, 2))
         @test !overlaps(Rect2f(1, -0.1, 3, 0.2), Rect2f(-1, -1, 2, 2))
         @test !overlaps(Rect2f(2, -0.1, 3, 0.2), Rect2f(-1, -1, 2, 2))
+
+        @test overlaps(Rect3f(-0.3, -0.3, -0.5, 0.6, 0.6, 1.0), Rect3f(-0.2, -0.2, -0.6, 0.4, 0.4, 1.2))
     end
 end
 

From 61d29fa307ce3960fb41ee4bcc14752b0a3ffe5a Mon Sep 17 00:00:00 2001
From: ffreyer <frederic481994@hotmail.de>
Date: Fri, 6 Mar 2026 17:15:12 +0100
Subject: [PATCH 10/10] consider all vs any widths zero in isempty

---
 src/primitives/rectangles.jl | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

diff --git a/src/primitives/rectangles.jl b/src/primitives/rectangles.jl
index da5905de..f5a0105b 100644
--- a/src/primitives/rectangles.jl
+++ b/src/primitives/rectangles.jl
@@ -320,11 +320,19 @@ 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})
@@ -332,8 +340,8 @@ Base.isempty(h::Rect{N,T}) where {N,T} = any(<=(zero(T)), h.widths)
 Returns a new `Rect{N}` which contains both r1 and r2.
 """
 function Base.union(h1::Rect{N}, h2::Rect{N}) where {N}
-    isempty(h1) && return h2
-    isempty(h2) && return h1
+    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)