spatial_decomposition/
kong_mount_roscoe.rs

1use std::num::NonZeroUsize;
2
3use approx::RelativeEq;
4use num_traits::AsPrimitive;
5use simba::scalar::RealField;
6
7use crate::Rectangle;
8
9#[allow(non_camel_case_types)]
10enum Decomposition<F> {
11    row(F),
12    col(F),
13}
14
15#[allow(non_snake_case)]
16#[derive(Clone, Debug, RelativeEq, PartialEq)]
17struct KongMountRoscoeValues<F> {
18    p: F,
19    h1: F,
20    h2: F,
21    k1: F,
22    k2: F,
23    A: F,
24    B: F,
25    C: F,
26    S: F,
27}
28
29impl<F> KongMountRoscoeValues<F>
30where
31    F: 'static + Copy + RealField,
32    usize: num_traits::cast::AsPrimitive<F>,
33{
34    #[allow(non_snake_case)]
35    fn calculate(A: F, B: F, n_subdomains: usize) -> Self {
36        let p: F = n_subdomains.as_();
37
38        let h1 = (A * p / B).sqrt().floor();
39        let h2 = (A * p / B).sqrt().ceil();
40        let k1 = (B * p / A).sqrt().floor();
41        let k2 = (B * p / A).sqrt().ceil();
42
43        let mut C = p;
44        let mut S = F::zero();
45        for k in 1..n_subdomains {
46            for h in 1..n_subdomains {
47                let k = k.as_();
48                let h = h.as_();
49                if (h - F::one()) * (k - F::one()) < p {
50                    let ci = (A / h).max(B / k);
51                    C = C.min(ci);
52                }
53                if p < (h + F::one()) * (k + F::one()) {
54                    let si = (A / h).min(B / k);
55                    S = S.max(si);
56                }
57            }
58        }
59
60        #[rustfmt::skip]
61        return KongMountRoscoeValues {
62            p,
63            h1, h2, k1, k2,
64            A,  B,  C,  S,
65        };
66    }
67}
68
69impl<F> Decomposition<F> {
70    #[rustfmt::skip]
71    fn figure_out(values: &KongMountRoscoeValues<F>) -> Option<Decomposition<F>>
72    where
73        F: RealField + Copy,
74    {
75        use Decomposition::*;
76        let KongMountRoscoeValues {
77            p,
78            h1, h2, k1, k2,
79            A,  B,  C,  S,
80        } = values.clone();
81        let h0 = (p / k2).floor();
82        let h3 = (p / k1).ceil();
83        let k0 = (p / h2).floor();
84        let k3 = (p / h1).ceil();
85
86        if h1==h2 && k1==k2 {return Some(col(h1));}
87
88        //  Lemma 2.9
89        if h0 < h1 && B / k0 <= S {return Some(row(h2));}
90        if k3 > k2 && A / h0 >= C {return Some(row(h1));}
91        if h0 < h1 && A / h0 <= S {return Some(col(k2));}
92        if h3 > h2 && A / h3 >= C {return Some(col(k1));}
93
94        // Lemma 2.10
95        if h0 == h1 && A / h1 <= S && A / h2 >= C {return Some(col(k2));}
96        if h3 == h2 && A / h1 <= S && A / h2 >= C {return Some(col(k1));}
97        if k0 == k1 && B / k1 <= S && B / k2 >= C {return Some(row(h2));}
98        if k3 == k2 && B / k1 <= S && B / k2 >= C {return Some(row(h1));}
99
100        // TODO is this exhaustive?
101        // -> probably not!
102
103        None
104    }
105
106    #[allow(non_snake_case)]
107    fn generate_rectangles(
108        &self,
109        p: F,
110        rectangle: &Rectangle<F>,
111    ) -> impl IntoIterator<Item = Rectangle<F>> + use<F>
112    where
113        F: 'static + RealField + Copy,
114        F: num_traits::cast::AsPrimitive<usize>,
115        usize: num_traits::cast::AsPrimitive<F>,
116    {
117        let min = rectangle.min;
118        let B = rectangle.max[0] - rectangle.min[0];
119        let A = rectangle.max[1] - rectangle.min[1];
120
121        use Decomposition::*;
122        match self {
123            row(hrow) => {
124                let hrow: F = *hrow;
125                let dx_row = A / hrow;
126                let (n_rows1, n_rows2, n_cols1, n_cols2): (usize, usize, F, F) =
127                    if (p / hrow).round() == p / hrow {
128                        (hrow.as_(), 0, (p / hrow), F::zero())
129                    } else {
130                        (
131                            (p - hrow * (p / hrow).floor()).round().as_(),
132                            (hrow * (p / hrow).ceil() - p).round().as_(),
133                            (p / hrow).ceil(),
134                            (p / hrow).floor(),
135                        )
136                    };
137                let dx_col1 = B / n_cols1;
138                let dx_col2 = B / n_cols2;
139                let n_cols1: usize = n_cols1.as_();
140                let n_cols2: usize = n_cols2.as_();
141
142                let rects1 =
143                    create_rectangles(0..n_rows1, 0..n_cols1, dx_row, dx_col1, min[1], min[0]);
144                let rects2 = create_rectangles(
145                    n_rows1..n_rows1 + n_rows2,
146                    0..n_cols2,
147                    dx_row,
148                    dx_col2,
149                    min[1],
150                    min[0],
151                );
152                rects1.into_iter().chain(rects2)
153            }
154            col(kcol) => {
155                let kcol: F = *kcol;
156                let dx_col = B / kcol;
157                let (n_cols1, n_cols2, n_rows1, n_rows2): (usize, usize, F, F) =
158                    if (p / kcol).round() == p / kcol {
159                        (kcol.as_(), 0, p / kcol, F::zero())
160                    } else {
161                        (
162                            (p - kcol * (p / kcol).floor()).round().as_(),
163                            (kcol * (p / kcol).ceil() - p).round().as_(),
164                            (p / kcol).ceil(),
165                            (p / kcol).floor(),
166                        )
167                    };
168                let dx_row1 = A / n_rows1;
169                let dx_row2 = A / n_rows2;
170                let n_rows1: usize = n_rows1.as_();
171                let n_rows2: usize = n_rows2.as_();
172
173                let rects1 =
174                    create_rectangles(0..n_rows1, 0..n_cols1, dx_row1, dx_col, min[1], min[0]);
175                let rects2 = create_rectangles(
176                    0..n_rows2,
177                    n_cols1..n_cols1 + n_cols2,
178                    dx_row2,
179                    dx_col,
180                    min[1],
181                    min[0],
182                );
183                rects1.into_iter().chain(rects2)
184            }
185        }
186    }
187}
188
189fn create_rectangles<F>(
190    n_rows_range: std::ops::Range<usize>,
191    n_cols_range: std::ops::Range<usize>,
192    dx_row: F,
193    dx_col: F,
194    x0_row: F,
195    x0_col: F,
196) -> impl IntoIterator<Item = Rectangle<F>>
197where
198    F: 'static + RealField + Copy,
199    F: num_traits::cast::AsPrimitive<usize>,
200    usize: num_traits::cast::AsPrimitive<F>,
201{
202    (n_rows_range).flat_map(move |n| {
203        (n_cols_range).clone().map(move |m| {
204            let n: F = n.as_();
205            let m: F = m.as_();
206            let min = [x0_col + m * dx_col, x0_row + n * dx_row];
207            let max = [
208                x0_col + (m + F::one()) * dx_col,
209                x0_row + (n + F::one()) * dx_row,
210            ];
211            Rectangle { min, max }
212        })
213    })
214}
215
216/// Partitions a rectangle into multiple smaller rectangles
217///
218/// This algorithms follows the paper by [Kong, Mount and Roscoe](https://scispace.com/pdf/the-decomposition-of-a-rectangle-into-rectangles-of-minimal-3whu99wjdy.pdf)
219///
220/// This algorithms divides a given rectangle into multiple smaller rectangles and minimizes the
221/// maximum rectangle perimeter.
222///
223/// ```
224/// use spatial_decomposition::{kmr_decompose, Rectangle};
225///
226/// let domain = Rectangle {
227///     min: [0., 40.],
228///     max: [100., 240.],
229/// };
230///
231/// let n_subdomains = 9;
232/// let subdomains = kmr_decompose(&domain, n_subdomains.try_into().unwrap());
233///
234/// assert_eq!(subdomains.len(), 9);
235///
236/// for subdomain in subdomains {
237/// }
238/// ```
239///
240/// ## Examples
241///
242/// Cases of very long/wide rectangles satisfying n_subdomain < max(a/b, b/a)
243/// will be split intuitively along the long axis.
244///
245/// ```text
246///                    B = 90
247///        ┌─────────┬─────────┬─────────┐
248///        │         │         │         │
249/// A = 20 │         │         │         │
250///        │         │         │         │
251///        └─────────┴─────────┴─────────┘
252/// ```
253///
254/// ```
255/// # use spatial_decomposition::{kmr_decompose, Rectangle};
256/// # use approx::assert_abs_diff_eq;
257/// let domain = Rectangle {
258///     min: [0.0; 2],
259///     max: [90.0, 20.0],
260/// };
261/// let subdomains = kmr_decompose(
262///     &domain,
263///     3.try_into().unwrap()
264/// );
265/// // ┌─────────┬─
266/// // │         │
267/// // │         │
268/// // │         │
269/// // └─────────┴─
270/// assert_abs_diff_eq!(
271///     subdomains[0],
272///     &Rectangle {
273///         min: [0.0; 2],
274///         max: [30.0, 20.0],
275/// });
276///
277/// //           ┬─────────┬─
278/// //           │         │
279/// //           │         │
280/// //           │         │
281/// //           ┴─────────┴─
282/// assert_abs_diff_eq!(
283///     subdomains[1],
284///     &Rectangle {
285///         min: [30.0, 0.0],
286///         max: [60.0, 20.0],
287/// });
288///
289/// //                     ┬─────────┐
290/// //                     │         │
291/// //                     │         │
292/// //                     │         │
293/// //                     ┴─────────┘
294/// assert_abs_diff_eq!(
295///     subdomains[2],
296///     &Rectangle {
297///         min: [60.0, 0.0],
298///         max: [90.0, 20.0],
299/// });
300/// ```
301#[allow(non_snake_case)]
302pub fn kmr_decompose<F>(rectangle: &Rectangle<F>, n_subdomains: NonZeroUsize) -> DecomposedDomain<F>
303where
304    F: 'static + Copy + RealField,
305    F: num_traits::cast::AsPrimitive<usize>,
306    usize: num_traits::cast::AsPrimitive<F>,
307{
308    let n_subdomains = n_subdomains.get();
309    if n_subdomains == 1 {
310        return vec![rectangle.clone()];
311    }
312
313    // Cover the very wide/long cases where n_subdomain < max(a/b, b/a)
314    //                  a
315    //   ┌─────────┬─────────┬─────────┐
316    //   │         │         │         │
317    // b │         │         │         │
318    //   │         │         │         │
319    //   └─────────┴─────────┴─────────┘
320    let n_subdomains_float: F = n_subdomains.as_();
321    let B = rectangle.max[0] - rectangle.min[0];
322    let A = rectangle.max[1] - rectangle.min[1];
323    let ratio_max = (B / A).max(A / B);
324    if n_subdomains_float <= ratio_max {
325        if B >= A {
326            let dx = B / n_subdomains_float;
327            return (0..n_subdomains)
328                .map(|n| Rectangle {
329                    min: [rectangle.min[0] + n.as_() * dx, rectangle.min[1]],
330                    max: [rectangle.min[0] + (n + 1).as_() * dx, rectangle.max[1]],
331                })
332                .collect();
333        } else {
334            let dx = A / n_subdomains_float;
335            return (0..n_subdomains)
336                .map(|n| Rectangle {
337                    min: [rectangle.min[0], rectangle.min[1] + n.as_() * dx],
338                    max: [rectangle.max[0], rectangle.min[1] + (n + 1).as_() * dx],
339                })
340                .collect();
341        };
342    }
343
344    let kmr_values = KongMountRoscoeValues::calculate(A, B, n_subdomains);
345    let decomposition = Decomposition::figure_out(&kmr_values);
346
347    decomposition
348        .map(|d| d.generate_rectangles(n_subdomains_float, rectangle))
349        .unwrap()
350        .into_iter()
351        .collect()
352}
353
354/// Returned from decomposition methods.
355pub type DecomposedDomain<F> = Vec<Rectangle<F>>;
356
357/// Error variants of decomposition or digitization
358#[derive(thiserror::Error, Debug)]
359pub enum Error {
360    #[error("Could not find suitable subdomain for given digit")]
361    Decompose,
362    #[error("Could not assign index to subdomain")]
363    Indexing,
364}
365
366/// Returned from digitization methods
367pub type SortedDigits<F, I> = Vec<(Rectangle<F>, Vec<(I, Rectangle<F>)>)>;
368
369pub type Result<T> = std::result::Result<T, Error>;
370
371/// Assign given digits to a number of subdomains generated by the [kmr_decompose] method.
372///
373/// It is the first algorithm presented in their paper.
374/// This particular implementation is `O(n*p)` where `n` is the number of digits and `p` is the
375/// number of subdomains generated.
376pub fn kmr_digitize_1<F, I>(
377    rectangle: &Rectangle<F>,
378    n_subdomains: NonZeroUsize,
379    digits: impl IntoIterator<Item = (I, Rectangle<F>)>,
380) -> Result<SortedDigits<F, I>>
381where
382    F: 'static + Copy + RealField,
383    F: num_traits::cast::AsPrimitive<usize>,
384    usize: num_traits::cast::AsPrimitive<F>,
385    I: 'static,
386{
387    let subdomains = kmr_decompose(rectangle, n_subdomains);
388    let mut res: Vec<_> = subdomains
389        .iter()
390        .map(|subdomain| (subdomain.clone(), Vec::new()))
391        .collect();
392    for (digit, rect) in digits.into_iter() {
393        let index = kmr_digitize_1_single(&subdomains, &rect)?;
394        res[index].1.push((digit, rect));
395    }
396    Ok(res)
397}
398
399/// Assigns an index to a given subspace.
400///
401/// Singular version of [kmr_digitize_1].
402pub fn kmr_digitize_1_single<F>(
403    decomposed_domain: &DecomposedDomain<F>,
404    subspace: &Rectangle<F>,
405) -> Result<usize>
406where
407    F: 'static + Copy + RealField,
408    F: num_traits::cast::AsPrimitive<usize>,
409    usize: num_traits::cast::AsPrimitive<F>,
410{
411    for (n_subdomain, subdomain) in decomposed_domain.iter().enumerate() {
412        let middle = [
413            (subspace.min[0] + subspace.max[0]) / (F::one() + F::one()),
414            (subspace.min[1] + subspace.max[1]) / (F::one() + F::one()),
415        ];
416        if middle[0] <= subdomain.max[0]
417            && middle[1] <= subdomain.max[1]
418            && middle[0] >= subdomain.min[0]
419            && middle[1] >= subdomain.min[1]
420        {
421            return Ok(n_subdomain);
422        }
423    }
424    Err(Error::Indexing)
425}
426
427#[test]
428fn kmr_digitize_square_2x2_in_4() {
429    let domain = Rectangle {
430        min: [0.0; 2],
431        max: [100.0; 2],
432    };
433    let digits = kmr_decompose(&domain, 16.try_into().unwrap())
434        .into_iter()
435        .enumerate();
436    let sorted = kmr_digitize_1(&domain, 4.try_into().unwrap(), digits).unwrap();
437    for (_, voxels) in sorted.iter() {
438        assert_eq!(voxels.len(), 4);
439    }
440}
441
442#[test]
443fn kmr_decompose_identity() {
444    let rectangle = Rectangle {
445        min: [0.0; 2],
446        max: [10.0; 2],
447    };
448    let rects = kmr_decompose(&rectangle, 1.try_into().unwrap());
449    assert_eq!(rects.len(), 1);
450    assert_eq!(rects[0], rectangle);
451}
452
453#[test]
454fn kmr_decompose_very_wide() {
455    use approx::assert_abs_diff_eq;
456    let rectangle = Rectangle {
457        min: [0.0; 2],
458        max: [100.0, 10.0],
459    };
460    let rects = kmr_decompose(&rectangle, 4.try_into().unwrap());
461    assert_eq!(rects.len(), 4);
462    assert_abs_diff_eq!(
463        rects[0],
464        Rectangle {
465            min: [0.0; 2],
466            max: [25.0, 10.0],
467        }
468    );
469    assert_abs_diff_eq!(
470        rects[1],
471        Rectangle {
472            min: [25.0, 0.0],
473            max: [50.0, 10.0],
474        }
475    );
476    assert_abs_diff_eq!(
477        rects[2],
478        Rectangle {
479            min: [50.0, 0.0],
480            max: [75.0, 10.0],
481        }
482    );
483    assert_abs_diff_eq!(
484        rects[3],
485        Rectangle {
486            min: [75.0, 0.0],
487            max: [100.0, 10.0],
488        }
489    );
490}
491
492#[test]
493fn kmr_decompose_very_long() {
494    use approx::assert_abs_diff_eq;
495    let rectangle = Rectangle {
496        min: [-10f32, -200f32],
497        max: [10f32, 200f32],
498    };
499    let rects = kmr_decompose(&rectangle, 10.try_into().unwrap());
500    let dx = (rectangle.max[1] - rectangle.min[1]) / 10f32;
501    for (n, rect) in rects.into_iter().enumerate() {
502        assert_abs_diff_eq!(
503            rect,
504            Rectangle {
505                min: [rectangle.min[0], rectangle.min[1] + n as f32 * dx],
506                max: [rectangle.max[0], rectangle.min[1] + (n + 1) as f32 * dx],
507            }
508        );
509    }
510}
511
512#[test]
513fn kmr_decompose_5x3_in_7() {
514    let rectangle = Rectangle {
515        min: [0.0; 2],
516        max: [5.0, 3.0],
517    };
518    let rects = kmr_decompose(&rectangle, 7.try_into().unwrap());
519    assert_eq!(rects.len(), 7);
520
521    for i in 0..4 {
522        let dx = (rectangle.max[0] - rectangle.min[0]) / 4.;
523        let i = i as f64;
524        assert!(rects.contains(&Rectangle {
525            min: [i * dx, 0.],
526            max: [(i + 1.) * dx, 1.5],
527        }));
528    }
529    for i in 0..3 {
530        let dx = (rectangle.max[0] - rectangle.min[0]) / 3.;
531        let i = i as f64;
532        assert!(rects.contains(&Rectangle {
533            min: [i * dx, 1.5],
534            max: [(i + 1.) * dx, 3.]
535        }))
536    }
537}
538
539#[test]
540fn kmr_decompose_6x6_in_14() {
541    let rectangle = Rectangle {
542        min: [-60., 0.],
543        max: [0., 60.],
544    };
545    let rects = kmr_decompose(&rectangle, 14.try_into().unwrap());
546    assert_eq!(rects.len(), 14);
547    for i in 0..5 {
548        for j in 0..2 {
549            let dx = 60. / 5.;
550            let i = i as f64;
551            let j = j as f64;
552            let dy = 60. / 3.;
553            let r = Rectangle {
554                min: [i * dx - 60., j * dy],
555                max: [(i + 1.) * dx - 60., (j + 1.) * dy],
556            };
557            assert!(rects.contains(&r));
558        }
559    }
560    for i in 0..4 {
561        let i = i as f64;
562        let j = 2.;
563        let dx = 60. / 4.;
564        let dy = 60. / 3.;
565        let r = Rectangle {
566            min: [i * dx - 60., j * dy],
567            max: [(i + 1.) * dx - 60., (j + 1.) * dy],
568        };
569        assert!(rects.contains(&r));
570    }
571}
572
573#[test]
574fn kmr_decompose_square_into_4() {
575    let rectangle = Rectangle {
576        min: [-40.0; 2],
577        max: [-20.0; 2],
578    };
579    let rects = kmr_decompose(&rectangle, 4.try_into().unwrap());
580    assert_eq!(rects.len(), 4);
581    assert!(rects.contains(&Rectangle {
582        min: [-40.; 2],
583        max: [-30.; 2]
584    }));
585    assert!(rects.contains(&Rectangle {
586        min: [-40.0, -30.0],
587        max: [-30.0, -20.0]
588    }));
589    assert!(rects.contains(&Rectangle {
590        min: [-30.0, -40.0],
591        max: [-20.0, -30.0]
592    }));
593    assert!(rects.contains(&Rectangle {
594        min: [-30.0; 2],
595        max: [-20.0; 2]
596    }));
597}
598
599#[test]
600fn kmr_decompose_square_into_5() {
601    let rectangle = Rectangle {
602        min: [0.0; 2],
603        max: [100.0; 2],
604    };
605    let subdomains = kmr_decompose(&rectangle, 5.try_into().unwrap());
606    assert_eq!(subdomains.len(), 5);
607    for s in subdomains.iter() {
608        println!("{s:7.2?}");
609    }
610    assert!(subdomains.contains(&Rectangle {
611        min: [0.0; 2],
612        max: [50.0, 100.0 / 3.],
613    }));
614    assert!(subdomains.contains(&Rectangle {
615        min: [0.0, 100. / 3.],
616        max: [50., 100. / 3. * 2.]
617    }));
618    assert!(subdomains.contains(&Rectangle {
619        min: [0.0, 100. / 3. * 2.],
620        max: [50., 100.]
621    }));
622    assert!(subdomains.contains(&Rectangle {
623        min: [50.0, 0.],
624        max: [100., 50.]
625    }));
626    assert!(subdomains.contains(&Rectangle {
627        min: [50.0, 50.],
628        max: [100., 100.]
629    }));
630}
spatial_decomposition/kong_mount_roscoe.rs

spatial_decomposition/
kong_mount_roscoe.rs
