diff --git a/src/ConcreteOperations/isdisjoint.jl b/src/ConcreteOperations/isdisjoint.jl index 7d122783ac..91c4ae602c 100644 --- a/src/ConcreteOperations/isdisjoint.jl +++ b/src/ConcreteOperations/isdisjoint.jl @@ -613,6 +613,53 @@ end end end +""" +# Extended help + + isdisjoint(H::Hyperrectangle, E::Ellipsoid) + +### Algorithm + +The sets are disjoint if the ellipse center lies outside the Minkowski sum +(rectangle expanded by the ellipsoid). Otherwise, we check one corner using +a quadratic form derived from the ellipsoids support function. + +### Notes +It works only for 2D axis-aligned rectangles and ellipsoids. + +### Reference +David Eberly, “Distance Between a Point and an Ellipse, an Ellipsoid, or a Hyperellipsoid”, +Geometric Tools, 2015. https://www.geometrictools.com/Documentation/DistancePointEllipseEllipsoid.pdf +""" + +@commutative function isdisjoint(H::Hyperrectangle, E::Ellipsoid) + @assert dim(H) == dim(E) == 2 "$H and $E must both have 2 dimensions." + + # center to the origin + H_trans = translate(H, -H.center) + E_trans = translate(E, -H.center) + K = E_trans.center + + bbox = overapproximate(H_trans ⊕ E_trans, Hyperrectangle) + if any(abs.(K) .> bbox.radius) + return true + end + + # find the rectangle corner corresponding to K. + s = sign.(K) + P = s .* H_trans.radius + Δ = K - P + + # support vector on the boundary of E_trans in direction Δ. + v = σ(Δ, E_trans) + if any(s .* v .<= 0) + return false + end + + # check the ellipse condition + return (ρ(Δ, E_trans))^2 ≤ 1 +end + # ============== # # disambiguation # # ============== # diff --git a/test/ConcreteOperations/isdisjoint.jl b/test/ConcreteOperations/isdisjoint.jl index 8bf3aee2a2..6db8b80428 100644 --- a/test/ConcreteOperations/isdisjoint.jl +++ b/test/ConcreteOperations/isdisjoint.jl @@ -25,3 +25,25 @@ for N in [Float64, Float32, Rational{Int}] end end end + +for N in [Float64, Float32] + # case 1: ellipse completely outside + R1 = Hyperrectangle(N[1.0, -1.0], N[2, 1]) + El1 = Ellipsoid(N[5.0, 5.0], Matrix{N}(I, 2, 2)) + @test isdisjoint(R1, El1) + + # case 2: ellipse completely inside + R2 = Hyperrectangle(N[0.0, 0.0], N[4.0, 1.0]) + El2 = Ellipsoid(N[1.0, 0.0], Matrix{N}(I, 2, 2)) + @test !isdisjoint(R2, El2) + + # case 3: ellipse overlapping + R3 = Hyperrectangle(N[0.0, 0.0], N[1.0, 1.0]) + El3 = Ellipsoid(N[2.0, 0.0], Diagonal(N[2.0, 1.0])) + @test !isdisjoint(R3, El3) + + # case 4: ellipse tangent in one point + R4 = Hyperrectangle(N[0.0, 0.0], N[1.0, 1.0]) + El4 = Ellipsoid(N[2.0, 0.0], Diagonal(N[1.0, 1.0])) + @test !isdisjoint(R4, El4) +end \ No newline at end of file