Add tests, field writes, update folder structure, add group border assigner, complete noise decoder, XYB encoder, reference DCT

src/ImageSharp/Formats/Jxl/Cms/JxlOpsinConstants.cs

@@ -21,4 +21,6 @@ internal static class JxlOpsinConstants
     public const float OpsinAbsorbanceBias0 = 0.0037930732552754493f;
     public const float OpsinAbsorbanceBias1 = OpsinAbsorbanceBias0;
     public const float OpsinAbsorbanceBias2 = OpsinAbsorbanceBias0;
+
+    public static ReadOnlySpan<float> OpsinAbsorbanceBias => [OpsinAbsorbanceBias0, OpsinAbsorbanceBias1, OpsinAbsorbanceBias2];
 }

src/ImageSharp/Formats/Jxl/Fields/JxlBitsCoder.cs

@@ -31,6 +31,16 @@ internal static class JxlBitsCoder
         return true;
     }
 
+    public static void Write(int bits, uint value, JxlBitWriter writer)
+    {
+        if (value >= (1u << bits))
+        {
+            throw new InvalidOperationException("Value is too large to encode in " + bits + " bits");
+        }
+
+        writer.Write(bits, value);
+    }
+
     // NOTE: BitsCoder::Read (fields.cc:418) returns a uint32_t,
     // suggesting the input bit size does not exceed 32 bits.
     public static uint Read(uint bits, JxlBitReader reader) => reader.ReadBits32(bits);

src/ImageSharp/Formats/Jxl/Fields/JxlBundle.cs

@@ -1,7 +1,12 @@
 // Copyright (c) Six Labors.
 // Licensed under the Six Labors Split License.
 
+using SixLabors.ImageSharp.Formats.Jxl.IO.FrameHeader;
+using SixLabors.ImageSharp.Formats.Jxl.IO.Metadata;
+using SixLabors.ImageSharp.Formats.Jxl.Processing;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder.AuxiliaryOutput;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Quantization;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Fields;
 
@@ -74,6 +79,23 @@ internal static class JxlBundle
         return visitor.OK;
     }
 
+    public static bool CanEncode(IJxlFields fields, ref int extensionBits, ref long totalBits)
+    {
+        JxlCanEncodeVisitor canEncodeVisitor = new();
+
+        if (!canEncodeVisitor.Visit(fields))
+        {
+            return false;
+        }
+
+        if (!canEncodeVisitor.GetSizes(ref extensionBits, ref totalBits))
+        {
+            return false;
+        }
+
+        return true;
+    }
+
     /// <summary>
     /// Tries to read the fields from a bit-reader.
     /// </summary>
@@ -86,4 +108,66 @@ internal static class JxlBundle
         _ = visitor.Visit(fields);
         return visitor.OK;
     }
+
+    public static bool Write(IJxlFields fields, JxlBitWriter writer, JxlLayerType layer, JxlAuxiliaryOutput auxOutput)
+    {
+        int extensionBits = 0;
+        long totalBits = 0;
+        if (!CanEncode(fields, ref extensionBits, ref totalBits))
+        {
+            return false;
+        }
+
+        return writer.WithMaxBits(totalBits, layer, auxOutput, () =>
+        {
+            JxlWriteVisitor visitor = new(extensionBits, writer);
+
+            if (!visitor.Visit(fields))
+            {
+                return false;
+            }
+
+            return visitor.OK();
+        });
+    }
+
+    public static bool WriteCodestreamHeaders(JxlCodecMetadata metadata, JxlBitWriter writer, JxlAuxiliaryOutput auxOut)
+    {
+        // Marker/signature
+        if (!writer.WithMaxBits(16, JxlLayerType.Header, auxOut, () =>
+        {
+            writer.Write(8, 0xFF);
+            writer.Write(8, JxlShared.CodestreamMarker);
+            return true;
+        }))
+        {
+            return false;
+        }
+
+        if (!WriteSizeHeader(metadata.Size!, writer, JxlLayerType.Header, auxOut))
+        {
+            return false;
+        }
+
+        if (!WriteImageMetadata(metadata.ImageMetadata!, writer, JxlLayerType.Header, auxOut))
+        {
+            return false;
+        }
+
+        metadata.CustomTransformData!.NonserializedXybEncoded = metadata.ImageMetadata!.XybEncoded;
+
+        return Write(metadata.CustomTransformData!, writer, JxlLayerType.Header, auxOut);
+    }
+
+    public static bool WriteFrameHeader(JxlFrameHeader frame, JxlBitWriter writer, JxlAuxiliaryOutput auxOut)
+        => Write(frame, writer, JxlLayerType.Header, auxOut);
+
+    public static bool WriteImageMetadata(JxlImageMetadata metadata, JxlBitWriter writer, JxlLayerType layer, JxlAuxiliaryOutput auxOut)
+        => Write(metadata, writer, layer, auxOut);
+
+    public static bool WriteQuantizerParameters(JxlQuantizerParameters metadata, JxlBitWriter writer, JxlLayerType layer, JxlAuxiliaryOutput auxOut)
+        => Write(metadata, writer, layer, auxOut);
+
+    public static bool WriteSizeHeader(JxlSizeHeader header, JxlBitWriter writer, JxlLayerType layer, JxlAuxiliaryOutput auxOut)
+        => Write(header, writer, layer, auxOut);
 }

src/ImageSharp/Formats/Jxl/Fields/JxlF16Coder.cs

@@ -66,4 +66,70 @@ internal static class JxlF16Coder
 
         return true;
     }
+
+    public static bool Write(float value, JxlBitWriter writer)
+    {
+        uint bits32 = BitConverter.SingleToUInt32Bits(value);
+
+        uint sign = bits32 >> 31;
+        uint biasedExp32 = (bits32 >> 23) & 0xFF;
+        uint mantissa32 = bits32 & 0x7FFFFF;
+
+        int exp = (int)biasedExp32 - 127;
+        if (exp > 15)
+        {
+            throw new InvalidOperationException("Too big to encode, CanEncode should return false");
+        }
+
+        // Tiny or zero => zero.
+        if (exp < -24)
+        {
+            writer.Write(16, 0);
+            return true;
+        }
+
+        uint biasedExp16 = 0;
+        uint mantissa16 = 0;
+
+        if (exp < -14)
+        {
+            biasedExp16 = 0;
+
+            uint subExp = unchecked((uint)(-14 - exp));
+
+            if (subExp is not (>= 1 and < 11))
+            {
+                return false;
+            }
+
+            mantissa16 = (1u << (int)(10 - subExp)) + (mantissa32 >> (int)(13 + subExp));
+        }
+        else
+        {
+            // exp = [-14, 15]
+            biasedExp16 = unchecked((uint)(exp + 15));
+
+            if (biasedExp16 is not (>= 1 and < 31))
+            {
+                return false;
+            }
+
+            mantissa16 = mantissa32 >> 13;
+        }
+
+        if (mantissa16 >= 1024)
+        {
+            return false;
+        }
+
+        uint bits16 = (sign << 15) | (biasedExp16 << 10) | mantissa16;
+
+        if (bits16 >= 0x10000)
+        {
+            return false;
+        }
+
+        writer.Write(16, bits16);
+        return true;
+    }
 }

src/ImageSharp/Formats/Jxl/Fields/JxlU32Coder.cs

@@ -73,7 +73,7 @@ internal static class JxlU32Coder
     /// </summary>
     public static bool ChooseSelector(in JxlU32Enc enc, uint value, ref uint selector, ref int totalBits)
     {
-        int bitsRequired = 32 - JxlMath.Num0BitsAboveMS1Bit(value);
+        uint bitsRequired = 32 - JxlMath.Num0BitsAboveMS1Bit(value);
 
         if (bitsRequired > 32)
         {
@@ -116,11 +116,38 @@ internal static class JxlU32Coder
 
         if (totalBits == 64)
         {
-            DebugGuard.IsTrue(false, "No matching selector");
+            throw new InvalidOperationException("No matching selector");
+        }
 
+        return true;
+    }
+
+    public static bool Write(in JxlU32Enc enc, uint value, JxlBitWriter writer)
+    {
+        uint selector = 0;
+        int totalBits = 0;
+
+        if (!ChooseSelector(in enc, value, ref selector, ref totalBits))
+        {
             return false;
         }
 
+        writer.Write(2, selector);
+
+        JxlU32Distribution d = enc.GetDistribution((int)selector);
+
+        if (!d.IsDirect)
+        {
+            uint offset = d.Offset;
+
+            if (value < offset)
+            {
+                return false;
+            }
+
+            writer.WriteBits(totalBits - 2, value - offset);
+        }
+
         return true;
     }
 }

src/ImageSharp/Formats/Jxl/Fields/JxlU64Coder.cs

@@ -99,6 +99,59 @@ internal static class JxlU64Coder
         return true;
     }
 
+    public static bool Write(ulong value, JxlBitWriter writer)
+    {
+        if (value == 0)
+        {
+            // Selector: use 0 bits, value 0
+            writer.Write(2, 0);
+        }
+        else if (value <= 16)
+        {
+            // Selector: use 4 bits, value 1..16
+            writer.Write(2, 1);
+            writer.Write(4, value - 1);
+        }
+        else if (value <= 272)
+        {
+            // Selector: use 8 bits, value 17..272
+            writer.Write(2, 2);
+            writer.Write(8, value - 17);
+        }
+        else
+        {
+            // Selector: varint, first a 12-bit group, after that per 8-bit group.
+            writer.Write(2, 3);
+            writer.Write(12, value & 4095);
+            value >>= 12;
+            int shift = 12;
+            while (value > 0 && shift < 60)
+            {
+                // Indicate varint not done
+                writer.Write(1, 1);
+                writer.Write(8, value & 255);
+                value >>= 8;
+                shift += 8;
+            }
+
+            if (value > 0)
+            {
+                // This only could happen if shift == N - 4.
+                writer.Write(1, 1);
+                writer.Write(4, value & 15);
+
+                // Implicitly closed sequence, no extra stop bit is required.
+            }
+            else
+            {
+                // Indicate end of varint
+                writer.Write(1, 0);
+            }
+        }
+
+        return true;
+    }
+
     /// <summary>
     /// Always returns 73.
     /// </summary>

