diff --git a/HEIF_IMPLEMENTATION_PLAN.md b/HEIF_IMPLEMENTATION_PLAN.md
index e00557367f..2886165b5e 100644
--- a/HEIF_IMPLEMENTATION_PLAN.md
+++ b/HEIF_IMPLEMENTATION_PLAN.md
@@ -856,6 +856,7 @@ Encoder verification contract:
 - [~] Retained encoder palette state and production palette writing now mirror current libaom's 50-byte palette-mode contents, separate luma and shared-chroma sizes, three eight-color planes, above-and-left sorted cache, 64-sample above-cache boundary, mode contexts, palette colors, color-index maps, and syntax order. The current block keeps one inline value in the reusable superblock workspace; only the 4x4-granularity top and left picture edges retain copies for later blocks. For a 3840x2160 tile these edges occupy about 73.4 KiB instead of about 6.2 MiB for a 50-byte palette value attached to every 8x8 mode allocation. Luma and chroma index maps share one lazily allocated 32 KiB owner containing two 128x128 maps, so the palette-disabled production path retains no map owner. The compact final-block workspace falls from about 10.3 KiB to about 8.3 KiB. The writer caps map traversal to the coded plane count, writes maps before transform syntax, and publishes palette edges only after the current block has consumed preceding contexts. Eight focused size, alignment, ownership, cache-boundary, round-trip, map-consumption, and edge-publication cases pass; all 114 palette cases, all 1,942 entropy cases, and all 8,996 HEIF/AV1 cases pass direct foreground net11 Release VSTest. The exact net11 Release rebuild reports 1,050 solution warnings and zero errors; Roslynk reports zero compiler errors and no diagnostics in the touched files. The current-main reference is `a40ed1ea9e4ecc3df58a5bccb76623f2c94ae727`. Production payloads remain unchanged because palette candidate generation is still disabled; that live rate-distortion search is the next checkpoint.
 - [~] Luma palette clustering now follows current libaom's one-dimensional search primitive exactly: equal-interval midpoint initialization, first-color tie order, rounded centroid means, deterministic empty-cluster replacement, the 50-iteration limit, and retention of the preceding state when distortion increases. Nearest-color assignment improves on libaom's AVX2 implementation by dispatching Vector512, Vector256, Vector128, then scalar through ImageSharp's shared vector-count helpers. The primitive uses only bounded stack scratch and introduces no allocator rent, managed array, or per-row copy. Three independent tests cover exact centroid convergence, initialization order, 12-bit nearest-color distortion, destination bounds, and every hardware-intrinsic tier. The complete AVIF set passes 8,930 of 8,930 cases and the HEIF set passes 230 of 230 cases through direct foreground net11 Release VSTest. The exact net11 Release rebuild reports 1,050 solution warnings and zero errors, and Roslynk reports zero compiler errors. Candidate enumeration, palette-cache snapping, transform RD selection, and production activation remain in the open luma-palette checkpoint.
 - [~] Live luma palette selection now follows current libaom's dominant-color and one-dimensional K-means candidate families, cache-bias threshold, sorted duplicate removal, active-edge map extension, and strict winner tie order. It improves on speed-configured libaom by evaluating both candidate families at every legal 2-through-8 size without early pruning, then exhaustively evaluates every legal transform using the existing SIMD prediction, residual, transform, quantization, and reconstruction operators. Candidate storage remains bounded stack memory; the reusable 32 KiB map owner is allocated only when an eligible block enters palette search. A production tile test proves that full 8x8 and clipped 5x3 blocks at 8 and 12 bits select exact colors and indices, extend the visible edges through coded padding, reconstruct every sample without coefficients, and emit a nonempty tile. The complete 57-case intra-superblock set, 8,931-case AVIF set, and 230-case HEIF set pass direct foreground net11 Release VSTest. The exact Release test-project build reports 1,992 baseline warnings and zero errors; Roslynk reports zero compiler errors and no touched-file analyzer warnings. Production frame activation remains gated until chroma palette mode and its rate accounting are complete.
