Support for JPEG XL (JXL) images

[winscripter]

Prerequisites

• I have written a descriptive pull-request title
• I have verified that there are no overlapping pull-requests open
• I have verified that I am following the existing coding patterns and practice as demonstrated in the repository. These follow strict Stylecop rules 👮.
• I have provided test coverage for my change (where applicable)

Description

This is a work-in-progress PR whose goal is to introduce decoding and encoding of JPEG XL (*.jxl) images.

Reference software
I use libjxl as reference. See https://github.com/libjxl/libjxl.

Performance
I will begin by applying light optimizations as I implement parts of the JPEG XL codec. Once the codec seems complete enough to handle decoding and encoding of JPEG XL images, I will apply heavier optimizations. Examples include but are not limited to stack allocation, array pooling, and SIMD.

Implementations
The JPEG XL codec lives under src/ImageSharp/Formats/Jxl.

Testing
I will start adding tests whenever the codec is complete enough to handle decoding of JPEG XL images.

Additionally, JPEG XL reference software, libjxl, contains its own tests too, which I might also implement without modification.

Progress

• 🔴 - work didn't start yet
• 🟠 - exists, but incomplete
• 🟡 - either decoding side, encoding side, or tests are complete, but not the other
• 🟢 - complete

---

• 🟡 AC strategy

• 🟢 AC strategy image/row

• 🟢 AC context

• 🔴 AC strategy tests

• 🟢 Adaptive Quantization (encoder)

• 🟠 ANS Entropy

• 🟢 ANS Entropy: Common

• 🟢 ANS Entropy: Common (Tests)

• 🟠 ANS Entropy: Decoder (symbol reader is incomplete)

• 🟠 ANS Entropy: Encoder (SIMD bit-cost calculation only)

• 🔴 ANS Entropy: Tests

• 🟢 Alpha Blending

• 🟡 Bit I/O

• 🟢 Bit I/O: Bit reader

• 🟢 Bit I/O: Bit writer

• 🔴 Bit I/O: tests

• 🟢 Box Content Decoder

• 🟢 Box Content Decoder: Uncompressed boxes

• 🟢 Box Content Decoder: Brotli-compressed boxes

• 🟢 Butteraugli

• 🟢 Butteraugli: Shared methods

• 🟢 Butteraugli: Abstractions

• 🟢 Butteraugli: Default comparator

• 🟢 Butteraugli: Encoder comparator

• 🟡 Cache

• 🟠 Cache: Decoder

• 🔴 Cache: Encoder

• 🟠 Chroma From Luma

• 🟢 Chroma From Luma Abstractions

• 🔴 Chroma From Luma Encoder

• 🟡 Coefficient Order

• 🟢 Coefficient Order: Forward

• 🟢 Coefficient Order: Main

• 🟢 Coefficient Order: Encoder

• 🔴 Coefficient Order: Tests

• 🟢 Compressed DC

• 🟡 Context Map

• 🟢 Context Map: Abstractions

• 🟢 Context Map: Decoder

• 🔴 Context Map: Encoder

• 🟡 Convolution

• 🟢 Convolution: Symmetric

• 🟢 Convolution: Separable

• 🟢 Convolution: Slow

• 🟢 Convolution: SIMD

• 🔴 Convolution: Separable5 Encoder

• 🔴 Convolution: Tests

• 🟠 Decoder: Frame

• 🟢 Decoder: Group

• 🟢 Decoder: Group Border

• 🟠 Decoder: Main

• 🟢 Decoder: Main: Codestream parser

• 🟢 Decoder: Main: Container format parser

• 🟡 Discrete Cosine Transform

• 🟢 Discrete Cosine Transform: DCT scales

• 🟢 Discrete Cosine Transform: Block data wrapper

• 🟠 Discrete Cosine Transform: Block-based

• 🟢 Discrete Cosine Transform: Slow DCT for reference in tests

• 🔴 Discrete Cosine Transform: Tests

• 🔴 Encoder dot detection

• 🔴 Encoder dot dictionary

• 🟠 Entropy coding

• 🔴 Encoder entropy coding

• 🟢 Fields

• 🟢 Fields: Visitor abstractions

