nx_uv_stretch_candidate_t281.nx source
↩ module page · 87 lines · 6617 B
1// nx_uv_stretch_candidate_t281.nx -- Optimizes UV stretch energy using symmetric Dirichlet penalties and gradient-based updates.
2import "nx_uv_gram_solver_candidate_t280.nx"
3// Candidate refinement of the existing chart solve. Symmetric Dirichlet
4// penalizes both singular-value expansion and contraction. No mesh changes.
5struct UvStretch{
6 s:*UvGramSolve,storage:*u8,bytes:i64,area:*f64,g:*f64,offset:*i64,incident:*i64,
7 scale:f64,initial:f64,energy:f64,accepted:i64,rejected:i64,evaluations:i64,sweeps:i64,
8 h00:f64,h01:f64,h11:f64,gu:f64,gv:f64,base:f64
9}
10func ust_abs(x:f64)->f64{if x<0.0{return 0.0-x};return x}
11func ust_free(q:*UvStretch)->i64{if q==(0 as *UvStretch){return 0};if q.storage!=(0 as *u8){sys_munmap_direct(q.storage,q.bytes)};return sys_munmap_direct(q as *u8,__size_of(UvStretch))}
12// UV gradient and exact 2x2 Hessian with respect to one vertex.
13// During a single-vertex update determinant is affine, so a positive
14// endpoint plus the positive start preserves each incident face's sign.
15func ust_face(q:*UvStretch,t:i64,v:i64,u:f64,w:f64,deriv:i64)->f64{
16 let s:*UvGramSolve=q.s;var j00:f64=0.0;var j01:f64=0.0;var j10:f64=0.0;var j11:f64=0.0;var gx:f64=0.0;var gy:f64=0.0
17 var k:i64=0;while k<3{let z:i64=s.corners[t*3+k];let dx:f64=q.g[t*6+k];let dy:f64=q.g[t*6+3+k]
18 var x:f64=s.x[z];var y:f64=s.x[s.n+z];if z==v{x=u;y=w;gx=dx;gy=dy}
19 j00=j00+x*dx;j01=j01+x*dy;j10=j10+y*dx;j11=j11+y*dy;k=k+1
20 }
21 let d:f64=j00*j11-j01*j10;if d<=0.0{return 0.0-1.0};if uls_finite(d)==0{return 0.0-1.0}
22 let norm:f64=j00*j00+j01*j01+j10*j10+j11*j11
23 let inv:f64=1.0/d;let inv2:f64=inv*inv;let a:f64=q.area[t];let e:f64=a*norm*(1.0+inv2);if uls_finite(e)==0{return 0.0-1.0}
24 if deriv==1{
25 let vu:f64=j00*gx+j01*gy;let vv:f64=j10*gx+j11*gy
26 let hu:f64=j11*gx-j10*gy;let hv:f64=j00*gy-j01*gx
27 let c:f64=2.0*a*(1.0+inv2)*(gx*gx+gy*gy)
28 let b:f64=4.0*a*inv2*inv;let h:f64=6.0*a*norm*inv2*inv2
29 q.gu=q.gu+2.0*a*((1.0+inv2)*vu-norm*inv2*inv*hu)
30 q.gv=q.gv+2.0*a*((1.0+inv2)*vv-norm*inv2*inv*hv)
31 q.h00=q.h00+c-2.0*b*vu*hu+h*hu*hu
32 q.h01=q.h01-b*(vu*hv+vv*hu)+h*hu*hv
33 q.h11=q.h11+c-2.0*b*vv*hv+h*hv*hv;q.base=q.base+c
34 }
35 return e
36}
37func ust_energy(q:*UvStretch)->f64{var e:f64=0.0;var t:i64=0;while t<q.s.nt{let f:f64=ust_face(q,t,0-1,0.0,0.0,0);if f<0.0{return f};e=e+f;t=t+1};return e}
38func ust_begin(s:*UvGramSolve)->*UvStretch{
39 if s==(0 as *UvGramSolve){return 0 as *UvStretch};if s.result!=UV_OK{return 0 as *UvStretch}
40 let q:*UvStretch=sys_mmap_try(__size_of(UvStretch)) as *UvStretch;if (q as i64)<=0{return 0 as *UvStretch};q.s=s
41 if s.nt>(NCP_MAX/8-2*s.n-1)/10{ust_free(q);return 0 as *UvStretch}
42 q.bytes=(10*s.nt+2*s.n+1)*8;q.storage=sys_mmap_try(q.bytes);if (q.storage as i64)<=0{q.storage=0 as *u8;ust_free(q);return 0 as *UvStretch}
43 var off:i64=q.storage as i64;q.area=off as *f64;off=off+s.nt*8;q.g=off as *f64;off=off+s.nt*6*8;q.offset=off as *i64;off=off+(s.n+1)*8;q.incident=off as *i64;off=off+s.nt*3*8;let cursor:*i64=off as *i64
44 var totalUV:f64=0.0;var totalSource:f64=0.0;var t:i64=0;let first:i64=s.m[s.m[UV_O_CSTART]+s.c]
45 while t<s.nt{
46 let st:i64=s.m[s.m[UV_O_CORDER]+first+t];let ia:i64=uv_ti(s.m,st,0);let ib:i64=uv_ti(s.m,st,1);let ic:i64=uv_ti(s.m,st,2)
47 let ax:f64=__f64_from_i64(uv_vx(s.m,ib)-uv_vx(s.m,ia));let ay:f64=__f64_from_i64(uv_vy(s.m,ib)-uv_vy(s.m,ia));let az:f64=__f64_from_i64(uv_vz(s.m,ib)-uv_vz(s.m,ia))
48 let bx:f64=__f64_from_i64(uv_vx(s.m,ic)-uv_vx(s.m,ia));let by:f64=__f64_from_i64(uv_vy(s.m,ic)-uv_vy(s.m,ia));let bz:f64=__f64_from_i64(uv_vz(s.m,ic)-uv_vz(s.m,ia))
49 let cx:f64=ay*bz-az*by;let cy:f64=az*bx-ax*bz;let cz:f64=ax*by-ay*bx;let a:f64=uls_sqrt(cx*cx+cy*cy+cz*cz)/2.0;if a<=0.0{ust_free(q);return 0 as *UvStretch}
50 let v0:i64=s.corners[t*3];let v1:i64=s.corners[t*3+1];let v2:i64=s.corners[t*3+2]
51 let du:f64=s.x[v1]-s.x[v0];let dv:f64=s.x[s.n+v1]-s.x[s.n+v0];let eu:f64=s.x[v2]-s.x[v0];let ev:f64=s.x[s.n+v2]-s.x[s.n+v0];let areaUV:f64=(du*ev-dv*eu)/2.0
52 if areaUV<=0.0{ust_free(q);return 0 as *UvStretch};totalUV=totalUV+areaUV;totalSource=totalSource+a;q.area[t]=a
53 let root:f64=uls_sqrt(a);var k:i64=0;while k<6{q.g[t*6+k]=s.coeff[t*6+k]/root;k=k+1}
54 k=0;while k<3{let v:i64=s.corners[t*3+k];q.offset[v+1]=q.offset[v+1]+1;k=k+1};t=t+1
55 }
56 q.scale=uls_sqrt(totalUV/totalSource);if q.scale<=0.0{ust_free(q);return 0 as *UvStretch};if uls_finite(q.scale)==0{ust_free(q);return 0 as *UvStretch}
57 var i:i64=0;while i<s.n{q.offset[i+1]=q.offset[i+1]+q.offset[i];cursor[i]=q.offset[i];i=i+1}
58 t=0;while t<s.nt{var k:i64=0;while k<3{let v:i64=s.corners[t*3+k];q.incident[cursor[v]]=t;cursor[v]=cursor[v]+1;k=k+1};t=t+1}
59 i=0;while i<s.n*2{s.x[i]=s.x[i]/q.scale;i=i+1};q.initial=ust_energy(q);q.energy=q.initial;return q
60}
61func ust_step(q:*UvStretch,sweeps:i64,backtracks:i64,cancel:*i64)->i64{
62 if q==(0 as *UvStretch){return UV_E_ARGS};if sweeps<1{return UV_E_ARGS};if backtracks<1{return UV_E_ARGS};if q.initial<0.0{return ULS_NUMERIC}
63 let s:*UvGramSolve=q.s;var sweep:i64=0
64 while sweep<sweeps{
65 if cancel!=(0 as *i64){if cancel[0]!=0{return ULS_CANCELLED}}
66 var v:i64=0;while v<s.n{
67 if v!=s.pa{if v!=s.pb{
68 q.gu=0.0;q.gv=0.0;q.h00=0.0;q.h01=0.0;q.h11=0.0;q.base=0.0;var old:f64=0.0;var i:i64=q.offset[v]
69 while i<q.offset[v+1]{let e:f64=ust_face(q,q.incident[i],v,s.x[v],s.x[s.n+v],1);if e<0.0{return ULS_NUMERIC};old=old+e;i=i+1}
70 let det:f64=q.h00*q.h11-q.h01*q.h01
71 if q.h00<=0.0{q.h00=q.h00+q.base+ust_abs(q.h00)+ust_abs(q.h01);q.h11=q.h11+q.base+ust_abs(q.h11)+ust_abs(q.h01)}else{if det<=0.0{let shift:f64=q.base+ust_abs(q.h00)+ust_abs(q.h11)+2.0*ust_abs(q.h01);q.h00=q.h00+shift;q.h11=q.h11+shift}}
72 let hd:f64=q.h00*q.h11-q.h01*q.h01;if hd<=0.0{return ULS_NUMERIC}
73 let du:f64=(q.h11*q.gu-q.h01*q.gv)/hd;let dv:f64=(q.h00*q.gv-q.h01*q.gu)/hd
74 if uls_finite(du)==0{return ULS_NUMERIC};if uls_finite(dv)==0{return ULS_NUMERIC}
75 var alpha:f64=1.0;var attempt:i64=0;var accepted:i64=0
76 while attempt<backtracks{
77 let u:f64=s.x[v]-alpha*du;let w:f64=s.x[s.n+v]-alpha*dv;var next:f64=0.0;var valid:i64=1;i=q.offset[v]
78 while i<q.offset[v+1]{let e:f64=ust_face(q,q.incident[i],v,u,w,0);q.evaluations=q.evaluations+1;if e<0.0{valid=0;break};next=next+e;i=i+1}
79 if valid==1{if next<old{s.x[v]=u;s.x[s.n+v]=w;accepted=1;q.accepted=q.accepted+1;break}}
80 alpha=alpha/2.0;attempt=attempt+1
81 };if accepted==0{q.rejected=q.rejected+1}
82 }};v=v+1
83 };q.sweeps=q.sweeps+1;sweep=sweep+1
84 }
85 q.energy=ust_energy(q);if q.energy<0.0{return ULS_NUMERIC};return UV_OK
86}
87func ust_restore_scale(q:*UvStretch)->i64{var i:i64=0;while i<q.s.n*2{q.s.x[i]=q.s.x[i]*q.scale;i=i+1};return UV_OK}