pdf-odin-llm/colorspace.odin

package pdf

import "core:math"

// Colour spaces (ISO 32000-2, 8.6).
//
// Every space is reduced to RGB for output. That is not colour management and
// does not pretend to be: an ICCBased space renders through its /Alternate (or
// through the device space its /N implies), which is what every
// non-colour-managed viewer does and what the format explicitly permits.

Colorspace_Kind :: enum u8 {
    Device_Gray,
    Device_RGB,
    Device_CMYK,
    Cal_Gray,
    Cal_RGB,
    Lab,
    ICC_Based,
    Indexed,
    Separation,
    Device_N,
    Pattern,
}

Colorspace :: struct {
    kind: Colorspace_Kind,
    // Components a colour value in this space carries. 1 for Indexed (the
    // index itself), not the base space's count.
    n:    int,
    name: string,

    // Indexed: `base` is the space the lookup table's entries live in, and
    // `lookup` holds (hival+1) * base.n bytes.
    base:   ^Colorspace,
    lookup: []byte,
    hival:  int,

    // Separation, DeviceN and ICCBased all render through `alt`. For the first
    // two, `tint` maps the colorant values into it.
    alt:            ^Colorspace,
    tint:           ^Function,
    colorant_names: []string,
    // The /None colorant marks nothing at all on the page (8.6.6.4).
    is_none:        bool,
    // The /All colorant marks every separation.
    is_all:         bool,

    // CalGray, CalRGB and Lab.
    white_point: [3]f64,
    black_point: [3]f64,
    gamma:       [3]f64,
    // CalRGB's /Matrix. For a Lab space this holds something else the caller
    // never sees: the 3x3 that takes XYZ under this space's own white point to
    // linear sRGB, chromatic adaptation folded in, computed once at load time
    // so that converting a colour is nine multiplies and no matrix building.
    mtx:         [9]f64,
    // Lab's a and b component ranges, defaulting to [-100 100 -100 100].
    range:       [4]f64,

    // Pattern's underlying space, present only for an uncoloured tiling
    // pattern, whose colour is supplied by the content that paints it.
    under: ^Colorspace,
}

// Nesting of colour spaces: an Indexed whose base is a Separation whose
// alternate is an ICCBased, and so on. Four is the deepest anything
// meaningful goes; the limit exists because /Resources can be made to name a
// space that names itself.
@(private)
MAX_COLORSPACE_DEPTH :: 8

// The most entries an Indexed lookup table can address. /HiVal is capped at
// 255 by the format (8.6.6.3), which is what makes an index a single byte.
@(private)
MAX_INDEXED_HIVAL :: 255

// ---------------------------------------------------------------------------
// The device spaces
//
// These are shared, immutable and outlive every Document: they hold no
// document memory and nothing ever writes to them, so handing the same pointer
// to two Documents (or to two threads rendering different files) is safe. It
// also means the overwhelmingly common `/DeviceRGB cs` costs nothing at all.
// ---------------------------------------------------------------------------

@(private)
cs_gray := Colorspace{kind = .Device_Gray, n = 1, name = "DeviceGray"}
@(private)
cs_rgb := Colorspace{kind = .Device_RGB, n = 3, name = "DeviceRGB"}
@(private)
cs_cmyk := Colorspace{kind = .Device_CMYK, n = 4, name = "DeviceCMYK"}
// A /Pattern space with no underlying space: the pattern carries its own
// colour, so a colour value in it has no components at all. `n` is therefore
// 0, not 1 -- `n` is what tells the content interpreter how many numeric
// operands `scn` takes, and a coloured pattern is set with `/P1 scn`, a name
// and nothing else. An uncoloured pattern is the case with components, and it
// gets them from its underlying space in cs_load_pattern.
@(private)
cs_pattern := Colorspace{kind = .Pattern, n = 0, name = "Pattern"}

colorspace_device_gray :: proc(doc: ^Document) -> ^Colorspace {
    return &cs_gray
}

colorspace_device_rgb :: proc(doc: ^Document) -> ^Colorspace {
    return &cs_rgb
}

colorspace_device_cmyk :: proc(doc: ^Document) -> ^Colorspace {
    return &cs_cmyk
}