• 🟢 Fields: Parser

• 🟢 Fields: Writer

• 🔴 Frame Encoder

• 🟢 Gaborish

• 🟢 Gaborish Encoder

• 🟢 Gaborish Tests

• 🔴 Group Encoder

• 🔴 Heuristics Encoder

• 🟢 Image Bundle

• 🟢 Image Bundle: Decoder/Common

• 🟢 Image Bundle: Encoder

• 🟡 LZ77 compression

• 🟢 LZ77: Fast Lossless Encoder

• 🔴 LZ77: Standard Encoder

• 🟢 Huffman compression

• 🟢 Huffman compression: Shared

• 🟢 Huffman compression: Decoder

• 🟢 Huffman compression: Encoder

• 🟡 Modular

• 🟢 Modular: Transforms

• 🟢 Modular: Transforms: Palette (Inverse)

• 🟢 Modular: Transforms: Palette (Forward)

• 🟢 Modular: Transforms: RCT (Inverse)

• 🟢 Modular: Transforms: RCT (Forward)

• 🟢 Modular: Transforms: Squeeze (Inverse)

• 🟢 Modular: Transforms: Squeeze (Forward)

• 🟢 Modular: Encoding

• 🟢 Modular: Encoding: Context Prediction

• 🟢 Modular: Encoding: MA decoder

• 🟢 Modular: Encoding: MA encoder

• 🟢 Modular: Encoding: Tree Samples

• 🟢 Modular: Encoding: Encoding decoder

• 🟢 Modular: Encoding: Encoding encoder

• 🔴 Modular: Decoder

• 🔴 Modular: Encoder

• 🔴 Modular: Encoder SIMD

• 🔴 Modular: Tests

• 🟠 Patch Dictionary: Decoder

• 🔴 Patch Dictionary: Encoder

• 🟠 Passes State: Decoder

• 🔴 Passes State: Encoder

• 🟢 Passes State: Shared

• 🔴 Encoder Main

• 🔴 Encoder Main

• 🔴 Encoder Internal

• 🔴 Encoder Tests

• 🟢 Encoder: Linear Algebra

• 🟢 Encoder: Linear Algebra Tests

• 🟢 Image Operations

• 🟢 Image Operations

• 🟢 Image Operations: Tests

• 🟢 Common I/O: Frame Header

• 🟢 Common I/O: Metadata

• 🟢 Common I/O: Container format

• 🟠 JPEG to JPEG XL lossless compression

• 🟢 JPEG to JPEG XL lossless compression: JPEG parser/writer

• 🟠 JPEG to JPEG XL lossless compression (decoder)

• 🔴 JPEG to JPEG XL lossless compression (encoder)

• 🟡 Splines

• 🟢 Splines

• 🔴 Splines Tests

• 🟢 Quantizer

• 🟢 Dequantizer matrices

• 🟢 Quantizer encoding

• 🟢 Quantizer weights

• 🟡 Noise

• 🟢 Noise: Shared

• 🟢 Noise: Decoder

• 🟠 Noise: Encoder

• 🟢 Noise: Simulation of Photon Noise

• 🟢 Noise: Simulation of Photon Noise (Tests)

• 🔴 Noise Tests

• 🔴 JPEG XL Testing Tools

• 🟠 Render Pipeline

• 🔴 Render Pipeline: Main

• 🔴 Render Pipeline: Low Memory Render Pipeline

• 🟢 Render Pipeline: Stages Abstractions

• 🔴 Render Pipeline: Stages: Blending

• 🟢 Render Pipeline: Stages: Chroma Upsampling

• 🔴 Render Pipeline: Stages: CMS

• 🟢 Render Pipeline: Stages: EPF

• 🟢 Render Pipeline: Stages: From Linear

• 🟢 Render Pipeline: Stages: Gaborish

• 🟢 Render Pipeline: Stages: Noise

• 🟢 Render Pipeline: Stages: Patches

• 🟢 Render Pipeline: Stages: Splines

• 🟢 Render Pipeline: Stages: Spot color

• 🟢 Render Pipeline: Stages: To Linear

• 🔴 Render Pipeline: Stages: Tone mapping

