Smallest enclosing circle problem
The smallest enclosing circle problem (aka minimum covering circle problem, bounding circle problem) is a mathematical problem of computing the smallest circle that contains all of a given set of points in the Euclidean plane.
Initially it was proposed by the English mathematician James Joseph Sylvester in 1857.
- Task
Find circle of smallest radius containing all of given points.
- Circle is defined by it's center and radius;
- Points are defined by their coordinates in n-dimensional space;
- Circle (sphere) contains point when distance between point and circle center <= circle radius.
Go
<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
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 Wezl'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