src/ImageSharp/Formats/Jxl/Processing/Blending/JxlBlending.cs

@@ -4,7 +4,7 @@
 using SixLabors.ImageSharp.Formats.Jxl.IO.FrameHeader;
 using SixLabors.ImageSharp.Formats.Jxl.IO.Metadata;
 using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes;
-using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Patch;
 using SixLabors.ImageSharp.Memory;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Blending;

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlAnsCode.cs → src/ImageSharp/Formats/Jxl/Processing/Decoder/Ans/JxlAnsCode.cs

@@ -4,7 +4,7 @@
 using System.Runtime.InteropServices;
 using SixLabors.ImageSharp.Formats.Jxl.IO.Entropy;
 
-namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
 
 internal sealed class JxlAnsCode
 {

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlAnsReader.cs → src/ImageSharp/Formats/Jxl/Processing/Decoder/Ans/JxlAnsReader.cs

@@ -5,7 +5,7 @@ using System.Buffers;
 using System.Runtime.CompilerServices;
 using SixLabors.ImageSharp.Formats.Jxl.IO.Entropy;
 
-namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
 
 internal static class JxlAnsReader
 {

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlAnsSymbolReader.cs → src/ImageSharp/Formats/Jxl/Processing/Decoder/Ans/JxlAnsSymbolReader.cs

@@ -3,7 +3,7 @@
 
 using SixLabors.ImageSharp.Formats.Jxl.IO.Entropy;
 
-namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
 
 internal sealed class JxlAnsSymbolReader
 {

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlBitReader.cs

@@ -8,7 +8,7 @@ namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
 /// <summary>
 /// Represents a bitstream reader.
 /// </summary>
-internal sealed class JxlBitReader(Stream stream)
+internal class JxlBitReader(Stream stream)
 {
     private ulong buffer;
     private uint bufferRemainingBits;
@@ -154,17 +154,17 @@ internal sealed class JxlBitReader(Stream stream)
         }
     }
 
-    public uint ReadBits32(uint bits) => this.ReadBits32Core(bits, peek: false);
+    public virtual uint ReadBits32(uint bits) => this.ReadBits32Core(bits, peek: false);
 
-    public uint PeekBits32(uint bits) => this.ReadBits32Core(bits, peek: true);
+    public virtual uint PeekBits32(uint bits) => this.ReadBits32Core(bits, peek: true);
 
-    public void SkipBits32(uint bits) => _ = this.ReadBits32(bits);
+    public virtual void SkipBits32(uint bits) => _ = this.ReadBits32(bits);
 
-    public ulong ReadBits64(ulong bits) => this.ReadBits64Core((uint)bits, peek: false);
+    public virtual ulong ReadBits64(ulong bits) => this.ReadBits64Core((uint)bits, peek: false);
 
-    public ulong PeekBits64(ulong bits) => this.ReadBits64Core((uint)bits, peek: true);
+    public virtual ulong PeekBits64(ulong bits) => this.ReadBits64Core((uint)bits, peek: true);
 
-    public void SkipBits64(ulong bits) => _ = this.ReadBits64(bits);
+    public virtual void SkipBits64(ulong bits) => _ = this.ReadBits64(bits);
 
-    public bool ReadBoolean() => this.ReadBits32Core(1, peek: false) == 1;
+    public virtual bool ReadBoolean() => this.ReadBits32Core(1, peek: false) == 1;
 }

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlGroupBorderAssigner.cs

@@ -0,0 +1,330 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+using SixLabors.ImageSharp.Memory;
+
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+
+internal sealed class JxlGroupBorderAssigner
+{
+    public const int MaxToFinalize = 3;
+
+    private readonly JxlFrameDimensions frameDim;
+    private readonly Buffer2D<int> counters;
+
+    private const byte TopLeft = 0x01;
+    private const byte TopRight = 0x02;
+    private const byte BottomRight = 0x04;
+    private const byte BottomLeft = 0x08;
+
+    public JxlGroupBorderAssigner(Configuration configuration, JxlFrameDimensions frameDim)
+    {
+        this.frameDim = frameDim;
+
+        int xSize = frameDim.XSizeGroups;
+        int ySize = frameDim.YSizeGroups;
+
+        int counterCount = (xSize + 1 + 7) / 8;
+        this.counters = configuration.MemoryAllocator.Allocate2D<int>(ySize + 1, counterCount);
+
+        void Set(int x, int y, uint corners)
+        {
+            int shift = 4 * (x & 7);
+            int bits = (int)(corners << shift);
+            this.counters[y, x / 8] |= bits;
+        }
+
+        for (int x = 0; x < xSize + 1; x++)
+        {
+            Set(x, 0, TopLeft | TopRight);
+            Set(x, ySize, BottomLeft | BottomRight);
+        }
+
+        for (int y = 0; y < ySize + 1; y++)
+        {
+            Set(0, y, TopLeft | BottomLeft);
+            Set(xSize, y, TopRight | BottomRight);
+        }
+    }
+
+    public void ClearDone(int groupId)
+    {
+        void Clear(int x, int y, uint corners)
+        {
+            int shift = 4 * (x & 7);
+            int mask = ~(int)(corners << shift);
+
+            ref int counter = ref this.counters[y, x / 8];
+            _ = Interlocked.And(ref counter, mask);
+        }
+
+        (int x, int y) = Math.DivRem(groupId, this.frameDim.XSizeGroups);
+
+        Clear(x, y, BottomRight);
+        Clear(x + 1, y, BottomLeft);
+        Clear(x, y + 1, TopRight);
+        Clear(x + 1, y + 1, TopLeft);
+    }
+
+    public void GroupDone(int groupId, int padx, int pady, Span<Rectangle> rectsToFinalize, out int numToFinalize)
+    {
+        (int x, int y) = Math.DivRem(groupId, this.frameDim.XSizeGroups);
+
+        Rectangle blockRect = RectangleUtils.CreateRectangle(
+            x * this.frameDim.GroupDimension / JxlFrameDimensions.BlockDimensions,
+            y * this.frameDim.GroupDimension / JxlFrameDimensions.BlockDimensions,
+            this.frameDim.GroupDimension / JxlFrameDimensions.BlockDimensions,
+            this.frameDim.GroupDimension / JxlFrameDimensions.BlockDimensions,
+            this.frameDim.XSizeBlocks,
+            this.frameDim.YSizeBlocks);
+
+        int FetchStatus(int cx, int cy, uint bit)
+        {
+            int shift = 4 * (cx & 7);
+            int bits = (int)(bit << shift);
+
+            ref int counter = ref this.counters[cy, cx / 8];
+
+            int status = Interlocked.Or(ref counter, bits);
+            status >>= shift;
+
+            return (int)((bit | (uint)status) & 0xF);
+        }
+
+        int topLeftStatus = FetchStatus(x, y, BottomRight);
+        int topRightStatus = FetchStatus(x + 1, y, BottomLeft);
+        int bottomRightStatus = FetchStatus(x + 1, y + 1, TopLeft);
+        int bottomLeftStatus = FetchStatus(x, y + 1, TopRight);
+
+        int x1 = blockRect.X + blockRect.Width;
+        int y1 = blockRect.Y + blockRect.Height;
+
+        bool isLastGroupX = this.frameDim.XSizeGroups == x + 1;
+        bool isLastGroupY = this.frameDim.YSizeGroups == y + 1;
+
+        Span<int> xpos =
+        [
+            blockRect.X == 0
+                ? 0
+                : (blockRect.X * JxlFrameDimensions.BlockDimensions) - padx,
+
+            blockRect.X == 0
+                ? 0
+                : Math.Min(
+                    this.frameDim.XSize,
+                    (blockRect.X * JxlFrameDimensions.BlockDimensions) + padx),
+
+            isLastGroupX
+                ? this.frameDim.XSize
+                : (x1 * JxlFrameDimensions.BlockDimensions) - padx,
+
+            Math.Min(
+                this.frameDim.XSize,
+                (x1 * JxlFrameDimensions.BlockDimensions) + padx)
+        ];
+
+        Span<int> ypos =
+        [
+            blockRect.Y == 0
+                ? 0
+                : (blockRect.Y * JxlFrameDimensions.BlockDimensions) - pady,
+
+            blockRect.Y == 0
+                ? 0
+                : Math.Min(
+                     this.frameDim.YSize,
+                     (blockRect.Y * JxlFrameDimensions.BlockDimensions) + pady),
+
+            isLastGroupY
+                ? this.frameDim.YSize
+                : (y1 * JxlFrameDimensions.BlockDimensions) - pady,
+
+            Math.Min(
+                 this.frameDim.YSize,
+                 (y1 * JxlFrameDimensions.BlockDimensions) + pady)
+        ];
+
+        numToFinalize = 0;
+
+        void AppendRect(int x0, int x1, int y0, int y1, ref int numToFinalize, Span<Rectangle> rectsToFinalize, Span<int> xpos, Span<int> ypos)
+        {
+            Rectangle rect = new(
+                xpos[x0],
+                ypos[y0],
+                xpos[x1] - xpos[x0],
+                ypos[y1] - ypos[y0]);
+
+            if (rect.Width == 0 || rect.Height == 0)
+            {
+                return;
+            }
+
+            rectsToFinalize[numToFinalize++] = rect;
+        }
+
+        bool[,] availablePartsMask = new bool[3, 3];
+
+        // Center
+        availablePartsMask[1, 1] = true;
+
+        // Corners
+        if (topLeftStatus == 0xF)
+        {
+            availablePartsMask[0, 0] = true;
+        }
+
+        if (topRightStatus == 0xF)
+        {
+            availablePartsMask[2, 0] = true;
+        }
+
+        if (bottomRightStatus == 0xF)
+        {
+            availablePartsMask[2, 2] = true;
+        }
+
+        if (bottomLeftStatus == 0xF)
+        {
+            availablePartsMask[0, 2] = true;
+        }
+
+        // Other borders
+        if ((topLeftStatus & TopRight) != 0)
+        {
+            availablePartsMask[1, 0] = true;
+        }
+
+        if ((topLeftStatus & BottomLeft) != 0)
+        {
+            availablePartsMask[0, 1] = true;
+        }
+
+        if ((topRightStatus & BottomRight) != 0)
+        {
+            availablePartsMask[2, 1] = true;
+        }
+
+        if ((bottomLeftStatus & BottomRight) != 0)
+        {
+            availablePartsMask[1, 2] = true;
+        }
+
+        const int noSegment = 3;
+
+        Span<(int First, int Last)> horizontalSegments =
+        [
+            (noSegment, noSegment),
+            (noSegment, noSegment),
+            (noSegment, noSegment)
+        ];
+
+        for (int py = 0; py < 3; py++)
+        {
+            for (int px = 0; px < 3; px++)
+            {
+                if (!availablePartsMask[px, py])
+                {
+                    continue;
+                }
+
+                if (horizontalSegments[py].First == noSegment)
+                {
+                    horizontalSegments[py] = (px, horizontalSegments[py].Last);
+                }
+
+                horizontalSegments[py] = (horizontalSegments[py].First, px + 1);
+            }
+        }
+
+        if (horizontalSegments[0] == horizontalSegments[1] &&
+            horizontalSegments[0] == horizontalSegments[2])
+        {
+            AppendRect(
+                horizontalSegments[0].First,
+                horizontalSegments[0].Last,
+                0,
+                3,
+                ref numToFinalize,
+                rectsToFinalize,
+                xpos,
+                ypos);
+        }
+        else if (horizontalSegments[0] == horizontalSegments[1])
+        {
+            AppendRect(
+                horizontalSegments[0].First,
+                horizontalSegments[0].Last,
+                0,
+                2,
+                ref numToFinalize,
+                rectsToFinalize,
+                xpos,
+                ypos);
+
+            AppendRect(
+                horizontalSegments[2].First,
+                horizontalSegments[2].Last,
+                2,
+                3,
+                ref numToFinalize,
+                rectsToFinalize,
+                xpos,
+                ypos);
+        }
+        else if (horizontalSegments[1] == horizontalSegments[2])
+        {
+            AppendRect(
+                horizontalSegments[0].First,
+                horizontalSegments[0].Last,
+                0,
+                1,
+                ref numToFinalize,
+                rectsToFinalize,
+                xpos,
+                ypos);
+
+            AppendRect(
+                horizontalSegments[1].First,
+                horizontalSegments[1].Last,
+                1,
+                3,
+                ref numToFinalize,
+                rectsToFinalize,
+                xpos,
+                ypos);
+        }
+        else
+        {
+            AppendRect(
+                horizontalSegments[0].First,
+                horizontalSegments[0].Last,
+                0,
+                1,
+                ref numToFinalize,
+                rectsToFinalize,
+                xpos,
+                ypos);
+
+            AppendRect(
+                horizontalSegments[1].First,
+                horizontalSegments[1].Last,
+                1,
+                2,
+                ref numToFinalize,
+                rectsToFinalize,
+                xpos,
+                ypos);
+
+            AppendRect(
+                horizontalSegments[2].First,
+                horizontalSegments[2].Last,
+                2,
+                3,
+                ref numToFinalize,
+                rectsToFinalize,
+                xpos,
+                ypos);
+        }
+    }
+}

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlNoiseDecoder.cs

@@ -1,11 +1,16 @@
 // Copyright (c) Six Labors.
 // Licensed under the Six Labors Split License.
 
+using System.Net.NetworkInformation;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
+using SixLabors.ImageSharp.Formats.Jxl.IO.FrameHeader;
 using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Noise;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+using SixLabors.ImageSharp.Formats.Tiff.Utils;
+using static System.Net.Mime.MediaTypeNames;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
 
@@ -61,4 +66,85 @@ internal static class JxlNoiseDecoder
             }
         }
     }