• 🔴 Render Pipeline: Stages: Upsampling

• 🟢 Render Pipeline: Stages: Write to Output

• 🔴 Render Pipeline: Stages: XYB

• 🟢 Render Pipeline: Stages: Y'Cb'Cr -> RGB

• 🟡 Color Management System (CMS)

• 🟢 CMS: Transfer Functions

• 🟢 CMS: Abstractions/Color Encoding

• 🔴 CMS: Tone Mapping

• 🔴 CMS: Interface

Other completed things:

• 🟢 Frame Dimensions
• 🟢 Field Encodings
• 🟢 Gamma correction (tests, encoder)
• 🟢 Shared math functions
• 🟢 Lehmer Codes
• 🟢 Lehmer Codes: Tests
• 🟢 Luminance helper
• 🟢 Opsin parameters
• 🟢 Passes shared state
• 🟢 SIMD utilities
• 🟢 TOC
• 🟢 XorShift128Plus
• 🟢 XorShift128Plus Tests

[CLAassistant]

All committers have signed the CLA.

[winscripter]

While I'm working on this, I'd like to note something important.

Libjxl is licensed under the BSD 3-Clause license, and since I'm using libjxl code as reference, that means the license must be included.

I'm not really sure what would be the proper way to include the license. I might place the LICENSE.txt file in the Jxl folder or add a README linking to the libjxl repo.

[gfoidl]

That's the same as above? Maybe extract it to a helper method?
And the * 0.25f on the 4 floats could be done vectorized in one pass then more easily.

[winscripter]

I don't actually know how to name this helper method, as this expression is just a core part of the Butterfly DCT algorithm.

As for vectorization, yes, it would probably bring performance benefits, but I'm not really sure how would one store the vector into the coefficients[0] coefficients[1] coefficients[8] coefficients[9] offsets instead of coefficient[0] through coefficient[3], it doesn't seem like vectors support something like this.

[gfoidl]

how to name this helper method

That's the hardest part 😉. Maybe just CoefficientMath?

it doesn't seem like vectors support something like this.

Gather is there on Avx2, but scatter (AVX512-F) hasn't landet and is tracked in dotnet/runtime#87097. So I'd just use the trivial coefficients[9] = vec128[3] idiom -- and hopefully a future JIT will recognise that pattern and emit a scatter store when profitable.

I tried a few trivial approaches, but none was really faster than the code as is. It's a gamer about ns, so not really worth it.
Maybe they tried it for Butterfly DCT too.

Benchmark results from trial

| Method  | Mean     | Error     | StdDev    | Ratio | RatioSD | Code Size |
|-------- |---------:|----------:|----------:|------:|--------:|----------:|
| Current | 5.000 ns | 0.0378 ns | 0.0316 ns |  1.00 |    0.01 |     203 B |
| A       | 6.246 ns | 0.0847 ns | 0.0751 ns |  1.25 |    0.02 |     206 B |
| B       | 5.977 ns | 0.0980 ns | 0.0869 ns |  1.20 |    0.02 |     255 B |
| C       | 5.527 ns | 0.0889 ns | 0.0832 ns |  1.11 |    0.02 |     214 B |

C# code

using System.Diagnostics;
using System.Runtime.InteropServices;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.X86;
using BenchmarkDotNet.Attributes;

#if !DEBUG
using BenchmarkDotNet.Running;
#endif

Bench bench = new();

bench.Setup();
bench.Current();
bench.Dump();

bench.Setup();
bench.A();
bench.Dump();

bench.Setup();
bench.B();
bench.Dump();

bench.Setup();
bench.C();
bench.Dump();

#if !DEBUG
BenchmarkRunner.Run<Bench>();
#endif

//[ShortRunJob]
[DisassemblyDiagnoser]
public class Bench
{
    private float[] _coefficients = null!;

    [GlobalSetup]
    public void Setup()
    {
        _coefficients = [.. Enumerable.Repeat(float.NaN, 10)];
        _coefficients[0] = 0f;
        _coefficients[1] = 1f;
        _coefficients[8] = 8f;
        _coefficients[9] = 9f;
    }

    [Benchmark(Baseline = true)]
    public void Current() => CurrentWorker(_coefficients);

