nx_geo_centroid.nx source

1// nx_geo_centroid.nx -- LIB: GEO-011 POLYGON CENTROID, two integer-exact variants. 2// 3// geo_centroid_vertex -- mean of the vertices (= turf.centroid). Exact at ANY scale; but sensitive 4// to how densely each edge is sampled with vertices. 5// geo_centroid_area -- the AREA-WEIGHTED centroid (= PostGIS ST_Centroid / turf.centerOfMass): 6// Cx = (1/(3*A2)) * sum (x_i+x_j)*cross, Cy likewise, cross = x_j*y_i-x_i*y_j. 7// 8// THE EXCEED ANGLE (measured): the area-weighted centroid is ROBUST -- it is unchanged by redundant 9// collinear edge vertices, whereas the vertex-mean shifts toward densely-sampled edges. The gate proves 10// they DIFFER on exactly such a polygon, and that the area-weighted one stays at the true center. Both 11// integer-exact (no float drift). 12// 13// OVERFLOW: geo_centroid_area origin-shifts to vertex 0 so the sums stay in delta space -- i64-exact 14// for regional polygons (sub-degree extent, modest vertex count = the geofence domain). Continental- 15// scale area-weighted centroids need the 128-bit path (nx_mul_wide); a documented future hardening. 16// Convention: poly = flat [lat0,lon0, ...], lat=Y, lon=X. license_tier: ORIGINAL 17import "nx_syscalls.nx" 18 19// mean of vertices -> out2 = [lat,lon]. 20func geo_centroid_vertex(poly: *i64, n: i64, out2: *i64) -> i64 { 21 var sy: i64 = 0 22 var sx: i64 = 0 23 var i: i64 = 0 24 while i < n { 25 sy = sy + poly[i * 2] 26 sx = sx + poly[i * 2 + 1] 27 i = i + 1 28 } 29 out2[0] = sy / n 30 out2[1] = sx / n 31 return 0 32} 33 34// area-weighted centroid -> out2 = [lat,lon]. Origin-shifted to vertex 0 for i64 headroom. 35func geo_centroid_area(poly: *i64, n: i64, out2: *i64) -> i64 { 36 let oy: i64 = poly[0] 37 let ox: i64 = poly[1] 38 var a2: i64 = 0 39 var sx: i64 = 0 40 var sy: i64 = 0 41 var i: i64 = 0 42 var j: i64 = n - 1 43 while i < n { 44 let yi: i64 = poly[i * 2] - oy 45 let xi: i64 = poly[i * 2 + 1] - ox 46 let yj: i64 = poly[j * 2] - oy 47 let xj: i64 = poly[j * 2 + 1] - ox 48 let cross: i64 = xj * yi - xi * yj 49 a2 = a2 + cross 50 sx = sx + (xj + xi) * cross 51 sy = sy + (yj + yi) * cross 52 j = i 53 i = i + 1 54 } 55 if a2 == 0 { out2[0] = oy; out2[1] = ox; return 0 } // degenerate (zero-area) -> vertex 0 56 out2[0] = oy + sy / (3 * a2) 57 out2[1] = ox + sx / (3 * a2) 58 return 0 59}