space_filling/solver/adf/
mod.rs

1//! Adaptive Distance Field, uses quadtree as an underlying data structire.
2//! Each node (bucket) stores several `Arc<dyn Fn(Point2D) -> {float}>`
3
4#![allow(clippy::mut_from_ref)]
5use {
6  crate::{
7    solver::LineSearch,
8    geometry::{Shape, shapes, P2, WorldSpace, BoundingBox},
9    sdf::SDF,
10  },
11  quadtree::{
12    Quadtree, TraverseCommand
13  },
14  std::{
15    sync::{
16      Arc, atomic::{AtomicBool, Ordering}
17    },
18    fmt::{Debug, Formatter}
19  },
20  euclid::{Point2D, Box2D, Rect},
21  num_traits::{Float, Signed}
22};
23
24#[cfg(test)] mod tests;
25pub(crate) mod quadtree;
26
27#[derive(Clone)]
28pub struct ADF<Float> {
29  pub tree: Quadtree<Vec<Arc<dyn Fn(P2<Float>) -> Float>>, Float>,
30  /// Gradient Descent lattice density, N^2
31  /// higher values improve precision
32  ipm_gd_lattice_density: u32,
33  ipm_line_config: LineSearch<Float>
34}
35
36unsafe impl<Float> Send for ADF<Float> {}
37unsafe impl<Float> Sync for ADF<Float> {}
38
39impl <_Float: Float> SDF<_Float> for &[Arc<dyn Fn(P2<_Float>) -> _Float>] {
40  fn sdf(&self, pixel: P2<_Float>) -> _Float {
41    self.iter()
42      .map(|f| f(pixel))
43      .reduce(|a, b| if a <= b { a } else { b })
44      .unwrap_or(_Float::max_value() / (_Float::one() + _Float::one()))
45  }
46}
47
48fn sdf_partialord<_Float: Float + Signed>(
49  f: impl Fn(P2<_Float>) -> _Float,
50  g: impl Fn(P2<_Float>) -> _Float,
51  domain: Rect<_Float, WorldSpace>,
52  lattice_density: u32,
53  line_search: LineSearch<_Float>
54) -> bool {
55  let boundary_constraint = |v| shapes::Rect { size: domain.size.to_vector().to_point() }
56    .translate(domain.center().to_vector())
57    .sdf(v); // IPM boundary
58
59  let control_points = |rect: Rect<_, _>| {
60    let p = (0..lattice_density).map(move |x| _Float::from(x).unwrap() / _Float::from(lattice_density - 1).unwrap());
61    itertools::iproduct!(p.clone(), p)
62      .map(move |p| rect.origin + rect.size.to_vector().component_mul(p.into()))
63  };
64
65  let test = |v| line_search.optimize_normal(
66    |v| if domain.contains(v) { g(v) - f(v) } else { -boundary_constraint(v) },
67    v
68  );
69
70  !match lattice_density {
71    1 => test(domain.center()),
72    2..=u32::MAX => control_points(domain).any(test),
73    _ => panic!("Invalid perturbation grid density: {lattice_density}")
74  }
75}
76
77impl <_Float: Float + Signed + Sync> ADF<_Float> {
78  /// Create a new ADF instance. `max_depth` specifies maximum number of quadtree subdivisions;
79  /// `init` specifies initial sdf primitives.
80  pub fn new(max_depth: u8, init: Vec<Arc<dyn Fn(P2<_Float>) -> _Float>>) -> Self {
81    Self {
82      tree: Quadtree::new(max_depth, init),
83      ipm_gd_lattice_density: 1,
84      ipm_line_config: LineSearch::default()
85    }
86  }
87  /// Controls precision of primitive pruning in a bucket.
88  pub fn with_gd_lattice_density(mut self, density: u32) -> Self {
89    self.ipm_gd_lattice_density = density;
90    self
91  }
92  /// Underlying GD settings for the interior point method (a part of primitive pruning).
93  pub fn with_ipm_line_config(mut self, line_config: LineSearch<_Float>) -> Self {
94    self.ipm_line_config = line_config;
95    self
96  }
97  /*
98    Upon insertion of a new SDF primitive (`f`), this function tests whether it does
99    change the distance field within a certain domain (remember that it is considered changed
100    if and only if at least one point of `f` is lower than combined distance field (`g`) within
101    a certain domain `D` (tree node).
102    If no change is present, therefore updating `g` within the domain may be safely
103    skipped. However, it is imperfect: the test is only performed within a static square grid of
104    25 control points, thus sometimes yielding incorrect result, and generally being very slow.
105
106    Proposition. Use gradient descent (see `sdf_partialord`) in order to pick the control points
107    more carefully, and answer following questions:
108    \begin{align*}
109&f(\overrightarrow{v}) < g(\overrightarrow{v}), \forall\, \overrightarrow{v}\epsilon \,\mathfrak{D} &(1)\\
110&f(\overrightarrow{v}) > g(\overrightarrow{v}), \forall\, \overrightarrow{v}\epsilon \,\mathfrak{D} &(2)\\
111&\exists \overrightarrow{v}\epsilon \,\mathfrak{D}: f(\overrightarrow{v}) < g(\overrightarrow{v}) &(3)
112\end{align}
113
114    Proposition 2. Use interior point method in order to specify the boundary constraint of `D`
115
116    Update: implemented in [adf::sdf_partialord]
117   */
118
119  /// f(v) > g(v) forall v e D
120  #[deprecated] #[allow(unused)]
121  fn higher_all(
122    f: &(dyn Fn(Point2D<f64, WorldSpace>) -> f64),
123    g: &(dyn Fn(Point2D<f64, WorldSpace>) -> f64),
124    d: Rect<f64, WorldSpace>
125  ) -> bool {
126    let control_points = |rect: Rect<_, _>| {
127      let n = 5;
128      let p = (0..n).map(move |x| x as f64 / (n - 1) as f64);
129      itertools::iproduct!(p.clone(), p)
130        .map(move |p| rect.origin + rect.size.to_vector().component_mul(p.into()))
131    };
132
133    !control_points(d)
134      .any(|v| g(v) > f(v))
135  }
136
137  /// Add a new sdf primitive function.
138  pub fn insert_sdf_domain(&mut self, domain: Rect<_Float, WorldSpace>, f: Arc<dyn Fn(P2<_Float>) -> _Float + Send + Sync>) -> bool {
139    let change_exists = AtomicBool::new(false);
140
141    self.tree.traverse_managed_parallel(|node| {
142      // no intersection with domain
143      if !node.rect.intersects(&domain) {
144        return TraverseCommand::Skip;
145      }
146
147      // not a leaf node
148      if node.children.is_some() {
149        return TraverseCommand::Ok;
150      }
151
152      // f(v) > g(v) forall v e D, no refinement is required
153      if sdf_partialord(
154        f.as_ref(),
155        |p| node.data.as_slice().sdf(p),
156        node.rect,
157        self.ipm_gd_lattice_density,
158        self.ipm_line_config
159      ) {
160        return TraverseCommand::Skip;
161      }
162
163      // f(v) <= g(v) forall v e D, a minor optimization
164      if sdf_partialord(
165        |p| node.data.as_slice().sdf(p),
166        f.as_ref(),
167        node.rect,
168        self.ipm_gd_lattice_density,
169        self.ipm_line_config
170      ) {
171        node.data = vec![f.clone()];
172        change_exists.store(true, Ordering::Relaxed);
173        return TraverseCommand::Skip;
174      };
175
176      change_exists.store(true, Ordering::Relaxed);
177      const BUCKET_SIZE: usize = 3;
178
179      // remove SDF primitives, that do not affect the field within `D`
180      let prune = |data: &[Arc<dyn Fn(P2<_Float>) -> _Float>], rect| {
181        let mut g = vec![];
182        for (i, f) in data.iter().enumerate() {
183          let sdf_old = |p|
184            data.iter().enumerate()
185              .filter_map(|(j, f)| if i != j {
186                Some(f(p))
187              } else { None })
188              .fold(_Float::max_value() / (_Float::one() + _Float::one()), |a, b| a.min(b));
189          // there exists v e D, such that f(v) < g(v)
190          if !sdf_partialord(
191            f.as_ref(),
192            sdf_old,
193            rect,
194            self.ipm_gd_lattice_density,
195            self.ipm_line_config
196          ) {
197            g.push(f.clone())
198          }
199        };
200        g
201      };
202
203      // max tree depth is reached, just append the primitive
204      if node.depth == node.max_depth || node.data.len() < BUCKET_SIZE {
205
206        node.data.push(f.clone());
207        //node.data = prune(node.data.as_slice(), node.rect);
208
209      } /*else if node.data.len() < BUCKET_SIZE {
210
211        //node.data = prune(node.data.as_slice(), node.rect);
212        //if !Self::higher_all(f.as_ref(), &node.sdf_vec(), node.rect) {
213        //  node.data.push(f.clone());
214        //}
215
216      }*/
217      else { // max bucket size is reached, subdivide
218
219        let mut g = node.data.clone();
220        g.push(f.clone());
221
222        node.subdivide(|rect_ch| prune(g.as_slice(), rect_ch));
223        /*node.subdivide(|rect_ch| prune(&g, rect_ch))
224          .as_deref_mut()
225          .unwrap()
226          .into_iter()
227          .for_each(|child| {
228            child.insert_sdf_domain(domain, f.clone());
229          });*/
230      }
231      //TODO: despite issuing Skip, newly divided node continues to be explored. Rework Traverse API
232      TraverseCommand::Skip
233    });
234
235    change_exists.load(Ordering::SeqCst)
236  }
237
238  /// # Safety
239  /// Nobody is safe
240  pub unsafe fn as_mut(&self) -> &mut Self {
241    let ptr = self as *const _ as usize;
242    &mut *(ptr as *const Self as *mut _)
243  }
244}
245
246impl <_Float: Float> SDF<_Float> for ADF<_Float> {
247  fn sdf(&self, pixel: P2<_Float>) -> _Float {
248    match self.tree.pt_to_node(pixel) {
249      Some(node) => node.data.as_slice().sdf(pixel),
250      None => self.tree.data.as_slice().sdf(pixel),
251    }}}
252
253impl <_Float: Float> BoundingBox<_Float> for ADF<_Float> {
254  fn bounding_box(&self) -> Box2D<_Float, WorldSpace> {
255    Box2D::new(
256      P2::splat(_Float::zero()),
257      P2::splat(_Float::one())
258    )}}
259
260impl <_Float: Float> Debug for ADF<_Float> {
261  fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
262    use humansize::{FileSize, file_size_opts as options};
263
264    let mut total_nodes = 0u64;
265    let mut total_size = 0usize;
266    let mut max_depth = 0u8;
267    self.tree.traverse(&mut |node| {
268      total_nodes += 1;
269      total_size += std::mem::size_of::<Self>()
270        + node.data.capacity() * std::mem::size_of::<Arc<dyn Fn(P2<f64>) -> f64>>();
271      max_depth = (max_depth).max(node.depth);
272      Ok(())
273    }).ok();
274    f.debug_struct("ADF")
275      .field("total_nodes", &total_nodes)
276      .field("max_depth", &max_depth)
277      .field("size", &total_size.file_size(options::BINARY).unwrap())
278      .finish()
279  }
280}
space_filling/solver/adf/mod.rs

space_filling/solver/adf/
mod.rs