    private static void CurrentWorker(Span<float> coefficients)
    {
        float block00 = coefficients[0];
        float block01 = coefficients[1];
        float block10 = coefficients[8];
        float block11 = coefficients[9];

        coefficients[0] = (block00 + block01 + block10 + block11) * 0.25f;
        coefficients[1] = (block00 + block01 - block10 - block11) * 0.25f;
        coefficients[8] = (block00 - block01 + block10 - block11) * 0.25f;
        coefficients[9] = (block00 - block01 - block10 + block11) * 0.25f;
    }

    [Benchmark]
    public void A() => AWorker(_coefficients);

    private static unsafe void AWorker(Span<float> coefficients)
    {
        float block00 = coefficients[0];
        float block01 = coefficients[1];
        float block10 = coefficients[8];
        float block11 = coefficients[9];

        Vector128<float> vec = Vector128.Create(
            block00 + block01 + block10 + block11,
            block00 + block01 - block10 - block11,
            block00 - block01 + block10 - block11,
            block00 - block01 - block10 + block11);

        vec *= 0.25f;

        coefficients[0] = vec[0];
        coefficients[1] = vec[1];
        coefficients[8] = vec[2];
        coefficients[9] = vec[3];
    }

    [Benchmark]
    public void B() => BWorker(_coefficients);

    private static unsafe void BWorker(Span<float> coefficients)
    {
        float c0 = coefficients[0];
        float c1 = coefficients[1];
        float c8 = coefficients[8];
        float c9 = coefficients[9];

        Vector128<float> vec = Vector128.Create(c0);
        vec += Vector128.Create(+c1, +c1, -c1, -c1);
        vec += Vector128.Create(+c8, -c8, +c8, -c8);
        vec += Vector128.Create(+c9, -c9, -c9, +c9);
        vec *= 0.25f;

        coefficients[0] = vec[0];
        coefficients[1] = vec[1];
        coefficients[8] = vec[2];
        coefficients[9] = vec[3];
    }

    [Benchmark]
    public void C() => CWorker(_coefficients);

    private static void CWorker(Span<float> coefficients)
    {
        Vector128<float> v0 = Vector128.Create(coefficients[0]);
        Vector128<float> v1 = Vector128.Create(coefficients[1]) * Vector128.Create(+1f, +1f, -1f, -1f);
        Vector128<float> v8 = Vector128.Create(coefficients[8]) * Vector128.Create(+1f, -1f, +1f, -1f);
        Vector128<float> v9 = Vector128.Create(coefficients[9]) * Vector128.Create(+1f, -1f, -1f, +1f);

        Vector128<float> res = (v0 + v1) + (v8 + v9);
        res *= 0.25f;

        coefficients[0] = res[0];
        coefficients[1] = res[1];
        coefficients[8] = res[2];
        coefficients[9] = res[3];
    }

    public void Dump()
    {
        Console.ForegroundColor =
               _coefficients[0] == 4.5f
            && _coefficients[1] == -4f
            && _coefficients[8] == -0.5f
            && _coefficients[9] == 0f
            ? ConsoleColor.Green
            : ConsoleColor.Red;

        Console.WriteLine($"{_coefficients[0]}\t{_coefficients[1]}\t{_coefficients[8]}\t{_coefficients[9]}");

        Console.ResetColor();
    }
}

dasm

; Bench.Current()
       sub       rsp,38
       vmovaps   [rsp+20],xmm6
       mov       rax,[rcx+8]
       test      rax,rax
       je        near ptr M00_L02
       lea       rcx,[rax+10]
       mov       eax,[rax+8]
M00_L00:
       cmp       eax,9
       jle       short M00_L03
       vmovss    xmm0,dword ptr [rcx]
       vmovss    xmm1,dword ptr [rcx+4]
       vmovss    xmm2,dword ptr [rcx+20]
       vmovss    xmm3,dword ptr [rcx+24]
