Browse Source

Implement SIMD-first HEVC inverse transforms

pull/2633/head
James Jackson-South 2 months ago
parent
commit
f512abf4a1
  1. 17
      HEIF_IMPLEMENTATION_PLAN.md
  2. 523
      src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.Operations.cs
  3. 141
      src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.OperatorImplementations.cs
  4. 270
      src/ImageSharp/Formats/Heif/Hevc/HevcInverseTransformer.cs
  5. 11
      src/ImageSharp/Formats/Heif/Hevc/HevcSequenceParameterSet.cs
  6. 90
      tests/ImageSharp.Benchmarks/Codecs/Heif/HevcInverseTransformBenchmarks.cs
  7. 248
      tests/ImageSharp.Tests/Formats/Heif/Hevc/HevcInverseTransformerTests.cs

17
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.

523
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
{
/// <summary>
/// Calculates the disjoint odd-frequency groups that seed the HEVC partial-butterfly reconstruction.
/// </summary>
/// <typeparam name="TOperator">The selected inverse-DCT operator.</typeparam>
/// <param name="source">The frequency rows followed by contiguous independent lines.</param>
/// <param name="groups">The destination group rows.</param>
/// <param name="lineCount">The number of independent lines transformed together.</param>
private static void PopulateButterflyGroups<TOperator>(ReadOnlySpan<int> source, Span<int> groups, int lineCount)
where TOperator : struct, IHevcInverseTransformOperator<TOperator>
{
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<TOperator>(
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<TOperator>(
source,
groups.Slice((groupOffset + position) * lineCount, lineCount),
lineCount,
0,
size >> 1,
position);
}
}
/// <summary>
/// Calculates one partial-butterfly group row across all independent lines.
/// </summary>
/// <typeparam name="TOperator">The selected inverse-DCT operator.</typeparam>
/// <param name="source">The complete frequency-row input.</param>
/// <param name="destination">The destination group row.</param>
/// <param name="lineCount">The number of independent lines.</param>
/// <param name="firstFrequency">The first frequency included by the group.</param>
/// <param name="frequencyStep">The distance between included frequencies.</param>
/// <param name="position">The group-relative spatial coordinate.</param>
private static void PopulateButterflyGroupRow<TOperator>(
ReadOnlySpan<int> source,
Span<int> destination,
int lineCount,
int firstFrequency,
int frequencyStep,
int position)
where TOperator : struct, IHevcInverseTransformOperator<TOperator>
{
ref int sourceBase = ref MemoryMarshal.GetReference(source);
ref int destinationBase = ref MemoryMarshal.GetReference(destination);
int x = 0;
if (Vector512.IsHardwareAccelerated)
{
int oneVectorFromEnd = lineCount - Vector512<int>.Count;
for (; x <= oneVectorFromEnd; x += Vector512<int>.Count)
{
Vector512<int> sum = Vector512<int>.Zero;
for (int frequency = firstFrequency; frequency < TOperator.Size; frequency += frequencyStep)
{
Vector512<int> 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<int>.Count;
for (; x <= oneVectorFromEnd; x += Vector256<int>.Count)
{
Vector256<int> sum = Vector256<int>.Zero;
for (int frequency = firstFrequency; frequency < TOperator.Size; frequency += frequencyStep)
{
Vector256<int> 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<int>.Count;
for (; x <= oneVectorFromEnd; x += Vector128<int>.Count)
{
Vector128<int> sum = Vector128<int>.Zero;
for (int frequency = firstFrequency; frequency < TOperator.Size; frequency += frequencyStep)
{
Vector128<int> 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;
}
}
/// <summary>
/// Expands the disjoint partial-butterfly groups into spatial rows.
/// </summary>
/// <typeparam name="TOperator">The selected inverse-DCT operator.</typeparam>
/// <param name="initial">The buffer containing every disjoint group.</param>
/// <param name="alternate">The alternate expansion buffer.</param>
/// <param name="lineCount">The number of independent lines transformed together.</param>
/// <param name="shift">The rounded right shift applied at the final hierarchy level.</param>
/// <param name="minimum">The inclusive output minimum.</param>
/// <param name="maximum">The inclusive output maximum.</param>
/// <returns>The buffer containing the completed spatial rows.</returns>
private static Span<int> CombineButterflyGroups<TOperator>(
Span<int> initial,
Span<int> alternate,
int lineCount,
int shift,
int minimum,
int maximum)
where TOperator : struct, IHevcInverseTransformOperator<TOperator>
{
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<int> even = currentIsInitial ? initial : alternate;
Span<int> destination = currentIsInitial ? alternate : initial;
for (int position = 0; position < combinedSize; position++)
{
ReadOnlySpan<int> evenRow = even.Slice((combinedStart + position) * lineCount, lineCount);
ReadOnlySpan<int> oddRow = initial.Slice((oddStart + position) * lineCount, lineCount);
Span<int> positiveRow = destination.Slice((oddStart + position) * lineCount, lineCount);
Span<int> 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;
}
/// <summary>
/// Combines one symmetric pair of partial-butterfly rows.
/// </summary>
/// <param name="even">The even-frequency contribution.</param>
/// <param name="odd">The odd-frequency contribution.</param>
/// <param name="positive">The destination receiving the added contribution.</param>
/// <param name="negative">The destination receiving the subtracted contribution.</param>
/// <param name="roundAndClip">Whether this is the final hierarchy level.</param>
/// <param name="shift">The rounded right shift applied at the final hierarchy level.</param>
/// <param name="minimum">The inclusive output minimum.</param>
/// <param name="maximum">The inclusive output maximum.</param>
private static void CombineButterflyRows(
ReadOnlySpan<int> even,
ReadOnlySpan<int> odd,
Span<int> positive,
Span<int> 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<int>.Count;
for (; x <= oneVectorFromEnd; x += Vector512<int>.Count)
{
Vector512<int> evenValues = Vector512.LoadUnsafe(ref evenBase, (nuint)x);
Vector512<int> oddValues = Vector512.LoadUnsafe(ref oddBase, (nuint)x);
Vector512<int> added = evenValues + oddValues;
Vector512<int> 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<int>.Count;
for (; x <= oneVectorFromEnd; x += Vector256<int>.Count)
{
Vector256<int> evenValues = Vector256.LoadUnsafe(ref evenBase, (nuint)x);
Vector256<int> oddValues = Vector256.LoadUnsafe(ref oddBase, (nuint)x);
Vector256<int> added = evenValues + oddValues;
Vector256<int> 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<int>.Count;
for (; x <= oneVectorFromEnd; x += Vector128<int>.Count)
{
Vector128<int> evenValues = Vector128.LoadUnsafe(ref evenBase, (nuint)x);
Vector128<int> oddValues = Vector128.LoadUnsafe(ref oddBase, (nuint)x);
Vector128<int> added = evenValues + oddValues;
Vector128<int> 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;
}
}
/// <summary>
/// Applies the four-point inverse-DST matrix across all independent lines.
/// </summary>
/// <typeparam name="TOperator">The selected dense transform operator.</typeparam>
/// <param name="source">The frequency rows followed by contiguous independent lines.</param>
/// <param name="destination">The spatial rows followed by contiguous independent lines.</param>
/// <param name="lineCount">The number of independent lines transformed together.</param>
/// <param name="shift">The rounded right shift applied to each result.</param>
/// <param name="minimum">The inclusive output minimum.</param>
/// <param name="maximum">The inclusive output maximum.</param>
private static void TransformDense<TOperator>(
ReadOnlySpan<int> source,
Span<int> destination,
int lineCount,
int shift,
int minimum,
int maximum)
where TOperator : struct, IHevcInverseTransformOperator<TOperator>
{
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<int>.Count;
for (; x <= oneVectorFromEnd; x += Vector512<int>.Count)
{
Vector512<int> sum = Vector512<int>.Zero;
for (int frequency = 0; frequency < TOperator.Size; frequency++)
{
Vector512<int> 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<int>.Count;
for (; x <= oneVectorFromEnd; x += Vector256<int>.Count)
{
Vector256<int> sum = Vector256<int>.Zero;
for (int frequency = 0; frequency < TOperator.Size; frequency++)
{
Vector256<int> 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<int>.Count;
for (; x <= oneVectorFromEnd; x += Vector128<int>.Count)
{
Vector128<int> sum = Vector128<int>.Zero;
for (int frequency = 0; frequency < TOperator.Size; frequency++)
{
Vector128<int> 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);
}
}
}
/// <summary>
/// Transposes one rectangular transform block into a separate full-block buffer.
/// </summary>
/// <param name="source">The source block in raster order.</param>
/// <param name="destination">The transposed destination block.</param>
/// <param name="sourceHeight">The source row count.</param>
/// <param name="sourceWidth">The source column count.</param>
private static void Transpose(ReadOnlySpan<int> source, Span<int> 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<int> row0 = Vector128.LoadUnsafe(ref sourceBase, (nuint)((y * sourceWidth) + x));
Vector128<int> row1 = Vector128.LoadUnsafe(ref sourceBase, (nuint)(((y + 1) * sourceWidth) + x));
Vector128<int> row2 = Vector128.LoadUnsafe(ref sourceBase, (nuint)(((y + 2) * sourceWidth) + x));
Vector128<int> 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];
}
}
}
/// <summary>
/// Adds a complete signed residual block to the predicted samples and clips to the component precision.
/// </summary>
/// <param name="residual">The signed residual block in raster order.</param>
/// <param name="destination">The predicted samples beginning at the block origin.</param>
/// <param name="destinationStride">The destination row stride in samples.</param>
/// <param name="width">The transform-block width.</param>
/// <param name="height">The transform-block height.</param>
/// <param name="bitDepth">The reconstructed component precision.</param>
private static void AddResidual(ReadOnlySpan<int> residual, Span<ushort> destination, int destinationStride, int width, int height, int bitDepth)
{
int maximum = (1 << bitDepth) - 1;
for (int y = 0; y < height; y++)
{
ReadOnlySpan<int> residualRow = residual.Slice(y * width, width);
Span<ushort> 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<int>.Count;
for (; x <= oneVectorFromEnd; x += Vector256<int>.Count)
{
Vector128<ushort> predicted16 = Vector128.LoadUnsafe(ref destinationBase, (nuint)x);
Vector256<int> predicted = Vector256.Create(Vector128.WidenLower(predicted16), Vector128.WidenUpper(predicted16)).AsInt32();
Vector256<int> reconstructed = Vector256.Clamp(
predicted + Vector256.LoadUnsafe(ref residualBase, (nuint)x),
Vector256<int>.Zero,
Vector256.Create(maximum));
Vector128.Narrow(reconstructed.GetLower().AsUInt32(), reconstructed.GetUpper().AsUInt32()).StoreUnsafe(ref destinationBase, (nuint)x);
}
}
if (Vector128.IsHardwareAccelerated)
{
int oneVectorFromEnd = width - Vector128<int>.Count;
for (; x <= oneVectorFromEnd; x += Vector128<int>.Count)
{
ulong packed = Unsafe.ReadUnaligned<ulong>(ref Unsafe.As<ushort, byte>(ref Unsafe.Add(ref destinationBase, x)));
Vector128<int> predicted = Vector128.WidenLower(Vector128.CreateScalarUnsafe(packed).AsUInt16()).AsInt32();
Vector128<int> reconstructed = Vector128.Clamp(
predicted + Vector128.LoadUnsafe(ref residualBase, (nuint)x),
Vector128<int>.Zero,
Vector128.Create(maximum));
Vector128<ushort> narrowed = Vector128.Narrow(reconstructed.AsUInt32(), Vector128<uint>.Zero);
Unsafe.WriteUnaligned(ref Unsafe.As<ushort, byte>(ref Unsafe.Add(ref destinationBase, x)), narrowed.AsUInt64().ToScalar());
}
}
for (; x < width; x++)
{
destinationRow[x] = (ushort)Math.Clamp(destinationRow[x] + residualRow[x], 0, maximum);
}
}
}
/// <summary>
/// Applies HEVC's rounded right shift and inclusive clipping to a 512-bit vector.
/// </summary>
/// <param name="value">The unnormalized transform values.</param>
/// <param name="shift">The right-shift count.</param>
/// <param name="minimum">The inclusive result minimum.</param>
/// <param name="maximum">The inclusive result maximum.</param>
/// <returns>The normalized and clipped values.</returns>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static Vector512<int> RoundShiftAndClamp(Vector512<int> value, int shift, int minimum, int maximum)
=> Vector512.Clamp((value + Vector512.Create(1 << (shift - 1))) >> shift, Vector512.Create(minimum), Vector512.Create(maximum));
/// <summary>
/// Applies HEVC's rounded right shift and inclusive clipping to a 256-bit vector.
/// </summary>
/// <param name="value">The unnormalized transform values.</param>
/// <param name="shift">The right-shift count.</param>
/// <param name="minimum">The inclusive result minimum.</param>
/// <param name="maximum">The inclusive result maximum.</param>
/// <returns>The normalized and clipped values.</returns>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static Vector256<int> RoundShiftAndClamp(Vector256<int> value, int shift, int minimum, int maximum)
=> Vector256.Clamp((value + Vector256.Create(1 << (shift - 1))) >> shift, Vector256.Create(minimum), Vector256.Create(maximum));
/// <summary>
/// Applies HEVC's rounded right shift and inclusive clipping to a 128-bit vector.
/// </summary>
/// <param name="value">The unnormalized transform values.</param>
/// <param name="shift">The right-shift count.</param>
/// <param name="minimum">The inclusive result minimum.</param>
/// <param name="maximum">The inclusive result maximum.</param>
/// <returns>The normalized and clipped values.</returns>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static Vector128<int> RoundShiftAndClamp(Vector128<int> value, int shift, int minimum, int maximum)
=> Vector128.Clamp((value + Vector128.Create(1 << (shift - 1))) >> shift, Vector128.Create(minimum), Vector128.Create(maximum));
/// <summary>
/// Applies HEVC's rounded right shift and inclusive clipping to one scalar value.
/// </summary>
/// <param name="value">The unnormalized transform value.</param>
/// <param name="shift">The right-shift count.</param>
/// <param name="minimum">The inclusive result minimum.</param>
/// <param name="maximum">The inclusive result maximum.</param>
/// <returns>The normalized and clipped value.</returns>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static int RoundShiftAndClamp(int value, int shift, int minimum, int maximum)
=> Math.Clamp((value + (1 << (shift - 1))) >> shift, minimum, maximum);
}

141
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
{
/// <summary>
/// Gets the common HEVC inverse-DCT magnitudes ordered on the pi-over-sixty-four angle grid.
/// </summary>
private static ReadOnlySpan<sbyte> 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
];
/// <summary>
/// Gets the four-point HEVC inverse-DST matrix in frequency-major order.
/// </summary>
private static ReadOnlySpan<sbyte> DiscreteSine4Coefficients =>
[
29, 55, 74, 84,
74, 74, 0, -74,
84, -29, -74, 55,
55, -84, 74, -29
];
/// <summary>
/// Reconstructs one normative inverse-DCT coefficient from the common HEVC angle grid.
/// </summary>
/// <param name="size">The transform side.</param>
/// <param name="frequency">The frequency-domain coordinate.</param>
/// <param name="position">The spatial-domain coordinate.</param>
/// <returns>The signed transform coefficient.</returns>
[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];
}
/// <summary>
/// Implements the four-point inverse discrete cosine transform.
/// </summary>
private readonly struct DiscreteCosine4Operator : IHevcInverseTransformOperator<DiscreteCosine4Operator>
{
/// <inheritdoc/>
public static int Size => 4;
/// <inheritdoc/>
public static bool UsesButterfly => true;
/// <inheritdoc/>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static int GetCoefficient(int frequency, int position) => GetDiscreteCosineCoefficient(Size, frequency, position);
}
/// <summary>
/// Implements the four-point inverse discrete sine transform.
/// </summary>
private readonly struct DiscreteSine4Operator : IHevcInverseTransformOperator<DiscreteSine4Operator>
{
/// <inheritdoc/>
public static int Size => 4;
/// <inheritdoc/>
public static bool UsesButterfly => false;
/// <inheritdoc/>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static int GetCoefficient(int frequency, int position) => DiscreteSine4Coefficients[(frequency * Size) + position];
}
/// <summary>
/// Implements the eight-point inverse discrete cosine transform.
/// </summary>
private readonly struct DiscreteCosine8Operator : IHevcInverseTransformOperator<DiscreteCosine8Operator>
{
/// <inheritdoc/>
public static int Size => 8;
/// <inheritdoc/>
public static bool UsesButterfly => true;
/// <inheritdoc/>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static int GetCoefficient(int frequency, int position) => GetDiscreteCosineCoefficient(Size, frequency, position);
}
/// <summary>
/// Implements the sixteen-point inverse discrete cosine transform.
/// </summary>
private readonly struct DiscreteCosine16Operator : IHevcInverseTransformOperator<DiscreteCosine16Operator>
{
/// <inheritdoc/>
public static int Size => 16;
/// <inheritdoc/>
public static bool UsesButterfly => true;
/// <inheritdoc/>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static int GetCoefficient(int frequency, int position) => GetDiscreteCosineCoefficient(Size, frequency, position);
}
/// <summary>
/// Implements the thirty-two-point inverse discrete cosine transform.
/// </summary>
private readonly struct DiscreteCosine32Operator : IHevcInverseTransformOperator<DiscreteCosine32Operator>
{
/// <inheritdoc/>
public static int Size => 32;
/// <inheritdoc/>
public static bool UsesButterfly => true;
/// <inheritdoc/>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static int GetCoefficient(int frequency, int position) => GetDiscreteCosineCoefficient(Size, frequency, position);
}
}

270
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;
/// <summary>
/// Applies HEVC inverse transforms and reconstructs predicted samples.
/// </summary>
internal static partial class HevcInverseTransformer
{
/// <summary>
/// The signed residual precision used after the second inverse-transform pass.
/// </summary>
private const int ResidualPrecision = 16;
/// <summary>
/// Defines one closed inverse-transform operation selected by the transform-unit syntax.
/// </summary>
/// <typeparam name="TOperator">The implementing operator type.</typeparam>
private interface IHevcInverseTransformOperator<TOperator>
where TOperator : struct, IHevcInverseTransformOperator<TOperator>
{
/// <summary>
/// Gets the transform side in samples.
/// </summary>
public static abstract int Size { get; }
/// <summary>
/// Gets a value indicating whether the transform uses the partial-butterfly factorization.
/// </summary>
public static abstract bool UsesButterfly { get; }
/// <summary>
/// Gets one inverse-transform matrix coefficient.
/// </summary>
/// <param name="frequency">The frequency-domain coordinate.</param>
/// <param name="position">The spatial-domain coordinate.</param>
/// <returns>The signed transform coefficient.</returns>
public static abstract int GetCoefficient(int frequency, int position);
}
/// <summary>
/// Gets the scratch length required for the specified rectangular transform block.
/// </summary>
/// <param name="log2Width">The base-two logarithm of the transform-block width.</param>
/// <param name="log2Height">The base-two logarithm of the transform-block height.</param>
/// <returns>The required number of signed thirty-two-bit elements.</returns>
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);
}
/// <summary>
/// Reconstructs one transform block by adding its inverse-transformed residual to the predicted samples.
/// </summary>
/// <param name="coefficients">The dequantized transform coefficients in raster order.</param>
/// <param name="destination">The predicted samples beginning at the transform-block origin.</param>
/// <param name="destinationStride">The destination row stride in samples.</param>
/// <param name="log2Width">The base-two logarithm of the transform-block width.</param>
/// <param name="log2Height">The base-two logarithm of the transform-block height.</param>
/// <param name="bitDepth">The reconstructed component precision.</param>
/// <param name="maxTransformDynamicRange">The transform dynamic range excluding its sign bit.</param>
/// <param name="useDiscreteSineTransform">Whether the four-by-four luma intra block uses the discrete sine transform.</param>
/// <param name="scratch">The caller-owned scratch returned by <see cref="GetScratchLength(int, int)"/>.</param>
public static void TransformAdd(
ReadOnlySpan<int> coefficients,
Span<ushort> destination,
int destinationStride,
int log2Width,
int log2Height,
int bitDepth,
int maxTransformDynamicRange,
bool useDiscreteSineTransform,
Span<int> 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<int> first = scratch[..sampleCount];
Span<int> second = scratch.Slice(sampleCount, sampleCount);
Span<int> residual = TransformCore(
coefficients[..sampleCount],
first,
second,
width,
height,
bitDepth,
maxTransformDynamicRange,
useDiscreteSineTransform);
AddResidual(residual, destination, destinationStride, width, height, bitDepth);
}
/// <summary>
/// Applies one two-dimensional inverse transform and writes signed residual samples in raster order.
/// </summary>
/// <param name="coefficients">The dequantized transform coefficients in raster order.</param>
/// <param name="residual">The destination residual samples in raster order.</param>
/// <param name="log2Width">The base-two logarithm of the transform-block width.</param>
/// <param name="log2Height">The base-two logarithm of the transform-block height.</param>
/// <param name="bitDepth">The reconstructed component precision.</param>
/// <param name="maxTransformDynamicRange">The transform dynamic range excluding its sign bit.</param>
/// <param name="useDiscreteSineTransform">Whether the four-by-four luma intra block uses the discrete sine transform.</param>
/// <param name="scratch">The caller-owned scratch returned by <see cref="GetScratchLength(int, int)"/>.</param>
public static void Transform(
ReadOnlySpan<int> coefficients,
Span<int> residual,
int log2Width,
int log2Height,
int bitDepth,
int maxTransformDynamicRange,
bool useDiscreteSineTransform,
Span<int> 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<int> transformed = TransformCore(
coefficients[..sampleCount],
scratch[..sampleCount],
scratch.Slice(sampleCount, sampleCount),
width,
height,
bitDepth,
maxTransformDynamicRange,
useDiscreteSineTransform);
transformed.CopyTo(residual);
}
/// <summary>
/// Dispatches both separable transform passes through closed operators selected from the block dimensions.
/// </summary>
/// <param name="coefficients">The complete dequantized coefficient block.</param>
/// <param name="first">The first full-block scratch buffer.</param>
/// <param name="second">The second full-block scratch buffer.</param>
/// <param name="width">The transform-block width.</param>
/// <param name="height">The transform-block height.</param>
/// <param name="bitDepth">The reconstructed component precision.</param>
/// <param name="maxTransformDynamicRange">The transform dynamic range excluding its sign bit.</param>
/// <param name="useDiscreteSineTransform">Whether both four-point passes use the discrete sine transform.</param>
/// <returns>The scratch buffer containing the raster-ordered residual.</returns>
private static Span<int> TransformCore(
ReadOnlySpan<int> coefficients,
Span<int> first,
Span<int> 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<int> vertical = TransformDimension(
coefficients,
first,
second,
height,
width,
7,
dynamicMinimum,
dynamicMaximum,
useDiscreteSineTransform);
Span<int> horizontalInput = vertical.Overlaps(first) ? second : first;
Transpose(vertical, horizontalInput, height, width);
Span<int> horizontalWorkspace = horizontalInput.Overlaps(first) ? second : first;
int secondShift = maxTransformDynamicRange + 5 - bitDepth;
Span<int> horizontal = TransformDimension(
horizontalInput,
horizontalWorkspace,
horizontalInput,
width,
height,
secondShift,
-(1 << (ResidualPrecision - 1)),
(1 << (ResidualPrecision - 1)) - 1,
useDiscreteSineTransform);
Span<int> residual = horizontal.Overlaps(first) ? second : first;
Transpose(horizontal, residual, width, height);
return residual;
}
/// <summary>
/// Selects the statically specialized operator for one transform dimension.
/// </summary>
/// <param name="source">The frequency rows followed by contiguous independent lines.</param>
/// <param name="initial">The initial operator output buffer.</param>
/// <param name="alternate">The alternate combination buffer.</param>
/// <param name="size">The transform dimension.</param>
/// <param name="lineCount">The number of independent lines transformed together.</param>
/// <param name="shift">The rounded right shift applied to the spatial results.</param>
/// <param name="minimum">The inclusive output minimum.</param>
/// <param name="maximum">The inclusive output maximum.</param>
/// <param name="useDiscreteSineTransform">Whether the four-point pass uses the discrete sine transform.</param>
/// <returns>The buffer containing spatial rows followed by contiguous independent lines.</returns>
private static Span<int> TransformDimension(
ReadOnlySpan<int> source,
Span<int> initial,
Span<int> alternate,
int size,
int lineCount,
int shift,
int minimum,
int maximum,
bool useDiscreteSineTransform)
=> (size, useDiscreteSineTransform) switch
{
(4, true) => TransformDimension<DiscreteSine4Operator>(source, initial, alternate, lineCount, shift, minimum, maximum),
(4, false) => TransformDimension<DiscreteCosine4Operator>(source, initial, alternate, lineCount, shift, minimum, maximum),
(8, _) => TransformDimension<DiscreteCosine8Operator>(source, initial, alternate, lineCount, shift, minimum, maximum),
(16, _) => TransformDimension<DiscreteCosine16Operator>(source, initial, alternate, lineCount, shift, minimum, maximum),
_ => TransformDimension<DiscreteCosine32Operator>(source, initial, alternate, lineCount, shift, minimum, maximum)
};
/// <summary>
/// Invokes one statically selected inverse-transform operator.
/// </summary>
/// <typeparam name="TOperator">The selected inverse-transform operator.</typeparam>
/// <param name="source">The frequency rows followed by contiguous independent lines.</param>
/// <param name="initial">The initial operator output buffer.</param>
/// <param name="alternate">The alternate combination buffer.</param>
/// <param name="lineCount">The number of independent lines transformed together.</param>
/// <param name="shift">The rounded right shift applied to the spatial results.</param>
/// <param name="minimum">The inclusive output minimum.</param>
/// <param name="maximum">The inclusive output maximum.</param>
/// <returns>The buffer containing spatial rows followed by contiguous independent lines.</returns>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static Span<int> TransformDimension<TOperator>(
ReadOnlySpan<int> source,
Span<int> initial,
Span<int> alternate,
int lineCount,
int shift,
int minimum,
int maximum)
where TOperator : struct, IHevcInverseTransformOperator<TOperator>
{
if (!TOperator.UsesButterfly)
{
TransformDense<TOperator>(source, initial, lineCount, shift, minimum, maximum);
return initial;
}
PopulateButterflyGroups<TOperator>(source, initial, lineCount);
return CombineButterflyGroups<TOperator>(initial, alternate, lineCount, shift, minimum, maximum);
}
}

11
src/ImageSharp/Formats/Heif/Hevc/HevcSequenceParameterSet.cs

@ -525,6 +525,17 @@ internal sealed class HevcSequenceParameterSet
/// </summary>
public bool CabacBypassAlignmentEnabled { get; }
/// <summary>
/// Gets the base-two logarithm of the transform dynamic range for the specified reconstructed plane.
/// </summary>
/// <param name="plane">The reconstructed plane.</param>
/// <returns>The transform dynamic range excluding its sign bit.</returns>
public int GetMaxTransformDynamicRange(HevcPlane plane)
{
int bitDepth = plane == HevcPlane.Y ? this.BitDepthLuma : this.BitDepthChroma;
return this.ExtendedPrecisionProcessingEnabled ? Math.Max(15, bitDepth + 6) : 15;
}
/// <summary>
/// Reads a conformance-window offset and converts it to luma-sample units.
/// </summary>

90
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;
/// <summary>
/// Measures complete coded-frame traversal through HEVC inverse-transform reconstruction.
/// </summary>
[MemoryDiagnoser(displayGenColumns: false)]
public class HevcInverseTransformBenchmarks
{
/// <summary>
/// The coded frame width, which is an exact multiple of the maximum transform-block side.
/// </summary>
private const int Width = 1920;
/// <summary>
/// The coded frame height including the final padded coding-tree row for a 1080-line presentation.
/// </summary>
private const int Height = 1088;
/// <summary>
/// The base-two logarithm of the benchmark transform-block side.
/// </summary>
private const int BlockLog2 = 5;
/// <summary>
/// The transform-block side in samples.
/// </summary>
private const int BlockSize = 1 << BlockLog2;
/// <summary>
/// The deterministic dequantized coefficients reused by each benchmark block.
/// </summary>
private readonly int[] coefficients = new int[BlockSize * BlockSize];
/// <summary>
/// The frame-wide predicted and reconstructed samples.
/// </summary>
private readonly ushort[] destination = new ushort[Width * Height];
/// <summary>
/// The maximum-block inverse-transform scratch reused throughout each coded frame.
/// </summary>
private readonly int[] scratch = new int[HevcInverseTransformer.GetScratchLength(BlockLog2, BlockLog2)];
/// <summary>
/// Populates a dense, deterministic twelve-bit transform workload outside the measured frame traversal.
/// </summary>
[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);
}
/// <summary>
/// Measures factorized inverse DCT, both transpositions, and saturated prediction addition for a coded frame.
/// </summary>
/// <returns>The final reconstructed sample, keeping the frame output observable.</returns>
[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];
}
}

248
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;
/// <summary>
/// Verifies HEVC inverse-DCT, inverse-DST, rectangular-block, dynamic-range, and reconstruction behavior.
/// </summary>
[Trait("Format", "Heic")]
public class HevcInverseTransformerTests
{
/// <summary>
/// Verifies the two normative normalization stages with a pure DC coefficient and a strided prediction block.
/// </summary>
[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]);
}
}
}
/// <summary>
/// Verifies a hand-derived four-by-four inverse-DST result from a pure lowest-frequency coefficient.
/// </summary>
[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));
}
/// <summary>
/// Compares factorized SIMD transforms with a dense specification-shaped oracle across supported dimensions and precisions.
/// </summary>
/// <param name="log2Width">The base-two logarithm of the tested block width.</param>
/// <param name="log2Height">The base-two logarithm of the tested block height.</param>
/// <param name="bitDepth">The reconstructed component precision.</param>
/// <param name="maxTransformDynamicRange">The transform dynamic range excluding its sign bit.</param>
/// <param name="useDiscreteSineTransform">Whether the tested four-by-four block uses the inverse DST.</param>
[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.");
}
/// <summary>
/// Verifies that residual addition clips both negative and positive reconstruction overflow at the component range.
/// </summary>
[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));
}
/// <summary>
/// Applies both inverse-transform dimensions using direct matrix products and normative rounding points.
/// </summary>
/// <param name="coefficients">The dequantized coefficient block.</param>
/// <param name="residual">The destination signed residual block.</param>
/// <param name="intermediate">The full-block intermediate buffer.</param>
/// <param name="width">The transform-block width.</param>
/// <param name="height">The transform-block height.</param>
/// <param name="bitDepth">The reconstructed component precision.</param>
/// <param name="maxTransformDynamicRange">The transform dynamic range excluding its sign bit.</param>
/// <param name="useDiscreteSineTransform">Whether both dimensions use the four-point inverse DST.</param>
private static void TransformDenseOracle(
ReadOnlySpan<int> coefficients,
Span<int> residual,
Span<int> 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);
}
}
}
/// <summary>
/// Gets one coefficient from the normative HEVC inverse-DCT or four-point inverse-DST matrix.
/// </summary>
/// <param name="size">The transform dimension.</param>
/// <param name="frequency">The frequency-domain coordinate.</param>
/// <param name="position">The spatial-domain coordinate.</param>
/// <param name="useDiscreteSineTransform">Whether the four-point inverse DST is selected.</param>
/// <returns>The signed matrix coefficient.</returns>
private static int GetInverseCoefficient(int size, int frequency, int position, bool useDiscreteSineTransform)
{
if (useDiscreteSineTransform)
{
ReadOnlySpan<sbyte> 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<sbyte> 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];
}
}
Loading…
Cancel
Save