Reuse prepared AV1 residuals

HEIF_IMPLEMENTATION_PLAN.md

@@ -874,6 +874,8 @@ Encoder verification contract:
 
 - [x] Uniform four-by-four luma transform search now uses two compact reconstruction and coefficient views already available in the aligned mode-decision workspace. Legal transform trials write into the non-winning view and exchange span ownership only on strict rate-distortion improvement; the selected coefficients and strided reconstruction mosaic are published once after the type search so the next raster transform sees the required decoded edge. This removes reconstruction and coefficient copies on every improving transform without adding storage, changing tie order, or repeating a transform. All four representative palette, filter-intra, and effort-eight output hashes are unchanged, the focused Release set passes 13 of 13 cases, and current-main `aomdec` accepts every checked stream.
 
+- [x] Candidate distortion now follows the reference separation between immutable source residuals and reconstructed-pixel error. The forward-transform boundary accepts a read-only residual, so all transform types for one prediction reuse that block directly instead of copying it into transform scratch before every trial. The shared residual API now measures strided source-versus-reconstruction squared error with documented Vector512, Vector256, Vector128, and scalar traversal through ImageSharp's vector-count helpers, eliminating the former residual destination write and second reduction pass. The same path covers ordinary intra, filter intra, palette, chroma-from-luma, split transforms, and intra-block copy at 8, 10, and 12 bits without an allocation. Independent scalar, stride, tail, intrinsic-tier, and zero-allocation coverage passes with the 140-case focused encoder set; the complete non-HEVC HEIF/AV1 namespace passes 9,301 of 9,301. Representative palette, filter-intra, and effort-eight output hashes remain byte-identical, and current-main `aomdec` accepts every checked stream.
+
 ### 7. Write complete AVIF output
 
 - [~] The encoder-side AV1 codec configuration is now derived directly from the encoded sequence header and writes the fixed four-byte `av1C` record with empty `configOBUs`. The image payload retains the required sequence header, so the property introduces no sequence-header allocation, retention, or copy. Four production-header cases cover main, high, and professional profiles; 8-, 10-, and 12-bit precision; monochrome, 4:2:0, 4:2:2, and 4:4:4 sampling; exact fixed bytes; decoder reparsing; and header/property equivalence through direct net11 Release VSTest. Property-container emission and public AVIF activation remain open.

src/ImageSharp/Formats/Heif/Av1/Pipeline/Av1ResidualBuilder.cs

@@ -8,7 +8,7 @@ using System.Runtime.Intrinsics;
 namespace SixLabors.ImageSharp.Formats.Heif.Av1.Pipeline;
 
 /// <summary>
-/// Builds signed AV1 residual planes from source and prediction samples.
+/// Builds signed AV1 residual planes and measures sample-domain error.
 /// </summary>
 internal static partial class Av1ResidualBuilder
 {
@@ -56,6 +56,44 @@ internal static partial class Av1ResidualBuilder
         int height)
         => Subtract<ushort, UInt16Operator>(source, sourceStride, prediction, predictionStride, residual, residualStride, width, height);
 
+    /// <summary>
+    /// Calculates the exact squared error between strided 8-bit sample planes.
+    /// </summary>
+    /// <param name="source">The source samples.</param>
+    /// <param name="sourceStride">The source row stride.</param>
+    /// <param name="prediction">The prediction or reconstruction samples.</param>
+    /// <param name="predictionStride">The prediction row stride.</param>
+    /// <param name="width">The number of samples per row.</param>
+    /// <param name="height">The number of rows.</param>
+    /// <returns>The sum of squared sample differences.</returns>
+    public static long SumSquaredError(
+        ReadOnlySpan<byte> source,
+        int sourceStride,
+        ReadOnlySpan<byte> prediction,
+        int predictionStride,
+        int width,
+        int height)
+        => SumSquaredError<byte, ByteOperator>(source, sourceStride, prediction, predictionStride, width, height);
+
+    /// <summary>
+    /// Calculates the exact squared error between strided high-bit-depth sample planes.
+    /// </summary>
+    /// <param name="source">The source samples.</param>
+    /// <param name="sourceStride">The source row stride.</param>
+    /// <param name="prediction">The prediction or reconstruction samples.</param>
+    /// <param name="predictionStride">The prediction row stride.</param>
+    /// <param name="width">The number of samples per row.</param>
+    /// <param name="height">The number of rows.</param>
+    /// <returns>The sum of squared sample differences.</returns>
+    public static long SumSquaredError(
+        ReadOnlySpan<ushort> source,
+        int sourceStride,
+        ReadOnlySpan<ushort> prediction,
+        int predictionStride,
+        int width,
+        int height)
+        => SumSquaredError<ushort, UInt16Operator>(source, sourceStride, prediction, predictionStride, width, height);
+
     /// <summary>
     /// Sums the squares of a contiguous signed residual block.
     /// </summary>