M00_L01:
       vaddss    xmm4,xmm0,xmm1
       vaddss    xmm5,xmm4,xmm2
       vaddss    xmm5,xmm5,xmm3
       vmovss    xmm6,dword ptr [7FFEED85AEF0]
       vmulss    xmm5,xmm5,xmm6
       vmovss    dword ptr [rcx],xmm5
       vsubss    xmm4,xmm4,xmm2
       vsubss    xmm4,xmm4,xmm3
       vmulss    xmm4,xmm4,xmm6
       vmovss    dword ptr [rcx+4],xmm4
       vsubss    xmm0,xmm0,xmm1
       vaddss    xmm1,xmm0,xmm2
       vsubss    xmm1,xmm1,xmm3
       vmulss    xmm1,xmm1,xmm6
       vmovss    dword ptr [rcx+20],xmm1
       vsubss    xmm0,xmm0,xmm2
       vaddss    xmm0,xmm0,xmm3
       vmulss    xmm0,xmm0,xmm6
       vmovss    dword ptr [rcx+24],xmm0
       vmovaps   xmm6,[rsp+20]
       add       rsp,38
       ret
M00_L02:
       xor       ecx,ecx
       xor       eax,eax
       jmp       short M00_L00
M00_L03:
       test      eax,eax
       je        short M00_L04
       vmovss    xmm0,dword ptr [rcx]
       cmp       eax,1
       jbe       short M00_L04
       vmovss    xmm1,dword ptr [rcx+4]
       cmp       eax,8
       jbe       short M00_L04
       vmovss    xmm2,dword ptr [rcx+20]
       cmp       eax,9
       jbe       short M00_L04
       vmovss    xmm3,dword ptr [rcx+24]
       jmp       near ptr M00_L01
M00_L04:
       call      CORINFO_HELP_RNGCHKFAIL
       int       3
; Total bytes of code 203

; Bench.A()
       sub       rsp,28
       mov       rax,[rcx+8]
       test      rax,rax
       je        near ptr M00_L02
       lea       rcx,[rax+10]
       mov       eax,[rax+8]
M00_L00:
       cmp       eax,9
       jle       near ptr M00_L03
       vmovss    xmm0,dword ptr [rcx]
       vmovss    xmm1,dword ptr [rcx+4]
       vmovss    xmm2,dword ptr [rcx+20]
       vmovss    xmm3,dword ptr [rcx+24]
M00_L01:
       vaddss    xmm4,xmm0,xmm1
       vaddss    xmm5,xmm4,xmm2
       vaddss    xmm5,xmm5,xmm3
       vsubss    xmm4,xmm4,xmm2
       vsubss    xmm4,xmm4,xmm3
       vinsertps xmm4,xmm5,xmm4,10
       vsubss    xmm0,xmm0,xmm1
       vaddss    xmm1,xmm0,xmm2
       vsubss    xmm1,xmm1,xmm3
       vinsertps xmm1,xmm4,xmm1,20
       vsubss    xmm0,xmm0,xmm2
       vaddss    xmm0,xmm0,xmm3
       vinsertps xmm0,xmm1,xmm0,30
       vmulps    xmm0,xmm0,[7FFEED85B0E0]
       vmovss    dword ptr [rcx],xmm0
       vextractps dword ptr [rcx+4],xmm0,1
       vextractps dword ptr [rcx+20],xmm0,2
       vextractps dword ptr [rcx+24],xmm0,3
       add       rsp,28
       ret
M00_L02:
       xor       ecx,ecx
       xor       eax,eax
       jmp       near ptr M00_L00
M00_L03:
       test      eax,eax
       je        short M00_L04
       vmovss    xmm0,dword ptr [rcx]
       cmp       eax,1
       jbe       short M00_L04
       vmovss    xmm1,dword ptr [rcx+4]
       cmp       eax,8
       jbe       short M00_L04
       vmovss    xmm2,dword ptr [rcx+20]
       cmp       eax,9
       jbe       short M00_L04
       vmovss    xmm3,dword ptr [rcx+24]
       jmp       near ptr M00_L01
M00_L04:
       call      CORINFO_HELP_RNGCHKFAIL
       int       3
; Total bytes of code 206

; Bench.B()
       sub       rsp,28
       mov       rax,[rcx+8]
       test      rax,rax
       je        near ptr M00_L02
       lea       rcx,[rax+10]
       mov       eax,[rax+8]