+- [~] Paired chroma palette clustering now preserves current libaom's squared two-component distance, first-centroid tie order, independently rounded U/V means, paired deterministic empty-cluster replacement, preceding-state retention on increased distortion, and 50-iteration limit. Keeping the source planes separate avoids interleave/deinterleave copies and improves on libaom's AVX2 ceiling with Vector512, Vector256, Vector128, then scalar dispatch through ImageSharp's shared vector-count helpers. Three independent tests cover exact paired convergence, midpoint initialization, 12-bit distance and index parity, untouched destination bounds, and every intrinsic tier. The exact Release test-project build reports 1,992 baseline warnings and zero errors; the focused three-case set, complete 8,934-case AVIF set, and complete 230-case HEIF set pass direct foreground net11 Release VSTest. Roslynk reports zero compiler errors and no touched-file analyzer warnings. Candidate integration and production activation remain in the open chroma-palette checkpoint.
 - [x] The expanded checkpoint exposed a pre-existing transform-block test that asserted uninitialized pooled padding was zero. The test now initializes the complete physical luma plane with a sentinel and proves the block operation leaves both adjacent padding samples unchanged. The exact net11 Release rebuild remains at 1,005 baseline warnings and zero errors, the focused allocator-order set passes 30 of 30 cases, and the complete HEIF/AV1 namespace passes 8,859 of 8,859 direct VSTest cases with zero failures or skips.
 - [x] Combined-frame OBU output now counts the byte-aligned frame and tile-group headers, non-final tile-size fields, and owned tile payloads before emitting the OBU size. It retains only the small allocator-owned header scratch and writes each entropy-coded tile span directly from its detached owner, removing the second file-sized allocator rent and complete-payload copy. A 64 KiB regression proves exactly one sub-payload-sized byte rent with a balanced return and verifies the exact streamed tile tail; the existing two-tile round trip proves size-prefix and ordering parity. The focused writer and production-frame set passes 32 of 32 direct net11 VSTest cases, current-main `aomdec` accepts all 29 generated native-format payloads, and the complete HEIF/AV1 namespace passes 8,860 of 8,860 cases with zero failures or skips.
 - [x] Finalized fixed-block decisions now set the block-level transform-skip flag only when every retained luma and coded chroma transform has zero EOB, matching current libaom's conjunction of per-plane skip state. The previous always-false flag produced legal but redundant non-skip and zero-coefficient syntax. Monochrome and 4:2:0 regressions prove both branches from actual coefficient state; the focused decision and production-frame set passes 32 of 32 direct net11 VSTest cases. Current-main `aomdec` accepts all 29 regenerated payloads, the recorded decoded-frame MD5s are unchanged, and affected 16x16 constant 8-bit and 10-bit payloads are one byte smaller. The complete HEIF/AV1 namespace passes 8,862 of 8,862 cases with zero failures or skips.
