Smallest enclosing circle problem: Difference between revisions
Line 40: | Line 40: | ||
P = eltype(pts) |
P = eltype(pts) |
||
F = eltype(P) |
F = eltype(P) |
||
projector = ProjectorStack(P[]) |
projector, centers, square_radii = ProjectorStack(P[]), P[], F[] |
||
centers = P[] |
|||
square_radii = F[] |
|||
empty_center = F(NaN) * first(pts) |
empty_center = F(NaN) * first(pts) |
||
GärtnerBoundary(centers, square_radii, projector, empty_center) |
GärtnerBoundary(centers, square_radii, projector, empty_center) |
||
end |
end |
||
⚫ | |||
⚫ | |||
⚫ | |||
function push_if_stable!(b::GärtnerBoundary, pt) |
function push_if_stable!(b::GärtnerBoundary, pt) |
||
Line 58: | Line 52: | ||
return true |
return true |
||
end |
end |
||
⚫ | |||
⚫ | |||
⚫ | |||
Qm, M = pt - q0, b.projector |
|||
C = center - q0 |
|||
r2 = b.square_radii[end] |
|||
Qm = pt - q0 |
|||
M = b.projector |
|||
Qm_bar = M*Qm |
Qm_bar = M*Qm |
||
residue = Qm - Qm_bar |
residue, e = Qm - Qm_bar, sqdist(Qm, C) - r2 |
||
⚫ | |||
e = sqdist(Qm, C) - r2 |
|||
z = 2*sqnorm(residue) |
|||
⚫ | |||
isstable = abs(z) > tol |
isstable = abs(z) > tol |
||
if isstable |
if isstable |
||
Line 91: | Line 76: | ||
pop!(b.projector) |
pop!(b.projector) |
||
end |
end |
||
b |
return b |
||
end |
end |
||
Line 100: | Line 85: | ||
function isinside(pt, nsphere::NSphere; atol=0, rtol=0) |
function isinside(pt, nsphere::NSphere; atol=0, rtol=0) |
||
r2 = sqdist(pt, center(nsphere)) |
r2, R2 = sqdist(pt, center(nsphere)), sqradius(nsphere) |
||
⚫ | |||
R2 = sqradius(nsphere) |
|||
⚫ | |||
end |
|||
function allinside(pts, nsphere; kw...) |
|||
for pt in pts |
|||
isinside(pt, nsphere; kw...) || return false |
|||
end |
|||
true |
|||
end |
end |
||
allinside(pts, nsphere; kw...) = all(pt -> isinside(pt, nsphere; kw...), pts) |
|||
function move_to_front!(pts, i) |
function move_to_front!(pts, i) |
||
Line 121: | Line 99: | ||
end |
end |
||
⚫ | |||
⚫ | |||
⚫ | |||
center(b::NSphere) = b.center |
center(b::NSphere) = b.center |
||
radius(b::NSphere) = sqrt(b.sqradius) |
radius(b::NSphere) = sqrt(b.sqradius) |
||
Line 127: | Line 108: | ||
sqdist(p1::AbstractVector, p2::AbstractVector) = sqnorm(p1-p2) |
sqdist(p1::AbstractVector, p2::AbstractVector) = sqnorm(p1-p2) |
||
sqnorm(p) = sum(abs2,p) |
sqnorm(p) = sum(abs2,p) |
||
ismaxlength(b::GärtnerBoundary) = length(b) == length(b.empty_center) + 1 |
|||
function NSphere(b::GärtnerBoundary) |
function NSphere(b::GärtnerBoundary) |
||
return isempty(b) ? NSphere(b.empty_center, 0.0) : |
|||
NSphere(b.centers[end], b.square_radii[end]) |
|||
r2 = zero(eltype(c)) |
|||
else |
|||
c = b.centers[end] |
|||
⚫ | |||
end |
|||
NSphere(c,r2) |
|||
end |
|||
function ismaxlength(b::GärtnerBoundary) |
|||
dim = length(b.empty_center) |
|||
length(b) == dim + 1 |
|||
end |
end |
||
function welzl!(pts, bdry) |
function welzl!(pts, bdry) |
||
bdry_len = length(bdry) |
bdry_len, support_count, nsphere = length(bdry), 0, NSphere(bdry) |
||
support_count = 0 |
|||
nsphere = NSphere(bdry) |
|||
ismaxlength(bdry) && return nsphere, 0 |
ismaxlength(bdry) && return nsphere, 0 |
||
for i in eachindex(pts) |
for i in eachindex(pts) |
||
Line 167: | Line 136: | ||
function find_max_excess(nsphere, pts, k1) |
function find_max_excess(nsphere, pts, k1) |
||
T = eltype(first(pts)) |
T = eltype(first(pts)) |
||
e_max = T(-Inf) |
e_max, k_max = T(-Inf), k1 -1 |
||
k_max = k1 -1 |
|||
for k in k1:length(pts) |
for k in k1:length(pts) |
||
pt = pts[k] |
pt = pts[k] |
||
Line 182: | Line 150: | ||
function welzl(points, maxiterations=2000) |
function welzl(points, maxiterations=2000) |
||
pts = deepcopy(points) |
pts = deepcopy(points) |
||
bdry = GärtnerBoundary(pts) |
bdry, t = GärtnerBoundary(pts), 1 |
||
t = 1 |
|||
nsphere, s = welzl!(prefix(pts, t), bdry) |
nsphere, s = welzl!(prefix(pts, t), bdry) |
||
for i in 1:maxiterations |
for i in 1:maxiterations |
||
Line 189: | Line 156: | ||
P = eltype(pts) |
P = eltype(pts) |
||
F = eltype(P) |
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 |
|||
s = s_new + 1 |
|||
else |
|||
return nsphere |
|||
end |
|||
end |
end |
||
return nsphere |
return nsphere |
||
end |
end |
||
function testwelzl() |
function testwelzl() |
||
testdata =[ |
testdata =[ |
Revision as of 05:53, 4 November 2020
The smallest enclosing circle problem (aka minimum covering circle problem, bounding circle problem) is a mathematical problem of computing the smallest circle that contains all of a given set of points in the Euclidean plane.
Initially it was proposed by the English mathematician James Joseph Sylvester in 1857.
- Task
Find circle of smallest radius containing all of given points.
- Circle is defined by it's center and radius;
- Points are defined by their coordinates in n-dimensional space;
- Circle (sphere) contains point when distance between point and circle center <= circle radius.
Julia
N-dimensional solution using the Welzl algorithm. Derived from the BoundingSphere.jl module at https://github.com/JuliaFEM. See also Bernd Gärtner's paper at people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf. <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]], [randn(5) for _ in 1:8] ] for test in testdata nsphere = welzl(test) println("For points: ", test) println(" Center is at: ", nsphere.center) println(" Radius is: ", radius(nsphere), "\n") end
end
testwelzl()
</lang>
- Output:
For points: [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]] Center is at: [0.5, 0.5] Radius is: 0.7071067811865476 For points: [[5.0, -2.0], [-3.0, -2.0], [-2.0, 5.0], [1.0, 6.0], [0.0, 2.0]] Center is at: [1.0, 1.0] Radius is: 5.0 For points: [[-0.6400900432782123, 2.643703255134232, 0.4016122094706093, 1.8959601399652273, -1.1624046608380516], [0.5632393652621324, -0.015981105190064373, -2.193725940351997, -0.9094586577358262, 0.7165036664470906], [-1.7704367632976061, 0.2518283698686299, -1.3810444289625348, -0.597516704360172, 1.089645656962647], [1.3448578652803103, -0.18140877132223002, -0.4288734015080915, 1.53271973321691, 0.41896461833399573], [0.2132336397213029, 0.07659442168765788, 0.148636431531099, 2.3386893481333795, -2.3000459709300927], [0.6023153188617328, 0.3051735340025314, 1.0732600437151525, 1.1148388039984898, 0.047605838564167786], [1.3655523661720959, 0.5461416420929995, -0.09321951900362889, -1.004771137760985, 1.6532914656050117], [-0.14974382165751837, -0.5375672589202939, -0.15845596754003466, -0.2750720340454088, -0.441247015836271]] Center is at: [0.0405866077655439, 0.5556683897981481, -0.2313678300507679, 0.6074066023194586, -0.2003463948612026] Radius is: 2.7946020036995454