Rebol3 Code Examplex
Smallest enclosing circle problem
Find the unique smallest circle that contains a given set of points in the plane.
Rebol [
title: "Rosetta code: Smallest enclosing circle problem"
file: %Smallest_enclosing_circle_problem.r3
url: https://rosettacode.org/wiki/Smallest_enclosing_circle_problem
]
welzl: function/with [
"Return the smallest enclosing circle of a block of pair! points using Welzl's algorithm"
points [block!]
][
;; shuffle points to ensure expected O(n) performance
welzl-rec random copy points []
][
make-circle: func [c r] [
reduce/no-set [centre: c radius: r]
]
encloses: func [point circle] [
;; true when point lies within or on the circle boundary
circle/radius >= distance point circle/centre
]
circle-from-two-points: func [a b] [
;; circle with diameter a-b
make-circle
as-pair (a/x + b/x) / 2.0 (a/y + b/y) / 2.0
(distance a b) / 2.0
]
circle-from-three-points: function [a b c][
;; circumscribed circle via perpendicular bisector intersection
ba: as-pair (b/x - a/x) (b/y - a/y)
ca: as-pair (c/x - a/x) (c/y - a/y)
bb: (ba/x * ba/x) + (ba/y * ba/y)
cc: (ca/x * ca/x) + (ca/y * ca/y)
dd: ((ba/x * ca/y) - (ba/y * ca/x)) * 2.0 ;; 2× cross product
centre: as-pair
(((ca/y * bb) - (ba/y * cc)) / dd) + a/x
(((ba/x * cc) - (ca/x * bb)) / dd) + a/y
make-circle centre distance a centre
]
circle-from-points: func [points][
;; base-case dispatcher: 0-3 boundary points determine a unique circle
switch length? points [
0 [make-circle 0.0x0.0 0.0]
1 [make-circle points/1 0.0]
2 [circle-from-two-points points/1 points/2]
3 [circle-from-three-points points/1 points/2 points/3]
]
]
welzl-rec: function [pts boundary][
;; base case: no points left or boundary fully determines the circle
if any [empty? pts 3 = length? boundary] [
return circle-from-points boundary
]
point: take/last pts: copy pts ;; pick next candidate (input is pre-shuffled)
candidate: welzl-rec pts boundary
if encloses point candidate [return candidate] ;; already inside — done
;; point must lie on the boundary; recurse with it added to boundary
welzl-rec pts append copy boundary point
]
]
tests: [
[0.0x0.0 0.0x1.0 1.0x0.0]
[5.0x-2.0 -3.0x-2.0 -2.0x5.0 1.0x6.0 0.0x2.0]
[0.0x0.0 -2.0x-1.0 3.0x-4.0 2.0x8.0 3.0x11.0 -8.0x-2.0 -14.0x-6.0 7.0x3.0 10.0x4.0 -1.0x4.0]
]
foreach test tests [
circle: welzl test
print [
"Centre: (" circle/centre/x "," circle/centre/y "), Radius:" circle/radius
]
]Output:
Centre: ( 0.5 , 0.5 ), Radius: 0.707106781186548
Centre: ( 1.0 , 1.0 ), Radius: 5.0
Centre: ( -2.0 , -1.0 ), Radius: 13.0