playpen/src/oklab.zig

//! sRGB in, sRGB out, with a perceptual space in the middle.
//!
//! This exists for `tint.zig`, which builds a whole palette out of one colour
//! and therefore has to answer questions like "the same hue, a third as bright"
//! and "one step lighter than that". Those questions have no good answer in RGB
//! — halving the channels of a mid blue gives a navy, halving the channels of a
//! mid yellow gives an olive, and the two do not look like they moved by the
//! same amount. HSL is no better: its lightness is the midpoint of the largest
//! and smallest channel, so a pure yellow and a pure blue are both "50% light"
//! when one of them is nearly white and the other is nearly black.
//!
//! OKLab is a fit to what people actually report seeing, and its lightness axis
//! is even enough that a ramp built by stepping it reads as an even ramp. The
//! polar form — lightness, chroma, hue — is what a palette is really made of:
//! hold the hue, walk the lightness, and the surfaces of a theme fall out.
//!
//! The one thing OKLCh will happily do that sRGB will not is name a colour that
//! does not exist on a monitor — a fully saturated yellow at the lightness of a
//! midtone, say. `toRgb` deals with that by giving up chroma rather than
//! lightness or hue: the result is the most colourful version of the colour that
//! can actually be shown, which is what "as close as the screen gets" should
//! mean for a palette. Clipping the channels instead would shift the hue, and a
//! generated theme whose reds drift orange as they darken looks broken in a way
//! that a slightly duller red does not.
//!
//! The matrices are Björn Ottosson's, unchanged.

const std = @import("std");

const palette = @import("palette.zig");

const Rgb = palette.Rgb;

/// A colour in OKLCh: how light, how colourful, and which colour.
pub const Lch = struct {
    /// Perceptual lightness. 0 is black, 1 is white, and 0.5 looks like a
    /// midtone rather than merely being one arithmetically.
    l: f32,

    /// Colourfulness. 0 is a grey; how high it can go before leaving sRGB
    /// depends on both the lightness and the hue, and is around 0.32 at best.
    c: f32,

    /// Hue angle in degrees. Roughly: 30 red, 100 yellow, 145 green, 195 cyan,
    /// 260 blue, 330 magenta.
    h: f32,
};

pub fn fromRgb(rgb: Rgb) Lch {
    const lin: Linear = .{
        .r = decode(rgb.r),
        .g = decode(rgb.g),
        .b = decode(rgb.b),
    };

    const lab = lin.toLab();
    const c = std.math.hypot(lab.a, lab.b);

    return .{
        .l = lab.l,
        .c = c,
        // A grey has no hue to report, and `atan2` on two zeroes is entitled to
        // say anything. Zero is as good as any other answer and is at least the
        // same one every time, which matters: the palette generator reads a hue
        // off the base colour and gives it to forty other colours, and a grey
        // base that produced a different hue on each launch would be a theme
        // that changed colour when you restarted the app.
        .h = if (c < 1e-6) 0 else std.math.radiansToDegrees(std.math.atan2(lab.b, lab.a)),
    };
}

/// The nearest colour a screen can show, giving up chroma before anything else.
pub fn toRgb(lch: Lch) Rgb {
    const l = std.math.clamp(lch.l, 0, 1);
    const rad = std.math.degreesToRadians(lch.h);
    const c = fit(l, @max(lch.c, 0), rad);

    const lin = Lab.at(l, c, rad).toLinear();
    return .{
        .r = encode(lin.r),
        .g = encode(lin.g),
        .b = encode(lin.b),
    };
}

/// The largest chroma at or below `c` that stays inside sRGB, to within a
/// rounding error of the 8-bit channels this is on its way to.
///
/// A bisection rather than a formula because the sRGB gamut boundary in OKLab
/// is not one: it is the image of a cube through a cube root, and the closed
/// forms for it are approximations with their own error. Twelve halvings of a
/// range that is at most 1.0 wide lands well inside a 1/255 step, and this runs
/// forty times when someone drags a colour picker.
fn fit(l: f32, c: f32, rad: f32) f32 {
    if (Lab.at(l, c, rad).inGamut()) return c;

    var lo: f32 = 0;
    var hi: f32 = c;
    for (0..12) |_| {
        const mid = (lo + hi) / 2;
        if (Lab.at(l, mid, rad).inGamut()) lo = mid else hi = mid;
    }
    return lo;
}

// -------------------------------------------------------------------------
// The two conversions, and the linear-light stage between them.

/// Light as the eye's cone responses model it: perceptual lightness, and two
/// opponent axes that carry the hue and how much of it there is.
const Lab = struct {
    l: f32,
    a: f32,
    b: f32,

    fn at(l: f32, c: f32, rad: f32) Lab {
        return .{ .l = l, .a = c * @cos(rad), .b = c * @sin(rad) };
    }

    fn toLinear(self: Lab) Linear {
        const l_ = self.l + 0.3963377774 * self.a + 0.2158037573 * self.b;
        const m_ = self.l - 0.1055613458 * self.a - 0.0638541728 * self.b;
        const s_ = self.l - 0.0894841775 * self.a - 1.2914855480 * self.b;

        const l = l_ * l_ * l_;
        const m = m_ * m_ * m_;
        const s = s_ * s_ * s_;

        return .{
            .r = 4.0767416621 * l - 3.3077115913 * m + 0.2309699292 * s,
            .g = -1.2684380046 * l + 2.6097574011 * m - 0.3413193965 * s,
            .b = -0.0041960863 * l - 0.7034186147 * m + 1.7076147010 * s,
        };
    }

    /// Whether this colour is one a screen can show. The tolerance is a hair
    /// over half of an 8-bit step in linear light near the top of the range, so
    /// a colour that is out of gamut only by the arithmetic isn't hunted down
    /// by the bisection above for no visible gain.
    fn inGamut(self: Lab) bool {
        const lin = self.toLinear();
        const tolerance = 1e-4;
        for ([_]f32{ lin.r, lin.g, lin.b }) |channel| {
            if (channel < -tolerance or channel > 1 + tolerance) return false;
        }
        return true;
    }
};

