salva2d/
z_order.rs

1use crate::math::{Point, Real};
2use num_traits::float::FloatCore;
3use std::cmp::Ordering;
4
5pub fn apply_permutation<T: Clone>(permutation: &[usize], data: &[T]) -> Vec<T> {
6    permutation.iter().map(|i| data[*i].clone()).collect()
7}
8
9pub fn compute_points_z_order(points: &[Point<Real>]) -> Vec<usize> {
10    let mut indices: Vec<_> = (0..points.len()).collect();
11    indices.sort_unstable_by(|i, j| {
12        z_order_floats(points[*i].coords.as_slice(), points[*j].coords.as_slice())
13            .unwrap_or(std::cmp::Ordering::Equal)
14    });
15    indices
16}
17
18// Fast construction of k-Nearest Neighbor Graphs for Point Clouds
19// Michael Connor, Piyush Kumar
20// Algorithm 1
21//
22// http://compgeom.com/~piyush/papers/tvcg_stann.pdf
23pub fn z_order_floats(p1: &[Real], p2: &[Real]) -> Option<Ordering> {
24    assert_eq!(
25        p1.len(),
26        p2.len(),
27        "Cannot compare array with different lengths."
28    );
29    let mut x = 0;
30    let mut dim = 0;
31
32    for j in 0..p1.len() {
33        let y = xor_msb_float(p1[j], p2[j]);
34        if x < y {
35            x = y;
36            dim = j;
37        }
38    }
39
40    p1[dim].partial_cmp(&p2[dim])
41}
42
43fn xor_msb_float(a: Real, b: Real) -> i16 {
44    let fa = na::try_convert::<_, f64>(a).unwrap();
45    let fb = na::try_convert::<_, f64>(b).unwrap();
46    let (mantissa1, exponent1, _sign1) = fa.integer_decode();
47    let (mantissa2, exponent2, _sign2) = fb.integer_decode();
48    let x = exponent1; // To keep the same notation as the paper.
49    let y = exponent2; // To keep the same notation as the paper.
50
51    if x == y {
52        let z = msdb(mantissa1, mantissa2);
53        x - z
54    } else if y < x {
55        x
56    } else {
57        y
58    }
59}
60
61fn msdb(x: u64, y: u64) -> i16 {
62    64i16 - (x ^ y).leading_zeros() as i16
63}