Smallest enclosing circle problem

Smallest enclosing circle problem is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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

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