Smallest enclosing circle problem: Difference between revisions
m
→{{header|Julia}}
Line 40:
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
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)
Line 58 ⟶ 52:
return true
end
q0, center = first(b.centers), b.centers[end]▼
▲ q0 = first(b.centers)
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))▼
▲ tol = eps(eltype(pt)) * max(r2, one(r2))
isstable = abs(z) > tol
if isstable
Line 91 ⟶ 76:
pop!(b.projector)
end
return b
end
Line 100 ⟶ 85:
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)▼
▲ 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)
Line 121 ⟶ 99:
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)
Line 127 ⟶ 108:
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)
▲ r2 = 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)
Line 167 ⟶ 136:
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]
Line 182 ⟶ 150:
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
Line 189 ⟶ 156:
P = eltype(pts)
F = eltype(P)
t, s = s
end
return nsphere
end
function testwelzl()
testdata =[
| |||