// colorspace_device_n returns the device space with `n` components, which is
// how an ICCBased space with no /Alternate and an inline image with no space
// at all are resolved (8.6.5.5).
@(private)
colorspace_device_n :: proc(n: int) -> ^Colorspace {
    switch n {
    case 1: return &cs_gray
    case 4: return &cs_cmyk
    }
    return &cs_rgb
}

// ---------------------------------------------------------------------------
// Loading
// ---------------------------------------------------------------------------

// colorspace_load builds the colour space `obj` names or describes.
//
// `res` is the /Resources dictionary of the page or form whose content stream
// named it, and is where a bare name that is not a device space is looked up
// (through /Resources /ColorSpace); it may be nil, and a caller that has only
// the /ColorSpace sub-dictionary may pass that instead.
//
// The result is owned by the Document and must not be freed. Spaces are cached
// by the object they were built from, because a page that names the same space
// in four hundred XObjects must build it once.
colorspace_load :: proc(doc: ^Document, obj: Object, res: ^Dict) -> (^Colorspace, Error) {
    return cs_load(doc, obj, res, 0)
}

@(private)
cs_load :: proc(doc: ^Document, obj: Object, res: ^Dict, depth: int) -> (^Colorspace, Error) {
    if doc == nil {
        return nil, .Invalid_Argument
    }
    if depth > MAX_COLORSPACE_DEPTH {
        add_repair(doc, .Recursion_Limit_Hit, 0)
        return nil, .Recursion_Limit
    }

    o := resolve(doc, obj)
    key := object_identity(o)
    if key != 0 {
        if cached, ok := doc.cs_cache[key]; ok {
            return cached, .None
        }
    }

    cs: ^Colorspace
    err := Error.Bad_Colorspace
    #partial switch v in o {
    case Name:
        cs, err = cs_load_name(doc, string(v), res, depth)
    case ^Array:
        cs, err = cs_load_array(doc, v, res, depth)
    }
    if err != .None || cs == nil {
        return nil, err == .None ? .Bad_Colorspace : err
    }

    if key != 0 {
        if doc.cs_cache == nil {
            doc.cs_cache = make(map[uintptr]^Colorspace, 16, doc.allocator)
        }
        doc.cs_cache[key] = cs
    }
    return cs, .None
}

@(private)
cs_load_name :: proc(doc: ^Document, name: string, res: ^Dict, depth: int) -> (^Colorspace, Error) {
    // The device spaces cannot be redefined by a resource dictionary
    // (8.6.3), so they are answered before it is consulted.
    switch name {
    case "DeviceGray": return &cs_gray, .None
    case "DeviceRGB":  return &cs_rgb, .None
    case "DeviceCMYK": return &cs_cmyk, .None
    case "Pattern":    return &cs_pattern, .None
    }

    if res != nil {
        // `res` is the /Resources dictionary, so a name is looked up in its
        // /ColorSpace sub-dictionary (8.6.3). A caller holding only that
        // sub-dictionary may pass it instead: it is tried second, so it cannot
        // shadow a genuine /ColorSpace entry.
        if csdict := dict_dict(doc, res, "ColorSpace"); csdict != nil {
            if entry := dict_get_raw(csdict, name); entry != nil {
                return cs_load(doc, entry, res, depth + 1)
            }
        }
        if entry := dict_get_raw(res, name); entry != nil {
            return cs_load(doc, entry, res, depth + 1)
        }
    }

    // The inline-image abbreviations (Table 93). They are only legal inside an
    // inline image dictionary, so they are tried after the resource lookup: a
    // resource genuinely called /G has to win.
    switch name {
    case "G":    return &cs_gray, .None
    case "RGB":  return &cs_rgb, .None
    case "CMYK": return &cs_cmyk, .None
    }
    return nil, .Bad_Colorspace
}