@@ -76,10 +114,7 @@ internal static partial class Av1ResidualBuilder
             for (; vectorCount > 0; vectorCount--, offset += Vector512<short>.Count)
             {
                 Vector512<short> values = Unsafe.As<short, Vector512<short>>(ref Unsafe.Add(ref residualBase, offset));
-                Vector512<int> lower = Vector512.WidenLower(values);
-                Vector512<int> upper = Vector512.WidenUpper(values);
-                sum += Vector512.Sum(lower * lower);
-                sum += Vector512.Sum(upper * upper);
+                sum += SumSquares(values);
             }
         }
 
@@ -90,10 +125,7 @@ internal static partial class Av1ResidualBuilder
             for (; vectorCount > 0; vectorCount--, offset += Vector256<short>.Count)
             {
                 Vector256<short> values = Unsafe.As<short, Vector256<short>>(ref Unsafe.Add(ref residualBase, offset));
-                Vector256<int> lower = Vector256.WidenLower(values);
-                Vector256<int> upper = Vector256.WidenUpper(values);
-                sum += Vector256.Sum(lower * lower);
-                sum += Vector256.Sum(upper * upper);
+                sum += SumSquares(values);
             }
         }
 
@@ -104,10 +136,7 @@ internal static partial class Av1ResidualBuilder
             for (; vectorCount > 0; vectorCount--, offset += Vector128<short>.Count)
             {
                 Vector128<short> values = Unsafe.As<short, Vector128<short>>(ref Unsafe.Add(ref residualBase, offset));
-                Vector128<int> lower = Vector128.WidenLower(values);
-                Vector128<int> upper = Vector128.WidenUpper(values);
-                sum += Vector128.Sum(lower * lower);
-                sum += Vector128.Sum(upper * upper);
+                sum += SumSquares(values);
             }
         }
 
@@ -120,6 +149,103 @@ internal static partial class Av1ResidualBuilder
         return sum;
     }
 
+    private static long SumSquaredError<TSample, TOperator>(
+        ReadOnlySpan<TSample> source,
+        int sourceStride,
+        ReadOnlySpan<TSample> prediction,
+        int predictionStride,
+        int width,
+        int height)
+        where TSample : unmanaged
+        where TOperator : struct, IResidualOperator<TSample>
+    {
+        long sum = 0;
+        for (int y = 0; y < height; y++)
+        {
+            ReadOnlySpan<TSample> sourceRow = source.Slice(y * sourceStride, width);
+            ReadOnlySpan<TSample> predictionRow = prediction.Slice(y * predictionStride, width);
+            ref TSample sourceBase = ref MemoryMarshal.GetReference(sourceRow);
+            ref TSample predictionBase = ref MemoryMarshal.GetReference(predictionRow);
+            int x = 0;
+
+            // Descending hardware widths consume every complete vector before the scalar tail. Byte
+            // subtraction produces two widened residual vectors; high-bit-depth subtraction produces one.
+            if (Vector512.IsHardwareAccelerated)
+            {
+                nuint vectorCount = sourceRow.Vector512Count<TSample>();
+
+                for (; vectorCount > 0; vectorCount--, x += Vector512<TSample>.Count)
+                {
+                    Vector512<TSample> sourceVector =
+                        Unsafe.As<TSample, Vector512<TSample>>(ref Unsafe.Add(ref sourceBase, x));
+
+                    Vector512<TSample> predictionVector =
+                        Unsafe.As<TSample, Vector512<TSample>>(ref Unsafe.Add(ref predictionBase, x));
+
+                    Vector512<short> lower = TOperator.Subtract(sourceVector, predictionVector, out Vector512<short> upper);
+                    sum += SumSquares(lower);
+                    if (Vector512<TSample>.Count != Vector512<short>.Count)
+                    {
+                        sum += SumSquares(upper);
+                    }
+                }
+            }
+
+            if (Vector256.IsHardwareAccelerated)
+            {
+                nuint vectorCount = sourceRow[x..].Vector256Count<TSample>();
+
+                for (; vectorCount > 0; vectorCount--, x += Vector256<TSample>.Count)
+                {
+                    Vector256<TSample> sourceVector =
+                        Unsafe.As<TSample, Vector256<TSample>>(ref Unsafe.Add(ref sourceBase, x));
+
+                    Vector256<TSample> predictionVector =
+                        Unsafe.As<TSample, Vector256<TSample>>(ref Unsafe.Add(ref predictionBase, x));
+
+                    Vector256<short> lower = TOperator.Subtract(sourceVector, predictionVector, out Vector256<short> upper);
+                    sum += SumSquares(lower);
+                    if (Vector256<TSample>.Count != Vector256<short>.Count)
+                    {
+                        sum += SumSquares(upper);
+                    }
+                }
+            }
+
+            if (Vector128.IsHardwareAccelerated)
+            {
+                nuint vectorCount = sourceRow[x..].Vector128Count<TSample>();
+
+                for (; vectorCount > 0; vectorCount--, x += Vector128<TSample>.Count)
+                {
+                    Vector128<TSample> sourceVector =
+                        Unsafe.As<TSample, Vector128<TSample>>(ref Unsafe.Add(ref sourceBase, x));
+
+                    Vector128<TSample> predictionVector =
+                        Unsafe.As<TSample, Vector128<TSample>>(ref Unsafe.Add(ref predictionBase, x));
+
+                    Vector128<short> lower = TOperator.Subtract(sourceVector, predictionVector, out Vector128<short> upper);
+                    sum += SumSquares(lower);
+                    if (Vector128<TSample>.Count != Vector128<short>.Count)
+                    {
+                        sum += SumSquares(upper);
+                    }
+                }
+            }
+
+            for (; x < width; x++)
+            {
+                int difference = TOperator.Subtract(
+                    Unsafe.Add(ref sourceBase, x),
+                    Unsafe.Add(ref predictionBase, x));
+
+                sum += difference * difference;
+            }
+        }
+
+        return sum;
+    }
+
     private static void Subtract<TSample, TOperator>(
         ReadOnlySpan<TSample> source,
         int sourceStride,
@@ -211,4 +337,28 @@ internal static partial class Av1ResidualBuilder
             }
         }
     }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static long SumSquares(Vector128<short> values)
