Smallest enclosing circle problem: Difference between revisions
SqrtNegInf (talk | contribs) m (→{{header|Raku}}: various idiomatifications) |
(Realize in F#) |
||
(8 intermediate revisions by 6 users not shown) | |||
Line 12: | Line 12: | ||
=={{header|11l}}== |
|||
{{trans|Nim}} |
|||
{{trans|Go}} |
|||
<syntaxhighlight lang="11l">T Point = (Float, Float) |
|||
T Circle |
|||
Point c |
|||
Float r |
|||
F (c, r) |
|||
.c = c |
|||
.r = r |
|||
F String() |
|||
R ‘Center ’(.c)‘, Radius ’(.r) |
|||
F getCircleCenter(Float bx, by, cx, cy) |
|||
V b = bx * bx + by * by |
|||
V c = cx * cx + cy * cy |
|||
V d = bx * cy - by * cx |
|||
R Point((cy * b - by * c) / (2 * d), (bx * c - cx * b) / (2 * d)) |
|||
F contains(Circle ci, Point p) |
|||
R sqlen(ci.c - p) <= ci.r ^ 2 |
|||
F encloses(Circle ci, [Point] ps) |
|||
L(p) ps |
|||
I !contains(ci, p) |
|||
R 0B |
|||
R 1B |
|||
F circleFrom3(Point a, b, c) |
|||
V i = getCircleCenter(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y) |
|||
i.x += a.x |
|||
i.y += a.y |
|||
R Circle(i, distance(i, a)) |
|||
F circleFrom2(Point a, b) |
|||
V c = Point((a.x + b.x) * 0.5, (a.y + b.y) * 0.5) |
|||
R Circle(c, distance(a, b) / 2) |
|||
F secTrivial([Point] rs) |
|||
I rs.empty |
|||
R Circle(Point(0, 0), 0) |
|||
I rs.len == 1 |
|||
R Circle(rs[0], 0) |
|||
I rs.len == 2 |
|||
R circleFrom2(rs[0], rs[1]) |
|||
assert(rs.len == 3, ‘There shouldn't be more than three points.’) |
|||
L(i) 2 |
|||
L(j) i + 1 .< 3 |
|||
V c = circleFrom2(rs[i], rs[j]) |
|||
I encloses(c, rs) |
|||
R c |
|||
R circleFrom3(rs[0], rs[1], rs[2]) |
|||
F welzlHelper([Point] &ps, [Point] rs, Int n) |
|||
V rc = copy(rs) |
|||
I n == 0 | rc.len == 3 |
|||
R secTrivial(rc) |
|||
V idx = random:(n) |
|||
V p = ps[idx] |
|||
swap(&ps[idx], &ps[n - 1]) |
|||
V d = welzlHelper(&ps, rc, n - 1) |
|||
I contains(d, p) |
|||
R d |
|||
rc.append(p) |
|||
R welzlHelper(&ps, rc, n - 1) |
|||
F welzl([Point] ps) |
|||
V pc = copy(ps) |
|||
random:shuffle(&pc) |
|||
R welzlHelper(&pc, [Point](), pc.len) |
|||
V Tests = [[Point(0.0, 0.0), Point(0.0, 1.0), Point(1.0, 0.0)], |
|||
[Point(5.0, -2.0), Point(-3.0, -2.0), Point(-2.0, 5.0), Point(1.0, 6.0), Point(0.0, 2.0)]] |
|||
L(test) Tests |
|||
print(welzl(test))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Center (0.5, 0.5), Radius 0.707106781 |
|||
Center (1, 1), Radius 5 |
|||
</pre> |
|||
=={{header|F_Sharp|F#}}== |
|||
This task uses [https://rosettacode.org/wiki/Centre_and_radius_of_a_circle_passing_through_3_points_in_a_plane#F# Centre and radius of a circle passing through 3 points in a plane (F#)] |
|||
<syntaxhighlight lang="fsharp"> |
|||
// Smallest enclosing circle problem. Nigel Galloway: February 22nd., 2024 |
|||
let mcc points= |
|||
let fN g=points|>List.map(fun n->(n,d g n))|>List.maxBy snd |
|||
let P1,P2=Seq.allPairs points points|>Seq.maxBy(fun(n,g)->d n g) |
|||
let c1,r1=let ((nx,ny) as n),((gx,gy) as g)=P1,P2 in (((nx+gx)/2.0,(ny+gy)/2.0),(d n g)/2.0) |
|||
match fN c1 with (_,d) when d=r1->(c1,r1) |
|||
|(P3,_)->let c2,r2=circ P1 P2 P3 |
|||
match fN c2 with (_,d) when d=r2->(c2,r2) |
|||
|(P4,_)->[circ P1 P3 P4;circ P2 P3 P4]|>List.minBy snd |
|||
let testMCC n=let centre,radius=mcc n in printfn "Minimum Covering Circle is centred at %A with a radius of %f" centre radius |
|||
testMCC [(0.0, 0.0);(0.0, 1.0);(1.0, 0.0)] |
|||
testMCC [(5.0, -2.0);(-3.0, -2.0);(-2.0, 5.0);(1.0, 6.0);(0.0, 2.0)] |
|||
testMCC [(15.85,6.05);(29.82,22.11);(4.84,20.32);(25.47,29.46);(33.65,17.31);(7.30,10.43);(28.15,6.67);(20.85,11.47);(22.83,2.07);(26.28,23.15);(14.39,30.24);(30.24,25.23)] |
|||
testMCC [(15.85,6.05);(29.82,22.11);(4.84,20.32);(25.47,29.46);(33.65,17.31);(7.30,10.43);(28.15,6.67);(20.85,11.47);(22.83,2.08);(26.28,23.15);(14.39,30.24);(30.24,25.23)] |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Minimum Covering Circle is centred at (0.5, 0.5) with a radius of 0.707107 |
|||
Minimum Covering Circle is centred at (1.0, 1.0) with a radius of 5.000000 |
|||
Minimum Covering Circle is centred at (18.97851566, 16.2654108) with a radius of 14.708624 |
|||
Minimum Covering Circle is centred at (18.97952365, 16.27401209) with a radius of 14.707010 |
|||
</pre> |
|||
=={{header|Go}}== |
=={{header|Go}}== |
||
{{trans|Wren}} |
{{trans|Wren}} |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 142: | Line 256: | ||
fmt.Println(welzl(test)) |
fmt.Println(welzl(test)) |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 148: | Line 262: | ||
{(0.500000, 0.500000), 0.707107} |
{(0.500000, 0.500000), 0.707107} |
||
{(1.000000, 1.000000), 5.000000} |
{(1.000000, 1.000000), 5.000000} |
||
</pre> |
|||
=={{header|jq}}== |
|||
{{works with|jq}} |
|||
'''Works with gojq, the Go implementation of jq''' |
|||
'''Adapted from [[#Wren|Wren]]''' |
|||
In the following, a Point (x,y) is represented by a JSON array [x,y], |
|||
and a Circle by a JSON object of the form {center, radius}. |
|||
<syntaxhighlight lang="jq"> |
|||
# Vector sum of the input array of Points |
|||
def sum: transpose | map(add); |
|||
# The square of the distance between two points |
|||
def distSq: |
|||
def sq: . * .; |
|||
. as [$a, $b] |
|||
| (($a[0] - $b[0]) | sq) +(($a[1] - $b[1]) | sq) ; |
|||
# The distance between two points |
|||
def distance: |
|||
distSq | sqrt; |
|||
# Does the input Circle contain the point $p? |
|||
def contains($p): |
|||
([.center, $p] | distSq) <= .radius * .radius; |
|||
# Does the input Circle contain the given array of points? |
|||
def encloses($ps): |
|||
. as $c |
|||
| all($ps[]; . as $p | $c | contains($p)); |
|||
# The center of a circle defined by 3 points |
|||
def getCircleCenter(bx; by; cx; cy): |
|||
(bx*bx + by*by) as $b |
|||
| (cx*cx + cy*cy) as $c |
|||
| (bx*cy - by*cx) as $d |
|||
| [(cy*$b - by*$c) / (2 * $d), (bx*$c - cx*$b) / (2 * $d)]; |
|||
# The (smallest) circle that intersects 1, 2 or 3 distinct points |
|||
def Circle: |
|||
if length == 1 then {center: .[0], radius: 0} |
|||
elif length == 2 |
|||
then {center: (sum | map(./2)), |
|||
radius: (distance/2) } |
|||
elif length == 3 |
|||
then . as [$a, $b, $c] |
|||
| ([getCircleCenter($b[0] - $a[0]; |
|||
$b[1] - $a[1]; |
|||
$c[0] - $a[0]; |
|||
$c[1] - $a[1]), $a] | sum) as $i |
|||
| {center: $i, radius: ([$i, $a]|distance)} |
|||
else "Circle/0 expects an input array of from 1 to 3 Points" | error |
|||
end; |
|||
# The smallest enclosing circle for n <= 3 |
|||
def secTrivial: |
|||
. as $rs |
|||
| length as $length |
|||
| if $length > 3 then "Internal error: secTrivial given more than 3 points." | error |
|||
elif $length == 0 then [[0, 0]] | Circle |
|||
elif $length <= 2 then Circle |
|||
else first( |
|||
range(0; 2) as $i |
|||
| range($i+1; 3) as $j |
|||
| ([$rs[$i], $rs[$j]] | Circle) |
|||
| select(encloses($rs)) ) |
|||
// Circle |
|||
end; |
|||
# Use the Welzl algorithm to find the SEC of the incoming array of points |
|||
def welzl: |
|||
# Helper function: |
|||
# $rs is an array of points on the boundary; |
|||
# $n keeps track of how many points remain: |
|||
def welzl_($rs; $n): |
|||
if $n == 0 or ($rs|length) == 3 |
|||
then $rs | secTrivial |
|||
else 0 as $idx # or Rand($n) |
|||
| .[$idx] as $p |
|||
| .[$idx] = .[$n-1] |
|||
| .[$n-1] = $p |
|||
| welzl_($rs; $n-1) as $d |
|||
| if ($d | contains($p)) then $d |
|||
else welzl_($rs + [$p]; $n-1) |
|||
end |
|||
end ; |
|||
welzl_([]; length); |
|||
def tests: |
|||
[ [0, 0], [0, 1], [1, 0] ], |
|||
[ [5, -2], [-3, -2], [-2, 5], [1, 6], [0, 2] ] |
|||
; |
|||
tests |
|||
| welzl |
|||
| "Circle(\(.center), \(.radius))" |
|||
</syntaxhighlight> |
|||
{{output}} |
|||
<pre> |
|||
Circle([0.5,0.5], 0.7071067811865476) |
|||
Circle([1,1], 5) |
|||
</pre> |
</pre> |
||
=={{header|Julia}}== |
=={{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. |
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. |
||
< |
<syntaxhighlight lang="julia">import Base.pop!, Base.push!, Base.length, Base.* |
||
using LinearAlgebra, Random |
using LinearAlgebra, Random |
||
Line 323: | Line 542: | ||
testwelzl() |
testwelzl() |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
For points: [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]] |
For points: [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]] |
||
Line 341: | Line 560: | ||
Radius is: 2.7946020036995454 |
Radius is: 2.7946020036995454 |
||
</pre> |
</pre> |
||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
|||
<syntaxhighlight lang="mathematica">f=BoundingRegion[#,"MinBall"]&; |
|||
f[{{0,0},{0,1},{1,0}}] |
|||
f[{{5,-2},{-3,-2},{-2,5},{1,6},{0,2}}] |
|||
f[{{2,4,-1},{1,5,-3},{8,-4,1},{3,9,-5}}] |
|||
f[{{-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`}}]</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Ball[{0.5,0.5},0.707107] |
|||
Ball[{1.,1.},5.] |
|||
Disk[{11/2,5/2,-2},Sqrt[115/2]] |
|||
Ball[{0.0405866,0.555668,-0.231368,0.607407,-0.200346},2.7946]</pre> |
|||
=={{header|Nim}}== |
=={{header|Nim}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
{{trans|Wren}} |
{{trans|Wren}} |
||
< |
<syntaxhighlight lang="nim">import math, random, strutils |
||
type |
type |
||
Line 446: | Line 677: | ||
for test in Tests: |
for test in Tests: |
||
echo welzl(test)</ |
echo welzl(test)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 454: | Line 685: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
Based on the same code as Wren, and likewise limited to 2D circles - I believe (but cannot prove) the main barrier to coping with more than two dimensions lies wholly within the circle_from3() routine. |
Based on the same code as Wren, and likewise limited to 2D circles - I believe (but cannot prove) the main barrier to coping with more than two dimensions lies wholly within the circle_from3() routine. |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>type point(sequence) return true end type |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
enum CENTRE, RADIUS |
|||
<span style="color: #008080;">type</span> <span style="color: #000000;">point</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #004600;">true</span> <span style="color: #008080;">end</span> <span style="color: #008080;">type</span> |
|||
type circle(sequence) return true end type |
|||
<span style="color: #008080;">enum</span> <span style="color: #000000;">CENTRE</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">RADIUS</span> |
|||
<span style="color: #008080;">type</span> <span style="color: #000000;">circle</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #004600;">true</span> <span style="color: #008080;">end</span> <span style="color: #008080;">type</span> |
|||
function distance(point a, b) |
|||
return sqrt(sum(sq_power(sq_sub(a,b),2))) |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">distance</span><span style="color: #0000FF;">(</span><span style="color: #000000;">point</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">),</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)))</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
function in_circle(point p, circle c) |
|||
return distance(p,c[CENTRE]) <= c[RADIUS] |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">in_circle</span><span style="color: #0000FF;">(</span><span style="color: #000000;">point</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">circle</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">)</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">distance</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">CENTRE</span><span style="color: #0000FF;">])</span> <span style="color: #0000FF;"><=</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">RADIUS</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
function circle_from2(point a, b) |
|||
-- return the smallest circle that intersects 2 points: |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">circle_from2</span><span style="color: #0000FF;">(</span><span style="color: #000000;">point</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> |
|||
point midpoint = sq_div(sq_add(a,b),2) |
|||
<span style="color: #000080;font-style:italic;">-- return the smallest circle that intersects 2 points:</span> |
|||
atom halfdiameter = distance(a,b)/2 |
|||
<span style="color: #000000;">point</span> <span style="color: #000000;">midpoint</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">),</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> |
|||
circle res = { midpoint, halfdiameter } |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">halfdiameter</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">distance</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">2</span> |
|||
return res |
|||
<span style="color: #000000;">circle</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">midpoint</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">halfdiameter</span> <span style="color: #0000FF;">}</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
function circle_from3(point a, b, c) |
|||
-- return a unique circle that intersects three points |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">circle_from3</span><span style="color: #0000FF;">(</span><span style="color: #000000;">point</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">)</span> |
|||
point bma = sq_sub(b,a), |
|||
<span style="color: #000080;font-style:italic;">-- return a unique circle that intersects three points </span> |
|||
cma = sq_sub(c,a) |
|||
<span style="color: #000000;">point</span> <span style="color: #000000;">bma</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">),</span> |
|||
atom {{aX,aY},{bX,bY},{cX,cY}} = {a,bma,cma}, |
|||
<span style="color: #000000;">cma</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">c</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> |
|||
B = sum(sq_power(bma,2)), |
|||
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">aX</span><span style="color: #0000FF;">,</span><span style="color: #000000;">aY</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">bX</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bY</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">cX</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cY</span><span style="color: #0000FF;">}}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bma</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cma</span><span style="color: #0000FF;">},</span> |
|||
C = sum(sq_power(cma,2)), |
|||
<span style="color: #000000;">B</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">bma</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)),</span> |
|||
D = (bX*cY - bY*cX)*2 |
|||
<span style="color: #000000;">C</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">cma</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)),</span> |
|||
point centre = {(cY*B - bY*C)/D + aX, |
|||
<span style="color: #000000;">D</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">bX</span><span style="color: #0000FF;">*</span><span style="color: #000000;">cY</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">bY</span><span style="color: #0000FF;">*</span><span style="color: #000000;">cX</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">2</span> |
|||
(bX*C - cX*B)/D + aY } |
|||
<span style="color: #000000;">point</span> <span style="color: #000000;">centre</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{(</span><span style="color: #000000;">cY</span><span style="color: #0000FF;">*</span><span style="color: #000000;">B</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">bY</span><span style="color: #0000FF;">*</span><span style="color: #000000;">C</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">D</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">aX</span><span style="color: #0000FF;">,</span> |
|||
atom radius = distance(centre,a) -- (=== b,c) |
|||
<span style="color: #0000FF;">(</span><span style="color: #000000;">bX</span><span style="color: #0000FF;">*</span><span style="color: #000000;">C</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">cX</span><span style="color: #0000FF;">*</span><span style="color: #000000;">B</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">D</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">aY</span> <span style="color: #0000FF;">}</span> |
|||
circle res = { centre, radius } |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">radius</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">distance</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centre</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (=== b,c)</span> |
|||
return res |
|||
<span style="color: #000000;">circle</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">centre</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">radius</span> <span style="color: #0000FF;">}</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
function trivial(sequence r) |
|||
integer l = length(r) |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">trivial</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> |
|||
switch l do |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> |
|||
case 0: return {{0,0},0} |
|||
<span style="color: #008080;">switch</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span> |
|||
case 1: return {r[1],0} |
|||
<span style="color: #008080;">case</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">:</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">},</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> |
|||
case 2: return circle_from2(r[1],r[2]) |
|||
<span style="color: #008080;">case</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">:</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> |
|||
case 3: return circle_from3(r[1],r[2],r[3]) |
|||
<span style="color: #008080;">case</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">:</span> <span style="color: #008080;">return</span> <span style="color: #000000;">circle_from2</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">])</span> |
|||
end switch |
|||
<span style="color: #008080;">case</span> <span style="color: #000000;">3</span><span style="color: #0000FF;">:</span> <span style="color: #008080;">return</span> <span style="color: #000000;">circle_from3</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">],</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">3</span><span style="color: #0000FF;">])</span> |
|||
end function |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">switch</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
function welzlr(sequence p, r) |
|||
if p={} or length(r)=3 then return trivial(r) end if |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">welzlr</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> |
|||
point p1 = p[1] |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">={}</span> <span style="color: #008080;">or</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">3</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">trivial</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
p = p[2..$] |
|||
<span style="color: #000000;">point</span> <span style="color: #000000;">p1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> |
|||
circle d = welzlr(p, r) |
|||
<span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">..$]</span> |
|||
if in_circle(p1,d) then return d end if |
|||
<span style="color: #000000;">circle</span> <span style="color: #000000;">d</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">welzlr</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> |
|||
return welzlr(p, append(r,p1)) |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">in_circle</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">d</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">welzlr</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">),</span><span style="color: #000000;">p1</span><span style="color: #0000FF;">))</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
procedure welzl(sequence p) |
|||
circle c = welzlr(shuffle(p),{}) |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">welzl</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
string s = sprintf("centre %v radius %.14g",c) |
|||
<span style="color: #000000;">circle</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">welzlr</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">shuffle</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">),{})</span> |
|||
printf(1,"Points %v ==> %s\n",{p,s}) |
|||
<span style="color: #004080;">string</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"centre %v radius %.14g"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">)</span> |
|||
end procedure |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Points %v ==> %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
constant tests = {{{0, 0},{ 0, 1},{ 1,0}}, |
|||
{{5,-2},{-3,-2},{-2,5},{1,6},{0,2}}} |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}},</span> |
|||
papply(tests,welzl)</lang> |
|||
<span style="color: #0000FF;">{{</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},{-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},{-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">6</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">}}}</span> |
|||
<span style="color: #7060A8;">papply</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">,</span><span style="color: #000000;">welzl</span><span style="color: #0000FF;">)</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 526: | Line 760: | ||
{{trans|Python}} |
{{trans|Python}} |
||
Uses the last test from Julia, however, since that's given more accurately. |
Uses the last test from Julia, however, since that's given more accurately. |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>constant inf = 1e300*1e300, |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
nan = -(inf/inf), |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">inf</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1e300</span><span style="color: #0000FF;">*</span><span style="color: #000000;">1e300</span><span style="color: #0000FF;">,</span> |
|||
eps = 1e-8 |
|||
<span style="color: #000000;">nan</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-(</span><span style="color: #000000;">inf</span><span style="color: #0000FF;">/</span><span style="color: #000000;">inf</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">eps</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1e-8</span> |
|||
function sqnorm(sequence p) |
|||
return sum(sq_power(p,2)) |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">sqnorm</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">))</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
function sqdist(sequence p1, p2) |
|||
return sqnorm(sq_sub(p1,p2)) |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">sqdist</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p2</span><span style="color: #0000FF;">)</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">sqnorm</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p2</span><span style="color: #0000FF;">))</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
class ProjectorStack |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">projector_mult</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">vs</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> |
|||
-- Stack of points that are shifted / projected to put first one at origin. |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">v</span><span style="color: #0000FF;">))</span> |
|||
-- |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">vs</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
sequence vs = {} |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">vi</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">vs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">vi</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">vi</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">))))</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">s</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #000080;font-style:italic;">-- Gaertner Boundary:</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">empty_center</span><span style="color: #0000FF;">,</span> |
|||
<span style="color: #000000;">centers</span><span style="color: #0000FF;">,</span> |
|||
<span style="color: #000000;">square_radii</span><span style="color: #0000FF;">,</span> |
|||
<span style="color: #000080;font-style:italic;">-- Stack of points that are shifted / projected to put first one at origin:</span> |
|||
<span style="color: #000000;">projector</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">push_if_stable</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">pt</span><span style="color: #0000FF;">)</span> |
|||
procedure push(sequence v) |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">==</span> <span style="color: #000000;">0</span> <span style="color: #008080;">then</span> |
|||
vs = append(vs, v) |
|||
<span style="color: #000000;">square_radii</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">square_radii</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> |
|||
end procedure |
|||
<span style="color: #000000;">centers</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #004600;">true</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">q0</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">centers</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span> |
|||
<span style="color: #000000;">center</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">centers</span><span style="color: #0000FF;">[$],</span> |
|||
<span style="color: #000000;">C</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q0</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">Qm</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q0</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">Qm_bar</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">projector_mult</span><span style="color: #0000FF;">(</span><span style="color: #000000;">projector</span><span style="color: #0000FF;">,</span><span style="color: #000000;">Qm</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">residue</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">Qm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">Qm_bar</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">r2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">square_radii</span><span style="color: #0000FF;">[$],</span> |
|||
<span style="color: #000000;">e</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sqdist</span><span style="color: #0000FF;">(</span><span style="color: #000000;">Qm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">C</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">r2</span><span style="color: #0000FF;">,</span> |
|||
<span style="color: #000000;">z</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">2</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">sqnorm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">residue</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">tol</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">eps</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">isstable</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">></span> <span style="color: #000000;">tol</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">isstable</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">center_new</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">e</span><span style="color: #0000FF;">/</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">residue</span><span style="color: #0000FF;">))</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">r2new</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">r2</span> <span style="color: #0000FF;">+</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">e</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">e</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">2</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">projector</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">projector</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">residue</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sqnorm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">residue</span><span style="color: #0000FF;">))))</span> |
|||
<span style="color: #000000;">centers</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">center_new</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">square_radii</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">square_radii</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r2new</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">isstable</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">pop_gb</span><span style="color: #0000FF;">()</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">centers</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">centers</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..$-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #000000;">square_radii</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">square_radii</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..$-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">>=</span> <span style="color: #000000;">2</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #000000;">projector</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">projector</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..$-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #000080;font-style:italic;">-- (pop/discard)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">isinside</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">pt</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">r2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sqdist</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">r2</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">nan</span> <span style="color: #008080;">and</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">r2</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">sqradius</span> <span style="color: #008080;">or</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r2</span><span style="color: #0000FF;">-</span><span style="color: #000000;">sqradius</span><span style="color: #0000FF;">)<</span><span style="color: #000000;">eps</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">move_to_front</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pt</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">i</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #0000FF;">{</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">pt</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]}</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">pts</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">ismaxlength</span><span style="color: #0000FF;">()</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">==</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">empty_center</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">_welzl</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">pos</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">support_count</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">center</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">)?</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">[$]:</span><span style="color: #000000;">empty_center</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">sqradius</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">)?</span><span style="color: #000000;">square_radii</span><span style="color: #0000FF;">[$]:</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">ismaxlength</span><span style="color: #0000FF;">()</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">pos</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #008080;">not</span> <span style="color: #000000;">isinside</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #004080;">bool</span> <span style="color: #000000;">isstable</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">push_if_stable</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">isstable</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #0000FF;">{</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">_welzl</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">pop_gb</span><span style="color: #0000FF;">()</span> |
|||
<span style="color: #000000;">pts</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">move_to_front</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">support_count</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">support_count</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">find_max_excess</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">k1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">err_max</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">inf</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">k_max</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">k1</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">1</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">k_max</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pt</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">err</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sqdist</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">sqradius</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">err</span> <span style="color: #0000FF;">></span> <span style="color: #000000;">err_max</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #000000;">err_max</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">err</span> |
|||
<span style="color: #000000;">k_max</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">k</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">err_max</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">k_max</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">welzl</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">points</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">maxiterations</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2000</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pts</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">points</span> |
|||
<span style="color: #000000;">empty_center</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nan</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]))</span> |
|||
<span style="color: #000000;">centers</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> |
|||
<span style="color: #000000;">square_radii</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> |
|||
<span style="color: #000000;">projector</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s_new</span> |
|||
function pop() |
|||
<span style="color: #0000FF;">{</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">_welzl</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
|||
sequence v = vs[$] |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">maxiterations</span> <span style="color: #008080;">do</span> |
|||
vs = vs[1..$-1] |
|||
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">e</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">find_max_excess</span><span style="color: #0000FF;">(</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
return v |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">e</span> <span style="color: #0000FF;"><=</span> <span style="color: #000000;">eps</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
end function |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pt</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">push_if_stable</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">)</span> |
|||
function mult(sequence v) |
|||
<span style="color: #0000FF;">{</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s_new</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">_welzl</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> |
|||
sequence s = repeat(0,length(v)) |
|||
<span style="color: #000000;">pop_gb</span><span style="color: #0000FF;">()</span> |
|||
for i=1 to length(vs) do |
|||
<span style="color: #000000;">pts</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">move_to_front</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> |
|||
sequence vi = vs[i] |
|||
<span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> |
|||
s = sq_add(s,sq_mul(vi,sum(sq_mul(vi, v)))) |
|||
<span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s_new</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> |
|||
end for |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
return s |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"For points: "</span><span style="color: #0000FF;">)</span> <span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">points</span><span style="color: #0000FF;">,{</span><span style="color: #004600;">pp_Indent</span><span style="color: #0000FF;">,</span><span style="color: #000000;">12</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_FltFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%11.8f"</span><span style="color: #0000FF;">})</span> |
|||
end function |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" Center is at: %v\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">center</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" Radius is: %.16g\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sqradius</span><span style="color: #0000FF;">)})</span> |
|||
end class |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
class GaertnerBoundary |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}},</span> |
|||
public: |
|||
<span style="color: #0000FF;">{{</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">}},</span> |
|||
sequence empty_center, |
|||
<span style="color: #0000FF;">{{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">5</span><span style="color: #0000FF;">}},</span> |
|||
centers = {}, |
|||
<span style="color: #0000FF;">{{-</span><span style="color: #000000;">0.6400900432782123</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.643703255134232</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.4016122094706093</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.8959601399652273</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1.1624046608380516</span><span style="color: #0000FF;">},</span> |
|||
square_radii = {} |
|||
<span style="color: #0000FF;">{</span> <span style="color: #000000;">0.5632393652621324</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.015981105190064373</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2.193725940351997</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">0.9094586577358262</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.7165036664470906</span><span style="color: #0000FF;">},</span> |
|||
ProjectorStack projector = new() |
|||
<span style="color: #0000FF;">{-</span><span style="color: #000000;">1.7704367632976061</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.2518283698686299</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1.3810444289625348</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.597516704360172</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.089645656962647</span><span style="color: #0000FF;">},</span> |
|||
end class |
|||
<span style="color: #0000FF;">{</span> <span style="color: #000000;">1.3448578652803103</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.18140877132223002</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">0.4288734015080915</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.53271973321691</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.41896461833399573</span><span style="color: #0000FF;">},</span> |
|||
<span style="color: #0000FF;">{</span> <span style="color: #000000;">0.2132336397213029</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.07659442168765788</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.1486364315310990</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.3386893481333795</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2.3000459709300927</span><span style="color: #0000FF;">},</span> |
|||
function push_if_stable(GaertnerBoundary bound, sequence pt) |
|||
<span style="color: #0000FF;">{</span> <span style="color: #000000;">0.6023153188617328</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.3051735340025314</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0732600437151525</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.1148388039984898</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.047605838564167786</span><span style="color: #0000FF;">},</span> |
|||
if length(bound.centers) == 0 then |
|||
<span style="color: #0000FF;">{</span> <span style="color: #000000;">1.3655523661720959</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.5461416420929995</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">0.09321951900362889</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1.004771137760985</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.6532914656050117</span><span style="color: #0000FF;">},</span> |
|||
bound.square_radii = append(bound.square_radii, 0) |
|||
<span style="color: #0000FF;">{-</span><span style="color: #000000;">0.14974382165751837</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.5375672589202939</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.15845596754003466</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.2750720340454088</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.441247015836271</span><span style="color: #0000FF;">}}}</span> |
|||
bound.centers = {pt} |
|||
<span style="color: #7060A8;">papply</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">,</span><span style="color: #000000;">welzl</span><span style="color: #0000FF;">)</span> |
|||
return true |
|||
<!--</syntaxhighlight>--> |
|||
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}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 741: | Line 941: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
{{trans|Julia}} |
{{trans|Julia}} |
||
< |
<syntaxhighlight lang="python">import numpy as np |
||
class ProjectorStack: |
class ProjectorStack: |
||
Line 900: | Line 1,100: | ||
print(" Center is at: ", nsphere.center) |
print(" Center is at: ", nsphere.center) |
||
print(" Radius is: ", np.sqrt(nsphere.sqradius), "\n") |
print(" Radius is: ", np.sqrt(nsphere.sqradius), "\n") |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
For points: [[0. 0.] |
For points: [[0. 0.] |
||
Line 937: | Line 1,137: | ||
=={{header|Raku}}== |
=={{header|Raku}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
<lang |
<syntaxhighlight lang="raku" line>class P { has ($.x, $.y) is rw; # Point |
||
method gist { "({$.x~", "~$.y})" } |
method gist { "({$.x~", "~$.y})" } |
||
} |
} |
||
Line 1,003: | Line 1,203: | ||
} |
} |
||
say "Solution for smallest circle enclosing {$_.gist} :\n", welzl $_ for @tests;</ |
say "Solution for smallest circle enclosing {$_.gist} :\n", welzl $_ for @tests;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Solution for smallest circle enclosing ((0, 0) (0, 1) (1, 0)) : |
<pre>Solution for smallest circle enclosing ((0, 0) (0, 1) (1, 0)) : |
||
Line 1,015: | Line 1,215: | ||
It is based on Welzl's algorithm and follows closely the C++ code [https://www.geeksforgeeks.org/minimum-enclosing-circle-set-2-welzls-algorithm/?ref=rp here]. |
It is based on Welzl's algorithm and follows closely the C++ code [https://www.geeksforgeeks.org/minimum-enclosing-circle-set-2-welzls-algorithm/?ref=rp here]. |
||
< |
<syntaxhighlight lang="wren">import "random" for Random |
||
var Rand = Random.new() |
var Rand = Random.new() |
||
Line 1,131: | Line 1,331: | ||
] |
] |
||
for (test in tests) System.print(Circle.welzl(test))</ |
for (test in tests) System.print(Circle.welzl(test))</syntaxhighlight> |
||
{{out}} |
{{out}} |
Revision as of 17:02, 22 February 2024
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.
11l
T Point = (Float, Float)
T Circle
Point c
Float r
F (c, r)
.c = c
.r = r
F String()
R ‘Center ’(.c)‘, Radius ’(.r)
F getCircleCenter(Float bx, by, cx, cy)
V b = bx * bx + by * by
V c = cx * cx + cy * cy
V d = bx * cy - by * cx
R Point((cy * b - by * c) / (2 * d), (bx * c - cx * b) / (2 * d))
F contains(Circle ci, Point p)
R sqlen(ci.c - p) <= ci.r ^ 2
F encloses(Circle ci, [Point] ps)
L(p) ps
I !contains(ci, p)
R 0B
R 1B
F circleFrom3(Point a, b, c)
V i = getCircleCenter(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y)
i.x += a.x
i.y += a.y
R Circle(i, distance(i, a))
F circleFrom2(Point a, b)
V c = Point((a.x + b.x) * 0.5, (a.y + b.y) * 0.5)
R Circle(c, distance(a, b) / 2)
F secTrivial([Point] rs)
I rs.empty
R Circle(Point(0, 0), 0)
I rs.len == 1
R Circle(rs[0], 0)
I rs.len == 2
R circleFrom2(rs[0], rs[1])
assert(rs.len == 3, ‘There shouldn't be more than three points.’)
L(i) 2
L(j) i + 1 .< 3
V c = circleFrom2(rs[i], rs[j])
I encloses(c, rs)
R c
R circleFrom3(rs[0], rs[1], rs[2])
F welzlHelper([Point] &ps, [Point] rs, Int n)
V rc = copy(rs)
I n == 0 | rc.len == 3
R secTrivial(rc)
V idx = random:(n)
V p = ps[idx]
swap(&ps[idx], &ps[n - 1])
V d = welzlHelper(&ps, rc, n - 1)
I contains(d, p)
R d
rc.append(p)
R welzlHelper(&ps, rc, n - 1)
F welzl([Point] ps)
V pc = copy(ps)
random:shuffle(&pc)
R welzlHelper(&pc, [Point](), pc.len)
V Tests = [[Point(0.0, 0.0), Point(0.0, 1.0), Point(1.0, 0.0)],
[Point(5.0, -2.0), Point(-3.0, -2.0), Point(-2.0, 5.0), Point(1.0, 6.0), Point(0.0, 2.0)]]
L(test) Tests
print(welzl(test))
- Output:
Center (0.5, 0.5), Radius 0.707106781 Center (1, 1), Radius 5
F#
This task uses Centre and radius of a circle passing through 3 points in a plane (F#)
// Smallest enclosing circle problem. Nigel Galloway: February 22nd., 2024
let mcc points=
let fN g=points|>List.map(fun n->(n,d g n))|>List.maxBy snd
let P1,P2=Seq.allPairs points points|>Seq.maxBy(fun(n,g)->d n g)
let c1,r1=let ((nx,ny) as n),((gx,gy) as g)=P1,P2 in (((nx+gx)/2.0,(ny+gy)/2.0),(d n g)/2.0)
match fN c1 with (_,d) when d=r1->(c1,r1)
|(P3,_)->let c2,r2=circ P1 P2 P3
match fN c2 with (_,d) when d=r2->(c2,r2)
|(P4,_)->[circ P1 P3 P4;circ P2 P3 P4]|>List.minBy snd
let testMCC n=let centre,radius=mcc n in printfn "Minimum Covering Circle is centred at %A with a radius of %f" centre radius
testMCC [(0.0, 0.0);(0.0, 1.0);(1.0, 0.0)]
testMCC [(5.0, -2.0);(-3.0, -2.0);(-2.0, 5.0);(1.0, 6.0);(0.0, 2.0)]
testMCC [(15.85,6.05);(29.82,22.11);(4.84,20.32);(25.47,29.46);(33.65,17.31);(7.30,10.43);(28.15,6.67);(20.85,11.47);(22.83,2.07);(26.28,23.15);(14.39,30.24);(30.24,25.23)]
testMCC [(15.85,6.05);(29.82,22.11);(4.84,20.32);(25.47,29.46);(33.65,17.31);(7.30,10.43);(28.15,6.67);(20.85,11.47);(22.83,2.08);(26.28,23.15);(14.39,30.24);(30.24,25.23)]
- Output:
Minimum Covering Circle is centred at (0.5, 0.5) with a radius of 0.707107 Minimum Covering Circle is centred at (1.0, 1.0) with a radius of 5.000000 Minimum Covering Circle is centred at (18.97851566, 16.2654108) with a radius of 14.708624 Minimum Covering Circle is centred at (18.97952365, 16.27401209) with a radius of 14.707010
Go
package main
import (
"fmt"
"log"
"math"
"math/rand"
"time"
)
type Point struct{ x, y float64 }
type Circle struct {
c Point
r float64
}
// string representation of a Point
func (p Point) String() string { return fmt.Sprintf("(%f, %f)", p.x, p.y) }
// returns the square of the distance between two points
func distSq(a, b Point) float64 {
return (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y)
}
// returns the center of a circle defined by 3 points
func getCircleCenter(bx, by, cx, cy float64) Point {
b := bx*bx + by*by
c := cx*cx + cy*cy
d := bx*cy - by*cx
return Point{(cy*b - by*c) / (2 * d), (bx*c - cx*b) / (2 * d)}
}
// returns whether a circle contains the point 'p'
func (ci Circle) contains(p Point) bool {
return distSq(ci.c, p) <= ci.r*ci.r
}
// returns whether a circle contains a slice of points
func (ci Circle) encloses(ps []Point) bool {
for _, p := range ps {
if !ci.contains(p) {
return false
}
}
return true
}
// string representation of a Circle
func (ci Circle) String() string { return fmt.Sprintf("{%v, %f}", ci.c, ci.r) }
// returns a unique circle that intersects 3 points
func circleFrom3(a, b, c Point) Circle {
i := getCircleCenter(b.x-a.x, b.y-a.y, c.x-a.x, c.y-a.y)
i.x += a.x
i.y += a.y
return Circle{i, math.Sqrt(distSq(i, a))}
}
// returns smallest circle that intersects 2 points
func circleFrom2(a, b Point) Circle {
c := Point{(a.x + b.x) / 2, (a.y + b.y) / 2}
return Circle{c, math.Sqrt(distSq(a, b)) / 2}
}
// returns smallest enclosing circle for n <= 3
func secTrivial(rs []Point) Circle {
size := len(rs)
if size > 3 {
log.Fatal("There shouldn't be more than 3 points.")
}
if size == 0 {
return Circle{Point{0, 0}, 0}
}
if size == 1 {
return Circle{rs[0], 0}
}
if size == 2 {
return circleFrom2(rs[0], rs[1])
}
for i := 0; i < 2; i++ {
for j := i + 1; j < 3; j++ {
c := circleFrom2(rs[i], rs[j])
if c.encloses(rs) {
return c
}
}
}
return circleFrom3(rs[0], rs[1], rs[2])
}
// helper function for Welzl method
func welzlHelper(ps, rs []Point, n int) Circle {
rc := make([]Point, len(rs))
copy(rc, rs) // 'rs' passed by value so need to copy
if n == 0 || len(rc) == 3 {
return secTrivial(rc)
}
idx := rand.Intn(n)
p := ps[idx]
ps[idx], ps[n-1] = ps[n-1], p
d := welzlHelper(ps, rc, n-1)
if d.contains(p) {
return d
}
rc = append(rc, p)
return welzlHelper(ps, rc, n-1)
}
// applies the Welzl algorithm to find the SEC
func welzl(ps []Point) Circle {
var pc = make([]Point, len(ps))
copy(pc, ps)
rand.Shuffle(len(pc), func(i, j int) {
pc[i], pc[j] = pc[j], pc[i]
})
return welzlHelper(pc, []Point{}, len(pc))
}
func main() {
rand.Seed(time.Now().UnixNano())
tests := [][]Point{
{Point{0, 0}, Point{0, 1}, Point{1, 0}},
{Point{5, -2}, Point{-3, -2}, Point{-2, 5}, Point{1, 6}, Point{0, 2}},
}
for _, test := range tests {
fmt.Println(welzl(test))
}
}
- Output:
{(0.500000, 0.500000), 0.707107} {(1.000000, 1.000000), 5.000000}
jq
Works with gojq, the Go implementation of jq
Adapted from Wren
In the following, a Point (x,y) is represented by a JSON array [x,y], and a Circle by a JSON object of the form {center, radius}.
# Vector sum of the input array of Points
def sum: transpose | map(add);
# The square of the distance between two points
def distSq:
def sq: . * .;
. as [$a, $b]
| (($a[0] - $b[0]) | sq) +(($a[1] - $b[1]) | sq) ;
# The distance between two points
def distance:
distSq | sqrt;
# Does the input Circle contain the point $p?
def contains($p):
([.center, $p] | distSq) <= .radius * .radius;
# Does the input Circle contain the given array of points?
def encloses($ps):
. as $c
| all($ps[]; . as $p | $c | contains($p));
# The center of a circle defined by 3 points
def getCircleCenter(bx; by; cx; cy):
(bx*bx + by*by) as $b
| (cx*cx + cy*cy) as $c
| (bx*cy - by*cx) as $d
| [(cy*$b - by*$c) / (2 * $d), (bx*$c - cx*$b) / (2 * $d)];
# The (smallest) circle that intersects 1, 2 or 3 distinct points
def Circle:
if length == 1 then {center: .[0], radius: 0}
elif length == 2
then {center: (sum | map(./2)),
radius: (distance/2) }
elif length == 3
then . as [$a, $b, $c]
| ([getCircleCenter($b[0] - $a[0];
$b[1] - $a[1];
$c[0] - $a[0];
$c[1] - $a[1]), $a] | sum) as $i
| {center: $i, radius: ([$i, $a]|distance)}
else "Circle/0 expects an input array of from 1 to 3 Points" | error
end;
# The smallest enclosing circle for n <= 3
def secTrivial:
. as $rs
| length as $length
| if $length > 3 then "Internal error: secTrivial given more than 3 points." | error
elif $length == 0 then [[0, 0]] | Circle
elif $length <= 2 then Circle
else first(
range(0; 2) as $i
| range($i+1; 3) as $j
| ([$rs[$i], $rs[$j]] | Circle)
| select(encloses($rs)) )
// Circle
end;
# Use the Welzl algorithm to find the SEC of the incoming array of points
def welzl:
# Helper function:
# $rs is an array of points on the boundary;
# $n keeps track of how many points remain:
def welzl_($rs; $n):
if $n == 0 or ($rs|length) == 3
then $rs | secTrivial
else 0 as $idx # or Rand($n)
| .[$idx] as $p
| .[$idx] = .[$n-1]
| .[$n-1] = $p
| welzl_($rs; $n-1) as $d
| if ($d | contains($p)) then $d
else welzl_($rs + [$p]; $n-1)
end
end ;
welzl_([]; length);
def tests:
[ [0, 0], [0, 1], [1, 0] ],
[ [5, -2], [-3, -2], [-2, 5], [1, 6], [0, 2] ]
;
tests
| welzl
| "Circle(\(.center), \(.radius))"
- Output:
Circle([0.5,0.5], 0.7071067811865476) Circle([1,1], 5)
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, 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]],
[[2.0, 4.0, -1.0], [1.0, 5.0, -3.0], [8.0, -4.0, 1.0], [3.0, 9.0, -5.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()
- 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: [[2.0, 4.0, -1.0], [1.0, 5.0, -3.0], [8.0, -4.0, 1.0], [3.0, 9.0, -5.0]] Center is at: [5.5, 2.5, -2.0] Radius is: 7.582875444051551 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
Mathematica/Wolfram Language
f=BoundingRegion[#,"MinBall"]&;
f[{{0,0},{0,1},{1,0}}]
f[{{5,-2},{-3,-2},{-2,5},{1,6},{0,2}}]
f[{{2,4,-1},{1,5,-3},{8,-4,1},{3,9,-5}}]
f[{{-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`}}]
- Output:
Ball[{0.5,0.5},0.707107] Ball[{1.,1.},5.] Disk[{11/2,5/2,-2},Sqrt[115/2]] Ball[{0.0405866,0.555668,-0.231368,0.607407,-0.200346},2.7946]
Nim
import math, random, strutils
type
Point = tuple[x, y: float]
Circle = tuple[c: Point; r: float]
func `$`(p: Point): string =
## Return the string representation of a point.
"($1, $2)".format(p.x, p.y)
func dist(a, b: Point): float =
## Return the distance between two points.
hypot(a.x - b.x, a.y - b.y)
func getCircleCenter(bx, by, cx, cy: float): Point =
## Return the center of a circle defined by three points.
let
b = bx * bx + by * by
c = cx * cx + cy * cy
d = bx * cy - by * cx
result = ((cy * b - by * c) / (2 * d), (bx * c - cx * b) / (2 * d))
func contains(ci: Circle; p: Point): bool =
## Return whether a circle contains the point 'p'.
## Used by 'in' and 'notin'.
dist(ci.c, p) <= ci.r
func contains(ci: Circle; ps: seq[Point]): bool =
## Return whether a circle contains a slice of points.
## Used by 'in'.
for p in ps:
if p notin ci: return false
result = true
func `$`(ci: Circle): string =
## Return the string representation of a circle.
"Center $1, Radius $2".format(ci.c, ci.r)
func circleFrom3(a, b, c: Point): Circle =
## Return the smallest circle that intersects three points.
var i = getCircleCenter(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y)
i.x += a.x
i.y += a.y
result = (c: i, r: dist(i, a))
func circleFrom2(a, b: Point): Circle =
## Return the smallest circle that intersects two points.
let c = ((a.x + b.x) * 0.5, (a.y + b.y) * 0.5)
result = (c: c, r: dist(a, b) * 0.5)
func secTrivial(rs: seq[Point]): Circle =
## Return the smallest enclosing circle for n <= 3.
case rs.len
of 0: return ((0.0, 0.0), 0.0)
of 1: return (rs[0], 0.0)
of 2: return circleFrom2(rs[0], rs[1])
of 3: discard
else: raise newException(ValueError, "There shouldn't be more than three points.")
# Three points.
for i in 0..1:
for j in (i + 1)..2:
let c = circleFrom2(rs[i], rs[j])
if rs in c: return c
result = circleFrom3(rs[0], rs[1], rs[2])
proc welzl(ps: var seq[Point]; rs: seq[Point]; n: int): Circle =
## Helper function for Welzl method.
var rc = rs
if n == 0 or rc.len == 3: return secTrivial(rc)
let idx = rand(n-1)
let p = ps[idx]
swap ps[idx], ps[n-1]
let d = welzl(ps, rc, n-1)
if p in d: return d
rc.add p
result = welzl(ps, rc, n-1)
proc welzl(ps: seq[Point]): Circle =
# Applies the Welzl algorithm to find the SEC.
var pc = ps
pc.shuffle()
result = welzl(pc, @[], pc.len)
when isMainModule:
randomize()
const Tests = [@[(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)]]
for test in Tests:
echo welzl(test)
- Output:
Center (0.5, 0.5), Radius 0.7071067811865476 Center (1.0, 1.0), Radius 5.0
Phix
Based on the same code as Wren, and likewise limited to 2D circles - I believe (but cannot prove) the main barrier to coping with more than two dimensions lies wholly within the circle_from3() routine.
with javascript_semantics type point(sequence s) return true end type enum CENTRE, RADIUS type circle(sequence s) return true end type function distance(point a, b) return sqrt(sum(sq_power(sq_sub(a,b),2))) end function function in_circle(point p, circle c) return distance(p,c[CENTRE]) <= c[RADIUS] end function function circle_from2(point a, b) -- return the smallest circle that intersects 2 points: point midpoint = sq_div(sq_add(a,b),2) atom halfdiameter = distance(a,b)/2 circle res = { midpoint, halfdiameter } return res end function function circle_from3(point a, b, c) -- return a unique circle that intersects three points point bma = sq_sub(b,a), cma = sq_sub(c,a) atom {{aX,aY},{bX,bY},{cX,cY}} = {a,bma,cma}, B = sum(sq_power(bma,2)), C = sum(sq_power(cma,2)), D = (bX*cY - bY*cX)*2 point centre = {(cY*B - bY*C)/D + aX, (bX*C - cX*B)/D + aY } atom radius = distance(centre,a) -- (=== b,c) circle res = { centre, radius } return res end function function trivial(sequence r) integer l = length(r) switch l do case 0: return {{0,0},0} case 1: return {r[1],0} case 2: return circle_from2(r[1],r[2]) case 3: return circle_from3(r[1],r[2],r[3]) end switch end function function welzlr(sequence p, r) if p={} or length(r)=3 then return trivial(r) end if point p1 = p[1] p = p[2..$] circle d = welzlr(p, r) if in_circle(p1,d) then return d end if return welzlr(p, append(deep_copy(r),p1)) end function procedure welzl(sequence p) circle c = welzlr(shuffle(p),{}) string s = sprintf("centre %v radius %.14g",c) printf(1,"Points %v ==> %s\n",{p,s}) end procedure constant tests = {{{0, 0},{ 0, 1},{ 1,0}}, {{5,-2},{-3,-2},{-2,5},{1,6},{0,2}}} papply(tests,welzl)
- Output:
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
n-dimensional
Uses the last test from Julia, however, since that's given more accurately.
with javascript_semantics 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 function projector_mult(sequence vs, 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 -- Gaertner Boundary: sequence empty_center, centers, square_radii, -- Stack of points that are shifted / projected to put first one at origin: projector function push_if_stable(sequence pt) if length(centers) == 0 then square_radii = append(square_radii, 0) centers = {pt} return true end if sequence q0 = centers[1], center = centers[$], C = sq_sub(center,q0), Qm = sq_sub(pt,q0), Qm_bar = projector_mult(projector,Qm), residue = sq_sub(Qm,Qm_bar) atom r2 = 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) projector = append(projector,sq_div(residue,sqrt(sqnorm(residue)))) centers = append(centers, center_new) square_radii = append(square_radii, r2new) end if return isstable end function procedure pop_gb() integer n = length(centers) centers = centers[1..$-1] square_radii = square_radii[1..$-1] if n >= 2 then projector = projector[1..$-1] -- (pop/discard) end if end procedure function isinside(sequence pt, center, atom sqradius) atom r2 = sqdist(pt, center) return r2!=nan and (r2<=sqradius or abs(r2-sqradius)<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() return length(centers) == length(empty_center) + 1 end function function _welzl(sequence pts, integer pos) integer support_count = 0, s sequence center = iff(length(centers)?centers[$]:empty_center) atom sqradius = iff(length(centers)?square_radii[$]:0) if ismaxlength() then return {center, sqradius, 0} end if for i=1 to pos do if not isinside(pts[i], center, sqradius) then bool isstable = push_if_stable(pts[i]) if isstable then {center, sqradius, s} = _welzl(pts, i) pop_gb() pts = move_to_front(deep_copy(pts), i) support_count = s + 1 end if end if end for return {center, sqradius, support_count} end function function find_max_excess(sequence center, pts, atom sqradius, 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, center) - 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 empty_center = repeat(nan,length(pts[1])) centers = {} square_radii = {} projector = {} integer t = 1, s_new {sequence center, atom sqradius, integer s} = _welzl(pts, t) for i=1 to maxiterations do atom {e, k} = find_max_excess(center, pts, sqradius, t + 1) if e <= eps then exit end if sequence pt = pts[k] {} = push_if_stable(pt) {center, sqradius, s_new} = _welzl(pts, s) pop_gb() pts = move_to_front(deep_copy(pts), k) 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", {center}) printf(1," Radius is: %.16g\n", {sqrt(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)
- Output:
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
Python
import numpy as np
class ProjectorStack:
"""
Stack of points that are shifted / projected to put first one at origin.
"""
def __init__(self, vec):
self.vs = np.array(vec)
def push(self, v):
if len(self.vs) == 0:
self.vs = np.array([v])
else:
self.vs = np.append(self.vs, [v], axis=0)
return self
def pop(self):
if len(self.vs) > 0:
ret, self.vs = self.vs[-1], self.vs[:-1]
return ret
def __mul__(self, v):
s = np.zeros(len(v))
for vi in self.vs:
s = s + vi * np.dot(vi, v)
return s
class GaertnerBoundary:
"""
GärtnerBoundary
See the passage regarding M_B in Section 4 of Gärtner's paper.
"""
def __init__(self, pts):
self.projector = ProjectorStack([])
self.centers, self.square_radii = np.array([]), np.array([])
self.empty_center = np.array([np.NaN for _ in pts[0]])
def push_if_stable(bound, pt):
if len(bound.centers) == 0:
bound.square_radii = np.append(bound.square_radii, 0.0)
bound.centers = np.array([pt])
return True
q0, center = bound.centers[0], bound.centers[-1]
C, r2 = center - q0, bound.square_radii[-1]
Qm, M = pt - q0, bound.projector
Qm_bar = M * Qm
residue, e = Qm - Qm_bar, sqdist(Qm, C) - r2
z, tol = 2 * sqnorm(residue), np.finfo(float).eps * max(r2, 1.0)
isstable = np.abs(z) > tol
if isstable:
center_new = center + (e / z) * residue
r2new = r2 + (e * e) / (2 * z)
bound.projector.push(residue / np.linalg.norm(residue))
bound.centers = np.append(bound.centers, np.array([center_new]), axis=0)
bound.square_radii = np.append(bound.square_radii, r2new)
return isstable
def pop(bound):
n = len(bound.centers)
bound.centers = bound.centers[:-1]
bound.square_radii = bound.square_radii[:-1]
if n >= 2:
bound.projector.pop()
return bound
class NSphere:
def __init__(self, c, sqr):
self.center = np.array(c)
self.sqradius = sqr
def isinside(pt, nsphere, atol=1e-6, rtol=0.0):
r2, R2 = sqdist(pt, nsphere.center), nsphere.sqradius
return r2 <= R2 or np.isclose(r2, R2, atol=atol**2,rtol=rtol**2)
def allinside(pts, nsphere, atol=1e-6, rtol=0.0):
for p in pts:
if not isinside(p, nsphere, atol, rtol):
return False
return True
def move_to_front(pts, i):
pt = pts[i]
for j in range(len(pts)):
pts[j], pt = pt, np.array(pts[j])
if j == i:
break
return pts
def dist(p1, p2):
return np.linalg.norm(p1 - p2)
def sqdist(p1, p2):
return sqnorm(p1 - p2)
def sqnorm(p):
return np.sum(np.array([x * x for x in p]))
def ismaxlength(bound):
len(bound.centers) == len(bound.empty_center) + 1
def makeNSphere(bound):
if len(bound.centers) == 0:
return NSphere(bound.empty_center, 0.0)
return NSphere(bound.centers[-1], bound.square_radii[-1])
def _welzl(pts, pos, bdry):
support_count, nsphere = 0, makeNSphere(bdry)
if ismaxlength(bdry):
return nsphere, 0
for i in range(pos):
if not isinside(pts[i], nsphere):
isstable = push_if_stable(bdry, pts[i])
if isstable:
nsphere, s = _welzl(pts, i, bdry)
pop(bdry)
move_to_front(pts, i)
support_count = s + 1
return nsphere, support_count
def find_max_excess(nsphere, pts, k1):
err_max, k_max = -np.Inf, k1 - 1
for (k, pt) in enumerate(pts[k_max:]):
err = sqdist(pt, nsphere.center) - nsphere.sqradius
if err > err_max:
err_max, k_max = err, k + k1
return err_max, k_max - 1
def welzl(points, maxiterations=2000):
pts, eps = np.array(points, copy=True), np.finfo(float).eps
bdry, t = GaertnerBoundary(pts), 1
nsphere, s = _welzl(pts, t, bdry)
for i in range(maxiterations):
e, k = find_max_excess(nsphere, pts, t + 1)
if e <= eps:
break
pt = pts[k]
push_if_stable(bdry, pt)
nsphere_new, s_new = _welzl(pts, s, bdry)
pop(bdry)
move_to_front(pts, k)
nsphere = nsphere_new
t, s = s + 1, s_new + 1
return nsphere
if __name__ == '__main__':
TESTDATA =[
np.array([[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]]),
np.array([[5.0, -2.0], [-3.0, -2.0], [-2.0, 5.0], [1.0, 6.0], [0.0, 2.0]]),
np.array([[2.0, 4.0, -1.0], [1.0, 5.0, -3.0], [8.0, -4.0, 1.0], [3.0, 9.0, -5.0]]),
np.random.normal(size=(8, 5))
]
for test in TESTDATA:
nsphere = welzl(test)
print("For points: ", test)
print(" Center is at: ", nsphere.center)
print(" Radius is: ", np.sqrt(nsphere.sqradius), "\n")
- Output:
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.0 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.28550471 -0.12787614 -0.65666165 -0.57319826 -0.09253894] [ 0.08947276 0.64778564 -0.54953015 0.05332253 -0.93055218] [ 1.26721546 -1.26410984 -0.79937865 -0.45096792 -1.23946668] [-0.62653779 0.56561466 0.62403237 0.23226903 -0.95820264] [ 0.90973949 -0.9821474 0.41400032 -1.11268937 0.19568717] [-0.4931165 0.94778404 -0.10534124 0.96431358 0.36504087] [ 1.43895269 -0.69957774 -0.31486014 -1.14980913 0.42550193] [ 1.4722404 1.19711551 -0.18735466 1.69208268 1.04594225]] Center is at: [ 1.00038333 0.08534583 -0.22414667 0.50927606 -0.15936026] Radius is: 2.076492284522762
Raku
class P { has ($.x, $.y) is rw; # Point
method gist { "({$.x~", "~$.y})" }
}
class C { has (P $.c, Numeric $.r); # Circle
method gist { "Center = " ~$.c.gist ~ " and Radius = $.r\n" }
# tests whether a circle contains the point 'p'
method contains(P \p --> Bool) { distSq($.c, p) ≤ $.r² }
# tests whether a circle contains a slice of point
method encloses(@ps --> Bool) { [and] @ps.map: { $.contains($_) } }
}
sub infix:<−> (P \a, P \b) { a.x - b.x, a.y - b.y } # note: Unicode 'minus'
# returns the square of the distance between two points
sub distSq (P \a, P \b) { [+] (a − b)»² }
sub getCenter (\bx, \by, \cx, \cy --> P) {
my (\b,\c,\d) = bx²+by², cx²+cy², bx*cy - by*cx;
P.new: x => (cy*b - by*c) / (2 * d), y => (bx*c - cx*b) / (2 * d)
} # returns the center of a circle defined by 3 points
sub circleFrom3 (P \a, P \b, P \c --> C) {
my \k = $ = getCenter |(b − a), |(c − a);
k.x, k.y Z[+=] a.x, a.y;
C.new: c => k, r => distSq(k, a).sqrt
} # returns a unique circle that intersects 3 points
sub circleFrom2 (P \a, P \b --> C ) {
my \center = P.new: x => ((a.x + b.x) / 2), y => ((a.y + b.y) / 2) ;
C.new: c => center, r => (distSq(a, b).sqrt / 2)
} # returns smallest circle that intersects 2 points
sub secTrivial( @rs --> C ) {
given @rs {
when * == 0 { return C.new: c => (P.new: x => 0, y => 0), r => 0 }
when * == 1 { return C.new: c => @rs[0], r => 0 }
when * == 2 { return circleFrom2 |@rs }
when * == 3 { #`{ no-op } }
when * > 3 { die "There shouldn't be more than 3 points." }
}
for 0, 1 X 1, 2 -> ( \i, \j ) {
return $_ if .encloses(@rs) given circleFrom2 |@rs[i,j]
}
circleFrom3 |@rs
} # returns smallest enclosing circle for n ≤ 3
sub Welzl-helper( @ps is copy, @rs is copy , \n --> C ) {
return secTrivial(@rs) if n == 0 or @rs == 3;
my \p = @ps.shift;
return $_ if .contains(p) given Welzl-helper @ps, @rs, n-1;
Welzl-helper @ps, @rs.append(p), n-1
} # helper function for Welzl method
# applies the Welzl algorithm to find the SEC
sub welzl(@ps --> C) { Welzl-helper @ps.pick(*), [], @ps }
my @tests = (
[ (0,0), (0,1), (1,0) ],
[ (5,-2), (-3,-2), (-2,5), (1,6), (0,2) ]
).map: {
@_.map: { P.new: x => $_[0], y => $_[1] }
}
say "Solution for smallest circle enclosing {$_.gist} :\n", welzl $_ for @tests;
- Output:
Solution for smallest circle enclosing ((0, 0) (0, 1) (1, 0)) : Center = (0.5, 0.5) and Radius = 0.7071067811865476 Solution for smallest circle enclosing ((5, -2) (-3, -2) (-2, 5) (1, 6) (0, 2)) : Center = (1, 1) and Radius = 5
Wren
Well a circle is a two dimensional figure and so, despite any contradictory indications in the task description, that's what this solution provides.
It is based on Welzl's algorithm and follows closely the C++ code here.
import "random" for Random
var Rand = Random.new()
class Point {
// returns the square of the distance between two points
static distSq(a, b) { (a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y) }
// returns the center of a circle defined by 3 points
static getCircleCenter(bx, by, cx, cy) {
var b = bx*bx + by*by
var c = cx*cx + cy*cy
var d = bx*cy - by*cx
return Point.new((cy*b - by*c) / (2 * d), (bx*c - cx*b) / (2 * d))
}
// constructs a new Point object
construct new(x, y) {
_x = x
_y = y
}
// basic properties
x { _x }
x=(o) { _x = o }
y { _y }
y=(o) { _y = o }
// returns a copy of the current instance
copy() { Point.new(_x, _y) }
// string representation
toString { "(%(_x), %(_y))" }
}
class Circle {
// returns a unique circle that intersects 3 points
static from(a, b, c) {
var i = Point.getCircleCenter(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y)
i.x = i.x + a.x
i.y = i.y + a.y
return Circle.new(i, Point.distSq(i, a).sqrt)
}
// returns smallest circle that intersects 2 points
static from(a, b) {
var c = Point.new((a.x + b.x) / 2, (a.y + b.y) / 2)
return Circle.new(c, Point.distSq(a, b).sqrt/2)
}
// returns smallest enclosing circle for n <= 3
static secTrivial(rs) {
var size = rs.count
if (size > 3) Fiber.abort("There shouldn't be more than 3 points.")
if (size == 0) return Circle.new(Point.new(0, 0), 0)
if (size == 1) return Circle.new(rs[0], 0)
if (size == 2) return Circle.from(rs[0], rs[1])
for (i in 0..1) {
for (j in i+1..2) {
var c = Circle.from(rs[i], rs[j])
if (c.encloses(rs)) return c
}
}
return Circle.from(rs[0], rs[1], rs[2])
}
// private helper method for Welzl method
static welzl_(ps, rs, n) {
rs = rs.toList // passed by value so need to copy
if (n == 0 || rs.count == 3) return secTrivial(rs)
var idx = Rand.int(n)
var p = ps[idx]
ps[idx] = ps[n-1]
ps[n-1] = p
var d = welzl_(ps, rs, n-1)
if (d.contains(p)) return d
rs.add(p)
return welzl_(ps, rs, n-1)
}
// applies the Welzl algorithm to find the SEC
static welzl(ps) {
var pc = List.filled(ps.count, null)
for (i in 0...pc.count) pc[i] = ps[i].copy()
Rand.shuffle(pc)
return welzl_(pc, [], pc.count)
}
// constructs a new Circle object from its center point and its radius
construct new(c, r) {
_c = c.copy()
_r = r
}
// basic properties
c { _c }
r { _r }
// returns whether the current instance contains the point 'p'
contains(p) { Point.distSq(_c, p) <= _r * _r }
// returns whether the current instance contains a list of points
encloses(ps) {
for (p in ps) if (!contains(p)) return false
return true
}
// string representation
toString { "Center %(_c), Radius %(_r)" }
}
var tests = [
[Point.new(0, 0), Point.new(0, 1), Point.new(1, 0)],
[Point.new(5, -2), Point.new(-3, -2), Point.new(-2, 5), Point.new(1, 6), Point.new(0, 2)]
]
for (test in tests) System.print(Circle.welzl(test))
- Output:
Center (0.5, 0.5), Radius 0.70710678118655 Center (1, 1), Radius 5