# 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]

[[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.

• I'd suggest looking into SIMD and Google's SwissTable (HashBrown on crates.io). Apr 2 at 15:57
• Vec2::new could be const and you could just #[derive()] Eq and Hash for it. Apr 2 at 20:11
• Also there are a lot of magic numbers in your code. Apr 2 at 20:16
• 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! Apr 3 at 7:59
• Sorry, you're right. f32 does not implement neither Eq not Hash, so #[derive()] is indeed not possible. Apr 8 at 19:34