umsh_regiondb/
blob.rs

1//! The compiled geometry encoding.
2//!
3//! One blob holds one connected polygon component: an exterior ring and its
4//! holes, as 1e-6-degree integer coordinates, delta-encoded and zigzag-varint
5//! packed, with ring closure implicit.
6//!
7//! Working in fixed-point integers is what makes the boundary rule exact. A
8//! point on a boundary counts as inside, everywhere, and deciding that with
9//! floating-point coordinates would mean choosing an epsilon and hoping three
10//! implementations chose the same one.
11
12/// Version of the blob layout this reader understands.
13pub const GEOMETRY_FORMAT_VERSION: u8 = 1;
14
15/// Coordinate units per degree.
16pub const COORD_SCALE: f64 = 1_000_000.0;
17
18/// The role a ring plays in its polygon.
19pub const RING_EXTERIOR: u8 = 0;
20/// A hole cut out of the exterior.
21pub const RING_HOLE: u8 = 1;
22
23/// A malformed or unsupported geometry blob.
24#[derive(Clone, Debug, PartialEq, Eq)]
25pub enum GeometryError {
26    /// The blob was empty.
27    Empty,
28    /// The blob declares a version this reader does not implement.
29    UnsupportedVersion(u8),
30    /// The blob ended in the middle of a value.
31    Truncated,
32    /// A varint ran past 64 bits.
33    VarintTooLong,
34    /// Bytes remained after the last ring.
35    TrailingBytes,
36    /// The first ring was not an exterior.
37    MissingExterior,
38}
39
40impl core::fmt::Display for GeometryError {
41    fn fmt(&self, formatter: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
42        match self {
43            Self::Empty => write!(formatter, "empty geometry blob"),
44            Self::UnsupportedVersion(version) => {
45                write!(formatter, "unsupported geometry format version {version}")
46            }
47            Self::Truncated => write!(formatter, "truncated geometry blob"),
48            Self::VarintTooLong => write!(formatter, "varint too long"),
49            Self::TrailingBytes => write!(formatter, "trailing bytes after geometry blob"),
50            Self::MissingExterior => write!(formatter, "geometry part has no exterior ring"),
51        }
52    }
53}
54
55impl std::error::Error for GeometryError {}
56
57/// One closed ring, without its repeated closing vertex.
58#[derive(Clone, Debug, PartialEq, Eq)]
59pub struct Ring {
60    /// [`RING_EXTERIOR`] or [`RING_HOLE`].
61    pub role: u8,
62    /// Vertices as `(longitude, latitude)` in 1e-6 degrees.
63    pub points: Vec<(i32, i32)>,
64}
65
66/// Quantize a degree value onto the storage grid, rounding half away from zero.
67pub fn to_e6(degrees: f64) -> i32 {
68    let scaled = degrees * COORD_SCALE;
69    let rounded = if scaled >= 0.0 {
70        (scaled + 0.5).floor()
71    } else {
72        -((-scaled + 0.5).floor())
73    };
74    rounded.clamp(i32::MIN as f64, i32::MAX as f64) as i32
75}
76
77/// Convert a stored coordinate back to degrees.
78pub fn from_e6(value: i32) -> f64 {
79    f64::from(value) / COORD_SCALE
80}
81
82struct Reader<'a> {
83    data: &'a [u8],
84    offset: usize,
85}
86
87impl<'a> Reader<'a> {
88    fn byte(&mut self) -> Result<u8, GeometryError> {
89        let byte = *self.data.get(self.offset).ok_or(GeometryError::Truncated)?;
90        self.offset += 1;
91        Ok(byte)
92    }
93
94    fn varint(&mut self) -> Result<u64, GeometryError> {
95        let mut result: u64 = 0;
96        let mut shift = 0u32;
97        loop {
98            let byte = self.byte()?;
99            result |= u64::from(byte & 0x7F) << shift;
100            if byte & 0x80 == 0 {
101                return Ok(result);
102            }
103            shift += 7;
104            if shift > 63 {
105                return Err(GeometryError::VarintTooLong);
106            }
107        }
108    }
109
110    fn zigzag(&mut self) -> Result<i64, GeometryError> {
111        let raw = self.varint()?;
112        Ok(((raw >> 1) as i64) ^ -((raw & 1) as i64))
113    }
114}
115
116/// Decode one polygon component.
117pub fn decode(data: &[u8]) -> Result<Vec<Ring>, GeometryError> {
118    if data.is_empty() {
119        return Err(GeometryError::Empty);
120    }
121    let version = data[0];
122    if version != GEOMETRY_FORMAT_VERSION {
123        return Err(GeometryError::UnsupportedVersion(version));
124    }
125
126    let mut reader = Reader { data, offset: 1 };
127    let ring_count = reader.varint()?;
128    let mut rings = Vec::with_capacity(ring_count as usize);
129    for index in 0..ring_count {
130        let role = reader.byte()?;
131        if index == 0 && role != RING_EXTERIOR {
132            return Err(GeometryError::MissingExterior);
133        }
134        let point_count = reader.varint()?;
135        let mut points = Vec::with_capacity(point_count as usize);
136        let mut longitude: i64 = 0;
137        let mut latitude: i64 = 0;
138        for _ in 0..point_count {
139            longitude += reader.zigzag()?;
140            latitude += reader.zigzag()?;
141            points.push((longitude as i32, latitude as i32));
142        }
143        rings.push(Ring { role, points });
144    }
145    if reader.offset != data.len() {
146        return Err(GeometryError::TrailingBytes);
147    }
148    Ok(rings)
149}
150
151/// Exact integer point-in-polygon, boundary inclusive.
152///
153/// A point exactly on a boundary is inside. Two abutting regions therefore both
154/// claim their shared edge, which costs an operator one extra region in a list
155/// they are reviewing anyway; a gap between them would leave a position with no
156/// region at all.
157pub fn point_in_rings(rings: &[Ring], longitude_e6: i32, latitude_e6: i32) -> bool {
158    let mut exterior = None;
159    for ring in rings {
160        if on_ring(ring, longitude_e6, latitude_e6) {
161            return true;
162        }
163        if ring.role == RING_EXTERIOR && exterior.is_none() {
164            exterior = Some(ring);
165        }
166    }
167
168    let Some(exterior) = exterior else {
169        return false;
170    };
171    if !strictly_inside(exterior, longitude_e6, latitude_e6) {
172        return false;
173    }
174    !rings
175        .iter()
176        .filter(|ring| ring.role == RING_HOLE)
177        .any(|hole| strictly_inside(hole, longitude_e6, latitude_e6))
178}
179
180fn on_ring(ring: &Ring, x: i32, y: i32) -> bool {
181    let points = &ring.points;
182    let x = i64::from(x);
183    let y = i64::from(y);
184    for index in 0..points.len() {
185        let (x1, y1) = points[index];
186        let (x2, y2) = points[(index + 1) % points.len()];
187        let (x1, y1, x2, y2) = (i64::from(x1), i64::from(y1), i64::from(x2), i64::from(y2));
188        if (x2 - x1) * (y - y1) - (y2 - y1) * (x - x1) != 0 {
189            continue;
190        }
191        if x1.min(x2) <= x && x <= x1.max(x2) && y1.min(y2) <= y && y <= y1.max(y2) {
192            return true;
193        }
194    }
195    false
196}
197
198/// Crossing-number test with a half-open rule on each edge's vertical span, so
199/// that a ray through a vertex counts once rather than twice or not at all.
200fn strictly_inside(ring: &Ring, x: i32, y: i32) -> bool {
201    let points = &ring.points;
202    let x = i64::from(x);
203    let y = i64::from(y);
204    let mut inside = false;
205    for index in 0..points.len() {
206        let (x1, y1) = points[index];
207        let (x2, y2) = points[(index + 1) % points.len()];
208        let (x1, y1, x2, y2) = (i64::from(x1), i64::from(y1), i64::from(x2), i64::from(y2));
209        if (y1 > y) != (y2 > y) {
210            let side = (x2 - x1) * (y - y1) - (y2 - y1) * (x - x1);
211            if side != 0 && (side > 0) == (y2 > y1) {
212                inside = !inside;
213            }
214        }
215    }
216    inside
217}