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.

Go

Translation of: Wren

<lang go>package main

import (

   "fmt"
   "log"
   "math"
   "math/rand"
   "time"

)

type Point struct{ x, y float64 }

type Circle struct {

   c Point
   r float64

}

// string representation of a Point func (p Point) String() string { return fmt.Sprintf("(%f, %f)", p.x, p.y) }

// returns the square of the distance between two points func distSq(a, b Point) float64 {

   return (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y)

}

// returns the center of a circle defined by 3 points func getCircleCenter(bx, by, cx, cy float64) Point {

   b := bx*bx + by*by
   c := cx*cx + cy*cy
   d := bx*cy - by*cx
   return Point{(cy*b - by*c) / (2 * d), (bx*c - cx*b) / (2 * d)}

}

// returns whether a circle contains the point 'p' func (ci Circle) contains(p Point) bool {

   return distSq(ci.c, p) <= ci.r*ci.r

}

// returns whether a circle contains a slice of points func (ci Circle) encloses(ps []Point) bool {

   for _, p := range ps {
       if !ci.contains(p) {
           return false
       }
   }
   return true

}

// string representation of a Circle func (ci Circle) String() string { return fmt.Sprintf("{%v, %f}", ci.c, ci.r) }

// returns a unique circle that intersects 3 points func circleFrom3(a, b, c Point) Circle {

   i := getCircleCenter(b.x-a.x, b.y-a.y, c.x-a.x, c.y-a.y)
   i.x += a.x
   i.y += a.y
   return Circle{i, math.Sqrt(distSq(i, a))}

}

// returns smallest circle that intersects 2 points func circleFrom2(a, b Point) Circle {

   c := Point{(a.x + b.x) / 2, (a.y + b.y) / 2}
   return Circle{c, math.Sqrt(distSq(a, b)) / 2}

}

// returns smallest enclosing circle for n <= 3 func secTrivial(rs []Point) Circle {

   size := len(rs)
   if size > 3 {
       log.Fatal("There shouldn't be more than 3 points.")
   }
   if size == 0 {
       return Circle{Point{0, 0}, 0}
   }
   if size == 1 {
       return Circle{rs[0], 0}
   }
   if size == 2 {
       return circleFrom2(rs[0], rs[1])
   }
   for i := 0; i < 2; i++ {
       for j := i + 1; j < 3; j++ {
           c := circleFrom2(rs[i], rs[j])
           if c.encloses(rs) {
               return c
           }
       }
   }
   return circleFrom3(rs[0], rs[1], rs[2])

}

// helper function for Welzl method func welzlHelper(ps, rs []Point, n int) Circle {

   rc := make([]Point, len(rs))
   copy(rc, rs) // 'rs' passed by value so need to copy
   if n == 0 || len(rc) == 3 {
       return secTrivial(rc)
   }
   idx := rand.Intn(n)
   p := ps[idx]
   ps[idx], ps[n-1] = ps[n-1], p
   d := welzlHelper(ps, rc, n-1)
   if d.contains(p) {
       return d
   }
   rc = append(rc, p)
   return welzlHelper(ps, rc, n-1)

}

// applies the Welzl algorithm to find the SEC func welzl(ps []Point) Circle {

   var pc = make([]Point, len(ps))
   copy(pc, ps)
   rand.Shuffle(len(pc), func(i, j int) {
       pc[i], pc[j] = pc[j], pc[i]
   })
   return welzlHelper(pc, []Point{}, len(pc))

}

func main() {

   rand.Seed(time.Now().UnixNano())
   tests := [][]Point{
       {Point{0, 0}, Point{0, 1}, Point{1, 0}},
       {Point{5, -2}, Point{-3, -2}, Point{-2, 5}, Point{1, 6}, Point{0, 2}},
   }
   for _, test := range tests {
       fmt.Println(welzl(test))
   }

}</lang>

Output:

{(0.500000, 0.500000), 0.707107}
{(1.000000, 1.000000), 5.000000}

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, centers, square_radii = ProjectorStack(P[]), P[], F[]
   empty_center = F(NaN) * first(pts)
   GärtnerBoundary(centers, square_radii, projector, empty_center)

