4
\$\begingroup\$

Introduction

I wanted to implement the Kirkpatrick-Seidel(KPS) convex hull algorithm and chose Rust as my language. I referred to the original KPS paper and tried to implement as closely it as possible.

The benchmarking code, including the implementation

benchmark.rs

#[global_allocator]
static GLOBAL: tikv_jemallocator::Jemalloc = tikv_jemallocator::Jemalloc;

use criterion::{criterion_group, criterion_main, Criterion};

use rustc_hash::FxHashSet as HashSet;
use std::hash::Hash;

// Implementation of the Point vectors
#[derive(Clone, Copy, PartialEq, PartialOrd)]
struct Vec2 {
    x: f32,
    y: f32,
}

impl Vec2 {
    fn new(x: f32, y: f32) -> Self {
        Self { x, y }
    }
}

impl Eq for Vec2 {}

impl Hash for Vec2 {
    fn hash<H: std::hash::Hasher>(&self, state: &mut H) {
        self.x.to_bits().hash(state);
        self.y.to_bits().hash(state);
    }
}


// Implementation of Jarvis March
fn jarvis_march(points: Vec<Vec2>) -> Vec<Vec2> {
    let n = points.len();
    let mut hull = Vec::with_capacity(n);

    if n < 3 {
        return hull;
    }

    let mut l = 0;
    for i in 1..n {
        if points[i].x < points[l].x {
            l = i;
        }
    }

    let mut p = l;
    let mut q;
    loop {
        hull.push(points[p]);

        q = (p + 1) % n;
        for r in 0..n {
            if let Orientation::Counterclockwise = orientation(&points[p], &points[r], &points[q]) {
                q = r;
            }
        }

        p = q;

        if p == l {
            break;
        }
    }

    hull
}

enum Orientation {
    Clockwise,
    Counterclockwise,
    Colinear,
}

#[inline(always)]
fn orientation(p: &Vec2, q: &Vec2, r: &Vec2) -> Orientation {
    let val = (q.y - p.y) * (r.x - q.x) - (q.x - p.x) * (r.y - q.y);
    if val == 0.0 {
        return Orientation::Colinear;
    }
    if val > 0.0 {
        return Orientation::Clockwise;
    }
    Orientation::Counterclockwise
}

// Implementation of KPS and its children functions
fn kirk_patrick_seidel(points: &[Vec2]) -> Vec<Vec2> {
    let mut upper_hull_vec = upper_hull(&points);

    let mut lower_hull_vec = upper_hull(
        &points
            .iter()
            .map(|point| Vec2 {
                x: point.x,
                y: -point.y,
            })
            .collect::<Vec<_>>(),
    );

    lower_hull_vec = lower_hull_vec
        .iter()
        .map(|point| Vec2 {
            x: point.x,
            y: -point.y,
        })
        .collect();

    let (upper_hull_min, upper_hull_max) = upper_hull_vec.iter().fold(
        (upper_hull_vec[0], upper_hull_vec[0]),
        |(min, max), &point| {
            (
                if point.x < min.x { point } else { min },
                if point.x > max.x { point } else { max },
            )
        },
    );

    let (lower_hull_min, lower_hull_max) = lower_hull_vec.iter().fold(
        (lower_hull_vec[0], lower_hull_vec[0]),
        |(min, max), &point| {
            (
                if point.x < min.x { point } else { min },
                if point.x > max.x { point } else { max },
            )
        },
    );

    if upper_hull_max.x == lower_hull_max.x && upper_hull_min.y != lower_hull_min.y {
        upper_hull_vec.push(lower_hull_max);
    }

    if upper_hull_min.x == lower_hull_min.x && upper_hull_min.y != lower_hull_min.y {
        upper_hull_vec.push(lower_hull_min);
    }

    upper_hull_vec.extend(lower_hull_vec);
    upper_hull_vec
}

fn upper_hull(points: &[Vec2]) -> Vec<Vec2> {
    let mut min_point = Vec2 {
        x: f32::MAX,
        y: f32::MIN,
    };
    let mut max_point = Vec2 {
        x: f32::MIN,
        y: f32::MAX,
    };

    for i in points.iter() {
        if i.x < min_point.x || (i.x == min_point.x && i.y > min_point.y) {
            min_point = *i;
        }
        if i.x > max_point.x || (i.x == max_point.x && i.y < max_point.y) {
            max_point = *i;
        }
    }

    let mut temporary = Vec::with_capacity(points.len());
    temporary.push(min_point);

    if min_point == max_point {
        return temporary;
    }

    temporary.push(max_point);
    temporary.extend(
        points
            .iter()
            .filter(|p| p.x > min_point.x && p.x < max_point.x),
    );

    connect(&min_point, &max_point, &temporary)
}

