1use crate::{Access, Coord, LineString, Provenance, SegmentDraft, Terrain, WayKind};
2use rstar::{AABB, RTree, RTreeObject};
3use serde::{Deserialize, Serialize};
4
5const METRES_PER_LATITUDE_DEGREE: f64 = 111_320.0;
6
7#[derive(Clone, Debug, PartialEq)]
10pub struct NetworkStratum {
11 pub precedence: u16,
12 pub drafts: Vec<SegmentDraft>,
13}
14
15#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
16#[serde(default)]
17pub struct ConflationPolicy {
18 pub parallel_tolerance_m: f64,
19 pub min_parallel_cosine: f64,
20 pub max_reported_decisions: usize,
21}
22
23impl Default for ConflationPolicy {
24 fn default() -> Self {
25 Self {
26 parallel_tolerance_m: 8.0,
27 min_parallel_cosine: 0.94,
28 max_reported_decisions: 2_048,
29 }
30 }
31}
32
33#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
34pub struct ConflationDecision {
35 pub preferred: Option<Provenance>,
36 pub suppressed: Option<Provenance>,
37 pub separation_m: f64,
38 pub geometry: LineString,
39}
40
41#[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
42pub struct ConflationReport {
43 pub strata: usize,
44 pub input_drafts: usize,
45 pub output_drafts: usize,
46 pub suppressed_parallel_segments: usize,
47 pub decisions: Vec<ConflationDecision>,
48}
49
50#[derive(Clone, Copy, Debug, Default, Eq, PartialEq, Serialize, Deserialize)]
51pub struct ConflationStats {
52 pub strata: usize,
53 pub input_drafts: usize,
54 pub output_drafts: usize,
55 pub suppressed_parallel_segments: usize,
56}
57
58impl From<&ConflationReport> for ConflationStats {
59 fn from(report: &ConflationReport) -> Self {
60 Self {
61 strata: report.strata,
62 input_drafts: report.input_drafts,
63 output_drafts: report.output_drafts,
64 suppressed_parallel_segments: report.suppressed_parallel_segments,
65 }
66 }
67}
68
69#[derive(Clone, Debug, PartialEq)]
70pub struct ConflatedNetwork {
71 pub drafts: Vec<SegmentDraft>,
72 pub report: ConflationReport,
73}
74
75#[derive(Clone, Copy)]
76struct IndexedPrimitive {
77 draft: usize,
78 a: Coord,
79 b: Coord,
80}
81
82impl RTreeObject for IndexedPrimitive {
83 type Envelope = AABB<[f64; 2]>;
84
85 fn envelope(&self) -> Self::Envelope {
86 AABB::from_corners(
87 [self.a.lon.min(self.b.lon), self.a.lat.min(self.b.lat)],
88 [self.a.lon.max(self.b.lon), self.a.lat.max(self.b.lat)],
89 )
90 }
91}
92
93#[derive(Clone, Copy)]
94struct ParallelMatch {
95 draft: usize,
96 separation_m: f64,
97}
98
99#[must_use]
100pub fn conflate(mut strata: Vec<NetworkStratum>, policy: ConflationPolicy) -> ConflatedNetwork {
101 strata.sort_by_key(|stratum| stratum.precedence);
102 let mut canonical = Vec::<SegmentDraft>::new();
103 let mut index = RTree::<IndexedPrimitive>::new();
104 let mut report = ConflationReport {
105 strata: strata.len(),
106 input_drafts: strata.iter().map(|stratum| stratum.drafts.len()).sum(),
107 ..ConflationReport::default()
108 };
109
110 for mut stratum in strata {
111 stratum.drafts.sort_by(draft_order);
112 let higher_precedence_count = canonical.len();
113 let mut admitted = Vec::new();
114 for draft in stratum.drafts {
115 let mut run = Vec::<Coord>::new();
116 for segment in draft.geometry.points.windows(2) {
117 let [a, b] = [segment[0], segment[1]];
118 if let Some(parallel) = nearest_parallel(&draft, a, b, &canonical, &index, policy) {
119 seal_run(&draft, &mut run, &mut admitted);
120 report.suppressed_parallel_segments += 1;
121 corroborate(&mut canonical[parallel.draft], &draft);
122 if report.decisions.len() < policy.max_reported_decisions {
123 report.decisions.push(ConflationDecision {
124 preferred: canonical[parallel.draft].provenance.first().cloned(),
125 suppressed: draft.provenance.first().cloned(),
126 separation_m: parallel.separation_m,
127 geometry: LineString::unchecked(vec![a, b]),
128 });
129 }
130 } else {
131 append_segment(&mut run, a, b);
132 }
133 }
134 seal_run(&draft, &mut run, &mut admitted);
135 }
136
137 canonical.extend(admitted);
138 for (draft, line) in canonical.iter().enumerate().skip(higher_precedence_count) {
139 for points in line.geometry.points.windows(2) {
140 index.insert(IndexedPrimitive {
141 draft,
142 a: points[0],
143 b: points[1],
144 });
145 }
146 }
147 }
148 report.output_drafts = canonical.len();
149 ConflatedNetwork {
150 drafts: canonical,
151 report,
152 }
153}
154
155fn nearest_parallel(
156 draft: &SegmentDraft,
157 a: Coord,
158 b: Coord,
159 canonical: &[SegmentDraft],
160 index: &RTree<IndexedPrimitive>,
161 policy: ConflationPolicy,
162) -> Option<ParallelMatch> {
163 if policy.parallel_tolerance_m <= 0.0 {
164 return None;
165 }
166 let midpoint = a.lerp(b, 0.5);
167 let latitude_radius = policy.parallel_tolerance_m / METRES_PER_LATITUDE_DEGREE;
168 let longitude_radius = latitude_radius / midpoint.lat.to_radians().cos().abs().max(0.05);
169 let envelope = AABB::from_corners(
170 [
171 midpoint.lon - longitude_radius,
172 midpoint.lat - latitude_radius,
173 ],
174 [
175 midpoint.lon + longitude_radius,
176 midpoint.lat + latitude_radius,
177 ],
178 );
179 index
180 .locate_in_envelope_intersecting(&envelope)
181 .filter_map(|candidate| {
182 if !same_facility(draft, &canonical[candidate.draft]) {
183 return None;
184 }
185 let (separation_m, projection) =
186 point_segment_distance(midpoint, candidate.a, candidate.b);
187 (separation_m <= policy.parallel_tolerance_m
188 && (0.0..=1.0).contains(&projection)
189 && parallel_cosine(a, b, candidate.a, candidate.b) >= policy.min_parallel_cosine)
190 .then_some(ParallelMatch {
191 draft: candidate.draft,
192 separation_m,
193 })
194 })
195 .min_by(|left, right| left.separation_m.total_cmp(&right.separation_m))
196}
197
198fn same_facility(left: &SegmentDraft, right: &SegmentDraft) -> bool {
199 if left.geometry_claim != right.geometry_claim {
200 return false;
201 }
202 let protected = |kind| {
203 matches!(
204 kind,
205 WayKind::Sidewalk
206 | WayKind::Crossing
207 | WayKind::PedestrianStreet
208 | WayKind::Roadway
209 | WayKind::ServiceRoad
210 | WayKind::Cycleway
211 | WayKind::Bushwhack
212 )
213 };
214 !(protected(left.way_kind) || protected(right.way_kind)) || left.way_kind == right.way_kind
215}
216
217fn point_segment_distance(point: Coord, a: Coord, b: Coord) -> (f64, f64) {
218 let latitude = point.lat.to_radians().cos();
219 let scale_x = METRES_PER_LATITUDE_DEGREE * latitude;
220 let scale_y = METRES_PER_LATITUDE_DEGREE;
221 let vx = (b.lon - a.lon) * scale_x;
222 let vy = (b.lat - a.lat) * scale_y;
223 let wx = (point.lon - a.lon) * scale_x;
224 let wy = (point.lat - a.lat) * scale_y;
225 let length2 = vx.mul_add(vx, vy * vy);
226 if length2 <= f64::EPSILON {
227 return (point.haversine_m(a), 0.0);
228 }
229 let projection = vx.mul_add(wx, vy * wy) / length2;
230 let t = projection.clamp(0.0, 1.0);
231 let dx = wx - t * vx;
232 let dy = wy - t * vy;
233 (dx.mul_add(dx, dy * dy).sqrt(), projection)
234}
235
236fn parallel_cosine(a: Coord, b: Coord, c: Coord, d: Coord) -> f64 {
237 let latitude = ((a.lat + b.lat + c.lat + d.lat) * 0.25).to_radians().cos();
238 let ab = ((b.lon - a.lon) * latitude, b.lat - a.lat);
239 let cd = ((d.lon - c.lon) * latitude, d.lat - c.lat);
240 let denominator = ab.0.hypot(ab.1) * cd.0.hypot(cd.1);
241 if denominator <= f64::EPSILON {
242 0.0
243 } else {
244 ab.0.mul_add(cd.0, ab.1 * cd.1).abs() / denominator
245 }
246}
247
248fn append_segment(run: &mut Vec<Coord>, a: Coord, b: Coord) {
249 if run.last().is_none_or(|last| !same_location(*last, a)) {
250 run.clear();
251 run.push(a);
252 }
253 run.push(b);
254}
255
256fn seal_run(draft: &SegmentDraft, run: &mut Vec<Coord>, admitted: &mut Vec<SegmentDraft>) {
257 if run.len() < 2 {
258 run.clear();
259 return;
260 }
261 admitted.push(draft.fragment(LineString::unchecked(std::mem::take(run))));
262}
263
264fn draft_order(left: &SegmentDraft, right: &SegmentDraft) -> std::cmp::Ordering {
265 left.provenance
266 .cmp(&right.provenance)
267 .then_with(|| {
268 left.geometry
269 .start()
270 .lon
271 .total_cmp(&right.geometry.start().lon)
272 })
273 .then_with(|| {
274 left.geometry
275 .start()
276 .lat
277 .total_cmp(&right.geometry.start().lat)
278 })
279 .then_with(|| left.geometry.end().lon.total_cmp(&right.geometry.end().lon))
280 .then_with(|| left.geometry.end().lat.total_cmp(&right.geometry.end().lat))
281}
282
283fn corroborate(preferred: &mut SegmentDraft, suppressed: &SegmentDraft) {
284 for provenance in &suppressed.provenance {
285 if !preferred.provenance.contains(provenance) {
286 preferred.provenance.push(provenance.clone());
287 }
288 }
289 preferred.confidence = preferred.confidence.max(suppressed.confidence);
290 if preferred.way_kind == WayKind::Unknown {
291 preferred.way_kind = suppressed.way_kind;
292 }
293 if preferred.standing == crate::TrailStanding::Unknown {
294 preferred.standing = suppressed.standing;
295 }
296 if preferred.marking == crate::TrailMarking::Unknown {
297 preferred.marking = suppressed.marking;
298 }
299 if preferred.terrain == Terrain::Unknown {
300 preferred.terrain = suppressed.terrain;
301 preferred.terrain_confidence = suppressed.terrain_confidence;
302 }
303 if preferred.surface.is_none() {
304 preferred.surface.clone_from(&suppressed.surface);
305 }
306 if preferred.access == Access::Unknown {
307 preferred.access = suppressed.access;
308 }
309}
310
311const fn same_location(left: Coord, right: Coord) -> bool {
312 left.lon.to_bits() == right.lon.to_bits() && left.lat.to_bits() == right.lat.to_bits()
313}
314
315#[cfg(test)]
316mod tests {
317 use super::*;
318 use crate::{
319 CrossingControl, EdgeTravel, GeometryClaim, JunctionPolicy, TrailMarking, TrailStanding,
320 WayRealm,
321 };
322
323 fn draft(name: &str, latitude: f64, standing: TrailStanding) -> SegmentDraft {
324 SegmentDraft {
325 geometry: LineString::unchecked(vec![
326 Coord::new(-74.1, latitude),
327 Coord::new(-74.0, latitude),
328 ]),
329 junctions: JunctionPolicy::Planar,
330 turn_ref: None,
331 junction_keys: None,
332 turn_restrictions: Vec::new(),
333 way_kind: WayKind::Path,
334 realm: WayRealm::default(),
335 geometry_claim: GeometryClaim::default(),
336 crossing_control: CrossingControl::default(),
337 standing,
338 marking: TrailMarking::Unknown,
339 terrain: Terrain::Trail,
340 terrain_confidence: Some(0.8),
341 surface: None,
342 access: Access::Open,
343 travel: EdgeTravel::Both,
344 road_exposure: 0.0,
345 confidence: 0.8,
346 provenance: vec![Provenance::fixture(name)],
347 }
348 }
349
350 #[test]
351 fn higher_strata_suppress_parallel_duplicates_and_keep_provenance() {
352 let mut primary = draft("osm", 41.2, TrailStanding::Informal);
353 primary.provenance[0].source = "z-authority".to_owned();
354 let mut secondary = draft("usgs", 41.200_03, TrailStanding::Established);
355 secondary.provenance[0].source = "a-corroborator".to_owned();
356 secondary.marking = TrailMarking::Marked;
357 let network = conflate(
358 vec![
359 NetworkStratum {
360 precedence: 0,
361 drafts: vec![primary],
362 },
363 NetworkStratum {
364 precedence: 10,
365 drafts: vec![secondary],
366 },
367 ],
368 ConflationPolicy::default(),
369 );
370 assert_eq!(network.drafts.len(), 1);
371 assert_eq!(network.drafts[0].standing, TrailStanding::Informal);
372 assert_eq!(network.drafts[0].marking, TrailMarking::Marked);
373 assert_eq!(network.drafts[0].provenance.len(), 2);
374 assert_eq!(network.drafts[0].provenance[0].source, "z-authority");
375 assert_eq!(network.drafts[0].provenance[1].source, "a-corroborator");
376 assert_eq!(network.report.suppressed_parallel_segments, 1);
377 }
378
379 #[test]
380 fn nonparallel_and_distant_segments_survive() {
381 let primary = draft("osm", 41.2, TrailStanding::Established);
382 let mut crossing = draft("usgs-crossing", 41.2, TrailStanding::Established);
383 crossing.geometry =
384 LineString::unchecked(vec![Coord::new(-74.05, 41.19), Coord::new(-74.05, 41.21)]);
385 let distant = draft("usgs-distant", 41.201, TrailStanding::Established);
386 let network = conflate(
387 vec![
388 NetworkStratum {
389 precedence: 0,
390 drafts: vec![primary],
391 },
392 NetworkStratum {
393 precedence: 10,
394 drafts: vec![crossing, distant],
395 },
396 ],
397 ConflationPolicy::default(),
398 );
399 assert_eq!(network.drafts.len(), 3);
400 assert_eq!(network.report.suppressed_parallel_segments, 0);
401 }
402
403 #[test]
404 fn conflation_cannot_erase_distinct_pedestrian_facilities() {
405 let mut road = draft("road", 41.2, TrailStanding::Established);
406 road.way_kind = WayKind::Roadway;
407 road.realm = WayRealm::Urban;
408 road.terrain = Terrain::Road;
409 let mut sidewalk = draft("sidewalk", 41.2, TrailStanding::Established);
410 sidewalk.way_kind = WayKind::Sidewalk;
411 sidewalk.realm = WayRealm::Urban;
412 sidewalk.terrain = Terrain::Pavement;
413
414 let network = conflate(
415 vec![
416 NetworkStratum {
417 precedence: 0,
418 drafts: vec![road],
419 },
420 NetworkStratum {
421 precedence: 10,
422 drafts: vec![sidewalk],
423 },
424 ],
425 ConflationPolicy::default(),
426 );
427
428 assert_eq!(network.drafts.len(), 2);
429 assert_eq!(network.report.suppressed_parallel_segments, 0);
430 }
431
432 #[test]
433 fn output_is_stable_under_equal_stratum_permutation() {
434 let a = draft("a", 41.2, TrailStanding::Established);
435 let b = draft("b", 41.3, TrailStanding::Established);
436 let forward = conflate(
437 vec![NetworkStratum {
438 precedence: 0,
439 drafts: vec![a.clone(), b.clone()],
440 }],
441 ConflationPolicy::default(),
442 );
443 let reverse = conflate(
444 vec![NetworkStratum {
445 precedence: 0,
446 drafts: vec![b, a],
447 }],
448 ConflationPolicy::default(),
449 );
450 let mut forward_ids = forward
451 .drafts
452 .iter()
453 .filter_map(|draft| draft.provenance[0].source_id.clone())
454 .collect::<Vec<_>>();
455 let mut reverse_ids = reverse
456 .drafts
457 .iter()
458 .filter_map(|draft| draft.provenance[0].source_id.clone())
459 .collect::<Vec<_>>();
460 forward_ids.sort();
461 reverse_ids.sort();
462 assert_eq!(forward_ids, reverse_ids);
463 }
464}