@(private)
cs_load_array :: proc(doc: ^Document, a: ^Array, res: ^Dict, depth: int) -> (^Colorspace, Error) {
    if array_len(a) == 0 {
        return nil, .Bad_Colorspace
    }
    family := array_name(doc, a, 0)
    // `[/DeviceRGB]` and `[/Pattern]` are both legal spellings.
    if array_len(a) == 1 {
        return cs_load_name(doc, family, res, depth)
    }

    switch family {
    case "DeviceGray", "G":    return &cs_gray, .None
    case "DeviceRGB", "RGB":   return &cs_rgb, .None
    case "DeviceCMYK", "CMYK": return &cs_cmyk, .None
    case "CalGray":            return cs_load_calgray(doc, a)
    case "CalRGB":             return cs_load_calrgb(doc, a)
    case "Lab":                return cs_load_lab(doc, a)
    case "ICCBased":           return cs_load_icc(doc, a, res, depth)
    case "Indexed", "I":       return cs_load_indexed(doc, a, res, depth)
    case "Separation":         return cs_load_separation(doc, a, res, depth)
    case "DeviceN":            return cs_load_devicen(doc, a, res, depth)
    case "Pattern":            return cs_load_pattern(doc, a, res, depth)
    }
    return nil, .Bad_Colorspace
}

@(private)
cs_new :: proc(doc: ^Document, kind: Colorspace_Kind, n: int, name: string) -> ^Colorspace {
    cs := new(Colorspace, doc.allocator)
    cs.kind = kind
    cs.n = n
    cs.name = name
    cs.white_point = {0.9642, 1.0, 0.8249}   // D50, the format's own default
    cs.gamma = {1, 1, 1}
    cs.range = {-100, 100, -100, 100}
    return cs
}

// cs_read_white_point reads /WhitePoint, which must be a positive XYZ triple
// with Y = 1. A file that gives something else gets D50, because a zero or
// negative component divides by zero further down.
@(private)
cs_read_white_point :: proc(doc: ^Document, d: ^Dict, cs: ^Colorspace) {
    wp := dict_array(doc, d, "WhitePoint")
    if array_len(wp) < 3 {
        return
    }
    x := clamp_finite(array_num(doc, wp, 0, 0), 0)
    y := clamp_finite(array_num(doc, wp, 1, 0), 0)
    z := clamp_finite(array_num(doc, wp, 2, 0), 0)
    if x <= 0 || y <= 0 || z <= 0 {
        add_repair(doc, .Bad_Resource, 0)
        return
    }
    cs.white_point = {x, y, z}

    bp := dict_array(doc, d, "BlackPoint")
    if array_len(bp) >= 3 {
        cs.black_point = {
            clamp_finite(array_num(doc, bp, 0, 0), 0),
            clamp_finite(array_num(doc, bp, 1, 0), 0),
            clamp_finite(array_num(doc, bp, 2, 0), 0),
        }
    }
}

@(private)
cs_load_calgray :: proc(doc: ^Document, a: ^Array) -> (^Colorspace, Error) {
    d := array_dict(doc, a, 1)
    cs := cs_new(doc, .Cal_Gray, 1, "CalGray")
    cs_read_white_point(doc, d, cs)
    g := clamp_finite(dict_num(doc, d, "Gamma", 1), 1)
    if g <= 0 {
        g = 1
    }
    cs.gamma = {g, g, g}
    return cs, .None
}

@(private)
cs_load_calrgb :: proc(doc: ^Document, a: ^Array) -> (^Colorspace, Error) {
    d := array_dict(doc, a, 1)
    cs := cs_new(doc, .Cal_RGB, 3, "CalRGB")
    cs_read_white_point(doc, d, cs)
    if ga := dict_array(doc, d, "Gamma"); array_len(ga) >= 3 {
        for i in 0 ..< 3 {
            g := clamp_finite(array_num(doc, ga, i, 1), 1)
            cs.gamma[i] = g > 0 ? g : 1
        }
    }
    if ma := dict_array(doc, d, "Matrix"); array_len(ma) >= 9 {
        for i in 0 ..< 9 {
            cs.mtx[i] = clamp_finite(array_num(doc, ma, i, 0), 0)
        }
    }
    return cs, .None
}