fn connect(min: &Vec2, max: &Vec2, points: &[Vec2]) -> Vec<Vec2> {
    let median = median_of_medians(&points.iter().map(|point| point.x).collect::<Vec<_>>());

    let (left, right) = bridge(points, median);

    let mut left_points = Vec::with_capacity(points.len());
    left_points.push(left);
    left_points.extend(points.iter().filter(|p| p.x < left.x));

    let mut right_points = Vec::with_capacity(points.len());
    right_points.push(right);
    right_points.extend(points.iter().filter(|p| p.x > right.x));

    let mut output = Vec::with_capacity(points.len());
    if left == *min {
        output.push(left);
    } else {
        output.extend(connect(&min, &left, &left_points));
    }

    if right == *max {
        output.push(right);
    } else {
        output.extend(connect(&right, &max, &right_points));
    }

    output
}

fn bridge(points: &[Vec2], median: f32) -> (Vec2, Vec2) {
    let mut candidates: HashSet<&Vec2> = HashSet::default();
    if points.len() == 2 {
        return if points[0].x < points[1].x {
            (points[0], points[1])
        } else {
            (points[1], points[0])
        };
    }

    let mut pairs: Vec<(&Vec2, &Vec2)> = Vec::new();

    for chunk in points.chunks(2) {
        if chunk.len() == 2 {
            if chunk[0] > chunk[1] {
                pairs.push((&chunk[1], &chunk[0]));
            } else {
                pairs.push((&chunk[0], &chunk[1]));
            }
        } else {
            candidates.insert(&chunk[0]);
        }
    }

    let mut slopes = Vec::with_capacity(pairs.len());

    for (point_i, point_j) in pairs.iter() {
        if point_i.x == point_j.x {
            if point_i.y > point_j.y {
                candidates.insert(point_i);
            } else {
                candidates.insert(point_j);
            }
        } else {
            slopes.push((
                point_i,
                point_j,
                (point_i.y - point_j.y) / (point_i.x - point_j.x),
            ));
        }
    }

    let median_slope =
        median_of_medians(&slopes.iter().map(|(_, _, slope)| slope).collect::<Vec<_>>());

    let mut small = Vec::with_capacity(slopes.len());
    let mut equal = Vec::with_capacity(slopes.len());
    let mut large = Vec::with_capacity(slopes.len());

    for slope in slopes.iter() {
        match slope.2 {
            s if s < *median_slope => small.push(slope),
            s if s == *median_slope => equal.push(slope),
            s if s > *median_slope => large.push(slope),
            _ => unreachable!(),
        }
    }

    // set of points with maximum value of p.y - median_slope * p.x
    let mut max_value = f32::MIN;
    let mut max_points = Vec::with_capacity(points.len());
    let mut min_point = None;
    let mut max_point = None;

    for p in points.iter() {
        let value = p.y - median_slope * p.x;
        if value > max_value {
            max_value = value;
            max_points.clear();
            max_points.push(p);
            min_point = Some(p);
            max_point = Some(p);
        } else if (value - max_value).abs() < f32::EPSILON {
            max_points.push(p);
            if let Some(min_p) = min_point {
                if p.x < min_p.x {
                    min_point = Some(p);
                }
            }
            if let Some(max_p) = max_point {
                if p.x > max_p.x {
                    max_point = Some(p);
                }
            }
        }
    }

    let min_point = min_point.unwrap();
    let max_point = max_point.unwrap();

    if min_point.x <= median && max_point.x > median {
        return (*min_point, *max_point);
    } else if max_point.x <= median {
        for (_, point2, _) in large {
            candidates.insert(point2);
        }

        for (_, point2, _) in equal {
            candidates.insert(point2);
        }

        for (point1, point2, _) in small {
            candidates.insert(point2);
            candidates.insert(point1);
        }
    } else if min_point.x > median {
        for (point1, _, _) in small {
            candidates.insert(point1);
        }

        for (point1, _, _) in equal {
            candidates.insert(point1);
        }

        for (point1, point2, _) in large {
            candidates.insert(point2);
            candidates.insert(point1);
        }
    }

    bridge(
        &candidates.into_iter().cloned().collect::<Vec<Vec2>>(),
        median,
    )
}

