From a4e4db6bc52075d5c878248436be1037bb9d16c9 Mon Sep 17 00:00:00 2001 From: James Jackson-South Date: Tue, 25 Aug 2026 02:21:07 +1000 Subject: [PATCH] Implement AV1 Wiener restoration kernel --- HEIF_IMPLEMENTATION_PLAN.md | 1 + .../LoopRestoration/Av1WienerFilter.cs | 193 ++++++++++++++++++ 2 files changed, 194 insertions(+) create mode 100644 src/ImageSharp/Formats/Heif/Av1/Pipeline/LoopRestoration/Av1WienerFilter.cs diff --git a/HEIF_IMPLEMENTATION_PLAN.md b/HEIF_IMPLEMENTATION_PLAN.md index 915512787..8f160e72f 100644 --- a/HEIF_IMPLEMENTATION_PLAN.md +++ b/HEIF_IMPLEMENTATION_PLAN.md @@ -75,6 +75,7 @@ This snapshot pins or classifies the available references and failures; it does | `Av1CdefDecoder`, `Av1CdefKernels`, and CDEF-unit strength storage | AV1 sections 7.15.2 through 7.15.4 constrained directional enhancement filtering | libaom `av1/common/cdef.c`, `av1/common/cdef_block.c`, `av1/common/cdef.h`, and `av1/common/cdef_block.h` at `03087864cf4bea6abb0d28f95cf7843511413d8f` | Port the scalar direction search, variance adjustment, constrained primary/secondary taps, subsampling direction conversion, skipped-8x8 selection, and frame-edge sentinel behavior. Use a frame-owned source snapshot so filtering never consumes already modified samples. This is normative AV1 still-image reconstruction and introduces no ISO BMFF, track, timing, or sequence-playback surface. | | `Av1SuperResolutionDecoder`, `Av1SuperResolutionKernels`, frame-size derivation, and decoded-image dimensions | AV1 section 7.16 normative super-resolution upscaling | libaom `av1/common/resize.c`, `av1/common/resize.h`, `av1/common/convolve.c`, and `aom_dsp/aom_filter.h` at `03087864cf4bea6abb0d28f95cf7843511413d8f` | Port the fixed 64-phase, 8-tap horizontal filter, phase/step derivation, replicated frame edges, chroma width rounding, signed rounding, and 8/10/12-bit clipping. Reuse ImageSharp's existing cross-platform `Vector128_.MultiplyAddAdjacent` helper for the exact eight-coefficient dot product with a scalar fallback. Generic image resizing is not normative AV1 super-resolution. This adds no track, timing, fragment, animation, or generic ISO BMFF model. | | `Av1TileReader` loop-restoration unit syntax, `Av1SymbolDecoder` restoration distributions/subexponential codes, and `Av1FrameInfo` unit storage | AV1 section 5.11.57 `read_lr` and `read_lr_unit` syntax | libaom `av1/decoder/decodeframe.c`, `av1/common/restoration.c`, `av1/common/restoration.h`, `av1/common/entropymode.c`, `aom_dsp/binary_codes_reader.c`, and `aom_dsp/recenter.h` at `03087864cf4bea6abb0d28f95cf7843511413d8f` | Decode tile-local switchable/Wiener/self-guided selections, finite reference-subexponential coefficients, chroma Wiener windows, self-guided parameter sets, super-resolution-adjusted unit corners, and the AV1 nearest-unit-count rule into frame-owned per-plane grids. This is compressed still-image syntax and adds no movie, track, timing, fragment, audio, or sequence surface. | +| `Av1WienerFilter` | AV1 section 7.17.3 Wiener restoration filtering | libaom `av1/common/restoration.c`, `av1/common/restoration.h`, `av1/common/convolve.c`, and `av1/common/convolve.h` at `03087864cf4bea6abb0d28f95cf7843511413d8f` | Preserve the implicit center-sample contribution, separable horizontal/vertical rounding, bit-depth-dependent 16-bit intermediate range, and final 8/10/12-bit clipping. Reuse `Vector128_.MultiplyAddAdjacent` for the contiguous horizontal eight-tap product with an exact scalar fallback. Keep this scalar/SIMD oracle disabled until restoration stripe boundaries and self-guided filtering are both complete. | | `Av1FrameInfo`, `Av1TileReader`, and `Av1BlockDecoder` transform/coefficient storage | AV1 section 5.11.39 coefficient syntax and section 7.11.2 reconstruction | libaom `av1/decoder/decodetxb.c` and `av1/decoder/decoder.h` at `03087864cf4bea6abb0d28f95cf7843511413d8f` | Preserve separate luma and chroma transform coefficients at monotonically advancing per-plane offsets within each superblock so reconstruction consumes the same transform-block order produced by tile parsing. | | `Av1InverseQuantizer` and `Av1InverseQuantizationLookup` | AV1 section 7.12.3 inverse quantization | libaom `aom_dsp/aom_dsp_common.h`, `av1/common/quant_common.c`, and `av1/decoder/decodetxb.c` at `03087864cf4bea6abb0d28f95cf7843511413d8f` | Select the per-segment matrix level, alias 64-pixel transform dimensions to their adjusted matrices, retain a flat level-15 matrix, and apply the five-bit inverse-matrix weight scale. The large managed lookup remains a single process-wide table. | | `Av1Inverse2dTransformer` and `Av1InverseTransformerFactory` | AV1 section 7.11.2 inverse transform and reconstruction | libaom `av1/common/av1_inv_txfm1d.c`, `av1/common/av1_inv_txfm2d.c`, and `av1/common/idct.c` at `03087864cf4bea6abb0d28f95cf7843511413d8f` | Scalar transform oracle for coefficient-row traversal, intermediate layout, stage ranges, clipping, and high-bit-depth sample addition. The managed 16-bit overload is also used as a parity oracle for the byte overload. | diff --git a/src/ImageSharp/Formats/Heif/Av1/Pipeline/LoopRestoration/Av1WienerFilter.cs b/src/ImageSharp/Formats/Heif/Av1/Pipeline/LoopRestoration/Av1WienerFilter.cs new file mode 100644 index 000000000..cb51e8e8f --- /dev/null +++ b/src/ImageSharp/Formats/Heif/Av1/Pipeline/LoopRestoration/Av1WienerFilter.cs @@ -0,0 +1,193 @@ +// 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.Common.Helpers; + +namespace SixLabors.ImageSharp.Formats.Heif.Av1.Pipeline.LoopRestoration; + +/// +/// Applies the normative separable Wiener filter used by AV1 loop restoration. +/// +internal static class Av1WienerFilter +{ + /// + /// The number of coefficients in the padded Wiener convolution kernel. + /// + private const int FilterTapCount = 8; + + /// + /// The number of independent symmetric coefficients transmitted for each filter direction. + /// + private const int TransmittedCoefficientCount = 3; + + /// + /// The number of fractional bits in the Wiener filter coefficients. + /// + private const int FilterBits = 7; + + /// + /// The default first-pass rounding shift. + /// + private const int InitialHorizontalRoundBits = 3; + + /// + /// The number of intermediate rows required beyond the destination stripe height. + /// + private const int IntermediateRowExtension = FilterTapCount - 1; + + /// + /// Gets the number of intermediate samples required to filter a stripe. + /// + /// The destination stripe width. + /// The destination stripe height. + /// The required scratch-span length. + public static int GetScratchLength(int width, int height) => width * (height + IntermediateRowExtension); + + /// + /// Filters one restoration stripe from a source rectangle containing the required three-sample borders. + /// + /// The source rectangle beginning three samples above and left of the destination stripe. + /// The number of samples between source rows. + /// The destination span beginning at the restored stripe origin. + /// The number of samples between destination rows. + /// The stripe width in plane samples. + /// The stripe height in plane samples. + /// The encoded sample bit depth. + /// The three transmitted horizontal coefficients. + /// The three transmitted vertical coefficients. + /// Intermediate sample storage sized according to . + public static void FilterStripe( + ReadOnlySpan source, + int sourceStride, + Span destination, + int destinationStride, + int width, + int height, + int bitDepth, + ReadOnlySpan horizontalCoefficients, + ReadOnlySpan verticalCoefficients, + Span scratch) + { + Span horizontalFilter = stackalloc short[FilterTapCount]; + Span verticalFilter = stackalloc short[FilterTapCount]; + PopulateFilter(horizontalCoefficients, horizontalFilter); + PopulateFilter(verticalCoefficients, verticalFilter); + + int horizontalRoundBits = InitialHorizontalRoundBits; + int intermediateBitCount = bitDepth + FilterBits - horizontalRoundBits + 2; + if (intermediateBitCount > 16) + { + // Twelve-bit input would otherwise exceed the unsigned 16-bit intermediate used by + // the normative two-pass convolution, so AV1 transfers those excess bits to pass two. + horizontalRoundBits += intermediateBitCount - 16; + } + + int verticalRoundBits = (FilterBits * 2) - horizontalRoundBits; + int intermediateMaximum = (1 << (bitDepth + 1 + FilterBits - horizontalRoundBits)) - 1; + int intermediateHeight = height + IntermediateRowExtension; + int horizontalBias = 1 << (bitDepth + FilterBits - 1); + for (int row = 0; row < intermediateHeight; row++) + { + int sourceRowOffset = row * sourceStride; + int intermediateRowOffset = row * width; + for (int column = 0; column < width; column++) + { + int sourceOffset = sourceRowOffset + column; + int sum = DotProduct(source, sourceOffset, horizontalFilter); + + // The transmitted center coefficient excludes its implicit 128 contribution. + // Adding the unfiltered center sample here reconstructs the complete kernel. + sum += (source[sourceOffset + TransmittedCoefficientCount] << FilterBits) + horizontalBias; + int value = RoundPowerOfTwo(sum, horizontalRoundBits); + scratch[intermediateRowOffset + column] = (ushort)Av1Math.Clip3(0, intermediateMaximum, value); + } + } + + int maximumSample = (1 << bitDepth) - 1; + int verticalBias = 1 << (bitDepth + verticalRoundBits - 1); + for (int row = 0; row < height; row++) + { + int destinationRowOffset = row * destinationStride; + for (int column = 0; column < width; column++) + { + int sum = 0; + for (int tap = 0; tap < FilterTapCount; tap++) + { + sum += scratch[((row + tap) * width) + column] * verticalFilter[tap]; + } + + int center = scratch[((row + TransmittedCoefficientCount) * width) + column]; + sum += (center << FilterBits) - verticalBias; + destination[destinationRowOffset + column] = + (ushort)Av1Math.Clip3(0, maximumSample, RoundPowerOfTwo(sum, verticalRoundBits)); + } + } + } + + /// + /// Expands three transmitted symmetric coefficients into the padded eight-tap convolution kernel. + /// + /// The transmitted outer-to-inner coefficients. + /// The destination eight-tap kernel. + private static void PopulateFilter(ReadOnlySpan coefficients, Span filter) + { + int outer = coefficients[0]; + int middle = coefficients[1]; + int inner = coefficients[2]; + filter[0] = (short)outer; + filter[1] = (short)middle; + filter[2] = (short)inner; + filter[3] = (short)(-2 * (outer + middle + inner)); + filter[4] = (short)inner; + filter[5] = (short)middle; + filter[6] = (short)outer; + + // libaom stores a seven-tap Wiener kernel in the shared eight-tap interpolation shape. + filter[7] = 0; + } + + /// + /// Computes one signed eight-tap horizontal filter product. + /// + /// The source rectangle containing the requested samples. + /// The first source sample consumed by the filter. + /// The eight signed filter coefficients. + /// The unrounded signed filter sum. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int DotProduct(ReadOnlySpan source, int sourceOffset, ReadOnlySpan filter) + { + if (Vector128.IsHardwareAccelerated) + { + ref ushort sourceReference = ref MemoryMarshal.GetReference(source); + ref short filterReference = ref MemoryMarshal.GetReference(filter); + Vector128 samples = Vector128.LoadUnsafe(ref sourceReference, (nuint)sourceOffset).AsInt16(); + Vector128 coefficients = Vector128.LoadUnsafe(ref filterReference); + Vector128 pairSums = Vector128_.MultiplyAddAdjacent(samples, coefficients); + + // The shared helper provides the architecture-specific adjacent products; reducing its + // four 32-bit lanes scalarly avoids an additional platform-specific shuffle sequence. + return pairSums.GetElement(0) + pairSums.GetElement(1) + pairSums.GetElement(2) + pairSums.GetElement(3); + } + + int sum = 0; + for (int tap = 0; tap < FilterTapCount; tap++) + { + sum += source[sourceOffset + tap] * filter[tap]; + } + + return sum; + } + + /// + /// Rounds a signed fixed-point value to the requested lower precision. + /// + /// The signed fixed-point value. + /// The number of low bits to discard. + /// The rounded signed value. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int RoundPowerOfTwo(int value, int bitCount) + => (value + (1 << (bitCount - 1))) >> bitCount; +}