Smallest enclosing circle problem: Difference between revisions
Content added Content deleted
(→{{header|Phix}}: added n-dimensional translation) |
|||
Line 411: | Line 411: | ||
Points {{0,0},{0,1},{1,0}} ==> centre {0.5,0.5} radius 0.70710678118655 |
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 |
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> |
</pre> |
||