M00_L00:
       cmp       eax,9
       jle       near ptr M00_L03
       vmovss    xmm0,dword ptr [rcx]
       vmovss    xmm1,dword ptr [rcx+4]
       vmovss    xmm2,dword ptr [rcx+20]
       vmovss    xmm3,dword ptr [rcx+24]
M00_L01:
       vmovaps   xmm4,xmm1
       vinsertps xmm4,xmm4,xmm1,10
       vxorps    xmm1,xmm1,[7FFEED87B1E0]
       vinsertps xmm4,xmm4,xmm1,20
       vinsertps xmm1,xmm4,xmm1,30
       vbroadcastss xmm0,xmm0
       vaddps    xmm0,xmm1,xmm0
       vmovaps   xmm1,xmm2
       vxorps    xmm4,xmm2,[7FFEED87B1E0]
       vinsertps xmm1,xmm1,xmm4,10
       vinsertps xmm1,xmm1,xmm2,20
       vinsertps xmm1,xmm1,xmm4,30
       vaddps    xmm0,xmm1,xmm0
       vmovaps   xmm1,xmm3
       vxorps    xmm2,xmm3,[7FFEED87B1E0]
       vinsertps xmm1,xmm1,xmm2,10
       vinsertps xmm1,xmm1,xmm2,20
       vinsertps xmm1,xmm1,xmm3,30
       vaddps    xmm0,xmm1,xmm0
       vmulps    xmm0,xmm0,[7FFEED87B1F0]
       vmovss    dword ptr [rcx],xmm0
       vextractps dword ptr [rcx+4],xmm0,1
       vextractps dword ptr [rcx+20],xmm0,2
       vextractps dword ptr [rcx+24],xmm0,3
       add       rsp,28
       ret
M00_L02:
       xor       ecx,ecx
       xor       eax,eax
       jmp       near ptr M00_L00
M00_L03:
       test      eax,eax
       je        short M00_L04
       vmovss    xmm0,dword ptr [rcx]
       cmp       eax,1
       jbe       short M00_L04
       vmovss    xmm1,dword ptr [rcx+4]
       cmp       eax,8
       jbe       short M00_L04
       vmovss    xmm2,dword ptr [rcx+20]
       cmp       eax,9
       jbe       short M00_L04
       vmovss    xmm3,dword ptr [rcx+24]
       jmp       near ptr M00_L01
M00_L04:
       call      CORINFO_HELP_RNGCHKFAIL
       int       3
; Total bytes of code 255

; Bench.C()
       sub       rsp,28
       mov       rax,[rcx+8]
       test      rax,rax
       je        short M00_L02
       lea       rcx,[rax+10]
       mov       eax,[rax+8]
M00_L00:
       cmp       eax,9
       jle       short M00_L03
       vbroadcastss xmm0,dword ptr [rcx]
       vbroadcastss xmm1,dword ptr [rcx+4]
       vmulps    xmm1,xmm1,[7FFEED86B130]
       vaddps    xmm0,xmm1,xmm0
       vbroadcastss xmm1,dword ptr [rcx+20]
       vmulps    xmm1,xmm1,[7FFEED86B140]
       vbroadcastss xmm2,dword ptr [rcx+24]
       vmulps    xmm2,xmm2,[7FFEED86B150]
       vaddps    xmm1,xmm2,xmm1
       vaddps    xmm0,xmm1,xmm0
M00_L01:
       vmulps    xmm0,xmm0,[7FFEED86B160]
       vmovss    dword ptr [rcx],xmm0
       vextractps dword ptr [rcx+4],xmm0,1
       vextractps dword ptr [rcx+20],xmm0,2
       vextractps dword ptr [rcx+24],xmm0,3
       add       rsp,28
       ret
M00_L02:
       xor       ecx,ecx
       xor       eax,eax
       jmp       short M00_L00