#[inline(always)]
pub fn median_of_medians<T>(nums: &[T]) -> T
where
    T: Clone + Copy + PartialOrd,
{
    match nums.len() {
        0 => unreachable!(),
        1 => nums[0],
        2..=5 => {
            let mut nums = nums.to_owned();
            nums.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
            nums[nums.len() / 2]
        }
        _ => median_of_medians(
            &nums
                .to_vec()
                .chunks(5)
                .map(|chunk| {
                    let mut chunk = chunk.to_vec();
                    chunk.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
                    chunk[chunk.len() / 2]
                })
                .collect::<Vec<T>>(),
        ),
    }
}

pub fn comparison(c: &mut Criterion) {
    use rand::rngs::StdRng;
    use rand::{Rng, SeedableRng};

    let seed = [32; 32]; // A seed for the RNG. You can put any number here.
    let mut rng: StdRng = SeedableRng::from_seed(seed);

    let points: Vec<Vec2> = (0..100_000)
        .map(|_| {
            Vec2::new(
                rng.gen_range(-50_000..50_000) as f32,
                rng.gen_range(-50_000..50_000) as f32,
            )
        })
        .collect();

    let mut group = c.benchmark_group("Convex-hull Algorithms Comparison");
    group.bench_function("Kirk Patrick Seidel", |b| {
        b.iter(|| kirk_patrick_seidel(&points.clone()))
    });
    group.bench_function("Jarvis March", |b| {
        b.iter(|| jarvis_march(points.clone()))
    });

    group.finish();
}

criterion_group!(benches, comparison);
criterion_main!(benches);

cargo.toml

[package]
name = "profiling"
version = "0.1.0"
edition = "2021"

[dependencies]
criterion = "0.5.1"
rand = "0.8.5"
rustc-hash = "1.1.0"
tikv-jemallocator = "0.5.4"

[profile.release]
debug = 1
lto = false
panic = "abort"

[rust]
debuginfo-level = 1

[target.x86_64-unknown-linux-gnu]
linker = "/usr/bin/clang"
rustflags = ["-Clink-arg=-fuse-ld=lld", "-Clink-arg=-Wl,--no-rosegment"]

[[bench]]
name = "benchmark"
harness = false

Grievances

But, the benchmarking results are horrible, as Jarvis march has complexity \$O(nh)\$ and KPS has complexity \$O(n\log(h))\$, where \$n\$ is the total number of points and \$h\$ is the total number of points on the hull.

The results are

Convex-hull Algorithms Comparison/Kirk Patrick Seidel
                        time:   [57.024 ms 57.348 ms 57.716 ms]
Found 7 outliers among 100 measurements (7.00%)
  6 (6.00%) high mild
  1 (1.00%) high severe
Convex-hull Algorithms Comparison/Jarvis March
                        time:   [4.1446 ms 4.2010 ms 4.2637 ms]
Found 9 outliers among 100 measurements (9.00%)
  4 (4.00%) high mild
  5 (5.00%) high severe

I tried it for more points, but KPS showed no signs of becoming faster than Jarvis March. I am quite new to Rust and newer to profiling (I used tricks from The Rust Performance Book) and am using the criterion crate.

I would appreciate your help and suggestions. Thank you for your time <3

Edit: I believe a per this resource, KPS is slow in general, but I am not sure if that should be the case even against Jarvis March.

\$\endgroup\$
5
  • \$\begingroup\$ I'd suggest looking into SIMD and Google's SwissTable (HashBrown on crates.io). \$\endgroup\$ Apr 2 at 15:57
  • \$\begingroup\$ Vec2::new could be const and you could just #[derive()] Eq and Hash for it. \$\endgroup\$ Apr 2 at 20:11
  • \$\begingroup\$ Also there are a lot of magic numbers in your code. \$\endgroup\$ Apr 2 at 20:16
  • \$\begingroup\$ Hi @RichardNeumann, I didn't quite grasp the significance of const Vec::new, could you please elaborate on that. Furthermore, I did try Eq and Hash using #[derive()] but they are not implemented for f32 :(. I will also try to reduce the magic numbers. Thanks for the feedback! \$\endgroup\$
    – Saphereye
    Apr 3 at 7:59
  • \$\begingroup\$ Sorry, you're right. f32 does not implement neither Eq not Hash, so #[derive()] is indeed not possible. \$\endgroup\$ Apr 8 at 19:34

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.