visual_odometry_rs/core/multires.rs
1// This Source Code Form is subject to the terms of the Mozilla Public
2// License, v. 2.0. If a copy of the MPL was not distributed with this
3// file, You can obtain one at http://mozilla.org/MPL/2.0/.
4
5//! Helper functions for generation of multi-resolution data.
6
7use nalgebra::{DMatrix, Scalar};
8
9use crate::core::gradient;
10
11/// Recursively generate a pyramid of matrices where each following level
12/// is half the previous resolution, computed with the mean of each 2x2 block.
13///
14/// It consumes the original matrix to keep it as first level of the pyramid without copy.
15/// If you need it later, simply use something like `pyramid[0]`.
16///
17/// PS: since we are using 2x2 blocs,
18/// border information is lost for odd resolutions.
19/// Some precision is also left to keep the pyramid data as `u8`.
20#[allow(clippy::cast_possible_truncation)]
21pub fn mean_pyramid(max_levels: usize, mat: DMatrix<u8>) -> Vec<DMatrix<u8>> {
22 limited_sequence(max_levels, mat, |m| {
23 halve(m, |a, b, c, d| {
24 let a = u16::from(a);
25 let b = u16::from(b);
26 let c = u16::from(c);
27 let d = u16::from(d);
28 ((a + b + c + d) / 4) as u8
29 })
30 })
31}
32
33/// Recursively apply a function transforming an image
34/// until it's not possible anymore or the max length is reached.
35///
36/// Using `max_length = 0` has the same effect than `max_length = 1` since
37/// the result vector contains always at least one matrix (the init matrix).
38pub fn limited_sequence<F: Fn(&T) -> Option<T>, T>(max_length: usize, data: T, f: F) -> Vec<T> {
39 let mut length = 1;
40 let f_limited = |x: &T| {
41 if length < max_length {
42 length += 1;
43 f(x)
44 } else {
45 None
46 }
47 };
48 sequence(data, f_limited)
49}
50
51/// Recursively apply a function transforming data
52/// until it's not possible anymore.
53pub fn sequence<F: FnMut(&T) -> Option<T>, T>(data: T, mut f: F) -> Vec<T> {
54 let mut s = Vec::new();
55 s.push(data);
56 while let Some(new_data) = f(s.last().unwrap()) {
57 s.push(new_data);
58 }
59 s
60}
61
62/// Halve the resolution of a matrix by applying a function to each 2x2 block.
63///
64/// If one size of the matrix is < 2 then this function returns None.
65/// If one size is odd, its last line/column is dropped.
66#[allow(clippy::many_single_char_names)]
67pub fn halve<F, T, U>(mat: &DMatrix<T>, f: F) -> Option<DMatrix<U>>
68where
69 F: Fn(T, T, T, T) -> U,
70 T: Scalar,
71 U: Scalar,
72{
73 let (r, c) = mat.shape();
74 let half_r = r / 2;
75 let half_c = c / 2;
76 if half_r == 0 || half_c == 0 {
77 None
78 } else {
79 let half_mat = DMatrix::<U>::from_fn(half_r, half_c, |i, j| {
80 let a = mat[(2 * i, 2 * j)];
81 let b = mat[(2 * i + 1, 2 * j)];
82 let c = mat[(2 * i, 2 * j + 1)];
83 let d = mat[(2 * i + 1, 2 * j + 1)];
84 f(a, b, c, d)
85 });
86 Some(half_mat)
87 }
88}
89
90// Gradients stuff ###################################################
91
92/// Compute centered gradients norm at each resolution from
93/// the image at the higher resolution.
94///
95/// As a consequence there is one less level in the gradients pyramid.
96pub fn gradients_squared_norm(multires_mat: &[DMatrix<u8>]) -> Vec<DMatrix<u16>> {
97 let nb_levels = multires_mat.len();
98 multires_mat
99 .iter()
100 .take(nb_levels - 1)
101 .map(|mat| {
102 halve(mat, gradient::bloc_squared_norm)
103 .expect("There is an issue in gradients_squared_norm")
104 })
105 .collect()
106}
107
108/// Compute centered gradients at each resolution from
109/// the image at the higher resolution.
110///
111/// As a consequence there is one less level in the gradients pyramid.
112pub fn gradients_xy(multires_mat: &[DMatrix<u8>]) -> Vec<(DMatrix<i16>, DMatrix<i16>)> {
113 // TODO: maybe it would be better to return Vec<DMatrix<(i16,i16)>>,
114 // to colocate the x and y gradient and do only one "halve" call?
115 let nb_levels = multires_mat.len();
116 multires_mat
117 .iter()
118 .take(nb_levels - 1)
119 .map(|mat| {
120 (
121 halve(mat, gradient::bloc_x).expect("There is an issue in gradients_xy x."),
122 halve(mat, gradient::bloc_y).expect("There is an issue in gradients_xy y."),
123 )
124 })
125 .collect()
126}