+    {
+        Vector128<int> lower = Vector128.WidenLower(values);
+        Vector128<int> upper = Vector128.WidenUpper(values);
+        return (long)Vector128.Sum(lower * lower) + Vector128.Sum(upper * upper);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static long SumSquares(Vector256<short> values)
+    {
+        Vector256<int> lower = Vector256.WidenLower(values);
+        Vector256<int> upper = Vector256.WidenUpper(values);
+        return (long)Vector256.Sum(lower * lower) + Vector256.Sum(upper * upper);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static long SumSquares(Vector512<short> values)
+    {
+        Vector512<int> lower = Vector512.WidenLower(values);
+        Vector512<int> upper = Vector512.WidenUpper(values);
+        return (long)Vector512.Sum(lower * lower) + Vector512.Sum(upper * upper);
+    }
 }

src/ImageSharp/Formats/Heif/Av1/Pipeline/Av1TransformBlockEncoder.cs

@@ -145,17 +145,15 @@ internal static class Av1TransformBlockEncoder
             plane,
             ref state);
 
-        Av1ResidualBuilder.Subtract(
+        long distortion = Av1ResidualBuilder.SumSquaredError(
             sourceSamples,
             source.Stride,
             reconstruction,
             width,
-            workspace.Residual,
-            width,
             width,
             height);
 
-        return Av1ResidualBuilder.SumSquares(workspace.Residual[..transformSize.GetSize2d()]) << 4;
+        return distortion << 4;
     }
 
     /// <summary>
@@ -199,16 +197,16 @@ internal static class Av1TransformBlockEncoder
         int sampleCount = transformSize.GetSize2d();
         ReadOnlySpan<byte> sourceSamples = GetPlaneSpan(source, blockOrigin);
 
-        // Each transform trial mutates reconstruction and residual scratch, so restore both prepared inputs.
+        // Each transform trial overwrites reconstruction but consumes the prepared residual read-only.
         // Row copies preserve a larger candidate surface without materializing a second compact block.
         for (int row = 0; row < height; row++)
         {
             prediction.Slice(row * width, width).CopyTo(reconstruction.Slice(row * reconstructionStride, width));
         }
 
-        residual[..sampleCount].CopyTo(workspace.Residual);
         EncodeLossy(
             workspace,
+            residual[..sampleCount],
             quantizedCoefficients,
             transformSize,
             transformType,
@@ -232,17 +230,15 @@ internal static class Av1TransformBlockEncoder
                 workspace.TransformWorkspace);
         }
 