/// sRGB with the display transfer function taken off, which is the only form in
/// which the channels can be mixed arithmetically.
const Linear = struct {
    r: f32,
    g: f32,
    b: f32,

    fn toLab(self: Linear) Lab {
        const l = 0.4122214708 * self.r + 0.5363325363 * self.g + 0.0514459929 * self.b;
        const m = 0.2119034982 * self.r + 0.6806995451 * self.g + 0.1073969566 * self.b;
        const s = 0.0883024619 * self.r + 0.2817188376 * self.g + 0.6299787005 * self.b;

        const l_ = std.math.cbrt(l);
        const m_ = std.math.cbrt(m);
        const s_ = std.math.cbrt(s);

        return .{
            .l = 0.2104542553 * l_ + 0.7936177850 * m_ - 0.0040720468 * s_,
            .a = 1.9779984951 * l_ - 2.4285922050 * m_ + 0.4505937099 * s_,
            .b = 0.0259040371 * l_ + 0.7827717662 * m_ - 0.8086757660 * s_,
        };
    }
};

fn decode(channel: u8) f32 {
    const v = @as(f32, @floatFromInt(channel)) / 255.0;
    if (v <= 0.04045) return v / 12.92;
    return std.math.pow(f32, (v + 0.055) / 1.055, 2.4);
}

fn encode(channel: f32) u8 {
    const v = std.math.clamp(channel, 0, 1);
    const encoded = if (v <= 0.0031308)
        v * 12.92
    else
        1.055 * std.math.pow(f32, v, 1.0 / 2.4) - 0.055;
    return @intFromFloat(@round(std.math.clamp(encoded, 0, 1) * 255.0));
}

// -------------------------------------------------------------------------

test "a colour survives the round trip" {
    // Every hue family, plus the two ends and a grey, since those are the cases
    // where the hue is undefined and the transfer function is at its steepest.
    const cases = [_][]const u8{
        "#000000", "#ffffff", "#808080", "#3a76f0", "#d2691e",
        "#12805a", "#f2717b", "#070c15", "#dde6f4", "#8a5a00",
    };

    for (cases) |text| {
        const rgb = Rgb.parse(text).?;
        const back = toRgb(fromRgb(rgb));

        // One 8-bit step of slack: the trip is through a cube root and back.
        var buf: [7:0]u8 = undefined;
        const drift = @max(
            @abs(@as(i16, back.r) - @as(i16, rgb.r)),
            @max(
                @abs(@as(i16, back.g) - @as(i16, rgb.g)),
                @abs(@as(i16, back.b) - @as(i16, rgb.b)),
            ),
        );
        if (drift > 1) {
            std.debug.print("{s} came back as {s}\n", .{ text, back.hex(&buf) });
            return error.RoundTripDrifted;
        }
    }
}

test "black and white are where they should be" {
    try std.testing.expectApproxEqAbs(@as(f32, 0), fromRgb(.{ .r = 0, .g = 0, .b = 0 }).l, 1e-4);
    try std.testing.expectApproxEqAbs(@as(f32, 1), fromRgb(.{ .r = 255, .g = 255, .b = 255 }).l, 1e-4);

    // A grey has no hue, and says so rather than saying whatever `atan2` makes
    // of two zeroes.
    const grey = fromRgb(.{ .r = 128, .g = 128, .b = 128 });
    try std.testing.expectApproxEqAbs(@as(f32, 0), grey.c, 1e-3);
    try std.testing.expectEqual(@as(f32, 0), grey.h);
}

test "an impossible colour gives up chroma, not hue" {
    // A fully saturated yellow at the lightness of a midtone: nothing like it
    // exists in sRGB, and asking for it has to produce *something*.
    const asked: Lch = .{ .l = 0.5, .c = 0.3, .h = 100 };
    const got = fromRgb(toRgb(asked));

    try std.testing.expectApproxEqAbs(asked.l, got.l, 0.01);
    try std.testing.expectApproxEqAbs(asked.h, got.h, 1.5);
    try std.testing.expect(got.c < asked.c);

    // And it is still as colourful as sRGB allows, rather than having been
    // rounded down to something safe: pushing it back up leaves the gamut.
    try std.testing.expect(!Lab.at(asked.l, got.c + 0.01, std.math.degreesToRadians(asked.h)).inGamut());
}

test "lightness is even enough to build a ramp on" {
    // The point of the whole module: equal steps in `l` have to look like equal
    // steps, whatever the hue. What is checked here is the weaker property that
    // makes that possible — the steps come back out the size they went in, for
    // hues whose RGB representations are nothing alike.
    for ([_]f32{ 30, 100, 145, 260, 330 }) |hue| {
        var previous: f32 = 0;
        var step: f32 = 0.2;
        while (step <= 0.8) : (step += 0.2) {
            const back = fromRgb(toRgb(.{ .l = step, .c = 0.05, .h = hue })).l;
            try std.testing.expectApproxEqAbs(step, back, 0.01);
            try std.testing.expect(back > previous);
            previous = back;
        }
    }
}