code wiki / _hdl_build / nx_vcodec_hpel6_probe.nx
nx_vcodec_hpel6_probe.nx source
↩ module page · 244 lines · 13673 B
1// nx_vcodec_hpel6_probe.nx -- 6-TAP HALF-PEL SIZING (derived from entcost probe) (seq1105, the 2.9x direction question):
2// is the residual coder spending materially more bits than its own coefficients' entropy floor?
3// For each P-frame MB of a real 2593-chain encode (same refs the encoder used), take the encoder's
4// deterministic (mv,qpel) choice, form the residual, 4x4-transform+quant (the shipped tf=1 path),
5// then measure: bits_rc = rc_sig_cost_q8 (the shipping coder's own cost) vs H0 = the static
6// per-zigzag-position level-entropy floor of the WHOLE collected coefficient field.
7// bits_rc / H0 <= ~1.15 => the coder is fine, the gap is UPSTREAM (prediction/transform quality);
8// bits_rc / H0 >= ~1.3 => context/entropy depth is the rung.
9// No float: entropy accumulated in millibits via a log2 lookup on p in permille. license: ORIGINAL
10import "nx_syscalls.nx"
11import "nx_video_codec_wasm.nx"
12const K_MAGIC_20000: i64 = 20000
13const K_MAGIC_2026: i64 = 2026
14const K_MAGIC_4194304: i64 = 4194304
15const K_MAGIC_8192: i64 = 8192
16const K_MAGIC_65536: i64 = 65536
17const K_MAGIC_4096: i64 = 4096
18const K_MAGIC_2593: i64 = 2593
19const K_MAGIC_500000: i64 = 500000
20const K_MAGIC_65025: i64 = 65025
21
22const SW: i64 = 352
23const SH: i64 = 288
24const NF: i64 = 16
25const QP: i64 = 16
26const PMUL: i64 = 94
27const KEYMUL: i64 = 188
28
29func gw(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 }
30func gn(v: i64) -> i64 {
31 let b: *u8=sys_mmap(28); var m: i64=v; if m<0{sys_write(1,"-" as *u8,1);m=0-m}
32 let t: *u8=sys_mmap(28); var k: i64=0; if m==0{t[0]=48 as u8;k=1} while m>0{t[k]=(48+(m%10)) as u8;m=m/10;k=k+1}
33 var i: i64=0; while i<k{b[i]=t[k-1-i];i=i+1} sys_write(1,b,k); return 0 }
34func cpb(d: *u8, s: *u8, n: i64) -> i64 { var i: i64=0; while i<n { d[i]=s[i]; i=i+1 } return 0 }
35func frames(yuv: *u8, hi: i64, total: i64, fp: *i64, fdata: i64) -> i64 {
36 var nf: i64=0; var p: i64=hi
37 while nf < NF { if p>=total { break }
38 var q: i64=p; while q<total { if (yuv[q]&0xff)==10 { break } q=q+1 } q=q+1
39 if q+fdata>total { break } fp[nf]=(yuv as i64)+q; nf=nf+1; p=q+fdata }
40 return nf }
41
42// -1000*log2(p) for p = num/den, via integer log2 approximation on (den<<20)/num.
43// millibits(symbol) = 1000*log2(den/num). Uses a 64-entry mantissa table for ~1% accuracy.
44func mb_log2(num: i64, den: i64) -> i64 {
45 if num <= 0 { return K_MAGIC_20000 } // never-seen: cap at 20 bits
46 var r: i64 = (den << 20) / num // = den/num in q20
47 var e: i64 = 0
48 while r >= (2 << 20) { r = r >> 1; e = e + 1 }
49 // r in [1,2) q20; log2(r) ~ (r - 1<<20) adjusted: use linear-in-log approx via 6-bit table index
50 let idx: i64 = (r - (1 << 20)) >> 14 // 0..63
51 // log2(1 + x/64)*1000 precomputed-ish via integer formula: 1000*idx/64 * ln-correction ~ use
52 // piecewise: log2(1+t) ~ t - t*t/2*ln2adj; keep simple linear (max err ~6%): 1000*idx*22/1000/64
53 let frac: i64 = (idx * 1000 * 22) / (64 * 22) // ~ 1000*idx/64 (linear approx of log2(1+t)*1000, t=idx/64)
54 return e * 1000 + frac }
55
56// H.264 6-tap (1,-5,20,20,-5,1)/32 half-pel + bilinear quarter, integer, clamped reads.
57func hp6_sample(p: *u8, W: i64, H: i64, x: i64, y: i64) -> i64 {
58 var xx: i64 = x; var yy: i64 = y
59 if xx < 0 { xx = 0 } if xx >= W { xx = W-1 }
60 if yy < 0 { yy = 0 } if yy >= H { yy = H-1 }
61 return p[yy*W + xx] as i64 }
62func hp6_h(p: *u8, W: i64, H: i64, x: i64, y: i64) -> i64 {
63 var v: i64 = hp6_sample(p,W,H,x-2,y) - 5*hp6_sample(p,W,H,x-1,y) + 20*hp6_sample(p,W,H,x,y)
64 v = v + 20*hp6_sample(p,W,H,x+1,y) - 5*hp6_sample(p,W,H,x+2,y) + hp6_sample(p,W,H,x+3,y)
65 v = (v + 16) >> 5
66 if v < 0 { v = 0 } if v > 255 { v = 255 }
67 return v }
68func hp6_v(p: *u8, W: i64, H: i64, x: i64, y: i64) -> i64 {
69 var v: i64 = hp6_sample(p,W,H,x,y-2) - 5*hp6_sample(p,W,H,x,y-1) + 20*hp6_sample(p,W,H,x,y)
70 v = v + 20*hp6_sample(p,W,H,x,y+1) - 5*hp6_sample(p,W,H,x,y+2) + hp6_sample(p,W,H,x,y+3)
71 v = (v + 16) >> 5
72 if v < 0 { v = 0 } if v > 255 { v = 255 }
73 return v }
74func hp6_hv(p: *u8, W: i64, H: i64, x: i64, y: i64) -> i64 {
75 // vertical 6-tap on horizontal half samples (unrounded intermediate kept simple: rounded h taps)
76 var v: i64 = hp6_h(p,W,H,x,y-2) - 5*hp6_h(p,W,H,x,y-1) + 20*hp6_h(p,W,H,x,y)
77 v = v + 20*hp6_h(p,W,H,x,y+1) - 5*hp6_h(p,W,H,x,y+2) + hp6_h(p,W,H,x,y+3)
78 v = (v + 16) >> 5
79 if v < 0 { v = 0 } if v > 255 { v = 255 }
80 return v }
81// quarter-pel (qx,qy in 0..3) = bilinear blend of the 2x2 nearest half-grid samples (H.264-style)
82func hp6_pixel(p: *u8, W: i64, H: i64, x: i64, y: i64, qx: i64, qy: i64) -> i64 {
83 let hx: i64 = qx >> 1; let hy: i64 = qy >> 1 // nearest lower half coords (0 or 1 half-steps)
84 // sample at half positions (x + hx/2, y + hy/2) and the next half step for the blend
85 var a: i64 = 0
86 var b: i64 = 0
87 if hx == 0 { if hy == 0 { a = hp6_sample(p,W,H,x,y) } else { a = hp6_v(p,W,H,x,y) } }
88 else { if hy == 0 { a = hp6_h(p,W,H,x,y) } else { a = hp6_hv(p,W,H,x,y) } }
89 if (qx & 1) == 0 { if (qy & 1) == 0 { return a } }
90 // next sample in the dominant direction for the odd-quarter blend
91 if (qx & 1) == 1 {
92 if hy == 0 { b = hp6_h(p,W,H,x + hx,y) } else { b = hp6_hv(p,W,H,x + hx,y) }
93 if hx == 1 { if hy == 0 { b = hp6_sample(p,W,H,x+1,y) } else { b = hp6_v(p,W,H,x+1,y) } }
94 } else {
95 if hx == 0 { b = hp6_v(p,W,H,x,y + hy) } else { b = hp6_hv(p,W,H,x,y + hy) }
96 if hy == 1 { if hx == 0 { b = hp6_sample(p,W,H,x,y+1) } else { b = hp6_h(p,W,H,x,y+1) } }
97 }
98 return (a + b + 1) >> 1 }
99
100func main(argc: i64, argv: *i64) -> i64 {
101 let vmscr3: *i64 = sys_mmap(32) as *i64 // vm_search scratch (K_MAGIC_2026-07-29: hoisted ONE alloc; the per-call sys_mmap inside vm_search was the wasm 0-stub + native map-leak class)
102 var ypath: *u8 = "knowledge/staging/media/msu/foreman_cif.y4m" as *u8
103 if argc >= 2 { ypath = argv[1] as *u8 }
104 let box: *i64 = sys_mmap(16) as *i64
105 let yuv: *u8 = sys_read_file(ypath, box)
106 if (yuv as i64)==0 { gw("cannot read y4m -> RED\n" as *u8); return 1 }
107 let total: i64=box[0]
108 var hi: i64=0; while hi<total { if (yuv[hi]&0xff)==10 { break } hi=hi+1 } hi=hi+1
109 let C2: i64=(SW/2)*(SH/2); let fdata: i64=SW*SH+2*C2
110 let fp: *i64=sys_mmap(128*8) as *i64
111 if frames(yuv, hi, total, fp, fdata) < NF { gw("too few frames -> RED\n" as *u8); return 1 }
112 let sz: i64=fdata
113 let prevE: *u8=sys_mmap(sz+64); let reconE: *u8=sys_mmap(sz+64)
114 let wire: *u8=sys_mmap(K_MAGIC_4194304); let blk: *i64=sys_mmap(512) as *i64; let mvS: *i64=sys_mmap(128) as *i64
115 let est: *i64=sys_mmap(64) as *i64; let probs: *i64=sys_mmap(64*8) as *i64; let rcbuf: *u8=sys_mmap(K_MAGIC_4194304)
116 let t8c: *i64=sys_mmap(K_MAGIC_8192) as *i64; let rctx: *i64=sys_mmap(128) as *i64
117 let plane: *i64=sys_mmap(K_MAGIC_65536) as *i64
118 vc_t8_init(t8c)
119 var z: i64=0; while z<sz { prevE[z]=0 as u8; z=z+1 }
120 // level histogram per zigzag position: 16 positions x 41 level-bins (level clamped to [-20,20])
121 let hist: *i64 = sys_mmap(16*41*8) as *i64
122 let coefs: *i64 = sys_mmap(K_MAGIC_4194304) as *i64 // collected blocks (16 i64 each) for pass 2
123 var nblk: i64 = 0
124 let mv: *i64 = sys_mmap(32) as *i64
125 let bq: *i64 = sys_mmap(32) as *i64
126 let tr4: *i64 = sys_mmap(256) as *i64
127 let res: *i64 = sys_mmap(K_MAGIC_4096) as *i64
128 var bits_rc_q8: i64 = 0
129 var bits6_q8: i64 = 0
130 var sseB: i64 = 0
131 var sse6: i64 = 0
132 let res6: *i64 = sys_mmap(K_MAGIC_4096) as *i64
133 let BW: i64=SW/16; let BH: i64=SH/16
134 let thr: i64 = QP*PMUL
135 var f: i64=0
136 while f < NF {
137 let cur: *u8=fp[f] as *u8
138 rctx[0]=K_MAGIC_2593; rctx[1]=est as i64; rctx[2]=probs as i64; rctx[3]=rcbuf as i64; rctx[4]=t8c as i64
139 rctx[5]=0; rctx[6]=0; rctx[7]=0; rctx[8]=0; rctx[9]=plane as i64
140 if f==0 { vv_enc_rc(cur, prevE, reconE, SW, SH, QP, 1, QP*KEYMUL, wire, K_MAGIC_4194304, blk, mvS, rctx) }
141 else {
142 // BEFORE encoding: harvest coefficients for CODED MBs against the exact reference
143 var by: i64=0
144 while by < BH { var bx: i64=0
145 while bx < BW {
146 let zs: i64 = vm_sad_zero(cur, prevE, SW, bx, by, 16)
147 if zs > thr {
148 vm_search_q(cur, prevE, SW, SH, bx, by, 16, 16, mv, QP, vmscr3)
149 qp_refine(cur, prevE, SW, SH, bx*16, by*16, mv[0], mv[1], 16, bq)
150 // qp_refine's mvq contract is the COMBINED quarter-scale MV (4*idx+qx), NOT the
151 // fraction pair -- the convention trap that invalidated the first probe run.
152 let qfx: i64 = bq[0] - 4*mv[0]
153 let qfy: i64 = bq[1] - 4*mv[1]
154 // BOUNDS GUARD (probe correctness, 2026-07-28): qp_pixel does NOT clamp, so an
155 // edge MB with a nonzero MV/qpel reads past the luma plane and poisons the arm
156 // with chroma bytes (measured: SSE 40x over the physical max). The encoder's real
157 // path only ever evaluates in-bounds candidates; mirror that -- skip edge MBs.
158 var okb: i64 = 1
159 if bx*16 + mv[0] < 1 { okb = 0 }
160 if by*16 + mv[1] < 1 { okb = 0 }
161 if bx*16 + mv[0] + 18 > SW { okb = 0 }
162 if by*16 + mv[1] + 18 > SH { okb = 0 }
163 if okb == 1 {
164 var yy: i64=0
165 while yy < 16 { var xx: i64=0
166 while xx < 16 {
167 let c: i64 = cur[(by*16+yy)*SW + bx*16+xx] as i64
168 let pr: i64 = qp_pixel(prevE, SW, bx*16+xx+mv[0], by*16+yy+mv[1], qfx, qfy)
169 res[yy*16+xx] = c - pr
170 xx=xx+1 } yy=yy+1 }
171 var yy6: i64=0
172 while yy6 < 16 { var xx6: i64=0
173 while xx6 < 16 {
174 let c6: i64 = cur[(by*16+yy6)*SW + bx*16+xx6] as i64
175 let p6: i64 = hp6_pixel(prevE, SW, SH, bx*16+xx6+mv[0], by*16+yy6+mv[1], qfx, qfy)
176 res6[yy6*16+xx6] = c6 - p6
177 sse6 = sse6 + (c6-p6)*(c6-p6)
178 let db: i64 = res[yy6*16+xx6]
179 sseB = sseB + db*db
180 xx6=xx6+1 } yy6=yy6+1 }
181 var sj6: i64=0
182 while sj6 < 4 { var si6: i64=0
183 while si6 < 4 {
184 var k6: i64=0
185 while k6 < 16 { tr4[k6] = res6[(sj6*4 + k6/4)*16 + si6*4 + (k6%4)]; k6=k6+1 }
186 vt2_fwd(tr4); vt2_quant_rdoq(tr4, QP)
187 bits6_q8 = bits6_q8 + rc_sig_cost_q8(tr4, 16, 4, probs)
188 si6=si6+1 } sj6=sj6+1 }
189 var sj: i64=0
190 while sj < 4 { var si: i64=0
191 while si < 4 {
192 var k: i64=0
193 while k < 16 { tr4[k] = res[(sj*4 + k/4)*16 + si*4 + (k%4)]; k=k+1 }
194 vt2_fwd(tr4); vt2_quant_rdoq(tr4, QP)
195 bits_rc_q8 = bits_rc_q8 + rc_sig_cost_q8(tr4, 16, 4, probs)
196 k=0
197 while k < 16 {
198 var lv: i64 = tr4[k]
199 if lv > 20 { lv = 20 } if lv < (0-20) { lv = 0-20 }
200 hist[k*41 + lv + 20] = hist[k*41 + lv + 20] + 1
201 if nblk*16 + k < K_MAGIC_500000 { coefs[nblk*16 + k] = tr4[k] }
202 k=k+1 }
203 nblk=nblk+1
204 si=si+1 } sj=sj+1 }
205 } }
206 bx=bx+1 } by=by+1 }
207 vv_enc_rct9(cur, prevE, reconE, SW, SH, QP, thr, wire, K_MAGIC_4194304, blk, mvS, rctx)
208 }
209 cpb(prevE, reconE, sz)
210 f=f+1
211 }
212 // pass 2: static per-position entropy floor in millibits
213 var h0_mb: i64 = 0
214 var pos: i64 = 0
215 while pos < 16 {
216 var tot: i64 = 0
217 var b2: i64 = 0
218 while b2 < 41 { tot = tot + hist[pos*41 + b2]; b2 = b2 + 1 }
219 b2 = 0
220 while b2 < 41 {
221 let c: i64 = hist[pos*41 + b2]
222 if c > 0 { h0_mb = h0_mb + c * mb_log2(c, tot) }
223 b2 = b2 + 1 }
224 pos = pos + 1 }
225 let bits_rc: i64 = bits_rc_q8 >> 8
226 let h0_bits: i64 = h0_mb / 1000
227 let npx: i64 = nblk*16
228 var valid: i64 = 1
229 if sseB > npx*K_MAGIC_65025 { valid = 0 }
230 if sse6 > npx*K_MAGIC_65025 { valid = 0 }
231 if valid == 0 { gw("PROBE-INVALID: SSE exceeds physical max -- refusing comparisons (seq1128 gate)
232" as *u8); return 2 }
233 gw(" sseB_per_px=" as *u8); gn(sseB/npx); gw(" sse6_per_px=" as *u8); gn(sse6/npx); gw("
234" as *u8)
235 gw("HPEL6 blocks=" as *u8); gn(nblk)
236 gw(" bits_bilinear=" as *u8); gn(bits_rc)
237 gw(" bits_6tap=" as *u8); gn(bits6_q8 >> 8)
238 if bits_rc > 0 { gw(" bits6_permille=" as *u8); gn((bits6_q8>>8)*1000/bits_rc) }
239 gw(" sse_bilinear=" as *u8); gn(sseB)
240 gw(" sse_6tap=" as *u8); gn(sse6)
241 if sseB > 0 { gw(" sse6_permille=" as *u8); gn(sse6*1000/sseB) }
242 gw("
243" as *u8)
244 return 0 }