type Vec3 = [f64; 3];
type Matrix3 = [Vec3; 3];
fn multiply(m: &Matrix3, v: Vec3) -> Vec3 {
let [x, y, z] = v;
[
m[0][0] * x + m[0][1] * y + m[0][2] * z,
m[1][0] * x + m[1][1] * y + m[1][2] * z,
m[2][0] * x + m[2][1] * y + m[2][2] * z,
]
}
const LAB_TO_LMS: Matrix3 = [
[1.0, 0.3963377773761749, 0.2158037573099136],
[1.0, -0.1055613458156586, -0.0638541728258133],
[1.0, -0.0894841775298119, -1.2914855480194092],
];
const LMS_TO_LINEAR_SRGB: Matrix3 = [
[4.0767416360759583, -3.3077115392580629, 0.2309699031821043],
[-1.2684379732850315, 2.6097573492876882, -0.341319376002657],
[-0.0041960761386756, -0.7034186179359362, 1.7076146940746117],
];
const SATURATION_FIT: [([f64; 2], [f64; 5]); 3] = [
(
[-1.8817031, -0.80936501],
[1.19086277, 1.76576728, 0.59662641, 0.75515197, 0.56771245],
),
(
[1.8144408, -1.19445267],
[0.73956515, -0.45954404, 0.08285427, 0.12541073, -0.14503204],
),
(
[0.13110758, 1.81333971],
[1.35733652, -0.00915799, -1.1513021, -0.50559606, 0.00692167],
),
];
const K1: f64 = 0.206;
const K2: f64 = 0.03;
const K3: f64 = (1.0 + K1) / (1.0 + K2);
fn okhsl_to_oklab_lightness(x: f64) -> f64 {
(x * x + K1 * x) / (K3 * (x + K2))
}
fn linear_to_srgb(value: f64) -> f64 {
if value > 0.0031308 {
1.055 * value.powf(1.0 / 2.4) - 0.055
} else {
12.92 * value
}
}
fn oklab_to_linear_srgb(lab: Vec3) -> Vec3 {
let lms = multiply(&LAB_TO_LMS, lab);
let cubed = [lms[0].powi(3), lms[1].powi(3), lms[2].powi(3)];
multiply(&LMS_TO_LINEAR_SRGB, cubed)
}
fn linear_srgb_to_rgb(linear: Vec3) -> [u8; 3] {
linear.map(|value| {
let encoded = linear_to_srgb(value.clamp(0.0, 1.0));
(encoded * 255.0).round() as u8
})
}
fn lms_slopes(a: f64, b: f64) -> Vec3 {
LAB_TO_LMS.map(|row| row[1] * a + row[2] * b)
}
fn max_saturation(a: f64, b: f64) -> f64 {
let channel = SATURATION_FIT
.iter()
.position(|([x, y], _)| x * a + y * b > 1.0)
.unwrap_or(2);
let [k0, k1, k2, k3, k4] = SATURATION_FIT[channel].1;
let weights = LMS_TO_LINEAR_SRGB[channel];
let saturation = k0 + k1 * a + k2 * b + k3 * a * a + k4 * a * b;
let slopes = lms_slopes(a, b);
let base = slopes.map(|k| 1.0 + saturation * k);
let dot =
|values: Vec3| weights[0] * values[0] + weights[1] * values[1] + weights[2] * values[2];
let f = dot(base.map(|value| value.powi(3)));
let f1 = dot(std::array::from_fn(|i| 3.0 * slopes[i] * base[i] * base[i]));
let f2 = dot(std::array::from_fn(|i| {
6.0 * slopes[i] * slopes[i] * base[i]
}));
saturation - (f * f1) / (f1 * f1 - 0.5 * f * f2)
}
fn cusp(a: f64, b: f64) -> (f64, f64) {
let saturation = max_saturation(a, b);
let max_linear = oklab_to_linear_srgb([1.0, saturation * a, saturation * b])
.into_iter()
.fold(f64::MIN, f64::max);
let lightness = (1.0 / max_linear).cbrt();
(lightness, lightness * saturation)
}
fn max_chroma(a: f64, b: f64, lightness: f64, cusp: (f64, f64)) -> f64 {
let (cusp_l, cusp_c) = cusp;
if lightness <= cusp_l {
return (cusp_c * lightness) / cusp_l;
}
let t = (cusp_c * (lightness - 1.0)) / (cusp_l - 1.0);
let slopes = lms_slopes(a, b);
let lms = slopes.map(|k| lightness + t * k);
let cubes = lms.map(|value| value.powi(3));
let first = std::array::from_fn(|i| 3.0 * slopes[i] * lms[i] * lms[i]);
let second = std::array::from_fn(|i| 6.0 * slopes[i] * slopes[i] * lms[i]);
let dot =
|row: Vec3, values: Vec3| row[0] * values[0] + row[1] * values[1] + row[2] * values[2];
let steps: [f64; 3] = LMS_TO_LINEAR_SRGB.map(|row| {
let f = dot(row, cubes) - 1.0;
let f1 = dot(row, first);
let f2 = dot(row, second);
let u = f1 / (f1 * f1 - 0.5 * f * f2);
if u >= 0.0 {
-f * u
} else {
f64::MAX
}
});
t + steps.into_iter().fold(f64::MAX, f64::min)
}
fn chroma_stops(l: f64, a: f64, b: f64) -> (f64, f64, f64) {
let peak = cusp(a, b);
let c_max = max_chroma(a, b, l, peak);
let k = c_max / (l * (peak.1 / peak.0)).min((1.0 - l) * (peak.1 / (1.0 - peak.0)));
let mid_s = 0.11516993
+ 1.0
/ (7.4477897
+ 4.1590124 * b
+ a * (-2.19557347
+ 1.75198401 * b
+ a * (-2.13704948 - 10.02301043 * b
+ a * (-4.24894561 + 5.38770819 * b + 4.69891013 * a))));
let mid_t = 0.11239642
+ 1.0
/ (1.6132032 - 0.68124379 * b
+ a * (0.40370612
+ 0.90148123 * b
+ a * (-0.27087943
+ 0.6122399 * b
+ a * (0.00299215 - 0.45399568 * b - 0.14661872 * a))));
let c_mid = 0.9
* k
* ((1.0 / (1.0 / (l * mid_s).powi(4) + 1.0 / ((1.0 - l) * mid_t).powi(4)))
.sqrt()
.sqrt());
let c0 = (1.0 / (1.0 / (l * 0.4).powi(2) + 1.0 / ((1.0 - l) * 0.8).powi(2))).sqrt();
(c0, c_mid, c_max)
}
pub fn okhsl_to_rgb(hue: f64, saturation: f64, lightness: f64) -> [u8; 3] {
let l = okhsl_to_oklab_lightness(lightness);
let mut lab: Vec3 = [l, 0.0, 0.0];
if l > 0.0 && l < 1.0 && saturation > 0.0 {
let angle = 2.0 * std::f64::consts::PI * (hue.rem_euclid(360.0)) / 360.0;
let a = angle.cos();
let b = angle.sin();
let (c0, c_mid, c_max) = chroma_stops(l, a, b);
let chroma = if saturation < 0.8 {
let t = 1.25 * saturation;
let k1 = 0.8 * c0;
(t * k1) / (1.0 - (1.0 - k1 / c_mid) * t)
} else {
let t = 5.0 * (saturation - 0.8);
let k1 = (0.2 * c_mid * c_mid * 1.25 * 1.25) / c0;
c_mid + (t * k1) / (1.0 - (1.0 - k1 / (c_max - c_mid)) * t)
};
lab = [l, chroma * a, chroma * b];
}
linear_srgb_to_rgb(oklab_to_linear_srgb(lab))
}
pub fn oklch_to_rgb(l: f64, c: f64, h: f64) -> [u8; 3] {
let radians = h * std::f64::consts::PI / 180.0;
let cos = radians.cos();
let sin = radians.sin();
let at_chroma = |chroma: f64| oklab_to_linear_srgb([l, chroma * cos, chroma * sin]);
let direct = at_chroma(c);
if direct.iter().all(|&v| v >= -1e-7 && v <= 1.0 + 1e-7) {
return linear_srgb_to_rgb(direct);
}
let mut linear = at_chroma(0.0);
let mut low = 0.0;
let mut high = c;
for _ in 0..20 {
let chroma = (low + high) / 2.0;
let candidate = at_chroma(chroma);
if candidate.iter().all(|&v| v >= -1e-7 && v <= 1.0 + 1e-7) {
low = chroma;
linear = candidate;
} else {
high = chroma;
}
}
linear_srgb_to_rgb(linear)
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_okhsl_spread_exercises_all_branches_without_panic() {
for hue in [0.0, 37.5, 83.0, 159.0, 232.0, 295.0, 360.0] {
for sat in [0.1, 0.5, 0.8, 1.0] {
for light in [0.1, 0.3, 0.5, 0.7, 0.9] {
let rgb = okhsl_to_rgb(hue, sat, light);
let _ = rgb;
}
}
}
}
#[test]
fn test_oklch_achromatic_fallback() {
assert_eq!(oklch_to_rgb(1.0, 0.3, 150.0), [255, 255, 255]);
}
#[test]
fn test_okhsl_gray_degenerates_to_lightness() {
assert_eq!(okhsl_to_rgb(30.0, 0.0, 0.5), okhsl_to_rgb(300.0, 0.0, 0.5));
}
#[test]
fn test_golden_okhsl_oklch_channels() {
assert_eq!(okhsl_to_rgb(295.0, 0.5, 0.67), [167, 152, 215]); assert_eq!(okhsl_to_rgb(231.0, 0.68, 0.47), [47, 120, 153]); assert_eq!(oklch_to_rgb(0.6, 0.15, 250.0), [39, 132, 213]); }
}