code wiki / _hdl_build / nx_geo_hull.nx
nx_geo_hull.nx source
↩ module page · 102 lines · 3476 B
1// nx_geo_hull.nx -- LIB: GEO-012 CONVEX HULL (Andrew's monotone chain), INTEGER-EXACT.
2//
3// THE EXCEED ANGLE (measured, not asserted): the turn-direction test is a single INTEGER cross
4// product over microdegree coordinates -- ZERO floating point. So the orientation predicate is EXACT,
5// and collinear / duplicate points are classified deterministically (cross == 0 is decided, never a
6// float epsilon coin-flip). Those are exactly the inputs where Turf.js / Shapely / qhull (float) emit
7// non-deterministic hulls or spurious near-collinear vertices. We drop collinear edge points (pop on
8// cross <= 0) for a minimal, exact hull.
9//
10// Convention (matches nx_geo): pts = flat [lat0,lon0, ...], n points, lat = Y, lon = X. Writes the
11// hull vertex INDICES (into pts) in counter-clockwise order to out_idx; returns the hull size.
12// Overflow-safe for Earth coords in i64 (cross of two deltas <= 3.6e8^2 ~ 1.3e17, no accumulation).
13// Selection-sort O(n^2) for now (exactness first; an O(n log n) sort is a later perf rung). license_tier: ORIGINAL
14import "nx_syscalls.nx"
15
16// is point i < point j ordered by (lon asc, then lat asc)? -> 1 yes, 0 no.
17func geo_cmp_lonlat(pts: *i64, i: i64, j: i64) -> i64 {
18 let xi: i64 = pts[i * 2 + 1]
19 let xj: i64 = pts[j * 2 + 1]
20 if xi < xj { return 1 }
21 if xi > xj { return 0 }
22 let yi: i64 = pts[i * 2]
23 let yj: i64 = pts[j * 2]
24 if yi < yj { return 1 }
25 return 0
26}
27
28// 2D cross product (a-o) x (b-o), x=lon, y=lat. > 0 left turn (CCW), < 0 right, == 0 collinear. Exact.
29func geo_cross3(pts: *i64, o: i64, a: i64, b: i64) -> i64 {
30 let ox: i64 = pts[o * 2 + 1]
31 let oy: i64 = pts[o * 2]
32 let ax: i64 = pts[a * 2 + 1]
33 let ay: i64 = pts[a * 2]
34 let bx: i64 = pts[b * 2 + 1]
35 let by: i64 = pts[b * 2]
36 return (ax - ox) * (by - oy) - (ay - oy) * (bx - ox)
37}
38
39func geo_convex_hull(pts: *i64, n: i64, out_idx: *i64) -> i64 {
40 if n < 3 {
41 var t: i64 = 0
42 while t < n { out_idx[t] = t; t = t + 1 }
43 return n
44 }
45
46 // sorted index array (selection sort by lon,lat).
47 let srt: *i64 = sys_mmap(8 * n) as *i64
48 var i: i64 = 0
49 while i < n { srt[i] = i; i = i + 1 }
50 var a: i64 = 0
51 while a < n {
52 var mn: i64 = a
53 var b: i64 = a + 1
54 while b < n {
55 if geo_cmp_lonlat(pts, srt[b], srt[mn]) == 1 { mn = b }
56 b = b + 1
57 }
58 let tmp: i64 = srt[a]; srt[a] = srt[mn]; srt[mn] = tmp
59 a = a + 1
60 }
61
62 // monotone chain into H (size 2n). H stores ORIGINAL point indices.
63 let H: *i64 = sys_mmap(8 * (2 * n)) as *i64
64 var k: i64 = 0
65
66 // lower hull
67 i = 0
68 while i < n {
69 let pi: i64 = srt[i]
70 var go: i64 = 1
71 while go == 1 {
72 if k < 2 { go = 0 } else {
73 if geo_cross3(pts, H[k - 2], H[k - 1], pi) <= 0 { k = k - 1 } else { go = 0 }
74 }
75 }
76 H[k] = pi
77 k = k + 1
78 i = i + 1
79 }
80
81 // upper hull
82 let lo_t: i64 = k + 1
83 i = n - 2
84 while i >= 0 {
85 let pi: i64 = srt[i]
86 var go: i64 = 1
87 while go == 1 {
88 if k < lo_t { go = 0 } else {
89 if geo_cross3(pts, H[k - 2], H[k - 1], pi) <= 0 { k = k - 1 } else { go = 0 }
90 }
91 }
92 H[k] = pi
93 k = k + 1
94 i = i - 1
95 }
96
97 // last point == first; drop it.
98 let hsz: i64 = k - 1
99 i = 0
100 while i < hsz { out_idx[i] = H[i]; i = i + 1 }
101 return hsz
102}