@(private)
cs_load_lab :: proc(doc: ^Document, a: ^Array) -> (^Colorspace, Error) {
    d := array_dict(doc, a, 1)
    cs := cs_new(doc, .Lab, 3, "Lab")
    cs_read_white_point(doc, d, cs)
    if r := dict_array(doc, d, "Range"); array_len(r) >= 4 {
        for i in 0 ..< 4 {
            cs.range[i] = clamp_finite(array_num(doc, r, i, 0), 0)
        }
        // An inverted range would clamp every value to nothing.
        if cs.range[0] > cs.range[1] {
            cs.range[0], cs.range[1] = cs.range[1], cs.range[0]
        }
        if cs.range[2] > cs.range[3] {
            cs.range[2], cs.range[3] = cs.range[3], cs.range[2]
        }
    }
    cs.mtx = lab_to_srgb_matrix(cs.white_point)
    return cs, .None
}

@(private)
cs_load_icc :: proc(doc: ^Document, a: ^Array, res: ^Dict, depth: int) -> (^Colorspace, Error) {
    stm := to_stream(array_get(doc, a, 1)) or_else nil
    d, _ := to_dict(array_get(doc, a, 1))
    n := int(dict_int(doc, d, "N", 0))
    if n != 1 && n != 3 && n != 4 {
        n = 0
    }

    alt: ^Colorspace
    if entry := dict_get_raw(d, "Alternate"); entry != nil {
        loaded, err := cs_load(doc, entry, res, depth + 1)
        if err == .None && loaded != nil && loaded.kind != .Pattern {
            // An /Alternate that disagrees with /N is the file contradicting
            // itself; /N is what the sample data is laid out by.
            if n == 0 || loaded.n == n {
                alt = loaded
            }
        }
    }
    if alt == nil {
        if n == 0 {
            // Neither /N nor a usable /Alternate. The profile itself would say,
            // but reading ICC profiles is exactly what this library does not
            // do, so the commonest case is assumed.
            add_repair(doc, .Bad_Resource, 0, stm != nil ? stm.ref : Ref{})
            n = 3
        }
        alt = colorspace_device_n(n)
    }

    // Rendering happens entirely through the alternate: no ICC profile is ever
    // parsed and no colour management is performed. Every non-colour-managed
    // viewer does the same, and 8.6.5.5 allows it.
    cs := cs_new(doc, .ICC_Based, alt.n, "ICCBased")
    cs.alt = alt
    return cs, .None
}

@(private)
cs_load_indexed :: proc(doc: ^Document, a: ^Array, res: ^Dict, depth: int) -> (^Colorspace, Error) {
    base, berr := cs_load(doc, array_get_raw(a, 1), res, depth + 1)
    if berr != .None || base == nil {
        return nil, .Bad_Colorspace
    }
    // The base of an Indexed space may not itself be Indexed or Pattern
    // (8.6.6.3); either would make an index mean nothing.
    if base.kind == .Indexed || base.kind == .Pattern {
        return nil, .Bad_Colorspace
    }

    cs := cs_new(doc, .Indexed, 1, "Indexed")
    cs.base = base
    cs.hival = int(clamp(array_int(doc, a, 2, 0), 0, MAX_INDEXED_HIVAL))

    lookup := array_get(doc, a, 3)
    #partial switch v in lookup {
    case String:
        cs.lookup = transmute([]byte)v
    case ^Stream:
        data, derr := stream_data(doc, v)
        if derr == .None {
            cs.lookup = data
        }
    }
    if len(cs.lookup) < (cs.hival + 1) * base.n {
        // A short table is common in files that were subset or edited. The
        // entries past its end read as zero rather than past the buffer.
        add_repair(doc, .Bad_Resource, 0)
    }
    return cs, .None
}

@(private)
cs_load_separation :: proc(doc: ^Document, a: ^Array, res: ^Dict, depth: int) -> (^Colorspace, Error) {
    cs := cs_new(doc, .Separation, 1, "Separation")
    colorant := array_name(doc, a, 1)
    names := make([]string, 1, doc.allocator)
    names[0] = colorant
    cs.colorant_names = names
    // /None marks nothing on the page at all; /All marks every separation
    // (8.6.6.4).
    cs.is_none = colorant == "None"
    cs.is_all = colorant == "All"

    cs_load_tint(doc, cs, a, res, depth)
    return cs, .None
}

