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.


11l

Translation of: Nim
Translation of: Go
T Point = (Float, Float)

T Circle
   Point c
   Float r

   F (c, r)
      .c = c
      .r = r

   F String()
      R ‘Center ’(.c)‘, Radius ’(.r)

F getCircleCenter(Float bx, by, cx, cy)
   V b = bx * bx + by * by
   V c = cx * cx + cy * cy
   V d = bx * cy - by * cx
   R Point((cy * b - by * c) / (2 * d), (bx * c - cx * b) / (2 * d))

F contains(Circle ci, Point p)
   R sqlen(ci.c - p) <= ci.r ^ 2

F encloses(Circle ci, [Point] ps)
   L(p) ps
      I !contains(ci, p)
         R 0B
   R 1B

F circleFrom3(Point a, b, c)
   V 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
   R Circle(i, distance(i, a))

F circleFrom2(Point a, b)
   V c = Point((a.x + b.x) * 0.5, (a.y + b.y) * 0.5)
   R Circle(c, distance(a, b) / 2)

F secTrivial([Point] rs)
   I rs.empty
      R Circle(Point(0, 0), 0)
   I rs.len == 1
      R Circle(rs[0], 0)
   I rs.len == 2
      R circleFrom2(rs[0], rs[1])
   assert(rs.len == 3, ‘There shouldn't be more than three points.’)

   L(i) 2
      L(j) i + 1 .< 3
         V c = circleFrom2(rs[i], rs[j])
         I encloses(c, rs)
            R c
   R circleFrom3(rs[0], rs[1], rs[2])

F welzlHelper([Point] &ps, [Point] rs, Int n)
   V rc = copy(rs)
   I n == 0 | rc.len == 3
      R secTrivial(rc)
   V idx = random:(n)
   V p = ps[idx]
   swap(&ps[idx], &ps[n - 1])
   V d = welzlHelper(&ps, rc, n - 1)
   I contains(d, p)
      R d
   rc.append(p)
   R welzlHelper(&ps, rc, n - 1)

F welzl([Point] ps)
   V pc = copy(ps)
   random:shuffle(&pc)
   R welzlHelper(&pc, [Point](), pc.len)

V Tests = [[Point(0.0, 0.0), Point(0.0, 1.0), Point(1.0, 0.0)],
           [Point(5.0, -2.0), Point(-3.0, -2.0), Point(-2.0, 5.0), Point(1.0, 6.0), Point(0.0, 2.0)]]

L(test) Tests
   print(welzl(test))
Output:
Center (0.5, 0.5), Radius 0.707106781
Center (1, 1), Radius 5

F#

This task uses Centre and radius of a circle passing through 3 points in a plane (F#)

// Smallest enclosing circle problem. Nigel Galloway: February 22nd., 2024
let mcc points=
  let fN g=points|>List.map(fun n->(n,d g n))|>List.maxBy snd
  let P1,P2=Seq.allPairs points points|>Seq.maxBy(fun(n,g)->d n g)
  let c1,r1=let ((nx,ny) as n),((gx,gy) as g)=P1,P2 in (((nx+gx)/2.0,(ny+gy)/2.0),(d n g)/2.0)
  match fN c1 with (_,d) when d=r1->(c1,r1)
                  |(P3,_)->let c2,r2=circ P1 P2 P3
                           match fN c2 with (_,d) when d=r2->(c2,r2)
                                           |(P4,_)->[circ P1 P3 P4;circ P2 P3 P4]|>List.minBy snd

let testMCC n=let centre,radius=mcc n in printfn "Minimum Covering Circle is centred at %A with a radius of %f" centre radius
testMCC [(0.0, 0.0);(0.0, 1.0);(1.0, 0.0)]
testMCC [(5.0, -2.0);(-3.0, -2.0);(-2.0, 5.0);(1.0, 6.0);(0.0, 2.0)]
testMCC [(15.85,6.05);(29.82,22.11);(4.84,20.32);(25.47,29.46);(33.65,17.31);(7.30,10.43);(28.15,6.67);(20.85,11.47);(22.83,2.07);(26.28,23.15);(14.39,30.24);(30.24,25.23)]
testMCC [(15.85,6.05);(29.82,22.11);(4.84,20.32);(25.47,29.46);(33.65,17.31);(7.30,10.43);(28.15,6.67);(20.85,11.47);(22.83,2.08);(26.28,23.15);(14.39,30.24);(30.24,25.23)]
Output:
Minimum Covering Circle is centred at (0.5, 0.5) with a radius of 0.707107
Minimum Covering Circle is centred at (1.0, 1.0) with a radius of 5.000000
Minimum Covering Circle is centred at (18.97851566, 16.2654108) with a radius of 14.708624
Minimum Covering Circle is centred at (18.97952365, 16.27401209) with a radius of 14.707010

Go

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}

jq

Works with: jq

Works with gojq, the Go implementation of jq

Adapted from Wren

In the following, a Point (x,y) is represented by a JSON array [x,y], and a Circle by a JSON object of the form {center, radius}.

# Vector sum of the input array of Points
def sum: transpose | map(add);

# The square of the distance between two points
def distSq:
  def sq: . * .;
  . as [$a, $b]
  | (($a[0] - $b[0]) | sq) +(($a[1] - $b[1]) | sq) ;

# The distance between two points
def distance:
  distSq | sqrt;

# Does the input Circle contain the point $p?
def contains($p):
  ([.center, $p] | distSq) <= .radius * .radius;

# Does the input Circle contain the given array of points?
def encloses($ps):
  . as $c
  | all($ps[]; . as $p | $c | contains($p));

# The center of a circle defined by 3 points
def getCircleCenter(bx; by; cx; cy):
  (bx*bx + by*by) as $b
  | (cx*cx + cy*cy) as $c
  | (bx*cy - by*cx) as $d
  | [(cy*$b - by*$c) / (2 * $d), (bx*$c - cx*$b) / (2 * $d)];

# The (smallest) circle that intersects 1, 2 or 3 distinct points
def Circle:
  if length == 1 then {center: .[0], radius: 0}
  elif length == 2
  then {center: (sum | map(./2)),
        radius: (distance/2) } 

  elif length == 3
  then . as [$a, $b, $c]
  | ([getCircleCenter($b[0] - $a[0];
                      $b[1] - $a[1];
                      $c[0] - $a[0];
                      $c[1] - $a[1]), $a] | sum) as $i
  | {center: $i, radius: ([$i, $a]|distance)}
  else "Circle/0 expects an input array of from 1 to 3 Points" | error
  end;
  
# The smallest enclosing circle for n <= 3
def secTrivial:
  . as $rs
  | length as $length
  | if $length > 3 then "Internal error: secTrivial given more than 3 points." | error
    elif $length == 0 then [[0, 0]] | Circle
    elif $length <= 2 then Circle
    else first(
       range(0; 2) as $i
       | range($i+1; 3) as $j
       | ([$rs[$i], $rs[$j]] | Circle)
       | select(encloses($rs)) )
    // Circle
    end;

# Use the Welzl algorithm to find the SEC of the incoming array of points
def welzl:
  # Helper function:
  #   $rs is an array of points on the boundary;
  #   $n keeps track of how many points remain:
  def welzl_($rs; $n):
    if $n == 0 or ($rs|length) == 3
    then $rs | secTrivial
    else 0 as $idx     # or Rand($n)
    | .[$idx] as $p
    | .[$idx] = .[$n-1]
    | .[$n-1] = $p
    | welzl_($rs; $n-1) as $d
    | if ($d | contains($p)) then $d
      else welzl_($rs + [$p]; $n-1)
      end
    end ;
  
  welzl_([]; length);

def tests:
    [ [0, 0], [0, 1], [1, 0] ],
    [ [5, -2], [-3, -2], [-2, 5], [1, 6], [0, 2] ]
;

tests
| welzl
| "Circle(\(.center), \(.radius))"
Output:
Circle([0.5,0.5], 0.7071067811865476)
Circle([1,1], 5)

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

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

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

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.

with javascript_semantics
type point(sequence s) return true end type
enum CENTRE, RADIUS
type circle(sequence s) 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(deep_copy(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

Translation of: Python

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

with javascript_semantics
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
 
function projector_mult(sequence vs, 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
 
-- Gaertner Boundary:
sequence empty_center,
         centers,
         square_radii,
         -- Stack of points that are shifted / projected to put first one at origin:
         projector

function push_if_stable(sequence pt)
    if length(centers) == 0 then
        square_radii = append(square_radii, 0)
        centers = {pt}
        return true
    end if
    sequence q0 = centers[1],
             center = centers[$],
             C = sq_sub(center,q0),
             Qm = sq_sub(pt,q0),
             Qm_bar = projector_mult(projector,Qm),
             residue = sq_sub(Qm,Qm_bar)
    atom r2 = 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)
        projector = append(projector,sq_div(residue,sqrt(sqnorm(residue))))
        centers = append(centers, center_new)
        square_radii = append(square_radii, r2new)
    end if
    return isstable
end function
 
procedure pop_gb()
    integer n = length(centers)
    centers = centers[1..$-1]
    square_radii = square_radii[1..$-1]
    if n >= 2 then
        projector = projector[1..$-1] -- (pop/discard)
    end if
end procedure
 
function isinside(sequence pt, center, atom sqradius)
    atom r2 = sqdist(pt, center)
    return r2!=nan and (r2<=sqradius or abs(r2-sqradius)<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()
    return length(centers) == length(empty_center) + 1
end function
 
function _welzl(sequence pts, integer pos)
    integer support_count = 0, s
    sequence center = iff(length(centers)?centers[$]:empty_center)
    atom sqradius = iff(length(centers)?square_radii[$]:0)
    if ismaxlength() then
        return {center, sqradius, 0}
    end if
    for i=1 to pos do
        if not isinside(pts[i], center, sqradius) then
            bool isstable = push_if_stable(pts[i])
            if isstable then
                {center, sqradius, s} = _welzl(pts, i)
                pop_gb()
                pts = move_to_front(deep_copy(pts), i)
                support_count = s + 1
            end if
        end if
    end for
    return {center, sqradius, support_count}
end function
 
function find_max_excess(sequence center, pts, atom sqradius, 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, center) - 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
    empty_center = repeat(nan,length(pts[1]))
    centers = {}
    square_radii = {}
    projector = {}

    integer t = 1, s_new
    {sequence center, atom sqradius, integer s} = _welzl(pts, t)
    for i=1 to maxiterations do
        atom {e, k} = find_max_excess(center, pts, sqradius, t + 1)
        if e <= eps then exit end if
        sequence pt = pts[k]
        {} = push_if_stable(pt)
        {center, sqradius, s_new} = _welzl(pts, s)
        pop_gb()
        pts = move_to_front(deep_copy(pts), k)
        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", {center})
    printf(1,"    Radius is: %.16g\n", {sqrt(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

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

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 |(ba), |(ca);
   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

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