Skip to main content

projective_grid/hex/
rectify.rs

1//! Single global homography from hex grid corners.
2//!
3//! Computes a global projective mapping between a rectified coordinate system
4//! (uniform hex lattice spacing) and the original image. Suitable when lens
5//! distortion is negligible.
6
7use crate::float_helpers::lit;
8use crate::grid_index::GridIndex;
9use crate::homography::{estimate_homography, Homography};
10use crate::Float;
11use nalgebra::Point2;
12use std::collections::HashMap;
13
14fn sqrt3_half<F: Float>() -> F {
15    lit::<F>(3.0).sqrt() / lit::<F>(2.0)
16}
17
18#[non_exhaustive]
19#[derive(thiserror::Error, Debug)]
20pub enum HexRectifyError {
21    #[error("not enough grid corners with positions (need >= 4, got {got})")]
22    NotEnoughPoints { got: usize },
23    #[error("homography estimation failed")]
24    HomographyFailed,
25    #[error("homography not invertible")]
26    NonInvertible,
27}
28
29/// A global homography mapping between rectified hex grid space and image space.
30///
31/// Rectified coordinates map axial `(q, r)` to 2D using:
32/// ```text
33/// x = px_per_cell * (q + r * 0.5)
34/// y = px_per_cell * (r * sqrt(3) / 2)
35/// ```
36#[derive(Clone, Debug)]
37pub struct HexGridHomography<F: Float = f32> {
38    /// Maps rectified coordinates to image coordinates.
39    pub h_img_from_rect: Homography<F>,
40    /// Maps image coordinates to rectified coordinates.
41    pub h_rect_from_img: Homography<F>,
42    /// Axial bounding box (with margin).
43    pub min_q: i32,
44    pub min_r: i32,
45    pub max_q: i32,
46    pub max_r: i32,
47    /// Rectified pixels per grid cell edge.
48    pub px_per_cell: F,
49    /// Rectified image dimensions.
50    pub rect_width: usize,
51    pub rect_height: usize,
52
53    /// Pixel-space origin offset subtracted during construction.
54    /// Needed by [`axial_to_rect`](Self::axial_to_rect) to produce coordinates
55    /// in the same frame as the stored homography.
56    x_offset: F,
57    y_offset: F,
58}
59
60impl<F: Float> HexGridHomography<F> {
61    /// Compute a global homography from hex grid corners to a rectified coordinate system.
62    ///
63    /// - `corners`: map from axial grid index `(q=i, r=j)` to image position.
64    /// - `px_per_cell`: rectified pixels per grid cell edge.
65    /// - `margin_cells`: extra margin around the grid bounding box (in cell units).
66    pub fn from_corners(
67        corners: &HashMap<GridIndex, Point2<F>>,
68        px_per_cell: F,
69        margin_cells: F,
70    ) -> Result<Self, HexRectifyError> {
71        if corners.len() < 4 {
72            return Err(HexRectifyError::NotEnoughPoints { got: corners.len() });
73        }
74
75        // Find axial bounding box
76        let (mut min_q, mut min_r) = (i32::MAX, i32::MAX);
77        let (mut max_q, mut max_r) = (i32::MIN, i32::MIN);
78        for g in corners.keys() {
79            min_q = min_q.min(g.i);
80            min_r = min_r.min(g.j);
81            max_q = max_q.max(g.i);
82            max_r = max_r.max(g.j);
83        }
84
85        let s = px_per_cell;
86        let s3h: F = sqrt3_half();
87        let half: F = lit(0.5);
88
89        let mut x_min = F::max_value().unwrap_or_else(|| lit(1e30));
90        let mut x_max = -x_min;
91        let mut y_min = x_min;
92        let mut y_max = -y_min;
93
94        for g in corners.keys() {
95            let q: F = lit(g.i as f64);
96            let r: F = lit(g.j as f64);
97            let x = s * (q + r * half);
98            let y = s * (r * s3h);
99            x_min = if x < x_min { x } else { x_min };
100            x_max = if x > x_max { x } else { x_max };
101            y_min = if y < y_min { y } else { y_min };
102            y_max = if y > y_max { y } else { y_max };
103        }
104
105        let margin_px = margin_cells * s;
106        x_min -= margin_px;
107        y_min -= margin_px;
108        x_max += margin_px;
109        y_max += margin_px;
110
111        let rect_width = nalgebra::try_convert::<F, f64>((x_max - x_min).round().max(F::one()))
112            .unwrap_or(1.0) as usize;
113        let rect_height = nalgebra::try_convert::<F, f64>((y_max - y_min).round().max(F::one()))
114            .unwrap_or(1.0) as usize;
115
116        // Build correspondences: rectified positions vs image positions
117        let mut rect_pts = Vec::with_capacity(corners.len());
118        let mut img_pts = Vec::with_capacity(corners.len());
119        for (g, &pos) in corners {
120            let q: F = lit(g.i as f64);
121            let r: F = lit(g.j as f64);
122            let rx = s * (q + r * half) - x_min;
123            let ry = s * (r * s3h) - y_min;
124            rect_pts.push(Point2::new(rx, ry));
125            img_pts.push(pos);
126        }
127
128        let h_img_from_rect =
129            estimate_homography(&rect_pts, &img_pts).ok_or(HexRectifyError::HomographyFailed)?;
130
131        let h_rect_from_img = h_img_from_rect
132            .inverse()
133            .ok_or(HexRectifyError::NonInvertible)?;
134
135        // Apply margin to axial bounds
136        let margin_ceil =
137            nalgebra::try_convert::<F, f64>(margin_cells.ceil()).unwrap_or(0.0) as i32;
138        let mq = min_q - margin_ceil;
139        let mr = min_r - margin_ceil;
140        let aq = max_q + margin_ceil;
141        let ar = max_r + margin_ceil;
142
143        Ok(Self {
144            h_img_from_rect,
145            h_rect_from_img,
146            min_q: mq,
147            min_r: mr,
148            max_q: aq,
149            max_r: ar,
150            px_per_cell,
151            rect_width,
152            rect_height,
153            x_offset: x_min,
154            y_offset: y_min,
155        })
156    }
157
158    /// Map a point from rectified space to image space.
159    pub fn rect_to_img(&self, p_rect: Point2<F>) -> Point2<F> {
160        self.h_img_from_rect.apply(p_rect)
161    }
162
163    /// Map a point from image space to rectified space.
164    pub fn img_to_rect(&self, p_img: Point2<F>) -> Point2<F> {
165        self.h_rect_from_img.apply(p_img)
166    }
167
168    /// Convert axial coordinates `(q, r)` to rectified pixel coordinates.
169    ///
170    /// This does **not** apply the homography — it maps grid indices to the
171    /// rectified coordinate system directly. The result is in the same shifted
172    /// frame used by the stored homography, so it can be passed to
173    /// [`rect_to_img`](Self::rect_to_img).
174    pub fn axial_to_rect(&self, q: F, r: F) -> Point2<F> {
175        let s = self.px_per_cell;
176        let half: F = lit(0.5);
177        let x = s * (q + r * half) - self.x_offset;
178        let y = s * (r * sqrt3_half::<F>()) - self.y_offset;
179        Point2::new(x, y)
180    }
181}
182
183#[cfg(test)]
184mod tests {
185    use super::*;
186
187    fn make_hex_corners(radius: i32, spacing: f32) -> HashMap<GridIndex, Point2<f32>> {
188        let sqrt3 = 3.0f32.sqrt();
189        let mut map = HashMap::new();
190        for q in -radius..=radius {
191            for r in -radius..=radius {
192                if (q + r).abs() > radius {
193                    continue;
194                }
195                let x = spacing * (q as f32 + r as f32 * 0.5);
196                let y = spacing * (r as f32 * sqrt3 / 2.0);
197                map.insert(GridIndex { i: q, j: r }, Point2::new(x, y));
198            }
199        }
200        map
201    }
202
203    #[test]
204    fn round_trip_rect_to_img() {
205        let corners = make_hex_corners(3, 60.0);
206        let h = HexGridHomography::from_corners(&corners, 60.0, 1.0).unwrap();
207
208        for &pos in corners.values() {
209            let rect = h.img_to_rect(pos);
210            let recovered = h.rect_to_img(rect);
211            assert!(
212                (recovered.x - pos.x).abs() < 0.5,
213                "x: {} vs {}",
214                recovered.x,
215                pos.x
216            );
217            assert!(
218                (recovered.y - pos.y).abs() < 0.5,
219                "y: {} vs {}",
220                recovered.y,
221                pos.y
222            );
223        }
224    }
225
226    #[test]
227    fn identity_case_with_ideal_positions() {
228        let corners = make_hex_corners(2, 50.0);
229        let h = HexGridHomography::from_corners(&corners, 50.0, 0.0).unwrap();
230
231        assert!(h.rect_width > 0);
232        assert!(h.rect_height > 0);
233
234        for &img_pos in corners.values() {
235            let rect_pos = h.img_to_rect(img_pos);
236            let recovered = h.rect_to_img(rect_pos);
237            assert!((recovered.x - img_pos.x).abs() < 0.1);
238            assert!((recovered.y - img_pos.y).abs() < 0.1);
239        }
240    }
241
242    #[test]
243    fn axial_to_rect_then_rect_to_img_matches_corners() {
244        let corners = make_hex_corners(3, 60.0);
245        let h = HexGridHomography::from_corners(&corners, 60.0, 1.0).unwrap();
246
247        for (g, &img_pos) in &corners {
248            let rect_pt = h.axial_to_rect(g.i as f32, g.j as f32);
249            let recovered = h.rect_to_img(rect_pt);
250            assert!(
251                (recovered.x - img_pos.x).abs() < 0.5,
252                "x mismatch at ({},{}): {} vs {}",
253                g.i,
254                g.j,
255                recovered.x,
256                img_pos.x,
257            );
258            assert!(
259                (recovered.y - img_pos.y).abs() < 0.5,
260                "y mismatch at ({},{}): {} vs {}",
261                g.i,
262                g.j,
263                recovered.y,
264                img_pos.y,
265            );
266        }
267    }
268
269    #[test]
270    fn too_few_corners_errors() {
271        let mut corners = HashMap::new();
272        corners.insert(GridIndex { i: 0, j: 0 }, Point2::new(0.0, 0.0));
273        corners.insert(GridIndex { i: 1, j: 0 }, Point2::new(50.0, 0.0));
274
275        let result = HexGridHomography::from_corners(&corners, 50.0, 0.0);
276        assert!(result.is_err());
277    }
278}