@(private)
cs_load_devicen :: proc(doc: ^Document, a: ^Array, res: ^Dict, depth: int) -> (^Colorspace, Error) {
    list := to_array(array_get(doc, a, 1)) or_else nil
    n := array_len(list)
    if n < 1 {
        return nil, .Bad_Colorspace
    }
    // 8.6.6.5 caps a DeviceN space at 32 colorants, which is also the width of
    // a Color.
    n = min(n, MAX_COLOR_COMPONENTS)

    cs := cs_new(doc, .Device_N, n, "DeviceN")
    names := make([]string, n, doc.allocator)
    none_count := 0
    for i in 0 ..< n {
        names[i] = array_name(doc, list, i)
        if names[i] == "None" {
            none_count += 1
        }
    }
    cs.colorant_names = names
    // A DeviceN marks nothing only when every one of its colorants is /None.
    cs.is_none = none_count == n

    cs_load_tint(doc, cs, a, res, depth)
    return cs, .None
}

// cs_load_tint fills in the alternate space and tint transform shared by
// Separation and DeviceN. Both may be missing or broken, and neither failure
// is a reason to reject the space: `colorspace_to_rgb` falls back to treating
// the colorant as ink coverage, which keeps a spot colour dark instead of
// dropping the artwork that uses it.
@(private)
cs_load_tint :: proc(doc: ^Document, cs: ^Colorspace, a: ^Array, res: ^Dict, depth: int) {
    // 8.6.6.4 and 8.6.6.5: the alternate space "shall be any device or
    // CIE-based colour space but shall not be another special colour space
    // (Pattern, Indexed, Separation, or DeviceN)". All four are refused, and
    // for the same reason -- none of them can be handed the output of a tint
    // transform and give a colour back. The components of a tint transform's
    // output run 0..1; an Indexed component is a palette index, so following
    // one collapses the whole ramp onto two entries, and a Pattern has no
    // colour value at all. pdfium refuses the whole space at this point;
    // dropping to the ink-coverage fallback below is gentler and keeps the
    // artwork on the page.
    refused := false
    if alt, err := cs_load(doc, array_get_raw(a, 2), res, depth + 1); err == .None && alt != nil {
        #partial switch alt.kind {
        case .Pattern, .Indexed, .Separation, .Device_N:
            refused = true
        case:
            cs.alt = alt
        }
    }
    if tint, err := function_load(doc, array_get_raw(a, 3)); err == .None && tint != nil && tint.n_out >= 1 {
        if cs.alt == nil && !refused {
            // A tint transform with no alternate space still says how many
            // components it produces, and a device space with that many is a
            // far better guess than giving up on the colour entirely. That
            // does not extend to an alternate that was named and refused: its
            // outputs mean whatever that space's components mean, which is
            // exactly what cannot be read as a device colour.
            cs.alt = colorspace_device_n(tint.n_out)
        }
        if cs.alt != nil && tint.n_out == cs.alt.n {
            cs.tint = tint
        }
    }
    if cs.tint == nil || cs.alt == nil {
        add_repair(doc, .Bad_Resource, 0)
    }
}

@(private)
cs_load_pattern :: proc(doc: ^Document, a: ^Array, res: ^Dict, depth: int) -> (^Colorspace, Error) {
    under, err := cs_load(doc, array_get_raw(a, 1), res, depth + 1)
    if err != .None || under == nil || under.kind == .Pattern {
        return &cs_pattern, .None
    }
    cs := cs_new(doc, .Pattern, under.n, "Pattern")
    cs.under = under
    return cs, .None
}

// ---------------------------------------------------------------------------
// Conversion to RGB
// ---------------------------------------------------------------------------

// colorspace_to_rgb converts one colour value to RGB in 0..1. It allocates
// nothing and never fails: an unconvertible colour comes back as black, since
// a page that draws the wrong colour is still a page and a page that stops
// drawing is not.
colorspace_to_rgb :: proc(cs: ^Colorspace, c: Color) -> [3]f64 {
    return cs_to_rgb(cs, c, 0)
}