diff --git a/src/ImageSharp/Formats/Heif/Av1/Pipeline/Av1PaletteKMeans2D.cs b/src/ImageSharp/Formats/Heif/Av1/Pipeline/Av1PaletteKMeans2D.cs
new file mode 100644
index 0000000000..74ebb5af44
--- /dev/null
+++ b/src/ImageSharp/Formats/Heif/Av1/Pipeline/Av1PaletteKMeans2D.cs
@@ -0,0 +1,358 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Runtime.InteropServices;
+using System.Runtime.Intrinsics;
+
+namespace SixLabors.ImageSharp.Formats.Heif.Av1.Pipeline;
+
+/// <summary>
+/// Assigns paired chroma samples to AV1 palette colors and refines two-dimensional palette centroids.
+/// </summary>
+internal static class Av1PaletteKMeans2D
+{
+    /// <summary>
+    /// Assigns every chroma pair to its nearest palette color.
+    /// </summary>
+    /// <param name="firstSamples">The active first-plane samples.</param>
+    /// <param name="secondSamples">The active second-plane samples.</param>
+    /// <param name="firstCentroids">The first-plane palette colors.</param>
+    /// <param name="secondCentroids">The second-plane palette colors.</param>
+    /// <param name="indices">The destination palette indices.</param>
+    /// <returns>The sum of squared two-plane sample-to-centroid distances.</returns>
+    public static long AssignIndices(
+        ReadOnlySpan<short> firstSamples,
+        ReadOnlySpan<short> secondSamples,
+        ReadOnlySpan<short> firstCentroids,
+        ReadOnlySpan<short> secondCentroids,
+        Span<byte> indices)
+    {
+        Span<int> distanceScratch = stackalloc int[Vector512<short>.Count];
+        Span<int> indexScratch = stackalloc int[Vector512<short>.Count];
+        ref short firstSampleBase = ref MemoryMarshal.GetReference(firstSamples);
+        ref short secondSampleBase = ref MemoryMarshal.GetReference(secondSamples);
+        int offset = 0;
+        long distortion = 0;
+
+        if (Vector512.IsHardwareAccelerated)
+        {
+            nuint vectorCount = Numerics.Vector512Count(firstSamples);
+            for (; vectorCount > 0; vectorCount--, offset += Vector512<short>.Count)
+            {
+                Vector512<short> firstSample = Vector512.LoadUnsafe(ref firstSampleBase, (nuint)offset);
+                Vector512<short> secondSample = Vector512.LoadUnsafe(ref secondSampleBase, (nuint)offset);
+                Vector512<short> firstDifference = firstSample - Vector512.Create(firstCentroids[0]);
+                Vector512<short> secondDifference = secondSample - Vector512.Create(secondCentroids[0]);
+                (Vector512<int> firstLower, Vector512<int> firstUpper) = Vector512.Widen(firstDifference);
+                (Vector512<int> secondLower, Vector512<int> secondUpper) = Vector512.Widen(secondDifference);
+                Vector512<int> bestDistanceLower = (firstLower * firstLower) + (secondLower * secondLower);
+                Vector512<int> bestDistanceUpper = (firstUpper * firstUpper) + (secondUpper * secondUpper);
+                Vector512<int> bestIndexLower = Vector512<int>.Zero;
+                Vector512<int> bestIndexUpper = Vector512<int>.Zero;
+                for (int centroidIndex = 1; centroidIndex < firstCentroids.Length; centroidIndex++)
+                {
+                    firstDifference = firstSample - Vector512.Create(firstCentroids[centroidIndex]);
+                    secondDifference = secondSample - Vector512.Create(secondCentroids[centroidIndex]);
+                    (firstLower, firstUpper) = Vector512.Widen(firstDifference);
+                    (secondLower, secondUpper) = Vector512.Widen(secondDifference);
+                    Vector512<int> distanceLower = (firstLower * firstLower) + (secondLower * secondLower);
+                    Vector512<int> distanceUpper = (firstUpper * firstUpper) + (secondUpper * secondUpper);
+                    Vector512<int> replaceLower = Vector512.LessThan(distanceLower, bestDistanceLower);
+                    Vector512<int> replaceUpper = Vector512.LessThan(distanceUpper, bestDistanceUpper);
+                    bestDistanceLower = Vector512.ConditionalSelect(replaceLower, distanceLower, bestDistanceLower);
+                    bestDistanceUpper = Vector512.ConditionalSelect(replaceUpper, distanceUpper, bestDistanceUpper);
+                    bestIndexLower = Vector512.ConditionalSelect(
+                        replaceLower,
+                        Vector512.Create(centroidIndex),
+                        bestIndexLower);
+
+                    bestIndexUpper = Vector512.ConditionalSelect(
+                        replaceUpper,
+                        Vector512.Create(centroidIndex),
+                        bestIndexUpper);
+                }
+
+                bestDistanceLower.CopyTo(distanceScratch);
+                bestDistanceUpper.CopyTo(distanceScratch[Vector512<int>.Count..]);
+                bestIndexLower.CopyTo(indexScratch);
+                bestIndexUpper.CopyTo(indexScratch[Vector512<int>.Count..]);
+                for (int lane = 0; lane < Vector512<short>.Count; lane++)
+                {
+                    indices[offset + lane] = (byte)indexScratch[lane];
+                    distortion += distanceScratch[lane];
+                }
+            }
+        }
+
+        if (Vector256.IsHardwareAccelerated)
+        {
+            nuint vectorCount = Numerics.Vector256Count(firstSamples[offset..]);
+            for (; vectorCount > 0; vectorCount--, offset += Vector256<short>.Count)
+            {
+                Vector256<short> firstSample = Vector256.LoadUnsafe(ref firstSampleBase, (nuint)offset);
+                Vector256<short> secondSample = Vector256.LoadUnsafe(ref secondSampleBase, (nuint)offset);
+                Vector256<short> firstDifference = firstSample - Vector256.Create(firstCentroids[0]);
+                Vector256<short> secondDifference = secondSample - Vector256.Create(secondCentroids[0]);
+                (Vector256<int> firstLower, Vector256<int> firstUpper) = Vector256.Widen(firstDifference);
+                (Vector256<int> secondLower, Vector256<int> secondUpper) = Vector256.Widen(secondDifference);
+                Vector256<int> bestDistanceLower = (firstLower * firstLower) + (secondLower * secondLower);
+                Vector256<int> bestDistanceUpper = (firstUpper * firstUpper) + (secondUpper * secondUpper);
+                Vector256<int> bestIndexLower = Vector256<int>.Zero;
+                Vector256<int> bestIndexUpper = Vector256<int>.Zero;
+                for (int centroidIndex = 1; centroidIndex < firstCentroids.Length; centroidIndex++)
+                {
+                    firstDifference = firstSample - Vector256.Create(firstCentroids[centroidIndex]);
+                    secondDifference = secondSample - Vector256.Create(secondCentroids[centroidIndex]);
+                    (firstLower, firstUpper) = Vector256.Widen(firstDifference);
+                    (secondLower, secondUpper) = Vector256.Widen(secondDifference);
+                    Vector256<int> distanceLower = (firstLower * firstLower) + (secondLower * secondLower);
+                    Vector256<int> distanceUpper = (firstUpper * firstUpper) + (secondUpper * secondUpper);
+                    Vector256<int> replaceLower = Vector256.LessThan(distanceLower, bestDistanceLower);
+                    Vector256<int> replaceUpper = Vector256.LessThan(distanceUpper, bestDistanceUpper);
+                    bestDistanceLower = Vector256.ConditionalSelect(replaceLower, distanceLower, bestDistanceLower);
+                    bestDistanceUpper = Vector256.ConditionalSelect(replaceUpper, distanceUpper, bestDistanceUpper);
+                    bestIndexLower = Vector256.ConditionalSelect(
+                        replaceLower,
+                        Vector256.Create(centroidIndex),
+                        bestIndexLower);
+
+                    bestIndexUpper = Vector256.ConditionalSelect(
+                        replaceUpper,
+                        Vector256.Create(centroidIndex),
+                        bestIndexUpper);
+                }
+
+                bestDistanceLower.CopyTo(distanceScratch);
+                bestDistanceUpper.CopyTo(distanceScratch[Vector256<int>.Count..]);
+                bestIndexLower.CopyTo(indexScratch);
+                bestIndexUpper.CopyTo(indexScratch[Vector256<int>.Count..]);
+                for (int lane = 0; lane < Vector256<short>.Count; lane++)
+                {
+                    indices[offset + lane] = (byte)indexScratch[lane];
+                    distortion += distanceScratch[lane];
+                }
+            }
+        }
+
+        if (Vector128.IsHardwareAccelerated)
+        {
+            nuint vectorCount = Numerics.Vector128Count(firstSamples[offset..]);
+            for (; vectorCount > 0; vectorCount--, offset += Vector128<short>.Count)
+            {
+                Vector128<short> firstSample = Vector128.LoadUnsafe(ref firstSampleBase, (nuint)offset);
+                Vector128<short> secondSample = Vector128.LoadUnsafe(ref secondSampleBase, (nuint)offset);
+                Vector128<short> firstDifference = firstSample - Vector128.Create(firstCentroids[0]);
+                Vector128<short> secondDifference = secondSample - Vector128.Create(secondCentroids[0]);
+                (Vector128<int> firstLower, Vector128<int> firstUpper) = Vector128.Widen(firstDifference);
+                (Vector128<int> secondLower, Vector128<int> secondUpper) = Vector128.Widen(secondDifference);
+                Vector128<int> bestDistanceLower = (firstLower * firstLower) + (secondLower * secondLower);
+                Vector128<int> bestDistanceUpper = (firstUpper * firstUpper) + (secondUpper * secondUpper);
+                Vector128<int> bestIndexLower = Vector128<int>.Zero;
+                Vector128<int> bestIndexUpper = Vector128<int>.Zero;
+                for (int centroidIndex = 1; centroidIndex < firstCentroids.Length; centroidIndex++)
+                {
+                    firstDifference = firstSample - Vector128.Create(firstCentroids[centroidIndex]);
+                    secondDifference = secondSample - Vector128.Create(secondCentroids[centroidIndex]);
+                    (firstLower, firstUpper) = Vector128.Widen(firstDifference);
+                    (secondLower, secondUpper) = Vector128.Widen(secondDifference);
+                    Vector128<int> distanceLower = (firstLower * firstLower) + (secondLower * secondLower);
+                    Vector128<int> distanceUpper = (firstUpper * firstUpper) + (secondUpper * secondUpper);
+                    Vector128<int> replaceLower = Vector128.LessThan(distanceLower, bestDistanceLower);
+                    Vector128<int> replaceUpper = Vector128.LessThan(distanceUpper, bestDistanceUpper);
+                    bestDistanceLower = Vector128.ConditionalSelect(replaceLower, distanceLower, bestDistanceLower);
+                    bestDistanceUpper = Vector128.ConditionalSelect(replaceUpper, distanceUpper, bestDistanceUpper);
+                    bestIndexLower = Vector128.ConditionalSelect(
+                        replaceLower,
+                        Vector128.Create(centroidIndex),
+                        bestIndexLower);
+
+                    bestIndexUpper = Vector128.ConditionalSelect(
+                        replaceUpper,
+                        Vector128.Create(centroidIndex),
+                        bestIndexUpper);
+                }
+
+                bestDistanceLower.CopyTo(distanceScratch);
+                bestDistanceUpper.CopyTo(distanceScratch[Vector128<int>.Count..]);
+                bestIndexLower.CopyTo(indexScratch);
+                bestIndexUpper.CopyTo(indexScratch[Vector128<int>.Count..]);
+                for (int lane = 0; lane < Vector128<short>.Count; lane++)
+                {
+                    indices[offset + lane] = (byte)indexScratch[lane];
+                    distortion += distanceScratch[lane];
+                }
+            }
+        }
+
+        for (; offset < firstSamples.Length; offset++)
+        {
+            int firstDifference = firstSamples[offset] - firstCentroids[0];
+            int secondDifference = secondSamples[offset] - secondCentroids[0];
+            int bestDistance = (firstDifference * firstDifference) + (secondDifference * secondDifference);
+            int bestIndex = 0;
+            for (int centroidIndex = 1; centroidIndex < firstCentroids.Length; centroidIndex++)
+            {
+                firstDifference = firstSamples[offset] - firstCentroids[centroidIndex];
+                secondDifference = secondSamples[offset] - secondCentroids[centroidIndex];
+                int distance = (firstDifference * firstDifference) + (secondDifference * secondDifference);
+                if (distance < bestDistance)
+                {
+                    bestDistance = distance;
+                    bestIndex = centroidIndex;
+                }
+            }
+
+            indices[offset] = (byte)bestIndex;
+            distortion += bestDistance;
+        }
+
+        return distortion;
+    }
+
+    /// <summary>
+    /// Refines initialized paired colors through the reference encoder's deterministic clustering sequence.
+    /// </summary>
+    /// <param name="firstSamples">The active first-plane samples.</param>
+    /// <param name="secondSamples">The active second-plane samples.</param>
+    /// <param name="firstCentroids">The initialized first-plane colors.</param>
+    /// <param name="secondCentroids">The initialized second-plane colors.</param>
+    /// <param name="indices">The palette indices belonging to the retained colors.</param>
+    /// <returns>The retained sum of squared two-plane distances.</returns>
+    public static long Cluster(
+        ReadOnlySpan<short> firstSamples,
+        ReadOnlySpan<short> secondSamples,
+        Span<short> firstCentroids,
+        Span<short> secondCentroids,
+        Span<byte> indices)
+    {
+        Span<short> alternateFirstCentroids = stackalloc short[Av1Constants.PaletteMaxSize];
+        Span<short> alternateSecondCentroids = stackalloc short[Av1Constants.PaletteMaxSize];
+        Span<byte> alternateIndices = stackalloc byte[firstSamples.Length];
+        alternateFirstCentroids = alternateFirstCentroids[..firstCentroids.Length];
+        alternateSecondCentroids = alternateSecondCentroids[..secondCentroids.Length];
+        long distortion = AssignIndices(
+            firstSamples,
+            secondSamples,
+            firstCentroids,
+            secondCentroids,
+            indices);
+
+        bool currentIsAlternate = false;
+        for (int iteration = 0; iteration < Av1PaletteKMeans.MaximumIterations; iteration++)
+        {
+            ReadOnlySpan<short> currentFirstCentroids = currentIsAlternate
+                ? alternateFirstCentroids
+                : firstCentroids;
+
+            ReadOnlySpan<short> currentSecondCentroids = currentIsAlternate
+                ? alternateSecondCentroids
+                : secondCentroids;
+
+            ReadOnlySpan<byte> currentIndices = currentIsAlternate ? alternateIndices : indices;
+            Span<short> nextFirstCentroids = currentIsAlternate ? firstCentroids : alternateFirstCentroids;
+            Span<short> nextSecondCentroids = currentIsAlternate ? secondCentroids : alternateSecondCentroids;
+            Span<byte> nextIndices = currentIsAlternate ? indices : alternateIndices;
+            CalculateCentroids(
+                firstSamples,
+                secondSamples,
+                currentIndices,
+                nextFirstCentroids,
+                nextSecondCentroids);
+
+            if (nextFirstCentroids.SequenceEqual(currentFirstCentroids) &&
+                nextSecondCentroids.SequenceEqual(currentSecondCentroids))
+            {
+                break;
+            }
+
+            long nextDistortion = AssignIndices(
+                firstSamples,
+                secondSamples,
+                nextFirstCentroids,
+                nextSecondCentroids,
+                nextIndices);
+
+            if (nextDistortion > distortion)
+            {
+                break;
+            }
+
+            distortion = nextDistortion;
+            currentIsAlternate = !currentIsAlternate;
+        }
+
+        if (currentIsAlternate)
+        {
+            alternateFirstCentroids.CopyTo(firstCentroids);
+            alternateSecondCentroids.CopyTo(secondCentroids);
+            alternateIndices.CopyTo(indices);
+        }
+
+        return distortion;
+    }
+
+    /// <summary>
+    /// Places paired initial colors at the midpoints of equal intervals spanning each plane's sample range.
+    /// </summary>
+    public static void InitializeCentroids(
+        short firstMinimum,
+        short firstMaximum,
+        short secondMinimum,
+        short secondMaximum,
+        Span<short> firstCentroids,
+        Span<short> secondCentroids)
+    {
+        int firstRange = firstMaximum - firstMinimum;
+        int secondRange = secondMaximum - secondMinimum;
+        for (int index = 0; index < firstCentroids.Length; index++)
+        {
+            firstCentroids[index] = (short)(firstMinimum + (((2 * index) + 1) * firstRange / firstCentroids.Length / 2));
+            secondCentroids[index] = (short)(secondMinimum + (((2 * index) + 1) * secondRange / secondCentroids.Length / 2));
+        }
+    }
+
+    private static void CalculateCentroids(
+        ReadOnlySpan<short> firstSamples,
+        ReadOnlySpan<short> secondSamples,
+        ReadOnlySpan<byte> indices,
+        Span<short> firstCentroids,
+        Span<short> secondCentroids)
+    {
+        Span<int> counts = stackalloc int[Av1Constants.PaletteMaxSize];
+        Span<int> firstSums = stackalloc int[Av1Constants.PaletteMaxSize];
+        Span<int> secondSums = stackalloc int[Av1Constants.PaletteMaxSize];
+        counts = counts[..firstCentroids.Length];
+        firstSums = firstSums[..firstCentroids.Length];
+        secondSums = secondSums[..secondCentroids.Length];
+        counts.Clear();
+        firstSums.Clear();
+        secondSums.Clear();
+        for (int index = 0; index < firstSamples.Length; index++)
+        {
+            int centroidIndex = indices[index];
+            counts[centroidIndex]++;
+            firstSums[centroidIndex] += firstSamples[index];
+            secondSums[centroidIndex] += secondSamples[index];
+        }
+
+        uint randomState = (uint)firstSamples[0];
+        for (int centroidIndex = 0; centroidIndex < firstCentroids.Length; centroidIndex++)
+        {
+            int count = counts[centroidIndex];
+            if (count == 0)
+            {
+                // Empty paired clusters copy both components from the same deterministically selected sample.
+                randomState = unchecked((randomState * 1103515245U) + 12345U);
+                uint random = (randomState / 65536U) % 32768U;
+                int sampleIndex = (int)(random % (uint)firstSamples.Length);
+                firstCentroids[centroidIndex] = firstSamples[sampleIndex];
+                secondCentroids[centroidIndex] = secondSamples[sampleIndex];
+            }
+            else
+            {
+                firstCentroids[centroidIndex] = (short)((firstSums[centroidIndex] + (count / 2)) / count);
+                secondCentroids[centroidIndex] = (short)((secondSums[centroidIndex] + (count / 2)) / count);
+            }
+        }
+    }
+}
diff --git a/tests/ImageSharp.Tests/Formats/Heif/Av1/Av1PaletteKMeans2DTests.cs b/tests/ImageSharp.Tests/Formats/Heif/Av1/Av1PaletteKMeans2DTests.cs
new file mode 100644
index 0000000000..185837d54b
--- /dev/null
+++ b/tests/ImageSharp.Tests/Formats/Heif/Av1/Av1PaletteKMeans2DTests.cs
@@ -0,0 +1,125 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Heif.Av1.Pipeline;
+using SixLabors.ImageSharp.Tests.TestUtilities;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Heif.Av1;
+
+/// <summary>
+/// Verifies paired AV1 palette clustering against independent scalar results at every intrinsic tier.
+/// </summary>
+[Trait("Format", "Avif")]
+public class Av1PaletteKMeans2DTests
+{
+    private const HwIntrinsics Configurations =
+        HwIntrinsics.AllowAll | HwIntrinsics.DisableAVX512F | HwIntrinsics.DisableAVX | HwIntrinsics.DisableHWIntrinsic;
+
+    [Fact]
+    public void AssignIndicesMatchesScalarAtEveryIntrinsicTier()
+        => FeatureTestRunner.RunWithHwIntrinsicsFeature(ValidateAssignment, Configurations);
+
+    [Fact]
+    public void ClusterMatchesReferenceFixture()
+    {
+        short[] firstSamples = [0, 2, 0, 100, 102, 100, 200, 202, 200];
+        short[] secondSamples = [10, 10, 12, 110, 110, 112, 210, 210, 212];
+        short[] firstCentroids = [20, 100, 180];
+        short[] secondCentroids = [30, 110, 190];
+        byte[] indices = new byte[firstSamples.Length];
+
+        long distortion = Av1PaletteKMeans2D.Cluster(
+            firstSamples,
+            secondSamples,
+            firstCentroids,
+            secondCentroids,
+            indices);
+
+        Assert.Equal([1, 101, 201], firstCentroids);
+        Assert.Equal([11, 111, 211], secondCentroids);
+        Assert.Equal([0, 0, 0, 1, 1, 1, 2, 2, 2], indices);
+        Assert.Equal(18, distortion);
+    }
+
+    [Fact]
+    public void InitializeCentroidsMatchesReferenceIntegerOrder()
+    {
+        short[] firstCentroids = new short[3];
+        short[] secondCentroids = new short[3];
+
+        Av1PaletteKMeans2D.InitializeCentroids(10, 250, 20, 260, firstCentroids, secondCentroids);
+
+        Assert.Equal([50, 130, 210], firstCentroids);
+        Assert.Equal([60, 140, 220], secondCentroids);
+    }
+
+    private static void ValidateAssignment()
+    {
+        const int SampleCount = 95;
+        short[] firstCentroids = [0, 512, 1024, 2048, 3072, 4095];
+        short[] secondCentroids = [4094, 3072, 2048, 1024, 512, 0];
+        short[] firstSamples = new short[SampleCount];
+        short[] secondSamples = new short[SampleCount];
+        for (int index = 0; index < SampleCount; index++)
+        {
+            firstSamples[index] = (short)(((index * 977) + (index * index * 17)) & 4095);
+            secondSamples[index] = (short)(((index * 619) + (index * index * 29)) & 4095);
+        }
+
+        // The first sample is equidistant from the first two colors and must retain the first palette index.
+        firstSamples[0] = 256;
+        secondSamples[0] = 3583;
+        byte[] expected = new byte[SampleCount];
+        long expectedDistortion = AssignReference(
+            firstSamples,
+            secondSamples,
+            firstCentroids,
+            secondCentroids,
+            expected);
+
+        byte[] actual = Enumerable.Repeat(byte.MaxValue, SampleCount + 7).ToArray();
+        long actualDistortion = Av1PaletteKMeans2D.AssignIndices(
+            firstSamples,
+            secondSamples,
+            firstCentroids,
+            secondCentroids,
+            actual);
+
+        Assert.Equal(expectedDistortion, actualDistortion);
+        Assert.Equal(expected, actual.AsSpan(..SampleCount).ToArray());
+        Assert.All(actual[SampleCount..], value => Assert.Equal(byte.MaxValue, value));
+    }
+
+    private static long AssignReference(
+        ReadOnlySpan<short> firstSamples,
+        ReadOnlySpan<short> secondSamples,
+        ReadOnlySpan<short> firstCentroids,
+        ReadOnlySpan<short> secondCentroids,
+        Span<byte> indices)
+    {
+        long distortion = 0;
+        for (int sampleIndex = 0; sampleIndex < firstSamples.Length; sampleIndex++)
+        {
+            int firstDifference = firstSamples[sampleIndex] - firstCentroids[0];
+            int secondDifference = secondSamples[sampleIndex] - secondCentroids[0];
+            int bestDistance = (firstDifference * firstDifference) + (secondDifference * secondDifference);
+            int bestIndex = 0;
+            for (int centroidIndex = 1; centroidIndex < firstCentroids.Length; centroidIndex++)
+            {
+                firstDifference = firstSamples[sampleIndex] - firstCentroids[centroidIndex];
+                secondDifference = secondSamples[sampleIndex] - secondCentroids[centroidIndex];
+                int distance = (firstDifference * firstDifference) + (secondDifference * secondDifference);
+                if (distance < bestDistance)
+                {
+                    bestDistance = distance;
+                    bestIndex = centroidIndex;
+                }
+            }
+
+            indices[sampleIndex] = (byte)bestIndex;
+            distortion += bestDistance;
+        }
+
+        return distortion;
+    }
+}