Smallest enclosing circle problem: Difference between revisions
Smallest enclosing circle problem (view source)
Revision as of 20:48, 28 November 2020
, 3 years ago→{{header|Phix}}: added n-dimensional translation
(→{{header|Phix}}: added n-dimensional translation) |
|||
Line 411:
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
</pre>
===n-dimensional===
{{trans|Python}}
Uses the last test from Julia, however, since that's given more accurately.
<lang Phix>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
class ProjectorStack
--
-- Stack of points that are shifted / projected to put first one at origin.
--
sequence vs = {}
procedure push(sequence v)
vs = append(vs, v)
end procedure
function pop()
sequence v = vs[$]
vs = vs[1..$-1]
return v
end function
function mult(sequence 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
end class
class GaertnerBoundary
public:
sequence empty_center,
centers = {},
square_radii = {}
ProjectorStack projector = new()
end class
function push_if_stable(GaertnerBoundary bound, sequence pt)
if length(bound.centers) == 0 then
bound.square_radii = append(bound.square_radii, 0)
bound.centers = {pt}
return true
end if
ProjectorStack M = bound.projector
sequence q0 = bound.centers[1],
center = bound.centers[$],
C = sq_sub(center,q0),
Qm = sq_sub(pt,q0),
Qm_bar = M.mult(Qm),
residue = sq_sub(Qm,Qm_bar)
atom r2 = bound.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)
ProjectorStack p = bound.projector
p.push(sq_div(residue,sqrt(sqnorm(residue))))
bound.centers = append(bound.centers, center_new)
bound.square_radii = append(bound.square_radii, r2new)
end if
return isstable
end function
function pop(GaertnerBoundary bound)
integer n = length(bound.centers)
bound.centers = bound.centers[1..$-1]
bound.square_radii = bound.square_radii[1..$-1]
if n >= 2 then
ProjectorStack p = bound.projector
{} = p.pop()
end if
return bound
end function
class NSphere
public:
sequence center
atom sqradius
end class
function isinside(sequence pt, NSphere nsphere)
atom r2 = sqdist(pt, nsphere.center),
R2 = nsphere.sqradius
if r2=nan then return false end if
return r2<=R2 or abs(r2-R2)<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(GaertnerBoundary bound)
return length(bound.centers) == length(bound.empty_center) + 1
end function
function makeNSphere(GaertnerBoundary bound)
NSphere res
if length(bound.centers) == 0 then
res = new({bound.empty_center, 0.0})
else
res = new({bound.centers[$], bound.square_radii[$]})
end if
return res
end function
function _welzl(sequence pts, integer pos, GaertnerBoundary bdry)
integer support_count = 0
NSphere nsphere = makeNSphere(bdry)
if ismaxlength(bdry) then
return {nsphere, 0}
end if
for i=1 to pos do
if not isinside(pts[i], nsphere) then
bool isstable = push_if_stable(bdry, pts[i])
if isstable then
{nsphere, integer s} = _welzl(pts, i, bdry)
bdry = pop(bdry)
pts = move_to_front(pts, i)
support_count = s + 1
end if
end if
end for
return {nsphere, support_count}
end function
function find_max_excess(NSphere nsphere, sequence pts, 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, nsphere.center) - nsphere.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
GaertnerBoundary bdry = new({repeat(nan,length(pts[1]))})
integer t = 1
{NSphere nsphere, integer s} = _welzl(pts, t, bdry)
for i=1 to maxiterations do
atom {e, k} = find_max_excess(nsphere, pts, t + 1)
if e <= eps then exit end if
sequence pt = pts[k]
{} = push_if_stable(bdry, pt)
{NSphere nsphere_new, integer s_new} = _welzl(pts, s, bdry)
bdry = pop(bdry)
pts = move_to_front(pts, k)
nsphere = nsphere_new
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", {nsphere.center})
printf(1," Radius is: %.16g\n", {sqrt(nsphere.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)</lang>
{{out}}
<pre>
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
</pre>
| |||