M00_L03:
       test      eax,eax
       je        short M00_L04
       vbroadcastss xmm0,dword ptr [rcx]
       cmp       eax,1
       jbe       short M00_L04
       vbroadcastss xmm1,dword ptr [rcx+4]
       vmulps    xmm1,xmm1,[7FFEED86B130]
       vaddps    xmm0,xmm1,xmm0
       cmp       eax,8
       jbe       short M00_L04
       vbroadcastss xmm1,dword ptr [rcx+20]
       vmulps    xmm1,xmm1,[7FFEED86B140]
       cmp       eax,9
       jbe       short M00_L04
       vbroadcastss xmm2,dword ptr [rcx+24]
       vmulps    xmm2,xmm2,[7FFEED86B150]
       vaddps    xmm1,xmm2,xmm1
       vaddps    xmm0,xmm1,xmm0
       jmp       short M00_L01
M00_L04:
       call      CORINFO_HELP_RNGCHKFAIL
       int       3
; Total bytes of code 214

[!NOTE]
It's cool to see how the JIT clones the code to avoid the bound checks.

PS: I'll have a look at the new commits the next days.

[winscripter]

Gather is there on Avx2, but scatter (AVX512-F) hasn't

Yep, but... even if .NET did support AVX512F scatter support, there would still be a problem. We only need 4 floats, and that's 128 bits. But AVX512F is 512-bit, and that's 16 floats. So even with scatter support, there would be too many floats.

Maybe just CoefficientMath?

Yup, I'll go ahead and use this as the method name.

[gfoidl]

But AVX512F is 512-bit

That's not true. Vector512<T> (or the native __m512) is 512 bits.
But the instruction set named AVX512 (with all the sub-categories) brings in just the instructions, and there are some that target Vector128<T> (__m128).

Same as e.g. _mm_i32gather_epi32 is part of AVX2 instruction set, but the operands are Vector128<T>, and not 256 bits.

Here the instruction for the scatter is _mm_mask_i32scatter_epi32 and is part of AVX512F + AVX512VL.

[winscripter]

there are some that target Vector128 (__m128).

but the operands are Vector128, and not 256 bits.

Very good to know! 💻✨

[gfoidl]

Is tableValues.Length the correct check here?

Let's assume AVX2, then Vector<int>.Count = 8, so table has 8 ints.
In tableValues a int-span with length 32 is created. Then in L182 the vec is copied to span, so only the first 8 items of the 32 have actually values. All others are either 0 or any gargabe that from stack alloc is there.

[gfoidl]

This is selecting elements from table as given by indices, so kind of shuffling.
Could then this be written as shuffle which is supported on Vector128, et. al.?

[gfoidl]

Have it as Vector3 argument? Or inline array?
In order to avoid the bound checks needed below in L50.

[gfoidl]

Store the Vector3 here directly for luminance?
The values are used only once , so no need to unpack the vector in the ctor beforehands.

[gfoidl]

When it's a Vector3 then this is rgb *= ratio and done.

[gfoidl]

Hm, just reading the code and seeing that in L45 int i = -2 this makes me feel it crashes here w/ a index out of range exception. Can you add an assert or comment why it works? I.e. that Debug.Assert(xStart >= 2); holds, etc.

[winscripter]

I'll explain this for row0 but it is the same idea for every other variable accessed like this.

On L32 we can see that row0 is defined as such:

Span<float> row0 = this.GetInputRow(inputRows, c, -2);

which basically says, get a row from the input rows at channel c for row offset -2.

Let's take a look at the GetInputRow method.

public Span<float> GetInputRow(Buffer2D<Memory<float>> inputRows, int c, int offset)
    => inputRows[c, this.Settings.BorderY + offset].Span[RenderPipelineXOffset..];

This returns a Memory<float> stored for channel c at offset this.Settings.BorderY + offset. In this case, this.Settings.BorderY would be equal to 2, because of this constructor in the ConvolveNoiseStage:

public ConvolveNoiseStage(Configuration configuration, int firstC)
    : base(configuration)
{
    this.Settings = RenderPipelineStageConfiguration.CreateSymmetricBorderOnly(2);
    this.firstC = firstC;
}

Specifically, this line:

this.Settings = RenderPipelineStageConfiguration.CreateSymmetricBorderOnly(2);

That means, create RenderPipelineStageConfiguration whose ShiftX, ShiftY = 0 and BorderX, BorderY = input parameter (2). So BorderY=2.

This gives us a Memory<float>, but we also get a Span to it, and:

