`nx_geo_clip.nx` source

1// nx_geo_clip.nx -- LIB: GEO-015 POLYGON INTERSECTION (Sutherland-Hodgman convex clip).
2//
3// THE EXCEED ANGLE (measured, not asserted): the inside/outside classification of every vertex against
4// every clip edge is the SIGN of an INTEGER cross product -- ZERO floating point, so the topology of
5// the result (which vertices survive, where edges cross) is EXACT and DETERMINISTIC. That is the part
6// float clippers (Turf/PostGIS) get wrong on touching/collinear edges. HONEST SCOPE: the intersection
7// VERTEX itself is S + (E-S)*dS/(dS-dE) computed in integer arithmetic and ROUNDED to microdegrees --
8// deterministic-integer, NOT exact-rational geometry. To stay overflow-safe at any scale the fraction
9// (num,den) is normalized (shifted down in lockstep) before the multiply, since (E-S)*dS would
10// otherwise exceed i64 at continental scale; the ratio is preserved to sub-microdegree.
11//
12// Convention (matches nx_geo): polygons = flat [lat0,lon0, ...], lat=Y, lon=X. The CLIP polygon must be
13// CONVEX and CCW. Writes the clipped polygon to out (caller sizes it >= (ns+nc+8) vertices); returns
14// the output vertex count (0 = no overlap). license_tier: ORIGINAL
15import "nx_syscalls.nx"
16const K_MAGIC_17179869184: i64 = 17179869184
17
18func geo_cabs(x: i64) -> i64 { if x < 0 { return 0 - x } return x }
19
20// is P left-of-or-on the directed clip edge A->B? (inside, for a CCW convex clip polygon). Exact.
21func geo_clip_inside(ax: i64, ay: i64, bx: i64, by: i64, px: i64, py: i64) -> i64 {
22    let cross: i64 = (bx - ax) * (py - ay) - (by - ay) * (px - ax)
23    if cross >= 0 { return 1 }
24    return 0
25}
26
27// intersection of segment S->E with the infinite line A->B, written as [lat,lon] to buf[pos]. Integer,
28// microdeg-rounded; fraction normalized to avoid i64 overflow at any coordinate scale.
29func geo_intersect_into(buf: *i64, pos: i64, sx: i64, sy: i64, ex: i64, ey: i64, ax: i64, ay: i64, bx: i64, by: i64) -> i64 {
30    let ds: i64 = (bx - ax) * (sy - ay) - (by - ay) * (sx - ax)
31    let de: i64 = (bx - ax) * (ey - ay) - (by - ay) * (ex - ax)
32    var num: i64 = ds
33    var den: i64 = ds - de
34    // keep |num|,|den| below 2^34 so (E-S)*num (|E-S| <= 3.6e8) stays under i64's 9.2e18.
35    var go: i64 = 1
36    while go == 1 {
37        if geo_cabs(num) > K_MAGIC_17179869184 { num = num / 2; den = den / 2 } else {
38            if geo_cabs(den) > K_MAGIC_17179869184 { num = num / 2; den = den / 2 } else { go = 0 }
39        }
40    }
41    if den == 0 { buf[pos * 2] = sy; buf[pos * 2 + 1] = sx; return 0 }
42    let ix: i64 = sx + (ex - sx) * num / den
43    let iy: i64 = sy + (ey - sy) * num / den
44    buf[pos * 2] = iy
45    buf[pos * 2 + 1] = ix
46    return 0
47}
48
49// clip the subject polygon by the convex CCW clip polygon (Sutherland-Hodgman). Returns vertex count.
50func geo_poly_clip(subj: *i64, ns: i64, clip: *i64, nc: i64, out: *i64) -> i64 {
51    let cap: i64 = ns + nc + 8
52    let bufA: *i64 = sys_mmap(8 * cap * 2) as *i64
53    let bufB: *i64 = sys_mmap(8 * cap * 2) as *i64
54    var cur: *i64 = bufA
55    var nxt: *i64 = bufB
56    var i: i64 = 0
57    while i < ns { cur[i * 2] = subj[i * 2]; cur[i * 2 + 1] = subj[i * 2 + 1]; i = i + 1 }
58    var curn: i64 = ns
59
60    var ci: i64 = 0
61    while ci < nc {
62        if curn == 0 { ci = nc } else {
63            let ax: i64 = clip[ci * 2 + 1]
64            let ay: i64 = clip[ci * 2]
65            let cj: i64 = (ci + 1) % nc
66            let bx: i64 = clip[cj * 2 + 1]
67            let by: i64 = clip[cj * 2]
68            var nxtn: i64 = 0
69            var k: i64 = 0
70            while k < curn {
71                let sx: i64 = cur[k * 2 + 1]
72                let sy: i64 = cur[k * 2]
73                let ek: i64 = (k + 1) % curn
74                let ex: i64 = cur[ek * 2 + 1]
75                let ey: i64 = cur[ek * 2]
76                let sIn: i64 = geo_clip_inside(ax, ay, bx, by, sx, sy)
77                let eIn: i64 = geo_clip_inside(ax, ay, bx, by, ex, ey)
78                if eIn == 1 {
79                    if sIn == 0 { geo_intersect_into(nxt, nxtn, sx, sy, ex, ey, ax, ay, bx, by); nxtn = nxtn + 1 }
80                    nxt[nxtn * 2] = ey; nxt[nxtn * 2 + 1] = ex; nxtn = nxtn + 1
81                } else {
82                    if sIn == 1 { geo_intersect_into(nxt, nxtn, sx, sy, ex, ey, ax, ay, bx, by); nxtn = nxtn + 1 }
83                }
84                k = k + 1
85            }
86            let tmp: *i64 = cur
87            cur = nxt
88            nxt = tmp
89            curn = nxtn
90            ci = ci + 1
91        }
92    }
93    i = 0
94    while i < curn { out[i * 2] = cur[i * 2]; out[i * 2 + 1] = cur[i * 2 + 1]; i = i + 1 }
95    return curn
96}
nx_geo_clip.nx source

`nx_geo_clip.nx` source
