`nx_mesh_qem.nx` source

↩ module page · 187 lines · 7866 B
1// nx_mesh_qem.nx -- mesh decimation by QUADRIC ERROR METRICS (Garland-Heckbert 1997), pure
2// integer.  Each vertex accumulates a quadric Q (sum of outer products of its incident face
3// planes); collapsing an edge costs p^T Q p (the squared distance of the new point to those
4// planes).  Collapsing LOWEST-cost edges first => flat regions decimate hard while curved/feature
5// regions keep their vertices -- the quality win over uniform vertex clustering (nx_mesh_decimate)
6// on CURVED geometry.  Sovereign + no-float: the mesh is scaled into a bounded [0,1024] cube so
7// the integer quadrics never overflow i64; collapse uses MIDPOINT placement (no matrix inverse).
8// Reuses nx_mesh + nx_isqrt.  license_tier: ORIGINAL
9import "nx_syscalls.nx"
10import "nx_mesh.nx"
11import "nx_mesh_print_check.nx"
12import "nx_isqrt.nx"
13
14const QEM_S: i64 = 1024      // normalized cube side
15const QEM_NQ: i64 = 256      // normal fixed-point scale
16
17// add face plane (a,b,c,d) outer-product quadric into vertex v's 10-entry quadric.
18func qem_add(quad: *i64, v: i64, a: i64, b: i64, c: i64, d: i64) -> i64 {
  let base: i64 = v * 10
  quad[base+0] = quad[base+0] + a*a
  quad[base+1] = quad[base+1] + a*b
  quad[base+2] = quad[base+2] + a*c
  quad[base+3] = quad[base+3] + a*d
  quad[base+4] = quad[base+4] + b*b
  quad[base+5] = quad[base+5] + b*c
  quad[base+6] = quad[base+6] + b*d
  quad[base+7] = quad[base+7] + c*c
  quad[base+8] = quad[base+8] + c*d
  quad[base+9] = quad[base+9] + d*d
  return 0
31}
32
33// p^T (Q[va]+Q[vb]) p  -- quadric error of placing the merged vertex at (px,py,pz).
34func qem_cost(quad: *i64, va: i64, vb: i64, px: i64, py: i64, pz: i64) -> i64 {
  let a: i64 = va*10; let b: i64 = vb*10
  let q0: i64=quad[a+0]+quad[b+0]; let q1: i64=quad[a+1]+quad[b+1]; let q2: i64=quad[a+2]+quad[b+2]
  let q3: i64=quad[a+3]+quad[b+3]; let q4: i64=quad[a+4]+quad[b+4]; let q5: i64=quad[a+5]+quad[b+5]
  let q6: i64=quad[a+6]+quad[b+6]; let q7: i64=quad[a+7]+quad[b+7]; let q8: i64=quad[a+8]+quad[b+8]
  let q9: i64=quad[a+9]+quad[b+9]
  let r1: i64 = px*px*q0 + 2*px*py*q1 + 2*px*pz*q2 + 2*px*q3
  let r2: i64 = py*py*q4 + 2*py*pz*q5 + 2*py*q6
  let r3: i64 = pz*pz*q7 + 2*pz*q8 + q9
  return r1 + r2 + r3
44}
45
46// Decimate `m` down to about `target_tris` triangles.  Returns a new NxMesh.
47func nx_mesh_qem(m: *NxMesh, target_tris: i64) -> *NxMesh {
  if (m as i64) == 0 { return 0 as *NxMesh }
  let nv: i64 = m.n_verts
  let nt: i64 = m.n_tris
  if nv <= 0 { return 0 as *NxMesh }
  if nt <= 0 { return 0 as *NxMesh }
53
  let bb: *NxMeshBBox = nx_mesh_bbox_compute(m)
  let minx: i64 = bb.min_x; let miny: i64 = bb.min_y; let minz: i64 = bb.min_z
  var ext: i64 = bb.max_x - minx
  if bb.max_y - miny > ext { ext = bb.max_y - miny }
  if bb.max_z - minz > ext { ext = bb.max_z - minz }
  if ext <= 0 { return 0 as *NxMesh }
60
  // working arrays
  let vx: *i64 = (sys_mmap(nv*8)) as *i64; let vy: *i64 = (sys_mmap(nv*8)) as *i64; let vz: *i64 = (sys_mmap(nv*8)) as *i64
  let valive: *i64 = (sys_mmap(nv*8)) as *i64
  let quad: *i64 = (sys_mmap(nv*10*8)) as *i64
  let tri: *i64 = (sys_mmap(nt*3*8)) as *i64
  let talive: *i64 = (sys_mmap(nt*8)) as *i64
67
  var v: i64 = 0
  while v < nv {
      vx[v] = (nx_mesh_get_vertex_x(m,v) - minx) * QEM_S / ext
      vy[v] = (nx_mesh_get_vertex_y(m,v) - miny) * QEM_S / ext
      vz[v] = (nx_mesh_get_vertex_z(m,v) - minz) * QEM_S / ext
      valive[v] = 1
      var k: i64 = 0; while k < 10 { quad[v*10+k] = 0; k = k + 1 }
      v = v + 1
  }
  var t: i64 = 0
  while t < nt {
      tri[t*3+0] = m.indices[t*3+0]; tri[t*3+1] = m.indices[t*3+1]; tri[t*3+2] = m.indices[t*3+2]
      talive[t] = 1
      t = t + 1
  }
83
  // accumulate face-plane quadrics
  t = 0
  while t < nt {
      let v0: i64 = tri[t*3+0]; let v1: i64 = tri[t*3+1]; let v2: i64 = tri[t*3+2]
      let e1x: i64=vx[v1]-vx[v0]; let e1y: i64=vy[v1]-vy[v0]; let e1z: i64=vz[v1]-vz[v0]
      let e2x: i64=vx[v2]-vx[v0]; let e2y: i64=vy[v2]-vy[v0]; let e2z: i64=vz[v2]-vz[v0]
      let nx2: i64 = e1y*e2z - e1z*e2y
      let ny2: i64 = e1z*e2x - e1x*e2z
      let nz2: i64 = e1x*e2y - e1y*e2x
      let ln: i64 = nx_isqrt(nx2*nx2 + ny2*ny2 + nz2*nz2)
      if ln > 0 {
          let a: i64 = nx2*QEM_NQ/ln; let b: i64 = ny2*QEM_NQ/ln; let c: i64 = nz2*QEM_NQ/ln
          let d: i64 = 0 - (a*vx[v0] + b*vy[v0] + c*vz[v0])
          qem_add(quad, v0, a,b,c,d); qem_add(quad, v1, a,b,c,d); qem_add(quad, v2, a,b,c,d)
      }
      t = t + 1
  }
101
  // greedy collapse: repeatedly merge the lowest-cost live edge (midpoint) until target.
  var live_t: i64 = nt
  while live_t > target_tris {
      var best: i64 = 0 - 1
      var bva: i64 = 0 - 1; var bvb: i64 = 0 - 1
      var bpx: i64 = 0; var bpy: i64 = 0; var bpz: i64 = 0
      var found: i64 = 0
      t = 0
      while t < nt {
          if talive[t] == 1 {
              var e: i64 = 0
              while e < 3 {
                  var va: i64 = tri[t*3+e]
                  var vb: i64 = tri[t*3+((e+1)%3)]
                  if valive[va] == 1 { if valive[vb] == 1 { if va != vb {
                      let px: i64 = (vx[va]+vx[vb])/2; let py: i64 = (vy[va]+vy[vb])/2; let pz: i64 = (vz[va]+vz[vb])/2
                      let c2: i64 = qem_cost(quad, va, vb, px, py, pz)
                      if found == 0 { found=1; best=c2; bva=va; bvb=vb; bpx=px; bpy=py; bpz=pz }
                      else { if c2 < best { best=c2; bva=va; bvb=vb; bpx=px; bpy=py; bpz=pz } }
                  } } }
                  e = e + 1
              }
          }
          t = t + 1
      }
      if found == 0 { live_t = target_tris }     // nothing collapsible -> stop
      else {
          // collapse vb -> va at midpoint
          vx[bva]=bpx; vy[bva]=bpy; vz[bva]=bpz
          var k: i64 = 0; while k < 10 { quad[bva*10+k] = quad[bva*10+k] + quad[bvb*10+k]; k = k + 1 }
          valive[bvb] = 0
          // retarget triangles + drop degenerates
          t = 0
          while t < nt {
              if talive[t] == 1 {
                  var i0: i64 = tri[t*3+0]; var i1: i64 = tri[t*3+1]; var i2: i64 = tri[t*3+2]
                  if i0 == bvb { i0 = bva }
                  if i1 == bvb { i1 = bva }
                  if i2 == bvb { i2 = bva }
                  tri[t*3+0]=i0; tri[t*3+1]=i1; tri[t*3+2]=i2
                  if i0==i1 { talive[t]=0; live_t=live_t-1 }
                  else { if i1==i2 { talive[t]=0; live_t=live_t-1 }
                  else { if i0==i2 { talive[t]=0; live_t=live_t-1 } } }
              }
              t = t + 1
          }
      }
  }
150
  // reindex live vertices
  let vmap: *i64 = (sys_mmap(nv*8)) as *i64
  var outv: i64 = 0
  v = 0
  while v < nv { if valive[v]==1 { vmap[v]=outv; outv=outv+1 } else { vmap[v]=0-1 } v=v+1 }
  // count live tris
  var outt: i64 = 0
  t = 0
  while t < nt { if talive[t]==1 { outt=outt+1 } t=t+1 }
  if outv <= 0 { return 0 as *NxMesh }
  if outt <= 0 { return 0 as *NxMesh }
162
  let out: *NxMesh = nx_mesh_alloc(outv, outt, 0)
  if (out as i64)==0 { return out }
  // emit verts scaled back to original space
  v = 0
  while v < nv {
      if valive[v]==1 {
          let ox: i64 = minx + vx[v]*ext/QEM_S
          let oy: i64 = miny + vy[v]*ext/QEM_S
          let oz: i64 = minz + vz[v]*ext/QEM_S
          nx_mesh_set_vertex(out, vmap[v], ox, oy, oz, 0)
      }
      v = v + 1
  }
  // emit live tris remapped
  var ti: i64 = 0
  t = 0
  while t < nt {
      if talive[t]==1 {
          nx_mesh_set_triangle(out, ti, vmap[tri[t*3+0]], vmap[tri[t*3+1]], vmap[tri[t*3+2]])
          ti = ti + 1
      }
      t = t + 1
  }
  return out
187}
nx_mesh_qem.nx source

`nx_mesh_qem.nx` source