end

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, center = first(b.centers), b.centers[end]
   C, r2  = center - q0, b.square_radii[end]
   Qm, M = pt - q0, b.projector
   Qm_bar = M*Qm
   residue, e = Qm - Qm_bar, sqdist(Qm, C) - r2
   z, tol = 2*sqnorm(residue), 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
   return b

end

struct NSphere{P,F}

   center::P
   sqradius::F

end

function isinside(pt, nsphere::NSphere; atol=0, rtol=0)

   r2, R2 = sqdist(pt, center(nsphere)), sqradius(nsphere)
   return r2 <= R2 || isapprox(r2, R2;atol=atol^2,rtol=rtol^2)

end allinside(pts, nsphere; kw...) = all(pt -> isinside(pt, nsphere; kw...), pts)

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

prefix(pts, i) = view(pts, 1:i) Base.length(b::GärtnerBoundary) = length(b.centers) Base.isempty(b::GärtnerBoundary) = length(b) == 0 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) ismaxlength(b::GärtnerBoundary) = length(b) == length(b.empty_center) + 1

function NSphere(b::GärtnerBoundary)

   return isempty(b) ? NSphere(b.empty_center, 0.0) :
       NSphere(b.centers[end], b.square_radii[end])

end

function welzl!(pts, bdry)

   bdry_len, support_count, nsphere = length(bdry), 0, 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, k_max = T(-Inf), 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, t = GärtnerBoundary(pts), 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)
       e <= eps(F) && break
       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 = s + 1, s_new + 1
   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]],
       [[2.0, 4.0, -1.0], [1.0, 5.0, -3.0], [8.0, -4.0, 1.0], [3.0, 9.0, -5.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: [[2.0, 4.0, -1.0], [1.0, 5.0, -3.0], [8.0, -4.0, 1.0], [3.0, 9.0, -5.0]]
    Center is at: [5.5, 2.5, -2.0]
    Radius is: 7.582875444051551

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

Phix

Based on the same code as Wren, and likewise limited to 2D circles - I believe (but cannot prove) the main barrier to coping with more than two dimensions lies wholly within the circle_from3() routine. <lang Phix>type point(sequence) return true end type enum CENTRE, RADIUS type circle(sequence) return true end type

function distance(point a, b)

   return sqrt(sum(sq_power(sq_sub(a,b),2)))

end function

function in_circle(point p, circle c)

   return distance(p,c[CENTRE]) <= c[RADIUS]

end function

function circle_from2(point a, b)

   -- return the smallest circle that intersects 2 points:
   point centre = sq_div(sq_add(a,b),2) -- (ie midpoint)
   atom radius = distance(a,b)/2
   circle res = { centre, radius }
   return res

end function

function circle_from3(point a, b, c)

   -- return a unique circle that intersects three points 
   atom {{aX,aY},{bX,bY},{cX,cY}} = {a,sq_sub(b,a),sq_sub(c,a)},
        B = bX * bX + bY * bY,
        C = cX * cX + cY * cY,
        D = bX * cY - bY * cX 
   point centre = { (cY * B - bY * C) / (2 * D) + aX, 
                    (bX * C - cX * B) / (2 * D) + aY }
   circle res = { centre, distance(centre, a) }
   return res

end function

function trivial(sequence r)

   integer l = length(r)
   switch l do
       case 0: return {{0,0},0}
       case 1: return {r[1],0}
       case 2: return circle_from2(r[1],r[2])
       case 3: return circle_from3(r[1],r[2],r[3])
   end switch

end function

function welzlr(sequence p, r)

   if p={} or length(r)=3 then return trivial(r) end if
   point p1 = p[1]
   p = p[2..$]
   circle d = welzlr(p, r)
   if in_circle(p1,d) then return d end if
   return welzlr(p, append(r,p1))

end function

procedure welzl(sequence p)

   printf(1,"Points %v ==> centre %v radius %.14g\n",prepend(welzlr(shuffle(p),{}),p))

end procedure

constant tests = {{ { 0, 0 }, { 0, 1 }, { 1, 0 } },

                 { { 5, -2 }, { -3, -2 }, { -2, 5 }, { 1, 6 }, { 0, 2 } }}

papply(tests,welzl)</lang>

Output:

Points {{0,0},{0,1},{1,0}} ==> centre {0.5,0.5} radius 0.70710678118655
Points {{5,-2},{-3,-2},{-2,5},{1,6},{0,2}} ==> centre {1,1} radius 5

Wren

Well a circle is a two dimensional figure and so, despite any contradictory indications in the task description, that's what this solution provides.

It is based on Welzl's algorithm and follows closely the C++ code here. <lang ecmascript>import "random" for Random

var Rand = Random.new()

class Point {

   // returns the square of the distance between two points
   static distSq(a, b) { (a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y) }

   // returns the center of a circle defined by 3 points
   static getCircleCenter(bx, by, cx, cy) {
       var b = bx*bx + by*by
       var c = cx*cx + cy*cy
       var d = bx*cy - by*cx
       return Point.new((cy*b - by*c) / (2 * d), (bx*c - cx*b) / (2 * d))
   }

   // constructs a new Point object
   construct new(x, y) {
       _x = x
       _y = y
   }

   // basic properties
   x { _x }
   x=(o) { _x = o }
   y { _y }
   y=(o) { _y = o }

   // returns a copy of the current instance
   copy() { Point.new(_x, _y) }

   // string representation
   toString { "(%(_x), %(_y))" }

}

class Circle {

   // returns a unique circle that intersects 3 points
   static from(a, b, c) {
       var i = Point.getCircleCenter(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y)
       i.x = i.x + a.x
       i.y = i.y + a.y
       return Circle.new(i, Point.distSq(i, a).sqrt)
   }

   // returns smallest circle that intersects 2 points
   static from(a, b) {
       var c = Point.new((a.x + b.x) / 2, (a.y + b.y) / 2)
       return Circle.new(c, Point.distSq(a, b).sqrt/2)
   }

   // returns smallest enclosing circle for n <= 3
   static secTrivial(rs) {
       var size = rs.count
       if (size > 3) Fiber.abort("There shouldn't be more than 3 points.")
       if (size == 0) return Circle.new(Point.new(0, 0), 0)
       if (size == 1) return Circle.new(rs[0], 0)
       if (size == 2) return Circle.from(rs[0], rs[1])
       for (i in 0..1) {
           for (j in i+1..2) {
               var c = Circle.from(rs[i], rs[j])
               if (c.encloses(rs)) return c
           }
       }
       return Circle.from(rs[0], rs[1], rs[2])
   }

   // private helper method for Welzl method
   static welzl_(ps, rs, n) {
       rs = rs.toList // passed by value so need to copy
       if (n == 0 || rs.count == 3) return secTrivial(rs)
       var idx = Rand.int(n)
       var p = ps[idx]
       ps[idx] = ps[n-1]
       ps[n-1] = p
       var d = welzl_(ps, rs, n-1)
       if (d.contains(p)) return d
       rs.add(p)
       return welzl_(ps, rs, n-1)
   }

   // applies the Welzl algorithm to find the SEC
   static welzl(ps) {
       var pc = List.filled(ps.count, null)
       for (i in 0...pc.count) pc[i] = ps[i].copy()
       Rand.shuffle(pc)
       return welzl_(pc, [], pc.count)
   }

   // constructs a new Circle object from its center point and its radius
   construct new(c, r) {
       _c = c.copy()
       _r = r
   }

   // basic properties
   c { _c }
   r { _r }

   // returns whether the current instance contains the point 'p'
   contains(p) { Point.distSq(_c, p) <= _r * _r }

   // returns whether the current instance contains a list of points
   encloses(ps) {
       for (p in ps) if (!contains(p)) return false
       return true
   }

   // string representation
   toString { "Center %(_c), Radius %(_r)" }

}

var tests = [

   [Point.new(0, 0), Point.new(0, 1), Point.new(1, 0)],
   [Point.new(5, -2), Point.new(-3, -2), Point.new(-2, 5), Point.new(1, 6), Point.new(0, 2)]

]

for (test in tests) System.print(Circle.welzl(test))</lang>

Output:

Center (0.5, 0.5), Radius 0.70710678118655
Center (1, 1), Radius 5