@(private)
cs_to_rgb :: proc(cs: ^Colorspace, c: Color, depth: int) -> [3]f64 {
    if cs == nil || depth > MAX_COLORSPACE_DEPTH {
        return {0, 0, 0}
    }

    switch cs.kind {
    case .Device_Gray:
        g := cs_comp(c, 0)
        return {g, g, g}

    case .Device_RGB:
        return {cs_comp(c, 0), cs_comp(c, 1), cs_comp(c, 2)}

    case .Device_CMYK:
        // The naive conversion, subtracting each ink and the black from full
        // intensity. It is not colorimetric and no non-colour-managed viewer
        // does anything else; a real CMYK profile is the only way to do
        // better and this library reads none.
        cy, m, y, k := cs_comp(c, 0), cs_comp(c, 1), cs_comp(c, 2), cs_comp(c, 3)
        return {
            1 - min(1, cy + k),
            1 - min(1, m + k),
            1 - min(1, y + k),
        }

    case .Cal_Gray:
        g := cs_cal_component(cs_comp(c, 0), cs.gamma[0])
        return {g, g, g}

    case .Cal_RGB:
        return {
            cs_cal_component(cs_comp(c, 0), cs.gamma[0]),
            cs_cal_component(cs_comp(c, 1), cs.gamma[1]),
            cs_cal_component(cs_comp(c, 2), cs.gamma[2]),
        }

    case .Lab:
        return cs_lab_to_rgb(cs, c)

    case .ICC_Based:
        return cs_to_rgb(cs.alt, c, depth + 1)

    case .Indexed:
        return cs_indexed_to_rgb(cs, c, depth)

    case .Separation, .Device_N:
        return cs_tint_to_rgb(cs, c, depth)

    case .Pattern:
        // An uncoloured pattern carries the colour it is painted in, in the
        // underlying space. A coloured one supplies its own, and this value is
        // never asked for.
        if cs.under != nil {
            return cs_to_rgb(cs.under, c, depth + 1)
        }
        return {0, 0, 0}
    }
    return {0, 0, 0}
}

// cs_comp reads component i, clamped to 0..1. Used for the spaces whose
// components are all 0..1, which is every one of them but Lab and Indexed.
@(private)
cs_comp :: proc(c: Color, i: int) -> f64 {
    if i < 0 || i >= c.n || i >= MAX_COLOR_COMPONENTS {
        return 0
    }
    return clamp(clamp_finite(c.v[i], 0), 0, 1)
}

// cs_cal_component decodes one CalGray or CalRGB component.
//
// The full colorimetry -- white point, /Matrix, adaptation -- is deliberately
// not implemented: the format allows a reader to treat these as DeviceGray and
// DeviceRGB and viewers do. /Gamma is honoured, because it is one pow: the
// component is decoded to linear light and re-encoded for an sRGB output
// device. At the usual /Gamma 2.2 that round trip is within half a level of
// leaving the value alone, so the common case still matches every other
// viewer; a gamma of exactly 1 is the default and is skipped entirely.
@(private)
cs_cal_component :: proc(v: f64, gamma: f64) -> f64 {
    x := clamp(clamp_finite(v, 0), 0, 1)
    if gamma == 1 {
        return x
    }
    return cs_srgb_encode(math.pow(x, gamma))
}

// cs_srgb_encode is the sRGB transfer function (IEC 61966-2-1), taking linear
// light to a display value.
@(private)
cs_srgb_encode :: proc(v: f64) -> f64 {
    x := clamp(clamp_finite(v, 0), 0, 1)
    if x <= 0.0031308 {
        return x * 12.92
    }
    return 1.055 * math.pow(x, 1.0 / 2.4) - 0.055
}

// ---------------------------------------------------------------------------
// Lab
// ---------------------------------------------------------------------------

// lab_to_srgb_matrix builds XYZ-under-`wp` to linear sRGB.
//
// sRGB is defined against D65 and PDF's Lab spaces are almost always D50, so
// the two are joined by a Bradford chromatic adaptation. Without it a D50
// white -- L* 100, the commonest colour in any Lab image -- renders as a
// yellow cast rather than as white.
@(private)
lab_to_srgb_matrix :: proc(wp: [3]f64) -> [9]f64 {
    // Bradford cone response and its inverse.
    ma := [9]f64{
         0.8951,  0.2664, -0.1614,
        -0.7502,  1.7135,  0.0367,
         0.0389, -0.0685,  1.0296,
    }
    ma_inv := [9]f64{
         0.9869929, -0.1470543, 0.1599627,
         0.4323053,  0.5183603, 0.0492912,
        -0.0085287,  0.0400428, 0.9684867,
    }
    // XYZ (D65) to linear sRGB.
    to_rgb := [9]f64{
         3.2404542, -1.5371385, -0.4985314,
        -0.9692660,  1.8760108,  0.0415560,
         0.0556434, -0.2040259,  1.0572252,
    }
    d65 := [3]f64{0.95047, 1.0, 1.08883}

    src := mat3_mul_vec(ma, wp)
    dst := mat3_mul_vec(ma, d65)
    scale: [9]f64
    for i in 0 ..< 3 {
        // A white point component of zero was rejected at load time; this is
        // belt and braces, since a division by zero here would poison every
        // colour the space produces.
        scale[i * 3 + i] = src[i] != 0 ? dst[i] / src[i] : 1
    }
    adapt := mat3_mul(ma_inv, mat3_mul(scale, ma))
    return mat3_mul(to_rgb, adapt)
}

