Skip to content

Commit

Permalink
generic sort trait
Browse files Browse the repository at this point in the history
  • Loading branch information
kylebarron committed Jan 2, 2024
1 parent 5e5cdbf commit a5265e1
Show file tree
Hide file tree
Showing 7 changed files with 210 additions and 160 deletions.
3 changes: 2 additions & 1 deletion benches/flatbush.rs
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
use bytemuck::cast_slice;
use criterion::{criterion_group, criterion_main, Criterion};
use flatbush::flatbush::HilbertSort;
use flatbush::{FlatbushBuilder, FlatbushIndex, OwnedFlatbush};
use rstar::primitives::{GeomWithData, Rectangle};
use rstar::{RTree, AABB};
Expand All @@ -19,7 +20,7 @@ fn construct_flatbush(boxes_buf: &[f64]) -> OwnedFlatbush {
let max_y = box_[3];
builder.add(min_x, min_y, max_x, max_y);
}
builder.finish()
builder.finish::<HilbertSort>()
}

fn construct_rstar(
Expand Down
166 changes: 11 additions & 155 deletions src/flatbush/builder.rs
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ use bytemuck::cast_slice_mut;

use crate::flatbush::constants::VERSION;
use crate::flatbush::index::OwnedFlatbush;
use crate::flatbush::sort::{Sort, SortParams};
use crate::flatbush::util::compute_num_nodes;
use crate::indices::MutableIndices;

Expand Down Expand Up @@ -109,7 +110,7 @@ impl FlatbushBuilder {
index
}

pub fn finish(mut self) -> OwnedFlatbush {
pub fn finish<S: Sort>(mut self) -> OwnedFlatbush {
assert_eq!(
self.pos >> 2,
self.num_items,
Expand Down Expand Up @@ -145,39 +146,15 @@ impl FlatbushBuilder {
};
}

let width = self.max_x - self.min_x; // || 1.0;
let height = self.max_y - self.min_y; // || 1.0;
let mut hilbert_values: Vec<u32> = Vec::with_capacity(self.num_items);
let hilbert_max = ((1 << 16) - 1) as f64;

{
// map item centers into Hilbert coordinate space and calculate Hilbert values
let mut pos = 0;
for _ in 0..self.num_items {
let min_x = boxes[pos];
pos += 1;
let min_y = boxes[pos];
pos += 1;
let max_x = boxes[pos];
pos += 1;
let max_y = boxes[pos];
pos += 1;

let x = (hilbert_max * ((min_x + max_x) / 2. - self.min_x) / width).floor() as u32;
let y = (hilbert_max * ((min_y + max_y) / 2. - self.min_y) / height).floor() as u32;
hilbert_values.push(hilbert(x, y));
}
}

// sort items by their Hilbert value (for packing later)
sort(
&mut hilbert_values,
boxes,
&mut indices,
0,
self.num_items - 1,
self.node_size,
);
let mut sort_params = SortParams {
num_items: self.num_items,
node_size: self.node_size,
min_x: self.min_x,
min_y: self.min_y,
max_x: self.max_x,
max_y: self.max_y,
};
S::sort(&mut sort_params, boxes, &mut indices);

{
// generate nodes at each tree level, bottom-up
Expand Down Expand Up @@ -248,124 +225,3 @@ fn split_data_borrow(
let indices = MutableIndices::new(indices_buf, num_nodes);
(boxes, indices)
}

/// Custom quicksort that partially sorts bbox data alongside the hilbert values.
// Partially taken from static_aabb2d_index under the MIT/Apache license
fn sort(
values: &mut [u32],
boxes: &mut [f64],
indices: &mut MutableIndices,
left: usize,
right: usize,
node_size: usize,
) {
debug_assert!(left <= right);

if left / node_size >= right / node_size {
return;
}

let midpoint = (left + right) / 2;
let pivot = values[midpoint];
let mut i = left.wrapping_sub(1);
let mut j = right.wrapping_add(1);

loop {
loop {
i = i.wrapping_add(1);
if values[i] >= pivot {
break;
}
}

loop {
j = j.wrapping_sub(1);
if values[j] <= pivot {
break;
}
}

if i >= j {
break;
}

swap(values, boxes, indices, i, j);
}

sort(values, boxes, indices, left, j, node_size);
sort(values, boxes, indices, j.wrapping_add(1), right, node_size);
}

/// Swap two values and two corresponding boxes.
#[inline]
fn swap(values: &mut [u32], boxes: &mut [f64], indices: &mut MutableIndices, i: usize, j: usize) {
values.swap(i, j);

let k = 4 * i;
let m = 4 * j;
boxes.swap(k, m);
boxes.swap(k + 1, m + 1);
boxes.swap(k + 2, m + 2);
boxes.swap(k + 3, m + 3);

indices.swap(i, j);
}

// Taken from static_aabb2d_index under the mit/apache license
// https://github.com/jbuckmccready/static_aabb2d_index/blob/9e6add59d77b74d4de0ac32159db47fbcb3acc28/src/static_aabb2d_index.rs#L486C1-L544C2
fn hilbert(x: u32, y: u32) -> u32 {
// Fast Hilbert curve algorithm by http://threadlocalmutex.com/
// Ported from C++ https://github.com/rawrunprotected/hilbert_curves (public domain)
let mut a_1 = x ^ y;
let mut b_1 = 0xFFFF ^ a_1;
let mut c_1 = 0xFFFF ^ (x | y);
let mut d_1 = x & (y ^ 0xFFFF);

let mut a_2 = a_1 | (b_1 >> 1);
let mut b_2 = (a_1 >> 1) ^ a_1;
let mut c_2 = ((c_1 >> 1) ^ (b_1 & (d_1 >> 1))) ^ c_1;
let mut d_2 = ((a_1 & (c_1 >> 1)) ^ (d_1 >> 1)) ^ d_1;

a_1 = a_2;
b_1 = b_2;
c_1 = c_2;
d_1 = d_2;
a_2 = (a_1 & (a_1 >> 2)) ^ (b_1 & (b_1 >> 2));
b_2 = (a_1 & (b_1 >> 2)) ^ (b_1 & ((a_1 ^ b_1) >> 2));
c_2 ^= (a_1 & (c_1 >> 2)) ^ (b_1 & (d_1 >> 2));
d_2 ^= (b_1 & (c_1 >> 2)) ^ ((a_1 ^ b_1) & (d_1 >> 2));

a_1 = a_2;
b_1 = b_2;
c_1 = c_2;
d_1 = d_2;
a_2 = (a_1 & (a_1 >> 4)) ^ (b_1 & (b_1 >> 4));
b_2 = (a_1 & (b_1 >> 4)) ^ (b_1 & ((a_1 ^ b_1) >> 4));
c_2 ^= (a_1 & (c_1 >> 4)) ^ (b_1 & (d_1 >> 4));
d_2 ^= (b_1 & (c_1 >> 4)) ^ ((a_1 ^ b_1) & (d_1 >> 4));

a_1 = a_2;
b_1 = b_2;
c_1 = c_2;
d_1 = d_2;
c_2 ^= (a_1 & (c_1 >> 8)) ^ (b_1 & (d_1 >> 8));
d_2 ^= (b_1 & (c_1 >> 8)) ^ ((a_1 ^ b_1) & (d_1 >> 8));

a_1 = c_2 ^ (c_2 >> 1);
b_1 = d_2 ^ (d_2 >> 1);

let mut i0 = x ^ y;
let mut i1 = b_1 | (0xFFFF ^ (i0 | a_1));

i0 = (i0 | (i0 << 8)) & 0x00FF00FF;
i0 = (i0 | (i0 << 4)) & 0x0F0F0F0F;
i0 = (i0 | (i0 << 2)) & 0x33333333;
i0 = (i0 | (i0 << 1)) & 0x55555555;

i1 = (i1 | (i1 << 8)) & 0x00FF00FF;
i1 = (i1 | (i1 << 4)) & 0x0F0F0F0F;
i1 = (i1 | (i1 << 2)) & 0x33333333;
i1 = (i1 | (i1 << 1)) & 0x55555555;

(i1 << 1) | i0
}
2 changes: 2 additions & 0 deletions src/flatbush/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,11 @@ pub mod builder;
pub mod constants;
pub mod error;
pub mod index;
pub mod sort;
pub mod r#trait;
pub mod util;

pub use builder::FlatbushBuilder;
pub use index::{FlatbushRef, OwnedFlatbush};
pub use r#trait::FlatbushIndex;
pub use sort::HilbertSort;
167 changes: 167 additions & 0 deletions src/flatbush/sort/hilbert.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,167 @@
use crate::flatbush::sort::{Sort, SortParams};
use crate::indices::MutableIndices;

#[derive(Debug, Clone, Copy)]
pub struct HilbertSort;

impl Sort for HilbertSort {
fn sort(params: &mut SortParams, boxes: &mut [f64], indices: &mut MutableIndices) {
let width = params.max_x - params.min_x; // || 1.0;
let height = params.max_y - params.min_y; // || 1.0;
let mut hilbert_values: Vec<u32> = Vec::with_capacity(params.num_items);
let hilbert_max = ((1 << 16) - 1) as f64;

{
// map item centers into Hilbert coordinate space and calculate Hilbert values
let mut pos = 0;
for _ in 0..params.num_items {
let min_x = boxes[pos];
pos += 1;
let min_y = boxes[pos];
pos += 1;
let max_x = boxes[pos];
pos += 1;
let max_y = boxes[pos];
pos += 1;

let x =
(hilbert_max * ((min_x + max_x) / 2. - params.min_x) / width).floor() as u32;
let y =
(hilbert_max * ((min_y + max_y) / 2. - params.min_y) / height).floor() as u32;
hilbert_values.push(hilbert(x, y));
}
}

// sort items by their Hilbert value (for packing later)
sort(
&mut hilbert_values,
boxes,
indices,
0,
params.num_items - 1,
params.node_size,
);
}
}

/// Custom quicksort that partially sorts bbox data alongside the hilbert values.
// Partially taken from static_aabb2d_index under the MIT/Apache license
fn sort(
values: &mut [u32],
boxes: &mut [f64],
indices: &mut MutableIndices,
left: usize,
right: usize,
node_size: usize,
) {
debug_assert!(left <= right);

if left / node_size >= right / node_size {
return;
}

let midpoint = (left + right) / 2;
let pivot = values[midpoint];
let mut i = left.wrapping_sub(1);
let mut j = right.wrapping_add(1);

loop {
loop {
i = i.wrapping_add(1);
if values[i] >= pivot {
break;
}
}

loop {
j = j.wrapping_sub(1);
if values[j] <= pivot {
break;
}
}

if i >= j {
break;
}

swap(values, boxes, indices, i, j);
}

sort(values, boxes, indices, left, j, node_size);
sort(values, boxes, indices, j.wrapping_add(1), right, node_size);
}

/// Swap two values and two corresponding boxes.
#[inline]
fn swap(values: &mut [u32], boxes: &mut [f64], indices: &mut MutableIndices, i: usize, j: usize) {
values.swap(i, j);

let k = 4 * i;
let m = 4 * j;
boxes.swap(k, m);
boxes.swap(k + 1, m + 1);
boxes.swap(k + 2, m + 2);
boxes.swap(k + 3, m + 3);

indices.swap(i, j);
}

// Taken from static_aabb2d_index under the mit/apache license
// https://github.com/jbuckmccready/static_aabb2d_index/blob/9e6add59d77b74d4de0ac32159db47fbcb3acc28/src/static_aabb2d_index.rs#L486C1-L544C2
#[inline]
fn hilbert(x: u32, y: u32) -> u32 {
// Fast Hilbert curve algorithm by http://threadlocalmutex.com/
// Ported from C++ https://github.com/rawrunprotected/hilbert_curves (public domain)
let mut a_1 = x ^ y;
let mut b_1 = 0xFFFF ^ a_1;
let mut c_1 = 0xFFFF ^ (x | y);
let mut d_1 = x & (y ^ 0xFFFF);

let mut a_2 = a_1 | (b_1 >> 1);
let mut b_2 = (a_1 >> 1) ^ a_1;
let mut c_2 = ((c_1 >> 1) ^ (b_1 & (d_1 >> 1))) ^ c_1;
let mut d_2 = ((a_1 & (c_1 >> 1)) ^ (d_1 >> 1)) ^ d_1;

a_1 = a_2;
b_1 = b_2;
c_1 = c_2;
d_1 = d_2;
a_2 = (a_1 & (a_1 >> 2)) ^ (b_1 & (b_1 >> 2));
b_2 = (a_1 & (b_1 >> 2)) ^ (b_1 & ((a_1 ^ b_1) >> 2));
c_2 ^= (a_1 & (c_1 >> 2)) ^ (b_1 & (d_1 >> 2));
d_2 ^= (b_1 & (c_1 >> 2)) ^ ((a_1 ^ b_1) & (d_1 >> 2));

a_1 = a_2;
b_1 = b_2;
c_1 = c_2;
d_1 = d_2;
a_2 = (a_1 & (a_1 >> 4)) ^ (b_1 & (b_1 >> 4));
b_2 = (a_1 & (b_1 >> 4)) ^ (b_1 & ((a_1 ^ b_1) >> 4));
c_2 ^= (a_1 & (c_1 >> 4)) ^ (b_1 & (d_1 >> 4));
d_2 ^= (b_1 & (c_1 >> 4)) ^ ((a_1 ^ b_1) & (d_1 >> 4));

a_1 = a_2;
b_1 = b_2;
c_1 = c_2;
d_1 = d_2;
c_2 ^= (a_1 & (c_1 >> 8)) ^ (b_1 & (d_1 >> 8));
d_2 ^= (b_1 & (c_1 >> 8)) ^ ((a_1 ^ b_1) & (d_1 >> 8));

a_1 = c_2 ^ (c_2 >> 1);
b_1 = d_2 ^ (d_2 >> 1);

let mut i0 = x ^ y;
let mut i1 = b_1 | (0xFFFF ^ (i0 | a_1));

i0 = (i0 | (i0 << 8)) & 0x00FF00FF;
i0 = (i0 | (i0 << 4)) & 0x0F0F0F0F;
i0 = (i0 | (i0 << 2)) & 0x33333333;
i0 = (i0 | (i0 << 1)) & 0x55555555;

i1 = (i1 | (i1 << 8)) & 0x00FF00FF;
i1 = (i1 | (i1 << 4)) & 0x0F0F0F0F;
i1 = (i1 | (i1 << 2)) & 0x33333333;
i1 = (i1 | (i1 << 1)) & 0x55555555;

(i1 << 1) | i0
}
5 changes: 5 additions & 0 deletions src/flatbush/sort/mod.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
pub mod hilbert;
pub mod r#trait;

pub use hilbert::HilbertSort;
pub use r#trait::{Sort, SortParams};
Loading

0 comments on commit a5265e1

Please sign in to comment.