nx_polygon.nx source
↩ module page · 512 lines · 19105 B
1// nx_polygon.nx -- 2D polygon primitive with foundational ops +
2// constant-distance Minkowski offset. Foundational substrate for
3// the slicer pipeline (per NISHI_3D_PRINT_ROADMAP §2-3): every
4// slice contour, every perimeter loop, every infill region is an
5// NxPolygon manipulated through this primitive.
6//
7// Distinct from:
8// nx_poly.nx -- polynomial algebra (p[0] + p[1]*x + ...)
9// nx_geom.nx -- Hough transform (vision-pipeline geometry)
10// nx_polyline -- (does not exist; if we add open polylines they
11// go in a separate primitive)
12//
13// Coordinate convention: Q14 fixed-point xy. 1 Q14 unit = 1/16384.
14// Matches nx_mesh + nx_stl pipeline so slice contours derived from
15// triangle-plane intersections compose without unit conversion.
16//
17// Winding convention: CCW = solid (positive signed area), CW = hole
18// (negative signed area). Same convention as nx_mesh (face front),
19// most CAD systems, and the SVG / PDF / PostScript "fill rule".
20//
21// Per cardinal feedback-bits-up-exceed-never-match: composes the
22// substrate's nx_isqrt (Q14 sqrt for unit-normal computation) +
23// nx_i128 (overflow protection in line-line intersection). No
24// CGAL, no libigl, no Boost.Geometry.
25//
26// Numerical robustness:
27// - Foundational ops (area, point-in-poly, segment-intersect) use
28// pure i64 arithmetic with cross-multiplication to avoid divides.
29// - Coordinate magnitudes ≤ 325 mm * 16384 ≈ 2^23 (Qidi X-Max 3
30// build volume). Products of two coord-deltas ≤ 2^46. Sums
31// of cross-products ≤ ~2^47. All fit in i64.
32// - Offset uses nx_i128 for the t_num * edge_delta product in
33// line-line intersection (would be 2^70 in i64, overflows).
34//
35// Per cardinal feedback-real-playtest-loop-not-just-logs: smoke
36// passes are necessary but not sufficient. Offset on degenerate
37// polygons (180-degree vertices, near-180 acute corners, near-
38// collinear edges) needs print-evidence validation before claiming
39// S-class.
40//
41// license_tier: ORIGINAL
42
43import "nx_syscalls.nx"
44import "nx_isqrt.nx"
45import "nx_i128.nx"
46
47const NX_POLYGON_VERT_STRIDE: i64 = 16 // 2 × i64 per vertex (x, y)
48const NX_POLYGON_Q14_ONE: i64 = 16384
49
50// ===== sealed-enum verdicts ========================================
51
52const NX_POLYGON_OK: i64 = 0
53const NX_POLYGON_ERR_TOO_FEW_VERTS: i64 = 1
54const NX_POLYGON_ERR_DEGENERATE_EDGE: i64 = 2
55const NX_POLYGON_ERR_OFFSET_COLLAPSED: i64 = 3
56const NX_POLYGON_ERR_CAPACITY: i64 = 4
57const NX_POLYGON_N_VERDICTS: i64 = 5
58
59func nx_polygon_verdict_is_valid(v: i64) -> i64 {
60 if v < 0 { return 0 }
61 if v >= NX_POLYGON_N_VERDICTS { return 0 }
62 return 1
63}
64
65func nx_polygon_verdict_name(v: i64) -> *u8 {
66 if v == NX_POLYGON_OK { return "OK" }
67 if v == NX_POLYGON_ERR_TOO_FEW_VERTS { return "TOO_FEW_VERTS" }
68 if v == NX_POLYGON_ERR_DEGENERATE_EDGE { return "DEGENERATE_EDGE" }
69 if v == NX_POLYGON_ERR_OFFSET_COLLAPSED { return "OFFSET_COLLAPSED" }
70 if v == NX_POLYGON_ERR_CAPACITY { return "CAPACITY" }
71 return "UNKNOWN"
72}
73
74// ===== struct ======================================================
75
76struct NxPolygon {
77 verts: *u8, // packed Q14 (x, y) pairs
78 n_verts: i64,
79 capacity: i64,
80}
81
82const NX_POLYGON_BYTES: i64 = 24
83
84// ===== alloc + setters + getters ===================================
85
86func nx_polygon_alloc(capacity: i64) -> *NxPolygon {
87 if capacity <= 0 { return 0 as *NxPolygon }
88 let p: *NxPolygon = (sys_mmap(NX_POLYGON_BYTES)) as *NxPolygon
89 p.verts = sys_mmap(capacity * NX_POLYGON_VERT_STRIDE)
90 p.n_verts = 0
91 p.capacity = capacity
92 return p
93}
94
95func nx_polygon_vert_ptr(p: *NxPolygon, i: i64) -> *i64 {
96 return ((p.verts as i64) + i * NX_POLYGON_VERT_STRIDE) as *i64
97}
98
99func nx_polygon_get_x(p: *NxPolygon, i: i64) -> i64 {
100 let v: *i64 = nx_polygon_vert_ptr(p, i)
101 return v[0]
102}
103
104func nx_polygon_get_y(p: *NxPolygon, i: i64) -> i64 {
105 let v: *i64 = nx_polygon_vert_ptr(p, i)
106 return v[1]
107}
108
109func nx_polygon_set_vert(p: *NxPolygon, i: i64, x: i64, y: i64) -> i64 {
110 if i < 0 { return NX_POLYGON_ERR_CAPACITY }
111 if i >= p.capacity { return NX_POLYGON_ERR_CAPACITY }
112 let v: *i64 = nx_polygon_vert_ptr(p, i)
113 v[0] = x
114 v[1] = y
115 return NX_POLYGON_OK
116}
117
118func nx_polygon_add_vert(p: *NxPolygon, x: i64, y: i64) -> i64 {
119 if p.n_verts >= p.capacity { return -1 }
120 let idx: i64 = p.n_verts
121 nx_polygon_set_vert(p, idx, x, y)
122 p.n_verts = idx + 1
123 return idx
124}
125
126// ===== signed area (Shoelace) =====================================
127//
128// Returns 2 * signed area, in Q28 units (product of two Q14 coords).
129// Sign: positive for CCW (solid), negative for CW (hole). Caller
130// halves + divides by 16384*16384 to get real mm² if needed.
131//
132// Uses nx_i128 accumulator to avoid overflow on large polygons.
133
134func nx_polygon_signed_2area_q28_hi(p: *NxPolygon) -> i64 {
135 // For polygons of slicer-typical size (≤ a few thousand verts on
136 // a 325mm-class print), the sum fits in i64. Returns the high
137 // (overflow) word from i128 accumulator -- callers can check
138 // this is zero (or sign-extended) for safe i64 fallback.
139 if p.n_verts < 3 { return 0 }
140 let acc: *I128 = nx_i128_alloc()
141 nx_i128_set_i64(acc, 0)
142 let term: *I128 = nx_i128_alloc()
143 var i: i64 = 0
144 while i < p.n_verts {
145 let j: i64 = (i + 1)
146 var jj: i64 = j
147 if jj >= p.n_verts { jj = 0 }
148 let xi: i64 = nx_polygon_get_x(p, i)
149 let yi: i64 = nx_polygon_get_y(p, i)
150 let xj: i64 = nx_polygon_get_x(p, jj)
151 let yj: i64 = nx_polygon_get_y(p, jj)
152 // term1 = xi * yj
153 nx_i128_mul_i64(xi, yj, term)
154 nx_i128_add(acc, term)
155 // term2 = -(xj * yi)
156 nx_i128_mul_i64(xj, yi, term)
157 nx_i128_neg(term)
158 nx_i128_add(acc, term)
159 i = i + 1
160 }
161 return acc.hi
162}
163
164func nx_polygon_signed_2area_q28(p: *NxPolygon) -> i64 {
165 // i64 fast path for the common case (polygon size where the
166 // shoelace sum doesn't overflow). Slicer-typical contours
167 // fit comfortably.
168 if p.n_verts < 3 { return 0 }
169 var acc: i64 = 0
170 var i: i64 = 0
171 while i < p.n_verts {
172 var jj: i64 = i + 1
173 if jj >= p.n_verts { jj = 0 }
174 let xi: i64 = nx_polygon_get_x(p, i)
175 let yi: i64 = nx_polygon_get_y(p, i)
176 let xj: i64 = nx_polygon_get_x(p, jj)
177 let yj: i64 = nx_polygon_get_y(p, jj)
178 acc = acc + xi * yj - xj * yi
179 i = i + 1
180 }
181 return acc
182}
183
184func nx_polygon_is_ccw(p: *NxPolygon) -> i64 {
185 if nx_polygon_signed_2area_q28(p) > 0 { return 1 }
186 return 0
187}
188
189// ===== bbox ========================================================
190
191func nx_polygon_bbox(p: *NxPolygon, out: *i64) -> i64 {
192 if p.n_verts <= 0 { return NX_POLYGON_ERR_TOO_FEW_VERTS }
193 var mnx: i64 = nx_polygon_get_x(p, 0)
194 var mny: i64 = nx_polygon_get_y(p, 0)
195 var mxx: i64 = mnx
196 var mxy: i64 = mny
197 var i: i64 = 1
198 while i < p.n_verts {
199 let x: i64 = nx_polygon_get_x(p, i)
200 let y: i64 = nx_polygon_get_y(p, i)
201 if x < mnx { mnx = x }
202 if y < mny { mny = y }
203 if x > mxx { mxx = x }
204 if y > mxy { mxy = y }
205 i = i + 1
206 }
207 out[0] = mnx; out[1] = mny; out[2] = mxx; out[3] = mxy
208 return NX_POLYGON_OK
209}
210
211// ===== point-in-polygon (ray-cast crossing) =======================
212//
213// Cast a horizontal ray from (x, y) in the +x direction; count
214// edges that strictly straddle y. Odd count = inside.
215//
216// Uses cross-multiplication to avoid division. All comparisons
217// exact in i64.
218
219func nx_polygon_point_inside(p: *NxPolygon, x: i64, y: i64) -> i64 {
220 if p.n_verts < 3 { return 0 }
221 var crossings: i64 = 0
222 var i: i64 = 0
223 while i < p.n_verts {
224 var jj: i64 = i + 1
225 if jj >= p.n_verts { jj = 0 }
226 let xi: i64 = nx_polygon_get_x(p, i)
227 let yi: i64 = nx_polygon_get_y(p, i)
228 let xj: i64 = nx_polygon_get_x(p, jj)
229 let yj: i64 = nx_polygon_get_y(p, jj)
230 // Edge straddles ray if (yi > y) != (yj > y).
231 var yi_above: i64 = 0
232 if yi > y { yi_above = 1 }
233 var yj_above: i64 = 0
234 if yj > y { yj_above = 1 }
235 if yi_above != yj_above {
236 // Crossing x-coord = xi + (y - yi) * (xj - xi) / (yj - yi).
237 // We want: x < crossing. Multiply through by (yj - yi),
238 // flipping the inequality if (yj - yi) < 0.
239 let dy: i64 = yj - yi
240 let dx: i64 = xj - xi
241 let lhs: i64 = (x - xi) * dy
242 let rhs: i64 = (y - yi) * dx
243 if dy > 0 {
244 if lhs < rhs { crossings = crossings + 1 }
245 }
246 if dy < 0 {
247 if lhs > rhs { crossings = crossings + 1 }
248 }
249 }
250 i = i + 1
251 }
252 return crossings & 1
253}
254
255// ===== orientation predicate =======================================
256//
257// orient(a, b, c) = sign of cross((b - a), (c - a)).
258// Returns +1 (CCW), -1 (CW), 0 (collinear). Pure i64.
259
260func nx_polygon_orient(ax: i64, ay: i64, bx: i64, by: i64,
261 cx: i64, cy: i64) -> i64 {
262 let cross: i64 = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax)
263 if cross > 0 { return 1 }
264 if cross < 0 { return -1 }
265 return 0
266}
267
268// ===== segment-segment proper intersection ========================
269//
270// Returns 1 iff segments (a,b) and (c,d) cross in their interiors.
271// Collinear-overlapping case returns 0 (slicer treats this as no
272// intersection; ambiguous geometry is the caller's problem).
273
274func nx_polygon_segments_cross(ax: i64, ay: i64, bx: i64, by: i64,
275 cx: i64, cy: i64, dx: i64, dy: i64) -> i64 {
276 let o1: i64 = nx_polygon_orient(ax, ay, bx, by, cx, cy)
277 let o2: i64 = nx_polygon_orient(ax, ay, bx, by, dx, dy)
278 let o3: i64 = nx_polygon_orient(cx, cy, dx, dy, ax, ay)
279 let o4: i64 = nx_polygon_orient(cx, cy, dx, dy, bx, by)
280 if o1 == 0 { return 0 }
281 if o2 == 0 { return 0 }
282 if o3 == 0 { return 0 }
283 if o4 == 0 { return 0 }
284 if o1 == o2 { return 0 }
285 if o3 == o4 { return 0 }
286 return 1
287}
288
289// ===== constant-distance offset ===================================
290//
291// Minkowski offset of `in_poly` by signed distance `d_q14`.
292// Positive d for CCW polygons = outward (polygon grows).
293// Negative d for CCW polygons = inward (polygon shrinks).
294//
295// Algorithm: for each vertex, compute the line-line intersection
296// of the two adjacent offset edges. Each edge is shifted by its
297// outward unit normal scaled by d. Sharp convex corners can
298// stretch out (no bevel / mitre limit for v1; queued for future
299// "sharp corner clamp" enhancement).
300//
301// Composition: nx_isqrt for unit-normal length, nx_i128 for the
302// t_num * edge_delta intermediate product (would overflow i64 for
303// large polygons).
304
305func nx_polygon_offset(in_poly: *NxPolygon, d_q14: i64) -> *NxPolygon {
306 if in_poly.n_verts < 3 { return 0 as *NxPolygon }
307 let n: i64 = in_poly.n_verts
308 let out: *NxPolygon = nx_polygon_alloc(n)
309 if (out as i64) == 0 { return out }
310 out.n_verts = n
311
312 var i: i64 = 0
313 while i < n {
314 var i_prev: i64 = i - 1
315 if i_prev < 0 { i_prev = n - 1 }
316 var i_next: i64 = i + 1
317 if i_next >= n { i_next = 0 }
318
319 let px: i64 = nx_polygon_get_x(in_poly, i_prev)
320 let py: i64 = nx_polygon_get_y(in_poly, i_prev)
321 let cx: i64 = nx_polygon_get_x(in_poly, i)
322 let cy: i64 = nx_polygon_get_y(in_poly, i)
323 let nx_: i64 = nx_polygon_get_x(in_poly, i_next)
324 let ny_: i64 = nx_polygon_get_y(in_poly, i_next)
325
326 // Edge vectors
327 let e_prev_dx: i64 = cx - px
328 let e_prev_dy: i64 = cy - py
329 let e_next_dx: i64 = nx_ - cx
330 let e_next_dy: i64 = ny_ - cy
331
332 // Outward unit normal of each edge in Q14 (CCW polygon ->
333 // rotate edge direction 90° clockwise = (edge.y, -edge.x)).
334 // length_sq is in Q28, nx_isqrt returns Q14 length.
335 let lsq_prev: i64 = e_prev_dx * e_prev_dx + e_prev_dy * e_prev_dy
336 let lsq_next: i64 = e_next_dx * e_next_dx + e_next_dy * e_next_dy
337 let len_prev_q14: i64 = nx_isqrt(lsq_prev)
338 let len_next_q14: i64 = nx_isqrt(lsq_next)
339
340 if len_prev_q14 == 0 { return 0 as *NxPolygon }
341 if len_next_q14 == 0 { return 0 as *NxPolygon }
342
343 // Q14 unit normals = edge_perp * Q14_ONE / length
344 let n_prev_x: i64 = (e_prev_dy * NX_POLYGON_Q14_ONE) / len_prev_q14
345 let n_prev_y: i64 = (0 - e_prev_dx * NX_POLYGON_Q14_ONE) / len_prev_q14
346 let n_next_x: i64 = (e_next_dy * NX_POLYGON_Q14_ONE) / len_next_q14
347 let n_next_y: i64 = (0 - e_next_dx * NX_POLYGON_Q14_ONE) / len_next_q14
348
349 // Offset displacement vectors (Q14)
350 let off_prev_x: i64 = (n_prev_x * d_q14) / NX_POLYGON_Q14_ONE
351 let off_prev_y: i64 = (n_prev_y * d_q14) / NX_POLYGON_Q14_ONE
352 let off_next_x: i64 = (n_next_x * d_q14) / NX_POLYGON_Q14_ONE
353 let off_next_y: i64 = (n_next_y * d_q14) / NX_POLYGON_Q14_ONE
354
355 // Points on each offset edge (anchor at the previous vertex
356 // of each edge)
357 let p0_prev_x: i64 = px + off_prev_x
358 let p0_prev_y: i64 = py + off_prev_y
359 let p0_next_x: i64 = cx + off_next_x
360 let p0_next_y: i64 = cy + off_next_y
361
362 // Solve P0_prev + t * (e_prev) = P0_next + s * (e_next)
363 // det = e_next_dx * e_prev_dy - e_prev_dx * e_next_dy
364 let det: i64 = e_next_dx * e_prev_dy - e_prev_dx * e_next_dy
365 let diff_x: i64 = p0_next_x - p0_prev_x
366 let diff_y: i64 = p0_next_y - p0_prev_y
367
368 var new_x: i64 = 0
369 var new_y: i64 = 0
370 if det == 0 {
371 // Parallel offset edges (straight or 180° vertex) -- just
372 // displace the vertex by the prev-edge offset normal.
373 new_x = cx + off_prev_x
374 new_y = cy + off_prev_y
375 }
376 if det != 0 {
377 // t_num = e_next_dx * diff_y - diff_x * e_next_dy
378 let t_num: i64 = e_next_dx * diff_y - diff_x * e_next_dy
379 // displacement = (t_num * e_prev_*) / det, via i128
380 // muldiv helper to avoid overflow on the t_num * e_prev_*
381 // product (would be ~2^70 for slicer-scale coords).
382 let disp_x: i64 = nx_muldiv_i64(t_num, e_prev_dx, det)
383 let disp_y: i64 = nx_muldiv_i64(t_num, e_prev_dy, det)
384 new_x = p0_prev_x + disp_x
385 new_y = p0_prev_y + disp_y
386 }
387 nx_polygon_set_vert(out, i, new_x, new_y)
388 i = i + 1
389 }
390
391 return out
392}
393
394// ===== factory: CCW square ========================================
395//
396// Convenience for smokes + slicer-substrate sanity tests. Builds
397// a CCW unit-area square spanning [x_min, x_max] × [y_min, y_max]
398// in Q14 coords.
399
400func nx_polygon_make_square(x_min_q14: i64, y_min_q14: i64,
401 x_max_q14: i64, y_max_q14: i64) -> *NxPolygon {
402 let p: *NxPolygon = nx_polygon_alloc(4)
403 if (p as i64) == 0 { return p }
404 nx_polygon_add_vert(p, x_min_q14, y_min_q14)
405 nx_polygon_add_vert(p, x_max_q14, y_min_q14)
406 nx_polygon_add_vert(p, x_max_q14, y_max_q14)
407 nx_polygon_add_vert(p, x_min_q14, y_max_q14)
408 return p
409}
410
411// ===== 3-point circle (arc fitting core) ==========================
412//
413// Given three Q14 points P1, P2, P3 (presumed non-collinear and on
414// a common circle), compute the circle's center (Q14 mm) and the
415// radius squared (Q28 mm²). Caller pre-allocates the out_* arrays.
416//
417// Formula (from analytic geometry):
418// D = 2 * (x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))
419// center_x = ((x1²+y1²)*(y2-y3) + (x2²+y2²)*(y3-y1) + (x3²+y3²)*(y1-y2)) / D
420// center_y = ((x1²+y1²)*(x3-x2) + (x2²+y2²)*(x1-x3) + (x3²+y3²)*(x2-x1)) / D
421//
422// Magnitudes: with Q14 coords up to ~2^23, squared = ~2^46. Times
423// delta-coord = ~2^67. Overflow without i128. Use nx_muldiv_i64
424// to handle the (sq * delta / D) intermediate safely.
425//
426// Returns 1 on success, 0 if D == 0 (points collinear -> no circle).
427
428func nx_polygon_3point_circle(x1: i64, y1: i64,
429 x2: i64, y2: i64,
430 x3: i64, y3: i64,
431 out_cx: *i64, out_cy: *i64,
432 out_r_sq_q28: *i64) -> i64 {
433 let dy23: i64 = y2 - y3
434 let dy31: i64 = y3 - y1
435 let dy12: i64 = y1 - y2
436 let dx32: i64 = x3 - x2
437 let dx13: i64 = x1 - x3
438 let dx21: i64 = x2 - x1
439
440 let d_q28: i64 = 2 * (x1 * dy23 + x2 * dy31 + x3 * dy12)
441 if d_q28 == 0 { return 0 } // collinear
442
443 let sq1: i64 = x1 * x1 + y1 * y1
444 let sq2: i64 = x2 * x2 + y2 * y2
445 let sq3: i64 = x3 * x3 + y3 * y3
446
447 // center_x = (sq1*dy23 + sq2*dy31 + sq3*dy12) / D
448 // Each term sq*delta/D = Q28*Q14/Q28 = Q14.
449 let cx1: i64 = nx_muldiv_i64(sq1, dy23, d_q28)
450 let cx2: i64 = nx_muldiv_i64(sq2, dy31, d_q28)
451 let cx3: i64 = nx_muldiv_i64(sq3, dy12, d_q28)
452 let cx: i64 = cx1 + cx2 + cx3
453
454 let cy1: i64 = nx_muldiv_i64(sq1, dx32, d_q28)
455 let cy2: i64 = nx_muldiv_i64(sq2, dx13, d_q28)
456 let cy3: i64 = nx_muldiv_i64(sq3, dx21, d_q28)
457 let cy: i64 = cy1 + cy2 + cy3
458
459 out_cx[0] = cx
460 out_cy[0] = cy
461
462 // r² = (x1 - cx)² + (y1 - cy)² in Q28
463 let dx1: i64 = x1 - cx
464 let dy1: i64 = y1 - cy
465 out_r_sq_q28[0] = dx1 * dx1 + dy1 * dy1
466 return 1
467}
468
469// ===== point-on-circle tolerance check ============================
470//
471// Returns 1 if (px, py) lies within tol_q14 of the circle defined
472// by (cx, cy, r_sq_q28). Tolerance compared in Q28 (squared) to
473// avoid one sqrt: |distance - r| <= tol equivalent to
474// |distance² - r²| <= 2*r*tol + tol²; for small tol vs r, ~ 2*r*tol.
475//
476// We use the cleaner exact form: compare distance² to (r-tol)² and
477// (r+tol)² — slightly larger comparison band but exact integer.
478
479func nx_polygon_point_on_circle(cx: i64, cy: i64, r_sq_q28: i64,
480 px: i64, py: i64, tol_q14: i64) -> i64 {
481 let dx: i64 = px - cx
482 let dy: i64 = py - cy
483 let d_sq_q28: i64 = dx * dx + dy * dy
484 let r_q14: i64 = nx_isqrt(r_sq_q28)
485 let lo: i64 = r_q14 - tol_q14
486 let hi: i64 = r_q14 + tol_q14
487 let lo_sq: i64 = lo * lo
488 let hi_sq: i64 = hi * hi
489 if d_sq_q28 < lo_sq { return 0 }
490 if d_sq_q28 > hi_sq { return 0 }
491 return 1
492}
493
494// ===== arc winding direction ======================================
495//
496// Given a circle center + two consecutive points on the circle,
497// returns +1 if going from P1 to P2 is counter-clockwise (G3),
498// -1 if clockwise (G2), 0 if degenerate. Cross product of
499// (P1 - center) × (P2 - center).
500
501func nx_polygon_arc_winding(cx: i64, cy: i64,
502 p1x: i64, p1y: i64,
503 p2x: i64, p2y: i64) -> i64 {
504 let ax: i64 = p1x - cx
505 let ay: i64 = p1y - cy
506 let bx: i64 = p2x - cx
507 let by: i64 = p2y - cy
508 let cross: i64 = ax * by - ay * bx
509 if cross > 0 { return 1 } // CCW = G3
510 if cross < 0 { return -1 } // CW = G2
511 return 0
512}