nx_geo_union.nx source

1// nx_geo_union.nx -- LIB: GEO-016 POLYGON UNION + INTERSECTION MEASURE (set algebra on areas). 2// 3// The union of two convex polygons is generally NON-convex, so it is not produced by a single 4// Sutherland-Hodgman clip. Instead we compute the exact set MEASURE by inclusion-exclusion: 5// area(A u B) = area(A) + area(B) - area(A n B) 6// composing GEO-010 (geo_area2, integer-exact shoelace) and GEO-015 (geo_poly_clip, Sutherland-Hodgman 7// intersection). HONEST SCOPE: this returns the union/intersection COVERAGE MEASURE (doubled area in 8// microdeg^2), not the union BOUNDARY polygon -- which is exactly what geofence coverage / overlap 9// scoring needs ("how much of region A is also covered by region B?"). 10// 11// EXCEED ANGLE: every term is integer arithmetic over microdegree coords -> deterministic, no float 12// drift in the inclusion-exclusion (where chained float set-ops accumulate error). geo_iou_permil adds 13// a sovereign Jaccard overlap score (intersection/union, per-mille) for geofence dedup/matching. 14// Convention: polygons = flat [lat0,lon0,...]; the SECOND polygon (b) is the clip and must be CONVEX+CCW. 15// license_tier: ORIGINAL 16import "nx_geo_clip.nx" 17import "nx_geo_area.nx" 18import "nx_syscalls.nx" 19 20// doubled area of A n B (microdeg^2); 0 if they do not overlap. 21func geo_inter_area2(a: *i64, na: i64, b: *i64, nb: i64) -> i64 { 22 let cap: i64 = na + nb + 8 23 let inter: *i64 = sys_mmap(8 * cap * 2) as *i64 24 let ni: i64 = geo_poly_clip(a, na, b, nb, inter) 25 if ni < 3 { return 0 } 26 return geo_area2(inter, ni) 27} 28 29// doubled area of A u B (microdeg^2) via inclusion-exclusion. 30func geo_union_area2(a: *i64, na: i64, b: *i64, nb: i64) -> i64 { 31 let aa: i64 = geo_area2(a, na) 32 let bb: i64 = geo_area2(b, nb) 33 let ai: i64 = geo_inter_area2(a, na, b, nb) 34 return aa + bb - ai 35} 36 37// Jaccard overlap (intersection / union) in per-mille [0,1000]; 0 if union is empty. The sovereign 38// geofence-similarity / dedup score. 39func geo_iou_permil(a: *i64, na: i64, b: *i64, nb: i64) -> i64 { 40 let inter: i64 = geo_inter_area2(a, na, b, nb) 41 let aa: i64 = geo_area2(a, na) 42 let bb: i64 = geo_area2(b, nb) 43 let uni: i64 = aa + bb - inter 44 if uni <= 0 { return 0 } 45 return inter * 1000 / uni 46}