nx_partid.nx source

1// nx_partid.nx -- geometric PART IDENTIFICATION (cadtwin P5): recognize a mesh part's shape CATEGORY by a 2// ROTATION- and SCALE-INVARIANT descriptor = the eigenvalue ratios of the vertex covariance (principal-moment 3// signature). A rod -> one dominant axis (long,thin,thin); a plate -> two dominant (flat); a cube -> three equal. 4// Invariant BY CONSTRUCTION: covariance rotates as R C R^T (eigenvalues unchanged) and scales by s^2 (ratios 5// unchanged). Match a query descriptor to the nearest catalog entry. Composes nx_mesh3 (verts) + nx_epose 6// (power-iteration 3x3 eigensolver). HONEST: coarse shape CLASS (fastener/plate/bracket/rod), NOT fine OEM 7// part-number ID -- that needs the EPC catalog corpus + finer/neural matching (the fuller P5 rung, benchmark-free 8// per the research). Deterministic. license_tier: ORIGINAL 9import "nx_mesh3.nx" 10import "nx_epose.nx" 11 12// vertex covariance (3x3 symmetric, row-major) = Sum (v-centroid)(v-centroid)^T (no /n; ratios are scale-free) 13func pid_cov(base: i64, C: *i64) -> i64 { 14 let h: *i64 = m3_hdr(base) 15 let n: i64 = h[0] 16 var cx: i64 = 0 17 var cy: i64 = 0 18 var cz: i64 = 0 19 var i: i64 = 0 20 while i < n { 21 let v: *i64 = m3_vert(base, i) 22 cx = cx + v[0]; cy = cy + v[1]; cz = cz + v[2] 23 i = i + 1 24 } 25 cx = cx / n; cy = cy / n; cz = cz / n 26 var i2: i64 = 0 27 while i2 < 9 { C[i2] = 0; i2 = i2 + 1 } 28 i = 0 29 while i < n { 30 let v: *i64 = m3_vert(base, i) 31 let dx: i64 = (v[0] - cx) / 16 // /16 to keep products in range (ratios unaffected) 32 let dy: i64 = (v[1] - cy) / 16 33 let dz: i64 = (v[2] - cz) / 16 34 C[0] = C[0] + dx * dx; C[1] = C[1] + dx * dy; C[2] = C[2] + dx * dz 35 C[4] = C[4] + dy * dy; C[5] = C[5] + dy * dz 36 C[8] = C[8] + dz * dz 37 i = i + 1 38 } 39 C[3] = C[1]; C[6] = C[2]; C[7] = C[5] // symmetric 40 return 0 41} 42// Rayleigh quotient v^T C v (v unit Q14, C raw) -> C-scale 43func pid_rayleigh(C: *i64, v: *i64) -> i64 { 44 let cv: *i64 = sys_mmap(32) as *i64 45 ep_mv(C, v[0], v[1], v[2], cv) // NOTE ep_mv divides by EP_Q -> cv in C-scale/Q14... see below 46 return (v[0] * cv[0] + v[1] * cv[1] + v[2] * cv[2]) / EP_Q 47} 48// C' = C - lambda * v v^T (deflate the eigenpair) -> out (3x3) 49func pid_deflate(C: *i64, l: i64, v: *i64, out: *i64) -> i64 { 50 var r: i64 = 0 51 while r < 3 { 52 var c: i64 = 0 53 while c < 3 { 54 let outer: i64 = (v[r] * v[c]) / EP_Q 55 out[r * 3 + c] = C[r * 3 + c] - (l * outer) / EP_Q 56 c = c + 1 57 } 58 r = r + 1 59 } 60 return 0 61} 62// three eigenvalues (sorted DESC, non-negative) of symmetric PSD C into out3, via THREE power-iter dominant + 63// deflation passes. Each lambda = Rayleigh of the dominant eigenvector of the successively-deflated matrix 64// (always the positive remaining eigenvalue -- NO trace-subtraction, which cancels catastrophically for a 65// degenerate spectrum like a plate's lambda1==lambda2). 66func pid_eigvals(C: *i64, out3: *i64) -> i64 { 67 let Cs: *i64 = sys_mmap(128) as *i64 68 var maxa: i64 = 1 69 var i: i64 = 0 70 while i < 9 { var aa: i64 = C[i]; if aa < 0 { aa = 0 - aa } if aa > maxa { maxa = aa } i = i + 1 } 71 let scl: i64 = (maxa / EP_Q) + 1 72 i = 0 73 while i < 9 { Cs[i] = C[i] / scl; i = i + 1 } 74 let v1: *i64 = sys_mmap(32) as *i64 75 let v2: *i64 = sys_mmap(32) as *i64 76 let v3: *i64 = sys_mmap(32) as *i64 77 let Cs2: *i64 = sys_mmap(128) as *i64 78 let Cs3: *i64 = sys_mmap(128) as *i64 79 ep_eig_dom(Cs, v1) 80 var l1: i64 = pid_rayleigh(Cs, v1) 81 pid_deflate(Cs, l1, v1, Cs2) 82 ep_eig_dom(Cs2, v2) 83 var l2: i64 = pid_rayleigh(Cs2, v2) 84 pid_deflate(Cs2, l2, v2, Cs3) 85 ep_eig_dom(Cs3, v3) 86 var l3: i64 = pid_rayleigh(Cs3, v3) 87 if l1 < 0 { l1 = 0 } 88 if l2 < 0 { l2 = 0 } 89 if l3 < 0 { l3 = 0 } 90 var a: i64 = l1 91 var b: i64 = l2 92 var cc: i64 = l3 93 if b < cc { let t: i64 = b; b = cc; cc = t } 94 if a < b { let t: i64 = a; a = b; b = t } 95 if b < cc { let t: i64 = b; b = cc; cc = t } 96 out3[0] = a; out3[1] = b; out3[2] = cc 97 return 0 98} 99// descriptor: (lambda2*1000/lambda1, lambda3*1000/lambda1) permille -> out2 (rotation+scale invariant) 100func pid_descriptor(base: i64, out2: *i64) -> i64 { 101 let C: *i64 = sys_mmap(128) as *i64 102 pid_cov(base, C) 103 let ev: *i64 = sys_mmap(32) as *i64 104 pid_eigvals(C, ev) 105 if ev[0] <= 0 { out2[0] = 0; out2[1] = 0; return 0 } 106 out2[0] = (ev[1] * 1000) / ev[0] 107 out2[1] = (ev[2] * 1000) / ev[0] 108 return 0 109} 110// match descriptor d(2) to catalog (stride2 descriptors, ncat entries) -> nearest index; dist into out_dist 111func pid_match(d: *i64, cat: *i64, ncat: i64, out_dist: *i64) -> i64 { 112 var best: i64 = 0 - 1 113 var bd: i64 = 0 - 1 114 var i: i64 = 0 115 while i < ncat { 116 var e0: i64 = d[0] - cat[i * 2] 117 if e0 < 0 { e0 = 0 - e0 } 118 var e1: i64 = d[1] - cat[i * 2 + 1] 119 if e1 < 0 { e1 = 0 - e1 } 120 let dd: i64 = e0 + e1 121 if best < 0 { best = i; bd = dd } else { if dd < bd { best = i; bd = dd } } 122 i = i + 1 123 } 124 out_dist[0] = bd 125 return best 126} 127// rotate all verts of src mesh into dst by Q14 rotation R (for invariance testing) 128func pid_rotate_mesh(src: i64, dst: i64, R: *i64) -> i64 { 129 m3_init(dst) 130 let h: *i64 = m3_hdr(src) 131 var i: i64 = 0 132 while i < h[0] { 133 let v: *i64 = m3_vert(src, i) 134 let x: i64 = (R[0] * v[0] + R[1] * v[1] + R[2] * v[2]) / EP_Q 135 let y: i64 = (R[3] * v[0] + R[4] * v[1] + R[5] * v[2]) / EP_Q 136 let z: i64 = (R[6] * v[0] + R[7] * v[1] + R[8] * v[2]) / EP_Q 137 m3_add_vert(dst, x, y, z) 138 i = i + 1 139 } 140 // copy tris 141 var t: i64 = 0 142 while t < h[1] { 143 let tr: *i64 = m3_tri(src, t) 144 m3_add_tri(dst, tr[0], tr[1], tr[2]) 145 t = t + 1 146 } 147 return 0 148}