.Span[RenderPipelineXOffset..]

we advance it by RenderPipelineStageXOffset.

That means it must start at 32 which is what the constant equals to.

Now, let's take a look at reference software, which is in C++. Let's take a look at their implementation of this GetInputRow method.

Note. RowInfo = std::vector<std::vector<float*>>.

// Returns a pointer to the input row of channel `c` with offset `y`.
// `y` must be in [-settings_.border_y, settings_.border_y]. `c` must be such
// that `GetChannelMode(c) != kIgnored`. The returned pointer points to the
// offset-ed row (i.e. kRenderPipelineXOffset has been applied).
float* GetInputRow(const RowInfo& input_rows, size_t c, int offset) const {
  JXL_DASSERT(GetChannelMode(c) != RenderPipelineChannelMode::kIgnored);
  JXL_DASSERT(-offset <= static_cast<int>(settings_.border_y));
  JXL_DASSERT(offset <= static_cast<int>(settings_.border_y));
  return input_rows[c][settings_.border_y + offset] + kRenderPipelineXOffset;
}

When we simplify it:

float* GetInputRow(const RowInfo& input_rows, size_t c, int offset) const {
  return input_rows[c][settings_.border_y + offset] + kRenderPipelineXOffset;
}

It's the same idea... or is it? That's because this C++ method returns a pointer to a float, while a Span would have bounds checks against values below 0. In this C++ version, negative accesses are fine here because the offset is kRenderPipelineXOffset (= 32) and the lowest offset would be 30, so it wouldn't lead to program termination. We should not assume that xStart >= 2.

That means we need to perhaps introduce a new method like this:

public Span<float> GetInputRow(Buffer2D<Memory<float>> inputRows, int c, int offset, int spanOffset)
    => inputRows[c, this.Settings.BorderY + offset].Span[(spanOffset + RenderPipelineXOffset)..];

Then we have to pass -2 as the spanOffset and pass +2 to every Span access offset, meaning the base offset is now 30, the 2 + -2 access would just be 0 and the 2+3 access would be 5, giving us the access to the 35th offset.

With that being said, we'd have to update this:

others += Vector.Create<float>(row0[(x + i)..]);

to:

others += Vector.Create<float>(row0[(x + 2 + i)..]);

Also. xStart cannot be negative.

[gfoidl]

Thanks for the explanation.

[gfoidl]

Use an abstract class here?
Makes it easier for the JIT to perform de-abstraction. Interfaces are harder (as a type can have multiple interface, but exactly one base class (even when it's the object type)).

[gfoidl]

All of these loops should be

for (int i = 0; i <= length - Vector<float>.Count; i += Vector<float>.Count)

to have the proper ranges, and to don't read too less.

Let's assume Vector<float>.Count = 4, length = 3 then for i = 0 the test i < length evaluates to true, and it would read a fulle vector (i.e. 4 elements) out of 3, which is not correct.

With the suggested loop declaration it's i <= 3 - 4 and false which is correct.

So that's a off-by-one error.

[gfoidl]

Please have a unit test for this with different types for T and see what happens 😉.

[gfoidl]

considered slow

It could be written as

if (Vector<T>.Count == Vector64<T>.Count)
{
    throw new NotSupportedException();
}
else if (Vector<T>.Count == Vector512<T>.Count)
{
    return Vector512.Create(vec).AsVector();
}
else if (Vector<T>.Count == Vector256<T>.Count)
{
    return Vector256.Create(vec).AsVector();
}
else if (Vector<T>.Count == Vector128<T>.Count)
{
    return vec.AsVector();
}
else
{
    Span<T> value = stackalloc T[Vector<T>.Count];
    for (int i = 0; i <= Vector<T>.Count - Vector128<T>.Count; i += Vector128<T>.Count)
    {
        vec.CopyTo(value[i..]);
    }

    return new Vector<T>(value);
}

in order to make it faster for the ${128, 256, 512}$ bits.
As ATM and in the near future IMO only ARM will have wider vectors, maybe there's ome special intrinsics to use for this (which I'm not aware of). Or .NET will extend Vector.Create<T> to support a similar duplication pattern as the VectorXYZ do.