+
+    public static void Random3Planes(
+        int visibleFrameIndex,
+        int nonVisibleFrameIndex,
+        int x0,
+        int y0,
+        (JxlImageF Image, Rectangle Bounds) plane0,
+        (JxlImageF Image, Rectangle Bounds) plane1,
+        (JxlImageF Image, Rectangle Bounds) plane2)
+    {
+        JxlXorShift rng = new((uint)visibleFrameIndex, (uint)nonVisibleFrameIndex, (uint)x0, (uint)y0);
+
+        GenerateRandomImage(rng, plane0.Bounds, plane0.Image);
+        GenerateRandomImage(rng, plane1.Bounds, plane1.Image);
+        GenerateRandomImage(rng, plane2.Bounds, plane2.Image);
+    }
+
+    private static void PrepareNoiseInput(
+        JxlPassesDecoderState decState,
+        JxlFrameDimensions frameDim,
+        JxlFrameHeader frameHeader,
+        int groupIndex,
+        int thread)
+    {
+        int groupDim = frameDim.GroupDimension;
+
+        int gx = groupIndex % frameDim.XSizeGroups;
+        int gy = groupIndex / frameDim.XSizeGroups;
+
+        JxlRenderPipelineInput input = decState.RenderPipeline.GetInputBuffers(groupIndex, thread);
+        int noiseCStart = 3 + frameHeader.Metadata!.ImageMetadata!.ExtraChannelCount;
+
+        InlineArray3<(JxlImageF Image, Rectangle Rect)> rects = default;
+
+        for (int iy = 0; iy < frameHeader.Upsampling; iy++)
+        {
+            for (int ix = 0; ix < frameHeader.Upsampling; ix++)
+            {
+                for (int c = 0; c < 3; c++)
+                {
+                    var r = input.GetBuffer(noiseCStart + c);
+
+                    int x1 = r.Rect.X0 + r.Rect.XSize;
+                    int y1 = r.Rect.Y0 + r.Rect.YSize;
+
+                    rects[c] = (
+                        r.Image,
+                        RectangleUtils.CreateRectangle(
+                            r.Rect.X0 + (ix * groupDim),
+                            r.Rect.Y0 + (iy * groupDim),
+                            groupDim,
+                            groupDim,
+                            x1,
+                            y1));
+                }
+
+                Random3Planes(
+                    decState.VisibleFrameIndex,
+                    decState.NonvisibleFrameIndex,
+                    (int)((gx * frameHeader.Upsampling) + ix) * groupDim,
+                    (int)((gy * frameHeader.Upsampling) + iy) * groupDim,
+                    rects[0],
+                    rects[1],
+                    rects[2]);
+            }
+        }
+    }
+
+    public static float DecodeFloatParam(float precision, JxlBitReader br)
+    {
+        int absValQuant = (int)br.ReadBits32(10);
+        return absValQuant / precision;
+    }
+
+    public static void DecodeNoise(JxlBitReader br, JxlNoiseParameters noiseParams)
+    {
+        for (int i = 0; i < noiseParams.Lookup.Length; i++)
+        {
+            noiseParams.Lookup[i] = DecodeFloatParam(JxlNoiseConstants.Precision, br);
+        }
+    }
 }

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlPassesDecoderState.cs

@@ -4,6 +4,7 @@
 using SixLabors.ImageSharp.Formats.Jxl.IO.FrameHeader;
 using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Dct;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Image;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;

src/ImageSharp/Formats/Jxl/Processing/Decoder/Modular/JxlModularFrameDecoder.cs

@@ -0,0 +1,25 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Modular;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Modular.Encoding;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Modular.Transforms;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Modular;
+
+internal sealed class JxlModularFrameDecoder
+{
+    private readonly Configuration configuration;
+    private readonly JxlModularImage fullImage;
+    private readonly List<JxlTransform> globalTransform;
+    private JxlFrameDimensions frameDimensions;
+    private bool doColor;
+    private bool haveSomething;
+    private bool useFullImage = true;
+    private bool allSameShift;
+    private JxlAnsCode code;
+    private readonly List<byte> contextMap;
+    private JxlGroupHeader groupHeader;
+}

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlPatchBlendMode.cs → src/ImageSharp/Formats/Jxl/Processing/Decoder/Patch/JxlPatchBlendMode.cs

@@ -1,7 +1,7 @@
 // Copyright (c) Six Labors.
 // Licensed under the Six Labors Split License.
 
-namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Patch;
 
 internal enum JxlPatchBlendMode : byte
 {

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlPatchBlending.cs → src/ImageSharp/Formats/Jxl/Processing/Decoder/Patch/JxlPatchBlending.cs

@@ -1,7 +1,7 @@
 // Copyright (c) Six Labors.
 // Licensed under the Six Labors Split License.
 
-namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Patch;
 
 internal struct JxlPatchBlending
 {

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlPatchDictionary.cs → src/ImageSharp/Formats/Jxl/Processing/Decoder/Patch/JxlPatchDictionary.cs

@@ -1,10 +1,11 @@
 // Copyright (c) Six Labors.
 // Licensed under the Six Labors Split License.
 
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Image;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
 
-namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Patch;
 
 internal sealed class JxlPatchDictionary
 {

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlPatchPosition.cs → src/ImageSharp/Formats/Jxl/Processing/Decoder/Patch/JxlPatchPosition.cs

@@ -1,7 +1,7 @@
 // Copyright (c) Six Labors.
 // Licensed under the Six Labors Split License.
 
-namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Patch;
 
 internal struct JxlPatchPosition
 {

src/ImageSharp/Formats/Jxl/Processing/Decoder/JxlPatchReferencePosition.cs → src/ImageSharp/Formats/Jxl/Processing/Decoder/Patch/JxlPatchReferencePosition.cs

@@ -1,7 +1,7 @@
 // Copyright (c) Six Labors.
 // Licensed under the Six Labors Split License.
 
-namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Patch;
 
 internal struct JxlPatchReferencePosition
 {

src/ImageSharp/Formats/Jxl/Processing/Encoder/JxlLinearAlgebra.cs

@@ -11,7 +11,7 @@ namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder;
 /// </summary>
 internal static class JxlLinearAlgebra
 {
-    public static void ConvertToDiagonal(Matrix2x2 a, Vector2 diag, Matrix2x2 u)
+    public static void ConvertToDiagonal(Matrix2x2 a, ref Vector2 diag, ref Matrix2x2 u)
     {
         DebugGuard.MustBeLessThan(Math.Abs(a[0][1] - a[1][0]), 1e-15, nameof(a));
 

src/ImageSharp/Formats/Jxl/Processing/Encoder/JxlXybEncoder.cs

@@ -0,0 +1,466 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Numerics;
+using System.Numerics.Tensors;
+using System.Runtime.InteropServices;
+using System.Runtime.Intrinsics;
+using SixLabors.ImageSharp.Formats.Jxl.Cms;
+using SixLabors.ImageSharp.Formats.Jxl.Cms.TransferFunctions;
+using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Image;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder;
+
+internal static class JxlXybEncoder
+{
+    private static void OpsinAbsorbance(
+        Vector<float> r,
+        Vector<float> g,
+        Vector<float> b,
+        ReadOnlySpan<float> premulAbsorb,
+        out Vector<float> mixed0,
+        out Vector<float> mixed1,
+        out Vector<float> mixed2)
+    {
+        ReadOnlySpan<float> bias = JxlOpsinConstants.OpsinAbsorbanceBias;
+
+        mixed0 = Vector.FusedMultiplyAdd(
+            Vector.Create(premulAbsorb),
+            r,
+            Vector.FusedMultiplyAdd(
+                Vector.Create(premulAbsorb[8..]),
+                g,
+                Vector.FusedMultiplyAdd(
+                    Vector.Create(premulAbsorb[16..]),
+                    b,
+                    Vector.Create(bias[0]))));
+
+        mixed1 = Vector.FusedMultiplyAdd(
+            Vector.Create(premulAbsorb[24..]),
+            r,
+            Vector.FusedMultiplyAdd(
+                Vector.Create(premulAbsorb[32..]),
+                g,
+                Vector.FusedMultiplyAdd(
+                    Vector.Create(premulAbsorb[40..]),
+                    b,
+                    Vector.Create(bias[1]))));
+
+        mixed2 = Vector.FusedMultiplyAdd(
+            Vector.Create(premulAbsorb[48..]),
+            r,
+            Vector.FusedMultiplyAdd(
+                Vector.Create(premulAbsorb[56..]),
+                g,
+                Vector.FusedMultiplyAdd(
+                    Vector.Create(premulAbsorb[64..]),
+                    b,
+                    Vector.Create(bias[2]))));
+    }
+
+    private static void StoreXYB(
+        Vector<float> r,
+        Vector<float> g,
+        Vector<float> b,
+        Span<float> valX,
+        Span<float> valY,
+        Span<float> valZ)
+    {
+        Vector<float> half = Vector.Create(0.5f);
+
+        (half * (r - g)).CopyTo(valX);
+        (half * (r + g)).CopyTo(valY);
+        b.CopyTo(valZ);
+    }
+
+    private static void LinearRgbToXyb(
+        Vector<float> r,
+        Vector<float> g,
+        Vector<float> b,
+        ReadOnlySpan<float> premulAbsorb,
+        Span<float> valX,
+        Span<float> valY,
+        Span<float> valZ)
+    {
+        OpsinAbsorbance(
+            r,
+            g,
+            b,
+            premulAbsorb,
+            out Vector<float> mixed0,
+            out Vector<float> mixed1,
+            out Vector<float> mixed2);
+
+        mixed0 = ZeroIfNegative(mixed0);
+        mixed1 = ZeroIfNegative(mixed1);
+        mixed2 = ZeroIfNegative(mixed2);
+
+        int n = Vector<float>.Count;
+
+        mixed0 = CubeRootAndAdd(
+            mixed0,
+            Vector.Create(premulAbsorb[(9 * n)..]));
+
+        mixed1 = CubeRootAndAdd(
+            mixed1,
+            Vector.Create(premulAbsorb[(10 * n)..]));
+
+        mixed2 = CubeRootAndAdd(
+            mixed2,
+            Vector.Create(premulAbsorb[(11 * n)..]));
+
+        StoreXYB(mixed0, mixed1, mixed2, valX, valY, valZ);
+    }
+
+    public static void LinearRgbRowToXyb(
+        Span<float> row0,
+        Span<float> row1,
+        Span<float> row2,
+        ReadOnlySpan<float> premulAbsorb,
+        int xSize)
+    {
+        int n = Vector<float>.Count;
+
+        for (int x = 0; x < xSize; x += n)
+        {
+            Vector<float> r = Vector.Create<float>(row0[x..]);
+            Vector<float> g = Vector.Create<float>(row1[x..]);
+            Vector<float> b = Vector.Create<float>(row2[x..]);
+
+            LinearRgbToXyb(
+                r,
+                g,
+                b,
+                premulAbsorb,
+                row0[x..],
+                row1[x..],
+                row2[x..]);
+        }
+    }
+
+    private static Vector<float> CubeRootAndAdd(Vector<float> v, Vector<float> add)
+    {
+        // HACK: Vector<T> doesn't support cuberoot
+        Span<float> span = stackalloc float[Vector<float>.Count];
+        v.StoreUnsafe(ref MemoryMarshal.GetReference(span));
+        TensorPrimitives.Cbrt(span, span);
+        return Vector.LoadUnsafe(ref MemoryMarshal.GetReference(span)) + add;
+    }
+
+    private static Vector<float> ZeroIfNegative(Vector<float> v) =>
+        Vector.ConditionalSelect(
+            Vector.LessThan(v, Vector<float>.Zero),
+            Vector<float>.Zero,
+            v);
+
+    private static Vector<float> LinearFromSRgb(Vector<float> encoded)
+        => JxlSRgbTransferFunction.DisplayFromEncoded(encoded);
+
+    private static void LinearSrgbToXyb(
+        Configuration configuration,
+        Memory<float> premulAbsorb,
+        JxlImage3F image)
+    {
+        int xSize = image.Width;
+
+        void ProcessRow(int y, ReadOnlySpan<float> premulAbsorb)
+        {
+            Span<float> row0 = image.PlaneRow(0, y);
+            Span<float> row1 = image.PlaneRow(1, y);
+            Span<float> row2 = image.PlaneRow(2, y);
+
+            int lanes = Vector<float>.Count;
+
+            for (int x = 0; x < xSize; x += lanes)
+            {
+                Vector<float> r = new(row0[x..]);
+                Vector<float> g = new(row1[x..]);
+                Vector<float> b = new(row2[x..]);
+
+                LinearRgbToXyb(
+                    r,
+                    g,
+                    b,
+                    premulAbsorb,
+                    row0[x..],
+                    row1[x..],
+                    row2[x..]);
+            }
+        }
+
+        _ = Parallel.For(0, image.Height, configuration.GetParallelOptions(), x =>
+        {
+            ProcessRow(x, premulAbsorb.Span);
+        });
+    }
+
+    private static void SRgbToXyb(
+        Configuration configuration,
+        Memory<float> premulAbsorb,
+        JxlImage3F image)
+    {
+        int xSize = image.Width;
+
+        void ProcessRow(int y, ReadOnlySpan<float> premulAbsorb)
+        {
+            Span<float> row0 = image.PlaneRow(0, y);
+            Span<float> row1 = image.PlaneRow(1, y);
+            Span<float> row2 = image.PlaneRow(2, y);
+
+            int lanes = Vector<float>.Count;
+
+            for (int x = 0; x < xSize; x += lanes)
+            {
+                Vector<float> r = LinearFromSRgb(new Vector<float>(row0[x..]));
+                Vector<float> g = LinearFromSRgb(new Vector<float>(row1[x..]));
+                Vector<float> b = LinearFromSRgb(new Vector<float>(row2[x..]));
+
+                LinearRgbToXyb(
+                    r,
+                    g,
+                    b,
+                    premulAbsorb,
+                    row0[x..],
+                    row1[x..],
+                    row2[x..]);
+            }
+        }
+
+        _ = Parallel.For(0, image.Height, configuration.GetParallelOptions(), x =>
+        {
+            ProcessRow(x, premulAbsorb.Span);
+        });
+    }
+
+    private static void SrgbToXybAndLinear(
+        Configuration configuration,
+        Memory<float> premulAbsorb,
+        JxlImage3F image,
+        JxlImage3F linear)
+    {
+        int xSize = image.Width;
+
+        void ProcessRow(int y, ReadOnlySpan<float> premulAbsorb)
+        {
+            Span<float> rowImage0 = image.PlaneRow(0, y);
+            Span<float> rowImage1 = image.PlaneRow(1, y);
+            Span<float> rowImage2 = image.PlaneRow(2, y);
+
+            Span<float> rowLinear0 = linear.PlaneRow(0, y);
+            Span<float> rowLinear1 = linear.PlaneRow(1, y);
+            Span<float> rowLinear2 = linear.PlaneRow(2, y);
+
+            int lanes = Vector<float>.Count;
+
+            for (int x = 0; x < xSize; x += lanes)
+            {
+                Vector<float> r = LinearFromSRgb(new Vector<float>(rowImage0[x..]));
+                Vector<float> g = LinearFromSRgb(new Vector<float>(rowImage1[x..]));
+                Vector<float> b = LinearFromSRgb(new Vector<float>(rowImage2[x..]));
+
+                r.CopyTo(rowLinear0[x..]);
+                g.CopyTo(rowLinear1[x..]);
+                b.CopyTo(rowLinear2[x..]);
+
+                LinearRgbToXyb(
+                    r,
+                    g,
+                    b,
+                    premulAbsorb,
+                    rowImage0[x..],
+                    rowImage1[x..],
+                    rowImage2[x..]);
+            }
+        }
+
+        _ = Parallel.For(0, image.Height, configuration.GetParallelOptions(), x =>
+        {
+            ProcessRow(x, premulAbsorb.Span);
+        });
+    }
+
+    private static void ComputePremulAbsorb(float intensityTarget, Span<float> premulAbsorb)
+    {
+        float mul = intensityTarget / 255.0f;
+        int lanes = Vector<float>.Count;
+
+        for (int j = 0; j < 3; j++)
+        {
+            for (int i = 0; i < 3; i++)
+            {
+                float absorb = JxlOpsinAbsorbanceMatrix[j][i] * mul;
+
+                premulAbsorb
+                    .Slice(((j * 3) + i) * lanes, lanes)
+                    .Fill(absorb);
+            }
+        }
+
+        for (int i = 0; i < 3; i++)
+        {
+            float negBiasCbrt = -MathF.Cbrt(JxlOpsinAbsorbanceBias[i]);
+
+            premulAbsorb
+                .Slice((9 + i) * lanes, lanes)
+                .Fill(negBiasCbrt);
+        }
+    }
+
+    public static bool ToXyb(
+        Configuration configuration,
+        JxlColorEncoding currentEncoding,
+        float intensityTarget,
+        JxlImageF? black,
+        JxlImage3F image,
+        JxlCmsInterface cms,
+        JxlImage3F? linear)
+    {
+        if (black is not null)
+        {
+            if (image.GetRectangle() != black.GetRectangle())
+            {
+                return false;
+            }
+        }
+
+        if (linear is not null)
+        {
+            if (image.GetRectangle() != linear.GetRectangle())
+            {
+                return false;
+            }
+        }
+
+        int lanes = Vector<float>.Count;
+
+        Memory<float> premulAbsorb = new float[lanes * 12];
+        ComputePremulAbsorb(intensityTarget, premulAbsorb.Span);
+
+        bool wantLinear = linear is not null;
+
+        JxlColorEncoding linearSrgb = JxlColorEncoding.LinearSRgb(currentEncoding.IsGray);
+
+        if (linearSrgb.SameColorEncoding(currentEncoding))
+        {
+            if (wantLinear)
+            {
+                JxlImageOperations.CopyImageTo(image, linear!);
+            }
+
+            LinearSrgbToXyb(configuration, premulAbsorb, image);
+            return true;
+        }
+
+        if (currentEncoding.IsSRgb)
+        {
+            if (wantLinear)
+            {
+                SrgbToXybAndLinear(configuration, premulAbsorb, image, linear!);
+                return true;
+            }
+
+            SRgbToXyb(configuration, premulAbsorb, image);
+            return true;
+        }
+
+        // TODO
+        ApplyColorTransform(
+            configuration,
+            currentEncoding,
+            intensityTarget,
+            image,
+            black,
+            image.GetRectangle(),
+            linearSrgb,
+            cms,
+            wantLinear ? linear! : image);
+
+        if (wantLinear)
+        {
+            JxlImageOperations.CopyImageTo(linear!, image);
+        }
+
+        LinearSrgbToXyb(configuration, premulAbsorb, image);
+        return true;
+    }
+
+    public static void RgbToYcbcr(
+        Configuration configuration,
+        JxlImageF rPlane,
+        JxlImageF gPlane,
+        JxlImageF bPlane,
+        JxlImageF yPlane,
+        JxlImageF cbPlane,
+        JxlImageF crPlane)
+    {
+        int xSize = rPlane.Width;
+        int ySize = rPlane.Height;
+
+        if (xSize == 0 || ySize == 0)
+        {
+            return;
+        }
+
+        Vector<float> k128 = new(128.0f / 255);
+        Vector<float> kR = new(0.299f);
+        Vector<float> kG = new(0.587f);
+        Vector<float> kB = new(0.114f);
+        Vector<float> kAmpR = new(0.701f);
+        Vector<float> kAmpB = new(0.886f);
+
+        Vector<float> kDiffR = kAmpR + kR;
+        Vector<float> kDiffB = kAmpB + kB;
+
+        Vector<float> kNormR = Vector<float>.One / (kAmpR + kG + kB);
+        Vector<float> kNormB = Vector<float>.One / (kR + kG + kAmpB);
+
+        int step = Vector<float>.Count;
+
+        int groupArea = JxlFrameDimensions.GroupDimensions * JxlFrameDimensions.GroupDimensions;
+        int linesPerGroup = JxlMath.DivCeil(groupArea, xSize);
+        int numStripes = JxlMath.DivCeil(ySize, linesPerGroup);
+
+        void Transform(int index)
+        {
+            int y0 = index * linesPerGroup;
+            int y1 = Math.Min(y0 + linesPerGroup, ySize);
+
+            for (int y = y0; y < y1; y++)
+            {
+                ReadOnlySpan<float> rRow = rPlane.GetRow(y);
+                ReadOnlySpan<float> gRow = gPlane.GetRow(y);
+                ReadOnlySpan<float> bRow = bPlane.GetRow(y);
+
+                Span<float> yRow = yPlane.GetRow(y);
+                Span<float> cbRow = cbPlane.GetRow(y);
+                Span<float> crRow = crPlane.GetRow(y);
+
+                for (int x = 0; x < xSize; x += step)
+                {
+                    Vector<float> r = new(rRow[x..]);
+                    Vector<float> g = new(gRow[x..]);
+                    Vector<float> b = new(bRow[x..]);
+
+                    Vector<float> rBase = r * kR;
+                    Vector<float> rDiff = r * kDiffR;
+                    Vector<float> gBase = g * kG;
+                    Vector<float> bBase = b * kB;
+                    Vector<float> bDiff = b * kDiffB;
+
+                    Vector<float> yBase = rBase + gBase + bBase;
+                    Vector<float> yVec = yBase - k128;
+                    Vector<float> cbVec = (bDiff - yBase) * kNormB;
+                    Vector<float> crVec = (rDiff - yBase) * kNormR;
+
+                    yVec.CopyTo(yRow[x..]);
+                    cbVec.CopyTo(cbRow[x..]);
+                    crVec.CopyTo(crRow[x..]);
+                }
+            }
+        }
+
+        _ = Parallel.For(0, numStripes, configuration.GetParallelOptions(), Transform);
+    }
+}

src/ImageSharp/Formats/Jxl/Processing/Image/JxlImageFeatures.cs

@@ -1,7 +1,7 @@
 // Copyright (c) Six Labors.
 // Licensed under the Six Labors Split License.
 
-using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Patch;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Noise;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Image;

src/ImageSharp/Formats/Jxl/Processing/Image/JxlImageOperations.cs

@@ -364,6 +364,20 @@ internal static class JxlImageOperations
         }
     }
 
+    /// <summary>
+    /// Fills every value in every plane of the image with <paramref name="value"/>.
+    /// </summary>
+    /// <typeparam name="T">The type of the image.</typeparam>
+    /// <param name="value">The value to fill every plane with.</param>
+    /// <param name="image">The image to fill.</param>
+    public static void FillImage<T>(T value, JxlImage3<T> image)
+        where T : unmanaged
+    {
+        FillImage(value, image.Plane(0));
+        FillImage(value, image.Plane(1));
+        FillImage(value, image.Plane(2));
+    }
+
     /// <summary>
     /// Sets every value in the plane to 0. See also <seealso cref="JxlPlane{T}.Clear()"/> .
     /// </summary>

src/ImageSharp/Formats/Jxl/Processing/JxlCoefficientOrder.cs

@@ -4,6 +4,7 @@
 using System.Runtime.CompilerServices;
 using SixLabors.ImageSharp.Formats.Jxl.IO.Entropy;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing;

src/ImageSharp/Formats/Jxl/Processing/JxlConvolve.cs

@@ -263,10 +263,10 @@ internal static class JxlConvolve
     }
 
     public static void SlowSymmetric3(
-        JxlImageF input,
+        JxlPlane<float> input,
         Rectangle rect,
         JxlWeightsSymmetric3 weights,
-        JxlImageF output)
+        JxlPlane<float> output)
     {
         int width = rect.Width;
         int height = rect.Height;

src/ImageSharp/Formats/Jxl/Processing/Modular/Encoding/JxlMaDecoder.cs

@@ -2,6 +2,7 @@
 // Licensed under the Six Labors Split License.
 
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Modular.Encoding.ContextPrediction;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
 using Tree = System.Collections.Generic.List<SixLabors.ImageSharp.Formats.Jxl.Processing.Modular.Encoding.JxlPropertyDecisionNode>;

src/ImageSharp/Formats/Jxl/Processing/Modular/Encoding/JxlModularEncoding.cs

@@ -5,6 +5,7 @@ using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
 using SixLabors.ImageSharp.Formats.Jxl.Fields;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Modular.Encoding.ContextPrediction;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Modular.Transforms;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;

src/ImageSharp/Formats/Jxl/Processing/Primitives/JxlXorShift.cs

@@ -23,7 +23,7 @@ internal sealed class JxlXorShift
         }
     }
 
-    public void XorShift128Plus(uint seed1, uint seed2, uint seed3, uint seed4)
+    public JxlXorShift(uint seed1, uint seed2, uint seed3, uint seed4)
     {
         this.s0[0] = SplitMix64((((ulong)seed1 << 32) + seed2) + 0x9E3779B97F4A7C15uL);
         this.s1[0] = SplitMix64((((ulong)seed3 << 32) + seed4) + 0x9E3779B97F4A7C15uL);

src/ImageSharp/Formats/Jxl/Processing/Primitives/RectangleUtils.cs

@@ -5,13 +5,13 @@ namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
 
 internal static class RectangleUtils
 {
-    public static int X0(in Rectangle rect) => rect.X;
+    public static int X0(this in Rectangle rect) => rect.X;
 
-    public static int Y0(in Rectangle rect) => rect.Y;
+    public static int Y0(this in Rectangle rect) => rect.Y;
 
-    public static int X1(in Rectangle rect) => rect.X + rect.Width;
+    public static int X1(this in Rectangle rect) => rect.X + rect.Width;
 
-    public static int Y1(in Rectangle rect) => rect.Y + rect.Height;
+    public static int Y1(this in Rectangle rect) => rect.Y + rect.Height;
 
     public static Rectangle Extend(Rectangle curr, int border, Rectangle parent)
     {
@@ -22,4 +22,18 @@ internal static class RectangleUtils
 
         return new(newX0, newY0, newX1 - newX0, newY1 - newY0);
     }
+
+    public static bool IsInside(this Rectangle a, Rectangle b) =>
+        a.X0() >= b.X0() &&
+        a.X1() <= b.X1() &&
+        a.Y0() >= b.Y0() &&
+        a.Y1() <= b.Y1();
+
+    public static Rectangle CreateRectangle(int xbegin, int ybegin, int xSizeMax, int ySizeMax, int xEnd, int yEnd)
+        => new(xbegin, ybegin, ClampedSize(xbegin, xSizeMax, xEnd), ClampedSize(ybegin, ySizeMax, yEnd));
+
+    private static int ClampedSize(int begin, int sizeMax, int end)
+        => begin + sizeMax <= end
+            ? sizeMax
+            : (end > begin ? end - begin : 0);
 }

src/ImageSharp/Formats/Jxl/Processing/RenderPipeline/PatchDictionaryStage.cs

@@ -2,7 +2,7 @@
 // Licensed under the Six Labors Split License.
 
 using SixLabors.ImageSharp.Formats.Jxl.IO.Metadata;
-using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Patch;
 using SixLabors.ImageSharp.Memory;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.RenderPipeline;

src/ImageSharp/Formats/Jxl/Processing/Splines/JxlQuantizedSpline.cs

@@ -5,6 +5,7 @@ using System.Buffers;
 using System.Runtime.CompilerServices;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Dct;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Splines;

src/ImageSharp/Formats/Jxl/Processing/Splines/JxlSplines.cs

@@ -5,6 +5,7 @@ using System.Buffers;
 using System.Runtime.InteropServices;
 using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Splines;

tests/ImageSharp.Tests/Formats/Jxl/Processing/AlphaBlendingTests.cs

@@ -0,0 +1,332 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.ColorProfiles;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Blending;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl.Processing;
+
+public class AlphaBlendingTests
+{
+    private static float[] RgbToArray(Rgb rgb) => [rgb.R, rgb.G, rgb.B];
+
+    private static readonly float[] Expected = [77.2f, 83.0f, 90.6f];
+
+    [Fact]
+    public void BlendingWithNonPremultiplied()
+    {
+        Rgb bgRgb = new(100, 110, 120);
+        float bgA = 180f / 255;
+
+        Rgb fgRgb = new(25, 21, 23);
+        float fgA = 15420f / 65535;
+        float fgA2 = 2.0f;
+
+        float[] outR = [0];
+        float[] outG = [0];
+        float[] outB = [0];
+        float[] outAlpha = [0];
+
+        JxlAlphaBlendingInputLayer inputBg = new()
+        {
+            R = [bgRgb.R],
+            G = [bgRgb.G],
+            B = [bgRgb.B],
+            A = [bgA],
+        };
+
+        JxlAlphaBlendingInputLayer inputFg = new()
+        {
+            R = [fgRgb.R],
+            G = [fgRgb.G],
+            B = [fgRgb.B],
+            A = [fgA],
+        };
+
+        JxlAlphaBlendingOutput output = new()
+        {
+            R = outR,
+            G = outG,
+            B = outB,
+            A = outAlpha,
+        };
+
+        JxlAlphaHelper.PerformAlphaBlending(
+            inputBg,
+            inputFg,
+            output,
+            1,
+            alphaIsPremultiplied: false,
+            clamp: false);
+
+        Assert.Equal(
+            Expected,
+            [outR[0], outG[0], outB[0]],
+            new ApproximateFloatComparer(0.05f));
+
+        Assert.True(MathF.Abs(outAlpha[0] - (3174f / 4095)) <= 1e-5f);
+
+        inputFg = new()
+        {
+            R = [fgRgb.R],
+            G = [fgRgb.G],
+            B = [fgRgb.B],
+            A = [fgA2],
+        };
+
+        JxlAlphaHelper.PerformAlphaBlending(
+            inputBg,
+            inputFg,
+            output,
+            1,
+            alphaIsPremultiplied: false,
+            clamp: true);
+
+        Assert.Equal(
+            RgbToArray(fgRgb),
+            [outR[0], outG[0], outB[0]],
+            new ApproximateFloatComparer(0.05f));
+
+        Assert.True(MathF.Abs(outAlpha[0] - 1.0f) <= 1e-5f);
+    }
+
+    [Fact]
+    public void BlendingWithPremultiplied()
+    {
+        Rgb bgRgb = new(100, 110, 120);
+        float bgA = 180f / 255;
+
+        Rgb fgRgb = new(25, 21, 23);
+        float fgA = 15420f / 65535;
+        float fgA2 = 2.0f;
+
+        float[] outR = [0];
+        float[] outG = [0];
+        float[] outB = [0];
+        float[] outAlpha = [0];
+
+        JxlAlphaBlendingInputLayer inputBg = new()
+        {
+            R = [bgRgb.R],
+            G = [bgRgb.G],
+            B = [bgRgb.B],
+            A = [bgA],
+        };
+
+        JxlAlphaBlendingInputLayer inputFg = new()
+        {
+            R = [fgRgb.R],
+            G = [fgRgb.G],
+            B = [fgRgb.B],
+            A = [fgA],
+        };
+
+        JxlAlphaBlendingOutput output = new()
+        {
+            R = outR,
+            G = outG,
+            B = outB,
+            A = outAlpha,
+        };
+
+        JxlAlphaHelper.PerformAlphaBlending(
+            inputBg,
+            inputFg,
+            output,
+            1,
+            alphaIsPremultiplied: true,
+            clamp: false);
+
+        Assert.Equal(
+            new float[] { 101.5f, 105.1f, 114.8f },
+            [outR[0], outG[0], outB[0]],
+            new ApproximateFloatComparer(0.05f));
+
+        Assert.True(MathF.Abs(outAlpha[0] - (3174f / 4095)) <= 1e-5f);
+
+        inputFg = new()
+        {
+            R = [fgRgb.R],
+            G = [fgRgb.G],
+            B = [fgRgb.B],
+            A = [fgA2],
+        };
+
+        JxlAlphaHelper.PerformAlphaBlending(
+            inputBg,
+            inputFg,
+            output,
+            1,
+            alphaIsPremultiplied: true,
+            clamp: true);
+
+        Assert.Equal(
+            RgbToArray(fgRgb),
+            [outR[0], outG[0], outB[0]],
+            new ApproximateFloatComparer(0.05f));
+
+        Assert.True(MathF.Abs(outAlpha[0] - 1.0f) <= 1e-5f);
+    }
+
+    [Fact]
+    public void Mul()
+    {
+        Span<float> bg = [100];
+        Span<float> fg = [25];
+        Span<float> output = [0];
+
+        JxlAlphaHelper.PerformMultiplyBlending(
+            bg,
+            fg,
+            output,
+            1,
+            clamp: false);
+
+        Assert.True(MathF.Abs(output[0] - (fg[0] * bg[0])) <= 0.05f);
+
+        JxlAlphaHelper.PerformMultiplyBlending(
+            bg,
+            fg,
+            output,
+            1,
+            clamp: true);
+
+        Assert.True(MathF.Abs(output[0] - bg[0]) <= 0.05f);
+    }
+
+    [Fact]
+    public void PremultiplyAndUnpremultiply()
+    {
+        float[] alpha =
+        {
+            0f,
+            63f / 255,
+            127f / 255,
+            1f
+        };
+
+        float[] r = { 120, 130, 140, 150 };
+        float[] g = { 124, 134, 144, 154 };
+        float[] b = { 127, 137, 147, 157 };
+
+        JxlAlphaHelper.PremultiplyAlpha(r, g, b, alpha, alpha.Length);
+
+        Assert.Equal(
+            new[]
+            {
+                0f,
+                130 * 63f / 255,
+                140 * 127f / 255,
+                150f
+            },
+            r,
+            new ApproximateFloatComparer(1e-5f));
+
+        Assert.Equal(
+            new[]
+            {
+                0f,
+                134 * 63f / 255,
+                144 * 127f / 255,
+                154f
+            },
+            g,
+            new ApproximateFloatComparer(1e-5f));
+
+        Assert.Equal(
+            new[]
+            {
+                0f,
+                137 * 63f / 255,
+                147 * 127f / 255,
+                157f
+            },
+            b,
+            new ApproximateFloatComparer(1e-5f));
+
+        JxlAlphaHelper.UnpremultiplyAlpha(r, g, b, alpha, alpha.Length);
+
+        Assert.Equal(
+            new[] { 120f, 130f, 140f, 150f },
+            r,
+            new ApproximateFloatComparer(1e-4f));
+
+        Assert.Equal(
+            new[] { 124f, 134f, 144f, 154f },
+            g,
+            new ApproximateFloatComparer(1e-4f));
+
+        Assert.Equal(
+            new[] { 127f, 137f, 147f, 157f },
+            b,
+            new ApproximateFloatComparer(1e-4f));
+    }
+
+    [Fact]
+    public void UnpremultiplyAndPremultiply()
+    {
+        float[] alpha =
+        {
+            0f,
+            63f / 255,
+            127f / 255,
+            1f
+        };
+
+        float[] r = { 50, 60, 70, 80 };
+        float[] g = { 54, 64, 74, 84 };
+        float[] b = { 57, 67, 77, 87 };
+
+        JxlAlphaHelper.UnpremultiplyAlpha(r, g, b, alpha, alpha.Length);
+
+        Assert.Equal(
+            new[]
+            {
+                50f * (1 << 26),
+                60 * 255f / 63,
+                70 * 255f / 127,
+                80f
+            },
+            r,
+            new ApproximateFloatComparer(1e-4f));
+
+        Assert.Equal(
+            new[]
+            {
+                54f * (1 << 26),
+                64 * 255f / 63,
+                74 * 255f / 127,
+                84f
+            },
+            g,
+            new ApproximateFloatComparer(1e-4f));
+
+        Assert.Equal(
+            new[]
+            {
+                57f * (1 << 26),
+                67 * 255f / 63,
+                77 * 255f / 127,
+                87f
+            },
+            b,
+            new ApproximateFloatComparer(1e-4f));
+
+        JxlAlphaHelper.PremultiplyAlpha(r, g, b, alpha, alpha.Length);
+
+        Assert.Equal(
+            new[] { 50f, 60f, 70f, 80f },
+            r,
+            new ApproximateFloatComparer(1e-4f));
+
+        Assert.Equal(
+            new[] { 54f, 64f, 74f, 84f },
+            g,
+            new ApproximateFloatComparer(1e-4f));
+
+        Assert.Equal(
+            new[] { 57f, 67f, 77f, 87f },
+            b,
+            new ApproximateFloatComparer(1e-4f));
+    }
+}

tests/ImageSharp.Tests/Formats/Jxl/Processing/Encoder/LinearAlgebraTests.cs

@@ -0,0 +1,169 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder;
+using Matrix2x2 = System.Runtime.CompilerServices.InlineArray2<System.Runtime.CompilerServices.InlineArray2<double>>;
+using Vector2 = System.Runtime.CompilerServices.InlineArray2<double>;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl.Processing.Encoder;
+
+public class LinearAlgebraTests
+{
+    private static Vector2 CreateVector2(double a, double b)
+    {
+        Vector2 v2 = default;
+
+        v2[0] = a;
+        v2[1] = b;
+
+        return v2;
+    }
+
+    private static Matrix2x2 Diagonal(Vector2 d)
+    {
+        Matrix2x2 result = default;
+
+        result[0][0] = d[0];
+        result[0][1] = 0;
+        result[1][0] = 0;
+        result[1][1] = d[1];
+
+        return result;
+    }
+
+    private static Matrix2x2 Identity() => Diagonal(CreateVector2(1.0f, 1.0f));
+
+    private static Matrix2x2 MatrixMultiply(Matrix2x2 a, Matrix2x2 b)
+    {
+        Matrix2x2 result = default;
+
+        for (int y = 0; y < 2; y++)
+        {
+            for (int x = 0; x < 2; x++)
+            {
+                result[y][x] = (a[0][x] * b[y][0]) + (a[1][x] * b[y][1]);
+            }
+        }
+
+        return result;
+    }
+
+    private static Matrix2x2 Transpose(Matrix2x2 a)
+    {
+        Matrix2x2 result = default;
+
+        result[0] = CreateVector2(a[0][0], a[1][0]);
+        result[1] = CreateVector2(a[0][1], a[1][1]);
+
+        return result;
+    }
+
+    private static Matrix2x2 RandomSymmetricMatrix(
+        Rng rng,
+        float vmin,
+        float vmax)
+    {
+        Matrix2x2 result = default;
+
+        result[0][0] = rng.UniformF(vmin, vmax);
+        result[0][1] = rng.UniformF(vmin, vmax);
+        result[1][0] = result[0][1];
+        result[1][1] = rng.UniformF(vmin, vmax);
+
+        return result;
+    }
+
+    private static void VerifyMatrixEqual(Matrix2x2 a, Matrix2x2 b, float epsilon)
+    {
+        TolerantMath comparer = new(epsilon);
+
+        for (int y = 0; y < 2; y++)
+        {
+            for (int x = 0; x < 2; x++)
+            {
+                Assert.True(
+                    comparer.AreEqual(a[y][x], b[y][x]),
+                    $"Matrices differ at [{y}][{x}]: {a[y][x]} != {b[y][x]}");
+            }
+        }
+    }
+
+    private static void VerifyOrthogonal(Matrix2x2 a, float epsilon) => VerifyMatrixEqual(
+        Identity(),
+        MatrixMultiply(Transpose(a), a),
+        epsilon);
+
+    [Fact]
+    public void ConvertToDiagonal()
+    {
+        {
+            Matrix2x2 i = Identity();
+            Matrix2x2 u = default;
+            Vector2 d = default;
+
+            JxlLinearAlgebra.ConvertToDiagonal(i, ref d, ref u);
+
+            VerifyMatrixEqual(i, u, 1e-15f);
+
+            for (int k = 0; k < 2; k++)
+            {
+                Assert.True(
+                    Math.Abs(d[k] - 1.0f) <= 1e-15f,
+                    $"d[{k}] = {d[k]}");
+            }
+        }
+
+        {
+            Matrix2x2 a = Identity();
+            a[0][1] = 2.0f;
+            a[1][0] = 2.0f;
+
+            Matrix2x2 u = default;
+            Vector2 d = default;
+
+            JxlLinearAlgebra.ConvertToDiagonal(a, ref d, ref u);
+
+            VerifyOrthogonal(u, 1e-12f);
+            VerifyMatrixEqual(
+                a,
+                MatrixMultiply(u, MatrixMultiply(Diagonal(d), Transpose(u))),
+                1e-12f);
+        }
+
+        {
+            Matrix2x2 a = default;
+            a[0] = CreateVector2(0.000208649f, 1.13687e-12f);
+            a[1] = CreateVector2(1.13687e-12f, 0.000208649f);
+
+            Matrix2x2 u = default;
+            Vector2 d = default;
+
+            JxlLinearAlgebra.ConvertToDiagonal(a, ref d, ref u);
+
+            VerifyOrthogonal(u, 1e-12f);
+            VerifyMatrixEqual(
+                a,
+                MatrixMultiply(u, MatrixMultiply(Diagonal(d), Transpose(u))),
+                1e-11f);
+        }
+
+        {
+            Rng rng = new(0);
+
+            for (int i = 0; i < 1_000_000; i++)
+            {
+                Matrix2x2 a = RandomSymmetricMatrix(rng, -1.0f, 1.0f);
+                Matrix2x2 u = default;
+                Vector2 d = default;
+
+                JxlLinearAlgebra.ConvertToDiagonal(a, ref d, ref u);
+
+                VerifyOrthogonal(u, 1e-12f);
+                VerifyMatrixEqual(
+                    a,
+                    MatrixMultiply(u, MatrixMultiply(Diagonal(d), Transpose(u))),
+                    5e-10f);
+            }
+        }
+    }
+}

tests/ImageSharp.Tests/Formats/Jxl/Processing/EntropyCoderTests.cs

@@ -0,0 +1,105 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Jxl.IO.Entropy;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder.Ans;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl.Processing;
+
+public class EntropyCoderTests
+{
+    [Fact]
+    public void PackUnpack()
+    {
+        for (int i = -31; i < 32; i++)
+        {
+            uint packed = JxlPackSigned.PackUnsigned(i);
+
+            Assert.True(packed < 63u);
+
+            int unpacked = JxlPackSigned.UnpackSigned(packed);
+
+            Assert.Equal(i, unpacked);
+        }
+    }
+
+    private sealed class MockBitReader : JxlBitReader
+    {
+        public ulong NBits { get; set; }
+
+        public ulong Bits { get; set; }
+
+        public MockBitReader(Stream stream)
+            : base(stream)
+        {
+        }
+
+        public override uint PeekBits32(uint n)
+        {
+            Assert.Equal(this.NBits, n);
+            return (uint)this.Bits;
+        }
+
+        public override uint ReadBits32(uint n)
+        {
+            Assert.Equal(this.NBits, n);
+            return (uint)this.Bits;
+        }
+    }
+
+    private static void HybridUintRoundtrip(
+        JxlAnsHybridUIntConfiguration config,
+        uint limit = 1u << 24)
+    {
+        Rng rng = new(0);
+
+        const int numIntegers = 1 << 20;
+
+        uint[] integers = new uint[numIntegers];
+        uint[] token = new uint[numIntegers];
+        uint[] nbits = new uint[numIntegers];
+        uint[] bits = new uint[numIntegers];
+
+        for (int i = 0; i < numIntegers; i++)
+        {
+            integers[i] = (uint)rng.UniformU(0, limit + 1);
+            config.Encode(
+                integers[i],
+                out token[i],
+                out nbits[i],
+                out bits[i]);
+        }
+
+        for (int i = 0; i < numIntegers; i++)
+        {
+            MockBitReader br = new(Stream.Null)
+            {
+                NBits = nbits[i],
+                Bits = bits[i]
+            };
+
+            Assert.Equal(
+                integers[i],
+                JxlAnsSymbolReader.ReadHybridUintConfig(
+                    config,
+                    token[i],
+                    ref br));
+        }
+    }
+
+    [Fact]
+    public void Test000() => HybridUintRoundtrip(new JxlAnsHybridUIntConfiguration(0, 0, 0));
+
+    [Fact]
+    public void Test411() => HybridUintRoundtrip(new JxlAnsHybridUIntConfiguration(4, 1, 1));
+
+    [Fact]
+    public void Test420() => HybridUintRoundtrip(new JxlAnsHybridUIntConfiguration(4, 2, 0));
+
+    [Fact]
+    public void Test421() => HybridUintRoundtrip(
+            new JxlAnsHybridUIntConfiguration(4, 2, 1),
+            256);
+}

tests/ImageSharp.Tests/Formats/Jxl/Processing/GaborishTests.cs

@@ -0,0 +1,113 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Runtime.CompilerServices;
+using System.Runtime.Intrinsics;
+using SixLabors.ImageSharp.Formats.Jxl.Memory;
+using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes;
+using SixLabors.ImageSharp.Formats.Jxl.Processing;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Image;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl.Processing;
+
+public class GaborishTests
+{
+    private static JxlWeightsSymmetric3 GaborishKernel(float weight1, float weight2)
+    {
+        const float weight0 = 1.0f;
+
+        // Normalize.
+        float mul = 1.0f / (weight0 + (4.0f * (weight1 + weight2)));
+        float w0 = weight0 * mul;
+        float w1 = weight1 * mul;
+        float w2 = weight2 * mul;
+
+        JxlWeightsSymmetric3 wgt = new();
+        wgt.SetC(Vector128.Create(w0));
+        wgt.SetD(Vector128.Create(w1));
+        wgt.SetR(Vector128.Create(w2));
+
+        return wgt;
+    }
+
+    private static void ConvolveGaborish(
+        JxlPlane<float> input,
+        float weight1,
+        float weight2,
+        JxlPlane<float> output)
+    {
+        Assert.True(JxlImageOperations.SameSize(input, output));
+
+        // The test will fail if this throws.
+        JxlConvolve.SlowSymmetric3(
+            input,
+            input.GetRectangle(),
+            GaborishKernel(weight1, weight2),
+            output);
+    }
+
+    private static void TestRoundTrip(
+        JxlImage3F input,
+        float maxL1)
+    {
+        JxlImage3F fwd = new(TestEnvironment.Configuration, input.XSize, input.YSize);
+
+        ConvolveGaborish(
+            input.Plane(0),
+            0,
+            0,
+            fwd.Plane(0));
+
+        ConvolveGaborish(
+            input.Plane(1),
+            0,
+            0,
+            fwd.Plane(1));
+
+        ConvolveGaborish(
+            input.Plane(2),
+            0,
+            0,
+            fwd.Plane(2));
+
+        const float w = 0.92718927264540152f;
+
+        InlineArray3<float> weights = default;
+        weights[0] = w;
+        weights[1] = w;
+        weights[2] = w;
+
+        // The test will fail if this throws.
+        JxlGaborish.InverseGaborish(
+            TestEnvironment.Configuration,
+            fwd,
+            fwd.GetRectangle(),
+            weights);
+
+        // TODO: VerifyRelativeError
+        Assert.True(
+            VerifyRelativeError(
+                input,
+                fwd,
+                maxL1,
+                1e-4f));
+    }
+
+    [Fact]
+    public void TestZero()
+    {
+        JxlImage3F input = new(TestEnvironment.Configuration, 20, 20);
+        JxlImageOperations.ZeroFillImage(input);
+        TestRoundTrip(input, 0.0f);
+    }
+
+    [Fact]
+    public void TestFlat()
+    {
+        JxlImage3F input = new(TestEnvironment.Configuration, 20, 20);
+        JxlImageOperations.FillImage(1.0f, input);
+        TestRoundTrip(input, 1e-5f);
+    }
+}

tests/ImageSharp.Tests/Formats/Jxl/Processing/ImageOperationsTests.cs

@@ -0,0 +1,146 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Jxl.Memory;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Image;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl.Processing;
+
+public class ImageOperationsTests
+{
+    private static void TestFillImpl<T>(JxlImage3<T> image, string layout)
+        where T : unmanaged
+    {
+        JxlImageOperations.FillImage((T)Convert.ChangeType(1, typeof(T), null), image);
+
+        for (int y = 0; y < image.YSize; y++)
+        {
+            for (int c = 0; c < 3; c++)
+            {
+                Span<T> row = image.PlaneRow(c, y);
+
+                for (int x = 0; x < image.XSize; x++)
+                {
+                    if (!EqualityComparer<T>.Default.Equals(
+                        row[x],
+                        (T)Convert.ChangeType(1, typeof(T), null)))
+                    {
+                        Assert.Fail(
+                            $"Not 1 at c={c} {x}, {y} " +
+                            $"({image.XSize} x {image.YSize}) ({layout})");
+                    }
+
+                    row[x] = (T)Convert.ChangeType(2, typeof(T), null);
+                }
+            }
+        }
+
+        JxlImageOperations.ZeroFillImage(image);
+
+        for (int c = 0; c < 3; c++)
+        {
+            for (int y = 0; y < image.YSize; y++)
+            {
+                Span<T> row = image.PlaneRow(c, y);
+
+                for (int x = 0; x < image.XSize; x++)
+                {
+                    if (!EqualityComparer<T>.Default.Equals(
+                        row[x],
+                        default))
+                    {
+                        Assert.Fail(
+                            $"Not 0 at c={c} {x}, {y} " +
+                            $"({image.XSize} x {image.YSize}) ({layout})");
+                    }
+
+                    row[x] = (T)Convert.ChangeType(3, typeof(T), null);
+                }
+            }
+        }
+    }
+
+    private static void TestFillT<T>()
+        where T : unmanaged
+    {
+        foreach (uint xSize in new uint[] { 0, 1, 15, 16, 31, 32 })
+        {
+            foreach (uint ySize in new uint[] { 0, 1, 15, 16, 31, 32 })
+            {
+                JxlImage3<T> image = JxlImage3<T>.Create(TestEnvironment.Configuration, (int)xSize, (int)ySize);
+
+                TestFillImpl(image, "size ctor");
+            }
+        }
+    }
+
+    [Fact]
+    public void TestFill()
+    {
+        TestFillT<byte>();
+        TestFillT<short>();
+        TestFillT<float>();
+    }
+
+    [Fact]
+    public void CopyImageToWithPaddingTest()
+    {
+        JxlPlane<uint> src = JxlPlane<uint>.Create(TestEnvironment.Configuration, 100, 61);
+
+        for (int y = 0; y < src.YSize; y++)
+        {
+            Span<uint> row = src.GetRow(y);
+
+            for (int x = 0; x < src.XSize; x++)
+            {
+                row[x] = (uint)((x * 1000) + y);
+            }
+        }
+
+        Rectangle srcRect = new(10, 20, 30, 40);
+        Assert.True(srcRect.IsInside(src.GetRectangle()));
+
+        JxlPlane<uint> dst = JxlPlane<uint>.Create(TestEnvironment.Configuration, 60, 50);
+
+        JxlImageOperations.FillImage(0u, dst);
+
+        Rectangle dstRect = new(20, 5, 30, 40);
+        Assert.True(dstRect.IsInside(dst.GetRectangle()));
+
+        Assert.True(
+            JxlImageOperations.CopyImageToWithPadding(
+                srcRect,
+                src,
+                padding: 2,
+                dstRect,
+                dst));
+
+        Rectangle paddedDstRect = new(
+            20 - 2,
+            5 - 2,
+            30 + 4,
+            40 + 3);
+
+        for (int y = 0; y < dst.YSize; y++)
+        {
+            Span<uint> row = dst.GetRow(y);
+
+            for (int x = 0; x < dst.XSize; x++)
+            {
+                if (new Rectangle(x, y, 1, 1).IsInside(paddedDstRect))
+                {
+                    Assert.Equal(
+                        (uint)(
+                            ((x - dstRect.X0() + srcRect.X0()) * 1000) +
+                            (y - dstRect.Y0() + srcRect.Y0())),
+                        row[x]);
+                }
+                else
+                {
+                    Assert.Equal(0u, row[x]);
+                }
+            }
+        }
+    }
+}

tests/ImageSharp.Tests/Formats/Jxl/TestUtils/TestDCT.cs

@@ -0,0 +1,108 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl.TestUtils;
+
+/// <summary>
+/// This is a naive Discrete Cosine Transform, without any optimization.
+/// It is a reference to test against the optimized DCT.
+/// </summary>
+internal static class TestDCT
+{
+    private static double Alpha(int u) => u == 0 ? 0.7071067811865475 : 1.0;
+
+    public static void Dct1D(ReadOnlySpan<double> block, Span<double> output, int n, int m)
+    {
+        Span<double> matrix = stackalloc double[n * n];
+        double scale = Math.Sqrt(2.0) / n;
+
+        for (int y = 0; y < n; y++)
+        {
+            for (int u = 0; u < n; u++)
+            {
+                matrix[(n * u) + y] =
+                    Alpha(u) *
+                    Math.Cos((y + 0.5) * u * Math.PI / n) *
+                    scale;
+            }
+        }
+
+        for (int x = 0; x < m; x++)
+        {
+            for (int u = 0; u < n; u++)
+            {
+                output[(m * u) + x] = 0;
+
+                for (int y = 0; y < n; y++)
+                {
+                    output[(m * u) + x] += matrix[(n * u) + y] * block[(m * y) + x];
+                }
+            }
+        }
+    }
+
+    public static void Idct1D(ReadOnlySpan<double> block, Span<double> output, int n, int m)
+    {
+        Span<double> matrix = stackalloc double[n * n];
+        double scale = Math.Sqrt(2.0);
+
+        for (int y = 0; y < n; y++)
+        {
+            for (int u = 0; u < n; u++)
+            {
+                matrix[(n * y) + u] =
+                    Alpha(u) *
+                    Math.Cos((y + 0.5) * u * Math.PI / n) *
+                    scale;
+            }
+        }
+
+        for (int x = 0; x < m; x++)
+        {
+            for (int u = 0; u < n; u++)
+            {
+                output[(m * u) + x] = 0;
+
+                for (int y = 0; y < n; y++)
+                {
+                    output[(m * u) + x] += matrix[(n * u) + y] * block[(m * y) + x];
+                }
+            }
+        }
+    }
+
+    public static void TransposeBlock(ReadOnlySpan<double> input, Span<double> output, int n, int m)
+    {
+        for (int x = 0; x < n; x++)
+        {
+            for (int y = 0; y < m; y++)
+            {
+                output[(y * n) + x] = input[(x * m) + y];
+            }
+        }
+    }
+
+    public static void DctSlow(Span<double> block, int n)
+    {
+        int blockSize = n * n;
+        Span<double> g = stackalloc double[blockSize];
+
+        Dct1D(block, g, n, n);
+        TransposeBlock(g, block, n, n);
+
+        Dct1D(block, g, n, n);
+        TransposeBlock(g, block, n, n);
+    }
+
+    public static void IdctSlow(Span<double> block, int n)
+    {
+        int blockSize = n * n;
+        Span<double> g = stackalloc double[blockSize];
+
+        Idct1D(block, g, n, n);
+        TransposeBlock(g, block, n, n);
+
+        Idct1D(block, g, n, n);
+        TransposeBlock(g, block, n, n);
+    }
+}