-        Av1ResidualBuilder.Subtract(
+        long distortion = Av1ResidualBuilder.SumSquaredError(
             sourceSamples,
             source.Stride,
             reconstruction,
             reconstructionStride,
-            workspace.Residual,
-            width,
             width,
             height);
 
-        return Av1ResidualBuilder.SumSquares(workspace.Residual[..sampleCount]) << 4;
+        return distortion << 4;
     }
 
     /// <summary>
@@ -322,17 +318,15 @@ internal static class Av1TransformBlockEncoder
         }
 
         // Final distortion is measured against the samples a decoder reconstructs, not the unquantized predictor.
-        Av1ResidualBuilder.Subtract(
+        long distortion = Av1ResidualBuilder.SumSquaredError(
             sourceSamples,
             source.Stride,
             reconstruction,
             width,
-            workspace.Residual,
-            width,
             width,
             height);
 
-        return Av1ResidualBuilder.SumSquares(workspace.Residual[..transformSize.GetSize2d()]) << 4;
+        return distortion << 4;
     }
 
     /// <summary>
@@ -470,17 +464,14 @@ internal static class Av1TransformBlockEncoder
             bitDepth,
             ref state);
 
-        Av1ResidualBuilder.Subtract(
+        long distortion = Av1ResidualBuilder.SumSquaredError(
             sourceSamples,
             source.Stride,
             reconstruction,
             width,
-            workspace.Residual,
-            width,
             width,
             height);
 
-        long distortion = Av1ResidualBuilder.SumSquares(workspace.Residual[..transformSize.GetSize2d()]);
         int shift = (bitDepth.GetBitCount() - 8) * 2;
         long normalizedDistortion = shift == 0
             ? distortion
@@ -532,16 +523,16 @@ internal static class Av1TransformBlockEncoder
         int sampleCount = transformSize.GetSize2d();
         ReadOnlySpan<ushort> sourceSamples = GetPlaneSpan(source, blockOrigin);
 
-        // Each transform trial mutates reconstruction and residual scratch, so restore both prepared inputs.
+        // Each transform trial overwrites reconstruction but consumes the prepared residual read-only.
         // Row copies preserve a larger candidate surface without materializing a second compact block.
         for (int row = 0; row < height; row++)
         {
             prediction.Slice(row * width, width).CopyTo(reconstruction.Slice(row * reconstructionStride, width));
         }
 
-        residual[..sampleCount].CopyTo(workspace.Residual);
         EncodeLossy(
             workspace,
+            residual[..sampleCount],
             quantizedCoefficients,
             transformSize,
             transformType,
@@ -566,17 +557,14 @@ internal static class Av1TransformBlockEncoder
                 workspace.TransformWorkspace);
         }
 
-        Av1ResidualBuilder.Subtract(
+        long distortion = Av1ResidualBuilder.SumSquaredError(
             sourceSamples,
             source.Stride,
             reconstruction,
             reconstructionStride,
-            workspace.Residual,
-            width,
             width,
             height);
 
-        long distortion = Av1ResidualBuilder.SumSquares(workspace.Residual[..sampleCount]);
         int shift = (bitDepth.GetBitCount() - 8) * 2;
         long normalizedDistortion = shift == 0
             ? distortion
@@ -672,17 +660,14 @@ internal static class Av1TransformBlockEncoder
                 workspace.TransformWorkspace);
         }
 
-        Av1ResidualBuilder.Subtract(
+        long distortion = Av1ResidualBuilder.SumSquaredError(
             sourceSamples,
             source.Stride,
             reconstruction,
             width,
-            workspace.Residual,
-            width,
             width,
             height);
 
-        long distortion = Av1ResidualBuilder.SumSquares(workspace.Residual[..transformSize.GetSize2d()]);
         int shift = (bitDepth.GetBitCount() - 8) * 2;
         long normalizedDistortion = shift == 0
             ? distortion
@@ -1047,16 +1032,39 @@ internal static class Av1TransformBlockEncoder
         int acDeltaQ,
         Av1BitDepth bitDepth,
         ref Av1EncoderTransformBlockState state)
