Oregami
Repositories/oxedyne/fe2o3

oxedyne/fe2o3/fe2o3_graphics/src/hevc/transform.rs

23.2 KiB, 65 runs

created by r1870400018:20494, which is this file's identity for as long as the history lasts, whatever it is later renamed to

download · who wrote it · its history

1//! Turning coded coefficients back into the residual that is added to a prediction.
2//!
3//! Two steps, and both are integer arithmetic specified to the bit -- a decoder that rounds
4//! differently from the specification does not produce a slightly different picture, it produces
5//! one that drifts further from the encoder's with every block that predicts from it.
6//!
7//! **Scaling** (§8.6.3) undoes the quantiser: each coefficient is multiplied by a factor chosen by
8//! the quantisation parameter, shifted left by that parameter divided by six, and rounded back down
9//! by a shift that depends on the bit depth and the block size. The six factors are the sixth roots
10//! of two to within a per cent, which is why the parameter's sixth part is a doubling.
11//!
12//! **The inverse transform** (§8.6.4) is a matrix multiplication down each column and then along
13//! each row, with a rounding shift between the two. Two matrices exist: a four-point sine transform
14//! used for the luma of four-by-four intra blocks, whose basis functions suit a residual that grows
15//! away from the predicted edge, and the cosine transform everything else uses. The cosine matrix is
16//! one thirty-two by thirty-two table for every size -- a sixteen-point transform reads every second
17//! column of it, an eight-point every fourth, and so on -- which is what makes a single
18//! implementation of the one-dimensional pass serve all four sizes.
19//!
20//! A block coded without a transform at all (`transform_skip_flag`) takes a rotation and a shift
21//! instead, and one coded without the quantiser either (`cu_transquant_bypass_flag`) is the
22//! residual already.
23//!
24//! [Written with AI entirely](https://need2know.ai/entirely-ai/code)\
25//! Anthropic Claude
26
27use oxedyne_fe2o3_core::prelude::*;
28
29pub const MAX_TB: usize = 32; // the largest transform this decoder will do, samples each way
30
31// What multiplies a coefficient before the shift, by the quantisation parameter's remainder on
32// division by six (§8.6.3). The six of them span one doubling: 40, 45, 51, 57, 64 and 72 are the
33// sixth powers of two to within a per cent, so six steps of the parameter double the step size.
34const LEVEL_SCALE: [i32; 6] = [40, 45, 51, 57, 64, 72];
35
36// The four-point sine transform, for the luma of an intra block of four (equation 8-316).
37const DST_4: [[i16; 4]; 4] = [
38 [29, 55, 74, 84],
39 [74, 74, 0, -74],
40 [84, -29, -74, 55],
41 [55, -84, 74, -29],
42];
43
44// Columns 0 to 15 of the transform matrix (§8.6.4.2, equation 8-319).
45//
46// The published table is transposed against the way the equation indexes it: a printed row is the
47// matrix's second subscript and a printed column its first, so transMatrix[m][n] is this array's
48// [n][m]. matrix() is the only place that knows it, and a check against the four-point inverse
49// everybody knows by heart is what settles that it is the right way round.
50const DCT_COL_0_15: [[i16; 16]; 32] = [
51 [64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64],
52 [90, 90, 88, 85, 82, 78, 73, 67, 61, 54, 46, 38, 31, 22, 13, 4],
53 [90, 87, 80, 70, 57, 43, 25, 9, -9, -25, -43, -57, -70, -80, -87, -90],
54 [90, 82, 67, 46, 22, -4, -31, -54, -73, -85, -90, -88, -78, -61, -38, -13],
55 [89, 75, 50, 18, -18, -50, -75, -89, -89, -75, -50, -18, 18, 50, 75, 89],
56 [88, 67, 31, -13, -54, -82, -90, -78, -46, -4, 38, 73, 90, 85, 61, 22],
57 [87, 57, 9, -43, -80, -90, -70, -25, 25, 70, 90, 80, 43, -9, -57, -87],
58 [85, 46, -13, -67, -90, -73, -22, 38, 82, 88, 54, -4, -61, -90, -78, -31],
59 [83, 36, -36, -83, -83, -36, 36, 83, 83, 36, -36, -83, -83, -36, 36, 83],
60 [82, 22, -54, -90, -61, 13, 78, 85, 31, -46, -90, -67, 4, 73, 88, 38],
61 [80, 9, -70, -87, -25, 57, 90, 43, -43, -90, -57, 25, 87, 70, -9, -80],
62 [78, -4, -82, -73, 13, 85, 67, -22, -88, -61, 31, 90, 54, -38, -90, -46],
63 [75, -18, -89, -50, 50, 89, 18, -75, -75, 18, 89, 50, -50, -89, -18, 75],
64 [73, -31, -90, -22, 78, 67, -38, -90, -13, 82, 61, -46, -88, -4, 85, 54],
65 [70, -43, -87, 9, 90, 25, -80, -57, 57, 80, -25, -90, -9, 87, 43, -70],
66 [67, -54, -78, 38, 85, -22, -90, 4, 90, 13, -88, -31, 82, 46, -73, -61],
67 [64, -64, -64, 64, 64, -64, -64, 64, 64, -64, -64, 64, 64, -64, -64, 64],
68 [61, -73, -46, 82, 31, -88, -13, 90, -4, -90, 22, 85, -38, -78, 54, 67],
69 [57, -80, -25, 90, -9, -87, 43, 70, -70, -43, 87, 9, -90, 25, 80, -57],
70 [54, -85, -4, 88, -46, -61, 82, 13, -90, 38, 67, -78, -22, 90, -31, -73],
71 [50, -89, 18, 75, -75, -18, 89, -50, -50, 89, -18, -75, 75, 18, -89, 50],
72 [46, -90, 38, 54, -90, 31, 61, -88, 22, 67, -85, 13, 73, -82, 4, 78],
73 [43, -90, 57, 25, -87, 70, 9, -80, 80, -9, -70, 87, -25, -57, 90, -43],
74 [38, -88, 73, -4, -67, 90, -46, -31, 85, -78, 13, 61, -90, 54, 22, -82],
75 [36, -83, 83, -36, -36, 83, -83, 36, 36, -83, 83, -36, -36, 83, -83, 36],
76 [31, -78, 90, -61, 4, 54, -88, 82, -38, -22, 73, -90, 67, -13, -46, 85],
77 [25, -70, 90, -80, 43, 9, -57, 87, -87, 57, -9, -43, 80, -90, 70, -25],
78 [22, -61, 85, -90, 73, -38, -4, 46, -78, 90, -82, 54, -13, -31, 67, -88],
79 [18, -50, 75, -89, 89, -75, 50, -18, -18, 50, -75, 89, -89, 75, -50, 18],
80 [13, -38, 61, -78, 88, -90, 85, -73, 54, -31, 4, 22, -46, 67, -82, 90],
81 [9, -25, 43, -57, 70, -80, 87, -90, 90, -87, 80, -70, 57, -43, 25, -9],
82 [4, -13, 22, -31, 38, -46, 54, -61, 67, -73, 78, -82, 85, -88, 90, -90],
83];
84
85// Columns 16 to 31 of the same matrix (equation 8-321), indexed the same way.
86const DCT_COL_16_31: [[i16; 16]; 32] = [
87 [64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64],
88 [-4, -13, -22, -31, -38, -46, -54, -61, -67, -73, -78, -82, -85, -88, -90, -90],
89 [-90, -87, -80, -70, -57, -43, -25, -9, 9, 25, 43, 57, 70, 80, 87, 90],
90 [13, 38, 61, 78, 88, 90, 85, 73, 54, 31, 4, -22, -46, -67, -82, -90],
91 [89, 75, 50, 18, -18, -50, -75, -89, -89, -75, -50, -18, 18, 50, 75, 89],
92 [-22, -61, -85, -90, -73, -38, 4, 46, 78, 90, 82, 54, 13, -31, -67, -88],
93 [-87, -57, -9, 43, 80, 90, 70, 25, -25, -70, -90, -80, -43, 9, 57, 87],
94 [31, 78, 90, 61, 4, -54, -88, -82, -38, 22, 73, 90, 67, 13, -46, -85],
95 [83, 36, -36, -83, -83, -36, 36, 83, 83, 36, -36, -83, -83, -36, 36, 83],
96 [-38, -88, -73, -4, 67, 90, 46, -31, -85, -78, -13, 61, 90, 54, -22, -82],
97 [-80, -9, 70, 87, 25, -57, -90, -43, 43, 90, 57, -25, -87, -70, 9, 80],
98 [46, 90, 38, -54, -90, -31, 61, 88, 22, -67, -85, -13, 73, 82, 4, -78],
99 [75, -18, -89, -50, 50, 89, 18, -75, -75, 18, 89, 50, -50, -89, -18, 75],
100 [-54, -85, 4, 88, 46, -61, -82, 13, 90, 38, -67, -78, 22, 90, 31, -73],
101 [-70, 43, 87, -9, -90, -25, 80, 57, -57, -80, 25, 90, 9, -87, -43, 70],
102 [61, 73, -46, -82, 31, 88, -13, -90, -4, 90, 22, -85, -38, 78, 54, -67],
103 [64, -64, -64, 64, 64, -64, -64, 64, 64, -64, -64, 64, 64, -64, -64, 64],
104 [-67, -54, 78, 38, -85, -22, 90, 4, -90, 13, 88, -31, -82, 46, 73, -61],
105 [-57, 80, 25, -90, 9, 87, -43, -70, 70, 43, -87, -9, 90, -25, -80, 57],
106 [73, 31, -90, 22, 78, -67, -38, 90, -13, -82, 61, 46, -88, 4, 85, -54],
107 [50, -89, 18, 75, -75, -18, 89, -50, -50, 89, -18, -75, 75, 18, -89, 50],
108 [-78, -4, 82, -73, -13, 85, -67, -22, 88, -61, -31, 90, -54, -38, 90, -46],
109 [-43, 90, -57, -25, 87, -70, -9, 80, -80, 9, 70, -87, 25, 57, -90, 43],
110 [82, -22, -54, 90, -61, -13, 78, -85, 31, 46, -90, 67, 4, -73, 88, -38],
111 [36, -83, 83, -36, -36, 83, -83, 36, 36, -83, 83, -36, -36, 83, -83, 36],
112 [-85, 46, 13, -67, 90, -73, 22, 38, -82, 88, -54, -4, 61, -90, 78, -31],
113 [-25, 70, -90, 80, -43, -9, 57, -87, 87, -57, 9, 43, -80, 90, -70, 25],
114 [88, -67, 31, 13, -54, 82, -90, 78, -46, 4, 38, -73, 90, -85, 61, -22],
115 [18, -50, 75, -89, 89, -75, 50, -18, -18, 50, -75, 89, -89, 75, -50, 18],
116 [-90, 82, -67, 46, -22, -4, 31, -54, 73, -85, 90, -88, 78, -61, 38, -13],
117 [-9, 25, -43, 57, -70, 80, -87, 90, -90, 87, -80, 70, -57, 43, -25, 9],
118 [90, -90, 88, -85, 82, -78, 73, -67, 61, -54, 46, -38, 31, -22, 13, -4],
119];
120/// One entry of the cosine transform matrix, undoing the transposition of the published tables.
121///
122/// `m` is the equation's first subscript and `n` its second, both nought to thirty-one.
123fn matrix(m: usize, n: usize) -> i32 {
124 if m < 16 {
125 DCT_COL_0_15[n][m] as i32
126 } else {
127 DCT_COL_16_31[n][m - 16] as i32
128 }
129}
130
131/// Which transform a block takes.
132#[derive(Clone, Copy, Debug, PartialEq, Eq)]
133pub enum Kind {
134 Sine, // four-point sine: intra luma, four by four, nothing else (§8.6.4.1)
135 Cosine, // the cosine transform, at whatever size the block is
136}
137
138impl Kind {
139
140 pub fn of(intra: bool, size: usize, chroma: bool) -> Self {
141 if intra && size == 4 && !chroma {
142 Self::Sine
143 } else {
144 Self::Cosine
145 }
146 }
147}
148
149/// Undoes the quantiser (§8.6.3), in place.
150///
151/// `coeffs` is the block in raster order, `size` its side, `qp` the quantisation parameter that
152/// applies to it, `depth` the bit depth of the component and `m` the scaling matrix, which is
153/// sixteen everywhere unless the sequence says otherwise.
154pub fn scale(coeffs: &mut [i32], size: usize, qp: i32, depth: u32, m: &[i32]) {
155 // Fifteen is the transform range every profile this decoder meets uses; the extended-precision
156 // flag that widens it belongs to profiles that do not appear in a photograph.
157 let shift = depth as i32 + log2(size) as i32 + 10 - 15;
158 let scale = LEVEL_SCALE[(qp % 6) as usize];
159 let up = qp / 6;
160 let round = 1i64 << (shift - 1);
161 for (i, c) in coeffs.iter_mut().take(size * size).enumerate() {
162 let factor = m.get(i).copied().unwrap_or(16) as i64;
163 // In sixty-four bits because the intermediate overflows thirty-two: a coefficient may be
164 // fifteen bits, the scale seven, and the shift up to eight more.
165 let v = (((*c as i64) * factor * (scale as i64)) << up) + round;
166 *c = ((v >> shift).clamp(-32_768, 32_767)) as i32;
167 }
168}
169
170// The default scaling lists (§7.4.5, Table 7-6), which is what a sequence that turns the lists on
171// without carrying any of its own means.
172//
173// Sixty-four values for an eight-by-eight block, in the diagonal scan's order; the sixteen and
174// thirty-two sample matrices are this one with each value covering two or four samples each way,
175// and the four-sample matrix is flat. The numbers climb away from the corner because the eye
176// notices an error in the coarse detail of a block more than in the fine, so the fine detail is
177// quantised harder.
178//
179// The first row is for a block predicted from within the picture and the second for one predicted
180// from another picture, which a still photograph never is -- it is here because the two are one
181// table in the document and splitting them would invite the wrong one being used.
182pub const DEFAULT_LIST: [[u8; 64]; 2] = [
183 [
184 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 17, 16, 17, 16, 17, 18,
185 17, 18, 18, 17, 18, 21, 19, 20, 21, 20, 19, 21, 24, 22, 22, 24,
186 24, 22, 22, 24, 25, 25, 27, 30, 27, 25, 25, 29, 31, 35, 35, 31,
187 29, 36, 41, 44, 41, 36, 47, 54, 54, 47, 65, 70, 65, 88, 88, 115,
188 ],
189 [
190 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 18,
191 18, 18, 18, 18, 18, 20, 20, 20, 20, 20, 20, 20, 24, 24, 24, 24,
192 24, 24, 24, 24, 25, 25, 25, 25, 25, 25, 25, 28, 28, 28, 28, 28,
193 28, 33, 33, 33, 33, 33, 41, 41, 41, 41, 54, 54, 54, 71, 71, 91,
194 ],
195];
196
197/// The default eight-by-eight scaling matrix, in raster order.
198///
199/// The published list is in the diagonal scan's order, so laying it out takes the same scan the
200/// coefficients themselves are read in.
201pub fn default_matrix(inter: bool) -> [u8; 64] {
202 let mut out = [16u8; 64];
203 let list = DEFAULT_LIST[inter as usize];
204 for (i, (x, y)) in crate::hevc::scan::positions(8, crate::hevc::scan::Order::Diagonal)
205 .iter()
206 .enumerate()
207 {
208 out[*y as usize * 8 + *x as usize] = list[i];
209 }
210 out
211}
212
213/// The base-two logarithm of a power of two.
214fn log2(n: usize) -> u32 {
215 n.trailing_zeros()
216}
217
218/// One dimension of the inverse transform (§8.6.4.2).
219///
220/// `src` holds `size` coefficients and `dst` takes `size` samples. The cosine matrix is read with a
221/// stride, which is what lets one table of thirty-two serve every size.
222fn one_way(src: &[i32], dst: &mut [i32], size: usize, kind: Kind) {
223 match kind {
224 Kind::Sine => for i in 0..4 {
225 let mut sum = 0i64;
226 for j in 0..4 {
227 sum += (DST_4[j][i] as i64) * (src[j] as i64);
228 }
229 dst[i] = sum as i32;
230 },
231 Kind::Cosine => {
232 let stride = 32 / size;
233 for i in 0..size {
234 let mut sum = 0i64;
235 for j in 0..size {
236 sum += (matrix(i, j * stride) as i64) * (src[j] as i64);
237 }
238 dst[i] = sum as i32;
239 }
240 },
241 }
242}
243
244/// The whole inverse transform: columns, a rounding shift, then rows (§8.6.4.1).
245///
246/// `block` arrives holding scaled coefficients in raster order and leaves holding the transform's
247/// output, which is **not yet the residual** -- [`finish`] applies the last shift, and it applies it
248/// to a skipped block too, which is what keeps the two paths in the same units.
249///
250/// The clip between the two passes is the specification's and is load-bearing: doing both passes
251/// and rounding once at the end is arithmetically similar and produces a different picture.
252pub fn inverse(block: &mut [i32], size: usize, kind: Kind) -> Outcome<()> {
253 if size > MAX_TB || !size.is_power_of_two() || size < 4 {
254 return Err(err!("A transform of {} samples was asked for.", size; Invalid, Input));
255 }
256 if block.len() < size * size {
257 return Err(err!(
258 "A transform of {0} wants {1} coefficients and was given {2}.",
259 size, size * size, block.len(); Invalid, Input));
260 }
261 let mut col = [0i32; MAX_TB];
262 let mut out = [0i32; MAX_TB];
263 // Down each column.
264 for x in 0..size {
265 for y in 0..size {
266 col[y] = block[y * size + x];
267 }
268 one_way(&col[..size], &mut out[..size], size, kind);
269 for y in 0..size {
270 // Clipped to sixteen bits, which is what keeps the second pass inside the range its
271 // matrix was designed for.
272 block[y * size + x] = ((out[y] + 64) >> 7).clamp(-32_768, 32_767);
273 }
274 }
275 // Then along each row, whose output the specification does not clip.
276 for y in 0..size {
277 col[..size].copy_from_slice(&block[y * size..y * size + size]);
278 one_way(&col[..size], &mut out[..size], size, kind);
279 block[y * size..y * size + size].copy_from_slice(&out[..size]);
280 }
281 Ok(())
282}
283
284/// A block coded without its transform (§8.6.2, equation 8-298): one shift left.
285///
286/// It stands where [`inverse`] would, and [`finish`] follows it just the same.
287pub fn skipped(block: &mut [i32], size: usize) {
288 let shift = 5 + log2(size);
289 for v in block.iter_mut().take(size * size) {
290 *v <<= shift;
291 }
292}
293
294/// The last shift, which turns either path's output into residual samples (equation 8-299).
295pub fn finish(block: &mut [i32], size: usize, depth: u32) {
296 let shift = 20 - depth as i32;
297 if shift <= 0 {
298 return;
299 }
300 let round = 1i32 << (shift - 1);
301 for v in block.iter_mut().take(size * size) {
302 *v = (*v + round) >> shift;
303 }
304}
305
306#[cfg(test)]
307mod tests {
308 use super::*;
309
310 /// The four-point inverse transforms, written out as the arithmetic everybody who has
311 /// implemented one knows by heart, rather than as a table lookup.
312 ///
313 /// This is what settles the orientation of the published matrices, which are printed
314 /// transposed against the way the equation subscripts them. A decoder that reads them the other
315 /// way round still produces a picture -- a wrong one, in a way no amount of staring at the table
316 /// reveals.
317 fn known_dct_4(x: [i32; 4]) -> [i32; 4] {
318 [
319 64 * x[0] + 83 * x[1] + 64 * x[2] + 36 * x[3],
320 64 * x[0] + 36 * x[1] - 64 * x[2] - 83 * x[3],
321 64 * x[0] - 36 * x[1] - 64 * x[2] + 83 * x[3],
322 64 * x[0] - 83 * x[1] + 64 * x[2] - 36 * x[3],
323 ]
324 }
325
326 fn known_dst_4(x: [i32; 4]) -> [i32; 4] {
327 [
328 29 * x[0] + 74 * x[1] + 84 * x[2] + 55 * x[3],
329 55 * x[0] + 74 * x[1] - 29 * x[2] - 84 * x[3],
330 74 * x[0] - 74 * x[2] + 74 * x[3],
331 84 * x[0] - 74 * x[1] + 55 * x[2] - 29 * x[3],
332 ]
333 }
334
335 #[test]
336 fn test_the_four_point_transforms_are_the_ones_everybody_knows_00() -> Outcome<()> {
337 let cases: [[i32; 4]; 5] = [
338 [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [17, -9, 240, -1000],
339 ];
340 for x in cases {
341 let mut got = [0i32; 4];
342 one_way(&x, &mut got, 4, Kind::Cosine);
343 req!(got, known_dct_4(x), "the cosine matrix is the wrong way round for {:?}", x);
344 let mut got = [0i32; 4];
345 one_way(&x, &mut got, 4, Kind::Sine);
346 req!(got, known_dst_4(x), "the sine matrix is the wrong way round for {:?}", x);
347 }
348 Ok(())
349 }
350
351 #[test]
352 fn test_every_size_of_the_cosine_matrix_is_orthogonal_01() -> Outcome<()> {
353 // The property that makes it a transform at all, and one this decoder's own arithmetic
354 // cannot fake: the rows are at right angles to each other and all the same length. A
355 // transposition, a wrong stride or a mistyped entry breaks it, and none of those is visible
356 // by reading the table.
357 //
358 // The integer matrix approximates an orthogonal one, so the off-diagonal terms are not
359 // quite nought. Measured against the diagonal, the worst of them is one part in 550 at
360 // sixteen points, 682 at thirty-two, 1,310 at eight and 2,340 at four, and the worst error
361 // in a row's own length is one part in 1,560. The bounds below sit just outside those. A
362 // single mistyped entry is worth far more than either: the smallest entry in the table is 4
363 // against a diagonal of 16,384, and getting one wrong moves a dot product by hundreds.
364 for size in [4usize, 8, 16, 32] {
365 let stride = 32 / size;
366 let mut worst = 0i64;
367 let diagonal = 64i64 * 64 * size as i64;
368 for i in 0..size {
369 for k in 0..size {
370 let mut dot = 0i64;
371 for j in 0..size {
372 dot += (matrix(i, j * stride) as i64) * (matrix(k, j * stride) as i64);
373 }
374 if i == k {
375 // Every row the same length, which is what keeps the picture's gain flat.
376 let off = (dot - diagonal).abs();
377 let near = off * 1_000 <= diagonal;
378 req!(near, true,
379 "row {} of the {}-point matrix has length squared {}, not {}",
380 i, size, dot, diagonal);
381 } else if dot.abs() > worst {
382 worst = dot.abs();
383 }
384 }
385 }
386 let square = worst * 500 <= diagonal;
387 req!(square, true,
388 "two rows of the {}-point matrix are {} out of square, against a diagonal of {}",
389 size, worst, diagonal);
390 }
391 Ok(())
392 }
393
394 #[test]
395 fn test_a_flat_block_comes_out_of_one_coefficient_02() -> Outcome<()> {
396 // What a decoder does far more often than anything else: a block whose only coefficient is
397 // the direct one, which must come back as a single level with no pattern in it. A stride
398 // worked out wrongly puts ripples in this, and ripples in a flat sky are exactly the
399 // artefact somebody reports.
400 for size in [4usize, 8, 16, 32] {
401 let mut block = vec![0i32; size * size];
402 block[0] = 512;
403 // Not the sine transform: that one is only ever used on four-by-four intra luma, where
404 // its basis is deliberately not flat.
405 res!(inverse(&mut block, size, Kind::Cosine));
406 finish(&mut block, size, 8);
407 let first = block[0];
408 for (i, v) in block.iter().enumerate() {
409 req!(*v, first, "sample {} of a {}-point flat block is not flat", i, size);
410 }
411 let something = first != 0;
412 req!(something, true, "a {}-point block of 512 came back empty", size);
413 }
414 Ok(())
415 }
416
417 #[test]
418 fn test_the_quantiser_steps_by_the_sixth_root_of_two_03() -> Outcome<()> {
419 // The published factors, against the arithmetic they approximate. Six steps of the
420 // quantisation parameter is meant to be one doubling of the step size, so consecutive
421 // factors should stand in the ratio of the sixth root of two -- and the sixth one, stepped
422 // once more, should land on twice the first. Nothing in this decoder is consulted; the
423 // oracle is the arithmetic the table was built from.
424 let root = 2f64.powf(1.0 / 6.0);
425 for k in 0..5 {
426 let ratio = LEVEL_SCALE[k + 1] as f64 / LEVEL_SCALE[k] as f64;
427 let close = (ratio - root).abs() < 0.02;
428 req!(close, true,
429 "levelScale {} to {} is a ratio of {:.4}, and the sixth root of two is {:.4}",
430 LEVEL_SCALE[k], LEVEL_SCALE[k + 1], ratio, root);
431 }
432 let wrapped = LEVEL_SCALE[5] as f64 * root;
433 let octave = (wrapped - 2.0 * LEVEL_SCALE[0] as f64).abs() < 2.0;
434 req!(octave, true,
435 "one step past the last factor is {:.1}, and twice the first is {}",
436 wrapped, 2 * LEVEL_SCALE[0]);
437
438 // And the scaling itself follows the table: six steps up doubles what comes out, to within
439 // the rounding the shift cannot avoid.
440 // A coefficient of one, so that the highest parameter still lands well inside the
441 // sixteen-bit range the scaled coefficients are clipped to -- a saturated value would
442 // double into itself and prove nothing.
443 for qp in 0..46i32 {
444 let mut low = [1i32; 16];
445 let mut high = [1i32; 16];
446 scale(&mut low, 4, qp, 8, &[16; 16]);
447 scale(&mut high, 4, qp + 6, 8, &[16; 16]);
448 let doubled = (high[0] - low[0] * 2).abs() <= 1;
449 req!(doubled, true,
450 "qp {} scales to {} and qp {} to {}", qp, low[0], qp + 6, high[0]);
451 }
452 Ok(())
453 }
454
455 #[test]
456 fn test_a_skipped_block_lands_in_the_same_units_04() -> Outcome<()> {
457 // A block coded without its transform takes a shift instead, and the two paths have to
458 // arrive in the same units or a picture that mixes them is a picture with a step in it.
459 // A flat block through the transform and the same block through the shift agree.
460 for size in [4usize, 8] {
461 let mut through = vec![0i32; size * size];
462 // The coefficient that yields a flat block of one at this size.
463 through[0] = 1;
464 res!(inverse(&mut through, size, Kind::Cosine));
465 finish(&mut through, size, 8);
466
467 let mut around = vec![0i32; size * size];
468 around[0] = 1;
469 skipped(&mut around, size);
470 finish(&mut around, size, 8);
471 // The transform spreads the one over the block and the shift leaves it in the corner,
472 // so what is compared is the level, not the position.
473 req!(around[0], through[0] * (size as i32),
474 "the two paths disagree at {} by more than the transform's own gain", size);
475 }
476 Ok(())
477 }
478
479 #[test]
480 fn test_the_matrix_is_the_published_one_05() -> Outcome<()> {
481 // The same discipline the context tables are held to: a thousand and twenty-four numbers
482 // copied out of a document, checked against the document.
483 //
484 // HEVC_SPEC_TEXT=~/.cache/specs/h265.txt cargo test -p oxedyne_fe2o3_graphics hevc
485 let path = match std::env::var("HEVC_SPEC_TEXT") {
486 Ok(p) => p,
487 Err(_) => {
488 println!(" skipped: set HEVC_SPEC_TEXT to a text rendering of Rec. ITU-T H.265");
489 return Ok(());
490 },
491 };
492 let text = match std::fs::read_to_string(&path) {
493 Ok(t) => t,
494 Err(e) => {
495 println!(" skipped: {} would not read ({})", path, e);
496 return Ok(());
497 },
498 };
499 for (marker, held) in [
500 ("transMatrixCol0to15 =", &DCT_COL_0_15),
501 ("transMatrixCol16to31 =", &DCT_COL_16_31),
502 ] {
503 let at = match text.find(marker) {
504 Some(at) => at,
505 None => return Err(err!("{} is not in {}.", marker, path; Test, Missing)),
506 };
507 let mut rows: Vec<Vec<i32>> = Vec::new();
508 for line in text[at..].lines() {
509 let t = line.trim();
510 if !t.starts_with('{') || !t.trim_end_matches(',').ends_with('}') {
511 continue;
512 }
513 let inner = t.trim_end_matches(',').trim_start_matches('{').trim_end_matches('}');
514 let mut row = Vec::new();
515 let mut ok = true;
516 for word in inner.split_whitespace() {
517 // The document writes a minus sign as U+2212, not as a hyphen.
518 match word.replace('\u{2212}', "-").parse::<i32>() {
519 Ok(n) => row.push(n),
520 Err(_) => { ok = false; break; },
521 }
522 }
523 if ok && row.len() == 16 {
524 rows.push(row);
525 }
526 if rows.len() == 32 {
527 break;
528 }
529 }
530 req!(rows.len(), 32usize, "{} does not have thirty-two rows under it", marker);
531 for (n, row) in rows.iter().enumerate() {
532 for (m, v) in row.iter().enumerate() {
533 req!(held[n][m] as i32, *v,
534 "{} row {} column {} is {} and the document says {}",
535 marker, n, m, held[n][m], v);
536 }
537 }
538 }
539 Ok(())
540 }
541}