From c653c8cda8f49d3bbe07190a6477367290ff7f04 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sun, 19 Apr 2020 16:28:10 -0400 Subject: Begin re-writing in Rust --- src/main.rs | 1 + 1 file changed, 1 insertion(+) create mode 100644 src/main.rs (limited to 'src') diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..f328e4d --- /dev/null +++ b/src/main.rs @@ -0,0 +1 @@ +fn main() {} -- cgit v1.2.3 From 8aa2911b4c978ad7131a430d07a580cedf6f8f65 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sun, 19 Apr 2020 16:37:16 -0400 Subject: metric: Add some general interfaces for metric spaces --- Cargo.lock | 104 ++++++++++++ Cargo.toml | 4 + src/main.rs | 2 + src/metric.rs | 531 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 641 insertions(+) create mode 100644 src/metric.rs (limited to 'src') diff --git a/Cargo.lock b/Cargo.lock index a63d677..45602e1 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1,5 +1,109 @@ # This file is automatically @generated by Cargo. # It is not intended for manual editing. +[[package]] +name = "autocfg" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f8aac770f1885fd7e387acedd76065302551364496e46b3dd00860b2f8359b9d" + +[[package]] +name = "cfg-if" +version = "0.1.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4785bdd1c96b2a846b2bd7cc02e86b6b3dbf14e7e53446c4f54c92a361040822" + +[[package]] +name = "getrandom" +version = "0.1.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7abc8dd8451921606d809ba32e95b6111925cd2906060d2dcc29c070220503eb" +dependencies = [ + "cfg-if", + "libc", + "wasi", +] + [[package]] name = "kd-forest" version = "2.0.0" +dependencies = [ + "ordered-float", + "rand", +] + +[[package]] +name = "libc" +version = "0.2.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "99e85c08494b21a9054e7fe1374a732aeadaff3980b6990b94bfd3a70f690005" + +[[package]] +name = "num-traits" +version = "0.2.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c62be47e61d1842b9170f0fdeec8eba98e60e90e5446449a0545e5152acd7096" +dependencies = [ + "autocfg", +] + +[[package]] +name = "ordered-float" +version = "1.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "18869315e81473c951eb56ad5558bbc56978562d3ecfb87abb7a1e944cea4518" +dependencies = [ + "num-traits", +] + +[[package]] +name = "ppv-lite86" +version = "0.2.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "74490b50b9fbe561ac330df47c08f3f33073d2d00c150f719147d7c54522fa1b" + +[[package]] +name = "rand" +version = "0.7.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6a6b1679d49b24bbfe0c803429aa1874472f50d9b363131f0e89fc356b544d03" +dependencies = [ + "getrandom", + "libc", + "rand_chacha", + "rand_core", + "rand_hc", +] + +[[package]] +name = "rand_chacha" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f4c8ed856279c9737206bf725bf36935d8666ead7aa69b52be55af369d193402" +dependencies = [ + "ppv-lite86", + "rand_core", +] + +[[package]] +name = "rand_core" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "90bde5296fc891b0cef12a6d03ddccc162ce7b2aff54160af9338f8d40df6d19" +dependencies = [ + "getrandom", +] + +[[package]] +name = "rand_hc" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ca3129af7b92a17112d59ad498c6f81eaf463253766b90396d39ea7a39d6613c" +dependencies = [ + "rand_core", +] + +[[package]] +name = "wasi" +version = "0.9.0+wasi-snapshot-preview1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cccddf32554fecc6acb585f82a32a72e28b48f8c4c1883ddfeeeaa96f7d8e519" diff --git a/Cargo.toml b/Cargo.toml index 65ffec8..5a93cf5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,3 +3,7 @@ name = "kd-forest" version = "2.0.0" authors = ["Tavian Barnes "] edition = "2018" + +[dependencies] +ordered-float = "1.0.2" +rand = "0.7.3" diff --git a/src/main.rs b/src/main.rs index f328e4d..0d7989b 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1 +1,3 @@ +pub mod metric; + fn main() {} diff --git a/src/metric.rs b/src/metric.rs new file mode 100644 index 0000000..e067771 --- /dev/null +++ b/src/metric.rs @@ -0,0 +1,531 @@ +//! [Metric spaces](https://en.wikipedia.org/wiki/Metric_space). + +use ordered_float::OrderedFloat; + +use std::cmp::Ordering; +use std::collections::BinaryHeap; +use std::iter::FromIterator; + +/// An [order embedding](https://en.wikipedia.org/wiki/Order_embedding) for distances. +/// +/// Implementations of this trait must satisfy, for all non-negative distances `x` and `y`: +/// +/// * `x == Self::from(x).into()` +/// * `x <= y` iff `Self::from(x) <= Self::from(y)` +/// +/// This trait exists to optimize the common case where distances can be compared more efficiently +/// than their exact values can be computed. For example, taking the square root can be avoided +/// when comparing Euclidean distances (see [SquaredDistance]). +pub trait Distance: Copy + From + Into + Ord {} + +/// A raw numerical distance. +#[derive(Debug, Clone, Copy, Eq, Ord, PartialEq, PartialOrd)] +pub struct RawDistance(OrderedFloat); + +impl From for RawDistance { + fn from(value: f64) -> Self { + Self(value.into()) + } +} + +impl From for f64 { + fn from(value: RawDistance) -> Self { + value.0.into_inner() + } +} + +impl Distance for RawDistance {} + +/// A squared distance, to avoid computing square roots unless absolutely necessary. +#[derive(Debug, Clone, Copy, Eq, Ord, PartialEq, PartialOrd)] +pub struct SquaredDistance(OrderedFloat); + +impl SquaredDistance { + /// Create a SquaredDistance from an already squared value. + pub fn from_squared(value: f64) -> Self { + Self(value.into()) + } +} + +impl From for SquaredDistance { + fn from(value: f64) -> Self { + Self::from_squared(value * value) + } +} + +impl From for f64 { + fn from(value: SquaredDistance) -> Self { + value.0.into_inner().sqrt() + } +} + +impl Distance for SquaredDistance {} + +/// A [metric space](https://en.wikipedia.org/wiki/Metric_space). +pub trait Metric { + /// The type used to represent distances. Use [RawDistance] to compare the actual values + /// directly, or another type if comparisons can be implemented more efficiently. + type Distance: Distance; + + /// Computes the distance between this point and another point. This function must satisfy + /// three conditions: + /// + /// * `x.distance(y) == 0` iff `x == y` (identity of indiscernibles) + /// * `x.distance(y) == y.distance(x)` (symmetry) + /// * `x.distance(z) <= x.distance(y) + y.distance(z)` (triangle inequality) + fn distance(&self, other: &T) -> Self::Distance; +} + +/// Blanket [Metric] implementation for references. +impl<'a, 'b, T, U: Metric> Metric<&'a T> for &'b U { + type Distance = U::Distance; + + fn distance(&self, other: &&'a T) -> Self::Distance { + (*self).distance(other) + } +} + +/// The standard [Euclidean distance](https://en.wikipedia.org/wiki/Euclidean_distance) metric. +impl Metric for [f64] { + type Distance = SquaredDistance; + + fn distance(&self, other: &Self) -> Self::Distance { + debug_assert!(self.len() == other.len()); + + let mut sum = 0.0; + for i in 0..self.len() { + let diff = self[i] - other[i]; + sum += diff * diff; + } + + Self::Distance::from_squared(sum) + } +} + +/// A nearest neighbor to a target. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Neighbor { + /// The found item. + pub item: T, + /// The distance from the target. + pub distance: f64, +} + +impl Neighbor { + /// Create a new Neighbor. + pub fn new(item: T, distance: f64) -> Self { + Self { item, distance } + } +} + +/// A candidate nearest neighbor found during a search. +#[derive(Debug)] +struct Candidate { + item: T, + distance: D, +} + +impl Candidate { + fn new(target: U, item: T) -> Self + where + U: Metric, + { + let distance = target.distance(&item); + Self { item, distance } + } + + fn into_neighbor(self) -> Neighbor { + Neighbor::new(self.item, self.distance.into()) + } +} + +impl PartialOrd for Candidate { + fn partial_cmp(&self, other: &Self) -> Option { + self.distance.partial_cmp(&other.distance) + } +} + +impl Ord for Candidate { + fn cmp(&self, other: &Self) -> Ordering { + self.distance.cmp(&other.distance) + } +} + +impl PartialEq for Candidate { + fn eq(&self, other: &Self) -> bool { + self.distance.eq(&other.distance) + } +} + +impl Eq for Candidate {} + +/// Accumulates nearest neighbor search results. +pub trait Neighborhood> { + /// Returns the target of the nearest neighbor search. + fn target(&self) -> U; + + /// Check whether a distance is within this neighborhood. + fn contains(&self, distance: f64) -> bool { + distance < 0.0 || self.contains_distance(distance.into()) + } + + /// Check whether a distance is within this neighborhood. + fn contains_distance(&self, distance: U::Distance) -> bool; + + /// Consider a new candidate neighbor. + fn consider(&mut self, item: T) -> U::Distance; +} + +/// A [Neighborhood] with at most one result. +#[derive(Debug)] +struct SingletonNeighborhood> { + /// The target of the nearest neighbor search. + target: U, + /// The current threshold distance to the farthest result. + threshold: Option, + /// The current nearest neighbor, if any. + candidate: Option>, +} + +impl SingletonNeighborhood +where + U: Copy + Metric, +{ + /// Create a new single metric result tracker. + /// + /// * `target`: The target fo the nearest neighbor search. + /// * `threshold`: The maximum allowable distance. + fn new(target: U, threshold: Option) -> Self { + Self { + target, + threshold: threshold.map(U::Distance::from), + candidate: None, + } + } + + /// Consider a candidate. + fn push(&mut self, candidate: Candidate) -> U::Distance { + let distance = candidate.distance; + + if self.contains_distance(distance) { + self.threshold = Some(distance); + self.candidate = Some(candidate); + } + + distance + } + + /// Convert this result into an optional neighbor. + fn into_option(self) -> Option> { + self.candidate.map(Candidate::into_neighbor) + } +} + +impl Neighborhood for SingletonNeighborhood +where + U: Copy + Metric, +{ + fn target(&self) -> U { + self.target + } + + fn contains_distance(&self, distance: U::Distance) -> bool { + self.threshold.map(|t| distance <= t).unwrap_or(true) + } + + fn consider(&mut self, item: T) -> U::Distance { + self.push(Candidate::new(self.target, item)) + } +} + +/// A [Neighborhood] of up to `k` results, using a binary heap. +#[derive(Debug)] +struct HeapNeighborhood> { + /// The target of the nearest neighbor search. + target: U, + /// The number of nearest neighbors to find. + k: usize, + /// The current threshold distance to the farthest result. + threshold: Option, + /// A max-heap of the best candidates found so far. + heap: BinaryHeap>, +} + +impl HeapNeighborhood +where + U: Copy + Metric, +{ + /// Create a new metric result tracker. + /// + /// * `target`: The target fo the nearest neighbor search. + /// * `k`: The number of nearest neighbors to find. + /// * `threshold`: The maximum allowable distance. + fn new(target: U, k: usize, threshold: Option) -> Self { + Self { + target, + k, + threshold: threshold.map(U::Distance::from), + heap: BinaryHeap::with_capacity(k), + } + } + + /// Consider a candidate. + fn push(&mut self, candidate: Candidate) -> U::Distance { + let distance = candidate.distance; + + if self.contains_distance(distance) { + let heap = &mut self.heap; + + if heap.len() == self.k { + heap.pop(); + } + + heap.push(candidate); + + if heap.len() == self.k { + self.threshold = self.heap.peek().map(|c| c.distance) + } + } + + distance + } + + /// Convert these results into a vector of neighbors. + fn into_vec(self) -> Vec> { + self.heap + .into_sorted_vec() + .into_iter() + .map(Candidate::into_neighbor) + .collect() + } +} + +impl Neighborhood for HeapNeighborhood +where + U: Copy + Metric, +{ + fn target(&self) -> U { + self.target + } + + fn contains_distance(&self, distance: U::Distance) -> bool { + self.k > 0 && self.threshold.map(|t| distance <= t).unwrap_or(true) + } + + fn consider(&mut self, item: T) -> U::Distance { + self.push(Candidate::new(self.target, item)) + } +} + +/// A [nearest neighbor search](https://en.wikipedia.org/wiki/Nearest_neighbor_search) index. +/// +/// Type parameters: +/// * `T`: The search result type. +/// * `U`: The query type. +pub trait NearestNeighbors = T> { + /// Returns the nearest neighbor to `target` (or `None` if this index is empty). + fn nearest(&self, target: &U) -> Option> { + self.search(SingletonNeighborhood::new(target, None)) + .into_option() + } + + /// Returns the nearest neighbor to `target` within the distance `threshold`, if one exists. + fn nearest_within(&self, target: &U, threshold: f64) -> Option> { + self.search(SingletonNeighborhood::new(target, Some(threshold))) + .into_option() + } + + /// Returns the up to `k` nearest neighbors to `target`. + fn k_nearest(&self, target: &U, k: usize) -> Vec> { + self.search(HeapNeighborhood::new(target, k, None)) + .into_vec() + } + + /// Returns the up to `k` nearest neighbors to `target` within the distance `threshold`. + fn k_nearest_within(&self, target: &U, k: usize, threshold: f64) -> Vec> { + self.search(HeapNeighborhood::new(target, k, Some(threshold))) + .into_vec() + } + + /// Search for nearest neighbors and add them to a neighborhood. + fn search<'a, 'b, N>(&'a self, neighborhood: N) -> N + where + T: 'a, + U: 'b, + N: Neighborhood<&'a T, &'b U>; +} + +/// A [NearestNeighbors] implementation that does exhaustive search. +#[derive(Debug)] +pub struct ExhaustiveSearch(Vec); + +impl ExhaustiveSearch { + /// Create an empty ExhaustiveSearch index. + pub fn new() -> Self { + Self(Vec::new()) + } + + /// Add a new item to the index. + pub fn push(&mut self, item: T) { + self.0.push(item); + } +} + +impl FromIterator for ExhaustiveSearch { + fn from_iter>(items: I) -> Self { + Self(items.into_iter().collect()) + } +} + +impl IntoIterator for ExhaustiveSearch { + type Item = T; + type IntoIter = std::vec::IntoIter; + + fn into_iter(self) -> Self::IntoIter { + self.0.into_iter() + } +} + +impl Extend for ExhaustiveSearch { + fn extend>(&mut self, iter: I) { + for value in iter { + self.push(value); + } + } +} + +impl> NearestNeighbors for ExhaustiveSearch { + fn search<'a, 'b, N>(&'a self, mut neighborhood: N) -> N + where + T: 'a, + U: 'b, + N: Neighborhood<&'a T, &'b U>, + { + for e in &self.0 { + neighborhood.consider(e); + } + neighborhood + } +} + +#[cfg(test)] +pub mod tests { + use super::*; + + use rand::prelude::*; + + #[derive(Clone, Copy, Debug, PartialEq)] + pub struct Point(pub [f64; 3]); + + impl Metric for Point { + type Distance = SquaredDistance; + + fn distance(&self, other: &Self) -> Self::Distance { + self.0.distance(&other.0) + } + } + + /// Test a [NearestNeighbors] impl. + pub fn test_nearest_neighbors(from_iter: F) + where + T: NearestNeighbors, + F: Fn(Vec) -> T, + { + test_empty(&from_iter); + test_pythagorean(&from_iter); + test_random_points(&from_iter); + } + + fn test_empty(from_iter: &F) + where + T: NearestNeighbors, + F: Fn(Vec) -> T, + { + let points = Vec::new(); + let index = from_iter(points); + let target = Point([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(from_iter: &F) + where + T: NearestNeighbors, + F: Fn(Vec) -> T, + { + let points = vec![ + Point([3.0, 4.0, 0.0]), + Point([5.0, 0.0, 12.0]), + Point([0.0, 8.0, 15.0]), + Point([1.0, 2.0, 2.0]), + Point([2.0, 3.0, 6.0]), + Point([4.0, 4.0, 7.0]), + ]; + let index = from_iter(points); + let target = Point([0.0, 0.0, 0.0]); + + assert_eq!( + index.nearest(&target), + Some(Neighbor::new(&Point([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), + Some(Neighbor::new(&Point([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(&Point([1.0, 2.0, 2.0]), 3.0), + Neighbor::new(&Point([3.0, 4.0, 0.0]), 5.0), + Neighbor::new(&Point([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(&Point([1.0, 2.0, 2.0]), 3.0), + Neighbor::new(&Point([3.0, 4.0, 0.0]), 5.0), + ] + ); + assert_eq!( + index.k_nearest_within(&target, 3, 8.0), + vec![ + Neighbor::new(&Point([1.0, 2.0, 2.0]), 3.0), + Neighbor::new(&Point([3.0, 4.0, 0.0]), 5.0), + Neighbor::new(&Point([2.0, 3.0, 6.0]), 7.0), + ] + ); + } + + fn test_random_points(from_iter: &F) + where + T: NearestNeighbors, + F: Fn(Vec) -> T, + { + let mut points = Vec::new(); + for _ in 0..255 { + points.push(Point([random(), random(), random()])); + } + let target = Point([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] + fn test_exhaustive_index() { + test_nearest_neighbors(ExhaustiveSearch::from_iter); + } +} -- cgit v1.2.3 From 8d9de0e1028daed981246174182a39dd917b72bc Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sun, 19 Apr 2020 16:38:24 -0400 Subject: metric/vp: Implement vantage-point trees --- src/metric.rs | 2 + src/metric/vp.rs | 168 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 170 insertions(+) create mode 100644 src/metric/vp.rs (limited to 'src') diff --git a/src/metric.rs b/src/metric.rs index e067771..549db67 100644 --- a/src/metric.rs +++ b/src/metric.rs @@ -1,5 +1,7 @@ //! [Metric spaces](https://en.wikipedia.org/wiki/Metric_space). +pub mod vp; + use ordered_float::OrderedFloat; use std::cmp::Ordering; diff --git a/src/metric/vp.rs b/src/metric/vp.rs new file mode 100644 index 0000000..8d5b091 --- /dev/null +++ b/src/metric/vp.rs @@ -0,0 +1,168 @@ +//! [Vantage-point trees](https://en.wikipedia.org/wiki/Vantage-point_tree). + +use super::{Metric, NearestNeighbors, Neighborhood}; + +use std::iter::FromIterator; + +/// A node in a VP tree. +#[derive(Debug)] +struct VpNode { + /// The vantage point itself. + item: T, + /// The radius of this node. + radius: f64, + /// The subtree inside the radius, if any. + inside: Option>, + /// The subtree outside the radius, if any. + outside: Option>, +} + +impl VpNode { + /// Create a new VpNode. + fn new(mut items: Vec) -> Option> { + if items.is_empty() { + return None; + } + + let item = items.pop().unwrap(); + + items.sort_by_cached_key(|a| item.distance(a)); + + let mid = items.len() / 2; + let outside: Vec = items.drain(mid..).collect(); + + let radius = items.last().map(|l| item.distance(l).into()).unwrap_or(0.0); + + Some(Box::new(Self { + item, + radius, + inside: Self::new(items), + outside: Self::new(outside), + })) + } +} + +trait VpSearch<'a, T, U, N> { + /// Recursively search for nearest neighbors. + fn search(&'a self, neighborhood: &mut N); + + /// Search the inside subtree. + fn search_inside(&'a self, distance: f64, neighborhood: &mut N); + + /// Search the outside subtree. + fn search_outside(&'a self, distance: f64, neighborhood: &mut N); +} + +impl<'a, T, U, N> VpSearch<'a, T, U, N> for VpNode +where + T: 'a, + U: Metric<&'a T>, + N: Neighborhood<&'a T, U>, +{ + fn search(&'a self, neighborhood: &mut N) { + let distance = neighborhood.consider(&self.item).into(); + + if distance <= self.radius { + self.search_inside(distance, neighborhood); + self.search_outside(distance, neighborhood); + } else { + self.search_outside(distance, neighborhood); + self.search_inside(distance, neighborhood); + } + } + + fn search_inside(&'a self, distance: f64, neighborhood: &mut N) { + if let Some(inside) = &self.inside { + if neighborhood.contains(distance - self.radius) { + inside.search(neighborhood); + } + } + } + + fn search_outside(&'a self, distance: f64, neighborhood: &mut N) { + if let Some(outside) = &self.outside { + if neighborhood.contains(self.radius - distance) { + outside.search(neighborhood); + } + } + } +} + +/// A [vantage-point tree](https://en.wikipedia.org/wiki/Vantage-point_tree). +#[derive(Debug)] +pub struct VpTree { + root: Option>>, +} + +impl FromIterator for VpTree { + fn from_iter>(items: I) -> Self { + Self { + root: VpNode::new(items.into_iter().collect::>()), + } + } +} + +impl NearestNeighbors for VpTree +where + T: Metric, + U: Metric, +{ + fn search<'a, 'b, N>(&'a self, mut neighborhood: N) -> N + where + T: 'a, + U: 'b, + N: Neighborhood<&'a T, &'b U>, + { + if let Some(root) = &self.root { + root.search(&mut neighborhood); + } + neighborhood + } +} + +/// An iterator that moves values out of a VP tree. +#[derive(Debug)] +pub struct IntoIter { + stack: Vec>>, +} + +impl IntoIter { + fn new(node: Option>>) -> Self { + Self { + stack: node.into_iter().collect(), + } + } +} + +impl Iterator for IntoIter { + type Item = T; + + fn next(&mut self) -> Option { + self.stack.pop().map(|node| { + self.stack.extend(node.inside); + self.stack.extend(node.outside); + node.item + }) + } +} + +impl IntoIterator for VpTree { + type Item = T; + type IntoIter = IntoIter; + + fn into_iter(self) -> Self::IntoIter { + IntoIter::new(self.root) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + use crate::metric::tests::test_nearest_neighbors; + + #[test] + fn test_vp_tree() { + test_nearest_neighbors(VpTree::from_iter); + } +} -- cgit v1.2.3 From 1c560791902a4ef72efa671106d8f6d97fea50c1 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sun, 19 Apr 2020 16:40:29 -0400 Subject: metric/kd: Implement k-d trees --- src/metric.rs | 1 + src/metric/kd.rs | 217 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 218 insertions(+) create mode 100644 src/metric/kd.rs (limited to 'src') diff --git a/src/metric.rs b/src/metric.rs index 549db67..95191da 100644 --- a/src/metric.rs +++ b/src/metric.rs @@ -1,5 +1,6 @@ //! [Metric spaces](https://en.wikipedia.org/wiki/Metric_space). +pub mod kd; pub mod vp; use ordered_float::OrderedFloat; diff --git a/src/metric/kd.rs b/src/metric/kd.rs new file mode 100644 index 0000000..ab0cd2e --- /dev/null +++ b/src/metric/kd.rs @@ -0,0 +1,217 @@ +//! [k-d trees](https://en.wikipedia.org/wiki/K-d_tree). + +use super::{Metric, NearestNeighbors, Neighborhood}; + +use ordered_float::OrderedFloat; + +use std::iter::FromIterator; + +/// A point in Cartesian space. +pub trait Cartesian { + /// Returns the number of dimensions necessary to describe this point. + fn dimensions(&self) -> usize; + + /// Returns the value of the `i`th coordinate of this point (`i < self.dimensions()`). + fn coordinate(&self, i: usize) -> f64; +} + +/// Blanket [Cartesian] implementation for references. +impl<'a, T: Cartesian> Cartesian for &'a T { + fn dimensions(&self) -> usize { + (*self).dimensions() + } + + fn coordinate(&self, i: usize) -> f64 { + (*self).coordinate(i) + } +} + +/// Standard cartesian space. +impl Cartesian for [f64] { + fn dimensions(&self) -> usize { + self.len() + } + + fn coordinate(&self, i: usize) -> f64 { + self[i] + } +} + +/// A node in a k-d tree. +#[derive(Debug)] +struct KdNode { + /// The value stored in this node. + item: T, + /// The left subtree, if any. + left: Option>, + /// The right subtree, if any. + right: Option>, +} + +trait KdSearch<'a, T, U, N> { + /// Recursively search for nearest neighbors. + fn search(&'a self, i: usize, neighborhood: &mut N); + + /// Search the left subtree. + fn search_left(&'a self, i: usize, distance: f64, neighborhood: &mut N); + + /// Search the right subtree. + fn search_right(&'a self, i: usize, distance: f64, neighborhood: &mut N); +} + +impl<'a, T, U, N> KdSearch<'a, T, U, N> for KdNode +where + T: 'a + Cartesian, + U: Cartesian + Metric<&'a T>, + N: Neighborhood<&'a T, U>, +{ + fn search(&'a self, i: usize, neighborhood: &mut N) { + neighborhood.consider(&self.item); + + let distance = neighborhood.target().coordinate(i) - self.item.coordinate(i); + let j = (i + 1) % self.item.dimensions(); + if distance <= 0.0 { + self.search_left(j, distance, neighborhood); + self.search_right(j, -distance, neighborhood); + } else { + self.search_right(j, -distance, neighborhood); + self.search_left(j, distance, neighborhood); + } + } + + fn search_left(&'a self, i: usize, distance: f64, neighborhood: &mut N) { + if let Some(left) = &self.left { + if neighborhood.contains(distance) { + left.search(i, neighborhood); + } + } + } + + fn search_right(&'a self, i: usize, distance: f64, neighborhood: &mut N) { + if let Some(right) = &self.right { + if neighborhood.contains(distance) { + right.search(i, neighborhood); + } + } + } +} + +impl KdNode { + /// Create a new KdNode. + fn new(i: usize, mut items: Vec) -> Option> { + if items.is_empty() { + return None; + } + + items.sort_unstable_by_key(|x| OrderedFloat::from(x.coordinate(i))); + + let mid = items.len() / 2; + let right: Vec = items.drain((mid + 1)..).collect(); + let item = items.pop().unwrap(); + let j = (i + 1) % item.dimensions(); + Some(Box::new(Self { + item, + left: Self::new(j, items), + right: Self::new(j, right), + })) + } +} + +/// A [k-d tree](https://en.wikipedia.org/wiki/K-d_tree). +#[derive(Debug)] +pub struct KdTree { + root: Option>>, +} + +impl FromIterator for KdTree { + /// Create a new k-d tree from a set of points. + fn from_iter>(items: I) -> Self { + Self { + root: KdNode::new(0, items.into_iter().collect()), + } + } +} + +impl NearestNeighbors for KdTree +where + T: Cartesian, + U: Cartesian + Metric, +{ + fn search<'a, 'b, N>(&'a self, mut neighborhood: N) -> N + where + T: 'a, + U: 'b, + N: Neighborhood<&'a T, &'b U>, + { + if let Some(root) = &self.root { + root.search(0, &mut neighborhood); + } + neighborhood + } +} + +/// An iterator that the moves values out of a k-d tree. +#[derive(Debug)] +pub struct IntoIter { + stack: Vec>>, +} + +impl IntoIter { + fn new(node: Option>>) -> Self { + Self { + stack: node.into_iter().collect(), + } + } +} + +impl Iterator for IntoIter { + type Item = T; + + fn next(&mut self) -> Option { + self.stack.pop().map(|node| { + self.stack.extend(node.left); + self.stack.extend(node.right); + node.item + }) + } +} + +impl IntoIterator for KdTree { + type Item = T; + type IntoIter = IntoIter; + + fn into_iter(self) -> Self::IntoIter { + IntoIter::new(self.root) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + use crate::metric::tests::{test_nearest_neighbors, Point}; + use crate::metric::SquaredDistance; + + impl Metric<[f64]> for Point { + type Distance = SquaredDistance; + + fn distance(&self, other: &[f64]) -> Self::Distance { + self.0.distance(other) + } + } + + impl Cartesian for Point { + fn dimensions(&self) -> usize { + self.0.dimensions() + } + + fn coordinate(&self, i: usize) -> f64 { + self.0.coordinate(i) + } + } + + #[test] + fn test_kd_tree() { + test_nearest_neighbors(KdTree::from_iter); + } +} -- cgit v1.2.3 From 5de377b2b00a927a4f6463c1c5a5fd18606ad006 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Thu, 30 Apr 2020 22:51:06 -0400 Subject: metric/kd: Prune k-d tree searches more aggressively --- src/metric/kd.rs | 116 +++++++++++++++++++++++++++++++------------------------ 1 file changed, 65 insertions(+), 51 deletions(-) (limited to 'src') diff --git a/src/metric/kd.rs b/src/metric/kd.rs index ab0cd2e..db1b2bd 100644 --- a/src/metric/kd.rs +++ b/src/metric/kd.rs @@ -7,7 +7,7 @@ use ordered_float::OrderedFloat; use std::iter::FromIterator; /// A point in Cartesian space. -pub trait Cartesian { +pub trait Cartesian: Metric<[f64]> { /// Returns the number of dimensions necessary to describe this point. fn dimensions(&self) -> usize; @@ -26,6 +26,15 @@ impl<'a, T: Cartesian> Cartesian for &'a T { } } +/// Blanket [Metric<[f64]>](Metric) implementation for [Cartesian] references. +impl<'a, T: Cartesian> Metric<[f64]> for &'a T { + type Distance = T::Distance; + + fn distance(&self, other: &[f64]) -> Self::Distance { + (*self).distance(other) + } +} + /// Standard cartesian space. impl Cartesian for [f64] { fn dimensions(&self) -> usize { @@ -37,6 +46,21 @@ impl Cartesian for [f64] { } } +/// Marker trait for cartesian metric spaces. +pub trait CartesianMetric: + Cartesian + Metric>::Distance> +{ +} + +/// Blanket [CartesianMetric] implementation for cartesian spaces with compatible metric distance +/// types. +impl CartesianMetric for U +where + T: ?Sized, + U: ?Sized + Cartesian + Metric>::Distance>, +{ +} + /// A node in a k-d tree. #[derive(Debug)] struct KdNode { @@ -48,54 +72,6 @@ struct KdNode { right: Option>, } -trait KdSearch<'a, T, U, N> { - /// Recursively search for nearest neighbors. - fn search(&'a self, i: usize, neighborhood: &mut N); - - /// Search the left subtree. - fn search_left(&'a self, i: usize, distance: f64, neighborhood: &mut N); - - /// Search the right subtree. - fn search_right(&'a self, i: usize, distance: f64, neighborhood: &mut N); -} - -impl<'a, T, U, N> KdSearch<'a, T, U, N> for KdNode -where - T: 'a + Cartesian, - U: Cartesian + Metric<&'a T>, - N: Neighborhood<&'a T, U>, -{ - fn search(&'a self, i: usize, neighborhood: &mut N) { - neighborhood.consider(&self.item); - - let distance = neighborhood.target().coordinate(i) - self.item.coordinate(i); - let j = (i + 1) % self.item.dimensions(); - if distance <= 0.0 { - self.search_left(j, distance, neighborhood); - self.search_right(j, -distance, neighborhood); - } else { - self.search_right(j, -distance, neighborhood); - self.search_left(j, distance, neighborhood); - } - } - - fn search_left(&'a self, i: usize, distance: f64, neighborhood: &mut N) { - if let Some(left) = &self.left { - if neighborhood.contains(distance) { - left.search(i, neighborhood); - } - } - } - - fn search_right(&'a self, i: usize, distance: f64, neighborhood: &mut N) { - if let Some(right) = &self.right { - if neighborhood.contains(distance) { - right.search(i, neighborhood); - } - } - } -} - impl KdNode { /// Create a new KdNode. fn new(i: usize, mut items: Vec) -> Option> { @@ -115,6 +91,40 @@ impl KdNode { right: Self::new(j, right), })) } + + /// Recursively search for nearest neighbors. + fn search<'a, U, N>(&'a self, i: usize, closest: &mut [f64], neighborhood: &mut N) + where + T: 'a, + U: CartesianMetric<&'a T>, + N: Neighborhood<&'a T, U>, + { + neighborhood.consider(&self.item); + + let target = neighborhood.target(); + let ti = target.coordinate(i); + let si = self.item.coordinate(i); + let j = (i + 1) % self.item.dimensions(); + + let (near, far) = if ti <= si { + (&self.left, &self.right) + } else { + (&self.right, &self.left) + }; + + if let Some(near) = near { + near.search(j, closest, neighborhood); + } + + if let Some(far) = far { + let saved = closest[i]; + closest[i] = si; + if neighborhood.contains_distance(target.distance(closest)) { + far.search(j, closest, neighborhood); + } + closest[i] = saved; + } + } } /// A [k-d tree](https://en.wikipedia.org/wiki/K-d_tree). @@ -135,7 +145,7 @@ impl FromIterator for KdTree { impl NearestNeighbors for KdTree where T: Cartesian, - U: Cartesian + Metric, + U: CartesianMetric, { fn search<'a, 'b, N>(&'a self, mut neighborhood: N) -> N where @@ -143,8 +153,12 @@ where U: 'b, N: Neighborhood<&'a T, &'b U>, { + let target = neighborhood.target(); + let dims = target.dimensions(); + let mut closest: Vec<_> = (0..dims).map(|i| target.coordinate(i)).collect(); + if let Some(root) = &self.root { - root.search(0, &mut neighborhood); + root.search(0, &mut closest, &mut neighborhood); } neighborhood } -- cgit v1.2.3 From 9699f4657ecaaf4361448f249e4f2e210a854af4 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Thu, 23 Apr 2020 09:55:26 -0400 Subject: metric/vp: Flatten the tree representation --- src/metric/vp.rs | 133 +++++++++++++++++++++---------------------------------- 1 file changed, 51 insertions(+), 82 deletions(-) (limited to 'src') diff --git a/src/metric/vp.rs b/src/metric/vp.rs index 8d5b091..fae62e5 100644 --- a/src/metric/vp.rs +++ b/src/metric/vp.rs @@ -11,78 +11,62 @@ struct VpNode { item: T, /// The radius of this node. radius: f64, - /// The subtree inside the radius, if any. - inside: Option>, - /// The subtree outside the radius, if any. - outside: Option>, + /// The size of the subtree inside the radius. + inside_len: usize, } impl VpNode { /// Create a new VpNode. - fn new(mut items: Vec) -> Option> { - if items.is_empty() { - return None; + fn new(item: T) -> Self { + Self { + item, + radius: 0.0, + inside_len: 0, } + } - let item = items.pop().unwrap(); - - items.sort_by_cached_key(|a| item.distance(a)); - - let mid = items.len() / 2; - let outside: Vec = items.drain(mid..).collect(); + /// Build a VP tree recursively. + fn build(slice: &mut [VpNode]) { + if let Some((node, children)) = slice.split_first_mut() { + let item = &node.item; + children.sort_by_cached_key(|n| item.distance(&n.item)); - let radius = items.last().map(|l| item.distance(l).into()).unwrap_or(0.0); + let (inside, outside) = children.split_at_mut(children.len() / 2); + if let Some(last) = inside.last() { + node.radius = item.distance(&last.item).into(); + } + node.inside_len = inside.len(); - Some(Box::new(Self { - item, - radius, - inside: Self::new(items), - outside: Self::new(outside), - })) + Self::build(inside); + Self::build(outside); + } } -} -trait VpSearch<'a, T, U, N> { /// Recursively search for nearest neighbors. - fn search(&'a self, neighborhood: &mut N); - - /// Search the inside subtree. - fn search_inside(&'a self, distance: f64, neighborhood: &mut N); - - /// Search the outside subtree. - fn search_outside(&'a self, distance: f64, neighborhood: &mut N); -} + fn recurse<'a, U, N>(slice: &'a [VpNode], neighborhood: &mut N) + where + T: 'a, + U: Metric<&'a T>, + N: Neighborhood<&'a T, U>, + { + let (node, children) = slice.split_first().unwrap(); + let (inside, outside) = children.split_at(node.inside_len); -impl<'a, T, U, N> VpSearch<'a, T, U, N> for VpNode -where - T: 'a, - U: Metric<&'a T>, - N: Neighborhood<&'a T, U>, -{ - fn search(&'a self, neighborhood: &mut N) { - let distance = neighborhood.consider(&self.item).into(); + let distance = neighborhood.consider(&node.item).into(); - if distance <= self.radius { - self.search_inside(distance, neighborhood); - self.search_outside(distance, neighborhood); + if distance <= node.radius { + if !inside.is_empty() && neighborhood.contains(distance - node.radius) { + Self::recurse(inside, neighborhood); + } + if !outside.is_empty() && neighborhood.contains(node.radius - distance) { + Self::recurse(outside, neighborhood); + } } else { - self.search_outside(distance, neighborhood); - self.search_inside(distance, neighborhood); - } - } - - fn search_inside(&'a self, distance: f64, neighborhood: &mut N) { - if let Some(inside) = &self.inside { - if neighborhood.contains(distance - self.radius) { - inside.search(neighborhood); + if !outside.is_empty() && neighborhood.contains(node.radius - distance) { + Self::recurse(outside, neighborhood); } - } - } - - fn search_outside(&'a self, distance: f64, neighborhood: &mut N) { - if let Some(outside) = &self.outside { - if neighborhood.contains(self.radius - distance) { - outside.search(neighborhood); + if !inside.is_empty() && neighborhood.contains(distance - node.radius) { + Self::recurse(inside, neighborhood); } } } @@ -90,15 +74,13 @@ where /// A [vantage-point tree](https://en.wikipedia.org/wiki/Vantage-point_tree). #[derive(Debug)] -pub struct VpTree { - root: Option>>, -} +pub struct VpTree(Vec>); impl FromIterator for VpTree { fn from_iter>(items: I) -> Self { - Self { - root: VpNode::new(items.into_iter().collect::>()), - } + let mut nodes: Vec<_> = items.into_iter().map(VpNode::new).collect(); + VpNode::build(nodes.as_mut_slice()); + Self(nodes) } } @@ -113,36 +95,23 @@ where U: 'b, N: Neighborhood<&'a T, &'b U>, { - if let Some(root) = &self.root { - root.search(&mut neighborhood); + if !self.0.is_empty() { + VpNode::recurse(&self.0, &mut neighborhood); } + neighborhood } } /// An iterator that moves values out of a VP tree. #[derive(Debug)] -pub struct IntoIter { - stack: Vec>>, -} - -impl IntoIter { - fn new(node: Option>>) -> Self { - Self { - stack: node.into_iter().collect(), - } - } -} +pub struct IntoIter(std::vec::IntoIter>); impl Iterator for IntoIter { type Item = T; fn next(&mut self) -> Option { - self.stack.pop().map(|node| { - self.stack.extend(node.inside); - self.stack.extend(node.outside); - node.item - }) + self.0.next().map(|n| n.item) } } @@ -151,7 +120,7 @@ impl IntoIterator for VpTree { type IntoIter = IntoIter; fn into_iter(self) -> Self::IntoIter { - IntoIter::new(self.root) + IntoIter(self.0.into_iter()) } } -- cgit v1.2.3 From e9a81a6d0df149252164003975addf175d5c6f4b Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Thu, 23 Apr 2020 09:55:13 -0400 Subject: metric/kd: Flatten the tree representation --- src/metric/kd.rs | 113 ++++++++++++++++++++++++++----------------------------- 1 file changed, 54 insertions(+), 59 deletions(-) (limited to 'src') diff --git a/src/metric/kd.rs b/src/metric/kd.rs index db1b2bd..2caf4a3 100644 --- a/src/metric/kd.rs +++ b/src/metric/kd.rs @@ -66,61 +66,71 @@ where struct KdNode { /// The value stored in this node. item: T, - /// The left subtree, if any. - left: Option>, - /// The right subtree, if any. - right: Option>, + /// The size of the left subtree. + left_len: usize, } impl KdNode { /// Create a new KdNode. - fn new(i: usize, mut items: Vec) -> Option> { - if items.is_empty() { - return None; + fn new(item: T) -> Self { + Self { item, left_len: 0 } + } + + /// Build a k-d tree recursively. + fn build(slice: &mut [KdNode], i: usize) { + if slice.is_empty() { + return; } - items.sort_unstable_by_key(|x| OrderedFloat::from(x.coordinate(i))); - - let mid = items.len() / 2; - let right: Vec = items.drain((mid + 1)..).collect(); - let item = items.pop().unwrap(); - let j = (i + 1) % item.dimensions(); - Some(Box::new(Self { - item, - left: Self::new(j, items), - right: Self::new(j, right), - })) + slice.sort_unstable_by_key(|n| OrderedFloat::from(n.item.coordinate(i))); + + let mid = slice.len() / 2; + slice.swap(0, mid); + + let (node, children) = slice.split_first_mut().unwrap(); + let (left, right) = children.split_at_mut(mid); + node.left_len = left.len(); + + let j = (i + 1) % node.item.dimensions(); + Self::build(left, j); + Self::build(right, j); } /// Recursively search for nearest neighbors. - fn search<'a, U, N>(&'a self, i: usize, closest: &mut [f64], neighborhood: &mut N) - where + fn recurse<'a, U, N>( + slice: &'a [KdNode], + i: usize, + closest: &mut [f64], + neighborhood: &mut N, + ) where T: 'a, U: CartesianMetric<&'a T>, N: Neighborhood<&'a T, U>, { - neighborhood.consider(&self.item); + let (node, children) = slice.split_first().unwrap(); + neighborhood.consider(&node.item); let target = neighborhood.target(); let ti = target.coordinate(i); - let si = self.item.coordinate(i); - let j = (i + 1) % self.item.dimensions(); + let ni = node.item.coordinate(i); + let j = (i + 1) % node.item.dimensions(); - let (near, far) = if ti <= si { - (&self.left, &self.right) + let (left, right) = children.split_at(node.left_len); + let (near, far) = if ti <= ni { + (left, right) } else { - (&self.right, &self.left) + (right, left) }; - if let Some(near) = near { - near.search(j, closest, neighborhood); + if !near.is_empty() { + Self::recurse(near, j, closest, neighborhood); } - if let Some(far) = far { + if !far.is_empty() { let saved = closest[i]; - closest[i] = si; + closest[i] = ni; if neighborhood.contains_distance(target.distance(closest)) { - far.search(j, closest, neighborhood); + Self::recurse(far, j, closest, neighborhood); } closest[i] = saved; } @@ -129,16 +139,14 @@ impl KdNode { /// A [k-d tree](https://en.wikipedia.org/wiki/K-d_tree). #[derive(Debug)] -pub struct KdTree { - root: Option>>, -} +pub struct KdTree(Vec>); impl FromIterator for KdTree { /// Create a new k-d tree from a set of points. fn from_iter>(items: I) -> Self { - Self { - root: KdNode::new(0, items.into_iter().collect()), - } + let mut nodes: Vec<_> = items.into_iter().map(KdNode::new).collect(); + KdNode::build(nodes.as_mut_slice(), 0); + Self(nodes) } } @@ -153,40 +161,27 @@ where U: 'b, N: Neighborhood<&'a T, &'b U>, { - let target = neighborhood.target(); - let dims = target.dimensions(); - let mut closest: Vec<_> = (0..dims).map(|i| target.coordinate(i)).collect(); + if !self.0.is_empty() { + let target = neighborhood.target(); + let dims = target.dimensions(); + let mut closest: Vec<_> = (0..dims).map(|i| target.coordinate(i)).collect(); - if let Some(root) = &self.root { - root.search(0, &mut closest, &mut neighborhood); + KdNode::recurse(&self.0, 0, &mut closest, &mut neighborhood); } + neighborhood } } /// An iterator that the moves values out of a k-d tree. #[derive(Debug)] -pub struct IntoIter { - stack: Vec>>, -} - -impl IntoIter { - fn new(node: Option>>) -> Self { - Self { - stack: node.into_iter().collect(), - } - } -} +pub struct IntoIter(std::vec::IntoIter>); impl Iterator for IntoIter { type Item = T; fn next(&mut self) -> Option { - self.stack.pop().map(|node| { - self.stack.extend(node.left); - self.stack.extend(node.right); - node.item - }) + self.0.next().map(|n| n.item) } } @@ -195,7 +190,7 @@ impl IntoIterator for KdTree { type IntoIter = IntoIter; fn into_iter(self) -> Self::IntoIter { - IntoIter::new(self.root) + IntoIter(self.0.into_iter()) } } -- cgit v1.2.3 From a4a75059f302de2a00971f1f485fcf4389710628 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sun, 19 Apr 2020 16:55:17 -0400 Subject: metric/forest: Implement dynamized forests --- src/metric.rs | 1 + src/metric/forest.rs | 152 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 153 insertions(+) create mode 100644 src/metric/forest.rs (limited to 'src') diff --git a/src/metric.rs b/src/metric.rs index 95191da..7a5f5f7 100644 --- a/src/metric.rs +++ b/src/metric.rs @@ -1,5 +1,6 @@ //! [Metric spaces](https://en.wikipedia.org/wiki/Metric_space). +pub mod forest; pub mod kd; pub mod vp; diff --git a/src/metric/forest.rs b/src/metric/forest.rs new file mode 100644 index 0000000..f23c451 --- /dev/null +++ b/src/metric/forest.rs @@ -0,0 +1,152 @@ +//! [Dynamization](https://en.wikipedia.org/wiki/Dynamization) for nearest neighbor search. + +use super::kd::KdTree; +use super::vp::VpTree; +use super::{Metric, NearestNeighbors, Neighborhood}; + +use std::iter::{Extend, Flatten, FromIterator}; + +/// 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(Vec>); + +impl Forest +where + U: FromIterator + IntoIterator, +{ + /// Create a new empty forest. + pub fn new() -> Self { + Self(Vec::new()) + } + + /// Add a new item to the forest. + pub fn push(&mut self, item: T) { + let mut items = vec![item]; + + for slot in &mut self.0 { + match slot.take() { + // Collect the items from any trees we encounter... + Some(tree) => { + items.extend(tree); + } + // ... and put them all in the first empty slot + None => { + *slot = Some(items.into_iter().collect()); + return; + } + } + } + + self.0.push(Some(items.into_iter().collect())); + } + + /// Get the number of items in the forest. + pub fn len(&self) -> usize { + let mut len = 0; + for (i, slot) in self.0.iter().enumerate() { + if slot.is_some() { + len |= 1 << i; + } + } + len + } +} + +impl Extend for Forest +where + U: FromIterator + IntoIterator, +{ + fn extend>(&mut self, items: I) { + for item in items { + self.push(item); + } + } +} + +impl FromIterator for Forest +where + U: FromIterator + IntoIterator, +{ + fn from_iter>(items: I) -> Self { + let mut forest = Self::new(); + forest.extend(items); + forest + } +} + +type IntoIterImpl = Flatten>>>; + +/// An iterator that moves items out of a forest. +pub struct IntoIter(IntoIterImpl); + +impl Iterator for IntoIter { + type Item = T::Item; + + fn next(&mut self) -> Option { + self.0.next() + } +} + +impl IntoIterator for Forest { + type Item = T::Item; + type IntoIter = IntoIter; + + fn into_iter(self) -> Self::IntoIter { + IntoIter(self.0.into_iter().flatten().flatten()) + } +} + +impl NearestNeighbors for Forest +where + U: Metric, + V: NearestNeighbors, +{ + fn search<'a, 'b, N>(&'a self, neighborhood: N) -> N + where + T: 'a, + U: 'b, + N: Neighborhood<&'a T, &'b U>, + { + self.0 + .iter() + .flatten() + .fold(neighborhood, |n, t| t.search(n)) + } +} + +/// A forest of k-d trees. +pub type KdForest = Forest>; + +/// A forest of vantage-point trees. +pub type VpForest = Forest>; + +#[cfg(test)] +mod tests { + use super::*; + + use crate::metric::tests::test_nearest_neighbors; + use crate::metric::ExhaustiveSearch; + + #[test] + fn test_exhaustive_forest() { + test_nearest_neighbors(Forest::>::from_iter); + } + + #[test] + fn test_forest_forest() { + test_nearest_neighbors(Forest::>>::from_iter); + } + + #[test] + fn test_kd_forest() { + test_nearest_neighbors(KdForest::from_iter); + } + + #[test] + fn test_vp_forest() { + test_nearest_neighbors(VpForest::from_iter); + } +} -- cgit v1.2.3 From b4a39a3f22fac361f6a535d281eee5586078281b Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Thu, 23 Apr 2020 14:58:47 -0400 Subject: metric/forest: Optimize bulk insertion --- src/metric/forest.rs | 47 +++++++++++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 20 deletions(-) (limited to 'src') diff --git a/src/metric/forest.rs b/src/metric/forest.rs index f23c451..29b6f55 100644 --- a/src/metric/forest.rs +++ b/src/metric/forest.rs @@ -4,7 +4,7 @@ use super::kd::KdTree; use super::vp::VpTree; use super::{Metric, NearestNeighbors, Neighborhood}; -use std::iter::{Extend, Flatten, FromIterator}; +use std::iter::{self, Extend, Flatten, FromIterator}; /// A dynamic wrapper for a static nearest neighbor search data structure. /// @@ -24,23 +24,7 @@ where /// Add a new item to the forest. pub fn push(&mut self, item: T) { - let mut items = vec![item]; - - for slot in &mut self.0 { - match slot.take() { - // Collect the items from any trees we encounter... - Some(tree) => { - items.extend(tree); - } - // ... and put them all in the first empty slot - None => { - *slot = Some(items.into_iter().collect()); - return; - } - } - } - - self.0.push(Some(items.into_iter().collect())); + self.extend(iter::once(item)); } /// Get the number of items in the forest. @@ -60,9 +44,32 @@ where U: FromIterator + IntoIterator, { fn extend>(&mut self, items: I) { - for item in items { - self.push(item); + let mut vec: Vec<_> = items.into_iter().collect(); + let new_len = self.len() + vec.len(); + + for i in 0.. { + let bit = 1 << i; + + if bit > new_len { + break; + } + + if i >= self.0.len() { + self.0.push(None); + } + + if new_len & bit == 0 { + if let Some(tree) = self.0[i].take() { + vec.extend(tree); + } + } else if self.0[i].is_none() { + let offset = vec.len() - bit; + self.0[i] = Some(vec.drain(offset..).collect()); + } } + + debug_assert!(vec.is_empty()); + debug_assert!(self.len() == new_len); } } -- cgit v1.2.3 From e53ace3e69ed6bacedb7de345df10d3e575a291e Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sun, 19 Apr 2020 16:58:52 -0400 Subject: metric/soft: Implement soft deletes --- src/metric.rs | 1 + src/metric/soft.rs | 282 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 283 insertions(+) create mode 100644 src/metric/soft.rs (limited to 'src') diff --git a/src/metric.rs b/src/metric.rs index 7a5f5f7..b46c8da 100644 --- a/src/metric.rs +++ b/src/metric.rs @@ -2,6 +2,7 @@ pub mod forest; pub mod kd; +pub mod soft; pub mod vp; use ordered_float::OrderedFloat; diff --git a/src/metric/soft.rs b/src/metric/soft.rs new file mode 100644 index 0000000..0d7dcdb --- /dev/null +++ b/src/metric/soft.rs @@ -0,0 +1,282 @@ +//! [Soft deletion](https://en.wiktionary.org/wiki/soft_deletion) for nearest neighbor search. + +use super::forest::{KdForest, VpForest}; +use super::kd::KdTree; +use super::vp::VpTree; +use super::{Metric, 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); + +impl Neighborhood for SoftNeighborhood +where + T: SoftDelete, + U: Metric, + N: Neighborhood, +{ + fn target(&self) -> U { + self.0.target() + } + + fn contains(&self, distance: f64) -> bool { + self.0.contains(distance) + } + + fn contains_distance(&self, distance: U::Distance) -> bool { + self.0.contains_distance(distance) + } + + fn consider(&mut self, item: T) -> U::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); + +impl SoftSearch +where + T: SoftDelete, + U: FromIterator + IntoIterator, +{ + /// 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, + { + 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> Extend for SoftSearch { + fn extend>(&mut self, iter: I) { + self.0.extend(iter); + } +} + +impl> FromIterator for SoftSearch { + fn from_iter>(iter: I) -> Self { + Self(U::from_iter(iter)) + } +} + +impl IntoIterator for SoftSearch { + type Item = T::Item; + type IntoIter = T::IntoIter; + + fn into_iter(self) -> Self::IntoIter { + self.0.into_iter() + } +} + +impl NearestNeighbors for SoftSearch +where + T: SoftDelete, + U: Metric, + V: NearestNeighbors, +{ + fn search<'a, 'b, N>(&'a self, neighborhood: N) -> N + where + T: 'a, + U: 'b, + N: Neighborhood<&'a T, &'b U>, + { + self.0.search(SoftNeighborhood(neighborhood)).0 + } +} + +/// A k-d forest that supports soft deletes. +pub type SoftKdForest = SoftSearch>; + +/// A k-d tree that supports soft deletes. +pub type SoftKdTree = SoftSearch>; + +/// A VP forest that supports soft deletes. +pub type SoftVpForest = SoftSearch>; + +/// A VP tree that supports soft deletes. +pub type SoftVpTree = SoftSearch>; + +#[cfg(test)] +mod tests { + use super::*; + + use crate::metric::kd::Cartesian; + use crate::metric::tests::Point; + use crate::metric::Neighbor; + + #[derive(Debug, PartialEq)] + struct SoftPoint { + point: Point, + deleted: bool, + } + + impl SoftPoint { + fn new(x: f64, y: f64, z: f64) -> Self { + Self { + point: Point([x, y, z]), + deleted: false, + } + } + + fn deleted(x: f64, y: f64, z: f64) -> Self { + Self { + point: Point([x, y, z]), + deleted: true, + } + } + } + + impl SoftDelete for SoftPoint { + fn is_deleted(&self) -> bool { + self.deleted + } + } + + impl Metric for SoftPoint { + type Distance = ::Distance; + + fn distance(&self, other: &Self) -> Self::Distance { + self.point.distance(&other.point) + } + } + + impl Metric<[f64]> for SoftPoint { + type Distance = ::Distance; + + fn distance(&self, other: &[f64]) -> Self::Distance { + self.point.distance(other) + } + } + + impl Cartesian for SoftPoint { + fn dimensions(&self) -> usize { + self.point.dimensions() + } + + fn coordinate(&self, i: usize) -> f64 { + self.point.coordinate(i) + } + } + + impl Metric for Point { + type Distance = ::Distance; + + fn distance(&self, other: &SoftPoint) -> Self::Distance { + self.distance(&other.point) + } + } + + fn test_index(index: &T) + where + T: NearestNeighbors, + { + let target = Point([0.0, 0.0, 0.0]); + + assert_eq!( + index.nearest(&target), + Some(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), + Some(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(index: &mut SoftSearch) + where + T: Extend, + T: FromIterator, + T: IntoIterator, + T: NearestNeighbors, + { + 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()); + } +} -- cgit v1.2.3 From 462daaffd5ec720ed80a2e7b1f445a73cabf5833 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sun, 19 Apr 2020 16:59:24 -0400 Subject: metric/approx: Implement approximate nearest neighbor search --- src/metric.rs | 1 + src/metric/approx.rs | 131 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 132 insertions(+) create mode 100644 src/metric/approx.rs (limited to 'src') diff --git a/src/metric.rs b/src/metric.rs index b46c8da..268aefd 100644 --- a/src/metric.rs +++ b/src/metric.rs @@ -1,5 +1,6 @@ //! [Metric spaces](https://en.wikipedia.org/wiki/Metric_space). +pub mod approx; pub mod forest; pub mod kd; pub mod soft; diff --git a/src/metric/approx.rs b/src/metric/approx.rs new file mode 100644 index 0000000..c23f9c7 --- /dev/null +++ b/src/metric/approx.rs @@ -0,0 +1,131 @@ +//! [Approximate nearest neighbor search](https://en.wikipedia.org/wiki/Nearest_neighbor_search#Approximate_nearest_neighbor). + +use super::{Metric, NearestNeighbors, Neighborhood}; + +/// An approximate [Neighborhood], for approximate nearest neighbor searches. +#[derive(Debug)] +struct ApproximateNeighborhood { + inner: N, + ratio: f64, + limit: usize, +} + +impl ApproximateNeighborhood { + fn new(inner: N, ratio: f64, limit: usize) -> Self { + Self { + inner, + ratio, + limit, + } + } +} + +impl Neighborhood for ApproximateNeighborhood +where + U: Metric, + N: Neighborhood, +{ + fn target(&self) -> U { + self.inner.target() + } + + fn contains(&self, distance: f64) -> bool { + if self.limit > 0 { + self.inner.contains(self.ratio * distance) + } else { + false + } + } + + fn contains_distance(&self, distance: U::Distance) -> bool { + self.contains(self.ratio * distance.into()) + } + + fn consider(&mut self, item: T) -> U::Distance { + self.limit = self.limit.saturating_sub(1); + self.inner.consider(item) + } +} + +/// An [approximate nearest neighbor search](https://en.wikipedia.org/wiki/Nearest_neighbor_search#Approximate_nearest_neighbor) +/// index. +/// +/// This wrapper converts an exact nearest neighbor search algorithm into an approximate one by +/// modifying the behavior of [Neighborhood::contains]. The approximation is controlled by two +/// parameters: +/// +/// * `ratio`: The [nearest neighbor distance ratio](https://en.wikipedia.org/wiki/Nearest_neighbor_search#Nearest_neighbor_distance_ratio), +/// which controls how much closer new candidates must be to be considered. For example, a ratio +/// of 2.0 means that a neighbor must be less than half of the current threshold away to be +/// considered. A ratio of 1.0 means an exact search. +/// +/// * `limit`: A limit on the number of candidates to consider. Typical nearest neighbor algorithms +/// find a close match quickly, so setting a limit bounds the worst-case search time while keeping +/// good accuracy. +#[derive(Debug)] +pub struct ApproximateSearch { + inner: T, + ratio: f64, + limit: usize, +} + +impl ApproximateSearch { + /// Create a new ApproximateSearch index. + /// + /// * `inner`: The [NearestNeighbors] implementation to wrap. + /// * `ratio`: The nearest neighbor distance ratio. + /// * `limit`: The maximum number of results to consider. + pub fn new(inner: T, ratio: f64, limit: usize) -> Self { + Self { + inner, + ratio, + limit, + } + } +} + +impl NearestNeighbors for ApproximateSearch +where + U: Metric, + V: NearestNeighbors, +{ + fn search<'a, 'b, N>(&'a self, neighborhood: N) -> N + where + T: 'a, + U: 'b, + N: Neighborhood<&'a T, &'b U>, + { + self.inner + .search(ApproximateNeighborhood::new( + neighborhood, + self.ratio, + self.limit, + )) + .inner + } +} + +#[cfg(test)] +mod tests { + use super::*; + + use crate::metric::kd::KdTree; + use crate::metric::tests::test_nearest_neighbors; + use crate::metric::vp::VpTree; + + use std::iter::FromIterator; + + #[test] + fn test_approx_kd_tree() { + test_nearest_neighbors(|iter| { + ApproximateSearch::new(KdTree::from_iter(iter), 1.0, std::usize::MAX) + }); + } + + #[test] + fn test_approx_vp_tree() { + test_nearest_neighbors(|iter| { + ApproximateSearch::new(VpTree::from_iter(iter), 1.0, std::usize::MAX) + }); + } +} -- cgit v1.2.3 From 232995aecf809309848b864a77e9d968c6185a29 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sat, 2 May 2020 13:43:16 -0400 Subject: color: Implement color spaces --- Cargo.lock | 296 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Cargo.toml | 1 + src/color.rs | 282 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/main.rs | 1 + 4 files changed, 580 insertions(+) create mode 100644 src/color.rs (limited to 'src') diff --git a/Cargo.lock b/Cargo.lock index 45602e1..7f81f15 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1,17 +1,119 @@ # This file is automatically @generated by Cargo. # It is not intended for manual editing. +[[package]] +name = "adler32" +version = "1.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5d2e7343e7fc9de883d1b0341e0b13970f764c14101234857d2ddafa1cb1cac2" + [[package]] name = "autocfg" version = "1.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f8aac770f1885fd7e387acedd76065302551364496e46b3dd00860b2f8359b9d" +[[package]] +name = "bitflags" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cf1de2fe8c75bc145a2f577add951f8134889b4795d47466a54a5c846d691693" + +[[package]] +name = "bytemuck" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "37fa13df2292ecb479ec23aa06f4507928bef07839be9ef15281411076629431" + +[[package]] +name = "byteorder" +version = "1.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08c48aae112d48ed9f069b33538ea9e3e90aa263cfa3d1c24309612b1f7472de" + [[package]] name = "cfg-if" version = "0.1.10" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "4785bdd1c96b2a846b2bd7cc02e86b6b3dbf14e7e53446c4f54c92a361040822" +[[package]] +name = "color_quant" +version = "1.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0dbbb57365263e881e805dc77d94697c9118fd94d8da011240555aa7b23445bd" + +[[package]] +name = "crc32fast" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ba125de2af0df55319f41944744ad91c71113bf74a4646efff39afe1f6842db1" +dependencies = [ + "cfg-if", +] + +[[package]] +name = "crossbeam-deque" +version = "0.7.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f02af974daeee82218205558e51ec8768b48cf524bd01d550abe5573a608285" +dependencies = [ + "crossbeam-epoch", + "crossbeam-utils", + "maybe-uninit", +] + +[[package]] +name = "crossbeam-epoch" +version = "0.8.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "058ed274caafc1f60c4997b5fc07bf7dc7cca454af7c6e81edffe5f33f70dace" +dependencies = [ + "autocfg", + "cfg-if", + "crossbeam-utils", + "lazy_static", + "maybe-uninit", + "memoffset", + "scopeguard", +] + +[[package]] +name = "crossbeam-queue" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c695eeca1e7173472a32221542ae469b3e9aac3a4fc81f7696bcad82029493db" +dependencies = [ + "cfg-if", + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-utils" +version = "0.7.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c3c7c73a2d1e9fc0886a08b93e98eb643461230d5f1925e4036204d5f2e261a8" +dependencies = [ + "autocfg", + "cfg-if", + "lazy_static", +] + +[[package]] +name = "deflate" +version = "0.8.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e7e5d2a2273fed52a7f947ee55b092c4057025d7a3e04e5ecdbd25d6c3fb1bd7" +dependencies = [ + "adler32", + "byteorder", +] + +[[package]] +name = "either" +version = "1.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bb1f6b1ce1c140482ea30ddd3335fc0024ac7ee112895426e0a629a6c20adfe3" + [[package]] name = "getrandom" version = "0.1.14" @@ -23,20 +125,145 @@ dependencies = [ "wasi", ] +[[package]] +name = "gif" +version = "0.10.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "471d90201b3b223f3451cd4ad53e34295f16a1df17b1edf3736d47761c3981af" +dependencies = [ + "color_quant", + "lzw", +] + +[[package]] +name = "hermit-abi" +version = "0.1.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "61565ff7aaace3525556587bd2dc31d4a07071957be715e63ce7b1eccf51a8f4" +dependencies = [ + "libc", +] + +[[package]] +name = "image" +version = "0.23.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9117f4167a8f21fa2bb3f17a652a760acd7572645281c98e3b612a26242c96ee" +dependencies = [ + "bytemuck", + "byteorder", + "gif", + "jpeg-decoder", + "num-iter", + "num-rational", + "num-traits", + "png", + "scoped_threadpool", + "tiff", +] + +[[package]] +name = "inflate" +version = "0.4.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1cdb29978cc5797bd8dcc8e5bf7de604891df2a8dc576973d71a281e916db2ff" +dependencies = [ + "adler32", +] + +[[package]] +name = "jpeg-decoder" +version = "0.1.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5b47b4c4e017b01abdc5bcc126d2d1002e5a75bbe3ce73f9f4f311a916363704" +dependencies = [ + "byteorder", + "rayon", +] + [[package]] name = "kd-forest" version = "2.0.0" dependencies = [ + "image", "ordered-float", "rand", ] +[[package]] +name = "lazy_static" +version = "1.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" + [[package]] name = "libc" version = "0.2.69" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "99e85c08494b21a9054e7fe1374a732aeadaff3980b6990b94bfd3a70f690005" +[[package]] +name = "lzw" +version = "0.10.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7d947cbb889ed21c2a84be6ffbaebf5b4e0f4340638cba0444907e38b56be084" + +[[package]] +name = "maybe-uninit" +version = "2.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60302e4db3a61da70c0cb7991976248362f30319e88850c487b9b95bbf059e00" + +[[package]] +name = "memoffset" +version = "0.5.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b4fc2c02a7e374099d4ee95a193111f72d2110197fe200272371758f6c3643d8" +dependencies = [ + "autocfg", +] + +[[package]] +name = "miniz_oxide" +version = "0.3.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa679ff6578b1cddee93d7e82e263b94a575e0bfced07284eb0c037c1d2416a5" +dependencies = [ + "adler32", +] + +[[package]] +name = "num-integer" +version = "0.1.42" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f6ea62e9d81a77cd3ee9a2a5b9b609447857f3d358704331e4ef39eb247fcba" +dependencies = [ + "autocfg", + "num-traits", +] + +[[package]] +name = "num-iter" +version = "0.1.40" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dfb0800a0291891dd9f4fe7bd9c19384f98f7fbe0cd0f39a2c6b88b9868bbc00" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5c000134b5dbf44adc5cb772486d335293351644b801551abe8f75c84cfa4aef" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + [[package]] name = "num-traits" version = "0.2.11" @@ -46,6 +273,16 @@ dependencies = [ "autocfg", ] +[[package]] +name = "num_cpus" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "05499f3756671c15885fee9034446956fff3f243d6077b91e5767df161f766b3" +dependencies = [ + "hermit-abi", + "libc", +] + [[package]] name = "ordered-float" version = "1.0.2" @@ -55,6 +292,18 @@ dependencies = [ "num-traits", ] +[[package]] +name = "png" +version = "0.16.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2c68a431ed29933a4eb5709aca9800989758c97759345860fa5db3cfced0b65d" +dependencies = [ + "bitflags", + "crc32fast", + "deflate", + "inflate", +] + [[package]] name = "ppv-lite86" version = "0.2.6" @@ -102,6 +351,53 @@ dependencies = [ "rand_core", ] +[[package]] +name = "rayon" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "db6ce3297f9c85e16621bb8cca38a06779ffc31bb8184e1be4bed2be4678a098" +dependencies = [ + "crossbeam-deque", + "either", + "rayon-core", +] + +[[package]] +name = "rayon-core" +version = "1.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08a89b46efaf957e52b18062fb2f4660f8b8a4dde1807ca002690868ef2c85a9" +dependencies = [ + "crossbeam-deque", + "crossbeam-queue", + "crossbeam-utils", + "lazy_static", + "num_cpus", +] + +[[package]] +name = "scoped_threadpool" +version = "0.1.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1d51f5df5af43ab3f1360b429fa5e0152ac5ce8c0bd6485cae490332e96846a8" + +[[package]] +name = "scopeguard" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd" + +[[package]] +name = "tiff" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "002351e428db1eb1d8656d4ca61947c3519ac3191e1c804d4600cd32093b77ad" +dependencies = [ + "byteorder", + "lzw", + "miniz_oxide", +] + [[package]] name = "wasi" version = "0.9.0+wasi-snapshot-preview1" diff --git a/Cargo.toml b/Cargo.toml index 5a93cf5..902c3a3 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -5,5 +5,6 @@ authors = ["Tavian Barnes "] edition = "2018" [dependencies] +image = "0.23.4" ordered-float = "1.0.2" rand = "0.7.3" diff --git a/src/color.rs b/src/color.rs new file mode 100644 index 0000000..a062795 --- /dev/null +++ b/src/color.rs @@ -0,0 +1,282 @@ +//! Colors and color spaces. + +use crate::metric::kd::{Cartesian, CartesianMetric}; +use crate::metric::{Metric, SquaredDistance}; + +use image::Rgb; + +use std::ops::Index; + +/// An 8-bit RGB color. +pub type Rgb8 = Rgb; + +/// A [color space](https://en.wikipedia.org/wiki/Color_space). +pub trait ColorSpace: Copy + From + CartesianMetric { + /// Compute the average of the given colors. + fn average>(colors: I) -> Self; +} + +/// [sRGB](https://en.wikipedia.org/wiki/SRGB) space. +#[derive(Clone, Copy, Debug)] +pub struct RgbSpace([f64; 3]); + +impl Index for RgbSpace { + type Output = f64; + + fn index(&self, i: usize) -> &f64 { + &self.0[i] + } +} + +impl From 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 Metric<[f64]> for RgbSpace { + type Distance = SquaredDistance; + + fn distance(&self, other: &[f64]) -> Self::Distance { + self.0.distance(other) + } +} + +impl Metric for RgbSpace { + type Distance = SquaredDistance; + + fn distance(&self, other: &Self) -> Self::Distance { + self.0.distance(&other.0) + } +} + +impl Cartesian for RgbSpace { + fn dimensions(&self) -> usize { + self.0.dimensions() + } + + fn coordinate(&self, i: usize) -> f64 { + self.0.coordinate(i) + } +} + +impl ColorSpace for RgbSpace { + fn average>(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 i in 0..3 { + sum[i] /= 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 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 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 for LabSpace { + type Output = f64; + + fn index(&self, i: usize) -> &f64 { + &self.0[i] + } +} + +impl From 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 Metric<[f64]> for LabSpace { + type Distance = SquaredDistance; + + fn distance(&self, other: &[f64]) -> Self::Distance { + self.0.distance(other) + } +} + +impl Metric for LabSpace { + type Distance = SquaredDistance; + + fn distance(&self, other: &Self) -> Self::Distance { + self.0.distance(&other.0) + } +} + +impl Cartesian for LabSpace { + fn dimensions(&self) -> usize { + self.0.dimensions() + } + + fn coordinate(&self, i: usize) -> f64 { + self.0.coordinate(i) + } +} + +impl ColorSpace for LabSpace { + fn average>(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 i in 0..3 { + sum[i] /= 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 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 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 Metric<[f64]> for LuvSpace { + type Distance = SquaredDistance; + + fn distance(&self, other: &[f64]) -> Self::Distance { + self.0.distance(other) + } +} + +impl Metric for LuvSpace { + type Distance = SquaredDistance; + + fn distance(&self, other: &Self) -> Self::Distance { + self.0.distance(&other.0) + } +} + +impl Cartesian for LuvSpace { + fn dimensions(&self) -> usize { + self.0.dimensions() + } + + fn coordinate(&self, i: usize) -> f64 { + self.0.coordinate(i) + } +} + +impl ColorSpace for LuvSpace { + fn average>(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 i in 0..3 { + sum[i] /= len as f64; + } + Self(sum) + } +} diff --git a/src/main.rs b/src/main.rs index 0d7989b..a7bda67 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,3 +1,4 @@ +pub mod color; pub mod metric; fn main() {} -- cgit v1.2.3 From 5aea448bd2a7f6157cf33a6a2c044d043bbcd3ec Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sat, 2 May 2020 13:46:25 -0400 Subject: color/source: Implement color sources --- src/color.rs | 2 ++ src/color/source.rs | 76 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 78 insertions(+) create mode 100644 src/color/source.rs (limited to 'src') diff --git a/src/color.rs b/src/color.rs index a062795..e0a3399 100644 --- a/src/color.rs +++ b/src/color.rs @@ -1,5 +1,7 @@ //! Colors and color spaces. +pub mod source; + use crate::metric::kd::{Cartesian, CartesianMetric}; use crate::metric::{Metric, SquaredDistance}; diff --git a/src/color/source.rs b/src/color/source.rs new file mode 100644 index 0000000..bd752d9 --- /dev/null +++ b/src/color/source.rs @@ -0,0 +1,76 @@ +//! 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: [usize; 3], +} + +impl AllColors { + /// Create an AllColors source with the given bit depth. + pub fn new(bits: usize) -> Self { + // Allocate bits from most to least perceptually important + let gbits = (bits + 2) / 3; + let rbits = (bits + 1) / 3; + let bbits = bits / 3; + + Self { + dims: [1 << rbits, 1 << gbits, 1 << bbits], + shifts: [8 - rbits, 8 - gbits, 8 - bbits], + } + } +} + +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 for ImageColors { + fn from(image: RgbImage) -> Self { + Self { + dims: [image.width() as usize, image.height() as usize], + image: 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) + } +} -- cgit v1.2.3 From 62e0fec044b5727efa1841138f44d9a1d9537bcf Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sat, 2 May 2020 13:48:07 -0400 Subject: color/order: Implement color orderings --- src/color.rs | 1 + src/color/order.rs | 196 +++++++++++++++++++++++++++++++++++++++++++++++++++++ src/hilbert.rs | 136 +++++++++++++++++++++++++++++++++++++ src/main.rs | 1 + 4 files changed, 334 insertions(+) create mode 100644 src/color/order.rs create mode 100644 src/hilbert.rs (limited to 'src') diff --git a/src/color.rs b/src/color.rs index e0a3399..64fd82b 100644 --- a/src/color.rs +++ b/src/color.rs @@ -1,5 +1,6 @@ //! Colors and color spaces. +pub mod order; pub mod source; use crate::metric::kd::{Cartesian, CartesianMetric}; 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 { + source: S, + coords: Vec, +} + +impl From for ColorSourceIter { + fn from(source: S) -> Self { + let coords = vec![0; source.dimensions().len()]; + + Self { source, coords } + } +} + +impl Iterator for ColorSourceIter { + type Item = Rgb8; + + fn next(&mut self) -> Option { + 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 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 { + Some(self.cmp(other)) + } +} + +/// Iterate over colors sorted by their hue. +pub fn hue_sorted(source: S) -> Vec { + 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(source: S, rng: &mut R) -> Vec { + 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::() as u32; + nbits - (n - 1).leading_zeros() +} + +/// Iterate over colors in Morton order (Z-order). +pub fn morton(source: S) -> Vec { + 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(source: S) -> Vec { + 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) -> Vec { + 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/hilbert.rs b/src/hilbert.rs new file mode 100644 index 0000000..c0982d4 --- /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 = 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 index a7bda67..a59a0cf 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,4 +1,5 @@ pub mod color; +pub mod hilbert; pub mod metric; fn main() {} -- cgit v1.2.3 From 5ac571d8d16a307cea2922587185557bc773e8ed Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sat, 2 May 2020 13:52:28 -0400 Subject: frontier: New trait for choosing color locations --- src/frontier.rs | 18 ++++++++++++++++++ src/main.rs | 1 + 2 files changed, 19 insertions(+) create mode 100644 src/frontier.rs (limited to 'src') diff --git a/src/frontier.rs b/src/frontier.rs new file mode 100644 index 0000000..2c6f43a --- /dev/null +++ b/src/frontier.rs @@ -0,0 +1,18 @@ +//! Frontiers on which to place pixels. + +use crate::color::Rgb8; + +/// 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; + + /// Place the given color on the frontier, and return its position. + fn place(&mut self, rgb8: Rgb8) -> Option<(u32, u32)>; +} diff --git a/src/main.rs b/src/main.rs index a59a0cf..07e138a 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,4 +1,5 @@ pub mod color; +pub mod frontier; pub mod hilbert; pub mod metric; -- cgit v1.2.3 From da9b2ad1310e1db0ccb26a53181fa7f9b9033d04 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sat, 2 May 2020 13:56:09 -0400 Subject: frontier/image: Implement image frontiers --- src/frontier.rs | 71 +++++++++++++++++++++++++++++++++++++++++++++++- src/frontier/image.rs | 74 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 144 insertions(+), 1 deletion(-) create mode 100644 src/frontier/image.rs (limited to 'src') diff --git a/src/frontier.rs b/src/frontier.rs index 2c6f43a..1c5a173 100644 --- a/src/frontier.rs +++ b/src/frontier.rs @@ -1,6 +1,13 @@ //! Frontiers on which to place pixels. -use crate::color::Rgb8; +pub mod image; + +use crate::color::{ColorSpace, Rgb8}; +use crate::metric::kd::Cartesian; +use crate::metric::soft::SoftDelete; +use crate::metric::Metric; + +use std::cell::Cell; /// A frontier of pixels. pub trait Frontier { @@ -16,3 +23,65 @@ pub trait Frontier { /// 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 { + pos: (u32, u32), + color: C, + deleted: Cell, +} + +impl Pixel { + 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); + } +} + +impl Metric> for C { + type Distance = C::Distance; + + fn distance(&self, other: &Pixel) -> Self::Distance { + self.distance(&other.color) + } +} + +impl> Metric<[f64]> for Pixel { + type Distance = C::Distance; + + fn distance(&self, other: &[f64]) -> Self::Distance { + self.color.distance(other) + } +} + +impl Metric for Pixel { + type Distance = C::Distance; + + fn distance(&self, other: &Pixel) -> Self::Distance { + self.color.distance(&other.color) + } +} + +impl Cartesian for Pixel { + fn dimensions(&self) -> usize { + self.color.dimensions() + } + + fn coordinate(&self, i: usize) -> f64 { + self.color.coordinate(i) + } +} + +impl SoftDelete for Pixel { + fn is_deleted(&self) -> bool { + self.deleted.get() + } +} diff --git a/src/frontier/image.rs b/src/frontier/image.rs new file mode 100644 index 0000000..3655580 --- /dev/null +++ b/src/frontier/image.rs @@ -0,0 +1,74 @@ +//! Frontier that targets an image. + +use super::{Frontier, Pixel}; + +use crate::color::{ColorSpace, Rgb8}; +use crate::metric::soft::SoftKdTree; +use crate::metric::NearestNeighbors; + +use image::RgbImage; + +/// A [Frontier] that places colors on the closest pixel of a target image. +#[derive(Debug)] +pub struct ImageFrontier { + nodes: SoftKdTree>, + width: u32, + height: u32, + len: usize, + deleted: usize, +} + +impl ImageFrontier { + /// 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 Frontier for ImageFrontier { + fn width(&self) -> u32 { + self.width + } + + fn height(&self) -> u32 { + self.height + } + + fn len(&self) -> usize { + self.len - self.deleted + } + + fn place(&mut self, rgb8: Rgb8) -> Option<(u32, u32)> { + let color = C::from(rgb8); + + if let Some(node) = self.nodes.nearest(&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 + } + } +} -- cgit v1.2.3 From dd2336f76bfaff01ccb5b57f4453629ff62016b3 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sat, 2 May 2020 13:57:30 -0400 Subject: frontier/min: Implement min selection --- src/frontier.rs | 61 +++++++++++++++++++++ src/frontier/min.rs | 150 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 211 insertions(+) create mode 100644 src/frontier/min.rs (limited to 'src') diff --git a/src/frontier.rs b/src/frontier.rs index 1c5a173..e1b6528 100644 --- a/src/frontier.rs +++ b/src/frontier.rs @@ -1,6 +1,7 @@ //! Frontiers on which to place pixels. pub mod image; +pub mod min; use crate::color::{ColorSpace, Rgb8}; use crate::metric::kd::Cartesian; @@ -8,6 +9,7 @@ use crate::metric::soft::SoftDelete; use crate::metric::Metric; use std::cell::Cell; +use std::rc::Rc; /// A frontier of pixels. pub trait Frontier { @@ -85,3 +87,62 @@ impl SoftDelete for Pixel { self.deleted.get() } } + +impl> Metric<[f64]> for Rc> { + type Distance = C::Distance; + + fn distance(&self, other: &[f64]) -> Self::Distance { + (**self).distance(other) + } +} + +impl Metric>> for C { + type Distance = C::Distance; + + fn distance(&self, other: &Rc>) -> Self::Distance { + self.distance(&other.color) + } +} + +impl Metric for Rc> { + type Distance = C::Distance; + + fn distance(&self, other: &Self) -> Self::Distance { + (**self).distance(&**other) + } +} + +impl Cartesian for Rc> { + fn dimensions(&self) -> usize { + (**self).dimensions() + } + + fn coordinate(&self, i: usize) -> f64 { + (**self).coordinate(i) + } +} + +impl SoftDelete for Rc> { + fn is_deleted(&self) -> bool { + (**self).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/min.rs b/src/frontier/min.rs new file mode 100644 index 0000000..b22b290 --- /dev/null +++ b/src/frontier/min.rs @@ -0,0 +1,150 @@ +//! Minimum selection frontier. + +use super::{neighbors, Frontier, Pixel}; + +use crate::color::{ColorSpace, Rgb8}; +use crate::metric::soft::SoftKdForest; +use crate::metric::NearestNeighbors; + +use rand::Rng; + +use std::rc::Rc; + +/// A pixel on a min frontier. +#[derive(Debug)] +struct MinPixel { + pixel: Option>>, + filled: bool, +} + +impl MinPixel { + 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 { + rng: R, + pixels: Vec>, + forest: SoftKdForest>>, + width: u32, + height: u32, + x0: u32, + y0: u32, + len: usize, + deleted: usize, +} + +impl MinFrontier { + /// 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 = Rc::new(Pixel::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 Frontier for MinFrontier { + fn width(&self) -> u32 { + self.width + } + + fn height(&self) -> u32 { + self.height + } + + fn len(&self) -> usize { + self.len - self.deleted + } + + fn place(&mut self, rgb8: Rgb8) -> Option<(u32, u32)> { + let color = C::from(rgb8); + let (x, y) = self + .forest + .nearest(&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) + } +} -- cgit v1.2.3 From b58924b5a2534f845b996a2d58ab618301f1d450 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sat, 2 May 2020 13:58:14 -0400 Subject: frontier/mean: Implement mean selection --- src/frontier.rs | 1 + src/frontier/mean.rs | 140 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 141 insertions(+) create mode 100644 src/frontier/mean.rs (limited to 'src') diff --git a/src/frontier.rs b/src/frontier.rs index e1b6528..7143cb7 100644 --- a/src/frontier.rs +++ b/src/frontier.rs @@ -1,6 +1,7 @@ //! Frontiers on which to place pixels. pub mod image; +pub mod mean; pub mod min; use crate::color::{ColorSpace, Rgb8}; diff --git a/src/frontier/mean.rs b/src/frontier/mean.rs new file mode 100644 index 0000000..889c5ba --- /dev/null +++ b/src/frontier/mean.rs @@ -0,0 +1,140 @@ +//! Mean selection frontier. + +use super::{neighbors, Frontier, Pixel}; + +use crate::color::{ColorSpace, Rgb8}; +use crate::metric::soft::SoftKdForest; +use crate::metric::NearestNeighbors; + +use std::iter; +use std::rc::Rc; + +/// A pixel on a mean frontier. +#[derive(Debug)] +enum MeanPixel { + Empty, + Fillable(Rc>), + Filled(C), +} + +impl MeanPixel { + fn filled_color(&self) -> Option { + match self { + Self::Filled(color) => Some(*color), + _ => None, + } + } +} + +/// A [Frontier] that looks at the average color of each pixel's neighbors. +pub struct MeanFrontier { + pixels: Vec>, + forest: SoftKdForest>>, + width: u32, + height: u32, + len: usize, + deleted: usize, +} + +impl MeanFrontier { + /// 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 = Rc::new(Pixel::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 = Rc::new(Pixel::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 Frontier for MeanFrontier { + fn width(&self) -> u32 { + self.width + } + + fn height(&self) -> u32 { + self.height + } + + fn len(&self) -> usize { + self.len - self.deleted + } + + fn place(&mut self, rgb8: Rgb8) -> Option<(u32, u32)> { + let color = C::from(rgb8); + let (x, y) = self.forest.nearest(&color).map(|n| n.item.pos)?; + + self.fill(x, y, color); + + Some((x, y)) + } +} -- cgit v1.2.3 From 2984e8f93fe88d0ee7eb3c0561dcd2da44807429 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Sat, 2 May 2020 13:58:45 -0400 Subject: main: Implement the main binary --- Cargo.lock | 192 +++++++++++++++++++++++++++++ Cargo.toml | 3 + src/main.rs | 397 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 591 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/Cargo.lock b/Cargo.lock index 7f81f15..c066dd8 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -6,18 +6,67 @@ version = "1.0.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5d2e7343e7fc9de883d1b0341e0b13970f764c14101234857d2ddafa1cb1cac2" +[[package]] +name = "ansi_term" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ee49baf6cb617b853aa8d93bf420db2383fab46d314482ca2803b40d5fde979b" +dependencies = [ + "winapi", +] + +[[package]] +name = "arrayref" +version = "0.3.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a4c527152e37cf757a3f78aae5a06fbeefdb07ccc535c980a3208ee3060dd544" + +[[package]] +name = "arrayvec" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cff77d8686867eceff3105329d4698d96c2391c176d5d03adc90c7389162b5b8" + +[[package]] +name = "atty" +version = "0.2.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d9b39be18770d11421cdb1b9947a45dd3f37e93092cbf377614828a319d5fee8" +dependencies = [ + "hermit-abi", + "libc", + "winapi", +] + [[package]] name = "autocfg" version = "1.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f8aac770f1885fd7e387acedd76065302551364496e46b3dd00860b2f8359b9d" +[[package]] +name = "base64" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b41b7ea54a0c9d92199de89e20e58d49f02f8e699814ef3fdf266f6f748d15c7" + [[package]] name = "bitflags" version = "1.2.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "cf1de2fe8c75bc145a2f577add951f8134889b4795d47466a54a5c846d691693" +[[package]] +name = "blake2b_simd" +version = "0.5.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d8fb2d74254a3a0b5cac33ac9f8ed0e44aa50378d9dbb2e5d83bd21ed1dc2c8a" +dependencies = [ + "arrayref", + "arrayvec", + "constant_time_eq", +] + [[package]] name = "bytemuck" version = "1.2.0" @@ -36,12 +85,33 @@ version = "0.1.10" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "4785bdd1c96b2a846b2bd7cc02e86b6b3dbf14e7e53446c4f54c92a361040822" +[[package]] +name = "clap" +version = "2.33.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5067f5bb2d80ef5d68b4c87db81601f0b75bca627bc2ef76b141d7b846a3c6d9" +dependencies = [ + "ansi_term", + "atty", + "bitflags", + "strsim", + "textwrap", + "unicode-width", + "vec_map", +] + [[package]] name = "color_quant" version = "1.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0dbbb57365263e881e805dc77d94697c9118fd94d8da011240555aa7b23445bd" +[[package]] +name = "constant_time_eq" +version = "0.1.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "245097e9a4535ee1e3e3931fcfcd55a796a44c643e8596ff6566d68f09b87bbc" + [[package]] name = "crc32fast" version = "1.2.0" @@ -108,6 +178,28 @@ dependencies = [ "byteorder", ] +[[package]] +name = "dirs" +version = "2.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "13aea89a5c93364a98e9b37b2fa237effbb694d5cfe01c5b70941f7eb087d5e3" +dependencies = [ + "cfg-if", + "dirs-sys", +] + +[[package]] +name = "dirs-sys" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "afa0b23de8fd801745c471deffa6e12d248f962c9fd4b4c33787b055599bde7b" +dependencies = [ + "cfg-if", + "libc", + "redox_users", + "winapi", +] + [[package]] name = "either" version = "1.5.3" @@ -185,9 +277,12 @@ dependencies = [ name = "kd-forest" version = "2.0.0" dependencies = [ + "clap", "image", "ordered-float", "rand", + "rand_pcg", + "term", ] [[package]] @@ -351,6 +446,15 @@ dependencies = [ "rand_core", ] +[[package]] +name = "rand_pcg" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "16abd0c1b639e9eb4d7c50c0b8100b0d0f849be2349829c740fe8e6eb4816429" +dependencies = [ + "rand_core", +] + [[package]] name = "rayon" version = "1.3.0" @@ -375,6 +479,35 @@ dependencies = [ "num_cpus", ] +[[package]] +name = "redox_syscall" +version = "0.1.56" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2439c63f3f6139d1b57529d16bc3b8bb855230c8efcc5d3a896c8bea7c3b1e84" + +[[package]] +name = "redox_users" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "09b23093265f8d200fa7b4c2c76297f47e681c655f6f1285a8780d6a022f7431" +dependencies = [ + "getrandom", + "redox_syscall", + "rust-argon2", +] + +[[package]] +name = "rust-argon2" +version = "0.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2bc8af4bda8e1ff4932523b94d3dd20ee30a87232323eda55903ffd71d2fb017" +dependencies = [ + "base64", + "blake2b_simd", + "constant_time_eq", + "crossbeam-utils", +] + [[package]] name = "scoped_threadpool" version = "0.1.9" @@ -387,6 +520,31 @@ version = "1.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd" +[[package]] +name = "strsim" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8ea5119cdb4c55b55d432abb513a0429384878c15dde60cc77b1c99de1a95a6a" + +[[package]] +name = "term" +version = "0.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c0863a3345e70f61d613eab32ee046ccd1bcc5f9105fe402c61fcd0c13eeb8b5" +dependencies = [ + "dirs", + "winapi", +] + +[[package]] +name = "textwrap" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d326610f408c7a4eb6f51c37c330e496b08506c9457c9d34287ecc38809fb060" +dependencies = [ + "unicode-width", +] + [[package]] name = "tiff" version = "0.4.0" @@ -398,8 +556,42 @@ dependencies = [ "miniz_oxide", ] +[[package]] +name = "unicode-width" +version = "0.1.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "caaa9d531767d1ff2150b9332433f32a24622147e5ebb1f26409d5da67afd479" + +[[package]] +name = "vec_map" +version = "0.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "05c78687fb1a80548ae3250346c3db86a80a7cdd77bda190189f2d0a0987c81a" + [[package]] name = "wasi" version = "0.9.0+wasi-snapshot-preview1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "cccddf32554fecc6acb585f82a32a72e28b48f8c4c1883ddfeeeaa96f7d8e519" + +[[package]] +name = "winapi" +version = "0.3.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8093091eeb260906a183e6ae1abdba2ef5ef2257a21801128899c3fc699229c6" +dependencies = [ + "winapi-i686-pc-windows-gnu", + "winapi-x86_64-pc-windows-gnu", +] + +[[package]] +name = "winapi-i686-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ac3b87c63620426dd9b991e5ce0329eff545bccbbb34f3be09ff6fb6ab51b7b6" + +[[package]] +name = "winapi-x86_64-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" diff --git a/Cargo.toml b/Cargo.toml index 902c3a3..de65bb7 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -5,6 +5,9 @@ authors = ["Tavian Barnes "] edition = "2018" [dependencies] +clap = "2.33.0" image = "0.23.4" ordered-float = "1.0.2" rand = "0.7.3" +rand_pcg = "0.2.1" +term = "0.6.1" diff --git a/src/main.rs b/src/main.rs index 07e138a..f016b4c 100644 --- a/src/main.rs +++ b/src/main.rs @@ -3,4 +3,399 @@ pub mod frontier; pub mod hilbert; pub mod metric; -fn main() {} +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, AppSettings}; + +use image::{self, Rgba, RgbaImage}; + +use rand::SeedableRng; +use rand_pcg::Pcg64; + +use term; + +use std::cmp; +use std::error::Error; +use std::fs; +use std::io::{self, Write}; +use std::path::PathBuf; +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. + AllRgb(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, +} + +/// Parse an argument into the appropriate type. +fn parse_arg(arg: Option<&str>) -> clap::Result> +where + F: FromStr, + F::Err: Error, +{ + match arg.map(|s| s.parse()) { + Some(Ok(f)) => Ok(Some(f)), + Some(Err(e)) => Err(clap::Error::with_description( + &e.to_string(), + clap::ErrorKind::InvalidValue, + )), + None => Ok(None), + } +} + +/// The parsed command line arguments. +#[derive(Debug)] +struct Args { + source: SourceArg, + order: OrderArg, + stripe: bool, + frontier: FrontierArg, + space: ColorSpaceArg, + width: Option, + height: Option, + x0: Option, + y0: Option, + animate: bool, + output: PathBuf, + seed: u64, +} + +impl Args { + fn parse() -> clap::Result { + let args = clap_app!((crate_name!()) => + (version: crate_version!()) + (author: crate_authors!()) + (setting: AppSettings::ColoredHelp) + (@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 { + SourceArg::AllRgb(parse_arg(args.value_of("DEPTH"))?.unwrap_or(24)) + }; + + 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 mut path = args.value_of("PATH").unwrap(); + if animate && args.occurrences_of("PATH") == 0 { + path = "kd-frames"; + } + 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, + }) + } +} + +/// main() return type. +type MainResult = Result<(), Box>; + +/// The kd-forest application itself. +#[derive(Debug)] +struct App { + args: Args, + rng: Pcg64, + width: Option, + height: Option, + 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) -> MainResult { + let colors = match self.args.source { + SourceArg::AllRgb(depth) => { + self.width.get_or_insert(1u32 << ((depth + 1) / 2)); + self.height.get_or_insert(1u32 << (depth / 2)); + self.get_colors(AllColors::new(depth as usize)) + } + 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::(colors), + ColorSpaceArg::Lab => self.paint::(colors), + ColorSpaceArg::Luv => self.paint::(colors), + } + } + + fn get_colors(&mut self, source: S) -> Vec { + 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(&mut self, colors: Vec) -> MainResult { + let width = self.width.unwrap(); + let height = self.height.unwrap(); + let x0 = self.args.x0.unwrap_or(width / 2); + let y0 = self.args.x0.unwrap_or(height / 2); + + match &self.args.frontier { + FrontierArg::Image(ref path) => { + let img = image::open(path)?.into_rgb(); + self.paint_on(colors, ImageFrontier::::new(&img)) + } + FrontierArg::Min => { + let rng = Pcg64::from_rng(&mut self.rng)?; + self.paint_on(colors, MinFrontier::::new(rng, width, height, x0, y0)) + } + FrontierArg::Mean => { + self.paint_on(colors, MeanFrontier::::new(width, height, x0, y0)) + } + } + } + + fn paint_on(&mut self, colors: Vec, mut frontier: F) -> MainResult { + 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() -> MainResult { + let args = match Args::parse() { + Ok(args) => args, + Err(e) => e.exit(), + }; + + App::new(args).run() +} -- cgit v1.2.3