@(private)
mat3_mul :: proc(a, b: [9]f64) -> (m: [9]f64) {
    for r in 0 ..< 3 {
        for c in 0 ..< 3 {
            sum := 0.0
            for k in 0 ..< 3 {
                sum += a[r * 3 + k] * b[k * 3 + c]
            }
            m[r * 3 + c] = sum
        }
    }
    return
}

@(private)
mat3_mul_vec :: proc(a: [9]f64, v: [3]f64) -> (o: [3]f64) {
    for r in 0 ..< 3 {
        o[r] = a[r * 3] * v[0] + a[r * 3 + 1] * v[1] + a[r * 3 + 2] * v[2]
    }
    return
}

// cs_lab_g inverts the CIE L* nonlinearity.
@(private)
cs_lab_g :: proc(x: f64) -> f64 {
    // 6/29, below which the curve is linear so that its slope stays finite.
    if x >= 6.0 / 29.0 {
        return x * x * x
    }
    return (108.0 / 841.0) * (x - 4.0 / 29.0)
}

@(private)
cs_lab_to_rgb :: proc(cs: ^Colorspace, c: Color) -> [3]f64 {
    l := clamp(clamp_finite(c.n > 0 ? c.v[0] : 0, 0), 0, 100)
    a := clamp(clamp_finite(c.n > 1 ? c.v[1] : 0, 0), cs.range[0], cs.range[1])
    b := clamp(clamp_finite(c.n > 2 ? c.v[2] : 0, 0), cs.range[2], cs.range[3])

    m := (l + 16) / 116
    xw := cs.white_point[0] * cs_lab_g(m + a / 500)
    yw := cs.white_point[1] * cs_lab_g(m)
    zw := cs.white_point[2] * cs_lab_g(m - b / 200)

    out: [3]f64
    for i in 0 ..< 3 {
        lin := cs.mtx[i * 3] * xw + cs.mtx[i * 3 + 1] * yw + cs.mtx[i * 3 + 2] * zw
        // Colours outside the sRGB gamut -- most of the Lab space is -- are
        // clipped per channel, which is what every viewer without a gamut
        // mapping engine does.
        out[i] = cs_srgb_encode(lin)
    }
    return out
}

// ---------------------------------------------------------------------------
// Indexed
// ---------------------------------------------------------------------------

@(private)
cs_indexed_to_rgb :: proc(cs: ^Colorspace, c: Color, depth: int) -> [3]f64 {
    base := cs.base
    if base == nil {
        return {0, 0, 0}
    }
    // An index out of range is clamped, never followed: the table is exactly
    // (hival+1) entries and a file is free to name entry 300 of a 4-entry one.
    idx := int(clamp(clamp_finite(math.round(c.n > 0 ? c.v[0] : 0), 0), 0, f64(cs.hival)))

    comp: Color
    comp.n = min(base.n, MAX_COLOR_COMPONENTS)
    for j in 0 ..< comp.n {
        off := idx * base.n + j
        raw := off >= 0 && off < len(cs.lookup) ? f64(cs.lookup[off]) : 0
        lo, hi := cs_component_range(base, j)
        comp.v[j] = lo + raw * (hi - lo) / 255
    }
    return cs_to_rgb(base, comp, depth + 1)
}

