diff --git a/HEIF_IMPLEMENTATION_PLAN.md b/HEIF_IMPLEMENTATION_PLAN.md index a05065e906..2930d9312f 100644 --- a/HEIF_IMPLEMENTATION_PLAN.md +++ b/HEIF_IMPLEMENTATION_PLAN.md @@ -453,7 +453,7 @@ Exit gate: - [ ] Independently encoded, opaque, single-item AVIF files reconstruct correctly across all AVIF profiles, bit depths, subsampling modes, and normative still-image compression tools. Pixel comparisons are made after applying the same signaled color conversion in the reference path. -### Phase 4: complete scalar HEVC still-image reconstruction +### Phase 4: complete HEVC still-image reconstruction Implement and verify in dependency order: @@ -466,11 +466,15 @@ Implement and verify in dependency order: - [x] Implement reconstructed-plane reference collection, substitution, asymmetric chroma-unit handling, and reference-filter selection from caller-owned availability and scratch spans. - [ ] Build each availability flag from reconstructed-plane and coding state, and connect the predictor to transform-unit traversal. - [ ] Scaling lists, inverse quantization, transform skip, every required inverse transform, range-extension precision, and lossless reconstruction. + - [x] Implement allocation-free SIMD-first inverse DCT for every 4/8/16/32 width and height combination, the four-by-four intra-luma inverse DST, normative intermediate and residual clipping, and saturated prediction addition for 8/10/12-bit samples. Verify normal and forced-scalar execution against fixed results and a dense scalar oracle. + - [ ] Implement scaling-list selection and inverse quantization for every transform size, component, prediction mode, quantization parameter, and range-extension precision path. + - [ ] Implement transform skip, coefficient rotation, implicit and explicit residual DPCM, transquant bypass, and lossless reconstruction. + - [ ] Connect coefficient decoding, inverse quantization, transform selection, reusable scratch, and add/clip to transform-unit traversal. - [ ] Deblocking and sample-adaptive offset for every signaled luma/chroma and bit-depth path. - [ ] Tiles, wavefront entry points, dependent slices, and all other parallelization syntax permitted by the exposed still-image profiles. - [ ] Supplemental enhancement information that changes image presentation or metadata exposed by ImageSharp. -Each subsystem starts with scalar conformance vectors derived from the HEVC specification and the pinned implementation reference. No valid syntax in the exposed HEIC profiles may terminate in an unsupported branch or silently omit a normative reconstruction stage. +Each subsystem keeps a specification-shaped scalar oracle derived from the HEVC specification and the pinned implementation reference, while production hot paths are designed SIMD-first with scalar fallback. No valid syntax in the exposed HEIC profiles may terminate in an unsupported branch or silently omit a normative reconstruction stage. Exit gate: @@ -552,7 +556,7 @@ Exit gate: ### Phase 8: SIMD and allocation optimization -SIMD work begins after the corresponding scalar stage has independent correctness tests; it then proceeds alongside subsequent functional phases. +Production hot paths are designed SIMD-first. Each SIMD implementation is verified against an independent specification-shaped scalar oracle before it is connected to the decoder pipeline. Tasks: @@ -561,6 +565,7 @@ Tasks: - [ ] Benchmark codec-specific costs for CABAC/range decode, inverse transforms, still-image prediction, deblocking, SAO, CDEF, restoration, chroma upsampling, color conversion, alpha packing, and grid copies. - [x] Add a permanent frame-wide HEVC intra-prediction benchmark. On .NET 10, SIMD planar, vertical-angular, and horizontal-angular prediction measured 155.2, 150.7, and 255.0 microseconds per padded 1920x1088 frame, compared with forced-scalar timings of 1.890, 1.081, and 0.979 milliseconds: 12.2, 7.2, and 3.8 times faster with zero managed allocations. - [x] Add a permanent frame-wide HEVC reference-preparation benchmark. On .NET 10, complete and partially substituted borders measured 465.3 and 566.7 microseconds per 2,040-block padded 1920x1088 frame, compared with forced-scalar timings of 477.5 and 608.0 microseconds, with zero managed allocations. + - [x] Add a permanent frame-wide HEVC inverse-transform benchmark. On .NET 10, dense-coefficient 32x32 twelve-bit inverse DCT, transposition, and add/clip measured 3.643 milliseconds per padded 1920x1088 frame, compared with 45.23 milliseconds with hardware intrinsics disabled: 12.4 times faster with zero managed allocations. - [ ] Implement vector paths only for confirmed hot loops, using existing `Vector128`, `Vector256`, and `Vector512` helper and dispatch patterns where supported. - [ ] Prioritize shared color conversion and pixel packing, chroma upsampling, inverse-transform add-and-clip, intra predictors, HEVC deblock/SAO, AV1 loop filter/CDEF/restoration, and contiguous grid copies. - [ ] Benchmark the complete decode color pipeline on representative 8/10/12-bit AVIF and HEIC images with and without embedded ICC profiles. Report absolute end-to-end timings and allocations in addition to the isolated YUV/CICP and ICC stage costs. @@ -631,12 +636,12 @@ The critical path is: - [ ] Baseline and provenance. - [ ] Public HEIC/AVIF boundaries and the shared parsed HEIF container model. -- [ ] Scalar AV1 still decode, starting with 8-bit 4:2:0 and completing every AVIF profile/tool. -- [ ] Scalar HEVC still decode, starting with 8-bit 4:2:0 and completing every exposed HEIC profile/tool. +- [ ] SIMD-first AV1 still decode with specification-shaped scalar oracles, starting with 8-bit 4:2:0 and completing every AVIF profile/tool. +- [ ] SIMD-first HEVC still decode with specification-shaped scalar oracles, starting with 8-bit 4:2:0 and completing every exposed HEIC profile/tool. - [ ] Shared color, alpha, grids, metadata, and presentation transforms across all bit depths and chroma formats. - [ ] Real AV1 encoder and libavif/libaom cross-decode across the complete matrix. - [ ] Real HEVC encoder and independent HEIC cross-decode across the complete matrix. -- [ ] Measured SIMD/allocation work integrated after each scalar subsystem stabilizes. +- [ ] Measured SIMD/allocation work integrated with each production hot-path subsystem before that subsystem is considered complete. - [ ] Hardening and release gates. This sequence does not define partial PR completion. The early 8-bit 4:2:0 decoders and smallest-valid encoders are temporary vertical slices, but the PR remains incomplete until all nine steps and the complete still-image verification matrix pass. diff --git a/src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.Operations.cs b/src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.Operations.cs new file mode 100644 index 0000000000..573008bc5c --- /dev/null +++ b/src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.Operations.cs @@ -0,0 +1,523 @@ +// Copyright (c) Six Labors. +// Licensed under the Six Labors Split License. + +using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; +using System.Runtime.Intrinsics; +using SixLabors.ImageSharp.Formats.Heif.Av1.Transform; + +namespace SixLabors.ImageSharp.Formats.Heif.Hevc; + +internal static partial class HevcInverseTransformer +{ + /// + /// Calculates the disjoint odd-frequency groups that seed the HEVC partial-butterfly reconstruction. + /// + /// The selected inverse-DCT operator. + /// The frequency rows followed by contiguous independent lines. + /// The destination group rows. + /// The number of independent lines transformed together. + private static void PopulateButterflyGroups(ReadOnlySpan source, Span groups, int lineCount) + where TOperator : struct, IHevcInverseTransformOperator + { + int size = TOperator.Size; + int groupOffset = 0; + for (int frequencyStep = 2; frequencyStep < size; frequencyStep <<= 1) + { + int outputCount = size / frequencyStep; + int firstFrequency = frequencyStep >> 1; + for (int position = 0; position < outputCount; position++) + { + PopulateButterflyGroupRow( + source, + groups.Slice((groupOffset + position) * lineCount, lineCount), + lineCount, + firstFrequency, + frequencyStep, + position); + } + + groupOffset += outputCount; + } + + // The deepest even group contains the DC term and the transform's Nyquist-frequency term. It remains a + // two-element group for every supported DCT size and closes the recursive butterfly hierarchy. + for (int position = 0; position < 2; position++) + { + PopulateButterflyGroupRow( + source, + groups.Slice((groupOffset + position) * lineCount, lineCount), + lineCount, + 0, + size >> 1, + position); + } + } + + /// + /// Calculates one partial-butterfly group row across all independent lines. + /// + /// The selected inverse-DCT operator. + /// The complete frequency-row input. + /// The destination group row. + /// The number of independent lines. + /// The first frequency included by the group. + /// The distance between included frequencies. + /// The group-relative spatial coordinate. + private static void PopulateButterflyGroupRow( + ReadOnlySpan source, + Span destination, + int lineCount, + int firstFrequency, + int frequencyStep, + int position) + where TOperator : struct, IHevcInverseTransformOperator + { + ref int sourceBase = ref MemoryMarshal.GetReference(source); + ref int destinationBase = ref MemoryMarshal.GetReference(destination); + int x = 0; + + if (Vector512.IsHardwareAccelerated) + { + int oneVectorFromEnd = lineCount - Vector512.Count; + for (; x <= oneVectorFromEnd; x += Vector512.Count) + { + Vector512 sum = Vector512.Zero; + for (int frequency = firstFrequency; frequency < TOperator.Size; frequency += frequencyStep) + { + Vector512 values = Vector512.LoadUnsafe(ref sourceBase, (nuint)((frequency * lineCount) + x)); + sum += values * Vector512.Create(TOperator.GetCoefficient(frequency, position)); + } + + sum.StoreUnsafe(ref destinationBase, (nuint)x); + } + } + + if (Vector256.IsHardwareAccelerated) + { + int oneVectorFromEnd = lineCount - Vector256.Count; + for (; x <= oneVectorFromEnd; x += Vector256.Count) + { + Vector256 sum = Vector256.Zero; + for (int frequency = firstFrequency; frequency < TOperator.Size; frequency += frequencyStep) + { + Vector256 values = Vector256.LoadUnsafe(ref sourceBase, (nuint)((frequency * lineCount) + x)); + sum += values * Vector256.Create(TOperator.GetCoefficient(frequency, position)); + } + + sum.StoreUnsafe(ref destinationBase, (nuint)x); + } + } + + if (Vector128.IsHardwareAccelerated) + { + int oneVectorFromEnd = lineCount - Vector128.Count; + for (; x <= oneVectorFromEnd; x += Vector128.Count) + { + Vector128 sum = Vector128.Zero; + for (int frequency = firstFrequency; frequency < TOperator.Size; frequency += frequencyStep) + { + Vector128 values = Vector128.LoadUnsafe(ref sourceBase, (nuint)((frequency * lineCount) + x)); + sum += values * Vector128.Create(TOperator.GetCoefficient(frequency, position)); + } + + sum.StoreUnsafe(ref destinationBase, (nuint)x); + } + } + + for (; x < lineCount; x++) + { + int sum = 0; + for (int frequency = firstFrequency; frequency < TOperator.Size; frequency += frequencyStep) + { + sum += source[(frequency * lineCount) + x] * TOperator.GetCoefficient(frequency, position); + } + + destination[x] = sum; + } + } + + /// + /// Expands the disjoint partial-butterfly groups into spatial rows. + /// + /// The selected inverse-DCT operator. + /// The buffer containing every disjoint group. + /// The alternate expansion buffer. + /// The number of independent lines transformed together. + /// The rounded right shift applied at the final hierarchy level. + /// The inclusive output minimum. + /// The inclusive output maximum. + /// The buffer containing the completed spatial rows. + private static Span CombineButterflyGroups( + Span initial, + Span alternate, + int lineCount, + int shift, + int minimum, + int maximum) + where TOperator : struct, IHevcInverseTransformOperator + { + int size = TOperator.Size; + int combinedSize = 2; + int combinedStart = size - combinedSize; + bool currentIsInitial = true; + while (combinedSize < size) + { + int oddStart = combinedStart - combinedSize; + bool finalLevel = (combinedSize << 1) == size; + ReadOnlySpan even = currentIsInitial ? initial : alternate; + Span destination = currentIsInitial ? alternate : initial; + + for (int position = 0; position < combinedSize; position++) + { + ReadOnlySpan evenRow = even.Slice((combinedStart + position) * lineCount, lineCount); + ReadOnlySpan oddRow = initial.Slice((oddStart + position) * lineCount, lineCount); + Span positiveRow = destination.Slice((oddStart + position) * lineCount, lineCount); + Span negativeRow = destination.Slice((oddStart + (2 * combinedSize) - 1 - position) * lineCount, lineCount); + CombineButterflyRows(evenRow, oddRow, positiveRow, negativeRow, finalLevel, shift, minimum, maximum); + } + + combinedStart = oddStart; + combinedSize <<= 1; + currentIsInitial = !currentIsInitial; + } + + return currentIsInitial ? initial : alternate; + } + + /// + /// Combines one symmetric pair of partial-butterfly rows. + /// + /// The even-frequency contribution. + /// The odd-frequency contribution. + /// The destination receiving the added contribution. + /// The destination receiving the subtracted contribution. + /// Whether this is the final hierarchy level. + /// The rounded right shift applied at the final hierarchy level. + /// The inclusive output minimum. + /// The inclusive output maximum. + private static void CombineButterflyRows( + ReadOnlySpan even, + ReadOnlySpan odd, + Span positive, + Span negative, + bool roundAndClip, + int shift, + int minimum, + int maximum) + { + ref int evenBase = ref MemoryMarshal.GetReference(even); + ref int oddBase = ref MemoryMarshal.GetReference(odd); + ref int positiveBase = ref MemoryMarshal.GetReference(positive); + ref int negativeBase = ref MemoryMarshal.GetReference(negative); + int x = 0; + + if (Vector512.IsHardwareAccelerated) + { + int oneVectorFromEnd = even.Length - Vector512.Count; + for (; x <= oneVectorFromEnd; x += Vector512.Count) + { + Vector512 evenValues = Vector512.LoadUnsafe(ref evenBase, (nuint)x); + Vector512 oddValues = Vector512.LoadUnsafe(ref oddBase, (nuint)x); + Vector512 added = evenValues + oddValues; + Vector512 subtracted = evenValues - oddValues; + if (roundAndClip) + { + added = RoundShiftAndClamp(added, shift, minimum, maximum); + subtracted = RoundShiftAndClamp(subtracted, shift, minimum, maximum); + } + + added.StoreUnsafe(ref positiveBase, (nuint)x); + subtracted.StoreUnsafe(ref negativeBase, (nuint)x); + } + } + + if (Vector256.IsHardwareAccelerated) + { + int oneVectorFromEnd = even.Length - Vector256.Count; + for (; x <= oneVectorFromEnd; x += Vector256.Count) + { + Vector256 evenValues = Vector256.LoadUnsafe(ref evenBase, (nuint)x); + Vector256 oddValues = Vector256.LoadUnsafe(ref oddBase, (nuint)x); + Vector256 added = evenValues + oddValues; + Vector256 subtracted = evenValues - oddValues; + if (roundAndClip) + { + added = RoundShiftAndClamp(added, shift, minimum, maximum); + subtracted = RoundShiftAndClamp(subtracted, shift, minimum, maximum); + } + + added.StoreUnsafe(ref positiveBase, (nuint)x); + subtracted.StoreUnsafe(ref negativeBase, (nuint)x); + } + } + + if (Vector128.IsHardwareAccelerated) + { + int oneVectorFromEnd = even.Length - Vector128.Count; + for (; x <= oneVectorFromEnd; x += Vector128.Count) + { + Vector128 evenValues = Vector128.LoadUnsafe(ref evenBase, (nuint)x); + Vector128 oddValues = Vector128.LoadUnsafe(ref oddBase, (nuint)x); + Vector128 added = evenValues + oddValues; + Vector128 subtracted = evenValues - oddValues; + if (roundAndClip) + { + added = RoundShiftAndClamp(added, shift, minimum, maximum); + subtracted = RoundShiftAndClamp(subtracted, shift, minimum, maximum); + } + + added.StoreUnsafe(ref positiveBase, (nuint)x); + subtracted.StoreUnsafe(ref negativeBase, (nuint)x); + } + } + + for (; x < even.Length; x++) + { + int added = even[x] + odd[x]; + int subtracted = even[x] - odd[x]; + if (roundAndClip) + { + added = RoundShiftAndClamp(added, shift, minimum, maximum); + subtracted = RoundShiftAndClamp(subtracted, shift, minimum, maximum); + } + + positive[x] = added; + negative[x] = subtracted; + } + } + + /// + /// Applies the four-point inverse-DST matrix across all independent lines. + /// + /// The selected dense transform operator. + /// The frequency rows followed by contiguous independent lines. + /// The spatial rows followed by contiguous independent lines. + /// The number of independent lines transformed together. + /// The rounded right shift applied to each result. + /// The inclusive output minimum. + /// The inclusive output maximum. + private static void TransformDense( + ReadOnlySpan source, + Span destination, + int lineCount, + int shift, + int minimum, + int maximum) + where TOperator : struct, IHevcInverseTransformOperator + { + ref int sourceBase = ref MemoryMarshal.GetReference(source); + ref int destinationBase = ref MemoryMarshal.GetReference(destination); + for (int position = 0; position < TOperator.Size; position++) + { + int x = 0; + if (Vector512.IsHardwareAccelerated) + { + int oneVectorFromEnd = lineCount - Vector512.Count; + for (; x <= oneVectorFromEnd; x += Vector512.Count) + { + Vector512 sum = Vector512.Zero; + for (int frequency = 0; frequency < TOperator.Size; frequency++) + { + Vector512 values = Vector512.LoadUnsafe(ref sourceBase, (nuint)((frequency * lineCount) + x)); + sum += values * Vector512.Create(TOperator.GetCoefficient(frequency, position)); + } + + RoundShiftAndClamp(sum, shift, minimum, maximum).StoreUnsafe(ref destinationBase, (nuint)((position * lineCount) + x)); + } + } + + if (Vector256.IsHardwareAccelerated) + { + int oneVectorFromEnd = lineCount - Vector256.Count; + for (; x <= oneVectorFromEnd; x += Vector256.Count) + { + Vector256 sum = Vector256.Zero; + for (int frequency = 0; frequency < TOperator.Size; frequency++) + { + Vector256 values = Vector256.LoadUnsafe(ref sourceBase, (nuint)((frequency * lineCount) + x)); + sum += values * Vector256.Create(TOperator.GetCoefficient(frequency, position)); + } + + RoundShiftAndClamp(sum, shift, minimum, maximum).StoreUnsafe(ref destinationBase, (nuint)((position * lineCount) + x)); + } + } + + if (Vector128.IsHardwareAccelerated) + { + int oneVectorFromEnd = lineCount - Vector128.Count; + for (; x <= oneVectorFromEnd; x += Vector128.Count) + { + Vector128 sum = Vector128.Zero; + for (int frequency = 0; frequency < TOperator.Size; frequency++) + { + Vector128 values = Vector128.LoadUnsafe(ref sourceBase, (nuint)((frequency * lineCount) + x)); + sum += values * Vector128.Create(TOperator.GetCoefficient(frequency, position)); + } + + RoundShiftAndClamp(sum, shift, minimum, maximum).StoreUnsafe(ref destinationBase, (nuint)((position * lineCount) + x)); + } + } + + for (; x < lineCount; x++) + { + int sum = 0; + for (int frequency = 0; frequency < TOperator.Size; frequency++) + { + sum += source[(frequency * lineCount) + x] * TOperator.GetCoefficient(frequency, position); + } + + destination[(position * lineCount) + x] = RoundShiftAndClamp(sum, shift, minimum, maximum); + } + } + } + + /// + /// Transposes one rectangular transform block into a separate full-block buffer. + /// + /// The source block in raster order. + /// The transposed destination block. + /// The source row count. + /// The source column count. + private static void Transpose(ReadOnlySpan source, Span destination, int sourceHeight, int sourceWidth) + { + if (Vector128.IsHardwareAccelerated) + { + ref int sourceBase = ref MemoryMarshal.GetReference(source); + ref int destinationBase = ref MemoryMarshal.GetReference(destination); + for (int y = 0; y < sourceHeight; y += 4) + { + for (int x = 0; x < sourceWidth; x += 4) + { + Vector128 row0 = Vector128.LoadUnsafe(ref sourceBase, (nuint)((y * sourceWidth) + x)); + Vector128 row1 = Vector128.LoadUnsafe(ref sourceBase, (nuint)(((y + 1) * sourceWidth) + x)); + Vector128 row2 = Vector128.LoadUnsafe(ref sourceBase, (nuint)(((y + 2) * sourceWidth) + x)); + Vector128 row3 = Vector128.LoadUnsafe(ref sourceBase, (nuint)(((y + 3) * sourceWidth) + x)); + Av1Transform2dOperations.Transpose(ref row0, ref row1, ref row2, ref row3); + row0.StoreUnsafe(ref destinationBase, (nuint)((x * sourceHeight) + y)); + row1.StoreUnsafe(ref destinationBase, (nuint)(((x + 1) * sourceHeight) + y)); + row2.StoreUnsafe(ref destinationBase, (nuint)(((x + 2) * sourceHeight) + y)); + row3.StoreUnsafe(ref destinationBase, (nuint)(((x + 3) * sourceHeight) + y)); + } + } + + return; + } + + for (int y = 0; y < sourceHeight; y++) + { + for (int x = 0; x < sourceWidth; x++) + { + destination[(x * sourceHeight) + y] = source[(y * sourceWidth) + x]; + } + } + } + + /// + /// Adds a complete signed residual block to the predicted samples and clips to the component precision. + /// + /// The signed residual block in raster order. + /// The predicted samples beginning at the block origin. + /// The destination row stride in samples. + /// The transform-block width. + /// The transform-block height. + /// The reconstructed component precision. + private static void AddResidual(ReadOnlySpan residual, Span destination, int destinationStride, int width, int height, int bitDepth) + { + int maximum = (1 << bitDepth) - 1; + for (int y = 0; y < height; y++) + { + ReadOnlySpan residualRow = residual.Slice(y * width, width); + Span destinationRow = destination.Slice(y * destinationStride, width); + ref int residualBase = ref MemoryMarshal.GetReference(residualRow); + ref ushort destinationBase = ref MemoryMarshal.GetReference(destinationRow); + int x = 0; + + if (Vector256.IsHardwareAccelerated) + { + int oneVectorFromEnd = width - Vector256.Count; + for (; x <= oneVectorFromEnd; x += Vector256.Count) + { + Vector128 predicted16 = Vector128.LoadUnsafe(ref destinationBase, (nuint)x); + Vector256 predicted = Vector256.Create(Vector128.WidenLower(predicted16), Vector128.WidenUpper(predicted16)).AsInt32(); + Vector256 reconstructed = Vector256.Clamp( + predicted + Vector256.LoadUnsafe(ref residualBase, (nuint)x), + Vector256.Zero, + Vector256.Create(maximum)); + + Vector128.Narrow(reconstructed.GetLower().AsUInt32(), reconstructed.GetUpper().AsUInt32()).StoreUnsafe(ref destinationBase, (nuint)x); + } + } + + if (Vector128.IsHardwareAccelerated) + { + int oneVectorFromEnd = width - Vector128.Count; + for (; x <= oneVectorFromEnd; x += Vector128.Count) + { + ulong packed = Unsafe.ReadUnaligned(ref Unsafe.As(ref Unsafe.Add(ref destinationBase, x))); + Vector128 predicted = Vector128.WidenLower(Vector128.CreateScalarUnsafe(packed).AsUInt16()).AsInt32(); + Vector128 reconstructed = Vector128.Clamp( + predicted + Vector128.LoadUnsafe(ref residualBase, (nuint)x), + Vector128.Zero, + Vector128.Create(maximum)); + + Vector128 narrowed = Vector128.Narrow(reconstructed.AsUInt32(), Vector128.Zero); + Unsafe.WriteUnaligned(ref Unsafe.As(ref Unsafe.Add(ref destinationBase, x)), narrowed.AsUInt64().ToScalar()); + } + } + + for (; x < width; x++) + { + destinationRow[x] = (ushort)Math.Clamp(destinationRow[x] + residualRow[x], 0, maximum); + } + } + } + + /// + /// Applies HEVC's rounded right shift and inclusive clipping to a 512-bit vector. + /// + /// The unnormalized transform values. + /// The right-shift count. + /// The inclusive result minimum. + /// The inclusive result maximum. + /// The normalized and clipped values. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static Vector512 RoundShiftAndClamp(Vector512 value, int shift, int minimum, int maximum) + => Vector512.Clamp((value + Vector512.Create(1 << (shift - 1))) >> shift, Vector512.Create(minimum), Vector512.Create(maximum)); + + /// + /// Applies HEVC's rounded right shift and inclusive clipping to a 256-bit vector. + /// + /// The unnormalized transform values. + /// The right-shift count. + /// The inclusive result minimum. + /// The inclusive result maximum. + /// The normalized and clipped values. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static Vector256 RoundShiftAndClamp(Vector256 value, int shift, int minimum, int maximum) + => Vector256.Clamp((value + Vector256.Create(1 << (shift - 1))) >> shift, Vector256.Create(minimum), Vector256.Create(maximum)); + + /// + /// Applies HEVC's rounded right shift and inclusive clipping to a 128-bit vector. + /// + /// The unnormalized transform values. + /// The right-shift count. + /// The inclusive result minimum. + /// The inclusive result maximum. + /// The normalized and clipped values. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static Vector128 RoundShiftAndClamp(Vector128 value, int shift, int minimum, int maximum) + => Vector128.Clamp((value + Vector128.Create(1 << (shift - 1))) >> shift, Vector128.Create(minimum), Vector128.Create(maximum)); + + /// + /// Applies HEVC's rounded right shift and inclusive clipping to one scalar value. + /// + /// The unnormalized transform value. + /// The right-shift count. + /// The inclusive result minimum. + /// The inclusive result maximum. + /// The normalized and clipped value. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int RoundShiftAndClamp(int value, int shift, int minimum, int maximum) + => Math.Clamp((value + (1 << (shift - 1))) >> shift, minimum, maximum); +} diff --git a/src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.OperatorImplementations.cs b/src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.OperatorImplementations.cs new file mode 100644 index 0000000000..4f13162326 --- /dev/null +++ b/src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.OperatorImplementations.cs @@ -0,0 +1,141 @@ +// Copyright (c) Six Labors. +// Licensed under the Six Labors Split License. + +using System.Runtime.CompilerServices; + +namespace SixLabors.ImageSharp.Formats.Heif.Hevc; + +internal static partial class HevcInverseTransformer +{ + /// + /// Gets the common HEVC inverse-DCT magnitudes ordered on the pi-over-sixty-four angle grid. + /// + private static ReadOnlySpan DiscreteCosineMagnitudes => + [ + 90, 90, 90, 90, 89, 88, 87, 85, 83, 82, 80, 78, 75, 73, 70, 67, 64, + 61, 57, 54, 50, 46, 43, 38, 36, 31, 25, 22, 18, 13, 9, 4, 0 + ]; + + /// + /// Gets the four-point HEVC inverse-DST matrix in frequency-major order. + /// + private static ReadOnlySpan DiscreteSine4Coefficients => + [ + 29, 55, 74, 84, + 74, 74, 0, -74, + 84, -29, -74, 55, + 55, -84, 74, -29 + ]; + + /// + /// Reconstructs one normative inverse-DCT coefficient from the common HEVC angle grid. + /// + /// The transform side. + /// The frequency-domain coordinate. + /// The spatial-domain coordinate. + /// The signed transform coefficient. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetDiscreteCosineCoefficient(int size, int frequency, int position) + { + if (frequency == 0) + { + return 64; + } + + int angle = ((2 * position) + 1) * frequency * (32 / size); + angle &= 127; + if (angle > 64) + { + angle = 128 - angle; + } + + // Angles in the second quadrant reuse the first-quadrant magnitude with a negative sign. Keeping the + // normative integer magnitudes in one compile-time span avoids runtime trigonometry and duplicate matrices. + if (angle > 32) + { + return -DiscreteCosineMagnitudes[64 - angle]; + } + + return DiscreteCosineMagnitudes[angle]; + } + + /// + /// Implements the four-point inverse discrete cosine transform. + /// + private readonly struct DiscreteCosine4Operator : IHevcInverseTransformOperator + { + /// + public static int Size => 4; + + /// + public static bool UsesButterfly => true; + + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static int GetCoefficient(int frequency, int position) => GetDiscreteCosineCoefficient(Size, frequency, position); + } + + /// + /// Implements the four-point inverse discrete sine transform. + /// + private readonly struct DiscreteSine4Operator : IHevcInverseTransformOperator + { + /// + public static int Size => 4; + + /// + public static bool UsesButterfly => false; + + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static int GetCoefficient(int frequency, int position) => DiscreteSine4Coefficients[(frequency * Size) + position]; + } + + /// + /// Implements the eight-point inverse discrete cosine transform. + /// + private readonly struct DiscreteCosine8Operator : IHevcInverseTransformOperator + { + /// + public static int Size => 8; + + /// + public static bool UsesButterfly => true; + + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static int GetCoefficient(int frequency, int position) => GetDiscreteCosineCoefficient(Size, frequency, position); + } + + /// + /// Implements the sixteen-point inverse discrete cosine transform. + /// + private readonly struct DiscreteCosine16Operator : IHevcInverseTransformOperator + { + /// + public static int Size => 16; + + /// + public static bool UsesButterfly => true; + + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static int GetCoefficient(int frequency, int position) => GetDiscreteCosineCoefficient(Size, frequency, position); + } + + /// + /// Implements the thirty-two-point inverse discrete cosine transform. + /// + private readonly struct DiscreteCosine32Operator : IHevcInverseTransformOperator + { + /// + public static int Size => 32; + + /// + public static bool UsesButterfly => true; + + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static int GetCoefficient(int frequency, int position) => GetDiscreteCosineCoefficient(Size, frequency, position); + } +} diff --git a/src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.cs b/src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.cs new file mode 100644 index 0000000000..4a350fcb1c --- /dev/null +++ b/src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.cs @@ -0,0 +1,270 @@ +// Copyright (c) Six Labors. +// Licensed under the Six Labors Split License. + +using System.Runtime.CompilerServices; +using SixLabors.ImageSharp.Common.Helpers; + +namespace SixLabors.ImageSharp.Formats.Heif.Hevc; + +/// +/// Applies HEVC inverse transforms and reconstructs predicted samples. +/// +internal static partial class HevcInverseTransformer +{ + /// + /// The signed residual precision used after the second inverse-transform pass. + /// + private const int ResidualPrecision = 16; + + /// + /// Defines one closed inverse-transform operation selected by the transform-unit syntax. + /// + /// The implementing operator type. + private interface IHevcInverseTransformOperator + where TOperator : struct, IHevcInverseTransformOperator + { + /// + /// Gets the transform side in samples. + /// + public static abstract int Size { get; } + + /// + /// Gets a value indicating whether the transform uses the partial-butterfly factorization. + /// + public static abstract bool UsesButterfly { get; } + + /// + /// Gets one inverse-transform matrix coefficient. + /// + /// The frequency-domain coordinate. + /// The spatial-domain coordinate. + /// The signed transform coefficient. + public static abstract int GetCoefficient(int frequency, int position); + } + + /// + /// Gets the scratch length required for the specified rectangular transform block. + /// + /// The base-two logarithm of the transform-block width. + /// The base-two logarithm of the transform-block height. + /// The required number of signed thirty-two-bit elements. + public static int GetScratchLength(int log2Width, int log2Height) + { + DebugGuard.MustBeBetweenOrEqualTo(log2Width, 2, 5, nameof(log2Width)); + DebugGuard.MustBeBetweenOrEqualTo(log2Height, 2, 5, nameof(log2Height)); + return 2 << (log2Width + log2Height); + } + + /// + /// Reconstructs one transform block by adding its inverse-transformed residual to the predicted samples. + /// + /// The dequantized transform coefficients in raster order. + /// The predicted samples beginning at the transform-block origin. + /// The destination row stride in samples. + /// The base-two logarithm of the transform-block width. + /// The base-two logarithm of the transform-block height. + /// The reconstructed component precision. + /// The transform dynamic range excluding its sign bit. + /// Whether the four-by-four luma intra block uses the discrete sine transform. + /// The caller-owned scratch returned by . + public static void TransformAdd( + ReadOnlySpan coefficients, + Span destination, + int destinationStride, + int log2Width, + int log2Height, + int bitDepth, + int maxTransformDynamicRange, + bool useDiscreteSineTransform, + Span scratch) + { + DebugGuard.MustBeBetweenOrEqualTo(log2Width, 2, 5, nameof(log2Width)); + DebugGuard.MustBeBetweenOrEqualTo(log2Height, 2, 5, nameof(log2Height)); + DebugGuard.MustBeBetweenOrEqualTo(bitDepth, 8, 16, nameof(bitDepth)); + DebugGuard.IsTrue(!useDiscreteSineTransform || (log2Width == 2 && log2Height == 2), "The HEVC inverse DST is defined only for four-by-four blocks."); + + int width = 1 << log2Width; + int height = 1 << log2Height; + int sampleCount = width * height; + DebugGuard.IsTrue(coefficients.Length >= sampleCount, "The coefficient span is shorter than the transform block."); + DebugGuard.IsTrue(scratch.Length >= sampleCount * 2, "The scratch span is shorter than the inverse-transform requirement."); + + Span first = scratch[..sampleCount]; + Span second = scratch.Slice(sampleCount, sampleCount); + Span residual = TransformCore( + coefficients[..sampleCount], + first, + second, + width, + height, + bitDepth, + maxTransformDynamicRange, + useDiscreteSineTransform); + + AddResidual(residual, destination, destinationStride, width, height, bitDepth); + } + + /// + /// Applies one two-dimensional inverse transform and writes signed residual samples in raster order. + /// + /// The dequantized transform coefficients in raster order. + /// The destination residual samples in raster order. + /// The base-two logarithm of the transform-block width. + /// The base-two logarithm of the transform-block height. + /// The reconstructed component precision. + /// The transform dynamic range excluding its sign bit. + /// Whether the four-by-four luma intra block uses the discrete sine transform. + /// The caller-owned scratch returned by . + public static void Transform( + ReadOnlySpan coefficients, + Span residual, + int log2Width, + int log2Height, + int bitDepth, + int maxTransformDynamicRange, + bool useDiscreteSineTransform, + Span scratch) + { + int width = 1 << log2Width; + int height = 1 << log2Height; + int sampleCount = width * height; + DebugGuard.IsTrue(residual.Length >= sampleCount, "The residual span is shorter than the transform block."); + DebugGuard.IsTrue(scratch.Length >= sampleCount * 2, "The scratch span is shorter than the inverse-transform requirement."); + + Span transformed = TransformCore( + coefficients[..sampleCount], + scratch[..sampleCount], + scratch.Slice(sampleCount, sampleCount), + width, + height, + bitDepth, + maxTransformDynamicRange, + useDiscreteSineTransform); + + transformed.CopyTo(residual); + } + + /// + /// Dispatches both separable transform passes through closed operators selected from the block dimensions. + /// + /// The complete dequantized coefficient block. + /// The first full-block scratch buffer. + /// The second full-block scratch buffer. + /// The transform-block width. + /// The transform-block height. + /// The reconstructed component precision. + /// The transform dynamic range excluding its sign bit. + /// Whether both four-point passes use the discrete sine transform. + /// The scratch buffer containing the raster-ordered residual. + private static Span TransformCore( + ReadOnlySpan coefficients, + Span first, + Span second, + int width, + int height, + int bitDepth, + int maxTransformDynamicRange, + bool useDiscreteSineTransform) + { + int dynamicMinimum = -(1 << maxTransformDynamicRange); + int dynamicMaximum = (1 << maxTransformDynamicRange) - 1; + + // HEVC moves one normalization bit from the second pass to the first. The intermediate clip therefore + // belongs after the vertical pass and must not be combined with the final residual clipping operation. + Span vertical = TransformDimension( + coefficients, + first, + second, + height, + width, + 7, + dynamicMinimum, + dynamicMaximum, + useDiscreteSineTransform); + + Span horizontalInput = vertical.Overlaps(first) ? second : first; + Transpose(vertical, horizontalInput, height, width); + + Span horizontalWorkspace = horizontalInput.Overlaps(first) ? second : first; + int secondShift = maxTransformDynamicRange + 5 - bitDepth; + Span horizontal = TransformDimension( + horizontalInput, + horizontalWorkspace, + horizontalInput, + width, + height, + secondShift, + -(1 << (ResidualPrecision - 1)), + (1 << (ResidualPrecision - 1)) - 1, + useDiscreteSineTransform); + + Span residual = horizontal.Overlaps(first) ? second : first; + Transpose(horizontal, residual, width, height); + return residual; + } + + /// + /// Selects the statically specialized operator for one transform dimension. + /// + /// The frequency rows followed by contiguous independent lines. + /// The initial operator output buffer. + /// The alternate combination buffer. + /// The transform dimension. + /// The number of independent lines transformed together. + /// The rounded right shift applied to the spatial results. + /// The inclusive output minimum. + /// The inclusive output maximum. + /// Whether the four-point pass uses the discrete sine transform. + /// The buffer containing spatial rows followed by contiguous independent lines. + private static Span TransformDimension( + ReadOnlySpan source, + Span initial, + Span alternate, + int size, + int lineCount, + int shift, + int minimum, + int maximum, + bool useDiscreteSineTransform) + => (size, useDiscreteSineTransform) switch + { + (4, true) => TransformDimension(source, initial, alternate, lineCount, shift, minimum, maximum), + (4, false) => TransformDimension(source, initial, alternate, lineCount, shift, minimum, maximum), + (8, _) => TransformDimension(source, initial, alternate, lineCount, shift, minimum, maximum), + (16, _) => TransformDimension(source, initial, alternate, lineCount, shift, minimum, maximum), + _ => TransformDimension(source, initial, alternate, lineCount, shift, minimum, maximum) + }; + + /// + /// Invokes one statically selected inverse-transform operator. + /// + /// The selected inverse-transform operator. + /// The frequency rows followed by contiguous independent lines. + /// The initial operator output buffer. + /// The alternate combination buffer. + /// The number of independent lines transformed together. + /// The rounded right shift applied to the spatial results. + /// The inclusive output minimum. + /// The inclusive output maximum. + /// The buffer containing spatial rows followed by contiguous independent lines. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static Span TransformDimension( + ReadOnlySpan source, + Span initial, + Span alternate, + int lineCount, + int shift, + int minimum, + int maximum) + where TOperator : struct, IHevcInverseTransformOperator + { + if (!TOperator.UsesButterfly) + { + TransformDense(source, initial, lineCount, shift, minimum, maximum); + return initial; + } + + PopulateButterflyGroups(source, initial, lineCount); + return CombineButterflyGroups(initial, alternate, lineCount, shift, minimum, maximum); + } +} diff --git a/src/ImageSharp/Formats/Heif/Hevc/HevcSequenceParameterSet.cs b/src/ImageSharp/Formats/Heif/Hevc/HevcSequenceParameterSet.cs index 1c70351b82..1f8e5585db 100644 --- a/src/ImageSharp/Formats/Heif/Hevc/HevcSequenceParameterSet.cs +++ b/src/ImageSharp/Formats/Heif/Hevc/HevcSequenceParameterSet.cs @@ -525,6 +525,17 @@ internal sealed class HevcSequenceParameterSet /// public bool CabacBypassAlignmentEnabled { get; } + /// + /// Gets the base-two logarithm of the transform dynamic range for the specified reconstructed plane. + /// + /// The reconstructed plane. + /// The transform dynamic range excluding its sign bit. + public int GetMaxTransformDynamicRange(HevcPlane plane) + { + int bitDepth = plane == HevcPlane.Y ? this.BitDepthLuma : this.BitDepthChroma; + return this.ExtendedPrecisionProcessingEnabled ? Math.Max(15, bitDepth + 6) : 15; + } + /// /// Reads a conformance-window offset and converts it to luma-sample units. /// diff --git a/tests/ImageSharp.Benchmarks/Codecs/Heif/HevcInverseTransformBenchmarks.cs b/tests/ImageSharp.Benchmarks/Codecs/Heif/HevcInverseTransformBenchmarks.cs new file mode 100644 index 0000000000..2195b82c90 --- /dev/null +++ b/tests/ImageSharp.Benchmarks/Codecs/Heif/HevcInverseTransformBenchmarks.cs @@ -0,0 +1,90 @@ +// Copyright (c) Six Labors. +// Licensed under the Six Labors Split License. + +using BenchmarkDotNet.Attributes; +using SixLabors.ImageSharp.Formats.Heif.Hevc; + +namespace SixLabors.ImageSharp.Benchmarks.Codecs.Heif; + +/// +/// Measures complete coded-frame traversal through HEVC inverse-transform reconstruction. +/// +[MemoryDiagnoser(displayGenColumns: false)] +public class HevcInverseTransformBenchmarks +{ + /// + /// The coded frame width, which is an exact multiple of the maximum transform-block side. + /// + private const int Width = 1920; + + /// + /// The coded frame height including the final padded coding-tree row for a 1080-line presentation. + /// + private const int Height = 1088; + + /// + /// The base-two logarithm of the benchmark transform-block side. + /// + private const int BlockLog2 = 5; + + /// + /// The transform-block side in samples. + /// + private const int BlockSize = 1 << BlockLog2; + + /// + /// The deterministic dequantized coefficients reused by each benchmark block. + /// + private readonly int[] coefficients = new int[BlockSize * BlockSize]; + + /// + /// The frame-wide predicted and reconstructed samples. + /// + private readonly ushort[] destination = new ushort[Width * Height]; + + /// + /// The maximum-block inverse-transform scratch reused throughout each coded frame. + /// + private readonly int[] scratch = new int[HevcInverseTransformer.GetScratchLength(BlockLog2, BlockLog2)]; + + /// + /// Populates a dense, deterministic twelve-bit transform workload outside the measured frame traversal. + /// + [GlobalSetup] + public void Setup() + { + for (int i = 0; i < this.coefficients.Length; i++) + { + this.coefficients[i] = (((i * 37) + 11) % 127) - 63; + } + + this.destination.AsSpan().Fill(2048); + } + + /// + /// Measures factorized inverse DCT, both transpositions, and saturated prediction addition for a coded frame. + /// + /// The final reconstructed sample, keeping the frame output observable. + [Benchmark] + public ushort TransformFrame() + { + for (int y = 0; y < Height; y += BlockSize) + { + for (int x = 0; x < Width; x += BlockSize) + { + HevcInverseTransformer.TransformAdd( + this.coefficients, + this.destination.AsSpan((y * Width) + x), + Width, + BlockLog2, + BlockLog2, + 12, + 18, + false, + this.scratch); + } + } + + return this.destination[^1]; + } +} diff --git a/tests/ImageSharp.Tests/Formats/Heif/Hevc/HevcInverseTransformerTests.cs b/tests/ImageSharp.Tests/Formats/Heif/Hevc/HevcInverseTransformerTests.cs new file mode 100644 index 0000000000..f7374e7598 --- /dev/null +++ b/tests/ImageSharp.Tests/Formats/Heif/Hevc/HevcInverseTransformerTests.cs @@ -0,0 +1,248 @@ +// Copyright (c) Six Labors. +// Licensed under the Six Labors Split License. + +using SixLabors.ImageSharp.Formats.Heif.Hevc; + +namespace SixLabors.ImageSharp.Tests.Formats.Heif.Hevc; + +/// +/// Verifies HEVC inverse-DCT, inverse-DST, rectangular-block, dynamic-range, and reconstruction behavior. +/// +[Trait("Format", "Heic")] +public class HevcInverseTransformerTests +{ + /// + /// Verifies the two normative normalization stages with a pure DC coefficient and a strided prediction block. + /// + [Fact] + public void DctDcCoefficientProducesUniformReconstruction() + { + const int size = 4; + const int stride = 7; + int[] coefficients = new int[size * size]; + coefficients[0] = 1024; + ushort[] destination = new ushort[stride * size]; + destination.AsSpan().Fill(100); + int[] scratch = new int[HevcInverseTransformer.GetScratchLength(2, 2)]; + + HevcInverseTransformer.TransformAdd(coefficients, destination, stride, 2, 2, 8, 15, false, scratch); + + for (int y = 0; y < size; y++) + { + for (int x = 0; x < size; x++) + { + Assert.Equal((ushort)108, destination[(y * stride) + x]); + } + } + } + + /// + /// Verifies a hand-derived four-by-four inverse-DST result from a pure lowest-frequency coefficient. + /// + [Fact] + public void DstLowestFrequencyMatchesFixedResult() + { + int[] coefficients = new int[16]; + coefficients[0] = 1024; + int[] actual = new int[16]; + int[] scratch = new int[HevcInverseTransformer.GetScratchLength(2, 2)]; + int[] expected = + [ + 2, 3, 4, 5, + 3, 6, 8, 9, + 4, 8, 11, 12, + 5, 9, 12, 14 + ]; + + HevcInverseTransformer.Transform(coefficients, actual, 2, 2, 8, 15, true, scratch); + + Assert.True(expected.AsSpan().SequenceEqual(actual)); + } + + /// + /// Compares factorized SIMD transforms with a dense specification-shaped oracle across supported dimensions and precisions. + /// + /// The base-two logarithm of the tested block width. + /// The base-two logarithm of the tested block height. + /// The reconstructed component precision. + /// The transform dynamic range excluding its sign bit. + /// Whether the tested four-by-four block uses the inverse DST. + [Theory] + [InlineData(2, 2, 8, 15, false)] + [InlineData(2, 2, 8, 15, true)] + [InlineData(2, 3, 10, 15, false)] + [InlineData(3, 2, 10, 15, false)] + [InlineData(3, 3, 10, 15, false)] + [InlineData(4, 3, 12, 18, false)] + [InlineData(3, 4, 12, 18, false)] + [InlineData(5, 4, 12, 18, false)] + [InlineData(4, 5, 12, 18, false)] + [InlineData(5, 5, 12, 15, false)] + public void TransformMatchesDenseOracle( + int log2Width, + int log2Height, + int bitDepth, + int maxTransformDynamicRange, + bool useDiscreteSineTransform) + { + int width = 1 << log2Width; + int height = 1 << log2Height; + int sampleCount = width * height; + int[] coefficients = new int[sampleCount]; + int coefficientMaximum = (1 << maxTransformDynamicRange) - 1; + for (int i = 0; i < coefficients.Length; i++) + { + // Alternating values across the complete dequantized range exercise both intermediate clipping bounds + // as well as every frequency group used by the factorized transform. + int unit = (((i * 37) + (width * 11) + height) % 127) - 63; + coefficients[i] = (unit * coefficientMaximum) / 63; + } + + int[] expected = new int[sampleCount]; + int[] actual = new int[sampleCount]; + int[] oracleScratch = new int[sampleCount]; + int[] transformScratch = new int[HevcInverseTransformer.GetScratchLength(log2Width, log2Height)]; + TransformDenseOracle( + coefficients, + expected, + oracleScratch, + width, + height, + bitDepth, + maxTransformDynamicRange, + useDiscreteSineTransform); + + HevcInverseTransformer.Transform( + coefficients, + actual, + log2Width, + log2Height, + bitDepth, + maxTransformDynamicRange, + useDiscreteSineTransform, + transformScratch); + + Assert.True(expected.AsSpan().SequenceEqual(actual), $"The {width}x{height} transform at {bitDepth} bits did not match the dense HEVC oracle."); + } + + /// + /// Verifies that residual addition clips both negative and positive reconstruction overflow at the component range. + /// + [Fact] + public void TransformAddClipsToComponentRange() + { + const int size = 4; + int[] positiveCoefficients = new int[size * size]; + int[] negativeCoefficients = new int[size * size]; + positiveCoefficients[0] = 32767; + negativeCoefficients[0] = -32768; + ushort[] positive = new ushort[size * size]; + ushort[] negative = new ushort[size * size]; + positive.AsSpan().Fill(4090); + negative.AsSpan().Fill(5); + int[] scratch = new int[HevcInverseTransformer.GetScratchLength(2, 2)]; + + HevcInverseTransformer.TransformAdd(positiveCoefficients, positive, size, 2, 2, 12, 18, false, scratch); + HevcInverseTransformer.TransformAdd(negativeCoefficients, negative, size, 2, 2, 12, 18, false, scratch); + + Assert.All(positive, value => Assert.Equal((ushort)4095, value)); + Assert.All(negative, value => Assert.Equal((ushort)0, value)); + } + + /// + /// Applies both inverse-transform dimensions using direct matrix products and normative rounding points. + /// + /// The dequantized coefficient block. + /// The destination signed residual block. + /// The full-block intermediate buffer. + /// The transform-block width. + /// The transform-block height. + /// The reconstructed component precision. + /// The transform dynamic range excluding its sign bit. + /// Whether both dimensions use the four-point inverse DST. + private static void TransformDenseOracle( + ReadOnlySpan coefficients, + Span residual, + Span intermediate, + int width, + int height, + int bitDepth, + int maxTransformDynamicRange, + bool useDiscreteSineTransform) + { + int dynamicMinimum = -(1 << maxTransformDynamicRange); + int dynamicMaximum = (1 << maxTransformDynamicRange) - 1; + for (int y = 0; y < height; y++) + { + for (int xFrequency = 0; xFrequency < width; xFrequency++) + { + int sum = 0; + for (int yFrequency = 0; yFrequency < height; yFrequency++) + { + sum += coefficients[(yFrequency * width) + xFrequency] * GetInverseCoefficient(height, yFrequency, y, useDiscreteSineTransform); + } + + intermediate[(y * width) + xFrequency] = Math.Clamp((sum + 64) >> 7, dynamicMinimum, dynamicMaximum); + } + } + + int secondShift = maxTransformDynamicRange + 5 - bitDepth; + int secondRounding = 1 << (secondShift - 1); + for (int y = 0; y < height; y++) + { + for (int x = 0; x < width; x++) + { + int sum = 0; + for (int xFrequency = 0; xFrequency < width; xFrequency++) + { + sum += intermediate[(y * width) + xFrequency] * GetInverseCoefficient(width, xFrequency, x, useDiscreteSineTransform); + } + + residual[(y * width) + x] = Math.Clamp((sum + secondRounding) >> secondShift, short.MinValue, short.MaxValue); + } + } + } + + /// + /// Gets one coefficient from the normative HEVC inverse-DCT or four-point inverse-DST matrix. + /// + /// The transform dimension. + /// The frequency-domain coordinate. + /// The spatial-domain coordinate. + /// Whether the four-point inverse DST is selected. + /// The signed matrix coefficient. + private static int GetInverseCoefficient(int size, int frequency, int position, bool useDiscreteSineTransform) + { + if (useDiscreteSineTransform) + { + ReadOnlySpan sine = + [ + 29, 55, 74, 84, + 74, 74, 0, -74, + 84, -29, -74, 55, + 55, -84, 74, -29 + ]; + + return sine[(frequency * 4) + position]; + } + + if (frequency == 0) + { + return 64; + } + + ReadOnlySpan magnitudes = + [ + 90, 90, 90, 90, 89, 88, 87, 85, 83, 82, 80, 78, 75, 73, 70, 67, 64, + 61, 57, 54, 50, 46, 43, 38, 36, 31, 25, 22, 18, 13, 9, 4, 0 + ]; + + int angle = (((2 * position) + 1) * frequency * (32 / size)) & 127; + if (angle > 64) + { + angle = 128 - angle; + } + + return angle > 32 ? -magnitudes[64 - angle] : magnitudes[angle]; + } +}