Add minimal-area rotated hyperrectangle overapproximation (#3795)

* use DocumenterCitations * remove 'Any' * v3.1.0 * Format .jl files * collect references in docstrings * Add minimal-area rotated hyperrectangle overapproximation and tests * Update test/Approximations/overapproximate.jl --------- Co-authored-by: schillic <[email protected]>
JuliaReach · Jan 27, 2025 · 25a485a · 25a485a
1 parent 1a38d7f
commit 25a485a
Show file tree

Hide file tree

Showing 2 changed files with 67 additions and 0 deletions.
diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl
@@ -603,6 +603,63 @@ function overapproximate(QM::QuadraticMap{N,<:SparsePolynomialZonotope},
     return SparsePolynomialZonotope(cz, Ḡ[:, H], Gz, Ē[1:p, H], indexvector(PZ))
 end
 
+"""
+    overapproximate(P::VPolygon, ::Type{<:LinearMap{N,<:Hyperrectangle}}) where {N}
+
+Overapproximate polygon with minimal-area rotated hyperrectangle.
+
+### Input
+
+- `P`               -- polygon
+- `LinearMap{N,<:Hyperrectangle}` -- target type (a linear map of a hyperrectangle)
+
+### Output 
+
+A LinearMap of a Hyperrectangle that represents the minimal area rotated bounding rectangle
+of the polygon P.
+
+### Algorithm
+
+This method uses an approach described in [Finding Minimum Area Rectangle for Given Points](https://gis.stackexchange.com/questions/22895/finding-minimum-area-rectangle-for-given-points/22934)
+
+"""
+
+function overapproximate(P::VPolygon, ::Type{<:LinearMap{N,<:Hyperrectangle}}) where {N}
+    min_area = N(Inf)
+    vert_P = P.vertices
+    n = length(vert_P)
+
+    center = Vector{N}(undef, 2)
+    radius = Vector{N}(undef, 2)
+    R = Matrix{N}(undef, 2, 2) 
+
+    @inbounds for i in eachindex(vert_P)
+        a = vert_P[i]
+        next_idx = i % n + 1
+        b = vert_P[next_idx]
+        e = b - a
+        v′ = normalize(e)
+        w′ = [-v′[2], v′[1]]
+
+        min_x, max_x = extrema(dot(vertex, v′) for vertex in vert_P)
+        min_y, max_y = extrema(dot(vertex, w′) for vertex in vert_P)
+
+        current_area = (max_x - min_x) * (max_y - min_y)
+
+        if current_area < min_area
+            min_area = current_area
+            center .= ((max_x + min_x) / 2, (max_y + min_y) / 2)
+            radius .= ((max_x - min_x) / 2, (max_y - min_y) / 2)
+            R[:, 1] = v′
+            R[:, 2] = w′
+        end
+    end
+
+    min_rectangle = Hyperrectangle(center, radius)
+
+    return LinearMap(R, min_rectangle)
+end
+
 # function to be loaded by Requires
 function load_paving_overapproximation()
     return quote

diff --git a/test/Approximations/overapproximate.jl b/test/Approximations/overapproximate.jl
@@ -137,6 +137,16 @@ for N in [Float64, Rational{Int}, Float32]
     @test overapproximate(Rectification(Z), Zonotope) == convert(Zonotope, Singleton(N[0, 0]))
     Z = Zonotope(N[-2, 2], N[1 1; 1 0])
     @test overapproximate(Rectification(Z), Zonotope) == Zonotope(N[0, 2], hcat(N[0, 1]))
+
+    # overapproximate polygon with minimal-area rotated hyperrectangle
+    P = VPolygon([N[1, 0], N[0, 1], N[-1, 0], N[0, -1]])
+    R = overapproximate(P, LinearMap{N, Hyperrectangle{N}})
+    @test set(R).center == N[0, 0]
+    @test set(R).radius == N[1 / √2, 1 / √2] 
+    @test R.M == N[-1/sqrt(2) -1/sqrt(2); 1/sqrt(2) -1/sqrt(2)] 
+    if N <: AbstractFloat
+        @test isequivalent(P, R)
+    end
 end
 
 # tests that do not work with Rational{Int}