A working version of a non-recursive implementation of Welzl's algorithm in Rust. The code below is part of a working solution to a coding challenge problem. It's been run against over a dozen test cases with large data sets and met timing constraints. Consider it tested code.
The only thing missing is the implementation for random() which produces a random value within a range of indexes - easy enough to implement for anyone wanting to use this.
The approach to conversion from recursive to non-recursive is very similar to my other posted answer in Python, but this posted answer is simpler and more complete. It is also more memory efficient.
Below is the primary method of a Circle object that implements the algorithm non-recursively:
fn from_points_internal(p: &mut Vec<Point>,
r: Vec<Point>,
n: usize)
-> Self
{
enum Block { First, Second }
use Block::*;
let mut params = vec![(r, n, None, First)];
let mut results = vec![];
while let Some((mut r, n, p1, jump)) = params.pop() {
match jump {
First if n == 0 || r.len() == 3 => {
let c1 = match r.len() {
3 => Self::from_3_points(&r[0], &r[1], &r[2]),
2 => Self::from_2_points(&r[0], &r[1]),
1 => Self::new(r[0], 0.),
0 => Self::new(Point { x: 0., y: 0. }, 0.),
_ => panic!("Too many `r` points."),
};
results.push(c1);
},
First => {
let idx = random(n);
let p1 = p[idx];
p.swap(idx, n - 1);
params.push((r.clone(), n, Some(p1), Second));
params.push((r, n - 1, None, First));
},
Second => {
let c1 = results.pop().unwrap();
let p1 = p1.unwrap();
if c1.contains_point(&p1) {
results.push(c1);
}
else {
r.push(p1);
params.push((r, n - 1, None, First));
}
}
}
}
results.pop().unwrap()
}
For comparison, here is the recursive version of the same method:
fn from_points_internal( p: &mut Vec<Point>,
mut r: Vec<Point>,
n: usize)
-> Self
{
if n == 0 || r.len() == 3 {
match r.len() {
3 => Self::from_3_points(&r[0], &r[1], &r[2]),
2 => Self::from_2_points(&r[0], &r[1]),
1 => Self::new(r[0], 0.),
0 => Self::new(Point { x: 0., y: 0. }, 0.),
_ => panic!("Too many 'r' points passed to \
Circle::from_points_internal()."),
}
} else {
let idx = random(n); // Get random index.
let p1 = p[idx];
p.swap(idx, n - 1);
let c1 = Self::from_points_internal(p, r.clone(), n - 1);
if c1.contains_point(&p1) {
c1
}
else {
r.push(p1);
Self::from_points_internal(p, r, n - 1)
}
}
}
Here's the full code for Circle and Point (iterative solution):
/// Object representing a circle and providing methods to calculate the minimum
/// enclosing circle of a collection of points.
///
#[derive(Clone, Copy)]
pub struct Circle {
center: Point,
radius: f64,
}
impl Circle {
/// Produce a new circle object with center point and given radius.
///
pub fn new(center: Point, radius: f64) -> Self
{
Self { center, radius }
}
/// Produce an enclosing circle for the slice of points given.
///
pub fn from_points<P>(points: P) -> Self
where
P: IntoIterator<Item = Point>
{
let mut points = points.into_iter().collect::<Vec<_>>();
let len = points.len();
// Shuffle points.
for i in 0..len {
let j = random(len);
points.swap(i, j);
}
Self::from_points_internal(&mut points, vec![], len)
}
/// Returns the radius of the enclosing circle.
///
pub fn radius(&self) -> f64
{
self.radius
}
/// Internal non-recursive function that finds the minimum sized circle that can
/// enclose all points. The time-complexity should work out to `O(n)`.
///
fn from_points_internal(p: &mut Vec<Point>,
r: Vec<Point>,
n: usize)
-> Self
{
enum Block { First, Second }
use Block::*;
let mut params = vec![(r, n, None, First)];
let mut results = vec![];
while let Some((mut r, n, p1, jump)) = params.pop() {
match jump {
First if n == 0 || r.len() == 3 => {
let c1 = match r.len() {
3 => Self::from_3_points(&r[0], &r[1], &r[2]),
2 => Self::from_2_points(&r[0], &r[1]),
1 => Self::new(r[0], 0.),
0 => Self::new(Point { x: 0., y: 0. }, 0.),
_ => panic!("Too many `r` points."),
};
results.push(c1);
},
First => {
let idx = random(n);
let p1 = p[idx];
p.swap(idx, n - 1);
params.push((r.clone(), n, Some(p1), Second));
params.push((r, n - 1, None, First));
},
Second => {
let c1 = results.pop().unwrap();
let p1 = p1.unwrap();
if c1.contains_point(&p1) {
results.push(c1);
}
else {
r.push(p1);
params.push((r, n - 1, None, First));
}
}
}
}
results.pop().unwrap()
}
/// Create circle that intersects the 3 points.
///
fn from_3_points(p1: &Point, p2: &Point, p3: &Point) -> Self
{
let p12 = Point::new(p2.x - p1.x, p2.y - p1.y);
let p13 = Point::new(p3.x - p1.x, p3.y - p1.y);
let c12 = p12.x * p12.x + p12.y * p12.y;
let c13 = p13.x * p13.x + p13.y * p13.y;
let c123 = p12.x * p13.y - p12.y * p13.x;
let cx = (p13.y * c12 - p12.y * c13) / (2. * c123);
let cy = (p12.x * c13 - p13.x * c12) / (2. * c123);
let center = Point::new(cx + p1.x, cy + p1.y);
Self::new(center, center.distance(p1))
}
/// Creates circle that intersects the two points.
///
fn from_2_points(p1: &Point, p2: &Point) -> Self
{
let center = Point::new((p1.x + p2.x) / 2.,
(p1.y + p2.y) / 2.);
Self::new(center, p1.distance(p2) / 2.)
}
/// Returns true if the circile contains the point.
///
fn contains_point(&self, point: &Point) -> bool
{
self.center.distance(point) <= self.radius
}
}
/// Represents a point on a 2D plane.
///
#[derive(Clone, Copy)]
pub struct Point {
x: f64,
y: f64,
}
impl Point {
pub fn new(x: f64, y: f64) -> Self
{
Self { x, y }
}
/// Returns the distance between `self` and `other`.
///
pub fn distance(&self, other: &Point) -> f64
{
((other.x - self.x).powf(2.) +
(other.y - self.y).powf(2.)).sqrt()
}
}