// cs_component_range is the interval component j of a colour in this space
// runs over. Everything is 0..1 except Lab, whose L is 0..100 and whose a and
// b carry their own /Range.
@(private)
cs_component_range :: proc(cs: ^Colorspace, j: int) -> (lo, hi: f64) {
    if cs != nil && cs.kind == .Lab {
        switch j {
        case 0: return 0, 100
        case 1: return cs.range[0], cs.range[1]
        case 2: return cs.range[2], cs.range[3]
        }
    }
    return 0, 1
}

// ---------------------------------------------------------------------------
// Separation and DeviceN
// ---------------------------------------------------------------------------

@(private)
cs_tint_to_rgb :: proc(cs: ^Colorspace, c: Color, depth: int) -> [3]f64 {
    if cs.is_none {
        // Nothing is painted in a /None space at all. The interpreter checks
        // `is_none` before drawing; this value only decides what a caller that
        // does not sees, and white is the least destructive answer.
        return {1, 1, 1}
    }

    if cs.tint != nil && cs.alt != nil {
        in_:  [MAX_COLOR_COMPONENTS]f64
        out_: [MAX_COLOR_COMPONENTS]f64
        n := min(cs.n, MAX_COLOR_COMPONENTS)
        for i in 0 ..< n {
            in_[i] = clamp_finite(i < c.n ? c.v[i] : 0, 0)
        }
        alt_n := min(cs.alt.n, MAX_COLOR_COMPONENTS)
        function_eval(cs.tint, in_[:n], out_[:alt_n])

        alt_color: Color
        alt_color.n = alt_n
        for i in 0 ..< alt_n {
            alt_color.v[i] = out_[i]
        }
        return cs_to_rgb(cs.alt, alt_color, depth + 1)
    }

    // No usable tint transform. Treating the darkest colorant as ink coverage
    // keeps a spot colour dark and a light tint light, which is much closer
    // than either black or white and never loses the artwork.
    ink := 0.0
    for i in 0 ..< min(cs.n, c.n) {
        ink = max(ink, clamp(clamp_finite(c.v[i], 0), 0, 1))
    }
    g := 1 - ink
    return {g, g, g}
}

// ---------------------------------------------------------------------------
// Initial values and image decode arrays
// ---------------------------------------------------------------------------

// colorspace_default_color is the colour a space starts at when it is selected
// with `cs` or `CS` (ISO 32000-2, 8.6.8): black for the device and CIE spaces,
// the first entry for Indexed, full tint for Separation and DeviceN.
colorspace_default_color :: proc(cs: ^Colorspace) -> (c: Color) {
    if cs == nil {
        return color_make(0)
    }
    n := min(cs.n, MAX_COLOR_COMPONENTS)
    switch cs.kind {
    case .Device_CMYK:
        // 0 0 0 1: black, not the white that four zero components would be.
        c.n = 4
        c.v[3] = 1
    case .Separation, .Device_N:
        c.n = n
        for i in 0 ..< n {
            c.v[i] = 1
        }
    case .Pattern:
        // Until a pattern is selected there is nothing to paint with. An
        // uncoloured pattern's components come from the underlying space.
        if cs.under != nil {
            return colorspace_default_color(cs.under)
        }
        c.n = 0
    case .Device_Gray, .Device_RGB, .Cal_Gray, .Cal_RGB, .Lab, .ICC_Based, .Indexed:
        c.n = n
    }
    return
}

// colorspace_default_decode is the /Decode array an image in this space gets
// when it does not supply one (ISO 32000-2, Table 89). `bpc` matters only for
// an Indexed space, whose default maps the raw sample straight onto the
// lookup table.
//
// The result is allocated with `allocator` and belongs to the caller.
colorspace_default_decode :: proc(cs: ^Colorspace, bpc: int, allocator := context.allocator) -> []f64 {
    if cs == nil {
        return nil
    }
    n := min(cs.n, MAX_COLOR_COMPONENTS)
    if n < 1 {
        return nil
    }
    out := make([]f64, 2 * n, allocator)

    if cs.kind == .Indexed {
        bits := clamp(bpc, 1, 16)
        out[0] = 0
        out[1] = f64(int(1) << uint(bits)) - 1
        return out
    }
    for j in 0 ..< n {
        out[2 * j], out[2 * j + 1] = cs_component_range(cs, j)
    }
    return out
}