Smallest enclosing circle problem
The smallest enclosing circle problem (aka minimum covering circle problem, bounding circle problem) is a mathematical problem of computing the smallest circle that contains all of a given set of points in the Euclidean plane.
Initially it was proposed by the English mathematician James Joseph Sylvester in 1857.
- Task
Find circle of smallest radius containing all of given points.
- Circle is defined by it's center and radius;
- Points are defined by their coordinates in n-dimensional space;
- Circle (sphere) contains point when distance between point and circle center <= circle radius.
Julia
N-dimensional solution using the Welzl algorithm. Derived from the BoundingSphere.jl module at https://github.com/JuliaFEM See also Bernd Gärtner's paper at people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf <lang julia>import Base.pop!, Base.push!, Base.length, Base.* using LinearAlgebra, Random
struct ProjectorStack{P <: AbstractVector}
vs::Vector{P}
end Base.push!(p::ProjectorStack, v) = (push!(p.vs, v); p) Base.pop!(p::ProjectorStack) = (pop!(p.vs); p) Base.:*(p::ProjectorStack, v) = (s = zero(v); for vi in p.vs s += vi * dot(vi, v) end; s)
"""
GärtnerBoundary
See the passage regarding M_B in Section 4 of Gärtner's paper. """ mutable struct GärtnerBoundary{P<:AbstractVector, F<:AbstractFloat}
centers::Vector{P} square_radii::Vector{F} # projection onto of affine space spanned by points # shifted such that first point becomes origin projector::ProjectorStack{P} empty_center::P # center of nsphere spanned by empty boundary
end
function GärtnerBoundary(pts)
P = eltype(pts) F = eltype(P) projector = ProjectorStack(P[]) centers = P[] square_radii = F[] empty_center = F(NaN) * first(pts) GärtnerBoundary(centers, square_radii, projector, empty_center)
end
prefix(pts, i) = view(pts, 1:i) Base.length(b::GärtnerBoundary) = length(b.centers) Base.isempty(b::GärtnerBoundary) = length(b) == 0
function push_if_stable!(b::GärtnerBoundary, pt)
if isempty(b) push!(b.square_radii, zero(eltype(pt))) push!(b.centers, pt) dim = length(pt) return true end
q0 = first(b.centers) center = b.centers[end] C = center - q0 r2 = b.square_radii[end]
Qm = pt - q0 M = b.projector Qm_bar = M*Qm residue = Qm - Qm_bar
e = sqdist(Qm, C) - r2 z = 2*sqnorm(residue)
tol = eps(eltype(pt)) * max(r2, one(r2)) isstable = abs(z) > tol if isstable center_new = center + (e/z) * residue r2new = r2 + (e^2)/(2z) push!(b.projector, residue / norm(residue)) push!(b.centers, center_new) push!(b.square_radii, r2new) end return isstable
end
function Base.pop!(b::GärtnerBoundary)
n = length(b) pop!(b.centers) pop!(b.square_radii) if n >= 2 pop!(b.projector) end b
end
struct NSphere{P,F}
center::P sqradius::F
end
function isinside(pt, nsphere::NSphere; atol=0, rtol=0)
r2 = sqdist(pt, center(nsphere)) R2 = sqradius(nsphere) r2 <= R2 || isapprox(r2, R2;atol=atol^2,rtol=rtol^2)
end
function allinside(pts, nsphere; kw...)
for pt in pts isinside(pt, nsphere; kw...) || return false end true
end
function move_to_front!(pts, i)
pt = pts[i] for j in eachindex(pts) pts[j], pt = pt, pts[j] j == i && break end pts
end
center(b::NSphere) = b.center radius(b::NSphere) = sqrt(b.sqradius) sqradius(b::NSphere) = b.sqradius dist(p1,p2) = norm(p1-p2) sqdist(p1::AbstractVector, p2::AbstractVector) = sqnorm(p1-p2) sqnorm(p) = sum(abs2,p)
function NSphere(b::GärtnerBoundary)
if isempty(b) c = b.empty_center r2 = zero(eltype(c)) else c = b.centers[end] r2 = b.square_radii[end] end NSphere(c,r2)
end
function ismaxlength(b::GärtnerBoundary)
dim = length(b.empty_center) length(b) == dim + 1
end
function welzl!(pts, bdry)
bdry_len = length(bdry) support_count = 0 nsphere = NSphere(bdry) ismaxlength(bdry) && return nsphere, 0 for i in eachindex(pts) pt = pts[i] if !isinside(pt, nsphere) pts_i = prefix(pts, i-1) isstable = push_if_stable!(bdry, pt) if isstable nsphere, s = welzl!(pts_i, bdry) pop!(bdry) move_to_front!(pts, i) support_count = s + 1 end end end return nsphere, support_count
end
function find_max_excess(nsphere, pts, k1)
T = eltype(first(pts)) e_max = T(-Inf) k_max = k1 -1 for k in k1:length(pts) pt = pts[k] e = sqdist(pt, center(nsphere)) - sqradius(nsphere) if e > e_max e_max = e k_max = k end end return e_max, k_max
end
function welzl(points, maxiterations=2000)
pts = deepcopy(points) bdry = GärtnerBoundary(pts) t = 1 nsphere, s = welzl!(prefix(pts, t), bdry) for i in 1:maxiterations e, k = find_max_excess(nsphere, pts, t + 1) P = eltype(pts) F = eltype(P) if e > eps(F) pt = pts[k] push_if_stable!(bdry, pt) nsphere_new, s_new = welzl!(prefix(pts, s), bdry) pop!(bdry) move_to_front!(pts, k) nsphere = nsphere_new t = s + 1 s = s_new + 1 else return nsphere end end return nsphere
end
function testwelzl()
testdata =[ [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]], [[5.0, -2.0], [-3.0, -2.0], [-2.0, 5.0], [1.0, 6.0], [0.0, 2.0]], [randn(5) for _ in 1:8] ] for test in testdata nsphere = welzl(test) println("For points: ", test) println(" Center is at: ", nsphere.center) println(" Radius is: ", radius(nsphere), "\n") end
end
testwelzl()
</lang>
- Output:
For points: [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]] Center is at: [0.5, 0.5] Radius is: 0.7071067811865476 For points: [[5.0, -2.0], [-3.0, -2.0], [-2.0, 5.0], [1.0, 6.0], [0.0, 2.0]] Center is at: [1.0, 1.0] Radius is: 5.0 For points: [[-0.6400900432782123, 2.643703255134232, 0.4016122094706093, 1.8959601399652273, -1.1624046608380516], [0.5632393652621324, -0.015981105190064373, -2.193725940351997, -0.9094586577358262, 0.7165036664470906], [-1.7704367632976061, 0.2518283698686299, -1.3810444289625348, -0.597516704360172, 1.089645656962647], [1.3448578652803103, -0.18140877132223002, -0.4288734015080915, 1.53271973321691, 0.41896461833399573], [0.2132336397213029, 0.07659442168765788, 0.148636431531099, 2.3386893481333795, -2.3000459709300927], [0.6023153188617328, 0.3051735340025314, 1.0732600437151525, 1.1148388039984898, 0.047605838564167786], [1.3655523661720959, 0.5461416420929995, -0.09321951900362889, -1.004771137760985, 1.6532914656050117], [-0.14974382165751837, -0.5375672589202939, -0.15845596754003466, -0.2750720340454088, -0.441247015836271]] Center is at: [0.0405866077655439, 0.5556683897981481, -0.2313678300507679, 0.6074066023194586, -0.2003463948612026] Radius is: 2.7946020036995454