conspire/geometry/ntree/defeature/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::geometry::ntree::{Orthotree, node::Kind, subdivide::insert_bit};
5use std::{
6    cmp::Reverse,
7    collections::{BinaryHeap, HashMap},
8    hash::Hash,
9};
10
11impl<const D: usize, const L: usize, const M: usize, const N: usize, T, U, V>
12    Orthotree<D, L, M, N, T, U, V>
13where
14    T: Copy + Into<usize>,
15    U: Copy + From<usize> + Into<usize>,
16    V: Copy + Eq + Hash,
17{
18    pub fn defeature(&mut self, minimum: usize) {
19        let count = self.len();
20        let mut parent: Vec<usize> = (0..count).collect();
21        let leaves: Vec<usize> = (0..count)
22            .filter(|&i| self.nodes[i].is_leaf() && self.nodes[i].value.is_some())
23            .collect();
24        let mut edges: Vec<(usize, usize, usize)> = Vec::new();
25        for &leaf in &leaves {
26            let value = self.nodes[leaf].value.unwrap();
27            let length: usize = self.nodes[leaf].length.into();
28            for face in 0..M {
29                if let Some(neighbor) = self.nodes[leaf].facets[face] {
30                    let mut others = Vec::new();
31                    self.face_leaves(neighbor.into(), face ^ 1, &mut others);
32                    for other in others {
33                        if let Some(adjacent) = self.nodes[other].value {
34                            if adjacent == value {
35                                union(&mut parent, leaf, other);
36                            } else {
37                                let span: usize = self.nodes[other].length.into();
38                                edges.push((leaf, other, length.min(span).pow((D - 1) as u32)));
39                            }
40                        }
41                    }
42                }
43            }
44        }
45        let mut volume: HashMap<usize, usize> = HashMap::new();
46        let mut value: HashMap<usize, V> = HashMap::new();
47        for &leaf in &leaves {
48            let root = find(&mut parent, leaf);
49            let length: usize = self.nodes[leaf].length.into();
50            *volume.entry(root).or_default() += length.pow(D as u32);
51            value
52                .entry(root)
53                .or_insert_with(|| self.nodes[leaf].value.unwrap());
54        }
55        let mut adjacency: HashMap<usize, HashMap<usize, usize>> = HashMap::new();
56        for (leaf, other, span) in edges {
57            let (a, b) = (find(&mut parent, leaf), find(&mut parent, other));
58            if a != b {
59                *adjacency.entry(a).or_default().entry(b).or_default() += span;
60                *adjacency.entry(b).or_default().entry(a).or_default() += span;
61            }
62        }
63        let mut queue: BinaryHeap<Reverse<(usize, usize)>> = volume
64            .iter()
65            .filter(|&(_, &size)| size < minimum)
66            .map(|(&root, &size)| Reverse((size, root)))
67            .collect();
68        while let Some(Reverse((size, root))) = queue.pop() {
69            if value.get(&root).is_none_or(|_| volume[&root] != size) || size >= minimum {
70                continue;
71            }
72            let neighbors = match adjacency.get(&root) {
73                Some(map) if !map.is_empty() => map,
74                _ => continue,
75            };
76            let mut by_value: HashMap<V, usize> = HashMap::new();
77            for (other, &span) in neighbors {
78                *by_value.entry(value[other]).or_default() += span;
79            }
80            let into = *by_value.iter().max_by_key(|&(_, &span)| span).unwrap().0;
81            let group: Vec<usize> = neighbors
82                .keys()
83                .copied()
84                .filter(|other| value[other] == into)
85                .chain([root])
86                .collect();
87            group[..group.len() - 1]
88                .iter()
89                .for_each(|&other| union(&mut parent, root, other));
90            let root = find(&mut parent, root);
91            let merged_volume = group.iter().map(|node| volume[node]).sum();
92            let mut merged: HashMap<usize, usize> = HashMap::new();
93            for node in &group {
94                if let Some(map) = adjacency.remove(node) {
95                    for (other, span) in map {
96                        let other = find(&mut parent, other);
97                        if other != root {
98                            *merged.entry(other).or_default() += span;
99                        }
100                    }
101                }
102            }
103            group.iter().filter(|&&node| node != root).for_each(|node| {
104                volume.remove(node);
105                value.remove(node);
106            });
107            for other in merged.keys() {
108                let map = adjacency.get_mut(other).unwrap();
109                let span: usize = group.iter().filter_map(|node| map.remove(node)).sum();
110                map.insert(root, span);
111            }
112            volume.insert(root, merged_volume);
113            value.insert(root, into);
114            adjacency.insert(root, merged);
115            if merged_volume < minimum {
116                queue.push(Reverse((merged_volume, root)));
117            }
118        }
119        leaves.iter().for_each(|&leaf| {
120            self.nodes[leaf].value = Some(value[&find(&mut parent, leaf)]);
121        });
122    }
123    fn face_leaves(&self, index: usize, face: usize, out: &mut Vec<usize>) {
124        match &self.nodes[index].kind {
125            Kind::Leaf => out.push(index),
126            Kind::Tree(orthants) => {
127                let (axis, side) = (face >> 1, face & 1);
128                for i in 0..L {
129                    let child = orthants[insert_bit(i, axis, side)].into();
130                    self.face_leaves(child, face, out);
131                }
132            }
133        }
134    }
135}
136
137fn find(parent: &mut [usize], mut i: usize) -> usize {
138    while parent[i] != i {
139        parent[i] = parent[parent[i]];
140        i = parent[i];
141    }
142    i
143}
144
145fn union(parent: &mut [usize], a: usize, b: usize) {
146    let (a, b) = (find(parent, a), find(parent, b));
147    if a != b {
148        parent[a] = b;
149    }
150}
conspire/geometry/ntree/defeature/mod.rs

conspire/geometry/ntree/defeature/
mod.rs
