Smallest enclosing circle problem: Difference between revisions
julia example
No edit summary |
(julia example) |
||
Line 10:
* Points are defined by their coordinates in n-dimensional space;
* Circle (sphere) contains point when ''distance between point and circle center <= circle radius''.
=={{header|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 = ProjectorStack(P[])
centers = P[]
square_radii = F[]
empty_center = F(NaN) * first(pts)
GärtnerBoundary(centers, square_radii, projector, empty_center)
end
prefix(pts, i) = view(pts, 1:i)
Base.length(b::GärtnerBoundary) = length(b.centers)
Base.isempty(b::GärtnerBoundary) = length(b) == 0
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 = first(b.centers)
center = b.centers[end]
C = center - q0
r2 = b.square_radii[end]
Qm = pt - q0
M = b.projector
Qm_bar = M*Qm
residue = Qm - Qm_bar
e = sqdist(Qm, C) - r2
z = 2*sqnorm(residue)
tol = 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
b
end
struct NSphere{P,F}
center::P
sqradius::F
end
function isinside(pt, nsphere::NSphere; atol=0, rtol=0)
r2 = sqdist(pt, center(nsphere))
R2 = sqradius(nsphere)
r2 <= R2 || isapprox(r2, R2;atol=atol^2,rtol=rtol^2)
end
function allinside(pts, nsphere; kw...)
for pt in pts
isinside(pt, nsphere; kw...) || return false
end
true
end
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
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)
function NSphere(b::GärtnerBoundary)
if isempty(b)
c = b.empty_center
r2 = zero(eltype(c))
else
c = b.centers[end]
r2 = b.square_radii[end]
end
NSphere(c,r2)
end
function ismaxlength(b::GärtnerBoundary)
dim = length(b.empty_center)
length(b) == dim + 1
end
function welzl!(pts, bdry)
bdry_len = length(bdry)
support_count = 0
nsphere = 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 = T(-Inf)
k_max = 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 = GärtnerBoundary(pts)
t = 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)
if e > eps(F)
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 + 1
s = s_new + 1
else
return nsphere
end
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>{{out}}
<pre>
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
</pre>
| |||