code wiki / _hdl_build / nx_geo_segint.nx
nx_geo_segint.nx source
↩ module page · 62 lines · 2604 B
1// nx_geo_segint.nx -- LIB: GEO-015b SEGMENT-SEGMENT INTERSECTION primitive (the building block for
2// non-convex polygon clipping / Greiner-Hormann, hardening GEO-015/016 toward real boundary polygons).
3//
4// Parametric: P1 + t*(P2-P1) = P3 + u*(P4-P3). With d1=P2-P1, d2=P4-P3, e=P3-P1:
5// denom = cross(d1,d2); t = cross(e,d2)/denom; u = cross(e,d1)/denom.
6// All INTEGER cross products -> the topological classification is EXACT and deterministic:
7// 0 = no intersection (parallel-disjoint, or lines meet outside one of the segments)
8// 1 = proper interior crossing (single point, strictly inside both)
9// 2 = touch at an endpoint (t or u exactly 0 or 1)
10// 3 = collinear (overlapping supports)
11// The intersection POINT (codes 1,2) is written to out2 = [x,y], microdeg-rounded, with the t fraction
12// normalized (lockstep >>1) so P1 + t*d1 can't overflow i64 at any scale (same trick as GEO-015).
13// Generic 2D (caller maps lat/lon as it likes). license_tier: ORIGINAL
14import "nx_syscalls.nx"
15const K_MAGIC_17179869184: i64 = 17179869184
16
17func geo_si_abs(x: i64) -> i64 { if x < 0 { return 0 - x } return x }
18
19func geo_seg_intersect(p1x: i64, p1y: i64, p2x: i64, p2y: i64, p3x: i64, p3y: i64, p4x: i64, p4y: i64, out2: *i64) -> i64 {
20 let d1x: i64 = p2x - p1x
21 let d1y: i64 = p2y - p1y
22 let d2x: i64 = p4x - p3x
23 let d2y: i64 = p4y - p3y
24 let ex: i64 = p3x - p1x
25 let ey: i64 = p3y - p1y
26
27 let denom: i64 = d1x * d2y - d1y * d2x
28 if denom == 0 {
29 // parallel; collinear iff e is parallel to d1 too.
30 if (ex * d1y - ey * d1x) == 0 { return 3 }
31 return 0
32 }
33
34 var tn: i64 = ex * d2y - ey * d2x
35 var un: i64 = ex * d1y - ey * d1x
36 var dn: i64 = denom
37 if dn < 0 { dn = 0 - dn; tn = 0 - tn; un = 0 - un } // normalize sign so checks use dn>0
38
39 if tn < 0 { return 0 }
40 if tn > dn { return 0 }
41 if un < 0 { return 0 }
42 if un > dn { return 0 }
43
44 // intersection point = P1 + (tn/dn)*d1 ; normalize the fraction to avoid overflow on tn*d1.
45 var n2: i64 = tn
46 var q2: i64 = dn
47 var go: i64 = 1
48 while go == 1 {
49 if geo_si_abs(n2) > K_MAGIC_17179869184 { n2 = n2 / 2; q2 = q2 / 2 } else {
50 if geo_si_abs(q2) > K_MAGIC_17179869184 { n2 = n2 / 2; q2 = q2 / 2 } else { go = 0 }
51 }
52 }
53 out2[0] = p1x + n2 * d1x / q2
54 out2[1] = p1y + n2 * d1y / q2
55
56 // classify: endpoint-touch if a parameter is exactly at a bound, else proper interior crossing.
57 if tn == 0 { return 2 }
58 if tn == dn { return 2 }
59 if un == 0 { return 2 }
60 if un == dn { return 2 }
61 return 1
62}