I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Smallest enclosing circle problem

From Rosetta Code
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[edit]

Translation of: Wren
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))
}
}
Output:
{(0.500000, 0.500000), 0.707107}
{(1.000000, 1.000000), 5.000000}

Julia[edit]

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, 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()
 
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

Mathematica/Wolfram Language[edit]

f=BoundingRegion[#,"MinBall"]&;
f[{{0,0},{0,1},{1,0}}]
f[{{5,-2},{-3,-2},{-2,5},{1,6},{0,2}}]
f[{{2,4,-1},{1,5,-3},{8,-4,1},{3,9,-5}}]
f[{{-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`}}]
Output:
Ball[{0.5,0.5},0.707107]
Ball[{1.,1.},5.]
Disk[{11/2,5/2,-2},Sqrt[115/2]]
Ball[{0.0405866,0.555668,-0.231368,0.607407,-0.200346},2.7946]

Nim[edit]

Translation of: Go
Translation of: Wren
import math, random, strutils
 
type
Point = tuple[x, y: float]
Circle = tuple[c: Point; r: float]
 
 
func `$`(p: Point): string =
## Return the string representation of a point.
"($1, $2)".format(p.x, p.y)
 
 
func dist(a, b: Point): float =
## Return the distance between two points.
hypot(a.x - b.x, a.y - b.y)
 
 
func getCircleCenter(bx, by, cx, cy: float): Point =
## Return the center of a circle defined by three points.
let
b = bx * bx + by * by
c = cx * cx + cy * cy
d = bx * cy - by * cx
result = ((cy * b - by * c) / (2 * d), (bx * c - cx * b) / (2 * d))
 
func contains(ci: Circle; p: Point): bool =
## Return whether a circle contains the point 'p'.
## Used by 'in' and 'notin'.
dist(ci.c, p) <= ci.r
 
 
func contains(ci: Circle; ps: seq[Point]): bool =
## Return whether a circle contains a slice of points.
## Used by 'in'.
for p in ps:
if p notin ci: return false
result = true
 
 
func `$`(ci: Circle): string =
## Return the string representation of a circle.
"Center $1, Radius $2".format(ci.c, ci.r)
 
 
func circleFrom3(a, b, c: Point): Circle =
## Return the smallest circle that intersects three points.
var 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
result = (c: i, r: dist(i, a))
 
 
func circleFrom2(a, b: Point): Circle =
## Return the smallest circle that intersects two points.
let c = ((a.x + b.x) * 0.5, (a.y + b.y) * 0.5)
result = (c: c, r: dist(a, b) * 0.5)
 
 
func secTrivial(rs: seq[Point]): Circle =
## Return the smallest enclosing circle for n <= 3.
case rs.len
of 0: return ((0.0, 0.0), 0.0)
of 1: return (rs[0], 0.0)
of 2: return circleFrom2(rs[0], rs[1])
of 3: discard
else: raise newException(ValueError, "There shouldn't be more than three points.")
 
# Three points.
for i in 0..1:
for j in (i + 1)..2:
let c = circleFrom2(rs[i], rs[j])
if rs in c: return c
result = circleFrom3(rs[0], rs[1], rs[2])
 
 
proc welzl(ps: var seq[Point]; rs: seq[Point]; n: int): Circle =
## Helper function for Welzl method.
var rc = rs
if n == 0 or rc.len == 3: return secTrivial(rc)
let idx = rand(n-1)
let p = ps[idx]
swap ps[idx], ps[n-1]
let d = welzl(ps, rc, n-1)
if p in d: return d
rc.add p
result = welzl(ps, rc, n-1)
 
 
proc welzl(ps: seq[Point]): Circle =
# Applies the Welzl algorithm to find the SEC.
var pc = ps
pc.shuffle()
result = welzl(pc, @[], pc.len)
 
 
when isMainModule:
randomize()
const Tests = [@[(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)]]
 
for test in Tests:
echo welzl(test)
Output:
Center (0.5, 0.5), Radius 0.7071067811865476
Center (1.0, 1.0), Radius 5.0

Phix[edit]

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.

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 midpoint = sq_div(sq_add(a,b),2)
atom halfdiameter = distance(a,b)/2
circle res = { midpoint, halfdiameter }
return res
end function
 
function circle_from3(point a, b, c)
-- return a unique circle that intersects three points
point bma = sq_sub(b,a),
cma = sq_sub(c,a)
atom {{aX,aY},{bX,bY},{cX,cY}} = {a,bma,cma},
B = sum(sq_power(bma,2)),
C = sum(sq_power(cma,2)),
D = (bX*cY - bY*cX)*2
point centre = {(cY*B - bY*C)/D + aX,
(bX*C - cX*B)/D + aY }
atom radius = distance(centre,a) -- (=== b,c)
circle res = { centre, radius }
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)
circle c = welzlr(shuffle(p),{})
string s = sprintf("centre %v radius %.14g",c)
printf(1,"Points %v ==> %s\n",{p,s})
end procedure
 
constant tests = {{{0, 0},{ 0, 1},{ 1,0}},
{{5,-2},{-3,-2},{-2,5},{1,6},{0,2}}}
papply(tests,welzl)
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

n-dimensional[edit]

Translation of: Python

Uses the last test from Julia, however, since that's given more accurately.

constant inf = 1e300*1e300,
nan = -(inf/inf),
eps = 1e-8
 
function sqnorm(sequence p)
return sum(sq_power(p,2))
end function
 
function sqdist(sequence p1, p2)
return sqnorm(sq_sub(p1,p2))
end function
 
class ProjectorStack
--
-- Stack of points that are shifted / projected to put first one at origin.
--
sequence vs = {}
 
procedure push(sequence v)
vs = append(vs, v)
end procedure
 
function pop()
sequence v = vs[$]
vs = vs[1..$-1]
return v
end function
 
function mult(sequence v)
sequence s = repeat(0,length(v))
for i=1 to length(vs) do
sequence vi = vs[i]
s = sq_add(s,sq_mul(vi,sum(sq_mul(vi, v))))
end for
return s
end function
 
end class
 
class GaertnerBoundary
public:
sequence empty_center,
centers = {},
square_radii = {}
ProjectorStack projector = new()
end class
 
function push_if_stable(GaertnerBoundary bound, sequence pt)
if length(bound.centers) == 0 then
bound.square_radii = append(bound.square_radii, 0)
bound.centers = {pt}
return true
end if
ProjectorStack M = bound.projector
sequence q0 = bound.centers[1],
center = bound.centers[$],
C = sq_sub(center,q0),
Qm = sq_sub(pt,q0),
Qm_bar = M.mult(Qm),
residue = sq_sub(Qm,Qm_bar)
atom r2 = bound.square_radii[$],
e = sqdist(Qm, C) - r2,
z = 2 * sqnorm(residue),
tol = eps * max(r2, 1),
isstable = abs(z) > tol
if isstable then
sequence center_new = sq_add(center,sq_mul(e/z,residue))
atom r2new = r2 + (e * e) / (2 * z)
ProjectorStack p = bound.projector
p.push(sq_div(residue,sqrt(sqnorm(residue))))
bound.centers = append(bound.centers, center_new)
bound.square_radii = append(bound.square_radii, r2new)
end if
return isstable
end function
 
function pop(GaertnerBoundary bound)
integer n = length(bound.centers)
bound.centers = bound.centers[1..$-1]
bound.square_radii = bound.square_radii[1..$-1]
if n >= 2 then
ProjectorStack p = bound.projector
{} = p.pop()
end if
return bound
end function
 
class NSphere
public:
sequence center
atom sqradius
end class
 
function isinside(sequence pt, NSphere nsphere)
atom r2 = sqdist(pt, nsphere.center),
R2 = nsphere.sqradius
if r2=nan then return false end if
return r2<=R2 or abs(r2-R2)<eps
end function
 
function move_to_front(sequence pts, integer i)
sequence pt = pts[i]
for j=1 to i do
{pts[j], pt} = {pt, pts[j]}
end for
return pts
end function
 
function ismaxlength(GaertnerBoundary bound)
return length(bound.centers) == length(bound.empty_center) + 1
end function
 
function makeNSphere(GaertnerBoundary bound)
NSphere res
if length(bound.centers) == 0 then
res = new({bound.empty_center, 0.0})
else
res = new({bound.centers[$], bound.square_radii[$]})
end if
return res
end function
 
function _welzl(sequence pts, integer pos, GaertnerBoundary bdry)
integer support_count = 0
NSphere nsphere = makeNSphere(bdry)
if ismaxlength(bdry) then
return {nsphere, 0}
end if
for i=1 to pos do
if not isinside(pts[i], nsphere) then
bool isstable = push_if_stable(bdry, pts[i])
if isstable then
{nsphere, integer s} = _welzl(pts, i, bdry)
bdry = pop(bdry)
pts = move_to_front(pts, i)
support_count = s + 1
end if
end if
end for
return {nsphere, support_count}
end function
 
function find_max_excess(NSphere nsphere, sequence pts, integer k1)
atom err_max = -inf
integer k_max = k1 - 1
for k=k_max to length(pts) do
sequence pt = pts[k]
atom err = sqdist(pt, nsphere.center) - nsphere.sqradius
if err > err_max then
err_max = err
k_max = k + 1
end if
end for
return {err_max, k_max - 1}
end function
 
procedure welzl(sequence points, integer maxiterations=2000)
sequence pts = points
GaertnerBoundary bdry = new({repeat(nan,length(pts[1]))})
integer t = 1
{NSphere nsphere, integer s} = _welzl(pts, t, bdry)
for i=1 to maxiterations do
atom {e, k} = find_max_excess(nsphere, pts, t + 1)
if e <= eps then exit end if
sequence pt = pts[k]
{} = push_if_stable(bdry, pt)
{NSphere nsphere_new, integer s_new} = _welzl(pts, s, bdry)
bdry = pop(bdry)
pts = move_to_front(pts, k)
nsphere = nsphere_new
t = s + 1
s = s_new + 1
end for
printf(1,"For points: ") pp(points,{pp_Indent,12,pp_FltFmt,"%11.8f"})
printf(1," Center is at: %v\n", {nsphere.center})
printf(1," Radius is: %.16g\n", {sqrt(nsphere.sqradius)})
end procedure
 
constant tests = {{{0, 0}, { 0, 1}, { 1, 0}},
{{5,-2}, {-3,-2}, {-2, 5}, {1, 6}, {0, 2}},
{{2, 4, -1}, {1, 5, -3}, {8, -4, 1}, {3, 9, -5}},
{{-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.1486364315310990, 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}}}
papply(tests,welzl)
Output:
For points: {{0,0}, {0,1}, {1,0}}
    Center is at: {0.5,0.5}
    Radius is: 0.7071067811865476
For points: {{5,-2}, {-3,-2}, {-2,5}, {1,6}, {0,2}}
    Center is at: {1,1}
    Radius is: 5
For points: {{2,4,-1}, {1,5,-3}, {8,-4,1}, {3,9,-5}}
    Center is at: {5.5,2.5,-2}
    Radius is: 7.582875444051551
For points: {{-0.64009004, 2.64370326, 0.40161221, 1.89596014,-1.16240466},
             { 0.56323937,-0.01598111,-2.19372594,-0.90945866, 0.71650367},
             {-1.77043676, 0.25182837,-1.38104443,-0.59751670, 1.08964566},
             { 1.34485787,-0.18140877,-0.42887340, 1.53271973, 0.41896462},
             { 0.21323364, 0.07659442, 0.14863643, 2.33868935,-2.30004597},
             { 0.60231532, 0.30517353, 1.07326004, 1.11483880, 0.04760584},
             { 1.36555237, 0.54614164,-0.09321952,-1.00477114, 1.65329147},
             {-0.14974382,-0.53756726,-0.15845597,-0.27507203,-0.44124702}}
    Center is at: {0.0405866078,0.5556683898,-0.2313678301,0.6074066023,-0.2003463949}
    Radius is: 2.794602003699545

Python[edit]

Translation of: Julia
import numpy as np
 
class ProjectorStack:
"""
Stack of points that are shifted / projected to put first one at origin.
"""

def __init__(self, vec):
self.vs = np.array(vec)
 
def push(self, v):
if len(self.vs) == 0:
self.vs = np.array([v])
else:
self.vs = np.append(self.vs, [v], axis=0)
return self
 
def pop(self):
if len(self.vs) > 0:
ret, self.vs = self.vs[-1], self.vs[:-1]
return ret
 
def __mul__(self, v):
s = np.zeros(len(v))
for vi in self.vs:
s = s + vi * np.dot(vi, v)
return s
 
class GaertnerBoundary:
"""
GärtnerBoundary
 
See the passage regarding M_B in Section 4 of Gärtner's paper.
"""

def __init__(self, pts):
self.projector = ProjectorStack([])
self.centers, self.square_radii = np.array([]), np.array([])
self.empty_center = np.array([np.NaN for _ in pts[0]])
 
 
def push_if_stable(bound, pt):
if len(bound.centers) == 0:
bound.square_radii = np.append(bound.square_radii, 0.0)
bound.centers = np.array([pt])
return True
q0, center = bound.centers[0], bound.centers[-1]
C, r2 = center - q0, bound.square_radii[-1]
Qm, M = pt - q0, bound.projector
Qm_bar = M * Qm
residue, e = Qm - Qm_bar, sqdist(Qm, C) - r2
z, tol = 2 * sqnorm(residue), np.finfo(float).eps * max(r2, 1.0)
isstable = np.abs(z) > tol
if isstable:
center_new = center + (e / z) * residue
r2new = r2 + (e * e) / (2 * z)
bound.projector.push(residue / np.linalg.norm(residue))
bound.centers = np.append(bound.centers, np.array([center_new]), axis=0)
bound.square_radii = np.append(bound.square_radii, r2new)
return isstable
 
def pop(bound):
n = len(bound.centers)
bound.centers = bound.centers[:-1]
bound.square_radii = bound.square_radii[:-1]
if n >= 2:
bound.projector.pop()
return bound
 
 
class NSphere:
def __init__(self, c, sqr):
self.center = np.array(c)
self.sqradius = sqr
 
def isinside(pt, nsphere, atol=1e-6, rtol=0.0):
r2, R2 = sqdist(pt, nsphere.center), nsphere.sqradius
return r2 <= R2 or np.isclose(r2, R2, atol=atol**2,rtol=rtol**2)
 
def allinside(pts, nsphere, atol=1e-6, rtol=0.0):
for p in pts:
if not isinside(p, nsphere, atol, rtol):
return False
return True
 
def move_to_front(pts, i):
pt = pts[i]
for j in range(len(pts)):
pts[j], pt = pt, np.array(pts[j])
if j == i:
break
return pts
 
def dist(p1, p2):
return np.linalg.norm(p1 - p2)
 
def sqdist(p1, p2):
return sqnorm(p1 - p2)
 
def sqnorm(p):
return np.sum(np.array([x * x for x in p]))
 
def ismaxlength(bound):
len(bound.centers) == len(bound.empty_center) + 1
 
def makeNSphere(bound):
if len(bound.centers) == 0:
return NSphere(bound.empty_center, 0.0)
return NSphere(bound.centers[-1], bound.square_radii[-1])
 
def _welzl(pts, pos, bdry):
support_count, nsphere = 0, makeNSphere(bdry)
if ismaxlength(bdry):
return nsphere, 0
for i in range(pos):
if not isinside(pts[i], nsphere):
isstable = push_if_stable(bdry, pts[i])
if isstable:
nsphere, s = _welzl(pts, i, bdry)
pop(bdry)
move_to_front(pts, i)
support_count = s + 1
return nsphere, support_count
 
def find_max_excess(nsphere, pts, k1):
err_max, k_max = -np.Inf, k1 - 1
for (k, pt) in enumerate(pts[k_max:]):
err = sqdist(pt, nsphere.center) - nsphere.sqradius
if err > err_max:
err_max, k_max = err, k + k1
return err_max, k_max - 1
 
def welzl(points, maxiterations=2000):
pts, eps = np.array(points, copy=True), np.finfo(float).eps
bdry, t = GaertnerBoundary(pts), 1
nsphere, s = _welzl(pts, t, bdry)
for i in range(maxiterations):
e, k = find_max_excess(nsphere, pts, t + 1)
if e <= eps:
break
pt = pts[k]
push_if_stable(bdry, pt)
nsphere_new, s_new = _welzl(pts, s, bdry)
pop(bdry)
move_to_front(pts, k)
nsphere = nsphere_new
t, s = s + 1, s_new + 1
return nsphere
 
if __name__ == '__main__':
TESTDATA =[
np.array([[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]]),
np.array([[5.0, -2.0], [-3.0, -2.0], [-2.0, 5.0], [1.0, 6.0], [0.0, 2.0]]),
np.array([[2.0, 4.0, -1.0], [1.0, 5.0, -3.0], [8.0, -4.0, 1.0], [3.0, 9.0, -5.0]]),
np.random.normal(size=(8, 5))
]
for test in TESTDATA:
nsphere = welzl(test)
print("For points: ", test)
print(" Center is at: ", nsphere.center)
print(" Radius is: ", np.sqrt(nsphere.sqradius), "\n")
 
Output:
For points:  [[0. 0.]
 [0. 1.]
 [1. 0.]]
    Center is at:  [0.5 0.5]
    Radius is:  0.7071067811865476 

For points:  [[ 5. -2.]
 [-3. -2.]
 [-2.  5.]
 [ 1.  6.]
 [ 0.  2.]]
    Center is at:  [1. 1.]
    Radius is:  5.0 

For points:  [[ 2.  4. -1.]
 [ 1.  5. -3.]
 [ 8. -4.  1.]
 [ 3.  9. -5.]]
    Center is at:  [ 5.5  2.5 -2. ]
    Radius is:  7.582875444051551 

For points:  [[-0.28550471 -0.12787614 -0.65666165 -0.57319826 -0.09253894]
 [ 0.08947276  0.64778564 -0.54953015  0.05332253 -0.93055218]
 [ 1.26721546 -1.26410984 -0.79937865 -0.45096792 -1.23946668]
 [-0.62653779  0.56561466  0.62403237  0.23226903 -0.95820264]
 [ 0.90973949 -0.9821474   0.41400032 -1.11268937  0.19568717]
 [-0.4931165   0.94778404 -0.10534124  0.96431358  0.36504087]
 [ 1.43895269 -0.69957774 -0.31486014 -1.14980913  0.42550193]
 [ 1.4722404   1.19711551 -0.18735466  1.69208268  1.04594225]]
    Center is at:  [ 1.00038333  0.08534583 -0.22414667  0.50927606 -0.15936026]
    Radius is:  2.076492284522762 

Raku[edit]

Translation of: Go
class P { has ($.x, $.y) is rw;                # Point
method gist { "({$.x~", "~$.y})" }
}
 
class C { has (P $.c, Numeric $.r); # Circle
method gist { "Center = " ~$.c.gist ~ " and Radius = $.r\n" }
 
# tests whether a circle contains the point 'p'
method contains(P \p --> Bool) { distSq($.c, p)$.r² }
 
# tests whether a circle contains a slice of point
method encloses(@ps --> Bool) { [and] @ps.map: { $.contains($_) } }
}
 
sub infix:<> (P \a, P \b) { a.x - b.x, a.y - b.y } # note: Unicode 'minus'
 
# returns the square of the distance between two points
sub distSq (P \a, P \b) { [+] (a − b)»² }
 
sub getCenter (\bx, \by, \cx, \cy --> P) {
my (\b,\c,\d) = bx²+by², cx²+cy², bx*cy - by*cx;
P.new: x => (cy*b - by*c) / (2 * d), y => (bx*c - cx*b) / (2 * d)
} # returns the center of a circle defined by 3 points
 
sub circleFrom3 (P \a, P \b, P \c --> C) {
my \k = $ = getCenter |(b − a), |(c − a);
k.x, k.y Z[+=] a.x, a.y;
C.new: c => k, r => distSq(k, a).sqrt
} # returns a unique circle that intersects 3 points
 
sub circleFrom2 (P \a, P \b --> C ) {
my \center = P.new: x => ((a.x + b.x) / 2), y => ((a.y + b.y) / 2) ;
C.new: c => center, r => (distSq(a, b).sqrt / 2)
} # returns smallest circle that intersects 2 points
 
sub secTrivial( @rs --> C ) {
given @rs {
when * == 0 { return C.new: c => (P.new: x => 0, y => 0), r => 0 }
when * == 1 { return C.new: c => @rs[0], r => 0 }
when * == 2 { return circleFrom2 |@rs }
when * == 3 { #`{ no-op } }
when * > 3 { die "There shouldn't be more than 3 points." }
}
for 0, 1 X 1, 2 -> ( \i, \j ) {
return $_ if .encloses(@rs) given circleFrom2 |@rs[i,j]
}
circleFrom3 |@rs
} # returns smallest enclosing circle for n ≤ 3
 
sub Welzl-helper( @ps is copy, @rs is copy , \n --> C ) {
return secTrivial(@rs) if n == 0 or @rs == 3;
my \p = @ps.shift;
return $_ if .contains(p) given Welzl-helper @ps, @rs, n-1;
Welzl-helper @ps, @rs.append(p), n-1
} # helper function for Welzl method
 
# applies the Welzl algorithm to find the SEC
sub welzl(@ps --> C) { Welzl-helper @ps.pick(*), [], @ps }
 
my @tests = (
[ (0,0), (0,1), (1,0) ],
[ (5,-2), (-3,-2), (-2,5), (1,6), (0,2) ]
).map: {
@_.map: { P.new: x => $_[0], y => $_[1] }
}
 
say "Solution for smallest circle enclosing {$_.gist} :\n", welzl $_ for @tests;
Output:
Solution for smallest circle enclosing ((0, 0) (0, 1) (1, 0)) :
Center = (0.5, 0.5) and Radius = 0.7071067811865476

Solution for smallest circle enclosing ((5, -2) (-3, -2) (-2, 5) (1, 6) (0, 2)) :
Center = (1, 1) and Radius = 5

Wren[edit]

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.

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))
Output:
Center (0.5, 0.5), Radius 0.70710678118655
Center (1, 1), Radius 5