summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/color.rs308
-rw-r--r--src/color/order.rs196
-rw-r--r--src/color/source.rs71
-rw-r--r--src/forest.rs292
-rw-r--r--src/frontier.rs194
-rw-r--r--src/frontier/image.rs79
-rw-r--r--src/frontier/mean.rs145
-rw-r--r--src/frontier/min.rs153
-rw-r--r--src/hilbert.rs136
-rw-r--r--src/main.rs492
-rw-r--r--src/soft.rs289
11 files changed, 2355 insertions, 0 deletions
diff --git a/src/color.rs b/src/color.rs
new file mode 100644
index 0000000..ff113a7
--- /dev/null
+++ b/src/color.rs
@@ -0,0 +1,308 @@
+//! Colors and color spaces.
+
+pub mod order;
+pub mod source;
+
+use acap::coords::{Coordinates, CoordinateMetric, CoordinateProximity};
+use acap::distance::{Metric, Proximity};
+use acap::euclid::{EuclideanDistance, euclidean_distance};
+
+use image::Rgb;
+
+use std::ops::Index;
+
+/// An 8-bit RGB color.
+pub type Rgb8 = Rgb<u8>;
+
+/// A [color space](https://en.wikipedia.org/wiki/Color_space).
+pub trait ColorSpace: Copy + From<Rgb8>
+ + Coordinates
+ + Metric
+ + CoordinateMetric<<Self as Coordinates>::Value, Distance = <Self as Proximity>::Distance>
+{
+ /// Compute the average of the given colors.
+ fn average<I: IntoIterator<Item = Self>>(colors: I) -> Self;
+}
+
+/// [sRGB](https://en.wikipedia.org/wiki/SRGB) space.
+#[derive(Clone, Copy, Debug)]
+pub struct RgbSpace([f64; 3]);
+
+impl Index<usize> for RgbSpace {
+ type Output = f64;
+
+ fn index(&self, i: usize) -> &f64 {
+ &self.0[i]
+ }
+}
+
+impl From<Rgb8> for RgbSpace {
+ fn from(rgb8: Rgb8) -> Self {
+ Self([
+ (rgb8[0] as f64) / 255.0,
+ (rgb8[1] as f64) / 255.0,
+ (rgb8[2] as f64) / 255.0,
+ ])
+ }
+}
+
+impl Coordinates for RgbSpace {
+ type Value = f64;
+
+ fn dims(&self) -> usize {
+ self.0.dims()
+ }
+
+ fn coord(&self, i: usize) -> f64 {
+ self.0.coord(i)
+ }
+}
+
+impl Proximity for RgbSpace {
+ type Distance = EuclideanDistance<f64>;
+
+ fn distance(&self, other: &Self) -> Self::Distance {
+ euclidean_distance(&self.0, &other.0)
+ }
+}
+
+impl Metric for RgbSpace {}
+
+impl CoordinateProximity<f64> for RgbSpace {
+ type Distance = EuclideanDistance<f64>;
+
+ fn distance_to_coords(&self, other: &[f64]) -> Self::Distance {
+ euclidean_distance(&self.0, other)
+ }
+}
+
+impl CoordinateMetric<f64> for RgbSpace {}
+
+impl ColorSpace for RgbSpace {
+ fn average<I: IntoIterator<Item = Self>>(colors: I) -> Self {
+ let mut sum = [0.0, 0.0, 0.0];
+ let mut len: usize = 0;
+ for color in colors.into_iter() {
+ for i in 0..3 {
+ sum[i] += color[i];
+ }
+ len += 1;
+ }
+ for s in &mut sum {
+ *s /= len as f64;
+ }
+ Self(sum)
+ }
+}
+
+/// [CIE XYZ](https://en.wikipedia.org/wiki/CIE_1931_color_space) space.
+#[derive(Clone, Copy, Debug)]
+struct XyzSpace([f64; 3]);
+
+impl Index<usize> for XyzSpace {
+ type Output = f64;
+
+ fn index(&self, i: usize) -> &f64 {
+ &self.0[i]
+ }
+}
+
+/// The inverse of the sRGB gamma function.
+fn srgb_inv_gamma(t: f64) -> f64 {
+ if t <= 0.040449936 {
+ t / 12.92
+ } else {
+ ((t + 0.055) / 1.055).powf(2.4)
+ }
+}
+
+impl From<Rgb8> for XyzSpace {
+ fn from(rgb8: Rgb8) -> Self {
+ let rgb = RgbSpace::from(rgb8);
+
+ let r = srgb_inv_gamma(rgb[0]);
+ let g = srgb_inv_gamma(rgb[1]);
+ let b = srgb_inv_gamma(rgb[2]);
+
+ Self([
+ 0.4123808838268995 * r + 0.3575728355732478 * g + 0.1804522977447919 * b,
+ 0.2126198631048975 * r + 0.7151387878413206 * g + 0.0721499433963131 * b,
+ 0.0193434956789248 * r + 0.1192121694056356 * g + 0.9505065664127130 * b,
+ ])
+ }
+}
+
+/// CIE D50 [white point](https://en.wikipedia.org/wiki/Standard_illuminant).
+const WHITE: XyzSpace = XyzSpace([0.9504060171449392, 0.9999085943425312, 1.089062231497274]);
+
+/// CIE L\*a\*b\* (and L\*u\*v\*) gamma
+fn lab_gamma(t: f64) -> f64 {
+ if t > 216.0 / 24389.0 {
+ t.cbrt()
+ } else {
+ 841.0 * t / 108.0 + 4.0 / 29.0
+ }
+}
+
+/// [CIE L\*a\*b\*](https://en.wikipedia.org/wiki/CIELAB_color_space) space.
+#[derive(Clone, Copy, Debug)]
+pub struct LabSpace([f64; 3]);
+
+impl Index<usize> for LabSpace {
+ type Output = f64;
+
+ fn index(&self, i: usize) -> &f64 {
+ &self.0[i]
+ }
+}
+
+impl From<Rgb8> for LabSpace {
+ fn from(rgb8: Rgb8) -> Self {
+ let xyz = XyzSpace::from(rgb8);
+
+ let x = lab_gamma(xyz[0] / WHITE[0]);
+ let y = lab_gamma(xyz[1] / WHITE[1]);
+ let z = lab_gamma(xyz[2] / WHITE[2]);
+
+ let l = 116.0 * y - 16.0;
+ let a = 500.0 * (x - y);
+ let b = 200.0 * (y - z);
+
+ Self([l, a, b])
+ }
+}
+
+impl Coordinates for LabSpace {
+ type Value = f64;
+
+ fn dims(&self) -> usize {
+ self.0.dims()
+ }
+
+ fn coord(&self, i: usize) -> f64 {
+ self.0.coord(i)
+ }
+}
+
+impl Proximity for LabSpace {
+ type Distance = EuclideanDistance<f64>;
+
+ fn distance(&self, other: &Self) -> Self::Distance {
+ euclidean_distance(self.0, other.0)
+ }
+}
+
+impl Metric for LabSpace {}
+
+impl CoordinateProximity<f64> for LabSpace {
+ type Distance = EuclideanDistance<f64>;
+
+ fn distance_to_coords(&self, other: &[f64]) -> Self::Distance {
+ euclidean_distance(&self.0, other)
+ }
+}
+
+impl CoordinateMetric<f64> for LabSpace {}
+
+impl ColorSpace for LabSpace {
+ fn average<I: IntoIterator<Item = Self>>(colors: I) -> Self {
+ let mut sum = [0.0, 0.0, 0.0];
+ let mut len: usize = 0;
+ for color in colors.into_iter() {
+ for i in 0..3 {
+ sum[i] += color[i];
+ }
+ len += 1;
+ }
+ for s in &mut sum {
+ *s /= len as f64;
+ }
+ Self(sum)
+ }
+}
+
+/// [CIE L\*u\*v\*](https://en.wikipedia.org/wiki/CIELUV) space.
+#[derive(Clone, Copy, Debug)]
+pub struct LuvSpace([f64; 3]);
+
+impl Index<usize> for LuvSpace {
+ type Output = f64;
+
+ fn index(&self, i: usize) -> &f64 {
+ &self.0[i]
+ }
+}
+
+/// Computes the u' and v' values for L\*u\*v\*.
+fn uv_prime(xyz: &XyzSpace) -> (f64, f64) {
+ let denom = xyz[0] + 15.0 * xyz[1] + 3.0 * xyz[2];
+ if denom == 0.0 {
+ (0.0, 0.0)
+ } else {
+ (4.0 * xyz[0] / denom, 9.0 * xyz[1] / denom)
+ }
+}
+
+impl From<Rgb8> for LuvSpace {
+ fn from(rgb8: Rgb8) -> Self {
+ let xyz = XyzSpace::from(rgb8);
+
+ let (uprime, vprime) = uv_prime(&xyz);
+ let (unprime, vnprime) = uv_prime(&WHITE);
+
+ let l = 116.0 * lab_gamma(xyz[1] / WHITE[1]) - 16.0;
+ let u = 13.0 * l * (uprime - unprime);
+ let v = 13.0 * l * (vprime - vnprime);
+
+ Self([l, u, v])
+ }
+}
+
+impl Coordinates for LuvSpace {
+ type Value = f64;
+
+ fn dims(&self) -> usize {
+ self.0.dims()
+ }
+
+ fn coord(&self, i: usize) -> f64 {
+ self.0.coord(i)
+ }
+}
+
+impl Proximity for LuvSpace {
+ type Distance = EuclideanDistance<f64>;
+
+ fn distance(&self, other: &Self) -> Self::Distance {
+ euclidean_distance(&self.0, &other.0)
+ }
+}
+
+impl Metric for LuvSpace {}
+
+impl CoordinateProximity<f64> for LuvSpace {
+ type Distance = EuclideanDistance<f64>;
+
+ fn distance_to_coords(&self, other: &[f64]) -> Self::Distance {
+ euclidean_distance(&self.0, other)
+ }
+}
+
+impl CoordinateMetric<f64> for LuvSpace {}
+
+impl ColorSpace for LuvSpace {
+ fn average<I: IntoIterator<Item = Self>>(colors: I) -> Self {
+ let mut sum = [0.0, 0.0, 0.0];
+ let mut len: usize = 0;
+ for color in colors.into_iter() {
+ for i in 0..3 {
+ sum[i] += color[i];
+ }
+ len += 1;
+ }
+ for s in &mut sum {
+ *s /= len as f64;
+ }
+ Self(sum)
+ }
+}
diff --git a/src/color/order.rs b/src/color/order.rs
new file mode 100644
index 0000000..300a556
--- /dev/null
+++ b/src/color/order.rs
@@ -0,0 +1,196 @@
+//! Linear orders for colors.
+
+use super::source::ColorSource;
+use super::Rgb8;
+
+use crate::hilbert::hilbert_point;
+
+use rand::seq::SliceRandom;
+use rand::Rng;
+
+use std::cmp::Ordering;
+
+/// An iterator over all colors from a source.
+#[derive(Debug)]
+struct ColorSourceIter<S> {
+ source: S,
+ coords: Vec<usize>,
+}
+
+impl<S: ColorSource> From<S> for ColorSourceIter<S> {
+ fn from(source: S) -> Self {
+ let coords = vec![0; source.dimensions().len()];
+
+ Self { source, coords }
+ }
+}
+
+impl<S: ColorSource> Iterator for ColorSourceIter<S> {
+ type Item = Rgb8;
+
+ fn next(&mut self) -> Option<Rgb8> {
+ if self.coords.is_empty() {
+ return None;
+ }
+
+ let color = self.source.get_color(&self.coords);
+
+ let dims = self.source.dimensions();
+ for i in 0..dims.len() {
+ self.coords[i] += 1;
+ if self.coords[i] < dims[i] {
+ break;
+ } else if i == dims.len() - 1 {
+ self.coords.clear();
+ } else {
+ self.coords[i] = 0;
+ }
+ }
+
+ Some(color)
+ }
+}
+
+/// Wrapper for sorting colors by hue.
+#[derive(Debug, Eq, PartialEq)]
+struct Hue {
+ /// The quadrant of the hue angle.
+ quad: i32,
+ /// The numerator of the hue calculation.
+ num: i32,
+ /// The denominator of the hue calculation.
+ denom: i32,
+}
+
+impl From<Rgb8> for Hue {
+ fn from(rgb8: Rgb8) -> Self {
+ // The hue angle is atan2(sqrt(3) * (G - B), 2 * R - G - B). We avoid actually computing
+ // the atan2() as an optimization.
+ let r = rgb8[0] as i32;
+ let g = rgb8[1] as i32;
+ let b = rgb8[2] as i32;
+
+ let num = g - b;
+ let mut denom = 2 * r - g - b;
+ if num == 0 && denom == 0 {
+ denom = 1;
+ }
+
+ let quad = match (num >= 0, denom >= 0) {
+ (true, true) => 0,
+ (true, false) => 1,
+ (false, false) => 2,
+ (false, true) => 3,
+ };
+
+ Self { quad, num, denom }
+ }
+}
+
+impl Ord for Hue {
+ fn cmp(&self, other: &Self) -> Ordering {
+ // Within the same quadrant,
+ //
+ // atan2(n1, d1) < atan2(n2, d2) iff
+ // n1 / d1 < n2 / d2 iff
+ // n1 * d2 < n2 * d1
+ self.quad
+ .cmp(&other.quad)
+ .then_with(|| (self.num * other.denom).cmp(&(other.num * self.denom)))
+ }
+}
+
+impl PartialOrd for Hue {
+ fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
+ Some(self.cmp(other))
+ }
+}
+
+/// Iterate over colors sorted by their hue.
+pub fn hue_sorted<S: ColorSource>(source: S) -> Vec<Rgb8> {
+ let mut colors: Vec<_> = ColorSourceIter::from(source).collect();
+ colors.sort_by_key(|c| Hue::from(*c));
+ colors
+}
+
+/// Iterate over colors in random order.
+pub fn shuffled<S: ColorSource, R: Rng>(source: S, rng: &mut R) -> Vec<Rgb8> {
+ let mut colors: Vec<_> = ColorSourceIter::from(source).collect();
+ colors.shuffle(rng);
+ colors
+}
+
+/// ceil(log_2(n)). for rounding up to powers of 2.
+fn log2(n: usize) -> u32 {
+ let nbits = 8 * std::mem::size_of::<usize>() as u32;
+ nbits - (n - 1).leading_zeros()
+}
+
+/// Iterate over colors in Morton order (Z-order).
+pub fn morton<S: ColorSource>(source: S) -> Vec<Rgb8> {
+ let mut colors = Vec::new();
+
+ let dims = source.dimensions();
+ let ndims = dims.len();
+
+ let nbits = ndims * dims.iter().map(|n| log2(*n) as usize).max().unwrap();
+
+ let size = 1usize << nbits;
+ let mut coords = vec![0; ndims];
+ for i in 0..size {
+ for x in &mut coords {
+ *x = 0;
+ }
+ for j in 0..nbits {
+ let bit = (i >> j) & 1;
+ coords[j % ndims] |= bit << (j / ndims);
+ }
+ if coords.iter().zip(dims.iter()).all(|(x, n)| x < n) {
+ colors.push(source.get_color(&coords));
+ }
+ }
+
+ colors
+}
+
+/// Iterate over colors in Hilbert curve order.
+pub fn hilbert<S: ColorSource>(source: S) -> Vec<Rgb8> {
+ let mut colors = Vec::new();
+
+ let dims = source.dimensions();
+ let ndims = dims.len();
+
+ let bits: Vec<_> = dims.iter().map(|n| log2(*n)).collect();
+ let nbits: u32 = bits.iter().sum();
+ let size = 1usize << nbits;
+
+ let mut coords = vec![0; ndims];
+
+ for i in 0..size {
+ hilbert_point(i, &bits, &mut coords);
+ if coords.iter().zip(dims.iter()).all(|(x, n)| x < n) {
+ colors.push(source.get_color(&coords));
+ }
+ }
+
+ colors
+}
+
+/// Stripe an ordered list of colors, to reduce artifacts in the generated image.
+///
+/// The striped ordering gives every other item first, then every other item from the remaining
+/// items, etc. For example, the striped form of `0..16` is
+/// `[0, 2, 4, 6, 8, 10, 12, 14, 1, 5, 9, 13, 3, 11, 7, 15]`.
+pub fn striped(colors: Vec<Rgb8>) -> Vec<Rgb8> {
+ let len = colors.len();
+ let mut result = Vec::with_capacity(len);
+ let mut stripe = 1;
+ while stripe <= len {
+ for i in ((stripe - 1)..len).step_by(2 * stripe) {
+ result.push(colors[i]);
+ }
+ stripe *= 2;
+ }
+
+ result
+}
diff --git a/src/color/source.rs b/src/color/source.rs
new file mode 100644
index 0000000..5cc9631
--- /dev/null
+++ b/src/color/source.rs
@@ -0,0 +1,71 @@
+//! Sources of colors.
+
+use super::Rgb8;
+
+use image::RgbImage;
+
+/// A source of colors in multidimensional space.
+pub trait ColorSource {
+ /// Get the size of each dimension in this space.
+ fn dimensions(&self) -> &[usize];
+
+ /// Get the color at some particular coordinates.
+ fn get_color(&self, coords: &[usize]) -> Rgb8;
+}
+
+/// The entire RGB space.
+#[derive(Debug)]
+pub struct AllColors {
+ dims: [usize; 3],
+ shifts: [u32; 3],
+}
+
+impl AllColors {
+ /// Create an AllColors source with the given bit depths.
+ pub fn new(r: u32, g: u32, b: u32) -> Self {
+ Self {
+ dims: [1 << r, 1 << g, 1 << b],
+ shifts: [8 - r, 8 - g, 8 - b],
+ }
+ }
+}
+
+impl ColorSource for AllColors {
+ fn dimensions(&self) -> &[usize] {
+ &self.dims
+ }
+
+ fn get_color(&self, coords: &[usize]) -> Rgb8 {
+ Rgb8::from([
+ (coords[0] << self.shifts[0]) as u8,
+ (coords[1] << self.shifts[1]) as u8,
+ (coords[2] << self.shifts[2]) as u8,
+ ])
+ }
+}
+
+/// Colors extracted from an image.
+#[derive(Debug)]
+pub struct ImageColors {
+ dims: [usize; 2],
+ image: RgbImage,
+}
+
+impl From<RgbImage> for ImageColors {
+ fn from(image: RgbImage) -> Self {
+ Self {
+ dims: [image.width() as usize, image.height() as usize],
+ image,
+ }
+ }
+}
+
+impl ColorSource for ImageColors {
+ fn dimensions(&self) -> &[usize] {
+ &self.dims
+ }
+
+ fn get_color(&self, coords: &[usize]) -> Rgb8 {
+ *self.image.get_pixel(coords[0] as u32, coords[1] as u32)
+ }
+}
diff --git a/src/forest.rs b/src/forest.rs
new file mode 100644
index 0000000..56dff7e
--- /dev/null
+++ b/src/forest.rs
@@ -0,0 +1,292 @@
+//! [Dynamization](https://en.wikipedia.org/wiki/Dynamization) for nearest neighbor search.
+
+use acap::distance::Proximity;
+use acap::kd::FlatKdTree;
+use acap::vp::FlatVpTree;
+use acap::{NearestNeighbors, Neighborhood};
+
+use std::iter::{self, Extend, FromIterator};
+
+/// The number of bits dedicated to the flat buffer.
+const BUFFER_BITS: usize = 6;
+/// The maximum size of the buffer.
+const BUFFER_SIZE: usize = 1 << BUFFER_BITS;
+
+/// A dynamic wrapper for a static nearest neighbor search data structure.
+///
+/// This type applies [dynamization](https://en.wikipedia.org/wiki/Dynamization) to an arbitrary
+/// nearest neighbor search structure `T`, allowing new items to be added dynamically.
+#[derive(Debug)]
+pub struct Forest<T: IntoIterator> {
+ /// A flat buffer used for the first few items, to avoid repeatedly rebuilding small trees.
+ buffer: Vec<T::Item>,
+ /// The trees of the forest, with sizes in geometric progression.
+ trees: Vec<Option<T>>,
+}
+
+impl<T, U> Forest<U>
+where
+ U: FromIterator<T> + IntoIterator<Item = T>,
+{
+ /// Create a new empty forest.
+ pub fn new() -> Self {
+ Self {
+ buffer: Vec::new(),
+ trees: Vec::new(),
+ }
+ }
+
+ /// Add a new item to the forest.
+ pub fn push(&mut self, item: T) {
+ self.extend(iter::once(item));
+ }
+
+ /// Get the number of items in the forest.
+ pub fn len(&self) -> usize {
+ let mut len = self.buffer.len();
+ for (i, slot) in self.trees.iter().enumerate() {
+ if slot.is_some() {
+ len += 1 << (i + BUFFER_BITS);
+ }
+ }
+ len
+ }
+
+ /// Check if this forest is empty.
+ pub fn is_empty(&self) -> bool {
+ if !self.buffer.is_empty() {
+ return false;
+ }
+
+ self.trees.iter().flatten().next().is_none()
+ }
+}
+
+impl<T, U> Default for Forest<U>
+where
+ U: FromIterator<T> + IntoIterator<Item = T>,
+{
+ fn default() -> Self {
+ Self::new()
+ }
+}
+
+impl<T, U> Extend<T> for Forest<U>
+where
+ U: FromIterator<T> + IntoIterator<Item = T>,
+{
+ fn extend<I: IntoIterator<Item = T>>(&mut self, items: I) {
+ self.buffer.extend(items);
+ if self.buffer.len() < BUFFER_SIZE {
+ return;
+ }
+
+ let len = self.len();
+
+ for i in 0.. {
+ let bit = 1 << (i + BUFFER_BITS);
+
+ if bit > len {
+ break;
+ }
+
+ if i >= self.trees.len() {
+ self.trees.push(None);
+ }
+
+ if len & bit == 0 {
+ if let Some(tree) = self.trees[i].take() {
+ self.buffer.extend(tree);
+ }
+ } else if self.trees[i].is_none() {
+ let offset = self.buffer.len() - bit;
+ self.trees[i] = Some(self.buffer.drain(offset..).collect());
+ }
+ }
+
+ debug_assert!(self.buffer.len() < BUFFER_SIZE);
+ debug_assert!(self.len() == len);
+ }
+}
+
+impl<T, U> FromIterator<T> for Forest<U>
+where
+ U: FromIterator<T> + IntoIterator<Item = T>,
+{
+ fn from_iter<I: IntoIterator<Item = T>>(items: I) -> Self {
+ let mut forest = Self::new();
+ forest.extend(items);
+ forest
+ }
+}
+
+impl<T: IntoIterator> IntoIterator for Forest<T> {
+ type Item = T::Item;
+ type IntoIter = std::vec::IntoIter<T::Item>;
+
+ fn into_iter(mut self) -> Self::IntoIter {
+ self.buffer.extend(self.trees.into_iter().flatten().flatten());
+ self.buffer.into_iter()
+ }
+}
+
+impl<K, V, T> NearestNeighbors<K, V> for Forest<T>
+where
+ K: Proximity<V>,
+ T: NearestNeighbors<K, V>,
+ T: IntoIterator<Item = V>,
+{
+ fn search<'k, 'v, N>(&'v self, mut neighborhood: N) -> N
+ where
+ K: 'k,
+ V: 'v,
+ N: Neighborhood<&'k K, &'v V>
+ {
+ for item in &self.buffer {
+ neighborhood.consider(item);
+ }
+
+ self.trees
+ .iter()
+ .flatten()
+ .fold(neighborhood, |n, t| t.search(n))
+ }
+}
+
+/// A forest of k-d trees.
+pub type KdForest<T> = Forest<FlatKdTree<T>>;
+
+/// A forest of vantage-point trees.
+pub type VpForest<T> = Forest<FlatVpTree<T>>;
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+
+ use acap::euclid::Euclidean;
+ use acap::exhaustive::ExhaustiveSearch;
+ use acap::{NearestNeighbors, Neighbor};
+
+ use rand::prelude::*;
+
+ type Point = Euclidean<[f32; 3]>;
+
+ fn test_empty<T, F>(from_iter: &F)
+ where
+ T: NearestNeighbors<Point>,
+ F: Fn(Vec<Point>) -> T,
+ {
+ let points = Vec::new();
+ let index = from_iter(points);
+ let target = Euclidean([0.0, 0.0, 0.0]);
+ assert_eq!(index.nearest(&target), None);
+ assert_eq!(index.nearest_within(&target, 1.0), None);
+ assert!(index.k_nearest(&target, 0).is_empty());
+ assert!(index.k_nearest(&target, 3).is_empty());
+ assert!(index.k_nearest_within(&target, 0, 1.0).is_empty());
+ assert!(index.k_nearest_within(&target, 3, 1.0).is_empty());
+ }
+
+ fn test_pythagorean<T, F>(from_iter: &F)
+ where
+ T: NearestNeighbors<Point>,
+ F: Fn(Vec<Point>) -> T,
+ {
+ let points = vec![
+ Euclidean([3.0, 4.0, 0.0]),
+ Euclidean([5.0, 0.0, 12.0]),
+ Euclidean([0.0, 8.0, 15.0]),
+ Euclidean([1.0, 2.0, 2.0]),
+ Euclidean([2.0, 3.0, 6.0]),
+ Euclidean([4.0, 4.0, 7.0]),
+ ];
+ let index = from_iter(points);
+ let target = Euclidean([0.0, 0.0, 0.0]);
+
+ assert_eq!(
+ index.nearest(&target).expect("No nearest neighbor found"),
+ Neighbor::new(&Euclidean([1.0, 2.0, 2.0]), 3.0)
+ );
+
+ assert_eq!(index.nearest_within(&target, 2.0), None);
+ assert_eq!(
+ index.nearest_within(&target, 4.0).expect("No nearest neighbor found within 4.0"),
+ Neighbor::new(&Euclidean([1.0, 2.0, 2.0]), 3.0)
+ );
+
+ assert!(index.k_nearest(&target, 0).is_empty());
+ assert_eq!(
+ index.k_nearest(&target, 3),
+ vec![
+ Neighbor::new(&Euclidean([1.0, 2.0, 2.0]), 3.0),
+ Neighbor::new(&Euclidean([3.0, 4.0, 0.0]), 5.0),
+ Neighbor::new(&Euclidean([2.0, 3.0, 6.0]), 7.0),
+ ]
+ );
+
+ assert!(index.k_nearest(&target, 0).is_empty());
+ assert_eq!(
+ index.k_nearest_within(&target, 3, 6.0),
+ vec![
+ Neighbor::new(&Euclidean([1.0, 2.0, 2.0]), 3.0),
+ Neighbor::new(&Euclidean([3.0, 4.0, 0.0]), 5.0),
+ ]
+ );
+ assert_eq!(
+ index.k_nearest_within(&target, 3, 8.0),
+ vec![
+ Neighbor::new(&Euclidean([1.0, 2.0, 2.0]), 3.0),
+ Neighbor::new(&Euclidean([3.0, 4.0, 0.0]), 5.0),
+ Neighbor::new(&Euclidean([2.0, 3.0, 6.0]), 7.0),
+ ]
+ );
+ }
+
+ fn test_random_points<T, F>(from_iter: &F)
+ where
+ T: NearestNeighbors<Point>,
+ F: Fn(Vec<Point>) -> T,
+ {
+ let mut points = Vec::new();
+ for _ in 0..255 {
+ points.push(Euclidean([random(), random(), random()]));
+ }
+ let target = Euclidean([random(), random(), random()]);
+
+ let eindex = ExhaustiveSearch::from_iter(points.clone());
+ let index = from_iter(points);
+
+ assert_eq!(index.k_nearest(&target, 3), eindex.k_nearest(&target, 3));
+ }
+
+ /// Test a [NearestNeighbors] impl.
+ fn test_nearest_neighbors<T, F>(from_iter: F)
+ where
+ T: NearestNeighbors<Point>,
+ F: Fn(Vec<Point>) -> T,
+ {
+ test_empty(&from_iter);
+ test_pythagorean(&from_iter);
+ test_random_points(&from_iter);
+ }
+
+ #[test]
+ fn test_exhaustive_forest() {
+ test_nearest_neighbors(Forest::<ExhaustiveSearch<_>>::from_iter);
+ }
+
+ #[test]
+ fn test_forest_forest() {
+ test_nearest_neighbors(Forest::<Forest<ExhaustiveSearch<_>>>::from_iter);
+ }
+
+ #[test]
+ fn test_kd_forest() {
+ test_nearest_neighbors(KdForest::from_iter);
+ }
+
+ #[test]
+ fn test_vp_forest() {
+ test_nearest_neighbors(VpForest::from_iter);
+ }
+}
diff --git a/src/frontier.rs b/src/frontier.rs
new file mode 100644
index 0000000..74c7398
--- /dev/null
+++ b/src/frontier.rs
@@ -0,0 +1,194 @@
+//! Frontiers on which to place pixels.
+
+pub mod image;
+pub mod mean;
+pub mod min;
+
+use crate::color::{ColorSpace, Rgb8};
+use crate::soft::SoftDelete;
+
+use acap::coords::{Coordinates, CoordinateMetric, CoordinateProximity};
+use acap::distance::{Proximity, Metric};
+
+use std::cell::Cell;
+use std::ops::Deref;
+use std::rc::Rc;
+
+/// A frontier of pixels.
+pub trait Frontier {
+ /// The width of the image.
+ fn width(&self) -> u32;
+
+ /// The height of the image.
+ fn height(&self) -> u32;
+
+ /// The number of pixels currently on the frontier.
+ fn len(&self) -> usize;
+
+ /// Whether the frontier is empty.
+ fn is_empty(&self) -> bool;
+
+ /// Place the given color on the frontier, and return its position.
+ fn place(&mut self, rgb8: Rgb8) -> Option<(u32, u32)>;
+}
+
+/// A pixel on a frontier.
+#[derive(Debug)]
+struct Pixel<C> {
+ pos: (u32, u32),
+ color: C,
+ deleted: Cell<bool>,
+}
+
+impl<C: ColorSpace> Pixel<C> {
+ fn new(x: u32, y: u32, color: C) -> Self {
+ Self {
+ pos: (x, y),
+ color,
+ deleted: Cell::new(false),
+ }
+ }
+
+ fn delete(&self) {
+ self.deleted.set(true);
+ }
+}
+
+/// A reference-counted pixel, to work around the coherence rules.
+#[derive(Clone, Debug)]
+struct RcPixel<C>(Rc<Pixel<C>>);
+
+impl<C: ColorSpace> RcPixel<C> {
+ fn new(x: u32, y: u32, color: C) -> Self {
+ Self(Rc::new(Pixel::new(x, y, color)))
+ }
+}
+
+impl<C> Deref for RcPixel<C> {
+ type Target = Pixel<C>;
+
+ fn deref(&self) -> &Self::Target {
+ self.0.deref()
+ }
+}
+
+/// A search target, to work around the coherence rules.
+#[derive(Debug)]
+struct Target<C>(C);
+
+impl<C: Proximity> Proximity<Pixel<C>> for Target<C> {
+ type Distance = C::Distance;
+
+ fn distance(&self, other: &Pixel<C>) -> Self::Distance {
+ self.0.distance(&other.color)
+ }
+}
+
+impl<C: Metric> Metric<Pixel<C>> for Target<C> {}
+
+impl<C: Proximity> Proximity for Pixel<C> {
+ type Distance = C::Distance;
+
+ fn distance(&self, other: &Pixel<C>) -> Self::Distance {
+ self.color.distance(&other.color)
+ }
+}
+
+impl<C: Metric> Metric for Pixel<C> {}
+
+impl<C: Coordinates> Coordinates for Pixel<C> {
+ type Value = C::Value;
+
+ fn dims(&self) -> usize {
+ self.color.dims()
+ }
+
+ fn coord(&self, i: usize) -> Self::Value {
+ self.color.coord(i)
+ }
+}
+
+impl<C> SoftDelete for Pixel<C> {
+ fn is_deleted(&self) -> bool {
+ self.deleted.get()
+ }
+}
+
+impl<C: Proximity> Proximity<RcPixel<C>> for Target<C> {
+ type Distance = C::Distance;
+
+ fn distance(&self, other: &RcPixel<C>) -> Self::Distance {
+ self.0.distance(&other.0.color)
+ }
+}
+
+impl<C: Metric> Metric<RcPixel<C>> for Target<C> {}
+
+impl<C: Coordinates> Coordinates for Target<C> {
+ type Value = C::Value;
+
+ fn dims(&self) -> usize {
+ self.0.dims()
+ }
+
+ fn coord(&self, i: usize) -> Self::Value {
+ self.0.coord(i)
+ }
+}
+
+impl<T, C: CoordinateProximity<T>> CoordinateProximity<T> for Target<C> {
+ type Distance = C::Distance;
+
+ fn distance_to_coords(&self, coords: &[T]) -> Self::Distance {
+ self.0.distance_to_coords(coords)
+ }
+}
+
+impl<T, C: CoordinateMetric<T>> CoordinateMetric<T> for Target<C> {}
+
+impl<C: Proximity> Proximity for RcPixel<C> {
+ type Distance = C::Distance;
+
+ fn distance(&self, other: &Self) -> Self::Distance {
+ (*self.0).distance(&*other.0)
+ }
+}
+
+impl<C: Metric> Metric for RcPixel<C> {}
+
+impl<C: Coordinates> Coordinates for RcPixel<C> {
+ type Value = C::Value;
+
+ fn dims(&self) -> usize {
+ (*self.0).dims()
+ }
+
+ fn coord(&self, i: usize) -> Self::Value {
+ (*self.0).coord(i)
+ }
+}
+
+impl<C> SoftDelete for RcPixel<C> {
+ fn is_deleted(&self) -> bool {
+ (*self.0).is_deleted()
+ }
+}
+
+/// Return all the neighbors of a pixel location.
+fn neighbors(x: u32, y: u32) -> [(u32, u32); 8] {
+ let xm1 = x.wrapping_sub(1);
+ let ym1 = y.wrapping_sub(1);
+ let xp1 = x + 1;
+ let yp1 = y + 1;
+
+ [
+ (xm1, ym1),
+ (xm1, y),
+ (xm1, yp1),
+ (x, ym1),
+ (x, yp1),
+ (xp1, ym1),
+ (xp1, y),
+ (xp1, yp1),
+ ]
+}
diff --git a/src/frontier/image.rs b/src/frontier/image.rs
new file mode 100644
index 0000000..18bf620
--- /dev/null
+++ b/src/frontier/image.rs
@@ -0,0 +1,79 @@
+//! Frontier that targets an image.
+
+use super::{Frontier, Pixel, Target};
+
+use crate::color::{ColorSpace, Rgb8};
+use crate::soft::SoftKdTree;
+
+use acap::NearestNeighbors;
+
+use image::RgbImage;
+
+/// A [Frontier] that places colors on the closest pixel of a target image.
+#[derive(Debug)]
+pub struct ImageFrontier<C> {
+ nodes: SoftKdTree<Pixel<C>>,
+ width: u32,
+ height: u32,
+ len: usize,
+ deleted: usize,
+}
+
+impl<C: ColorSpace> ImageFrontier<C> {
+ /// Create an ImageFrontier from an image.
+ pub fn new(img: &RgbImage) -> Self {
+ let width = img.width();
+ let height = img.height();
+ let len = (width as usize) * (height as usize);
+
+ Self {
+ nodes: img
+ .enumerate_pixels()
+ .map(|(x, y, p)| Pixel::new(x, y, C::from(*p)))
+ .collect(),
+ width,
+ height,
+ len,
+ deleted: 0,
+ }
+ }
+}
+
+impl<C: ColorSpace> Frontier for ImageFrontier<C> {
+ fn width(&self) -> u32 {
+ self.width
+ }
+
+ fn height(&self) -> u32 {
+ self.height
+ }
+
+ fn len(&self) -> usize {
+ self.len - self.deleted
+ }
+
+ fn is_empty(&self) -> bool {
+ self.len() == 0
+ }
+
+ fn place(&mut self, rgb8: Rgb8) -> Option<(u32, u32)> {
+ let color = C::from(rgb8);
+
+ if let Some(node) = self.nodes.nearest(&Target(color)).map(|n| n.item) {
+ let pos = node.pos;
+
+ node.delete();
+ self.deleted += 1;
+
+ if 32 * self.deleted >= self.len {
+ self.nodes.rebuild();
+ self.len -= self.deleted;
+ self.deleted = 0;
+ }
+
+ Some(pos)
+ } else {
+ None
+ }
+ }
+}
diff --git a/src/frontier/mean.rs b/src/frontier/mean.rs
new file mode 100644
index 0000000..3c441b8
--- /dev/null
+++ b/src/frontier/mean.rs
@@ -0,0 +1,145 @@
+//! Mean selection frontier.
+
+use super::{neighbors, Frontier, RcPixel, Target};
+
+use crate::color::{ColorSpace, Rgb8};
+use crate::soft::SoftKdForest;
+
+use acap::NearestNeighbors;
+
+use std::iter;
+
+/// A pixel on a mean frontier.
+#[derive(Debug)]
+enum MeanPixel<C> {
+ Empty,
+ Fillable(RcPixel<C>),
+ Filled(C),
+}
+
+impl<C: ColorSpace> MeanPixel<C> {
+ fn filled_color(&self) -> Option<C> {
+ match self {
+ Self::Filled(color) => Some(*color),
+ _ => None,
+ }
+ }
+}
+
+/// A [Frontier] that looks at the average color of each pixel's neighbors.
+#[derive(Debug)]
+pub struct MeanFrontier<C> {
+ pixels: Vec<MeanPixel<C>>,
+ forest: SoftKdForest<RcPixel<C>>,
+ width: u32,
+ height: u32,
+ len: usize,
+ deleted: usize,
+}
+
+impl<C: ColorSpace> MeanFrontier<C> {
+ /// Create a MeanFrontier with the given dimensions and initial pixel location.
+ pub fn new(width: u32, height: u32, x0: u32, y0: u32) -> Self {
+ let size = (width as usize) * (height as usize);
+ let mut pixels = Vec::with_capacity(size);
+ for _ in 0..size {
+ pixels.push(MeanPixel::Empty);
+ }
+
+ let pixel0 = RcPixel::new(x0, y0, C::from(Rgb8::from([0, 0, 0])));
+ let i = (x0 + y0 * width) as usize;
+ pixels[i] = MeanPixel::Fillable(pixel0.clone());
+
+ Self {
+ pixels,
+ forest: iter::once(pixel0).collect(),
+ width,
+ height,
+ len: 1,
+ deleted: 0,
+ }
+ }
+
+ fn pixel_index(&self, x: u32, y: u32) -> usize {
+ debug_assert!(x < self.width);
+ debug_assert!(y < self.height);
+
+ (x + y * self.width) as usize
+ }
+
+ fn fill(&mut self, x: u32, y: u32, color: C) {
+ let i = self.pixel_index(x, y);
+ match &self.pixels[i] {
+ MeanPixel::Empty => {}
+ MeanPixel::Fillable(pixel) => {
+ pixel.delete();
+ self.deleted += 1;
+ }
+ _ => unreachable!(),
+ }
+ self.pixels[i] = MeanPixel::Filled(color);
+
+ let mut pixels = Vec::new();
+ for &(x, y) in &neighbors(x, y) {
+ if x < self.width && y < self.height {
+ let i = self.pixel_index(x, y);
+ match &self.pixels[i] {
+ MeanPixel::Empty => {}
+ MeanPixel::Fillable(pixel) => {
+ pixel.delete();
+ self.deleted += 1;
+ }
+ MeanPixel::Filled(_) => continue,
+ }
+ let color = C::average(
+ neighbors(x, y)
+ .iter()
+ .filter(|(x, y)| *x < self.width && *y < self.height)
+ .map(|(x, y)| self.pixel_index(*x, *y))
+ .map(|i| &self.pixels[i])
+ .map(MeanPixel::filled_color)
+ .flatten(),
+ );
+ let pixel = RcPixel::new(x, y, color);
+ self.pixels[i] = MeanPixel::Fillable(pixel.clone());
+ pixels.push(pixel);
+ }
+ }
+
+ self.len += pixels.len();
+ self.forest.extend(pixels);
+
+ if 2 * self.deleted >= self.len {
+ self.forest.rebuild();
+ self.len -= self.deleted;
+ self.deleted = 0;
+ }
+ }
+}
+
+impl<C: ColorSpace> Frontier for MeanFrontier<C> {
+ fn width(&self) -> u32 {
+ self.width
+ }
+
+ fn height(&self) -> u32 {
+ self.height
+ }
+
+ fn len(&self) -> usize {
+ self.len - self.deleted
+ }
+
+ fn is_empty(&self) -> bool {
+ self.len() == 0
+ }
+
+ fn place(&mut self, rgb8: Rgb8) -> Option<(u32, u32)> {
+ let color = C::from(rgb8);
+ let (x, y) = self.forest.nearest(&Target(color)).map(|n| n.item.pos)?;
+
+ self.fill(x, y, color);
+
+ Some((x, y))
+ }
+}
diff --git a/src/frontier/min.rs b/src/frontier/min.rs
new file mode 100644
index 0000000..95b3321
--- /dev/null
+++ b/src/frontier/min.rs
@@ -0,0 +1,153 @@
+//! Minimum selection frontier.
+
+use super::{neighbors, Frontier, RcPixel, Target};
+
+use crate::color::{ColorSpace, Rgb8};
+use crate::soft::SoftKdForest;
+
+use acap::NearestNeighbors;
+
+use rand::Rng;
+
+/// A pixel on a min frontier.
+#[derive(Debug)]
+struct MinPixel<C> {
+ pixel: Option<RcPixel<C>>,
+ filled: bool,
+}
+
+impl<C: ColorSpace> MinPixel<C> {
+ fn new() -> Self {
+ Self {
+ pixel: None,
+ filled: false,
+ }
+ }
+}
+
+/// A [Frontier] that places colors on a neighbor of the closest pixel so far.
+#[derive(Debug)]
+pub struct MinFrontier<C, R> {
+ rng: R,
+ pixels: Vec<MinPixel<C>>,
+ forest: SoftKdForest<RcPixel<C>>,
+ width: u32,
+ height: u32,
+ x0: u32,
+ y0: u32,
+ len: usize,
+ deleted: usize,
+}
+
+impl<C: ColorSpace, R: Rng> MinFrontier<C, R> {
+ /// Create a MinFrontier with the given dimensions and initial pixel location.
+ pub fn new(rng: R, width: u32, height: u32, x0: u32, y0: u32) -> Self {
+ let size = (width as usize) * (height as usize);
+ let mut pixels = Vec::with_capacity(size);
+ for _ in 0..size {
+ pixels.push(MinPixel::new());
+ }
+
+ Self {
+ rng,
+ pixels,
+ forest: SoftKdForest::new(),
+ width,
+ height,
+ x0,
+ y0,
+ len: 0,
+ deleted: 0,
+ }
+ }
+
+ fn pixel_index(&self, x: u32, y: u32) -> usize {
+ debug_assert!(x < self.width);
+ debug_assert!(y < self.height);
+
+ (x + y * self.width) as usize
+ }
+
+ fn free_neighbor(&mut self, x: u32, y: u32) -> Option<(u32, u32)> {
+ // Pick a pseudo-random neighbor
+ let offset: usize = self.rng.gen();
+
+ let neighbors = neighbors(x, y);
+ for i in 0..8 {
+ let (x, y) = neighbors[(i + offset) % 8];
+ if x < self.width && y < self.height {
+ let i = self.pixel_index(x, y);
+ if !self.pixels[i].filled {
+ return Some((x, y));
+ }
+ }
+ }
+
+ None
+ }
+
+ fn fill(&mut self, x: u32, y: u32, color: C) -> Option<(u32, u32)> {
+ let i = self.pixel_index(x, y);
+ let pixel = &mut self.pixels[i];
+ if pixel.filled {
+ return None;
+ }
+
+ let rc = RcPixel::new(x, y, color);
+ pixel.pixel = Some(rc.clone());
+ pixel.filled = true;
+
+ if self.free_neighbor(x, y).is_some() {
+ self.forest.push(rc);
+ self.len += 1;
+ }
+
+ for &(x, y) in &neighbors(x, y) {
+ if x < self.width && y < self.height && self.free_neighbor(x, y).is_none() {
+ let i = self.pixel_index(x, y);
+ if let Some(pixel) = self.pixels[i].pixel.take() {
+ pixel.delete();
+ self.deleted += 1;
+ }
+ }
+ }
+
+ if 2 * self.deleted >= self.len {
+ self.forest.rebuild();
+ self.len -= self.deleted;
+ self.deleted = 0;
+ }
+
+ Some((x, y))
+ }
+}
+
+impl<C: ColorSpace, R: Rng> Frontier for MinFrontier<C, R> {
+ fn width(&self) -> u32 {
+ self.width
+ }
+
+ fn height(&self) -> u32 {
+ self.height
+ }
+
+ fn len(&self) -> usize {
+ self.len - self.deleted
+ }
+
+ fn is_empty(&self) -> bool {
+ self.len() == 0
+ }
+
+ fn place(&mut self, rgb8: Rgb8) -> Option<(u32, u32)> {
+ let color = C::from(rgb8);
+ let (x, y) = self
+ .forest
+ .nearest(&Target(color))
+ .map(|n| n.item.pos)
+ .map(|(x, y)| self.free_neighbor(x, y).unwrap())
+ .unwrap_or((self.x0, self.y0));
+
+ self.fill(x, y, color)
+ }
+}
diff --git a/src/hilbert.rs b/src/hilbert.rs
new file mode 100644
index 0000000..36464dc
--- /dev/null
+++ b/src/hilbert.rs
@@ -0,0 +1,136 @@
+//! Implementation of [Compact Hilbert Indices](https://dl.acm.org/doi/10.1109/CISIS.2007.16) by
+//! Chris Hamilton.
+
+/// Right rotation of x by b bits out of n.
+fn rotate_right(x: usize, b: u32, n: u32) -> usize {
+ let l = x & ((1 << b) - 1);
+ let r = x >> b;
+ (l << (n - b)) | r
+}
+
+/// Left rotation of x by b bits out of n.
+fn rotate_left(x: usize, b: u32, n: u32) -> usize {
+ rotate_right(x, n - b, n)
+}
+
+/// Binary reflected Gray code.
+fn gray_code(i: usize) -> usize {
+ i ^ (i >> 1)
+}
+
+/// e(i), the entry point for the ith sub-hypercube.
+fn entry_point(i: usize) -> usize {
+ if i == 0 {
+ 0
+ } else {
+ gray_code((i - 1) & !1)
+ }
+}
+
+/// g(i), the inter sub-hypercube direction.
+fn inter_direction(i: usize) -> u32 {
+ // g(i) counts the trailing set bits in i
+ (!i).trailing_zeros()
+}
+
+/// d(i), the intra sub-hypercube direction.
+fn intra_direction(i: usize) -> u32 {
+ if i & 1 != 0 {
+ inter_direction(i)
+ } else if i > 0 {
+ inter_direction(i - 1)
+ } else {
+ 0
+ }
+}
+
+/// T transformation inverse
+fn t_inverse(dims: u32, e: usize, d: u32, a: usize) -> usize {
+ rotate_left(a, d, dims) ^ e
+}
+
+/// GrayCodeRankInverse
+fn gray_code_rank_inverse(
+ dims: u32,
+ mu: usize,
+ pi: usize,
+ r: usize,
+ free_bits: u32,
+) -> (usize, usize) {
+ // The inverse rank of r
+ let mut i = 0;
+ // gray_code(i)
+ let mut g = 0;
+
+ let mut j = free_bits - 1;
+ for k in (0..dims).rev() {
+ if mu & (1 << k) == 0 {
+ g |= pi & (1 << k);
+ i |= (g ^ (i >> 1)) & (1 << k);
+ } else {
+ i |= ((r >> j) & 1) << k;
+ g |= (i ^ (i >> 1)) & (1 << k);
+ j = j.wrapping_sub(1);
+ }
+ }
+
+ (i, g)
+}
+
+/// ExtractMask.
+fn extract_mask(bits: &[u32], i: u32) -> (usize, u32) {
+ // The mask
+ let mut mu = 0;
+ // popcount(mu)
+ let mut free_bits = 0;
+
+ let dims = bits.len();
+ for j in (0..dims).rev() {
+ mu <<= 1;
+ if bits[j] > i {
+ mu |= 1;
+ free_bits += 1;
+ }
+ }
+
+ (mu, free_bits)
+}
+
+/// Compute the corresponding point for a Hilbert index (CompactHilbertIndexInverse).
+pub fn hilbert_point(index: usize, bits: &[u32], point: &mut [usize]) {
+ let dims = bits.len() as u32;
+ let max = *bits.iter().max().unwrap();
+ let sum: u32 = bits.iter().sum();
+
+ let mut e = 0;
+ let mut k = 0;
+
+ // Next direction; we use d instead of d + 1 everywhere
+ let mut d = 1;
+
+ for x in point.iter_mut() {
+ *x = 0;
+ }
+
+ for i in (0..max).rev() {
+ let (mut mu, free_bits) = extract_mask(bits, i);
+ mu = rotate_right(mu, d, dims);
+
+ let pi = rotate_right(e, d, dims) & !mu;
+
+ let r = (index >> (sum - k - free_bits)) & ((1 << free_bits) - 1);
+
+ k += free_bits;
+
+ let (w, mut l) = gray_code_rank_inverse(dims, mu, pi, r, free_bits);
+ l = t_inverse(dims, e, d, l);
+
+ for x in point.iter_mut() {
+ *x |= (l & 1) << i;
+ l >>= 1;
+ }
+
+ e ^= rotate_right(entry_point(w), d, dims);
+ d = (d + intra_direction(w) + 1) % dims;
+ }
+}
diff --git a/src/main.rs b/src/main.rs
new file mode 100644
index 0000000..ce54939
--- /dev/null
+++ b/src/main.rs
@@ -0,0 +1,492 @@
+pub mod color;
+pub mod forest;
+pub mod frontier;
+pub mod hilbert;
+pub mod soft;
+
+use crate::color::source::{AllColors, ColorSource, ImageColors};
+use crate::color::{order, ColorSpace, LabSpace, LuvSpace, Rgb8, RgbSpace};
+use crate::frontier::image::ImageFrontier;
+use crate::frontier::mean::MeanFrontier;
+use crate::frontier::min::MinFrontier;
+use crate::frontier::Frontier;
+
+use clap::{self, clap_app, crate_authors, crate_name, crate_version};
+
+use image::{self, ImageError, Rgba, RgbaImage};
+
+use rand::{self, SeedableRng};
+use rand_pcg::Pcg64;
+
+use std::cmp;
+use std::error::Error;
+use std::fs;
+use std::io::{self, Write};
+use std::path::PathBuf;
+use std::process::exit;
+use std::str::FromStr;
+use std::time::Instant;
+
+/// The color source specified on the command line.
+#[derive(Debug, Eq, PartialEq)]
+enum SourceArg {
+ /// All RGB colors of the given bit depth(s).
+ AllRgb(u32, u32, u32),
+ /// Take the colors from an image.
+ Image(PathBuf),
+}
+
+/// The order to process colors in.
+#[derive(Debug, Eq, PartialEq)]
+enum OrderArg {
+ /// Sorted by hue.
+ HueSort,
+ /// Shuffled randomly.
+ Random,
+ /// Morton/Z-order.
+ Morton,
+ /// Hilbert curve order.
+ Hilbert,
+}
+
+/// The frontier implementation.
+#[derive(Debug, Eq, PartialEq)]
+enum FrontierArg {
+ /// Pick a neighbor of the closest pixel so far.
+ Min,
+ /// Pick the pixel with the closest mean color of all its neighbors.
+ Mean,
+ /// Target the closest pixel on an image.
+ Image(PathBuf),
+}
+
+/// The color space to operate in.
+#[derive(Debug, Eq, PartialEq)]
+enum ColorSpaceArg {
+ /// sRGB space.
+ Rgb,
+ /// CIE L*a*b* space.
+ Lab,
+ /// CIE L*u*v* space.
+ Luv,
+}
+
+/// Error type for this app.
+#[derive(Debug)]
+enum AppError {
+ ArgError(clap::Error),
+ RuntimeError(Box<dyn Error>),
+}
+
+impl AppError {
+ /// Create an error for an invalid argument.
+ fn invalid_value(msg: &str) -> Self {
+ Self::ArgError(clap::Error::with_description(
+ msg,
+ clap::ErrorKind::InvalidValue,
+ ))
+ }
+
+ /// Exit the program with this error.
+ fn exit(&self) -> ! {
+ match self {
+ Self::ArgError(err) => err.exit(),
+ Self::RuntimeError(err) => {
+ eprintln!("{}", err);
+ exit(1)
+ }
+ }
+ }
+}
+
+impl From<clap::Error> for AppError {
+ fn from(err: clap::Error) -> Self {
+ Self::ArgError(err)
+ }
+}
+
+impl From<ImageError> for AppError {
+ fn from(err: ImageError) -> Self {
+ Self::RuntimeError(Box::new(err))
+ }
+}
+
+impl From<io::Error> for AppError {
+ fn from(err: io::Error) -> Self {
+ Self::RuntimeError(Box::new(err))
+ }
+}
+
+impl From<rand::Error> for AppError {
+ fn from(err: rand::Error) -> Self {
+ Self::RuntimeError(Box::new(err))
+ }
+}
+
+/// Result type for this app.
+type AppResult<T> = Result<T, AppError>;
+
+/// Parse an argument into the appropriate type.
+fn parse_arg<F>(arg: Option<&str>) -> AppResult<Option<F>>
+where
+ F: FromStr,
+ F::Err: Error,
+{
+ match arg.map(|s| s.parse()) {
+ Some(Ok(f)) => Ok(Some(f)),
+ Some(Err(e)) => Err(AppError::invalid_value(&e.to_string())),
+ None => Ok(None),
+ }
+}
+
+/// The parsed command line arguments.
+#[derive(Debug)]
+struct Args {
+ source: SourceArg,
+ order: OrderArg,
+ stripe: bool,
+ frontier: FrontierArg,
+ space: ColorSpaceArg,
+ width: Option<u32>,
+ height: Option<u32>,
+ x0: Option<u32>,
+ y0: Option<u32>,
+ animate: bool,
+ output: PathBuf,
+ seed: u64,
+}
+
+impl Args {
+ fn parse() -> AppResult<Self> {
+ let args = clap_app!((crate_name!()) =>
+ (version: crate_version!())
+ (author: crate_authors!())
+ (@setting ColoredHelp)
+ (@setting DeriveDisplayOrder)
+ (@setting UnifiedHelpMessage)
+ (@group source =>
+ (@arg DEPTH: -b --("bit-depth") +takes_value "Use all DEPTH-bit colors")
+ (@arg INPUT: -i --input +takes_value "Use colors from the INPUT image")
+ )
+ (@group order =>
+ (@arg HUE: -s --hue-sort "Sort colors by hue [default]")
+ (@arg RANDOM: -r --random "Randomize colors")
+ (@arg MORTON: -M --morton "Place colors in Morton order (Z-order)")
+ (@arg HILBERT: -H --hilbert "Place colors in Hilbert curve order")
+ )
+ (@group stripe =>
+ (@arg STRIPE: -t --stripe "Reduce artifacts by iterating through the colors in multiple stripes [default]")
+ (@arg NOSTRIPE: -T --("no-stripe") "Don't stripe")
+ )
+ (@group frontier =>
+ (@arg MODE: -l --selection +takes_value possible_value[min mean] "Specify the selection mode")
+ (@arg TARGET: -g --target +takes_value "Place colors on the closest pixels of the TARGET image")
+ )
+ (@arg SPACE: -c --("color-space") default_value("Lab") possible_value[RGB Lab Luv] "Use the given color space")
+ (@arg WIDTH: -w --width +takes_value "The width of the generated image")
+ (@arg HEIGHT: -h --height +takes_value "The height of the generated image")
+ (@arg X: -x +takes_value "The x coordinate of the first pixel")
+ (@arg Y: -y +takes_value "The y coordinate of the first pixel")
+ (@arg ANIMATE: -a --animate "Generate frames of an animation")
+ (@arg PATH: -o --output default_value("kd-forest.png") "Save the image to PATH")
+ (@arg SEED: -e --seed default_value("0") "Seed the random number generator")
+ )
+ .get_matches_safe()?;
+
+ let source = if let Some(input) = args.value_of("INPUT") {
+ SourceArg::Image(PathBuf::from(input))
+ } else {
+ let arg = args.value_of("DEPTH");
+ let depths: Vec<_> = arg
+ .iter()
+ .map(|s| s.split(','))
+ .flatten()
+ .map(|n| n.parse().ok())
+ .collect();
+
+ let (r, g, b) = match depths.as_slice() {
+ [] => (8, 8, 8),
+
+ // Allocate bits from most to least perceptually important
+ [Some(d)] => ((d + 1) / 3, (d + 2) / 3, d / 3),
+
+ [Some(r), Some(g), Some(b)] => (*r, *g, *b),
+
+ _ => {
+ return Err(AppError::invalid_value(
+ &format!("invalid bit depth {}", arg.unwrap()),
+ ));
+ }
+ };
+
+ if r > 8 || g > 8 || b > 8 {
+ return Err(AppError::invalid_value(
+ &format!("bit depth of {} is too deep!", arg.unwrap()),
+ ));
+ }
+
+ SourceArg::AllRgb(r, g, b)
+ };
+
+ let order = if args.is_present("RANDOM") {
+ OrderArg::Random
+ } else if args.is_present("MORTON") {
+ OrderArg::Morton
+ } else if args.is_present("HILBERT") {
+ OrderArg::Hilbert
+ } else {
+ OrderArg::HueSort
+ };
+
+ let stripe = !args.is_present("NOSTRIPE") && order != OrderArg::Random;
+
+ let frontier = if let Some(target) = args.value_of("TARGET") {
+ FrontierArg::Image(PathBuf::from(target))
+ } else {
+ match args.value_of("MODE") {
+ Some("min") | None => FrontierArg::Min,
+ Some("mean") => FrontierArg::Mean,
+ _ => unreachable!(),
+ }
+ };
+
+ let space = match args.value_of("SPACE").unwrap() {
+ "RGB" => ColorSpaceArg::Rgb,
+ "Lab" => ColorSpaceArg::Lab,
+ "Luv" => ColorSpaceArg::Luv,
+ _ => unreachable!(),
+ };
+
+ let width = parse_arg(args.value_of("WIDTH"))?;
+ let height = parse_arg(args.value_of("HEIGHT"))?;
+ let x0 = parse_arg(args.value_of("X"))?;
+ let y0 = parse_arg(args.value_of("Y"))?;
+
+ let animate = args.is_present("ANIMATE");
+
+ let path = if animate && args.occurrences_of("PATH") == 0 {
+ "kd-frames"
+ } else {
+ args.value_of("PATH").unwrap()
+ };
+ let output = PathBuf::from(path);
+
+ let seed = parse_arg(args.value_of("SEED"))?.unwrap_or(0);
+
+ Ok(Self {
+ source,
+ order,
+ stripe,
+ frontier,
+ space,
+ width,
+ height,
+ x0,
+ y0,
+ animate,
+ output,
+ seed,
+ })
+ }
+}
+
+/// The kd-forest application itself.
+#[derive(Debug)]
+struct App {
+ args: Args,
+ rng: Pcg64,
+ width: Option<u32>,
+ height: Option<u32>,
+ start_time: Instant,
+}
+
+impl App {
+ /// Make the App.
+ fn new(args: Args) -> Self {
+ let rng = Pcg64::seed_from_u64(args.seed);
+ let width = args.width;
+ let height = args.height;
+ let start_time = Instant::now();
+
+ Self {
+ args,
+ rng,
+ width,
+ height,
+ start_time,
+ }
+ }
+
+ fn run(&mut self) -> AppResult<()> {
+ let colors = match self.args.source {
+ SourceArg::AllRgb(r, g, b) => {
+ let total = r + g + b;
+ self.width.get_or_insert(1u32 << ((total + 1) / 2));
+ self.height.get_or_insert(1u32 << (total / 2));
+ self.get_colors(AllColors::new(r, g, b))
+ }
+ SourceArg::Image(ref path) => {
+ let img = image::open(path)?.into_rgb();
+ self.width.get_or_insert(img.width());
+ self.height.get_or_insert(img.height());
+ self.get_colors(ImageColors::from(img))
+ }
+ };
+
+ match self.args.space {
+ ColorSpaceArg::Rgb => self.paint::<RgbSpace>(colors),
+ ColorSpaceArg::Lab => self.paint::<LabSpace>(colors),
+ ColorSpaceArg::Luv => self.paint::<LuvSpace>(colors),
+ }
+ }
+
+ fn get_colors<S: ColorSource>(&mut self, source: S) -> Vec<Rgb8> {
+ let colors = match self.args.order {
+ OrderArg::HueSort => order::hue_sorted(source),
+ OrderArg::Random => order::shuffled(source, &mut self.rng),
+ OrderArg::Morton => order::morton(source),
+ OrderArg::Hilbert => order::hilbert(source),
+ };
+
+ if self.args.stripe {
+ order::striped(colors)
+ } else {
+ colors
+ }
+ }
+
+ fn paint<C: ColorSpace>(&mut self, colors: Vec<Rgb8>) -> AppResult<()> {
+ let width = self.width.unwrap();
+ let height = self.height.unwrap();
+ let x0 = self.args.x0.unwrap_or(width / 2);
+ let y0 = self.args.y0.unwrap_or(height / 2);
+
+ if x0 >= width || y0 >= height {
+ return Err(AppError::invalid_value(
+ &format!("Initial pixel ({}, {}) is out of bounds ({}, {})", x0, y0, width, height),
+ ));
+ }
+
+ match &self.args.frontier {
+ FrontierArg::Image(ref path) => {
+ let img = image::open(path)?.into_rgb();
+ self.paint_on(colors, ImageFrontier::<C>::new(&img))
+ }
+ FrontierArg::Min => {
+ let rng = Pcg64::from_rng(&mut self.rng)?;
+ self.paint_on(colors, MinFrontier::<C, _>::new(rng, width, height, x0, y0))
+ }
+ FrontierArg::Mean => {
+ self.paint_on(colors, MeanFrontier::<C>::new(width, height, x0, y0))
+ }
+ }
+ }
+
+ fn paint_on<F: Frontier>(&mut self, colors: Vec<Rgb8>, mut frontier: F) -> AppResult<()> {
+ let width = frontier.width();
+ let height = frontier.height();
+ let mut output = RgbaImage::new(width, height);
+
+ let size = cmp::min((width * height) as usize, colors.len());
+ println!("Generating a {}x{} image ({} pixels)", width, height, size);
+
+ if self.args.animate {
+ fs::create_dir_all(&self.args.output)?;
+ output.save(&self.args.output.join("0000.png"))?;
+ }
+
+ let interval = cmp::max(width, height) as usize;
+
+ let mut max_frontier = frontier.len();
+
+ for (i, color) in colors.into_iter().enumerate() {
+ let pos = frontier.place(color);
+ if pos.is_none() {
+ break;
+ }
+
+ let (x, y) = pos.unwrap();
+ let rgba = Rgba([color[0], color[1], color[2], 255]);
+ output.put_pixel(x, y, rgba);
+
+ max_frontier = cmp::max(max_frontier, frontier.len());
+
+ if (i + 1) % interval == 0 {
+ if self.args.animate {
+ let frame = (i + 1) / interval;
+ output.save(&self.args.output.join(format!("{:04}.png", frame)))?;
+ }
+
+ if i + 1 < size {
+ self.print_progress(i + 1, size, frontier.len())?;
+ }
+ }
+ }
+
+ if self.args.animate && size % interval != 0 {
+ let frame = size / interval;
+ output.save(&self.args.output.join(format!("{:04}.png", frame)))?;
+ }
+
+ self.print_progress(size, size, max_frontier)?;
+
+ if !self.args.animate {
+ output.save(&self.args.output)?;
+ }
+
+ Ok(())
+ }
+
+ fn print_progress(&self, i: usize, size: usize, frontier_len: usize) -> io::Result<()> {
+ let mut term = match term::stderr() {
+ Some(term) => term,
+ None => return Ok(()),
+ };
+
+ let progress = 100.0 * (i as f64) / (size as f64);
+ let mut rate = (i as f64) / self.start_time.elapsed().as_secs_f64();
+ let mut unit = "px/s";
+
+ if rate >= 10_000.0 {
+ rate /= 1_000.0;
+ unit = "Kpx/s";
+ }
+
+ if rate >= 10_000.0 {
+ rate /= 1_000.0;
+ unit = "Mpx/s";
+ }
+
+ if rate >= 10_000.0 {
+ rate /= 1_000.0;
+ unit = "Gpx/s";
+ }
+
+ let (frontier_label, newline) = if i == size {
+ ("max frontier size", "\n")
+ } else {
+ ("frontier size", "")
+ };
+
+ term.carriage_return()?;
+ term.delete_line()?;
+
+ write!(
+ term,
+ "{:>6.2}% | {:4.0} {:>5} | {}: {}{}",
+ progress, rate, unit, frontier_label, frontier_len, newline,
+ )
+ }
+}
+
+fn main() {
+ let args = match Args::parse() {
+ Ok(args) => args,
+ Err(e) => e.exit(),
+ };
+
+ match App::new(args).run() {
+ Ok(_) => {},
+ Err(e) => e.exit(),
+ }
+}
diff --git a/src/soft.rs b/src/soft.rs
new file mode 100644
index 0000000..3edaa10
--- /dev/null
+++ b/src/soft.rs
@@ -0,0 +1,289 @@
+//! [Soft deletion](https://en.wiktionary.org/wiki/soft_deletion) for nearest neighbor search.
+
+use super::forest::{KdForest, VpForest};
+
+use acap::distance::Proximity;
+use acap::kd::FlatKdTree;
+use acap::vp::FlatVpTree;
+use acap::{NearestNeighbors, Neighborhood};
+
+use std::iter;
+use std::iter::FromIterator;
+use std::mem;
+
+/// A trait for objects that can be soft-deleted.
+pub trait SoftDelete {
+ /// Check whether this item is deleted.
+ fn is_deleted(&self) -> bool;
+}
+
+/// Blanket [SoftDelete] implementation for references.
+impl<'a, T: SoftDelete> SoftDelete for &'a T {
+ fn is_deleted(&self) -> bool {
+ (*self).is_deleted()
+ }
+}
+
+/// [Neighborhood] wrapper that ignores soft-deleted items.
+#[derive(Debug)]
+struct SoftNeighborhood<N>(N);
+
+impl<K, V, N> Neighborhood<K, V> for SoftNeighborhood<N>
+where
+ V: SoftDelete,
+ K: Proximity<V>,
+ N: Neighborhood<K, V>,
+{
+ fn target(&self) -> K {
+ self.0.target()
+ }
+
+ fn contains<D>(&self, distance: D) -> bool
+ where
+ D: PartialOrd<K::Distance>
+ {
+ self.0.contains(distance)
+ }
+
+ fn consider(&mut self, item: V) -> K::Distance {
+ if item.is_deleted() {
+ self.target().distance(&item)
+ } else {
+ self.0.consider(item)
+ }
+ }
+}
+
+/// A [NearestNeighbors] implementation that supports [soft deletes](https://en.wiktionary.org/wiki/soft_deletion).
+#[derive(Debug)]
+pub struct SoftSearch<T>(T);
+
+impl<T, U> SoftSearch<U>
+where
+ T: SoftDelete,
+ U: FromIterator<T> + IntoIterator<Item = T>,
+{
+ /// Create a new empty soft index.
+ pub fn new() -> Self {
+ Self(iter::empty().collect())
+ }
+
+ /// Push a new item into this index.
+ pub fn push(&mut self, item: T)
+ where
+ U: Extend<T>,
+ {
+ self.0.extend(iter::once(item));
+ }
+
+ /// Rebuild this index, discarding deleted items.
+ pub fn rebuild(&mut self) {
+ let items = mem::replace(&mut self.0, iter::empty().collect());
+ self.0 = items.into_iter().filter(|e| !e.is_deleted()).collect();
+ }
+}
+
+impl<T, U> Default for SoftSearch<U>
+where
+ T: SoftDelete,
+ U: FromIterator<T> + IntoIterator<Item = T>,
+{
+ fn default() -> Self {
+ Self::new()
+ }
+}
+
+impl<T, U: Extend<T>> Extend<T> for SoftSearch<U> {
+ fn extend<I: IntoIterator<Item = T>>(&mut self, iter: I) {
+ self.0.extend(iter);
+ }
+}
+
+impl<T, U: FromIterator<T>> FromIterator<T> for SoftSearch<U> {
+ fn from_iter<I: IntoIterator<Item = T>>(iter: I) -> Self {
+ Self(U::from_iter(iter))
+ }
+}
+
+impl<T: IntoIterator> IntoIterator for SoftSearch<T> {
+ type Item = T::Item;
+ type IntoIter = T::IntoIter;
+
+ fn into_iter(self) -> Self::IntoIter {
+ self.0.into_iter()
+ }
+}
+
+impl<K, V, T> NearestNeighbors<K, V> for SoftSearch<T>
+where
+ K: Proximity<V>,
+ V: SoftDelete,
+ T: NearestNeighbors<K, V>,
+{
+ fn search<'k, 'v, N>(&'v self, neighborhood: N) -> N
+ where
+ K: 'k,
+ V: 'v,
+ N: Neighborhood<&'k K, &'v V>
+ {
+ self.0.search(SoftNeighborhood(neighborhood)).0
+ }
+}
+
+/// A k-d forest that supports soft deletes.
+pub type SoftKdForest<T> = SoftSearch<KdForest<T>>;
+
+/// A k-d tree that supports soft deletes.
+pub type SoftKdTree<T> = SoftSearch<FlatKdTree<T>>;
+
+/// A VP forest that supports soft deletes.
+pub type SoftVpForest<T> = SoftSearch<VpForest<T>>;
+
+/// A VP tree that supports soft deletes.
+pub type SoftVpTree<T> = SoftSearch<FlatVpTree<T>>;
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+
+ use acap::coords::Coordinates;
+ use acap::euclid::{euclidean_distance, Euclidean, EuclideanDistance};
+ use acap::Neighbor;
+
+ type Point = Euclidean<[f32; 3]>;
+
+ #[derive(Debug, PartialEq)]
+ struct SoftPoint {
+ point: [f32; 3],
+ deleted: bool,
+ }
+
+ impl SoftPoint {
+ fn new(x: f32, y: f32, z: f32) -> Self {
+ Self {
+ point: [x, y, z],
+ deleted: false,
+ }
+ }
+
+ fn deleted(x: f32, y: f32, z: f32) -> Self {
+ Self {
+ point: [x, y, z],
+ deleted: true,
+ }
+ }
+ }
+
+ impl SoftDelete for SoftPoint {
+ fn is_deleted(&self) -> bool {
+ self.deleted
+ }
+ }
+
+ impl Proximity for SoftPoint {
+ type Distance = EuclideanDistance<f32>;
+
+ fn distance(&self, other: &Self) -> Self::Distance {
+ euclidean_distance(&self.point, &other.point)
+ }
+ }
+
+ impl Coordinates for SoftPoint {
+ type Value = <Point as Coordinates>::Value;
+
+ fn dims(&self) -> usize {
+ self.point.dims()
+ }
+
+ fn coord(&self, i: usize) -> Self::Value {
+ self.point.coord(i)
+ }
+ }
+
+ impl Proximity<SoftPoint> for Point {
+ type Distance = EuclideanDistance<f32>;
+
+ fn distance(&self, other: &SoftPoint) -> Self::Distance {
+ euclidean_distance(&self, &other.point)
+ }
+ }
+
+ fn test_index<T>(index: &T)
+ where
+ T: NearestNeighbors<Point, SoftPoint>,
+ {
+ let target = Euclidean([0.0, 0.0, 0.0]);
+
+ assert_eq!(
+ index.nearest(&target).expect("No nearest neighbor found"),
+ Neighbor::new(&SoftPoint::new(1.0, 2.0, 2.0), 3.0)
+ );
+
+ assert_eq!(index.nearest_within(&target, 2.0), None);
+ assert_eq!(
+ index.nearest_within(&target, 4.0).expect("No nearest neighbor found within 4.0"),
+ Neighbor::new(&SoftPoint::new(1.0, 2.0, 2.0), 3.0)
+ );
+
+ assert_eq!(
+ index.k_nearest(&target, 3),
+ vec![
+ Neighbor::new(&SoftPoint::new(1.0, 2.0, 2.0), 3.0),
+ Neighbor::new(&SoftPoint::new(3.0, 4.0, 0.0), 5.0),
+ Neighbor::new(&SoftPoint::new(2.0, 3.0, 6.0), 7.0),
+ ]
+ );
+
+ assert_eq!(
+ index.k_nearest_within(&target, 3, 6.0),
+ vec![
+ Neighbor::new(&SoftPoint::new(1.0, 2.0, 2.0), 3.0),
+ Neighbor::new(&SoftPoint::new(3.0, 4.0, 0.0), 5.0),
+ ]
+ );
+ assert_eq!(
+ index.k_nearest_within(&target, 3, 8.0),
+ vec![
+ Neighbor::new(&SoftPoint::new(1.0, 2.0, 2.0), 3.0),
+ Neighbor::new(&SoftPoint::new(3.0, 4.0, 0.0), 5.0),
+ Neighbor::new(&SoftPoint::new(2.0, 3.0, 6.0), 7.0),
+ ]
+ );
+ }
+
+ fn test_soft_index<T>(index: &mut SoftSearch<T>)
+ where
+ T: Extend<SoftPoint>,
+ T: FromIterator<SoftPoint>,
+ T: IntoIterator<Item = SoftPoint>,
+ T: NearestNeighbors<Point, SoftPoint>,
+ {
+ let points = vec![
+ SoftPoint::deleted(0.0, 0.0, 0.0),
+ SoftPoint::new(3.0, 4.0, 0.0),
+ SoftPoint::new(5.0, 0.0, 12.0),
+ SoftPoint::new(0.0, 8.0, 15.0),
+ SoftPoint::new(1.0, 2.0, 2.0),
+ SoftPoint::new(2.0, 3.0, 6.0),
+ SoftPoint::new(4.0, 4.0, 7.0),
+ ];
+
+ for point in points {
+ index.push(point);
+ }
+ test_index(index);
+
+ index.rebuild();
+ test_index(index);
+ }
+
+ #[test]
+ fn test_soft_kd_forest() {
+ test_soft_index(&mut SoftKdForest::new());
+ }
+
+ #[test]
+ fn test_soft_vp_forest() {
+ test_soft_index(&mut SoftVpForest::new());
+ }
+}
generated by cgit v1.2.3 (git 2.39.1) at 2025-07-04 13:01:06 +0000