nx_geo_simplify.nx source

1// nx_geo_simplify.nx -- LIB: GEO-014 LINE SIMPLIFY (Douglas-Peucker), DETERMINISTIC INTEGER ARITH. 2// 3// THE EXCEED ANGLE (measured, not asserted): the perpendicular distance of a point to a segment is 4// computed as |cross| / isqrt(|AB|^2) using INTEGER cross products + nx_isqrt -- ZERO floating point. 5// So the keep/drop decision against the tolerance is DETERMINISTIC and reproducible to the bit, where 6// Turf.js / Mapbox simplify (float) can keep or drop a near-tolerance vertex differently across 7// platforms. (We divide rather than compare cross^2 vs eps^2*ab^2 because cross^2 overflows i64 at 8// Earth scale; |cross| <= ~2.6e17 and ab2 <= ~2.6e17 both fit, so the integer divide is overflow-safe.) 9// 10// Convention (matches nx_geo): pts = flat [lat0,lon0, ...], n points, lat=Y, lon=X. Writes the kept 11// vertex INDICES (in order, endpoints always kept) to out_idx; returns the simplified count. Iterative 12// (explicit stack, no recursion). license_tier: ORIGINAL 13import "nx_isqrt.nx" 14import "nx_syscalls.nx" 15 16func geo_iabs(x: i64) -> i64 { if x < 0 { return 0 - x } return x } 17 18// integer perpendicular distance (microdeg) of point p to the segment a->b (all indices into pts). 19func geo_perp_dist_idx(pts: *i64, a: i64, b: i64, p: i64) -> i64 { 20 let ax: i64 = pts[a * 2 + 1] 21 let ay: i64 = pts[a * 2] 22 let bx: i64 = pts[b * 2 + 1] 23 let by: i64 = pts[b * 2] 24 let px: i64 = pts[p * 2 + 1] 25 let py: i64 = pts[p * 2] 26 let dx: i64 = bx - ax 27 let dy: i64 = by - ay 28 if dx == 0 { 29 if dy == 0 { 30 // degenerate segment: distance is point-to-point. 31 let ex: i64 = px - ax 32 let ey: i64 = py - ay 33 return nx_isqrt(ex * ex + ey * ey) 34 } 35 } 36 var cross: i64 = dx * (py - ay) - dy * (px - ax) 37 cross = geo_iabs(cross) 38 let ab: i64 = nx_isqrt(dx * dx + dy * dy) 39 if ab == 0 { return 0 } 40 return cross / ab 41} 42 43// Douglas-Peucker: keep endpoints + every vertex whose perpendicular distance exceeds eps (microdeg). 44func geo_simplify(pts: *i64, n: i64, eps: i64, out_idx: *i64) -> i64 { 45 if n <= 2 { 46 var t: i64 = 0 47 while t < n { out_idx[t] = t; t = t + 1 } 48 return n 49 } 50 let keep: *i64 = sys_mmap(8 * n) as *i64 51 var i: i64 = 0 52 while i < n { keep[i] = 0; i = i + 1 } 53 keep[0] = 1 54 keep[n - 1] = 1 55 56 // explicit segment stack (lo,hi); max depth <= n. 57 let slo: *i64 = sys_mmap(8 * (n + 2)) as *i64 58 let shi: *i64 = sys_mmap(8 * (n + 2)) as *i64 59 var top: i64 = 0 60 slo[top] = 0; shi[top] = n - 1; top = top + 1 61 62 while top > 0 { 63 top = top - 1 64 let lo: i64 = slo[top] 65 let hi: i64 = shi[top] 66 var idx: i64 = 0 - 1 67 var maxd: i64 = 0 - 1 68 var k: i64 = lo + 1 69 while k < hi { 70 let d: i64 = geo_perp_dist_idx(pts, lo, hi, k) 71 if d > maxd { maxd = d; idx = k } 72 k = k + 1 73 } 74 if idx >= 0 { 75 if maxd > eps { 76 keep[idx] = 1 77 slo[top] = lo; shi[top] = idx; top = top + 1 78 slo[top] = idx; shi[top] = hi; top = top + 1 79 } 80 } 81 } 82 83 var c: i64 = 0 84 i = 0 85 while i < n { 86 if keep[i] == 1 { out_idx[c] = i; c = c + 1 } 87 i = i + 1 88 } 89 return c 90}