+        => EncodeLossy(
+            workspace,
+            workspace.Residual,
+            quantizedCoefficients,
+            transformSize,
+            transformType,
+            qIndex,
+            dcDeltaQ,
+            acDeltaQ,
+            bitDepth,
+            ref state);
+
+    private static void EncodeLossy(
+        Av1EncoderBlockWorkspace workspace,
+        ReadOnlySpan<short> residual,
+        Span<int> quantizedCoefficients,
+        Av1TransformSize transformSize,
+        Av1TransformType transformType,
+        int qIndex,
+        int dcDeltaQ,
+        int acDeltaQ,
+        Av1BitDepth bitDepth,
+        ref Av1EncoderTransformBlockState state)
     {
         int coefficientCount = transformSize.GetAdjusted().GetSize2d();
         Span<int> transformed = workspace.TransformCoefficients[..coefficientCount];
         Span<int> quantized = quantizedCoefficients[..coefficientCount];
         Span<int> dequantized = workspace.DequantizedCoefficients[..coefficientCount];
 
-        // The encoder keeps transformed, quantized, and reconstructed coefficients separate because mode decision
-        // consumes all three while only the quantized values survive in the frame coefficient owner.
+        // The forward transform reads the prepared residual without changing it, so every type candidate can
+        // reuse one source-minus-prediction block. Separate coefficient spans preserve each later representation.
         Av1ForwardTransformer.Transform2d(
-            workspace.Residual,
+            residual,
             transformed,
             (uint)transformSize.GetWidth(),
             transformType,

src/ImageSharp/Formats/Heif/Av1/Transform/Av1ForwardTransformer.cs

@@ -31,7 +31,7 @@ internal static partial class Av1ForwardTransformer
     /// <param name="bitDepth">The source sample bit depth.</param>
     /// <param name="workspace">The reusable workspace owned by the containing encode operation.</param>
     public static void Transform2d(
-        Span<short> input,
+        ReadOnlySpan<short> input,
         Span<int> coefficients,
         uint stride,
         Av1TransformType transformType,
@@ -55,7 +55,7 @@ internal static partial class Av1ForwardTransformer
     /// <param name="config">The resolved transform functions, shifts, and axis orientation.</param>
     /// <param name="workspace">The reusable transform workspace.</param>
     private static void DispatchColumn(
-        Span<short> input,
+        ReadOnlySpan<short> input,
         Span<int> coefficients,
         uint stride,
         int bitDepth,
@@ -116,7 +116,7 @@ internal static partial class Av1ForwardTransformer
     /// <param name="config">The resolved transform functions, shifts, and axis orientation.</param>
     /// <param name="workspace">The reusable transform workspace.</param>
     private static void DispatchRow<TColumnOperator>(
-        Span<short> input,
+        ReadOnlySpan<short> input,
         Span<int> coefficients,
         uint stride,
         int bitDepth,
@@ -179,7 +179,7 @@ internal static partial class Av1ForwardTransformer
     /// <param name="config">The resolved transform functions, shifts, and axis orientation.</param>
     /// <param name="workspace">The reusable transform workspace.</param>
     private static void Transform2d<TColumnOperator, TRowOperator>(
-        Span<short> input,
+        ReadOnlySpan<short> input,
         Span<int> coefficients,
         uint stride,
         int bitDepth,
@@ -210,7 +210,7 @@ internal static partial class Av1ForwardTransformer
     /// <param name="config">The resolved transform functions, shifts, and axis orientation.</param>
     /// <param name="workspace">The reusable transform workspace.</param>
     private static void TransformPacked<TColumnOperator, TRowOperator>(
-        Span<short> input,
+        ReadOnlySpan<short> input,
         Span<int> coefficients,
         uint stride,
         ref Av1Transform2dFlipConfiguration config,
@@ -324,7 +324,7 @@ internal static partial class Av1ForwardTransformer
     /// <param name="config">The resolved transform functions, shifts, and axis orientation.</param>
     /// <param name="workspace">The reusable transform workspace.</param>
     private static void TransformExpanded<TColumnOperator, TRowOperator>(
-        Span<short> input,
+        ReadOnlySpan<short> input,
         Span<int> coefficients,
         uint stride,
         ref Av1Transform2dFlipConfiguration config,
@@ -395,7 +395,7 @@ internal static partial class Av1ForwardTransformer
     /// <param name="flipUpsideDown">Whether to reverse the input row order.</param>
     /// <param name="flipLeftToRight">Whether to reverse the samples within each row.</param>
     private static void LoadPacked(
-        Span<short> input,
+        ReadOnlySpan<short> input,
         uint inputStride,
         ref short destination,
         int destinationStride,
@@ -476,7 +476,7 @@ internal static partial class Av1ForwardTransformer
     /// <param name="flipUpsideDown">Whether to reverse the input row order.</param>
     /// <param name="flipLeftToRight">Whether to reverse the samples within each row.</param>
     private static void LoadExpanded(
-        Span<short> input,
+        ReadOnlySpan<short> input,
         uint inputStride,
         ref int destination,
         int destinationStride,

tests/ImageSharp.Tests/Formats/Heif/Av1/Av1ResidualBuilderTests.cs

@@ -106,7 +106,7 @@ public class Av1ResidualBuilderTests
         => FeatureTestRunner.RunWithHwIntrinsicsFeature(ValidateSumSquares, ResidualConfigurations);
 
     /// <summary>
-    /// Verifies that repeated maximum-transform residual construction uses only caller-owned buffers.
+    /// Verifies that repeated maximum-transform residual and squared-error processing uses only caller-owned buffers.
     /// </summary>
     [Fact]
     public void ResidualConstructionDoesNotAllocate()
@@ -124,15 +124,34 @@ public class Av1ResidualBuilderTests
 
         Av1ResidualBuilder.Subtract(source, width, prediction, width, residual, width, width, height);
         Av1ResidualBuilder.Subtract(highBitDepthSource, width, highBitDepthPrediction, width, highBitDepthResidual, width, width, height);
-
+        _ = Av1ResidualBuilder.SumSquaredError(source, width, prediction, width, width, height);
+        _ = Av1ResidualBuilder.SumSquaredError(
+            highBitDepthSource,
+            width,
+            highBitDepthPrediction,
+            width,
+            width,
+            height);
+
+        long sum = 0;
         long before = GC.GetAllocatedBytesForCurrentThread();
         for (int iteration = 0; iteration < 32; iteration++)
         {
             Av1ResidualBuilder.Subtract(source, width, prediction, width, residual, width, width, height);
             Av1ResidualBuilder.Subtract(highBitDepthSource, width, highBitDepthPrediction, width, highBitDepthResidual, width, width, height);
+            sum += Av1ResidualBuilder.SumSquaredError(source, width, prediction, width, width, height);
+            sum += Av1ResidualBuilder.SumSquaredError(
+                highBitDepthSource,
+                width,
+                highBitDepthPrediction,
+                width,
+                width,
+                height);
         }
 
-        Assert.Equal(0, GC.GetAllocatedBytesForCurrentThread() - before);
+        long allocated = GC.GetAllocatedBytesForCurrentThread() - before;
+        Assert.True(sum > 0);
+        Assert.Equal(0, allocated);
     }
 
     private static void ValidateSumSquares()
@@ -186,6 +205,26 @@ public class Av1ResidualBuilderTests
         Av1ResidualBuilder.Subtract(sourcePlane, sourceStride, predictionPlane, predictionStride, actualPlane, residualStride, width, height);
 
         Assert.Equal(expected, actual);
+
+        long expectedSquaredError = 0;
+        for (int y = 0; y < height; y++)
+        {
+            for (int x = 0; x < width; x++)
+            {
+                int difference = sourcePlane[(y * sourceStride) + x] - predictionPlane[(y * predictionStride) + x];
+                expectedSquaredError += difference * difference;
+            }
+        }
+
+        Assert.Equal(
+            expectedSquaredError,
+            Av1ResidualBuilder.SumSquaredError(
+                sourcePlane,
+                sourceStride,
+                predictionPlane,
+                predictionStride,
+                width,
+                height));
     }
 
     private static void ValidateUInt16Residuals(int maximumSample)
@@ -215,6 +254,26 @@ public class Av1ResidualBuilderTests
         Av1ResidualBuilder.Subtract(sourcePlane, sourceStride, predictionPlane, predictionStride, actualPlane, residualStride, width, height);
 
         Assert.Equal(expected, actual);
+
+        long expectedSquaredError = 0;
+        for (int y = 0; y < height; y++)
+        {
+            for (int x = 0; x < width; x++)
+            {
+                int difference = sourcePlane[(y * sourceStride) + x] - predictionPlane[(y * predictionStride) + x];
+                expectedSquaredError += difference * difference;
+            }
+        }
+
+        Assert.Equal(
+            expectedSquaredError,
+            Av1ResidualBuilder.SumSquaredError(
+                sourcePlane,
+                sourceStride,
+                predictionPlane,
+                predictionStride,
+                width,
+                height));
     }
 
     private static void FillBytePlanes(Span<byte> source, int sourceStride, Span<byte> prediction, int predictionStride, int width, int height)