📄 Colors/src/OklchConverter.cs
using System.Numerics;
using SixLabors.ImageSharp;
using SixLabors.ImageSharp.PixelFormats;

namespace Mosaic.Colors;

public static class OklchConverter
{
    extension(Image<Rgba32> image)
    {
        public void ConvertToScaledOklch()
        {
            image.ProcessPixelRows(accessor =>
            {
                for (var y = 0; y < accessor.Height; y++)
                {
                    var row = accessor.GetRowSpan(y);
                    for (var x = 0; x < row.Length; x++)
                    {
                        row[x] = row[x].ToScaledOklch();
                    }
                }
            });
        }
    }

    extension(Rgba32 pixel)
    {
        public Rgba32 ToScaledOklch()
        {
            var oklch = pixel.ToOklch();
            return new(oklch.X, oklch.Y / 0.4f, oklch.Z / 360.0f);
        }

        public Vector3 ToOklch()
        {
            var lr = SrgbToLinear(pixel.R);
            var lg = SrgbToLinear(pixel.G);
            var lb = SrgbToLinear(pixel.B);

            // Linear sRGB → LMS (using Oklab M1 matrix)
            var l_ = 0.4122214708f * lr + 0.5363325363f * lg + 0.0514459929f * lb;
            var m_ = 0.2119034982f * lr + 0.6806995451f * lg + 0.1073969566f * lb;
            var s_ = 0.0883024619f * lr + 0.2817188376f * lg + 0.6299787005f * lb;

            // Cube root (LMS → Lab cone response)
            var l3 = float.Cbrt(l_);
            var m3 = float.Cbrt(m_);
            var s3 = float.Cbrt(s_);

            // LMS cone response → OKLab
            var L = 0.2104542553f * l3 + 0.7936177850f * m3 - 0.0040720468f * s3;
            var a = 1.9779984951f * l3 - 2.4285922050f * m3 + 0.4505937099f * s3;
            var bLab = 0.0259040371f * l3 + 0.7827717662f * m3 - 0.8086757660f * s3;

            // OKLab → OKLCH
            var C = float.Sqrt(a * a + bLab * bLab);
            var H = float.Atan2(bLab, a) * (180.0f / float.Pi);
            while (H < 0.0f)
            {
                H += 360.0f;
            }

            return new(L, C, H);
        }

        public Rgba32 FromScaledOklch() =>
            Rgba32.FromOklch(new(pixel.R / 255.0f, pixel.G / 255.0f * 0.4f, pixel.B / 255.0f * 360.0f));

        public static Rgba32 FromOklch(Vector3 oklch)
        {
            // OKLCH → OKLab
            var hRad = oklch.Z * (float.Pi / 180.0f);
            var a = oklch.Y * float.Cos(hRad);
            var bLab = oklch.Y * float.Sin(hRad);

            // OKLab → LMS cone response
            var l3 = oklch.X + 0.3963377774f * a + 0.2158037573f * bLab;
            var m3 = oklch.X - 0.1055613458f * a - 0.0638541728f * bLab;
            var s3 = oklch.X - 0.0894841775f * a - 1.2914855480f * bLab;

            // Cube (cone response → LMS)
            var l_ = l3 * l3 * l3;
            var m_ = m3 * m3 * m3;
            var s_ = s3 * s3 * s3;

            // LMS → linear sRGB (inverse of M1)
            var lr = +4.0767416621f * l_ - 3.3077115913f * m_ + 0.2309699292f * s_;
            var lg = -1.2684380046f * l_ + 2.6097574011f * m_ - 0.3413193965f * s_;
            var lb = -0.0041960863f * l_ - 0.7034186147f * m_ + 1.7076147010f * s_;

            return new(LinearToSrgb(lr), LinearToSrgb(lg), LinearToSrgb(lb));
        }
    }

    private static float SrgbToLinear(byte v)
    {
        var s = v / 255.0f;
        return s <= 0.04045f ? s / 12.92f : float.Pow((s + 0.055f) / 1.055f, 2.4f);
    }

    private static byte LinearToSrgb(float v)
    {
        var s = v <= 0.0031308f ? 12.92f * v : 1.055f * float.Pow(v, 1.0f / 2.4f) - 0.055f;
        return (byte)float.Round(float.Clamp(s * 255.0f, 0.0f, 255.0f));
    }
}