diff --git a/crates/export/src/lib.rs b/crates/export/src/lib.rs index 941a7723..c9196625 100644 --- a/crates/export/src/lib.rs +++ b/crates/export/src/lib.rs @@ -7,7 +7,7 @@ //! | PNG / JPEG / TIFF sequence | one file per frame | 8-bit | PNG, TIFF | — | //! | OpenEXR sequence | one file per frame | 32-bit float, linear light, premultiplied | yes | — | //! | Animated GIF | GIF89a | 256-colour palette per frame (NeuQuant) | 1-bit | — | -//! | WebM | WebM (Matroska) | VP9 profile 0 intra frames (`effectcraft-vp9enc`), 8-bit 4:2:0 | VP9 alpha (BlockAdditional) | Opus (`effectcraft-opusenc`), 48 kHz stereo | +//! | WebM | WebM (Matroska) | VP9 profile 0 (`effectcraft-vp9enc`: key + inter frames, a key frame every 2 s, loop filter), 8-bit 4:2:0 | VP9 alpha (BlockAdditional) | Opus (`effectcraft-opusenc`), 48 kHz stereo | //! | WAV / AIFF | RIFF WAVE / AIFF | — | — | 16-bit PCM stereo (audio only) | //! //! Frames are rendered in parallel batches (rayon; one batch ≈ one frame per core) and handed to diff --git a/crates/export/src/webm.rs b/crates/export/src/webm.rs index eb99cd22..65ca1189 100644 --- a/crates/export/src/webm.rs +++ b/crates/export/src/webm.rs @@ -5,8 +5,10 @@ //! (size patched at the end), Info, Tracks and one Cluster per render batch. Video frames are //! SimpleBlocks; with alpha (Channels: RGB + Alpha) each frame is a BlockGroup whose //! BlockAdditional (BlockAddID 1) carries a second VP9 frame coding the alpha channel as luma, -//! with `AlphaMode` 1 on the track, as WebM defines. VP9 comes from `effectcraft-vp9enc` (every -//! frame a key frame), Opus from `effectcraft-opusenc` (CELT, 48 kHz stereo, 20 ms packets, +//! with `AlphaMode` 1 on the track, as WebM defines. VP9 comes from `effectcraft-vp9enc` (key +//! frames every two seconds, inter frames between them: SimpleBlocks carry the key flag, inter +//! BlockGroups a ReferenceBlock; the alpha stream's key frames line up with the colour stream's), +//! Opus from `effectcraft-opusenc` (CELT, 48 kHz stereo, 20 ms packets, //! `CodecDelay` = pre-skip). WAV (RIFF, 16-bit PCM) and AIFF (big-endian 16-bit, 80-bit //! extended sample rate) follow their published file layouts. @@ -82,6 +84,10 @@ fn block(track: u64, rel: i16, flags: u8, data: &[u8]) -> Vec { struct Frame { ms: i64, track: u64, + /// A key frame (VP9) or any audio packet. + key: bool, + /// Timestamp of the previous video frame (ReferenceBlock of inter frames). + prev_ms: Option, data: Vec, alpha: Option>, } @@ -94,10 +100,15 @@ fn cluster(frames: &mut [Frame]) -> Vec { for f in frames.iter() { let rel = (f.ms - base).clamp(i16::MIN as i64, i16::MAX as i64) as i16; match &f.alpha { - None => el(&mut body, 0xA3, &block(f.track, rel, 0x80, &f.data)), + None => el(&mut body, 0xA3, &block(f.track, rel, if f.key { 0x80 } else { 0x00 }, &f.data)), Some(a) => { let mut g = vec![]; el(&mut g, 0xA1, &block(f.track, rel, 0x00, &f.data)); + if let (false, Some(p)) = (f.key, f.prev_ms) { + // ReferenceBlock: a signed timestamp relative to this block (an inter frame). + let d = (p - f.ms).clamp(i16::MIN as i64, -1) as i16; + el(&mut g, 0xFB, &d.to_be_bytes()); + } let mut more = vec![]; el_uint(&mut more, 0xEE, 1); el(&mut more, 0xA5, a); @@ -117,7 +128,19 @@ pub(crate) fn webm(job: &Job, comp: &Comp, w: u32, h: u32, st: &mut State) -> Re let rate = job.settings.rate(comp); let alpha = job.output.channels == Channels::Rgba; let quality = job.output.quality.clamp(1, 100); - let cfg = effectcraft_vp9enc::EncoderConfig { width: w, height: h, quality, full_range: false }; + let fps = rate.as_f64().max(1.0); + let cfg = effectcraft_vp9enc::EncoderConfig { + width: w, + height: h, + quality, + full_range: false, + keyframe_interval: (fps * 2.0).round().max(1.0) as u32, + target_kbps: None, + frame_rate: fps, + loop_filter: true, + }; + let mut color_enc = effectcraft_vp9enc::Vp9Encoder::new(cfg.clone()); + let mut alpha_enc = alpha.then(|| effectcraft_vp9enc::Vp9Encoder::new(cfg.clone())); let with_audio = wants_audio(job); let mut opus = with_audio.then(|| effectcraft_opusenc::OpusEncoder::new(2, OPUS_BITRATE)); let total = st.total; @@ -209,7 +232,7 @@ pub(crate) fn webm(job: &Job, comp: &Comp, w: u32, h: u32, st: &mut State) -> Re while pcm_left.len() >= fs { let pkt = enc.encode_float(&pcm_left[..fs]); pcm_left.drain(..fs); - frames.push(Frame { ms: packets as i64 * 20, track: 2, data: pkt, alpha: None }); + frames.push(Frame { ms: packets as i64 * 20, track: 2, key: true, prev_ms: None, data: pkt, alpha: None }); packets += 1; } }; @@ -218,23 +241,28 @@ pub(crate) fn webm(job: &Job, comp: &Comp, w: u32, h: u32, st: &mut State) -> Re let mut i = 0; while i < total { let end = (i + batch).min(total); - let encoded: Vec<(Vec, Option>)> = (i..end) - .into_par_iter() - .map(|k| { - let px = rgba8(&render_frame(job, comp, k), comp, if alpha { Channels::Rgba } else { Channels::Rgb }, w, h); - let mut enc = effectcraft_vp9enc::Vp9Encoder::new(cfg.clone()); - let (y, u, vv) = effectcraft_vp9enc::rgba_to_yuv420(&px, w, h, false); - let color = enc.encode_yuv420(&y, &u, &vv); - let a = alpha.then(|| { - let mut ea = effectcraft_vp9enc::Vp9Encoder::new(cfg.clone()); - let (y, u, vv) = effectcraft_vp9enc::alpha_to_yuv420(&px, w, h); - ea.encode_yuv420(&y, &u, &vv) - }); - (color, a) - }) - .collect(); - let mut frames: Vec = - encoded.into_iter().enumerate().map(|(j, (data, alpha))| Frame { ms: frame_ms(i + j as u64), track: 1, data, alpha }).collect(); + // Frames render in parallel; the VP9 streams (inter-coded) encode in order, colour and + // alpha side by side. + let pixels: Vec> = + (i..end).into_par_iter().map(|k| rgba8(&render_frame(job, comp, k), comp, if alpha { Channels::Rgba } else { Channels::Rgb }, w, h)).collect(); + let mut frames: Vec = Vec::with_capacity(pixels.len()); + for (j, px) in pixels.iter().enumerate() { + let (color, a) = rayon::join( + || { + let (y, u, vv) = effectcraft_vp9enc::rgba_to_yuv420(px, w, h, false); + color_enc.encode_frame(&y, &u, &vv, false) + }, + || { + alpha_enc.as_mut().map(|ea| { + let (y, u, vv) = effectcraft_vp9enc::alpha_to_yuv420(px, w, h); + ea.encode_frame(&y, &u, &vv, false).data + }) + }, + ); + let k = i + j as u64; + frames.push(Frame { ms: frame_ms(k), track: 1, key: color.key, prev_ms: (k > 0).then(|| frame_ms(k - 1)), data: color.data, alpha: a }); + } + drop(pixels); let t_end = if end >= total { span_end } else { job.settings.frame_time(comp, end) }; audio_until(t_end, &mut frames, end >= total); let c = cluster(&mut frames); diff --git a/crates/export/tests/vp9_roundtrip.rs b/crates/export/tests/vp9_roundtrip.rs index 8156d850..5bf46e0a 100644 --- a/crates/export/tests/vp9_roundtrip.rs +++ b/crates/export/tests/vp9_roundtrip.rs @@ -61,7 +61,7 @@ fn lossless_is_exact_for_odd_sizes() { for (w, h) in [(1, 1), (33, 17), (64, 64), (200, 120), (7, 70)] { for kind in 0..3 { let (y, u, v) = source(w, h, kind); - let mut e = Vp9Encoder::new(EncoderConfig { width: w as u32, height: h as u32, quality: 100, full_range: false }); + let mut e = Vp9Encoder::new(EncoderConfig { width: w as u32, height: h as u32, quality: 100, full_range: false, ..Default::default() }); let f = e.encode_yuv420(&y, &u, &v); let p = decode(&f); assert_eq!((p.width, p.height), (w as u32, h as u32)); @@ -80,7 +80,7 @@ fn lossy_quality_and_size() { let (y, u, v) = source(w, h, kind); let mut last = 0usize; for q in [30u8, 60, 80, 95] { - let mut e = Vp9Encoder::new(EncoderConfig { width: w as u32, height: h as u32, quality: q, full_range: false }); + let mut e = Vp9Encoder::new(EncoderConfig { width: w as u32, height: h as u32, quality: q, full_range: false, ..Default::default() }); let f = e.encode_yuv420(&y, &u, &v); let p = decode(&f); let got = crop(&plane_u8(&p.y), p.y_stride, w, h); @@ -103,7 +103,7 @@ fn wide_frames_use_tile_columns() { // 4160 px > 64 superblocks: two tile columns are required. let (w, h) = (4160, 16); let (y, u, v) = source(w, h, 0); - let mut e = Vp9Encoder::new(EncoderConfig { width: w as u32, height: h as u32, quality: 90, full_range: false }); + let mut e = Vp9Encoder::new(EncoderConfig { width: w as u32, height: h as u32, quality: 90, full_range: false, ..Default::default() }); let f = e.encode_yuv420(&y, &u, &v); let p = decode(&f); assert!(psnr(&crop(&plane_u8(&p.y), p.y_stride, w, h), &y) > 35.0); @@ -113,7 +113,7 @@ fn wide_frames_use_tile_columns() { fn ffmpeg_decodes_our_stream() { let (w, h) = (96, 64); let (y, u, v) = source(w, h, 0); - let mut e = Vp9Encoder::new(EncoderConfig { width: w as u32, height: h as u32, quality: 100, full_range: false }); + let mut e = Vp9Encoder::new(EncoderConfig { width: w as u32, height: h as u32, quality: 100, full_range: false, ..Default::default() }); let frame = e.encode_yuv420(&y, &u, &v); // IVF container (one frame). let mut ivf = Vec::new(); @@ -147,3 +147,176 @@ fn ffmpeg_decodes_our_stream() { assert_eq!(&out.stdout[..w * h], &y[..], "ffmpeg luma matches (lossless)"); assert_eq!(&out.stdout[w * h..w * h + u.len()], &u[..]); } + +// ---------------------------------------------------------------- inter frames + +/// Frame `t` of a moving sequence: a diagonal gradient with a sinusoidal texture and bright +/// squares, all translating by a fractional amount per frame (`kind` 1 adds noise that changes +/// every frame). +fn moving(w: usize, h: usize, t: usize, kind: u32) -> (Vec, Vec, Vec) { + let (cw, ch) = (w.div_ceil(2), h.div_ceil(2)); + let (dx, dy) = (1.75 * t as f64, 0.6 * t as f64); + let mut s = 1234u32 + t as u32 * 77; + let mut rnd = || { + s = s.wrapping_mul(1664525).wrapping_add(1013904223); + (s >> 28) as f64 + }; + let lum = |x: f64, y: f64| -> f64 { + let g = 40.0 + 120.0 * (x / w as f64) + 50.0 * (y / h as f64); + let tex = 18.0 * ((x * 0.31).sin() * (y * 0.23).cos()); + let sq = if (x - 20.0).rem_euclid(64.0) < 18.0 && (y - 10.0).rem_euclid(48.0) < 14.0 { 45.0 } else { 0.0 }; + g + tex + sq + }; + let mut y = vec![0u8; w * h]; + for yy in 0..h { + for x in 0..w { + let n = if kind == 1 { rnd() - 7.5 } else { 0.0 }; + y[yy * w + x] = (lum(x as f64 - dx, yy as f64 - dy) + n).round().clamp(0.0, 255.0) as u8; + } + } + let u: Vec = (0..cw * ch).map(|i| (100.0 + 40.0 * (((i % cw) as f64 * 2.0 - dx) / w as f64)).round() as u8).collect(); + let v: Vec = (0..cw * ch).map(|i| (150.0 - 30.0 * (((i / cw) as f64 * 2.0 - dy) / h as f64)).round() as u8).collect(); + (y, u, v) +} + +fn planes_of(p: &filmcraft_vp9::Picture, w: usize, h: usize) -> [Vec; 3] { + let (cw, ch) = (w.div_ceil(2), h.div_ceil(2)); + [crop(&plane_u8(&p.y), p.y_stride, w, h), crop(&plane_u8(&p.u), p.uv_stride, cw, ch), crop(&plane_u8(&p.v), p.uv_stride, cw, ch)] +} + +/// Encode `n` frames and decode them with one decoder: every decoded frame must equal the +/// encoder's reconstruction. Returns (total bytes, mean luma PSNR vs the source). +fn run_sequence(cfg: EncoderConfig, n: usize, kind: u32) -> (usize, f64) { + let (w, h) = (cfg.width as usize, cfg.height as usize); + let mut e = Vp9Encoder::new(cfg); + let mut d = Decoder::new(); + let (mut bytes, mut db) = (0usize, 0.0); + for t in 0..n { + let (y, u, v) = moving(w, h, t, kind); + let f = e.encode_frame(&y, &u, &v, false); + bytes += f.data.len(); + let kind_s = if f.key { "key" } else { "inter" }; + let mut pics = d.decode(&f.data, t as i64).unwrap_or_else(|err| panic!("{w}x{h} frame {t} ({kind_s}): {err:?}")); + assert_eq!(pics.len(), 1); + let got = planes_of(&pics.remove(0), w, h); + let rec = e.reconstruction().expect("reconstruction"); + for k in 0..3 { + if got[k] != rec[k] { + let first = got[k].iter().zip(&rec[k]).position(|(a, b)| a != b).unwrap_or(0); + let pw = if k == 0 { w } else { w.div_ceil(2) }; + panic!("{w}x{h} frame {t} ({kind_s}, q{} lf{}): plane {k} differs first at ({}, {})", f.q_idx, f.filter_level, first % pw, first / pw); + } + } + db += psnr(&got[0], &y); + } + (bytes, db / n as f64) +} + +#[test] +fn inter_frames_decode_exactly() { + for (w, h) in [(64, 64), (96, 72), (33, 17), (130, 66)] { + for (q, kind) in [(30u8, 0u32), (70, 1), (92, 0), (100, 0)] { + let cfg = EncoderConfig { width: w as u32, height: h as u32, quality: q, keyframe_interval: 6, ..Default::default() }; + run_sequence(cfg, 8, kind); + } + } +} + +#[test] +fn inter_frames_with_tile_columns_decode_exactly() { + let cfg = EncoderConfig { width: 4160, height: 24, quality: 80, keyframe_interval: 0, ..Default::default() }; + run_sequence(cfg, 3, 0); +} + +#[test] +fn inter_psnr_rises_with_bitrate() { + let (w, h) = (160, 96); + let mut last = (0usize, 0.0f64); + let mut rows = vec![]; + for kbps in [60u32, 150, 400, 1200] { + let cfg = EncoderConfig { width: w, height: h, quality: 60, target_kbps: Some(kbps), frame_rate: 30.0, keyframe_interval: 30, ..Default::default() }; + let (bytes, db) = run_sequence(cfg, 20, 0); + rows.push(format!("{kbps} kbps target: {:.0} kbps actual, {db:.2} dB", bytes as f64 * 8.0 * 30.0 / 20.0 / 1000.0)); + assert!(bytes > last.0 && db > last.1, "{rows:#?}"); + last = (bytes, db); + } + eprintln!("{rows:#?}"); +} + +#[test] +fn inter_frames_are_much_smaller_than_intra() { + let (w, h) = (192u32, 108u32); + let mut report = vec![]; + for q in [50u8, 75, 90] { + let intra = EncoderConfig { width: w, height: h, quality: q, keyframe_interval: 1, ..Default::default() }; + let inter = EncoderConfig { keyframe_interval: 0, ..intra.clone() }; + let (bi, pi) = run_sequence(intra, 24, 0); + let (bp, pp) = run_sequence(inter, 24, 0); + report.push(format!("q{q}: intra {bi} B {pi:.2} dB, inter {bp} B {pp:.2} dB, ratio {:.2}", bi as f64 / bp as f64)); + assert!(pp > pi - 1.0, "{report:#?}"); + assert!(bi as f64 >= 3.0 * bp as f64, "{report:#?}"); + } + eprintln!("{report:#?}"); +} + +#[test] +fn ffmpeg_decodes_inter_frames_like_the_encoder() { + let (w, h) = (96usize, 64usize); + let mut e = Vp9Encoder::new(EncoderConfig { width: w as u32, height: h as u32, quality: 75, keyframe_interval: 5, ..Default::default() }); + let n = 7; + let mut ivf = Vec::new(); + ivf.extend_from_slice(b"DKIF"); + ivf.extend_from_slice(&0u16.to_le_bytes()); + ivf.extend_from_slice(&32u16.to_le_bytes()); + ivf.extend_from_slice(b"VP90"); + ivf.extend_from_slice(&(w as u16).to_le_bytes()); + ivf.extend_from_slice(&(h as u16).to_le_bytes()); + ivf.extend_from_slice(&30u32.to_le_bytes()); + ivf.extend_from_slice(&1u32.to_le_bytes()); + ivf.extend_from_slice(&(n as u32).to_le_bytes()); + ivf.extend_from_slice(&0u32.to_le_bytes()); + let mut recon = vec![]; + for t in 0..n { + let (y, u, v) = moving(w, h, t, 0); + let f = e.encode_frame(&y, &u, &v, false); + ivf.extend_from_slice(&(f.data.len() as u32).to_le_bytes()); + ivf.extend_from_slice(&(t as u64).to_le_bytes()); + ivf.extend_from_slice(&f.data); + recon.push(e.reconstruction().expect("recon")); + } + let dir = std::env::temp_dir().join(format!("effectcraft-vp9-inter-{}", std::process::id())); + std::fs::create_dir_all(&dir).unwrap(); + let path = dir.join("a.ivf"); + std::fs::write(&path, &ivf).unwrap(); + let out = std::process::Command::new("ffmpeg") + .args(["-hide_banner", "-loglevel", "error", "-i"]) + .arg(&path) + .args(["-f", "rawvideo", "-pix_fmt", "yuv420p", "-"]) + .output(); + let Ok(out) = out else { + eprintln!("ffmpeg not found: skipping oracle"); + return; + }; + assert!(out.status.success(), "{}", String::from_utf8_lossy(&out.stderr)); + let fsz = w * h + 2 * w.div_ceil(2) * h.div_ceil(2); + assert_eq!(out.stdout.len(), fsz * n); + for (t, r) in recon.iter().enumerate() { + let frame = &out.stdout[t * fsz..(t + 1) * fsz]; + assert_eq!(&frame[..w * h], &r[0][..], "ffmpeg luma of frame {t}"); + assert_eq!(&frame[w * h..w * h + r[1].len()], &r[1][..], "ffmpeg U of frame {t}"); + } +} + +/// Encoder speed on 720p (run with `--ignored --nocapture`). +#[test] +#[ignore] +fn bench_720p() { + let (w, h) = (1280usize, 720usize); + let mut e = Vp9Encoder::new(EncoderConfig { width: w as u32, height: h as u32, quality: 80, keyframe_interval: 0, ..Default::default() }); + for t in 0..6 { + let (y, u, v) = moving(w, h, t, 0); + let t0 = std::time::Instant::now(); + let f = e.encode_frame(&y, &u, &v, false); + eprintln!("frame {t} ({}): {} B in {:.0} ms", if f.key { "key" } else { "inter" }, f.data.len(), t0.elapsed().as_secs_f64() * 1000.0); + } +} diff --git a/crates/vp9enc/Cargo.toml b/crates/vp9enc/Cargo.toml index 62144b25..fe2e767c 100644 --- a/crates/vp9enc/Cargo.toml +++ b/crates/vp9enc/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "effectcraft-vp9enc" -description = "VP9 intra-frame encoder (profile 0, 8-bit 4:2:0) for WebM export, from the VP9 Bitstream & Decoding Process Specification" +description = "VP9 encoder (profile 0, 8-bit 4:2:0; intra + inter frames, loop filter, rate control) for WebM export, from the VP9 Bitstream & Decoding Process Specification" version.workspace = true edition.workspace = true license.workspace = true diff --git a/crates/vp9enc/README.md b/crates/vp9enc/README.md index c8f781b8..b9f94fae 100644 --- a/crates/vp9enc/README.md +++ b/crates/vp9enc/README.md @@ -1,15 +1,47 @@ # effectcraft-vp9enc -A small VP9 encoder for WebM export: every frame is a key frame (profile 0, 8-bit 4:2:0, -BT.709), coded with 8×8 intra blocks (DC / V / H / TM prediction chosen per block), 4×4 -DCT/ADST transforms with default probabilities, or the Walsh–Hadamard transform at quality 100 -(mathematically lossless). +A small VP9 encoder for WebM export (profile 0, 8-bit 4:2:0, BT.709). + +* **Key frames**: 8×8 intra blocks (DC / V / H / TM prediction chosen per block), 4×4 DCT/ADST + transforms, or the Walsh–Hadamard transform at quality 100 (mathematically lossless). +* **Inter frames**: one reference (the previous frame, LAST), 16×16 or 8×8 blocks chosen by a + cost comparison; motion search from the predicted vectors (NearestMv / NearMv / last frame's + co-located vector) with a whole-sample cross/diamond search refined to half and quarter + samples with the regular 8-tap interpolation filters; NEARESTMV / NEARMV / ZEROMV / NEWMV + chosen by bit cost; blocks with no coefficients worth coding are skipped; 4×4 DCT residual + with a rate–distortion check per transform block. Inter frames are coded error resilient + (`error_resilient_mode` = 1: no dependence on earlier probability contexts or motion vectors), + so `UsePrevFrameMvs` is never needed. No high-precision (1/8) vectors, no compound prediction, + no intra blocks inside inter frames. +* **Loop filter**: the in-loop deblocking filter (§8.8) on every lossy frame; the level is the + best of three trials against the source (sharpness 0, no mode/reference deltas). +* **Probabilities**: per-frame forward updates in the compressed header (`diff_update_prob` / + `update_mv_prob`) of the coefficient, skip, inter-mode, is-inter, reference, partition and + motion-vector probabilities. Tiles are recorded as symbols first, the frame's own counts pick + the updates that pay for their signalling, then the symbols are entropy coded with the updated + probabilities. There is no backward adaptation (`refresh_frame_context` = 0). +* **Rate control**: constant quality (quality → `base_q_idx`) or a target bitrate: a per-frame + q index from a bits ≈ complexity / step model per frame type, key frames getting a larger share + of each key-frame interval, with the running over/undershoot paid back over about a second. +* **Key frame interval** is configurable (1 = every frame a key frame, the old intra-only mode). + +The encoder reconstructs exactly what a decoder reconstructs (inter prediction with the +decoder's vector clamping and edge-clamped reference fetches, integer inverse transforms, the +loop filter) and exposes it (`Vp9Encoder::reconstruction`); the tests require FilmCraft's VP9 +decoder (and ffmpeg, when installed) to output exactly that reconstruction, frame by frame, for +moving synthetic sequences at several sizes and qualities, with tile columns, and lossless. + +On a moving-gradient sequence (192×108, 24 frames) inter coding is 3.9–7× smaller than all-intra +at the same quality setting with equal or better PSNR. Written from the **VP9 Bitstream & Decoding Process Specification, version 0.6 (31 March 2016)**: the encoder emits the uncompressed and compressed headers (§6.2, §6.3), partition, mode -info and token syntax (§6.4) and reconstructs with the specification's intra prediction (§8.5.1) -and integer inverse transforms (§8.7) so it stays in step with any conformant decoder. The -boolean encoder is the inverse of the §9.2 decoder in the form given in RFC 6386 §7.3. Spec -tables (§10) were transcribed from FilmCraft's first-party decoder `filmcraft-vp9` (MIT OR -Apache-2.0), which is also the decoder the tests round-trip through; ffmpeg is used only as an -external test oracle. No libvpx or other third-party encoder/decoder source was consulted. +info, motion vector and token syntax (§6.4), derives motion vector candidates with the §6.4 +`find_mv_refs` / `find_best_ref_mvs` processes, and reconstructs with the specification's intra +and inter prediction (§8.5), integer inverse transforms (§8.7) and loop filter (§8.8) so it stays +in step with any conformant decoder. The boolean encoder is the inverse of the §9.2 decoder in +the form given in RFC 6386 §7.3. Spec tables (§10) are the specification's; the key-frame tables +were transcribed from FilmCraft's first-party decoder `filmcraft-vp9` (MIT OR Apache-2.0), which +is also the decoder the tests round-trip through and whose behaviour (loop filter edge rules, +vector candidate order) the encoder was checked against; ffmpeg is used only as an external test +oracle. No libvpx or other third-party encoder/decoder source was consulted. diff --git a/crates/vp9enc/src/bool.rs b/crates/vp9enc/src/bool.rs index 45db92dd..0cd91b0e 100644 --- a/crates/vp9enc/src/bool.rs +++ b/crates/vp9enc/src/bool.rs @@ -64,31 +64,6 @@ impl BoolEncoder { } } - /// Write `value` with a tree (leaves stored negated, as in the decoder's tables). - pub fn tree(&mut self, tree: &[i8], probs: &[u8], value: u8) { - // Find the path to the leaf. - fn path(tree: &[i8], node: usize, value: u8, out: &mut Vec<(usize, bool)>) -> bool { - for b in 0..2 { - let t = tree[node + b]; - out.push((node >> 1, b == 1)); - if t <= 0 { - if (-t) as u8 == value { - return true; - } - } else if path(tree, t as usize, value, out) { - return true; - } - out.pop(); - } - false - } - let mut p = Vec::with_capacity(8); - path(tree, 0, value, &mut p); - for (i, bit) in p { - self.write(bit, probs[i]); - } - } - /// Flush (RFC 6386 §7.3) and return the bytes. pub fn finish(mut self) -> Vec { let mut c = self.bit_count; diff --git a/crates/vp9enc/src/lf.rs b/crates/vp9enc/src/lf.rs new file mode 100644 index 00000000..d781144f --- /dev/null +++ b/crates/vp9enc/src/lf.rs @@ -0,0 +1,124 @@ +//! The in-loop deblocking filter (spec 8.8), as the decoder applies it to a frame coded with 4×4 +//! transforms, 8×8 / 16×16 blocks, no segmentation and no reference/mode deltas: superblocks in +//! raster order, for each plane the vertical edges (pass 0) then the horizontal edges (pass 1). +//! Transform edges inside a skipped inter block are left alone; every 4×4 edge uses the 4-tap +//! filter except edges on a 32-sample boundary, which get the 8-tap (flat) filter. + +use crate::{BLOCK_16X16, Mi, Plane}; + +#[derive(Clone, Copy)] +struct Limits { + limit: i32, + blimit: i32, + thresh: i32, +} + +/// Filter the reconstructed planes (Y, U, V, or a prefix of them) in place at `level` (0 = off; +/// sharpness 0). +pub fn filter_frame(planes: &mut [Plane], mi: &[Mi], mi_rows: usize, mi_cols: usize, level: u8) { + if level == 0 { + return; + } + let l = level as i32; + let limit = l.max(1); + let lim = Limits { limit, blimit: 2 * (l + 2) + limit, thresh: l >> 4 }; + for sbr in 0..mi_rows.div_ceil(8) { + for sbc in 0..mi_cols.div_ceil(8) { + for (plane, p) in planes.iter_mut().enumerate() { + for pass in 0..2 { + sb_plane(p, mi, mi_rows, mi_cols, plane, pass, sbr * 8, sbc * 8, lim); + } + } + } + } +} + +#[allow(clippy::too_many_arguments)] +fn sb_plane(p: &mut Plane, mi: &[Mi], mi_rows: usize, mi_cols: usize, plane: usize, pass: usize, row: usize, col: usize, lim: Limits) { + let (sx, sy) = if plane > 0 { (1, 1) } else { (0, 0) }; + let (sub, edge_len) = if pass == 0 { (sx, 64 >> sy) } else { (sy, 64 >> sx) }; + let stride = p.w; + for edge in 0..(16usize >> sub) { + let mut i = 0; + while i < edge_len { + let (x, y) = if pass == 0 { (col * 8 + edge * (4 << sx), row * 8 + (i << sy)) } else { (col * 8 + (i << sx), row * 8 + edge * (4 << sy)) }; + if x >= 8 * mi_cols || y >= 8 * mi_rows || (pass == 0 && x == 0) || (pass == 1 && y == 0) { + i += 8; + continue; + } + let lc = ((x >> 3) >> sx) << sx; + let lr = ((y >> 3) >> sy) << sy; + let m = mi[lr * mi_cols + lc]; + let size = if sub == 0 { m.size } else { m.size.max(BLOCK_16X16) }; + let n8 = if size >= BLOCK_16X16 { 2 } else { 1 }; + let is_block_edge = if pass == 0 { x & (8 * n8 - 1) == 0 } else { y & (8 * n8 - 1) == 0 }; + let tx_ok = !m.inter || !m.skip; + if !is_block_edge && !tx_ok { + i += 8; + continue; + } + let wide = edge % 8 == 0; + let n = if pass == 0 { (8 * mi_rows - y).div_ceil(1 << sy) } else { (8 * mi_cols - x).div_ceil(1 << sx) }.min(8).min(edge_len - i); + // A transform edge is not filtered in the last odd chroma column. + let odd_last = pass == 1 && sx == 1 && mi_cols & 1 == 1 && edge & 1 == 1; + let (px, py) = (x >> sx, y >> sy); + for l in 0..n { + let is_tx_edge = tx_ok && !(odd_last && x + (l << sx) + 8 >= mi_cols * 8); + if !(is_block_edge || is_tx_edge) { + continue; + } + let (pos, step) = if pass == 0 { ((py + l) * stride + px, 1) } else { (py * stride + px + l, stride) }; + filter_lane(&mut p.rec, pos, step, wide, lim); + } + i += 8; + } + } +} + +/// Filter mask, narrow (4-tap) and 8-tap flat filters of spec 8.8.5 across one edge position; +/// `pos` is q0, samples across the edge are `step` apart. +fn filter_lane(d: &mut [u8], pos: usize, step: usize, wide: bool, lim: Limits) { + let at = |k: isize| d[(pos as isize + k * step as isize) as usize] as i32; + let (p3, p2, p1, p0) = (at(-4), at(-3), at(-2), at(-1)); + let (q0, q1, q2, q3) = (at(0), at(1), at(2), at(3)); + let Limits { limit, blimit, thresh } = lim; + let mask = (p3 - p2).abs() <= limit + && (p2 - p1).abs() <= limit + && (p1 - p0).abs() <= limit + && (q1 - q0).abs() <= limit + && (q2 - q1).abs() <= limit + && (q3 - q2).abs() <= limit + && (p0 - q0).abs() * 2 + (p1 - q1).abs() / 2 <= blimit; + if !mask { + return; + } + let flat = + wide && (p1 - p0).abs() <= 1 && (q1 - q0).abs() <= 1 && (p2 - p0).abs() <= 1 && (q2 - q0).abs() <= 1 && (p3 - p0).abs() <= 1 && (q3 - q0).abs() <= 1; + let mut set = |k: isize, v: i32| d[(pos as isize + k * step as isize) as usize] = v as u8; + if flat { + // 8-tap: positions p2..q2, each the rounded mean of a 7-wide clamped window plus itself. + let s = [p3, p2, p1, p0, q0, q1, q2, q3]; + let px = |k: isize| s[(k.clamp(-4, 3) + 4) as usize]; + for k in -3isize..3 { + let mut sum = px(k); + for j in -3..=3 { + sum += px(k + j); + } + set(k, (sum + 4) >> 3); + } + return; + } + let c = |v: i32| v.clamp(-128, 127); + let hev = (p1 - p0).abs() > thresh || (q1 - q0).abs() > thresh; + let (ps1, ps0, qs0, qs1) = (p1 - 128, p0 - 128, q0 - 128, q1 - 128); + let filter = c(if hev { c(ps1 - qs1) } else { 0 } + 3 * (qs0 - ps0)); + let f1 = c(filter + 4) >> 3; + let f2 = c(filter + 3) >> 3; + set(0, c(qs0 - f1) + 128); + set(-1, c(ps0 + f2) + 128); + if !hev { + let f = (f1 + 1) >> 1; + set(1, c(qs1 - f) + 128); + set(-2, c(ps1 + f) + 128); + } +} diff --git a/crates/vp9enc/src/lib.rs b/crates/vp9enc/src/lib.rs index 02cbb0e2..8f72cb7b 100644 --- a/crates/vp9enc/src/lib.rs +++ b/crates/vp9enc/src/lib.rs @@ -1,290 +1,89 @@ -//! A VP9 intra-frame encoder (profile 0, 8-bit 4:2:0), written from the *VP9 Bitstream & -//! Decoding Process Specification* v0.6 by mirroring its decoding process in reverse (see -//! README.md). Every frame is a key frame: one tile row (tile columns as the frame width -//! requires), 64×64 superblocks split to 8×8 blocks, 4×4 transforms (the Walsh–Hadamard -//! transform at quality 100, which is lossless), intra modes DC / V / H / TM chosen per block by -//! prediction error, default probabilities (no probability updates). The encoder reconstructs -//! exactly what a decoder reconstructs (the specification's integer inverse transforms and -//! intra edge rules), so prediction never drifts. +//! A VP9 encoder (profile 0, 8-bit 4:2:0), written from the *VP9 Bitstream & Decoding Process +//! Specification* v0.6 by mirroring its decoding process in reverse (see README.md). +//! +//! * **Key frames**: 64×64 superblocks split to 8×8 intra blocks (DC / V / H / TM chosen per +//! block), 4×4 DCT/ADST transforms (the Walsh–Hadamard transform at quality 100: lossless). +//! * **Inter frames** (error resilient, single reference: the previous frame as LAST): 16×16 or +//! 8×8 blocks with motion vectors found by a candidate + diamond search refined to quarter +//! samples with the regular 8-tap filters; NEARESTMV / NEARMV / ZEROMV / NEWMV chosen by cost, +//! motion vector prediction exactly as the decoder derives it; skipped blocks; 4×4 DCT residual. +//! * **In-loop deblocking** (8.8) on every lossy frame; the level is picked per frame by trial. +//! * **Probabilities**: per-frame forward updates (compressed header) of the coefficient, skip, +//! inter mode, reference, partition and motion vector probabilities, chosen from the frame's +//! own symbol counts (no backward adaptation: `refresh_frame_context` = 0). +//! * **Rate control**: constant quality (quality → base_q_idx), or a target bitrate with a +//! per-frame q_idx model. +//! +//! The encoder reconstructs exactly what a decoder reconstructs (prediction, integer inverse +//! transforms, loop filter), so prediction never drifts; [`Vp9Encoder::reconstruction`] returns +//! it. mod bool; +mod lf; +mod mc; +mod probs; mod tables; +mod tile; +mod transform; -use bool::BoolEncoder; -use tables::*; +pub use mc::Mv; +use probs::Probs; +use tables::{AC_QLOOKUP, DC_QLOOKUP}; +use transform::Quant; -const DC_PRED: u8 = 0; -const V_PRED: u8 = 1; -const H_PRED: u8 = 2; -const TM_PRED: u8 = 9; -const MODES: [u8; 4] = [DC_PRED, V_PRED, H_PRED, TM_PRED]; - -const DCT_DCT: u8 = 0; -const ADST_DCT: u8 = 1; -const DCT_ADST: u8 = 2; -const ADST_ADST: u8 = 3; - -/// mode2txfm_map for the modes used here. -fn mode_tx_type(mode: u8) -> u8 { - match mode { - V_PRED => ADST_DCT, - H_PRED => DCT_ADST, - TM_PRED => ADST_ADST, - _ => DCT_DCT, - } -} - -const PARTITION_TREE: [i8; 6] = [0, 2, -1, 4, -2, -3]; -const INTRA_MODE_TREE: [i8; 18] = [0, 2, -9, 4, -1, 6, 8, 12, -2, 10, -4, -5, -3, 14, -8, 16, -6, -7]; - -/// extra_bits[token] = (cat, numExtra, base). -const EXTRA_BITS: [(u8, u8, u16); 11] = - [(0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 0, 3), (0, 0, 4), (1, 1, 5), (2, 2, 7), (3, 3, 11), (4, 4, 19), (5, 5, 35), (6, 14, 67)]; -const CAT_PROBS: [&[u8]; 7] = [ - &[], - &[159], - &[165, 145], - &[173, 148, 140], - &[176, 155, 140, 135], - &[180, 157, 141, 134, 130], - &[254, 254, 254, 252, 249, 243, 230, 196, 177, 153, 140, 133, 130, 129], -]; +pub(crate) const BLOCK_8X8: u8 = 3; +pub(crate) const BLOCK_16X16: u8 = 6; /// Encoder settings. -#[derive(Clone, Debug, PartialEq, Eq)] +#[derive(Clone, Debug, PartialEq)] pub struct EncoderConfig { pub width: u32, pub height: u32, - /// 1..=100; 100 is lossless. + /// 1..=100; 100 is lossless. With a target bitrate: the quality of the first frame. pub quality: u8, /// Full-range (0–255) rather than limited (16–235) YUV; only signalled. pub full_range: bool, + /// A key frame every this many frames (1 = every frame is a key frame; 0 = only the first). + pub keyframe_interval: u32, + /// Target bitrate in kbit/s (`None`: constant quality). + pub target_kbps: Option, + /// Frames per second (rate control). + pub frame_rate: f64, + /// The in-loop deblocking filter. + pub loop_filter: bool, } -/// Probabilities of the full token tree (9.3.2) for model probability `p`. -fn pareto(p: u8) -> [u8; 8] { - let p = p.max(1) as usize; - let x = (p - 1) / 2; - let mut t = [0u8; 8]; - for n in 0..8 { - t[n] = if p & 1 == 1 { PARETO_TABLE[x * 8 + n] } else { ((PARETO_TABLE[x * 8 + n] as u32 + PARETO_TABLE[(x + 1) * 8 + n] as u32) >> 1) as u8 }; - } - t -} - -fn scan(tx_type: u8) -> &'static [u16; 16] { - match tx_type { - ADST_DCT => &ROW_SCAN_4X4, - DCT_ADST => &COL_SCAN_4X4, - _ => &DEFAULT_SCAN_4X4, +impl Default for EncoderConfig { + fn default() -> Self { + EncoderConfig { width: 0, height: 0, quality: 80, full_range: false, keyframe_interval: 60, target_kbps: None, frame_rate: 30.0, loop_filter: true } } } -/// Coefficient context neighbours (9.3.2) of each scan position. -fn neighbors(tx_type: u8) -> [(usize, usize); 16] { - let sc = scan(tx_type); - let mut nb = [(0, 0); 16]; - for (c, &pos) in sc.iter().enumerate().skip(1) { - let pos = pos as usize; - let (i, j) = (pos / 4, pos % 4); - nb[c] = if i > 0 && j > 0 { - let a = (i - 1) * 4 + j; - let b = i * 4 + j - 1; - match tx_type { - DCT_ADST => (a, a), - ADST_DCT => (b, b), - _ => (a, b), - } - } else if i > 0 { - ((i - 1) * 4 + j, (i - 1) * 4 + j) - } else { - (j - 1, j - 1) - }; - } - nb +/// One coded frame. +#[derive(Clone, Debug)] +pub struct EncodedFrame { + pub data: Vec, + pub key: bool, + pub q_idx: u8, + pub filter_level: u8, } -// ---------------------------------------------------------------- transforms (8.7) - -fn r14(v: i32) -> i32 { - (v + (1 << 13)) >> 14 -} - -fn idct4(x: [i32; 4]) -> [i32; 4] { - let v0 = r14(x[0] * 11585 - x[2] * 11585); - let v1 = r14(x[0] * 11585 + x[2] * 11585); - let v2 = r14(x[1] * 6270 - x[3] * 15137); - let v3 = r14(x[1] * 15137 + x[3] * 6270); - [v1 + v3, v0 + v2, v0 - v2, v1 - v3] -} - -fn iadst4(x: [i32; 4]) -> [i32; 4] { - let (s19, s29, s39, s49) = (5283i64, 9929i64, 13377i64, 15212i64); - let x: [i64; 4] = x.map(|v| v as i64); - let s0 = s19 * x[0]; - let s1 = s29 * x[0]; - let s2 = s39 * x[1]; - let s3 = s49 * x[2]; - let s4 = s19 * x[2]; - let s5 = s29 * x[3]; - let s6 = s49 * x[3]; - let s7 = s39 * (x[0] - x[2] + x[3]); - let x0 = s0 + s3 + s5; - let x1 = s1 - s4 - s6; - let x2 = s7; - let x3 = s2; - let r = |v: i64| ((v + (1 << 13)) >> 14) as i32; - [r(x0 + x3), r(x1 + x3), r(x2), r(x0 + x1 - x3)] -} - -/// The decoder's 4×4 inverse transform of dequantised coefficients (row-major) to a residual. -fn inverse4(coefs: &[i32; 16], tx_type: u8) -> [i32; 16] { - let row_adst = matches!(tx_type, DCT_ADST | ADST_ADST); - let col_adst = matches!(tx_type, ADST_DCT | ADST_ADST); - let mut t = [[0i32; 4]; 4]; - for i in 0..4 { - let r = [coefs[i * 4], coefs[i * 4 + 1], coefs[i * 4 + 2], coefs[i * 4 + 3]]; - t[i] = if row_adst { iadst4(r) } else { idct4(r) }; - } - let mut out = [0i32; 16]; - for j in 0..4 { - let c = [t[0][j], t[1][j], t[2][j], t[3][j]]; - let c = if col_adst { iadst4(c) } else { idct4(c) }; - for i in 0..4 { - out[i * 4 + j] = (c[i] + 8) >> 4; - } - } - out -} - -/// Inverse Walsh–Hadamard (8.7.1.10) on one line. -fn iwht(x: [i32; 4], shift: u32) -> [i32; 4] { - let mut a = x[0] >> shift; - let mut c = x[1] >> shift; - let mut d = x[2] >> shift; - let mut b = x[3] >> shift; - a += c; - d -= b; - let e = (a - d) >> 1; - b = e - b; - c = e - c; - a -= b; - d += c; - [a, b, c, d] -} - -/// The exact inverse of [`iwht`] (with shift 0): inputs that produce outputs `y`. -fn fwht(y: [i32; 4]) -> [i32; 4] { - let (aa, bb, cc, dd) = (y[0], y[1], y[2], y[3]); - let d1 = dd - cc; - let a1 = aa + bb; - let e = (a1 - d1) >> 1; - let b0 = e - bb; - let c0 = e - cc; - let a0 = a1 - c0; - let d0 = d1 + b0; - [a0, c0, d0, b0] -} - -fn inverse_wht(coefs: &[i32; 16]) -> [i32; 16] { - let mut t = [[0i32; 4]; 4]; - for i in 0..4 { - t[i] = iwht([coefs[i * 4], coefs[i * 4 + 1], coefs[i * 4 + 2], coefs[i * 4 + 3]], 2); - } - let mut out = [0i32; 16]; - for j in 0..4 { - let r = iwht([t[0][j], t[1][j], t[2][j], t[3][j]], 0); - for i in 0..4 { - out[i * 4 + j] = r[i]; - } - } - out -} - -/// Lossless forward transform: levels (dequantised by 4) whose inverse is exactly `res`. -fn forward_wht(res: &[i32; 16]) -> [i32; 16] { - let mut t = [[0i32; 4]; 4]; - for j in 0..4 { - let col = fwht([res[j], res[4 + j], res[8 + j], res[12 + j]]); - for i in 0..4 { - t[i][j] = col[i]; - } - } - let mut lv = [0i32; 16]; - for i in 0..4 { - let r = fwht(t[i]); - lv[i * 4..i * 4 + 4].copy_from_slice(&r); - } - lv -} - -/// Forward transform matrices: the inverses of the decoder's (linear part of the) inverse -/// transforms, so `inverse4(forward(res)) ≈ res`. -fn forward_matrices() -> &'static [[[f64; 16]; 16]; 4] { - use std::sync::OnceLock; - static M: OnceLock<[[[f64; 16]; 16]; 4]> = OnceLock::new(); - M.get_or_init(|| { - let mut out = [[[0.0; 16]; 16]; 4]; - for (ty, o) in out.iter_mut().enumerate() { - // Columns of the inverse: response to a scaled unit coefficient. - let k = 1 << 16; - let mut inv = [[0.0f64; 16]; 16]; - for c in 0..16 { - let mut e = [0i32; 16]; - e[c] = k; - let r = inverse4(&e, ty as u8); - for p in 0..16 { - inv[p][c] = r[p] as f64 / k as f64; - } - } - *o = invert16(inv); - } - out - }) -} - -fn invert16(m: [[f64; 16]; 16]) -> [[f64; 16]; 16] { - let mut a = m; - let mut inv = [[0.0; 16]; 16]; - for (i, r) in inv.iter_mut().enumerate() { - r[i] = 1.0; - } - for col in 0..16 { - let piv = (col..16).max_by(|&x, &y| a[x][col].abs().total_cmp(&a[y][col].abs())).unwrap_or(col); - a.swap(col, piv); - inv.swap(col, piv); - let d = a[col][col]; - for j in 0..16 { - a[col][j] /= d; - inv[col][j] /= d; - } - for r in 0..16 { - if r != col { - let f = a[r][col]; - if f != 0.0 { - for j in 0..16 { - a[r][j] -= f * a[col][j]; - inv[r][j] -= f * inv[col][j]; - } - } - } - } - } - inv -} - -// ---------------------------------------------------------------- planes and prediction - -struct Plane { - w: usize, - /// Last valid column / row for prediction edges (the 8-aligned mode-info grid). - max_x: usize, - max_y: usize, - src: Vec, - rec: Vec, +/// A plane being coded: source and reconstruction at superblock-aligned size. +pub(crate) struct Plane { + /// Stride (superblock-aligned width). + pub w: usize, + /// Last valid column / row for intra prediction edges (the 8-aligned mode-info grid). + pub max_x: usize, + pub max_y: usize, + /// Visible size. + pub vis_w: usize, + pub vis_h: usize, + pub src: Vec, + pub rec: Vec, } impl Plane { + #[allow(clippy::too_many_arguments)] fn new(data: &[u8], vw: usize, vh: usize, w: usize, h: usize, max_x: usize, max_y: usize) -> Plane { let mut src = vec![0u8; w * h]; for y in 0..h { @@ -294,103 +93,141 @@ impl Plane { src[y * w + x] = data.get(sy * vw + sx).copied().unwrap_or(128); } } - Plane { w, max_x, max_y, src, rec: vec![0; w * h] } + Plane { w, max_x, max_y, vis_w: vw, vis_h: vh, src, rec: vec![0; w * h] } } } -/// Intra prediction (8.5.1) of a square block for DC / V / H / TM. -fn predict(p: &Plane, x: usize, y: usize, log2: usize, mode: u8, have_left: bool, have_above: bool) -> [i32; 64] { - let size = 1usize << log2; - let base = 128i32; - let buf = &p.rec; - let w = p.w; - let mut above = [0i32; 8]; - let mut left = [0i32; 8]; - if have_above { - for i in 0..size { - above[i] = buf[(y - 1) * w + (x + i).min(p.max_x)] as i32; - } - } else { - above[..size].fill(base - 1); - } - let top_left = if have_above && have_left { - buf[(y - 1) * w + (x - 1).min(p.max_x)] as i32 - } else if have_above { - base + 1 - } else { - base - 1 - }; - if have_left { - for i in 0..size { - left[i] = buf[(y + i).min(p.max_y) * w + x - 1] as i32; - } - } else { - left[..size].fill(base + 1); - } - let mut out = [0i32; 64]; - for i in 0..size { - for j in 0..size { - out[i * size + j] = match mode { - V_PRED => above[j], - H_PRED => left[i], - TM_PRED => (above[j] + left[i] - top_left).clamp(0, 255), - _ => 0, - }; - } - } - if mode == DC_PRED { - let v = match (have_left, have_above) { - (true, true) => (left[..size].iter().sum::() + above[..size].iter().sum::() + size as i32) >> (log2 + 1), - (true, false) => (left[..size].iter().sum::() + (1 << (log2 - 1))) >> log2, - (false, true) => (above[..size].iter().sum::() + (1 << (log2 - 1))) >> log2, - (false, false) => base, - }; - out[..size * size].fill(v); - } - out +/// Mode info of one 8×8 position. +#[derive(Clone, Copy, Debug, Default)] +pub(crate) struct Mi { + pub inter: bool, + /// Intra mode (0..=9) or inter mode (NEARESTMV..=NEWMV). + pub mode: u8, + pub mv: Mv, + /// Skip as the loop filter sees it (inter blocks without coefficients). + pub skip: bool, + pub size: u8, } -// ---------------------------------------------------------------- the encoder +/// Rate control state (target bitrate mode). +#[derive(Clone, Debug)] +struct RateControl { + /// bits ≈ complexity / step: per frame type, from the last frame of that type. + key_complexity: Option, + inter_complexity: Option, + /// Bits spent minus bits budgeted so far. + debt: f64, + last_inter_q: Option, +} pub struct Vp9Encoder { cfg: EncoderConfig, - qidx: u8, + base_qidx: u8, lossless: bool, + frames: u64, + since_key: u32, + /// The reconstruction of the last frame (the LAST reference): planes at their strides. + reference: Option<[Plane; 3]>, + mvs: Vec, + rc: RateControl, } -/// Per-tile coding state. -struct Tile<'a> { - enc: &'a Vp9Encoder, - bc: BoolEncoder, - planes: &'a mut [Plane; 3], - mi_cols: usize, - mi_rows: usize, - col_start: usize, - above_nz: [Vec; 3], - left_nz: [Vec; 3], - above_part: Vec, - left_part: Vec, - /// y mode of every 8×8 block (for mode contexts). - modes: &'a mut [u8], +fn quality_to_qidx(q: u8) -> u8 { + // 99 ≈ 3, 90 ≈ 30, 50 ≈ 150, 1 ≈ 255. + ((100 - q.clamp(1, 99) as i32) * 255 / 85).clamp(1, 255) as u8 } impl Vp9Encoder { pub fn new(cfg: EncoderConfig) -> Self { let q = cfg.quality.clamp(1, 100); - let lossless = q >= 100; - // Quality → base_q_idx: 99 ≈ 3, 90 ≈ 30, 50 ≈ 150, 1 ≈ 255. - let qidx = if lossless { 0 } else { ((100 - q as i32) * 255 / 85).clamp(1, 255) as u8 }; - Vp9Encoder { cfg, qidx, lossless } + let lossless = q >= 100 && cfg.target_kbps.is_none(); + let base_qidx = if lossless { 0 } else { quality_to_qidx(q) }; + Vp9Encoder { + cfg, + base_qidx, + lossless, + frames: 0, + since_key: 0, + reference: None, + mvs: vec![], + rc: RateControl { key_complexity: None, inter_complexity: None, debt: 0.0, last_inter_q: None }, + } } pub fn config(&self) -> &EncoderConfig { &self.cfg } - /// Encode one key frame from planar 8-bit 4:2:0 (`y`: width×height, `u`/`v`: + /// Encode one frame from planar 8-bit 4:2:0 (`y`: width×height, `u`/`v`: /// ((w+1)/2)×((h+1)/2), tightly packed). pub fn encode_yuv420(&mut self, y: &[u8], u: &[u8], v: &[u8]) -> Vec { + self.encode_frame(y, u, v, false).data + } + + /// The decoder's reconstruction of the last coded frame (visible area, tightly packed Y, U, + /// V). + pub fn reconstruction(&self) -> Option<[Vec; 3]> { + let r = self.reference.as_ref()?; + Some(std::array::from_fn(|i| { + let p = &r[i]; + (0..p.vis_h).flat_map(|y| p.rec[y * p.w..y * p.w + p.vis_w].iter().copied()).collect() + })) + } + + fn is_key_due(&self, force: bool) -> bool { + force || self.reference.is_none() || self.cfg.keyframe_interval == 1 || (self.cfg.keyframe_interval > 1 && self.since_key >= self.cfg.keyframe_interval) + } + + /// The q index of the next frame. + fn pick_qidx(&mut self, key: bool) -> u8 { + if self.lossless { + return 0; + } + let Some(kbps) = self.cfg.target_kbps else { return self.base_qidx }; + let per_frame = kbps as f64 * 1000.0 / self.cfg.frame_rate.max(1.0); + let interval = if self.cfg.keyframe_interval == 0 { 300.0 } else { self.cfg.keyframe_interval as f64 }; + // Key frames get a larger share of the group's budget. + let boost = (interval / 6.0).clamp(1.0, 5.0); + let group = per_frame * interval; + let key_share = group * boost / (boost + interval - 1.0); + let inter_share = (group - key_share) / (interval - 1.0).max(1.0); + let mut target = if key { key_share } else { inter_share }; + // Pay back (or spend) the running difference over about a second. + let horizon = self.cfg.frame_rate.max(1.0); + target = (target - self.rc.debt / horizon).max(target * 0.25); + let complexity = if key { self.rc.key_complexity } else { self.rc.inter_complexity.or(self.rc.key_complexity.map(|c| c / (boost * 2.0))) }; + let Some(cx) = complexity else { return self.base_qidx }; + let step = (cx / target).max(4.0); + let mut q = (0..=255u8).min_by_key(|&q| ((AC_QLOOKUP[q as usize] as f64 - step).abs() * 16.0) as i64).unwrap_or(self.base_qidx); + if !key && let Some(last) = self.rc.last_inter_q { + q = q.clamp(last.saturating_sub(24), last.saturating_add(24)); + } + q.max(1) + } + + fn rc_update(&mut self, key: bool, qidx: u8, bytes: usize) { + let Some(kbps) = self.cfg.target_kbps else { return }; + let bits = bytes as f64 * 8.0; + let cx = bits * AC_QLOOKUP[qidx as usize] as f64; + if key { + self.rc.key_complexity = Some(cx); + } else { + self.rc.inter_complexity = Some(match self.rc.inter_complexity { + Some(c) => c * 0.5 + cx * 0.5, + None => cx, + }); + self.rc.last_inter_q = Some(qidx); + } + let per_frame = kbps as f64 * 1000.0 / self.cfg.frame_rate.max(1.0); + self.rc.debt += bits - per_frame; + } + + /// Encode one frame; `force_key` starts a new key frame. + pub fn encode_frame(&mut self, y: &[u8], u: &[u8], v: &[u8], force_key: bool) -> EncodedFrame { let (w, h) = (self.cfg.width.max(1) as usize, self.cfg.height.max(1) as usize); + let key = self.is_key_due(force_key); + let qidx = self.pick_qidx(key); + let lossless = qidx == 0; let mi_cols = w.div_ceil(8); let mi_rows = h.div_ceil(8); let sb_cols = mi_cols.div_ceil(8); @@ -413,71 +250,77 @@ impl Vp9Encoder { } max_log2 -= 1; let tile_cols = 1usize << min_log2; - let mut modes = vec![0u8; mi_cols * mi_rows]; - let mut tiles = vec![]; + let quant = Quant { dcq: DC_QLOOKUP[qidx as usize], acq: AC_QLOOKUP[qidx as usize], lossless, round_dc: 0.5, round_ac: 0.38 }; + let step = AC_QLOOKUP[qidx as usize] as f64 / 8.0; + let lambda_sse = 0.35 * step * step; + let refs = if key { + None + } else { + self.reference.as_ref().map(|r| std::array::from_fn::<_, 3, _>(|i| mc::RefPlane { data: &r[i].rec, stride: r[i].w, w: r[i].vis_w, h: r[i].vis_h })) + }; + let prev = (!key && self.mvs.len() == mi_cols * mi_rows).then_some(&self.mvs[..]); + let phases = refs.as_ref().map(|r| mc::PhasePlanes::new(&r[0], 48)); + let fp = tile::FrameParams { key, quant, mi_cols, mi_rows, refs, phases, prev_mvs: prev, lambda_sad: lambda_sse.sqrt(), lambda_sse }; + let mut mi = vec![Mi::default(); mi_cols * mi_rows]; + let mut writers = vec![]; for t in 0..tile_cols { let start = (((t * sb_cols) >> min_log2) * 8).min(mi_cols); let end = ((((t + 1) * sb_cols) >> min_log2) * 8).min(mi_cols); - let mut tile = Tile { - enc: self, - bc: BoolEncoder::new(), - planes: &mut planes, - mi_cols, - mi_rows, - col_start: start, - above_nz: [vec![0; mi_cols * 2 + 16], vec![0; mi_cols * 2 + 16], vec![0; mi_cols * 2 + 16]], - left_nz: [vec![0; 16], vec![0; 16], vec![0; 16]], - above_part: vec![0; mi_cols + 8], - left_part: vec![0; 8], - modes: &mut modes, - }; - tile.bc.write(false, 128); - for r in (0..mi_rows).step_by(8) { - for p in 0..3 { - tile.left_nz[p].fill(0); - } - tile.left_part.fill(0); - for c in (start..end).step_by(8) { - tile.partition(r, c, 12); - } - } - tiles.push(tile.bc.finish()); + let mut tile = tile::Tile::new(&fp, &mut planes, &mut mi, start, end); + tile.code(); + writers.push(tile.w); } - // Compressed header. - let mut ch_enc = BoolEncoder::new(); - ch_enc.write(false, 128); - if !self.lossless { - ch_enc.literal(0, 2); // tx_mode ONLY_4X4 + // Forward probability updates from this frame's counts. + let defaults = Probs::defaults(); + let mut counts = vec![[0u32; 2]; probs::NPROBS]; + for wr in &writers { + wr.count(&mut counts); } - ch_enc.literal(0, 1); // no coefficient probability updates (TX_4X4) - for _ in 0..3 { - ch_enc.write(false, 252); // skip probabilities unchanged - } - let compressed = ch_enc.finish(); + let wanted = probs::plan_updates(&defaults, &counts, !key); + let (compressed, used) = probs::compressed_header(&defaults, &wanted, lossless, !key); + let tiles: Vec> = writers.iter().map(|wr| wr.emit(&used)).collect(); + drop(fp); + // Loop filter: the level with the least error against the source. + let level = if lossless || !self.cfg.loop_filter { 0 } else { choose_level(&mut planes, &mi, mi_rows, mi_cols, qidx, key) }; + lf::filter_frame(&mut planes, &mi, mi_rows, mi_cols, level); // Uncompressed header. let mut bw = BitWriter::default(); bw.put(2, 2); // frame_marker bw.put(0, 1); // profile low bit bw.put(0, 1); // profile high bit bw.put(0, 1); // show_existing_frame - bw.put(0, 1); // frame_type = KEY_FRAME + bw.put(!key as u32, 1); // frame_type bw.put(1, 1); // show_frame - bw.put(0, 1); // error_resilient_mode - bw.put(0x49, 8); - bw.put(0x83, 8); - bw.put(0x42, 8); - bw.put(2, 3); // color_space CS_BT_709 - bw.put(self.cfg.full_range as u32, 1); - bw.put(w as u32 - 1, 16); - bw.put(h as u32 - 1, 16); - bw.put(0, 1); // render_and_frame_size_different - bw.put(0, 1); // refresh_frame_context - bw.put(1, 1); // frame_parallel_decoding_mode + if key { + bw.put(0, 1); // error_resilient_mode + bw.put(0x49, 8); + bw.put(0x83, 8); + bw.put(0x42, 8); + bw.put(2, 3); // color_space CS_BT_709 + bw.put(self.cfg.full_range as u32, 1); + bw.put(w as u32 - 1, 16); + bw.put(h as u32 - 1, 16); + bw.put(0, 1); // render_and_frame_size_different + bw.put(0, 1); // refresh_frame_context + bw.put(1, 1); // frame_parallel_decoding_mode + } else { + bw.put(1, 1); // error_resilient_mode (no reliance on earlier contexts or vectors) + bw.put(0x01, 8); // refresh_frame_flags: slot 0 becomes the new LAST + for _ in 0..3 { + bw.put(0, 3); // ref_frame_idx: slot 0 + bw.put(0, 1); // sign bias + } + bw.put(1, 1); // found_ref: the size of LAST + bw.put(0, 1); // render_and_frame_size_different + bw.put(0, 1); // allow_high_precision_mv + bw.put(0, 1); // is_filter_switchable + bw.put(1, 2); // raw_interpolation_filter → EIGHTTAP (regular) + } bw.put(0, 2); // frame_context_idx - bw.put(0, 6); // loop_filter_level + bw.put(level as u32, 6); // loop_filter_level bw.put(0, 3); // sharpness bw.put(0, 1); // mode_ref_delta_enabled - bw.put(self.qidx as u32, 8); + bw.put(qidx as u32, 8); bw.put(0, 1); // delta_q_y_dc bw.put(0, 1); // delta_q_uv_dc bw.put(0, 1); // delta_q_uv_ac @@ -496,10 +339,46 @@ impl Vp9Encoder { } out.extend_from_slice(&t); } - out + // State for the next frame. + self.mvs = mi.iter().map(|m| if m.inter { m.mv } else { Mv::ZERO }).collect(); + self.reference = Some(planes); + self.frames += 1; + self.since_key = if key { 1 } else { self.since_key + 1 }; + self.rc_update(key, qidx, out.len()); + EncodedFrame { data: out, key, q_idx: qidx, filter_level: level } } } +/// Try a few loop filter levels on the luma reconstruction; the one closest to the source +/// (squared error over the visible area) wins (0 when filtering doesn't help). +fn choose_level(planes: &mut [Plane; 3], mi: &[Mi], mi_rows: usize, mi_cols: usize, qidx: u8, key: bool) -> u8 { + let guess = ((qidx as f64 * if key { 0.25 } else { 0.2 }).round() as u8).clamp(1, 63); + let sse = |p: &Plane| -> u64 { + let mut s = 0u64; + for y in 0..p.vis_h { + for x in 0..p.vis_w { + let d = p.src[y * p.w + x] as i64 - p.rec[y * p.w + x] as i64; + s += (d * d) as u64; + } + } + s + }; + let mut best = (sse(&planes[0]), 0u8); + let saved = planes[0].rec.clone(); + for level in [guess / 2, guess, (guess as u16 * 3 / 2).min(63) as u8] { + if level == 0 { + continue; + } + lf::filter_frame(&mut planes[..1], mi, mi_rows, mi_cols, level); + let s = sse(&planes[0]); + planes[0].rec.copy_from_slice(&saved); + if s < best.0 { + best = (s, level); + } + } + best.1 +} + #[derive(Default)] struct BitWriter { out: Vec, @@ -527,254 +406,6 @@ impl BitWriter { } } -const BLOCK_8X8: u8 = 3; -const NUM_8X8: [usize; 13] = [1, 1, 1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8]; -const MI_WIDTH_LOG2: [usize; 13] = [0, 0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3]; - -impl Tile<'_> { - fn partition(&mut self, r: usize, c: usize, bsize: u8) { - if r >= self.mi_rows || c >= self.mi_cols { - return; - } - let num8 = NUM_8X8[bsize as usize]; - let half = num8 >> 1; - let has_rows = (r + half) < self.mi_rows; - let has_cols = (c + half) < self.mi_cols; - let bsl = MI_WIDTH_LOG2[bsize as usize]; - let boffset = 3 - bsl; - let mut above = 0u8; - let mut left = 0u8; - for i in 0..num8 { - above |= self.above_part[c + i]; - left |= self.left_part[(r & 7) + i]; - } - let above = (above >> boffset) & 1; - let left = (left >> boffset) & 1; - let ctx = bsl * 4 + left as usize * 2 + above as usize; - let probs = &KF_PARTITION_PROBS[ctx * 3..ctx * 3 + 3]; - if bsize == BLOCK_8X8 { - // PARTITION_NONE. - self.bc.write(false, probs[0]); - self.block(r, c); - self.above_part[c] = 15 >> 1; - self.left_part[r & 7] = 15 >> 1; - return; - } - // PARTITION_SPLIT. - if has_rows && has_cols { - self.bc.tree(&PARTITION_TREE, probs, 3); - } else if has_cols { - self.bc.write(true, probs[1]); - } else if has_rows { - self.bc.write(true, probs[2]); - } - let sub = bsize - 3; - self.partition(r, c, sub); - self.partition(r, c + half, sub); - self.partition(r + half, c, sub); - self.partition(r + half, c + half, sub); - } - - fn block(&mut self, r: usize, c: usize) { - let avail_u = r > 0; - let avail_l = c > self.col_start; - // Mode decision on the 8×8 luma block and the 4×4 chroma blocks. - let y_mode = self.choose(0, c * 8, r * 8, 3, avail_l, avail_u); - let uv_mode = { - let mut best = (i64::MAX, DC_PRED); - for m in MODES { - let cost: i64 = (1..3).map(|p| self.cost(p, c * 4, r * 4, 2, m, avail_l, avail_u)).sum(); - if cost < best.0 { - best = (cost, m); - } - } - best.1 - }; - // skip = 0 (context: neighbours' skip, always 0 here); tx_size implied (4×4). - self.bc.write(false, DEFAULT_SKIP_PROB[0]); - let above_mode = if avail_u { self.modes[(r - 1) * self.mi_cols + c] } else { DC_PRED }; - let left_mode = if avail_l { self.modes[r * self.mi_cols + c - 1] } else { DC_PRED }; - let o = (above_mode as usize * 10 + left_mode as usize) * 9; - self.bc.tree(&INTRA_MODE_TREE, &KF_Y_MODE_PROBS[o..o + 9], y_mode); - let o = y_mode as usize * 9; - self.bc.tree(&INTRA_MODE_TREE, &KF_UV_MODE_PROBS[o..o + 9], uv_mode); - self.modes[r * self.mi_cols + c] = y_mode; - // Residual: luma 2×2 transform blocks, then one 4×4 block per chroma plane. - for (dy, dx) in [(0, 0), (0, 1), (1, 0), (1, 1)] { - let tx = if self.enc.lossless { DCT_DCT } else { mode_tx_type(y_mode) }; - self.tx_block(0, c * 8 + dx * 4, r * 8 + dy * 4, y_mode, tx, avail_l || dx > 0, avail_u || dy > 0); - } - for p in 1..3 { - self.tx_block(p, c * 4, r * 4, uv_mode, DCT_DCT, avail_l, avail_u); - } - } - - fn cost(&self, plane: usize, x: usize, y: usize, log2: usize, mode: u8, have_left: bool, have_above: bool) -> i64 { - let p = &self.planes[plane]; - let pred = predict(p, x, y, log2, mode, have_left, have_above); - let size = 1 << log2; - let mut s = 0i64; - for i in 0..size { - for j in 0..size { - s += (p.src[(y + i) * p.w + x + j] as i32 - pred[i * size + j]).abs() as i64; - } - } - s - } - - fn choose(&self, plane: usize, x: usize, y: usize, log2: usize, have_left: bool, have_above: bool) -> u8 { - let mut best = (i64::MAX, DC_PRED); - for m in MODES { - let cost = self.cost(plane, x, y, log2, m, have_left, have_above); - if cost < best.0 { - best = (cost, m); - } - } - best.1 - } - - #[allow(clippy::too_many_arguments)] - fn tx_block(&mut self, plane: usize, x: usize, y: usize, mode: u8, tx_type: u8, have_left: bool, have_above: bool) { - let ss = (plane > 0) as usize; - let max_x = (self.mi_cols * 8) >> ss; - let max_y = (self.mi_rows * 8) >> ss; - let (x4, y4) = (x >> 2, y >> 2); - let ly = y4 & (15 >> ss); - if x >= max_x || y >= max_y { - self.above_nz[plane][x4] = 0; - self.left_nz[plane][ly] = 0; - return; - } - let p = &self.planes[plane]; - let pred = predict(p, x, y, 2, mode, have_left, have_above); - let mut res = [0i32; 16]; - for i in 0..4 { - for j in 0..4 { - res[i * 4 + j] = p.src[(y + i) * p.w + x + j] as i32 - pred[i * 4 + j]; - } - } - let (dcq, acq) = (DC_QLOOKUP[self.enc.qidx as usize], AC_QLOOKUP[self.enc.qidx as usize]); - // Quantised levels (raster order). - let levels: [i32; 16] = if self.enc.lossless { - forward_wht(&res) - } else { - let f = &forward_matrices()[tx_type as usize]; - let mut lv = [0i32; 16]; - for k in 0..16 { - let coef: f64 = (0..16).map(|q| f[k][q] * res[q] as f64).sum(); - let qs = if k == 0 { dcq } else { acq } as f64; - let a = coef.abs() / qs; - // Rounding with a small dead zone for AC coefficients. - let l = if k == 0 { (a + 0.5).floor() } else { (a + 0.38).floor() }; - lv[k] = if coef < 0.0 { -(l as i32) } else { l as i32 }; - } - lv - }; - let sc = scan(tx_type); - let eob = (0..16).rev().find(|&i| levels[sc[i] as usize] != 0).map_or(0, |i| i + 1); - // Tokens. - let ctx0 = (self.above_nz[plane][x4] + self.left_nz[plane][ly]) as usize; - self.tokens(plane, tx_type, &levels, eob, ctx0); - // Reconstruct exactly as the decoder does. - let mut deq = [0i32; 16]; - for k in 0..16 { - deq[k] = levels[k] * if k == 0 { dcq } else { acq }; - } - let resid = if eob == 0 { - [0i32; 16] - } else if self.enc.lossless { - inverse_wht(&deq) - } else { - inverse4(&deq, tx_type) - }; - let p = &mut self.planes[plane]; - for i in 0..4 { - for j in 0..4 { - p.rec[(y + i) * p.w + x + j] = (pred[i * 4 + j] + resid[i * 4 + j]).clamp(0, 255) as u8; - } - } - let nz = (eob > 0) as u8; - self.above_nz[plane][x4] = nz; - self.left_nz[plane][ly] = nz; - } - - fn tokens(&mut self, plane: usize, tx_type: u8, levels: &[i32; 16], eob: usize, ctx0: usize) { - let sc = scan(tx_type); - let nb = neighbors(tx_type); - let ptype = (plane > 0) as usize; - // coef_probs[TX_4X4][ptype][intra]: 6 bands × 6 contexts × 3. - let base = (ptype * 2) * 6 * 6 * 3; - let mut tc = [0u8; 16]; - let mut check_eob = true; - let mut c = 0; - while c < 16 { - let pos = sc[c] as usize; - let band = COEFBAND_4X4[c] as usize; - let ctx = if c == 0 { ctx0 } else { ((1 + tc[nb[c].0] + tc[nb[c].1]) >> 1) as usize }; - let o = base + (band * 6 + ctx) * 3; - let p = [DEFAULT_COEF_PROBS[o], DEFAULT_COEF_PROBS[o + 1], DEFAULT_COEF_PROBS[o + 2]]; - if check_eob { - let more = c < eob; - self.bc.write(more, p[0]); - if !more { - break; - } - } - let v = levels[pos].unsigned_abs(); - if v == 0 { - self.bc.write(false, p[1]); - tc[pos] = 0; - check_eob = false; - c += 1; - continue; - } - self.bc.write(true, p[1]); - check_eob = true; - let token = match v { - 1 => 1, - 2 => 2, - 3 => 3, - 4 => 4, - 5..=6 => 5, - 7..=10 => 6, - 11..=18 => 7, - 19..=34 => 8, - 35..=66 => 9, - _ => 10, - }; - if token == 1 { - self.bc.write(false, p[2]); - } else { - self.bc.write(true, p[2]); - let pp = pareto(p[2]); - let path: &[(usize, bool)] = match token { - 2 => &[(0, false), (1, false)], - 3 => &[(0, false), (1, true), (2, false)], - 4 => &[(0, false), (1, true), (2, true)], - 5 => &[(0, true), (3, false), (4, false)], - 6 => &[(0, true), (3, false), (4, true)], - 7 => &[(0, true), (3, true), (5, false), (6, false)], - 8 => &[(0, true), (3, true), (5, false), (6, true)], - 9 => &[(0, true), (3, true), (5, true), (7, false)], - _ => &[(0, true), (3, true), (5, true), (7, true)], - }; - for &(n, b) in path { - self.bc.write(b, pp[n]); - } - let (cat, num_extra, base_v) = EXTRA_BITS[token]; - let extra = (v - base_v as u32).min((1 << num_extra) - 1); - let cp = CAT_PROBS[cat as usize]; - for e in 0..num_extra as usize { - self.bc.write((extra >> (num_extra as usize - 1 - e)) & 1 != 0, cp[e]); - } - } - tc[pos] = ENERGY_CLASS[token]; - self.bc.write(levels[pos] < 0, 128); - c += 1; - } - } -} - // ---------------------------------------------------------------- colour conversion /// BT.709 RGB(A) 8-bit → planar YUV 4:2:0 (limited or full range). @@ -822,36 +453,3 @@ pub fn alpha_to_yuv420(rgba: &[u8], w: u32, h: u32) -> (Vec, Vec, Vec = rgba.chunks_exact(4).take(w * h).map(|p| p[3]).collect(); (y, vec![128; cw * ch], vec![128; cw * ch]) } - -#[cfg(test)] -mod tests { - use super::*; - - #[test] - fn wht_is_exactly_invertible() { - let mut s = 7u32; - for _ in 0..2000 { - let mut r = [0i32; 16]; - for v in &mut r { - s = s.wrapping_mul(1664525).wrapping_add(1013904223); - *v = (s >> 16) as i32 % 511 - 255; - } - let lv = forward_wht(&r); - let deq = lv.map(|l| l * 4); - assert_eq!(inverse_wht(&deq), r); - } - } - - #[test] - fn forward_transforms_invert_the_decoder() { - for ty in 0..4u8 { - let f = &forward_matrices()[ty as usize]; - let r: [i32; 16] = std::array::from_fn(|i| (i as i32 * 37 % 61) - 30); - let c: [i32; 16] = std::array::from_fn(|k| (0..16).map(|q| f[k][q] * r[q] as f64).sum::().round() as i32); - let back = inverse4(&c, ty); - for i in 0..16 { - assert!((back[i] - r[i]).abs() <= 1, "type {ty}: {back:?} vs {r:?}"); - } - } - } -} diff --git a/crates/vp9enc/src/mc.rs b/crates/vp9enc/src/mc.rs new file mode 100644 index 00000000..779a1ea5 --- /dev/null +++ b/crates/vp9enc/src/mc.rs @@ -0,0 +1,207 @@ +//! Inter prediction exactly as the decoder forms it: motion vector clamping (spec 8.5.2.2) and +//! the unscaled block inter prediction process (8.5.2.4) with the regular 8-tap filters and +//! reference samples clamped to the reference frame's visible area. + +use crate::tables::SUBPEL_REGULAR; + +/// A motion vector in 1/8 luma samples. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub struct Mv { + pub row: i32, + pub col: i32, +} + +impl Mv { + pub const ZERO: Mv = Mv { row: 0, col: 0 }; + pub fn new(row: i32, col: i32) -> Mv { + Mv { row, col } + } + pub fn minus(self, o: Mv) -> Mv { + Mv::new(self.row - o.row, self.col - o.col) + } +} + +/// A reference plane: `w`×`h` visible samples at `stride`. +pub struct RefPlane<'a> { + pub data: &'a [u8], + pub stride: usize, + pub w: usize, + pub h: usize, +} + +/// The block's motion vector for `plane` in 1/16 samples of that plane, clamped as the +/// decoder clamps it. `n8`: block size in 8×8 units. +pub fn clamp_for_plane(mv: Mv, r: usize, c: usize, n8: usize, mi_rows: usize, mi_cols: usize, chroma: bool) -> (i32, i32) { + let s = chroma as i32; + let (r, c, b) = (r as i32, c as i32, n8 as i32); + let to_top = (-(r * 8 * 16)) >> s; + let to_bottom = (((mi_rows as i32 - b - r) * 8) * 16) >> s; + let to_left = (-(c * 8 * 16)) >> s; + let to_right = (((mi_cols as i32 - b - c) * 8) * 16) >> s; + let spel = (4 + ((b * 8) >> s)) << 4; + let row = ((2 * mv.row) >> s).clamp(to_top - spel, (to_bottom + spel - 16).max(to_top - spel)); + let col = ((2 * mv.col) >> s).clamp(to_left - spel, (to_right + spel - 16).max(to_left - spel)); + (row, col) +} + +/// Predict a `w`×`h` block whose top-left sample is at (`x16`, `y16`) (1/16 samples) into `out` +/// (stride `w`). +pub fn predict(r: &RefPlane, x16: i32, y16: i32, w: usize, h: usize, out: &mut [u8]) { + let (x0, y0) = (x16 >> 4, y16 >> 4); + let (fx, fy) = ((x16 & 15) as usize, (y16 & 15) as usize); + let (lx, ly) = (r.w as i32 - 1, r.h as i32 - 1); + let inside = x0 - 3 >= 0 && y0 - 3 >= 0 && x0 + w as i32 + 4 <= lx && y0 + h as i32 + 4 <= ly; + let at = |x: i32, y: i32| -> i32 { + if inside { r.data[y as usize * r.stride + x as usize] as i32 } else { r.data[y.clamp(0, ly) as usize * r.stride + x.clamp(0, lx) as usize] as i32 } + }; + let fh = &SUBPEL_REGULAR[fx]; + let fv = &SUBPEL_REGULAR[fy]; + if fx == 0 && fy == 0 { + for i in 0..h { + for j in 0..w { + out[i * w + j] = at(x0 + j as i32, y0 + i as i32) as u8; + } + } + return; + } + // Horizontal pass over h + 7 rows (from y0 - 3), then vertical. + let mut tmp = [0u8; 23 * 16]; + for i in 0..h + 7 { + let y = y0 - 3 + i as i32; + for j in 0..w { + let x = x0 + j as i32 - 3; + let mut acc = 64; + for t in 0..8 { + acc += fh[t] * at(x + t as i32, y); + } + tmp[i * w + j] = (acc >> 7).clamp(0, 255) as u8; + } + } + for i in 0..h { + for j in 0..w { + let mut acc = 64; + for t in 0..8 { + acc += fv[t] * tmp[(i + t) * w + j] as i32; + } + out[i * w + j] = (acc >> 7).clamp(0, 255) as u8; + } + } +} + +/// The reference luma pre-interpolated at the 16 quarter-sample phases the encoder uses (motion +/// vectors are even in 1/8 units: fractions 0, 4, 8, 12 in 1/16), over the visible area plus a +/// border. Each phase plane holds exactly what [`predict`] returns at that phase (the decoder's +/// separable filter with edge clamping), so motion search reads predictions directly. +pub struct PhasePlanes { + border: usize, + stride: usize, + rows: usize, + planes: Vec>, +} + +impl PhasePlanes { + pub fn new(r: &RefPlane, border: usize) -> PhasePlanes { + const M: usize = 4; + let (pw, ph) = (r.w + 2 * border, r.h + 2 * border); + // Edge-replicated source with M more samples on every side for the filter taps. + let rs = pw + 2 * M; + let rh = ph + 2 * M; + let off = (border + M) as i32; + let mut raw = vec![0u8; rs * rh]; + for y in 0..rh { + let sy = (y as i32 - off).clamp(0, r.h as i32 - 1) as usize; + let src = &r.data[sy * r.stride..sy * r.stride + r.w]; + for x in 0..rs { + raw[y * rs + x] = src[(x as i32 - off).clamp(0, r.w as i32 - 1) as usize]; + } + } + // Horizontal phases over every raw row (inner columns). + let hpass = |fx: usize| -> Vec { + let f = &SUBPEL_REGULAR[fx]; + let mut out = vec![0u8; pw * rh]; + for y in 0..rh { + let row = &raw[y * rs..y * rs + rs]; + for x in 0..pw { + let s = &row[x + M - 3..x + M + 5]; + let acc: i32 = 64 + (0..8).map(|t| f[t] * s[t] as i32).sum::(); + out[y * pw + x] = (acc >> 7).clamp(0, 255) as u8; + } + } + out + }; + let unfiltered: Vec = (0..rh).flat_map(|y| raw[y * rs + M..y * rs + M + pw].iter().copied()).collect(); + let h: Vec> = (0..4).map(|k| if k == 0 { unfiltered.clone() } else { hpass(k * 4) }).collect(); + let mut planes = Vec::with_capacity(16); + for fy in 0..4 { + for hp in &h { + if fy == 0 { + planes.push(hp[M * pw..(M + ph) * pw].to_vec()); + continue; + } + let f = &SUBPEL_REGULAR[fy * 4]; + let mut out = vec![0u8; pw * ph]; + for y in 0..ph { + for x in 0..pw { + let mut acc = 64; + for t in 0..8 { + acc += f[t] * hp[(y + M + t - 3) * pw + x] as i32; + } + out[y * pw + x] = (acc >> 7).clamp(0, 255) as u8; + } + } + planes.push(out); + } + } + PhasePlanes { border, stride: pw, rows: ph, planes } + } + + /// The prediction of the `n`×`n` block at (`x16`, `y16`) (fractions multiples of 4) as rows of + /// a phase plane: (plane, offset of the first sample, stride). `None` beyond the border. + pub fn block(&self, x16: i32, y16: i32, n: usize) -> Option<(&[u8], usize, usize)> { + let (x0, y0) = ((x16 >> 4) + self.border as i32, (y16 >> 4) + self.border as i32); + let (fx, fy) = ((x16 & 15) as usize, (y16 & 15) as usize); + if fx % 4 != 0 || fy % 4 != 0 || x0 < 0 || y0 < 0 || x0 as usize + n > self.stride || y0 as usize + n > self.rows { + return None; + } + Some((&self.planes[(fy / 4) * 4 + fx / 4], y0 as usize * self.stride + x0 as usize, self.stride)) + } +} + +#[cfg(test)] +mod tests { + #[test] + fn phase_planes_match_predict() { + let (w, h) = (37usize, 29usize); + let data: Vec = (0..w * h).map(|i| ((i * 7919) % 251) as u8).collect(); + let r = RefPlane { data: &data, stride: w, w, h }; + let pp = PhasePlanes::new(&r, 24); + let mut out = [0u8; 64]; + for y16 in (-20 * 16..(h as i32 + 10) * 16).step_by(4 * 7) { + for x16 in (-20 * 16..(w as i32 + 10) * 16).step_by(4 * 5) { + let Some((p, o, s)) = pp.block(x16, y16, 8) else { continue }; + predict(&r, x16, y16, 8, 8, &mut out); + for i in 0..8 { + assert_eq!(&p[o + i * s..o + i * s + 8], &out[i * 8..i * 8 + 8], "({x16}, {y16}) row {i}"); + } + } + } + } + + use super::*; + + #[test] + fn whole_sample_and_linear_ramps() { + let (w, h) = (40usize, 30usize); + let data: Vec = (0..w * h).map(|i| ((i % w) * 4 + (i / w)) as u8).collect(); + let r = RefPlane { data: &data, stride: w, w, h }; + let mut out = [0u8; 64]; + predict(&r, 10 << 4, 8 << 4, 8, 8, &mut out); + assert_eq!(out[0], data[8 * w + 10]); + // Half a sample right on a linear ramp: exact midpoint. + predict(&r, (10 << 4) + 8, 8 << 4, 8, 8, &mut out); + assert_eq!(out[0], data[8 * w + 10] + 2); + // Off the left edge: clamped to column 0. + predict(&r, -(20 << 4), 0, 4, 4, &mut out); + assert_eq!(out[0], data[0]); + } +} diff --git a/crates/vp9enc/src/probs.rs b/crates/vp9enc/src/probs.rs new file mode 100644 index 00000000..8f55bbba --- /dev/null +++ b/crates/vp9enc/src/probs.rs @@ -0,0 +1,426 @@ +//! Probabilities: the frame's adaptable probability table (one flat array laid out in the order +//! the compressed header updates them), a symbol recorder that codes every boolean decision +//! against either a fixed probability or an entry of that table, bit-cost estimates, and the +//! per-frame forward updates (spec 6.3: `diff_update_prob` with `decode_term_subexp` / +//! `inv_remap_prob`, and `update_mv_prob`). +//! +//! Tiles are first recorded, then the encoder counts how often each table entry coded a 0 or a +//! 1, picks updates that pay for themselves, writes them in the compressed header and finally +//! entropy codes the recorded symbols with the updated probabilities. + +use crate::bool::BoolEncoder; +use crate::tables::*; + +/// Coefficient probabilities `[ptype][ref][band][ctx][3]` for 4×4 transforms (band 0 uses +/// contexts 0..3 only). +pub const COEF: usize = 0; +pub const SKIP: usize = 432; +pub const INTER_MODE: usize = SKIP + 3; +pub const IS_INTER: usize = INTER_MODE + 21; +pub const SINGLE_REF: usize = IS_INTER + 4; +pub const PARTITION: usize = SINGLE_REF + 10; +pub const MV_JOINT: usize = PARTITION + 48; +pub const MV_SIGN: usize = MV_JOINT + 3; +pub const MV_CLASS: usize = MV_SIGN + 2; +pub const MV_CLASS0_BIT: usize = MV_CLASS + 20; +pub const MV_BITS: usize = MV_CLASS0_BIT + 2; +pub const MV_CLASS0_FR: usize = MV_BITS + 20; +pub const MV_FR: usize = MV_CLASS0_FR + 12; +pub const NPROBS: usize = MV_FR + 6; + +/// Index of the first of the three coefficient probabilities of a context. +pub fn coef_index(ptype: usize, inter: bool, band: usize, ctx: usize) -> usize { + COEF + (((ptype * 2 + inter as usize) * 6 + band) * 6 + ctx) * 3 +} + +#[derive(Clone)] +pub struct Probs(pub [u8; NPROBS]); + +impl Probs { + /// The default probabilities (setup_past_independence / key frames). + pub fn defaults() -> Probs { + let mut p = [128u8; NPROBS]; + p[COEF..COEF + 432].copy_from_slice(&DEFAULT_COEF_PROBS); + p[SKIP..SKIP + 3].copy_from_slice(&DEFAULT_SKIP_PROB); + p[INTER_MODE..INTER_MODE + 21].copy_from_slice(&DEFAULT_INTER_MODE_PROBS); + p[IS_INTER..IS_INTER + 4].copy_from_slice(&DEFAULT_IS_INTER_PROB); + p[SINGLE_REF..SINGLE_REF + 10].copy_from_slice(&DEFAULT_SINGLE_REF_PROB); + p[PARTITION..PARTITION + 48].copy_from_slice(&DEFAULT_PARTITION_PROBS); + p[MV_JOINT..MV_JOINT + 3].copy_from_slice(&DEFAULT_MV_JOINT_PROBS); + p[MV_SIGN..MV_SIGN + 2].copy_from_slice(&DEFAULT_MV_SIGN_PROB); + p[MV_CLASS..MV_CLASS + 20].copy_from_slice(&DEFAULT_MV_CLASS_PROBS); + p[MV_CLASS0_BIT..MV_CLASS0_BIT + 2].copy_from_slice(&DEFAULT_MV_CLASS0_BIT_PROB); + p[MV_BITS..MV_BITS + 20].copy_from_slice(&DEFAULT_MV_BITS_PROB); + p[MV_CLASS0_FR..MV_CLASS0_FR + 12].copy_from_slice(&DEFAULT_MV_CLASS0_FR_PROBS); + p[MV_FR..MV_FR + 6].copy_from_slice(&DEFAULT_MV_FR_PROBS); + Probs(p) + } +} + +/// Probabilities of the full token tree (9.3.2) for model probability `p`. +pub fn pareto(p: u8) -> [u8; 8] { + let p = p.max(1) as usize; + let x = (p - 1) / 2; + let mut t = [0u8; 8]; + for n in 0..8 { + t[n] = if p & 1 == 1 { PARETO_TABLE[x * 8 + n] } else { ((PARETO_TABLE[x * 8 + n] as u32 + PARETO_TABLE[(x + 1) * 8 + n] as u32) >> 1) as u8 }; + } + t +} + +/// The probability a recorded symbol is coded with. +#[derive(Clone, Copy, Debug)] +pub enum P { + Fixed(u8), + /// An entry of the frame's [`Probs`]. + Ctx(u16), + /// Node `n` of the Pareto token tree derived from the coefficient probability at the index. + Pareto(u16, u8), +} + +#[derive(Clone, Copy)] +struct Sym { + bit: bool, + p: P, +} + +/// Records the boolean decisions of a tile. +#[derive(Default)] +pub struct Writer { + syms: Vec, +} + +impl Writer { + #[inline] + pub fn put(&mut self, bit: bool, p: P) { + self.syms.push(Sym { bit, p }); + } + #[inline] + pub fn fixed(&mut self, bit: bool, prob: u8) { + self.put(bit, P::Fixed(prob)); + } + #[inline] + pub fn ctx(&mut self, bit: bool, idx: usize) { + self.put(bit, P::Ctx(idx as u16)); + } + /// Write `value` with `tree` (leaves stored negated); node `i` uses `prob(i)`. + pub fn tree(&mut self, tree: &[i8], value: u8, prob: impl Fn(usize) -> P) { + for (node, bit) in tree_path(tree, value) { + self.put(bit, prob(node)); + } + } + + /// Add the 0/1 counts of every table entry this tile used. + pub fn count(&self, counts: &mut [[u32; 2]]) { + for s in &self.syms { + if let P::Ctx(i) = s.p { + counts[i as usize][s.bit as usize] += 1; + } + } + } + + /// Entropy code the recorded symbols. + pub fn emit(&self, probs: &Probs) -> Vec { + let mut bc = BoolEncoder::new(); + bc.write(false, 128); // marker bit + for s in &self.syms { + let p = match s.p { + P::Fixed(p) => p, + P::Ctx(i) => probs.0[i as usize], + P::Pareto(i, n) => pareto(probs.0[i as usize])[n as usize], + }; + bc.write(s.bit, p); + } + bc.finish() + } +} + +/// The (node, bit) path to `value` in a tree (node = tree index / 2). +pub fn tree_path(tree: &[i8], value: u8) -> impl Iterator { + struct Path { + steps: [(u8, bool); 16], + len: usize, + } + fn walk(tree: &[i8], node: usize, value: u8, out: &mut Path) -> bool { + for b in 0..2 { + let t = tree[node + b]; + out.steps[out.len] = ((node >> 1) as u8, b == 1); + out.len += 1; + if t <= 0 { + if (-t) as u8 == value { + return true; + } + } else if walk(tree, t as usize, value, out) { + return true; + } + out.len -= 1; + } + false + } + let mut p = Path { steps: [(0, false); 16], len: 0 }; + walk(tree, 0, value, &mut p); + let Path { steps, len } = p; + steps.into_iter().take(len).map(|(n, b)| (n as usize, b)) +} + +// ---------------------------------------------------------------- costs + +/// Cost in 1/256 bit of coding a 0 (`[p][0]`) or a 1 (`[p][1]`) with probability `p`. +pub fn cost_table() -> &'static [[u32; 2]; 256] { + use std::sync::OnceLock; + static T: OnceLock<[[u32; 2]; 256]> = OnceLock::new(); + T.get_or_init(|| { + let mut t = [[0u32; 2]; 256]; + for (p, e) in t.iter_mut().enumerate() { + let q = (p.max(1) as f64) / 256.0; + e[0] = (-q.log2() * 256.0).round() as u32; + e[1] = (-(1.0 - q).log2() * 256.0).round() as u32; + } + t + }) +} + +#[inline] +pub fn cost(bit: bool, p: u8) -> u32 { + cost_table()[p as usize][bit as usize] +} + +/// Cost of `value` in a tree with probabilities `probs` (node i → probs[i]). +pub fn tree_cost(tree: &[i8], probs: &[u8], value: u8) -> u32 { + tree_path(tree, value).map(|(n, b)| cost(b, probs[n])).sum() +} + +// ---------------------------------------------------------------- forward updates + +fn inv_recenter_nonneg(v: i32, m: i32) -> i32 { + if v > 2 * m { + v + } else if v & 1 != 0 { + m - ((v + 1) >> 1) + } else { + m + (v >> 1) + } +} + +/// `inv_remap_prob` (6.3.5). +fn inv_remap_prob(delta: usize, prob: u8) -> u8 { + let v = inv_map_table()[delta.min(254)] as i32; + let m = prob as i32 - 1; + let r = if (m << 1) <= 255 { 1 + inv_recenter_nonneg(v, m) } else { 255 - inv_recenter_nonneg(v, 255 - 1 - m) }; + r.clamp(1, 255) as u8 +} + +/// The `decode_term_subexp` value that moves `old` to `new`. +fn delta_for(old: u8, new: u8) -> Option { + (0..255).find(|&d| inv_remap_prob(d, old) == new) +} + +/// Bits of `decode_term_subexp` for `d`. +fn subexp_bits(d: usize) -> u32 { + match d { + 0..16 => 5, + 16..32 => 6, + 32..64 => 8, + 64..129 => 10, + _ => 11, + } +} + +fn write_subexp(bc: &mut BoolEncoder, d: usize) { + let d = d as u32; + if d < 16 { + bc.literal(0, 1); + bc.literal(d, 4); + } else if d < 32 { + bc.literal(2, 2); + bc.literal(d - 16, 4); + } else if d < 64 { + bc.literal(6, 3); + bc.literal(d - 32, 5); + } else { + bc.literal(7, 3); + if d < 129 { + bc.literal(d - 64, 7); + } else { + let v = (d + 1) >> 1; + bc.literal(v, 7); + bc.literal((d + 1) & 1, 1); + } + } +} + +fn counts_cost(c: [u32; 2], p: u8) -> u64 { + c[0] as u64 * cost(false, p) as u64 + c[1] as u64 * cost(true, p) as u64 +} + +/// Whether table entry `i` is an MV probability (7-bit `update_mv_prob` rather than +/// `diff_update_prob`). +fn is_mv(i: usize) -> bool { + i >= MV_JOINT +} + +/// Choose updated probabilities: an entry changes when the bits saved on this frame's symbols +/// exceed the cost of signalling the change. `inter`: the inter-frame tables are updatable. +pub fn plan_updates(old: &Probs, counts: &[[u32; 2]], inter: bool) -> Probs { + let mut new = old.clone(); + let flag = cost(true, 252) as i64 - cost(false, 252) as i64; + let limit = if inter { NPROBS } else { INTER_MODE }; + for i in 0..limit { + let c = counts[i]; + let n = c[0] + c[1]; + if n == 0 { + continue; + } + let o = old.0[i]; + let ideal = ((c[0] as f64 * 256.0 / n as f64).round() as i32).clamp(1, 255) as u8; + let cands: Vec = if is_mv(i) { vec![ideal | 1, ideal.saturating_sub(1).max(1) | 1] } else { vec![ideal] }; + let base = counts_cost(c, o) as i64; + let mut best = (0i64, o); + for p in cands { + if p == o { + continue; + } + let signal = if is_mv(i) { + 7 * 256 + } else { + match delta_for(o, p) { + Some(d) => subexp_bits(d) as i64 * 256, + None => continue, + } + }; + let gain = base - counts_cost(c, p) as i64 - signal - flag; + if gain > best.0 { + best = (gain, p); + } + } + new.0[i] = best.1; + } + new +} + +/// `diff_update_prob` for one entry. +fn diff_update(bc: &mut BoolEncoder, old: &Probs, new: &Probs, i: usize) { + let (o, n) = (old.0[i], new.0[i]); + match (o != n).then(|| delta_for(o, n)).flatten() { + Some(d) => { + bc.write(true, 252); + write_subexp(bc, d); + } + None => bc.write(false, 252), + } +} + +/// `update_mv_prob` for one entry. +fn mv_update(bc: &mut BoolEncoder, old: &Probs, new: &Probs, i: usize) { + let n = new.0[i]; + if n != old.0[i] && n & 1 == 1 { + bc.write(true, 252); + bc.literal(n as u32 >> 1, 7); + } else { + bc.write(false, 252); + } +} + +/// The compressed header (6.3) of a frame coded with ONLY_4X4 transforms, single LAST +/// references and a fixed interpolation filter, moving the probabilities from `old` to `new`. +/// Entries `new` cannot signal are left at `old` (and must be coded with `old`): the caller +/// codes the tiles with the probabilities this returns. +pub fn compressed_header(old: &Probs, new: &Probs, lossless: bool, inter: bool) -> (Vec, Probs) { + let mut bc = BoolEncoder::new(); + bc.write(false, 128); + if !lossless { + bc.literal(0, 2); // tx_mode = ONLY_4X4 + } + // read_coef_probs: TX_4X4 only. + let coef_changed = (COEF..COEF + 432).any(|i| new.0[i] != old.0[i]); + bc.literal(coef_changed as u32, 1); + if coef_changed { + for i in 0..2 { + for j in 0..2 { + for k in 0..6 { + for l in 0..if k == 0 { 3 } else { 6 } { + for m in 0..3 { + diff_update(&mut bc, old, new, coef_index(i, j == 1, k, l) + m); + } + } + } + } + } + } + for i in 0..3 { + diff_update(&mut bc, old, new, SKIP + i); + } + if inter { + for i in 0..21 { + diff_update(&mut bc, old, new, INTER_MODE + i); + } + for i in 0..4 { + diff_update(&mut bc, old, new, IS_INTER + i); + } + // reference_mode is SINGLE_REFERENCE (no compound: equal sign biases). + for i in 0..10 { + diff_update(&mut bc, old, new, SINGLE_REF + i); + } + for _ in 0..36 { + bc.write(false, 252); // y mode probabilities (no intra blocks in inter frames) + } + for i in 0..48 { + diff_update(&mut bc, old, new, PARTITION + i); + } + for j in 0..3 { + mv_update(&mut bc, old, new, MV_JOINT + j); + } + for i in 0..2 { + mv_update(&mut bc, old, new, MV_SIGN + i); + for j in 0..10 { + mv_update(&mut bc, old, new, MV_CLASS + i * 10 + j); + } + mv_update(&mut bc, old, new, MV_CLASS0_BIT + i); + for j in 0..10 { + mv_update(&mut bc, old, new, MV_BITS + i * 10 + j); + } + } + for i in 0..2 { + for j in 0..6 { + mv_update(&mut bc, old, new, MV_CLASS0_FR + i * 6 + j); + } + for j in 0..3 { + mv_update(&mut bc, old, new, MV_FR + i * 3 + j); + } + } + } + // What the decoder ends up with. + let mut used = old.clone(); + for i in 0..NPROBS { + let ok = if is_mv(i) { new.0[i] & 1 == 1 } else { delta_for(old.0[i], new.0[i]).is_some() }; + if ok && (inter || i < INTER_MODE) { + used.0[i] = new.0[i]; + } + } + (bc.finish(), used) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn remap_reaches_every_probability() { + for old in 1..=255u8 { + for new in 1..=255u8 { + if new != old { + assert!(delta_for(old, new).is_some(), "{old} -> {new}"); + } + } + } + } + + #[test] + fn inv_map_table_shape() { + let t = inv_map_table(); + assert_eq!(&t[..4], &[7, 20, 33, 46]); + assert_eq!(t[19], 254); + assert_eq!(&t[20..27], &[1, 2, 3, 4, 5, 6, 8]); + assert_eq!(&t[252..], &[252, 253, 253]); + } +} diff --git a/crates/vp9enc/src/tables.rs b/crates/vp9enc/src/tables.rs index dc8027a6..ffb26641 100644 --- a/crates/vp9enc/src/tables.rs +++ b/crates/vp9enc/src/tables.rs @@ -128,3 +128,76 @@ pub const DEFAULT_COEF_PROBS: [u8; 432] = [ 79, 1, 20, 39, 99, 167, 239, 47, 141, 224, 10, 104, 178, 2, 73, 133, 1, 44, 85, 1, 22, 47, 127, 145, 243, 71, 129, 228, 17, 93, 177, 3, 61, 124, 1, 41, 84, 1, 21, 52, 157, 78, 244, 140, 72, 231, 69, 58, 184, 31, 44, 137, 14, 38, 105, 8, 23, 61, ]; + +// ---------------------------------------------------------------- inter frames (spec section 10.5) + +/// `inv_map_table` (used by `inv_remap_prob`): the 20 values 7 + 13k first, then every other value +/// 1..=253 in increasing order, then 253 again. +pub fn inv_map_table() -> &'static [u8; 255] { + use std::sync::OnceLock; + static T: OnceLock<[u8; 255]> = OnceLock::new(); + T.get_or_init(|| { + let mut t = [253u8; 255]; + let mut n = 0; + for k in 0..20 { + t[n] = 7 + 13 * k; + n += 1; + } + for v in 1..=253u8 { + if v < 7 || (v - 7) % 13 != 0 { + t[n] = v; + n += 1; + } + } + t + }) +} + +/// The regular 8-tap sub-sample interpolation filters (`subpel_filters[EIGHTTAP]`), 16 positions. +pub const SUBPEL_REGULAR: [[i32; 8]; 16] = [ + [0, 0, 0, 128, 0, 0, 0, 0], + [0, 1, -5, 126, 8, -3, 1, 0], + [-1, 3, -10, 122, 18, -6, 2, 0], + [-1, 4, -13, 118, 27, -9, 3, -1], + [-1, 4, -16, 112, 37, -11, 4, -1], + [-1, 5, -18, 105, 48, -14, 4, -1], + [-1, 5, -19, 97, 58, -16, 5, -1], + [-1, 6, -19, 88, 68, -18, 5, -1], + [-1, 6, -19, 78, 78, -19, 6, -1], + [-1, 5, -18, 68, 88, -19, 6, -1], + [-1, 5, -16, 58, 97, -19, 5, -1], + [-1, 4, -14, 48, 105, -18, 5, -1], + [-1, 4, -11, 37, 112, -16, 4, -1], + [-1, 3, -9, 27, 118, -13, 4, -1], + [0, 2, -6, 18, 122, -10, 3, -1], + [0, 1, -3, 8, 126, -5, 1, 0], +]; + +/// `mv_ref_blocks` (row, column offsets) for BLOCK_8X8 and BLOCK_16X16. +pub const MV_REF_8X8: [(i32, i32); 8] = [(-1, 0), (0, -1), (-1, -1), (-2, 0), (0, -2), (-2, -1), (-1, -2), (-2, -2)]; +pub const MV_REF_16X16: [(i32, i32); 8] = [(-1, 0), (0, -1), (-1, 1), (1, -1), (-1, -1), (-3, 0), (0, -3), (-3, -3)]; + +/// `mode_2_counter` indexed by y mode (intra modes 0..=9, then NEARESTMV, NEARMV, ZEROMV, NEWMV). +pub const MODE_2_COUNTER: [u8; 14] = [9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 0, 0, 3, 1]; + +/// `counter_to_context`; INVALID_CASE entries are 9. +pub const COUNTER_TO_CONTEXT: [u8; 19] = [2, 3, 4, 1, 3, 9, 0, 9, 9, 5, 5, 9, 5, 9, 9, 9, 9, 9, 6]; + +pub const DEFAULT_PARTITION_PROBS: [u8; 48] = [ + 199, 122, 141, 147, 63, 159, 148, 133, 118, 121, 104, 114, 174, 73, 87, 92, 41, 83, 82, 99, 50, 53, 39, 39, 177, 58, 59, 68, 26, 63, 52, 79, 25, 17, 14, + 12, 222, 34, 30, 72, 16, 44, 58, 32, 12, 10, 7, 6, +]; + +pub const DEFAULT_IS_INTER_PROB: [u8; 4] = [9, 102, 187, 225]; + +pub const DEFAULT_SINGLE_REF_PROB: [u8; 10] = [33, 16, 77, 74, 142, 142, 172, 170, 238, 247]; + +pub const DEFAULT_INTER_MODE_PROBS: [u8; 21] = [2, 173, 34, 7, 145, 85, 7, 166, 63, 7, 94, 66, 8, 64, 46, 17, 81, 31, 25, 29, 30]; + +pub const DEFAULT_MV_JOINT_PROBS: [u8; 3] = [32, 64, 96]; +pub const DEFAULT_MV_SIGN_PROB: [u8; 2] = [128, 128]; +pub const DEFAULT_MV_CLASS_PROBS: [u8; 20] = [224, 144, 192, 168, 192, 176, 192, 198, 198, 245, 216, 128, 176, 160, 176, 176, 192, 198, 198, 208]; +pub const DEFAULT_MV_CLASS0_BIT_PROB: [u8; 2] = [216, 208]; +pub const DEFAULT_MV_BITS_PROB: [u8; 20] = [136, 140, 148, 160, 176, 192, 224, 234, 234, 240, 136, 140, 148, 160, 176, 192, 224, 234, 234, 240]; +pub const DEFAULT_MV_CLASS0_FR_PROBS: [u8; 12] = [128, 128, 64, 96, 112, 64, 128, 128, 64, 96, 112, 64]; +pub const DEFAULT_MV_FR_PROBS: [u8; 6] = [64, 96, 64, 64, 96, 64]; diff --git a/crates/vp9enc/src/tile.rs b/crates/vp9enc/src/tile.rs new file mode 100644 index 00000000..fe3a1a49 --- /dev/null +++ b/crates/vp9enc/src/tile.rs @@ -0,0 +1,894 @@ +//! Tile coding: partitions, intra blocks (key frames) and inter blocks (motion search, mode +//! decision, motion vector prediction, residual), recorded as symbols (see [`crate::probs`]). +//! +//! Mode info, motion vector candidate lists and contexts follow the decoding process of the +//! VP9 specification (6.4, 6.4.x `find_mv_refs` / `find_best_ref_mvs`) for the subset of the +//! syntax this encoder emits: 8×8 and 16×16 blocks, single LAST references, no +//! `UsePrevFrameMvs` (inter frames are error resilient), no high-precision vectors. + +use crate::mc::{Mv, PhasePlanes, RefPlane, clamp_for_plane, predict}; +use crate::probs::{self, P, Writer, cost, tree_cost}; +use crate::tables::*; +use crate::transform::{ADST_ADST, ADST_DCT, DCT_ADST, DCT_DCT, Quant}; +use crate::{BLOCK_8X8, BLOCK_16X16, Mi, Plane}; + +pub const DC_PRED: u8 = 0; +pub const V_PRED: u8 = 1; +pub const H_PRED: u8 = 2; +pub const TM_PRED: u8 = 9; +const MODES: [u8; 4] = [DC_PRED, V_PRED, H_PRED, TM_PRED]; +pub const NEARESTMV: u8 = 10; +pub const NEARMV: u8 = 11; +pub const ZEROMV: u8 = 12; +pub const NEWMV: u8 = 13; + +const PARTITION_TREE: [i8; 6] = [0, 2, -1, 4, -2, -3]; +const INTRA_MODE_TREE: [i8; 18] = [0, 2, -9, 4, -1, 6, 8, 12, -2, 10, -4, -5, -3, 14, -8, 16, -6, -7]; +/// Leaves are `mode - NEARESTMV`. +const INTER_MODE_TREE: [i8; 6] = [-2, 2, 0, 4, -1, -3]; +const MV_JOINT_TREE: [i8; 6] = [0, 2, -1, 4, -2, -3]; +const MV_CLASS_TREE: [i8; 20] = [0, 2, -1, 4, 6, 8, -2, -3, 10, 12, -4, -5, -6, 14, 16, 18, -7, -8, -9, -10]; +const MV_FR_TREE: [i8; 6] = [0, 2, -1, 4, -2, -3]; + +/// extra_bits[token] = (cat, numExtra, base). +const EXTRA_BITS: [(u8, u8, u16); 11] = + [(0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 0, 3), (0, 0, 4), (1, 1, 5), (2, 2, 7), (3, 3, 11), (4, 4, 19), (5, 5, 35), (6, 14, 67)]; +const CAT_PROBS: [&[u8]; 7] = [ + &[], + &[159], + &[165, 145], + &[173, 148, 140], + &[176, 155, 140, 135], + &[180, 157, 141, 134, 130], + &[254, 254, 254, 252, 249, 243, 230, 196, 177, 153, 140, 133, 130, 129], +]; + +/// mode2txfm_map for the intra modes used here (inter blocks use DCT_DCT). +fn mode_tx_type(mode: u8) -> u8 { + match mode { + V_PRED => ADST_DCT, + H_PRED => DCT_ADST, + TM_PRED => ADST_ADST, + _ => DCT_DCT, + } +} + +fn scan(tx_type: u8) -> &'static [u16; 16] { + match tx_type { + ADST_DCT => &ROW_SCAN_4X4, + DCT_ADST => &COL_SCAN_4X4, + _ => &DEFAULT_SCAN_4X4, + } +} + +/// Coefficient context neighbours (9.3.2) of each scan position. +fn neighbors(tx_type: u8) -> [(usize, usize); 16] { + let sc = scan(tx_type); + let mut nb = [(0, 0); 16]; + for (c, &pos) in sc.iter().enumerate().skip(1) { + let pos = pos as usize; + let (i, j) = (pos / 4, pos % 4); + nb[c] = if i > 0 && j > 0 { + let a = (i - 1) * 4 + j; + let b = i * 4 + j - 1; + match tx_type { + DCT_ADST => (a, a), + ADST_DCT => (b, b), + _ => (a, b), + } + } else if i > 0 { + ((i - 1) * 4 + j, (i - 1) * 4 + j) + } else { + (j - 1, j - 1) + }; + } + nb +} + +fn eob_of(levels: &[i32; 16], tx_type: u8) -> usize { + let sc = scan(tx_type); + (0..16).rev().find(|&i| levels[sc[i] as usize] != 0).map_or(0, |i| i + 1) +} + +/// Rough coefficient bit cost (1/256 bit) for rate-distortion decisions. +fn levels_cost(levels: &[i32; 16], eob: usize) -> u32 { + if eob == 0 { + return 0; + } + let mut bits = 256; + for &l in levels { + let a = l.unsigned_abs(); + bits += if a == 0 { 160 } else { 640 + 384 * (32 - a.leading_zeros()) }; + } + bits +} + +/// Intra prediction (8.5.1) of a square block for DC / V / H / TM. +fn predict_intra(p: &Plane, x: usize, y: usize, log2: usize, mode: u8, have_left: bool, have_above: bool) -> [i32; 64] { + let size = 1usize << log2; + let base = 128i32; + let buf = &p.rec; + let w = p.w; + let mut above = [0i32; 8]; + let mut left = [0i32; 8]; + if have_above { + for i in 0..size { + above[i] = buf[(y - 1) * w + (x + i).min(p.max_x)] as i32; + } + } else { + above[..size].fill(base - 1); + } + let top_left = if have_above && have_left { + buf[(y - 1) * w + (x - 1).min(p.max_x)] as i32 + } else if have_above { + base + 1 + } else { + base - 1 + }; + if have_left { + for i in 0..size { + left[i] = buf[(y + i).min(p.max_y) * w + x - 1] as i32; + } + } else { + left[..size].fill(base + 1); + } + let mut out = [0i32; 64]; + for i in 0..size { + for j in 0..size { + out[i * size + j] = match mode { + V_PRED => above[j], + H_PRED => left[i], + TM_PRED => (above[j] + left[i] - top_left).clamp(0, 255), + _ => 0, + }; + } + } + if mode == DC_PRED { + let v = match (have_left, have_above) { + (true, true) => (left[..size].iter().sum::() + above[..size].iter().sum::() + size as i32) >> (log2 + 1), + (true, false) => (left[..size].iter().sum::() + (1 << (log2 - 1))) >> log2, + (false, true) => (above[..size].iter().sum::() + (1 << (log2 - 1))) >> log2, + (false, false) => base, + }; + out[..size * size].fill(v); + } + out +} + +// ---------------------------------------------------------------- motion vector coding + +/// `mv_class` of an offset (|component| - 1). +fn mv_class(offset: u32) -> u32 { + if offset < 16 { 0 } else { (31 - offset.leading_zeros()) - 3 } +} + +fn mv_component_cost(i: usize, v: i32) -> u32 { + let d = &probs_defaults().0; + let offset = v.unsigned_abs() - 1; + let c = mv_class(offset); + let mut bits = cost(v < 0, d[probs::MV_SIGN + i]); + bits += tree_cost(&MV_CLASS_TREE, &d[probs::MV_CLASS + i * 10..], c as u8); + let fr = ((offset >> 1) & 3) as u8; + if c == 0 { + let b = offset >> 3; + bits += cost(b == 1, d[probs::MV_CLASS0_BIT + i]); + bits += tree_cost(&MV_FR_TREE, &d[probs::MV_CLASS0_FR + i * 6 + b as usize * 3..], fr); + } else { + let dd = (offset - (1 << (c + 3))) >> 3; + for b in 0..c { + bits += cost((dd >> b) & 1 == 1, d[probs::MV_BITS + i * 10 + b as usize]); + } + bits += tree_cost(&MV_FR_TREE, &d[probs::MV_FR + i * 3..], fr); + } + bits +} + +const MV_TABLE: i32 = 8192; + +/// Component costs for |v| ≤ MV_TABLE (index v + MV_TABLE) and the four joint costs. +fn mv_tables() -> &'static ([Vec; 2], [u32; 4]) { + use std::sync::OnceLock; + static T: OnceLock<([Vec; 2], [u32; 4])> = OnceLock::new(); + T.get_or_init(|| { + let comp = std::array::from_fn(|i| (-MV_TABLE..=MV_TABLE).map(|v| if v == 0 { 0 } else { mv_component_cost(i, v) }).collect()); + let d = &probs_defaults().0; + let joints = std::array::from_fn(|j| tree_cost(&MV_JOINT_TREE, &d[probs::MV_JOINT..], j as u8)); + (comp, joints) + }) +} + +/// Bits (1/256) of a NEWMV difference with the default probabilities. +pub fn mv_cost(diff: Mv) -> u32 { + let (comp, joints) = mv_tables(); + let joint = (diff.row != 0) as usize * 2 + (diff.col != 0) as usize; + let c = |i: usize, v: i32| if v.abs() <= MV_TABLE { comp[i][(v + MV_TABLE) as usize] } else { mv_component_cost(i, v) }; + joints[joint] + if diff.row != 0 { c(0, diff.row) } else { 0 } + if diff.col != 0 { c(1, diff.col) } else { 0 } +} + +fn probs_defaults() -> &'static probs::Probs { + use std::sync::OnceLock; + static D: OnceLock = OnceLock::new(); + D.get_or_init(probs::Probs::defaults) +} + +// ---------------------------------------------------------------- the tile + +/// Per-frame coding parameters shared by the tiles. +pub struct FrameParams<'a> { + pub key: bool, + pub quant: Quant, + pub mi_cols: usize, + pub mi_rows: usize, + /// Reference planes (inter frames). + pub refs: Option<[RefPlane<'a>; 3]>, + /// The reference luma at every quarter-sample phase (motion search). + pub phases: Option, + /// Encoder-side motion hints from the previous frame (by MI position). + pub prev_mvs: Option<&'a [Mv]>, + /// Lagrangian multipliers: per bit, in SAD and in squared-error units. + pub lambda_sad: f64, + pub lambda_sse: f64, +} + +pub struct Tile<'a> { + pub f: &'a FrameParams<'a>, + pub w: Writer, + pub planes: &'a mut [Plane; 3], + pub mi: &'a mut [Mi], + pub col_start: usize, + pub col_end: usize, + above_nz: [Vec; 3], + left_nz: [Vec; 3], + above_part: Vec, + left_part: Vec, + /// Motion vectors of the four 8×8 blocks evaluated for the current 16×16 area. + pending: [(usize, usize, Mv); 4], +} + +impl<'a> Tile<'a> { + pub fn new(f: &'a FrameParams<'a>, planes: &'a mut [Plane; 3], mi: &'a mut [Mi], col_start: usize, col_end: usize) -> Tile<'a> { + let n = f.mi_cols * 2 + 16; + Tile { + f, + w: Writer::default(), + planes, + mi, + col_start, + col_end, + above_nz: [vec![0; n], vec![0; n], vec![0; n]], + left_nz: [vec![0; 16], vec![0; 16], vec![0; 16]], + above_part: vec![0; f.mi_cols + 8], + left_part: vec![0; 8], + pending: [(usize::MAX, 0, Mv::ZERO); 4], + } + } + + /// Code the tile's superblocks. + pub fn code(&mut self) { + for r in (0..self.f.mi_rows).step_by(8) { + for p in 0..3 { + self.left_nz[p].fill(0); + } + self.left_part.fill(0); + for c in (self.col_start..self.col_end).step_by(8) { + self.partition(r, c, 12); + } + } + } + + fn mi_at(&self, r: usize, c: usize) -> Mi { + self.mi[r * self.f.mi_cols + c] + } + + fn set_mi(&mut self, r: usize, c: usize, n8: usize, m: Mi) { + for dr in 0..n8 { + for dc in 0..n8 { + if r + dr < self.f.mi_rows && c + dc < self.f.mi_cols { + self.mi[(r + dr) * self.f.mi_cols + c + dc] = m; + } + } + } + } + + fn partition(&mut self, r: usize, c: usize, bsize: u8) { + let (mi_rows, mi_cols) = (self.f.mi_rows, self.f.mi_cols); + if r >= mi_rows || c >= mi_cols { + return; + } + let num8 = [1usize, 1, 1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8][bsize as usize]; + let half = num8 >> 1; + let has_rows = (r + half) < mi_rows; + let has_cols = (c + half) < mi_cols; + let bsl = [0usize, 0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3][bsize as usize]; + let boffset = 3 - bsl; + let mut above = 0u8; + let mut left = 0u8; + for i in 0..num8 { + above |= self.above_part[c + i]; + left |= self.left_part[(r & 7) + i]; + } + let ctx = bsl * 4 + (((left >> boffset) & 1) as usize) * 2 + ((above >> boffset) & 1) as usize; + let key = self.f.key; + let prob = |i: usize| if key { P::Fixed(KF_PARTITION_PROBS[ctx * 3 + i]) } else { P::Ctx((probs::PARTITION + ctx * 3 + i) as u16) }; + if bsize == BLOCK_8X8 { + self.w.put(false, prob(0)); // PARTITION_NONE + if key { + self.intra_block(r, c); + } else { + let mv = self.pending.iter().find(|p| p.0 == r && p.1 == c).map(|p| p.2); + let mv = mv.unwrap_or_else(|| self.search(r, c, BLOCK_8X8).0); + self.inter_block(r, c, BLOCK_8X8, mv); + } + self.above_part[c] = 15 >> 1; + self.left_part[r & 7] = 15 >> 1; + return; + } + if !key && bsize == BLOCK_16X16 && has_rows && has_cols { + if let Some(mv) = self.decide16(r, c) { + self.w.put(false, prob(0)); // PARTITION_NONE + self.inter_block(r, c, BLOCK_16X16, mv); + for i in 0..2 { + self.above_part[c + i] = 15 >> 2; + self.left_part[(r & 7) + i] = 15 >> 2; + } + return; + } + } else if !key && bsize == BLOCK_16X16 { + self.pending = [(usize::MAX, 0, Mv::ZERO); 4]; + } + // PARTITION_SPLIT. + if has_rows && has_cols { + self.w.tree(&PARTITION_TREE, 3, prob); + } else if has_cols { + self.w.put(true, prob(1)); + } else if has_rows { + self.w.put(true, prob(2)); + } + let sub = bsize - 3; + self.partition(r, c, sub); + self.partition(r, c + half, sub); + self.partition(r + half, c, sub); + self.partition(r + half, c + half, sub); + } + + // ------------------------------------------------------------ intra (key frames) + + fn intra_block(&mut self, r: usize, c: usize) { + let avail_u = r > 0; + let avail_l = c > self.col_start; + let y_mode = self.choose_intra(0, c * 8, r * 8, 3, avail_l, avail_u); + let uv_mode = { + let mut best = (i64::MAX, DC_PRED); + for m in MODES { + let cost: i64 = (1..3).map(|p| self.intra_cost(p, c * 4, r * 4, 2, m, avail_l, avail_u)).sum(); + if cost < best.0 { + best = (cost, m); + } + } + best.1 + }; + // skip = 0 (neighbours never skip in key frames: context 0); tx_size implied (4×4). + self.w.ctx(false, probs::SKIP); + let above_mode = if avail_u { self.mi_at(r - 1, c).mode } else { DC_PRED }; + let left_mode = if avail_l { self.mi_at(r, c - 1).mode } else { DC_PRED }; + let o = (above_mode as usize * 10 + left_mode as usize) * 9; + self.w.tree(&INTRA_MODE_TREE, y_mode, |i| P::Fixed(KF_Y_MODE_PROBS[o + i])); + let o = y_mode as usize * 9; + self.w.tree(&INTRA_MODE_TREE, uv_mode, |i| P::Fixed(KF_UV_MODE_PROBS[o + i])); + self.set_mi(r, c, 1, Mi { inter: false, mode: y_mode, mv: Mv::ZERO, skip: false, size: BLOCK_8X8 }); + // Residual: luma 2×2 transform blocks, then one 4×4 block per chroma plane. + for (dy, dx) in [(0, 0), (0, 1), (1, 0), (1, 1)] { + let tx = if self.f.quant.lossless { DCT_DCT } else { mode_tx_type(y_mode) }; + self.intra_tx(0, c * 8 + dx * 4, r * 8 + dy * 4, y_mode, tx, avail_l || dx > 0, avail_u || dy > 0); + } + for p in 1..3 { + self.intra_tx(p, c * 4, r * 4, uv_mode, DCT_DCT, avail_l, avail_u); + } + } + + #[allow(clippy::too_many_arguments)] + fn intra_cost(&self, plane: usize, x: usize, y: usize, log2: usize, mode: u8, have_left: bool, have_above: bool) -> i64 { + let p = &self.planes[plane]; + let pred = predict_intra(p, x, y, log2, mode, have_left, have_above); + let size = 1 << log2; + let mut s = 0i64; + for i in 0..size { + for j in 0..size { + s += (p.src[(y + i) * p.w + x + j] as i32 - pred[i * size + j]).abs() as i64; + } + } + s + } + + fn choose_intra(&self, plane: usize, x: usize, y: usize, log2: usize, have_left: bool, have_above: bool) -> u8 { + let mut best = (i64::MAX, DC_PRED); + for m in MODES { + let cost = self.intra_cost(plane, x, y, log2, m, have_left, have_above); + if cost < best.0 { + best = (cost, m); + } + } + best.1 + } + + #[allow(clippy::too_many_arguments)] + fn intra_tx(&mut self, plane: usize, x: usize, y: usize, mode: u8, tx_type: u8, have_left: bool, have_above: bool) { + let ss = (plane > 0) as usize; + let max_x = (self.f.mi_cols * 8) >> ss; + let max_y = (self.f.mi_rows * 8) >> ss; + let (x4, y4) = (x >> 2, y >> 2); + let ly = y4 & (15 >> ss); + if x >= max_x || y >= max_y { + self.above_nz[plane][x4] = 0; + self.left_nz[plane][ly] = 0; + return; + } + let p = &self.planes[plane]; + let pred = predict_intra(p, x, y, 2, mode, have_left, have_above); + let mut res = [0i32; 16]; + for i in 0..4 { + for j in 0..4 { + res[i * 4 + j] = p.src[(y + i) * p.w + x + j] as i32 - pred[i * 4 + j]; + } + } + let q = Quant { round_dc: 0.5, round_ac: 0.38, ..self.f.quant }; + let levels = q.levels(&res, tx_type); + let eob = eob_of(&levels, tx_type); + let ctx0 = (self.above_nz[plane][x4] + self.left_nz[plane][ly]) as usize; + self.tokens(plane, false, tx_type, &levels, eob, ctx0); + let resid = q.residual(&levels, tx_type, eob); + let p = &mut self.planes[plane]; + for i in 0..4 { + for j in 0..4 { + p.rec[(y + i) * p.w + x + j] = (pred[i * 4 + j] + resid[i * 4 + j]).clamp(0, 255) as u8; + } + } + let nz = (eob > 0) as u8; + self.above_nz[plane][x4] = nz; + self.left_nz[plane][ly] = nz; + } + + fn tokens(&mut self, plane: usize, inter: bool, tx_type: u8, levels: &[i32; 16], eob: usize, ctx0: usize) { + let sc = scan(tx_type); + let nb = neighbors(tx_type); + let ptype = (plane > 0) as usize; + let mut tc = [0u8; 16]; + let mut check_eob = true; + let mut c = 0; + while c < 16 { + let pos = sc[c] as usize; + let band = COEFBAND_4X4[c] as usize; + let ctx = if c == 0 { ctx0 } else { ((1 + tc[nb[c].0] + tc[nb[c].1]) >> 1) as usize }; + let o = probs::coef_index(ptype, inter, band, ctx); + if check_eob { + let more = c < eob; + self.w.ctx(more, o); + if !more { + break; + } + } + let v = levels[pos].unsigned_abs(); + if v == 0 { + self.w.ctx(false, o + 1); + tc[pos] = 0; + check_eob = false; + c += 1; + continue; + } + self.w.ctx(true, o + 1); + check_eob = true; + let token = match v { + 1 => 1, + 2 => 2, + 3 => 3, + 4 => 4, + 5..=6 => 5, + 7..=10 => 6, + 11..=18 => 7, + 19..=34 => 8, + 35..=66 => 9, + _ => 10, + }; + if token == 1 { + self.w.ctx(false, o + 2); + } else { + self.w.ctx(true, o + 2); + let path: &[(u8, bool)] = match token { + 2 => &[(0, false), (1, false)], + 3 => &[(0, false), (1, true), (2, false)], + 4 => &[(0, false), (1, true), (2, true)], + 5 => &[(0, true), (3, false), (4, false)], + 6 => &[(0, true), (3, false), (4, true)], + 7 => &[(0, true), (3, true), (5, false), (6, false)], + 8 => &[(0, true), (3, true), (5, false), (6, true)], + 9 => &[(0, true), (3, true), (5, true), (7, false)], + _ => &[(0, true), (3, true), (5, true), (7, true)], + }; + for &(n, b) in path { + self.w.put(b, P::Pareto((o + 2) as u16, n)); + } + let (cat, num_extra, base_v) = EXTRA_BITS[token]; + let extra = (v - base_v as u32).min((1 << num_extra) - 1); + let cp = CAT_PROBS[cat as usize]; + for e in 0..num_extra as usize { + self.w.fixed((extra >> (num_extra as usize - 1 - e)) & 1 != 0, cp[e]); + } + } + tc[pos] = ENERGY_CLASS[token]; + self.w.fixed(levels[pos] < 0, 128); + c += 1; + } + } + + // ------------------------------------------------------------ inter + + fn is_inside(&self, r: i32, c: i32) -> bool { + r >= 0 && (r as usize) < self.f.mi_rows && c >= self.col_start as i32 && (c as usize) < self.col_end + } + + /// `find_mv_refs` + `find_best_ref_mvs` for LAST_FRAME: (NearestMv, NearMv, mode context). + fn mv_refs(&self, r: usize, c: usize, bsize: u8) -> (Mv, Mv, usize) { + let search = if bsize == BLOCK_16X16 { &MV_REF_16X16 } else { &MV_REF_8X8 }; + let mut list = [Mv::ZERO; 2]; + let mut count = 0; + let add = |mv: Mv, list: &mut [Mv; 2], count: &mut usize| { + if *count >= 2 || (*count > 0 && mv == list[0]) { + return; + } + list[*count] = mv; + *count += 1; + }; + let mut counter = 0usize; + for (i, &(dr, dc)) in search.iter().enumerate() { + let (cr, cc) = (r as i32 + dr, c as i32 + dc); + if !self.is_inside(cr, cc) { + continue; + } + let m = self.mi_at(cr as usize, cc as usize); + if i < 2 { + counter += MODE_2_COUNTER[m.mode as usize] as usize; + } + if m.inter { + add(m.mv, &mut list, &mut count); + } + } + let ctx = COUNTER_TO_CONTEXT[counter.min(18)] as usize; + let n8 = if bsize == BLOCK_16X16 { 2 } else { 1 }; + let clamp = |mv: Mv, border: i32| -> Mv { + let (r, c, b) = (r as i32, c as i32, n8); + let to_top = -(r * 64); + let to_bottom = (self.f.mi_rows as i32 - b - r) * 64; + let to_left = -(c * 64); + let to_right = (self.f.mi_cols as i32 - b - c) * 64; + Mv::new( + mv.row.clamp(to_top - border, (to_bottom + border).max(to_top - border)), + mv.col.clamp(to_left - border, (to_right + border).max(to_left - border)), + ) + }; + for mv in list.iter_mut() { + *mv = clamp(*mv, 128); + // No high precision: odd components move toward zero. + let even = |v: i32| if v & 1 != 0 { v + if v > 0 { -1 } else { 1 } } else { v }; + *mv = clamp(Mv::new(even(mv.row), even(mv.col)), (160 - 4) << 3); + } + (list[0], list[1], ctx) + } + + /// The allowed motion vector range of a block (no decoder clamping inside it). + fn mv_range(&self, r: usize, c: usize, n: i32) -> (i32, i32, i32, i32) { + let n8 = n / 8; + let (r, c) = (r as i32, c as i32); + let lo_c = -(c * 64) - (4 + n) * 8; + let hi_c = (self.f.mi_cols as i32 - n8 - c) * 64 + (4 + n) * 8 - 8; + let lo_r = -(r * 64) - (4 + n) * 8; + let hi_r = (self.f.mi_rows as i32 - n8 - r) * 64 + (4 + n) * 8 - 8; + let lim = 8 * 1000; + (lo_r.max(-lim), hi_r.min(lim), lo_c.max(-lim), hi_c.min(lim)) + } + + /// SAD of the luma prediction at `mv` against the source. + fn sad(&self, r: usize, c: usize, n: usize, mv: Mv) -> u32 { + let p = &self.planes[0]; + let (x16, y16) = (((c * 8) << 4) as i32 + 2 * mv.col, ((r * 8) << 4) as i32 + 2 * mv.row); + let mut s = 0u32; + if let Some((pl, o, st)) = self.f.phases.as_ref().and_then(|ph| ph.block(x16, y16, n)) { + for i in 0..n { + let a = &p.src[(r * 8 + i) * p.w + c * 8..][..n]; + let b = &pl[o + i * st..o + i * st + n]; + s += a.iter().zip(b).map(|(a, b)| (*a as i32 - *b as i32).unsigned_abs()).sum::(); + } + return s; + } + let refs = self.f.refs.as_ref().expect("reference"); + let mut pred = [0u8; 256]; + predict(&refs[0], x16, y16, n, n, &mut pred); + for i in 0..n { + let a = &p.src[(r * 8 + i) * p.w + c * 8..][..n]; + s += a.iter().zip(&pred[i * n..i * n + n]).map(|(a, b)| (*a as i32 - *b as i32).unsigned_abs()).sum::(); + } + s + } + + /// Motion search for the block at (`r`, `c`): the best vector and its cost + /// (SAD + λ·bits). + fn search(&self, r: usize, c: usize, bsize: u8) -> (Mv, f64) { + let n = if bsize == BLOCK_16X16 { 16 } else { 8 }; + let (nearest, near, mctx) = self.mv_refs(r, c, bsize); + let (lo_r, hi_r, lo_c, hi_c) = self.mv_range(r, c, n as i32); + let inside = |mv: Mv| mv.row >= lo_r && mv.row <= hi_r && mv.col >= lo_c && mv.col <= hi_c; + let lambda = self.f.lambda_sad / 256.0; + let d = &probs_defaults().0; + let mb: [u32; 4] = std::array::from_fn(|m| tree_cost(&INTER_MODE_TREE, &d[probs::INTER_MODE + mctx * 3..], m as u8)); + let mode_bits = |m: u8| mb[(m - NEARESTMV) as usize]; + let cost = |mv: Mv| -> f64 { + let bits = if mv == nearest { + mode_bits(NEARESTMV) + } else if mv == Mv::ZERO { + mode_bits(ZEROMV) + } else if mv == near { + mode_bits(NEARMV) + } else { + mode_bits(NEWMV) + mv_cost(mv.minus(nearest)) + }; + self.sad(r, c, n, mv) as f64 + lambda * bits as f64 + }; + let mut best = (Mv::ZERO, cost(Mv::ZERO)); + let prev = self.f.prev_mvs.map(|p| p[r * self.f.mi_cols + c]); + for cand in [Some(nearest), Some(near), prev].into_iter().flatten() { + if cand != Mv::ZERO && inside(cand) { + let j = cost(cand); + if j < best.1 { + best = (cand, j); + } + } + } + // A perfect match needs no search. + let good_enough = lambda * 256.0 * 0.5; + if best.1 > good_enough { + // Whole-sample search from the best candidate: a coarse cross, then a small diamond. + let round8 = |v: i32| (v + 4).div_euclid(8) * 8; + let mut cur = Mv::new(round8(best.0.row), round8(best.0.col)); + let mut cur_j = if inside(cur) { cost(cur) } else { f64::MAX }; + if cur_j < best.1 { + best = (cur, cur_j); + } else { + cur = best.0; + cur_j = best.1; + } + for (step, max_iters) in [(16, 1), (8, 2), (4, 3), (2, 4), (1, 8)] { + let s = step * 8; + for _ in 0..max_iters { + let mut moved = false; + let center = cur; + for (dr, dc) in [(-s, 0), (s, 0), (0, -s), (0, s)] { + let cand = Mv::new(center.row + dr, center.col + dc); + if !inside(cand) { + continue; + } + let j = cost(cand); + if j < cur_j { + cur = cand; + cur_j = j; + moved = true; + } + } + if !moved { + break; + } + } + } + if cur_j < best.1 { + best = (cur, cur_j); + } + } + // Sub-sample refinement: half then quarter samples (1/8 units, even: no high precision). + for s in [4, 2] { + let center = best.0; + for (dr, dc) in [(-s, 0), (s, 0), (0, -s), (0, s), (-s, -s), (-s, s), (s, -s), (s, s)] { + let cand = Mv::new(center.row + dr, center.col + dc); + if !inside(cand) { + continue; + } + let j = cost(cand); + if j < best.1 { + best = (cand, j); + } + } + } + best + } + + /// Choose between one 16×16 block and four 8×8 blocks at (`r`, `c`): `Some(mv)` for 16×16. + fn decide16(&mut self, r: usize, c: usize) -> Option { + let (mv16, j16) = self.search(r, c, BLOCK_16X16); + let lambda = self.f.lambda_sad; + // A near-perfect 16×16 prediction needs no split. + if j16 < 2.0 * 256.0 / 256.0 * 16.0 { + return Some(mv16); + } + let mut j8 = lambda * 4.0; // partition split + three more blocks of mode info (≈ 4 bits) + for (k, (dr, dc)) in [(0, 0), (0, 1), (1, 0), (1, 1)].into_iter().enumerate() { + let (rr, cc) = (r + dr, c + dc); + let (mv, j) = self.search(rr, cc, BLOCK_8X8); + j8 += j; + self.pending[k] = (rr, cc, mv); + // Tentative mode info so the next 8×8 blocks see this one as a candidate. + self.set_mi(rr, cc, 1, Mi { inter: true, mode: NEWMV, mv, skip: false, size: BLOCK_8X8 }); + } + if j16 <= j8 { + self.pending = [(usize::MAX, 0, Mv::ZERO); 4]; + Some(mv16) + } else { + None + } + } + + /// Code an inter block with motion vector `mv`. + fn inter_block(&mut self, r: usize, c: usize, bsize: u8, mv: Mv) { + let (mi_rows, mi_cols) = (self.f.mi_rows, self.f.mi_cols); + let n8 = if bsize == BLOCK_16X16 { 2 } else { 1 }; + let n = 8 * n8; + let (nearest, near, mctx) = self.mv_refs(r, c, bsize); + let d = &probs_defaults().0; + let mode_bits = |m: u8| tree_cost(&INTER_MODE_TREE, &d[probs::INTER_MODE + mctx * 3..], m - NEARESTMV); + let mut mode = (NEWMV, mode_bits(NEWMV) + mv_cost(mv.minus(nearest))); + for (m, ok) in [(NEARESTMV, mv == nearest), (NEARMV, mv == near), (ZEROMV, mv == Mv::ZERO)] { + if ok && mode_bits(m) < mode.1 { + mode = (m, mode_bits(m)); + } + } + let mode = mode.0; + // Prediction of every plane (written straight into the reconstruction). + let refs = self.f.refs.as_ref().expect("reference"); + for plane in 0..3 { + let ch = plane > 0; + let (row16, col16) = clamp_for_plane(mv, r, c, n8, mi_rows, mi_cols, ch); + let bn = if ch { n / 2 } else { n }; + let (bx, by) = if ch { (c * 4, r * 4) } else { (c * 8, r * 8) }; + let mut pred = [0u8; 256]; + predict(&refs[plane], ((bx as i32) << 4) + col16, ((by as i32) << 4) + row16, bn, bn, &mut pred); + let p = &mut self.planes[plane]; + for i in 0..bn { + p.rec[(by + i) * p.w + bx..][..bn].copy_from_slice(&pred[i * bn..i * bn + bn]); + } + } + // Residual: quantise every 4×4, dropping blocks whose coefficients don't pay for themselves. + let q = Quant { round_dc: 0.5, round_ac: 0.42, ..self.f.quant }; + // Weight of the rate term when dropping a transform block's coefficients (low: keep quality on par + // with key frames at the same q index). + let zl = 0.07; + let tx = DCT_DCT; + let mut blocks: Vec<(usize, usize, usize, [i32; 16], usize)> = Vec::with_capacity(24); + let mut any = false; + for plane in 0..3 { + let ch = plane > 0; + let bn = if ch { n / 2 } else { n }; + let (bx, by) = if ch { (c * 4, r * 4) } else { (c * 8, r * 8) }; + let p = &self.planes[plane]; + for ty in (0..bn).step_by(4) { + for tx4 in (0..bn).step_by(4) { + let (x, y) = (bx + tx4, by + ty); + let mut res = [0i32; 16]; + for i in 0..4 { + for j in 0..4 { + let k = (y + i) * p.w + x + j; + res[i * 4 + j] = p.src[k] as i32 - p.rec[k] as i32; + } + } + let mut levels = q.levels(&res, tx); + let mut eob = eob_of(&levels, tx); + if eob > 0 && !q.lossless { + let back = q.residual(&levels, tx, eob); + let d_skip: i64 = res.iter().map(|v| (*v as i64) * (*v as i64)).sum(); + let d_code: i64 = res.iter().zip(&back).map(|(a, b)| ((a - b) as i64).pow(2)).sum(); + if d_skip as f64 <= d_code as f64 + zl * self.f.lambda_sse * levels_cost(&levels, eob) as f64 / 256.0 { + levels = [0; 16]; + eob = 0; + } + } + any |= eob > 0; + blocks.push((plane, x, y, levels, eob)); + } + } + } + let skip = !any; + // Mode info. + let avail_u = r > 0; + let avail_l = c > self.col_start; + let above = avail_u.then(|| self.mi_at(r - 1, c)); + let left = avail_l.then(|| self.mi_at(r, c - 1)); + let skip_ctx = above.is_some_and(|m| m.skip) as usize + left.is_some_and(|m| m.skip) as usize; + self.w.ctx(skip, probs::SKIP + skip_ctx); + // is_inter context: neighbours are inter (no intra blocks in inter frames). + let intra = |m: Option| m.is_some_and(|m| !m.inter); + let inter_ctx = match (above, left) { + (Some(_), Some(_)) => { + if intra(above) && intra(left) { + 3 + } else { + (intra(above) || intra(left)) as usize + } + } + (Some(_), None) => 2 * intra(above) as usize, + (None, Some(_)) => 2 * intra(left) as usize, + (None, None) => 0, + }; + self.w.ctx(true, probs::IS_INTER + inter_ctx); + // single_ref_p1 = 0 (LAST_FRAME); context for single-reference LAST neighbours. + let ref_ctx = match (above, left) { + (Some(a), Some(l)) => match (a.inter, l.inter) { + (false, false) => 2, + (false, true) => 4, + (true, false) => 4, + (true, true) => 4, + }, + (Some(m), None) | (None, Some(m)) => { + if m.inter { + 4 + } else { + 2 + } + } + (None, None) => 2, + }; + self.w.ctx(false, probs::SINGLE_REF + ref_ctx * 2); + self.w.tree(&INTER_MODE_TREE, mode - NEARESTMV, |i| P::Ctx((probs::INTER_MODE + mctx * 3 + i) as u16)); + if mode == NEWMV { + self.write_mv(mv.minus(nearest)); + } + self.set_mi(r, c, n8, Mi { inter: true, mode, mv, skip, size: bsize }); + // Tokens and reconstruction. + for (plane, x, y, levels, eob) in blocks { + let ss = (plane > 0) as usize; + let (x4, ly) = (x >> 2, (y >> 2) & (15 >> ss)); + if !skip { + let ctx0 = (self.above_nz[plane][x4] + self.left_nz[plane][ly]) as usize; + self.tokens(plane, true, tx, &levels, eob, ctx0); + if eob > 0 { + let resid = q.residual(&levels, tx, eob); + let p = &mut self.planes[plane]; + for i in 0..4 { + for j in 0..4 { + let k = (y + i) * p.w + x + j; + p.rec[k] = (p.rec[k] as i32 + resid[i * 4 + j]).clamp(0, 255) as u8; + } + } + } + } + let nz = (!skip && eob > 0) as u8; + self.above_nz[plane][x4] = nz; + self.left_nz[plane][ly] = nz; + } + } + + fn write_mv(&mut self, diff: Mv) { + let joint = (diff.row != 0) as u8 * 2 + (diff.col != 0) as u8; + self.w.tree(&MV_JOINT_TREE, joint, |i| P::Ctx((probs::MV_JOINT + i) as u16)); + if diff.row != 0 { + self.write_mv_component(0, diff.row); + } + if diff.col != 0 { + self.write_mv_component(1, diff.col); + } + } + + fn write_mv_component(&mut self, i: usize, v: i32) { + debug_assert!(v != 0 && v & 1 == 0, "vectors are even without high precision"); + let offset = v.unsigned_abs() - 1; + let c = mv_class(offset); + self.w.ctx(v < 0, probs::MV_SIGN + i); + self.w.tree(&MV_CLASS_TREE, c as u8, |k| P::Ctx((probs::MV_CLASS + i * 10 + k) as u16)); + let fr = ((offset >> 1) & 3) as u8; + if c == 0 { + let b = (offset >> 3) as usize; + self.w.ctx(b == 1, probs::MV_CLASS0_BIT + i); + self.w.tree(&MV_FR_TREE, fr, |k| P::Ctx((probs::MV_CLASS0_FR + i * 6 + b * 3 + k) as u16)); + } else { + let dd = (offset - (1 << (c + 3))) >> 3; + for b in 0..c { + self.w.ctx((dd >> b) & 1 == 1, probs::MV_BITS + i * 10 + b as usize); + } + self.w.tree(&MV_FR_TREE, fr, |k| P::Ctx((probs::MV_FR + i * 3 + k) as u16)); + } + } +} diff --git a/crates/vp9enc/src/transform.rs b/crates/vp9enc/src/transform.rs new file mode 100644 index 00000000..e6a685d5 --- /dev/null +++ b/crates/vp9enc/src/transform.rs @@ -0,0 +1,249 @@ +//! 4×4 transforms: the decoder's integer inverse DCT / ADST (spec 8.7.1) and Walsh–Hadamard +//! transform (8.7.1.10), and their forward counterparts (the exact inverse matrices of the +//! decoder's linear transforms, so `inverse4(forward(res)) ≈ res`; the WHT inverse is exact). + +pub const DCT_DCT: u8 = 0; +pub const ADST_DCT: u8 = 1; +pub const DCT_ADST: u8 = 2; +pub const ADST_ADST: u8 = 3; + +fn r14(v: i32) -> i32 { + (v + (1 << 13)) >> 14 +} + +fn idct4(x: [i32; 4]) -> [i32; 4] { + let v0 = r14(x[0] * 11585 - x[2] * 11585); + let v1 = r14(x[0] * 11585 + x[2] * 11585); + let v2 = r14(x[1] * 6270 - x[3] * 15137); + let v3 = r14(x[1] * 15137 + x[3] * 6270); + [v1 + v3, v0 + v2, v0 - v2, v1 - v3] +} + +fn iadst4(x: [i32; 4]) -> [i32; 4] { + let (s19, s29, s39, s49) = (5283i64, 9929i64, 13377i64, 15212i64); + let x: [i64; 4] = x.map(|v| v as i64); + let s0 = s19 * x[0]; + let s1 = s29 * x[0]; + let s2 = s39 * x[1]; + let s3 = s49 * x[2]; + let s4 = s19 * x[2]; + let s5 = s29 * x[3]; + let s6 = s49 * x[3]; + let s7 = s39 * (x[0] - x[2] + x[3]); + let x0 = s0 + s3 + s5; + let x1 = s1 - s4 - s6; + let x2 = s7; + let x3 = s2; + let r = |v: i64| ((v + (1 << 13)) >> 14) as i32; + [r(x0 + x3), r(x1 + x3), r(x2), r(x0 + x1 - x3)] +} + +/// The decoder's 4×4 inverse transform of dequantised coefficients (row-major) to a residual. +pub fn inverse4(coefs: &[i32; 16], tx_type: u8) -> [i32; 16] { + let row_adst = matches!(tx_type, DCT_ADST | ADST_ADST); + let col_adst = matches!(tx_type, ADST_DCT | ADST_ADST); + let mut t = [[0i32; 4]; 4]; + for i in 0..4 { + let r = [coefs[i * 4], coefs[i * 4 + 1], coefs[i * 4 + 2], coefs[i * 4 + 3]]; + t[i] = if row_adst { iadst4(r) } else { idct4(r) }; + } + let mut out = [0i32; 16]; + for j in 0..4 { + let c = [t[0][j], t[1][j], t[2][j], t[3][j]]; + let c = if col_adst { iadst4(c) } else { idct4(c) }; + for i in 0..4 { + out[i * 4 + j] = (c[i] + 8) >> 4; + } + } + out +} + +/// Inverse Walsh–Hadamard (8.7.1.10) on one line. +fn iwht(x: [i32; 4], shift: u32) -> [i32; 4] { + let mut a = x[0] >> shift; + let mut c = x[1] >> shift; + let mut d = x[2] >> shift; + let mut b = x[3] >> shift; + a += c; + d -= b; + let e = (a - d) >> 1; + b = e - b; + c = e - c; + a -= b; + d += c; + [a, b, c, d] +} + +/// The exact inverse of [`iwht`] (with shift 0): inputs that produce outputs `y`. +fn fwht(y: [i32; 4]) -> [i32; 4] { + let (aa, bb, cc, dd) = (y[0], y[1], y[2], y[3]); + let d1 = dd - cc; + let a1 = aa + bb; + let e = (a1 - d1) >> 1; + let b0 = e - bb; + let c0 = e - cc; + let a0 = a1 - c0; + let d0 = d1 + b0; + [a0, c0, d0, b0] +} + +pub fn inverse_wht(coefs: &[i32; 16]) -> [i32; 16] { + let mut t = [[0i32; 4]; 4]; + for i in 0..4 { + t[i] = iwht([coefs[i * 4], coefs[i * 4 + 1], coefs[i * 4 + 2], coefs[i * 4 + 3]], 2); + } + let mut out = [0i32; 16]; + for j in 0..4 { + let r = iwht([t[0][j], t[1][j], t[2][j], t[3][j]], 0); + for i in 0..4 { + out[i * 4 + j] = r[i]; + } + } + out +} + +/// Lossless forward transform: levels (dequantised by 4) whose inverse is exactly `res`. +pub fn forward_wht(res: &[i32; 16]) -> [i32; 16] { + let mut t = [[0i32; 4]; 4]; + for j in 0..4 { + let col = fwht([res[j], res[4 + j], res[8 + j], res[12 + j]]); + for i in 0..4 { + t[i][j] = col[i]; + } + } + let mut lv = [0i32; 16]; + for i in 0..4 { + let r = fwht(t[i]); + lv[i * 4..i * 4 + 4].copy_from_slice(&r); + } + lv +} + +/// Forward transform matrices: the inverses of the decoder's (linear part of the) inverse +/// transforms. +pub fn forward_matrices() -> &'static [[[f64; 16]; 16]; 4] { + use std::sync::OnceLock; + static M: OnceLock<[[[f64; 16]; 16]; 4]> = OnceLock::new(); + M.get_or_init(|| { + let mut out = [[[0.0; 16]; 16]; 4]; + for (ty, o) in out.iter_mut().enumerate() { + // Columns of the inverse: response to a scaled unit coefficient. + let k = 1 << 16; + let mut inv = [[0.0f64; 16]; 16]; + for c in 0..16 { + let mut e = [0i32; 16]; + e[c] = k; + let r = inverse4(&e, ty as u8); + for p in 0..16 { + inv[p][c] = r[p] as f64 / k as f64; + } + } + *o = invert16(inv); + } + out + }) +} + +fn invert16(m: [[f64; 16]; 16]) -> [[f64; 16]; 16] { + let mut a = m; + let mut inv = [[0.0; 16]; 16]; + for (i, r) in inv.iter_mut().enumerate() { + r[i] = 1.0; + } + for col in 0..16 { + let piv = (col..16).max_by(|&x, &y| a[x][col].abs().total_cmp(&a[y][col].abs())).unwrap_or(col); + a.swap(col, piv); + inv.swap(col, piv); + let d = a[col][col]; + for j in 0..16 { + a[col][j] /= d; + inv[col][j] /= d; + } + for r in 0..16 { + if r != col { + let f = a[r][col]; + if f != 0.0 { + for j in 0..16 { + a[r][j] -= f * a[col][j]; + inv[r][j] -= f * inv[col][j]; + } + } + } + } + } + inv +} + +/// Quantisation settings of one block. +#[derive(Clone, Copy)] +pub struct Quant { + pub dcq: i32, + pub acq: i32, + pub lossless: bool, + /// Rounding offsets (fractions of a step) for DC and AC; smaller = wider dead zone. + pub round_dc: f64, + pub round_ac: f64, +} + +impl Quant { + /// Quantised levels (raster order) of a residual. + pub fn levels(&self, res: &[i32; 16], tx_type: u8) -> [i32; 16] { + if self.lossless { + return forward_wht(res); + } + let f = &forward_matrices()[tx_type as usize]; + let mut lv = [0i32; 16]; + for k in 0..16 { + let coef: f64 = (0..16).map(|q| f[k][q] * res[q] as f64).sum(); + let (qs, rnd) = if k == 0 { (self.dcq, self.round_dc) } else { (self.acq, self.round_ac) }; + let l = (coef.abs() / qs as f64 + rnd).floor() as i32; + lv[k] = if coef < 0.0 { -l } else { l }; + } + lv + } + + /// The decoder's residual for `levels` (all zero when `eob` is 0). + pub fn residual(&self, levels: &[i32; 16], tx_type: u8, eob: usize) -> [i32; 16] { + if eob == 0 { + return [0; 16]; + } + let mut deq = [0i32; 16]; + for k in 0..16 { + deq[k] = levels[k] * if k == 0 { self.dcq } else { self.acq }; + } + if self.lossless { inverse_wht(&deq) } else { inverse4(&deq, tx_type) } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn wht_is_exactly_invertible() { + let mut s = 7u32; + for _ in 0..2000 { + let mut r = [0i32; 16]; + for v in &mut r { + s = s.wrapping_mul(1664525).wrapping_add(1013904223); + *v = (s >> 16) as i32 % 511 - 255; + } + let lv = forward_wht(&r); + let deq = lv.map(|l| l * 4); + assert_eq!(inverse_wht(&deq), r); + } + } + + #[test] + fn forward_transforms_invert_the_decoder() { + for ty in 0..4u8 { + let f = &forward_matrices()[ty as usize]; + let r: [i32; 16] = std::array::from_fn(|i| (i as i32 * 37 % 61) - 30); + let c: [i32; 16] = std::array::from_fn(|k| (0..16).map(|q| f[k][q] * r[q] as f64).sum::().round() as i32); + let back = inverse4(&c, ty); + for i in 0..16 { + assert!((back[i] - r[i]).abs() <= 1, "type {ty}: {back:?} vs {r:?}"); + } + } + } +} diff --git a/docs/parity.md b/docs/parity.md index 38df7882..8fc67b9a 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -110,7 +110,7 @@ the disk cache in the formats wave; the 3D Camera Tracker in M12.6.) | Area | Weighted parity | Remaining (agent-hours) | Biggest gaps | |---|---|---|---| | Layers | 88% | 4.8 | frame blending, collapse transformations, slip edit | -| Output | 90% | 3.0 | multiple output modules, pre-render (WebM with VP9 alpha + Opus and WAV/AIFF audio-only landed: VP9 is intra-only, Opus CELT-only) | +| Output | 90% | 3.0 | multiple output modules, pre-render (WebM with VP9 alpha + Opus and WAV/AIFF audio-only landed: VP9 key + inter frames with motion search, loop filter and rate control; Opus CELT-only) | | Audio | 85% | 0.5 | audio to keyframes | | Import | 88% | 2.0 | AI/EPS/PDF vector footage, PSD smart objects and 3D layers (PSD as footage/composition/retain layer sizes, SVG footage and Create Shapes from Vector Layer landed) | | Automation | ≈ 90% | 0.8 | scripting covers the documented core object model (AUT-2, M14.4: `app`, project items, comps, layers, properties and keyframes, text documents, markers, render queue, Script Console, `effectcraft-cli script`, MCP `run_script`); still missing: ScriptUI panels/dialogs, `.jsxbin`, sockets |