package pdf

import "core:math"

// Shadings (ISO 32000-2, 8.7.4.5). Types 1-3 are evaluated per pixel from a
// function; types 4-7 are meshes that are decomposed into triangles.
//
// HOW A SHADING REACHES THE PIXELS
//
// Everything here paints through `blend_span_rgba8`, one horizontal run at a
// time, in straight RGBA8 -- the same road an image takes, so a shading, an
// image and a fill all round identically and a page cannot show a seam where
// one meets another. Runs are assembled in `Shade_Ctx.row`, one scratch row as
// wide as the clip, and a pixel the shading does not cover is left with alpha
// zero and skipped rather than painted transparent, because "not painted" and
// "painted with alpha 0" differ under every blend mode but Normal.
//
// Types 1-3 walk the clip box row by row, mapping each pixel centre back
// through the inverse CTM into the shading's own space. Types 4-7 read
// triangles out of the mesh stream and scan-convert each one with a top-left
// fill rule at pixel centres. That rule is not decoration: adjacent triangles
// share an edge exactly, and a pixel centre landing on a shared edge must be
// claimed by exactly one of them or the mesh grows either a bright seam (the
// pixel painted twice, visible whenever alpha < 1) or a dark one (painted by
// neither). Triangles are drawn without anti-aliasing for the same reason --
// two anti-aliased edges meeting would composite twice and leave a lighter
// line down every shared edge. Triangles that genuinely overlap -- which only
// a self-intersecting mesh produces -- are composited in the order the stream
// lists them, so the later one shows, which is what the format asks for; with
// a constant alpha below 1 the overlap is painted twice, and correcting that
// would mean rendering the whole shading into a buffer of its own for a case
// no real file contains.
//
// WHAT IS NOT DONE HERE
//
// /Background is parsed but never painted. It applies only when a shading is
// used as a pattern's fill, and even there only outside the shading's own
// extent; no viewer this library was checked against paints it, and painting
// it would make `sh` disagree with `scn` for the same shading dictionary.
// /AntiAlias is likewise parsed and ignored -- it is a hint, and the mesh
// rasteriser's fill rule is chosen for seam-freedom rather than smoothness.

Shading_Kind :: enum u8 {
    Function_Based    = 1,
    Axial             = 2,
    Radial            = 3,
    Free_Triangle     = 4,
    Lattice_Triangle  = 5,
    Coons             = 6,
    Tensor            = 7,
}

Shading :: struct {
    kind: Shading_Kind,
    cs:   ^Colorspace,
    // One function producing cs.n outputs, or cs.n functions each producing
    // one. Both spellings are legal and both occur.
    funcs: []^Function,

    // /Background, painted outside the shading's own extent when the shading
    // is used as a pattern (it is ignored by the `sh` operator).
    background:     []f64,
    has_background: bool,
    // /BBox in the shading's own space.
    bbox:     Rect,
    has_bbox: bool,
    antialias: bool,

    // Type 1.
    domain: [4]f64,
    mtx:    Matrix,

    // Types 2 and 3: /Coords is 4 numbers for axial, 6 for radial.
    coords: [6]f64,
    // /Domain, defaulting to [0 1].
    t0, t1: f64,
    extend: [2]bool,

    // Types 4-7: the mesh stream and its sample layout.
    stream:            ^Stream,
    bits_per_coord:    int,
    bits_per_component: int,
    bits_per_flag:     int,
    decode:            []f64,
    // Type 5 only.
    vertices_per_row:  int,

    // Colour values per mesh vertex: 1 when /Function maps a parametric value,
    // otherwise one per component of `cs`. It fixes the size of every vertex
    // in the stream, so it is resolved once at load time.
    n_comp: int,
}

// ---------------------------------------------------------------------------
// Limits
//
// Every one of these bounds a number that came out of a file.
// ---------------------------------------------------------------------------

// SHADE_LUT_SIZE is how finely the colour function of an axial or radial
// shading is sampled. The function is 1-D there, so sampling it once per
// entry and looking the colour up per pixel turns a type 4 PostScript
// calculator running several million times into 256 runs. 256 entries put the
// error below the 8-bit quantisation the output has anyway.
@(private)
SHADE_LUT_SIZE :: 256

// Vertices one mesh stream may yield. A megabyte of 16-bit vertices is already
// far past any real gradient mesh.
@(private)
MAX_MESH_VERTICES :: 1 << 20

// /VerticesPerRow. Two is the minimum that forms a triangle; the cap is what
// stops a four-byte integer from sizing an allocation.
@(private)
MAX_MESH_ROW :: 1 << 14

@(private)
MAX_MESH_TRIANGLES :: 1 << 20

@(private)
MAX_MESH_PATCHES :: 1 << 16

// Steps a patch is subdivided into along each parameter. 32 gives 2048
// triangles for the largest patch, which is smooth well past 300 dpi.
@(private)
MAX_PATCH_SUBDIV :: 32

// Device pixels a patch aims to put between two subdivision steps.
@(private)
PATCH_STEP_PIXELS :: 3.0

// Pixels one shading may examine. This is not a quality knob -- a full page at
// 600 dpi is under 35 million -- it is the stop on a mesh built to make the
// rasteriser walk the clip a thousand times over.
@(private)
SHADE_PIXEL_FUEL :: 1 << 26

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

// shading_load builds the shading `obj` describes: a dictionary for types 1-3,
// a stream for types 4-7 (whose dictionary carries the same entries).
//
// `res` is the /Resources dictionary the shading was named from, needed only
// to resolve a colour space given as a bare name.
//
// The result is owned by the Document and must not be freed. Nothing is
// cached: unlike a colour space or a function, a shading is built once per
// `sh` operator or once per pattern, and the mesh stream it points at is
// decoded lazily by `stream_data`, which does its own caching.
shading_load :: proc(doc: ^Document, obj: Object, res: ^Dict) -> (out: ^Shading, err: Error) {
    if doc == nil {
        return nil, .Invalid_Argument
    }
    o := resolve(doc, obj)
    d, has_dict := to_dict(o)
    if !has_dict {
        return nil, .Bad_Shading
    }

    st := dict_int(doc, d, "ShadingType", 0)
    if st < 1 || st > 7 {
        return nil, .Bad_Shading
    }

    sh := new(Shading, doc.allocator)
    sh.kind = Shading_Kind(u8(st))

    cs, cs_err := colorspace_load(doc, dict_get(doc, d, "ColorSpace"), res)
    // A shading whose colour space will not load cannot be guessed at: for a
    // mesh the component count is what fixes the size of every vertex in the
    // stream, so a wrong guess does not degrade, it garbles.
    if cs_err != .None || cs == nil || cs.n < 1 || cs.n > MAX_COLOR_COMPONENTS {
        return nil, .Bad_Shading
    }
    // 8.7.4.5.1: the colour space of a shading shall not be a Pattern space.
    if cs.kind == .Pattern {
        return nil, .Bad_Shading
    }
    sh.cs = cs

    sh.funcs = shade_load_functions(doc, dict_get(doc, d, "Function"))
    // Types 1-3 have no colour of their own; without a function there is
    // nothing to paint with.
    if len(sh.funcs) == 0 && (sh.kind == .Function_Based || sh.kind == .Axial || sh.kind == .Radial) {
        return nil, .Bad_Shading
    }

    if ba := dict_array(doc, d, "Background"); array_len(ba) > 0 {
        n := min(array_len(ba), cs.n)
        bg := make([]f64, n, doc.allocator)
        for i in 0 ..< n {
            bg[i] = clamp_finite(array_num(doc, ba, i), 0)
        }
        sh.background = bg
        sh.has_background = true
    }

    if bb, bb_ok := object_to_rect(doc, dict_get(doc, d, "BBox")); bb_ok {
        sh.bbox = bb
        sh.has_bbox = true
    }
    sh.antialias = dict_bool(doc, d, "AntiAlias")

    switch sh.kind {
    case .Function_Based:
        shade_load_type1(doc, sh, d)
    case .Axial, .Radial:
        shade_load_axial_radial(doc, sh, d) or_return
    case .Free_Triangle, .Lattice_Triangle, .Coons, .Tensor:
        stm, is_stream := to_stream(o)
        if !is_stream {
            return nil, .Bad_Shading
        }
        sh.stream = stm
        shade_load_mesh(doc, sh, d) or_return
    }
    return sh, .None
}

// shade_load_functions reads /Function, which is either one function with as
// many outputs as the colour space has components or an array of that many
// one-output functions. Both spellings are legal and both occur, so the array
// is what is stored and the count is what tells them apart.
//
// The array form is INDEXED BY COMPONENT, so its length and the position of
// every entry in it are load-bearing: a sub-function that will not load holds
// its slot as nil rather than being dropped. An array in which nothing loads
// comes back empty, which for types 1-3 is what makes the shading fail.
@(private)
shade_load_functions :: proc(doc: ^Document, obj: Object) -> []^Function {
    if obj == nil {
        return nil
    }
    if a, is_array := to_array(obj); is_array {
        count := min(array_len(a), MAX_COLOR_COMPONENTS)
        if count == 0 {
            return nil
        }
        out, aerr := make([]^Function, count, doc.allocator)
        if aerr != nil {
            return nil
        }
        any := false
        for i in 0 ..< count {
            f, err := function_load(doc, array_get(doc, a, i))
            if err != .None || f == nil {
                // One broken sub-function must not lose the others, and it
                // must not move them either: the slot is left nil so that
                // every surviving function stays on the component the file
                // wrote it for, and only the component this one would have
                // supplied reads as zero. Dropping the entry instead would
                // shift every later function down one index and repaint the
                // whole shading in the wrong colours.
                add_repair(doc, .Bad_Resource, 0)
                continue
            }
            out[i] = f
            any = true
        }
        if !any {
            return nil
        }
        return out
    }
    f, err := function_load(doc, obj)
    if err != .None || f == nil {
        return nil
    }
    one := make([]^Function, 1, doc.allocator)
    one[0] = f
    return one
}

@(private)
shade_load_type1 :: proc(doc: ^Document, sh: ^Shading, d: ^Dict) {
    sh.domain = {0, 1, 0, 1}
    if da := dict_array(doc, d, "Domain"); array_len(da) >= 4 {
        dom: [4]f64
        ok := true
        for i in 0 ..< 4 {
            dom[i] = array_num(doc, da, i)
            if !is_finite(dom[i]) {
                ok = false
            }
        }
        // A reversed or non-finite /Domain is not worth failing the shading
        // over; normalising it paints the same region the file meant.
        if ok {
            sh.domain = {
                min(dom[0], dom[1]), max(dom[0], dom[1]),
                min(dom[2], dom[3]), max(dom[2], dom[3]),
            }
        }
    }
    sh.mtx = dict_matrix(doc, d, "Matrix")
}

@(private)
shade_load_axial_radial :: proc(doc: ^Document, sh: ^Shading, d: ^Dict) -> Error {
    need := sh.kind == .Axial ? 4 : 6
    ca := dict_array(doc, d, "Coords")
    if array_len(ca) < need {
        return .Bad_Shading
    }
    for i in 0 ..< need {
        sh.coords[i] = array_num(doc, ca, i)
        if !is_finite(sh.coords[i]) {
            return .Bad_Shading
        }
    }

    sh.t0, sh.t1 = 0, 1
    if da := dict_array(doc, d, "Domain"); array_len(da) >= 2 {
        a := array_num(doc, da, 0)
        b := array_num(doc, da, 1)
        if is_finite(a) && is_finite(b) {
            sh.t0, sh.t1 = a, b
        }
    }

    if ea := dict_array(doc, d, "Extend"); array_len(ea) >= 2 {
        sh.extend[0], _ = to_bool(array_get(doc, ea, 0))
        sh.extend[1], _ = to_bool(array_get(doc, ea, 1))
    }
    return .None
}

@(private)
shade_load_mesh :: proc(doc: ^Document, sh: ^Shading, d: ^Dict) -> Error {
    bpc := int(dict_int(doc, d, "BitsPerCoordinate", 0))
    switch bpc {
    case 1, 2, 4, 8, 12, 16, 24, 32:
        sh.bits_per_coord = bpc
    case:
        return .Bad_Shading
    }

    bpcomp := int(dict_int(doc, d, "BitsPerComponent", 0))
    switch bpcomp {
    case 1, 2, 4, 8, 12, 16:
        sh.bits_per_component = bpcomp
    case:
        return .Bad_Shading
    }

    if sh.kind != .Lattice_Triangle {
        bpf := int(dict_int(doc, d, "BitsPerFlag", 0))
        switch bpf {
        case 2, 4, 8:
            sh.bits_per_flag = bpf
        case:
            return .Bad_Shading
        }
    }

    // With a /Function the stream carries one parametric value per vertex, not
    // one value per colour component (8.7.4.5.5).
    sh.n_comp = len(sh.funcs) > 0 ? 1 : sh.cs.n
    need := 4 + 2 * sh.n_comp
    da := dict_array(doc, d, "Decode")
    if array_len(da) < need {
        return .Bad_Shading
    }
    dec := make([]f64, need, doc.allocator)
    for i in 0 ..< need {
        dec[i] = array_num(doc, da, i)
        if !is_finite(dec[i]) {
            return .Bad_Shading
        }
    }
    sh.decode = dec

    if sh.kind == .Lattice_Triangle {
        vpr := int(dict_int(doc, d, "VerticesPerRow", 0))
        if vpr < 2 || vpr > MAX_MESH_ROW {
            return .Bad_Shading
        }
        sh.vertices_per_row = vpr
    }
    return .None
}

// ---------------------------------------------------------------------------
// Colour
// ---------------------------------------------------------------------------

// shade_components_to_rgb converts raw colour components -- the ones a mesh
// vertex carries when there is no /Function -- through the shading's space.
@(private)
shade_components_to_rgb :: proc(sh: ^Shading, comps: []f64) -> [3]f64 {
    c := Color{n = sh.cs.n}
    for i in 0 ..< min(sh.cs.n, len(comps)) {
        c.v[i] = comps[i]
    }
    return colorspace_to_rgb(sh.cs, c)
}

// shade_eval runs the shading's function(s) at `in_` and converts the result
// to RGB. `in_` is one parametric value for types 2-7 and two coordinates for
// type 1.
@(private)
shade_eval :: proc(sh: ^Shading, in_: []f64) -> [3]f64 {
    comps: [MAX_COLOR_COMPONENTS]f64
    n := sh.cs.n
    switch len(sh.funcs) {
    case 0:
        // Only reachable for a mesh whose colours come from the stream; the
        // caller uses shade_components_to_rgb there.
        return [3]f64{0, 0, 0}
    case 1:
        // One function supplying every component: the array spelling with a
        // single entry means the same thing as the bare one, and both occur.
        function_eval(sh.funcs[0], in_, comps[:n])
    case:
        one: [1]f64
        for i in 0 ..< min(len(sh.funcs), n) {
            // A nil slot is a sub-function that would not load. Its component
            // stays zero; the ones around it keep their own indices.
            if sh.funcs[i] == nil {
                continue
            }
            function_eval(sh.funcs[i], in_, one[:])
            comps[i] = one[0]
        }
    }
    return shade_components_to_rgb(sh, comps[:n])
}

// Shade_Lut is the colour of an axial or radial shading sampled across its
// /Domain. Entry i is the colour at t0 + (i/(N-1)) * (t1 - t0).
@(private)
Shade_Lut :: struct {
    rgb: [SHADE_LUT_SIZE][3]f64,
}

@(private)
shade_build_lut :: proc(sh: ^Shading, lut: ^Shade_Lut) {
    in_: [1]f64
    for i in 0 ..< SHADE_LUT_SIZE {
        s := f64(i) / f64(SHADE_LUT_SIZE - 1)
        in_[0] = sh.t0 + s * (sh.t1 - sh.t0)
        lut.rgb[i] = shade_eval(sh, in_[:])
    }
}

// shade_lut_at looks a parametric position 0..1 up, rounding to the nearest
// sample so the error is half a step rather than a whole one.
@(private)
shade_lut_at :: proc(lut: ^Shade_Lut, s: f64) -> [3]f64 {
    if !(s > 0) {
        return lut.rgb[0]
    }
    if s >= 1 {
        return lut.rgb[SHADE_LUT_SIZE - 1]
    }
    i := int(s * f64(SHADE_LUT_SIZE - 1) + 0.5)
    return lut.rgb[clamp(i, 0, SHADE_LUT_SIZE - 1)]
}

// ---------------------------------------------------------------------------
// Bounds
// ---------------------------------------------------------------------------

// shading_bounds is the device-space region the shading can paint under `ctm`.
//
// RECT_INFINITE means "unbounded", which is the honest answer for an axial
// shading (the strip between its two ends is infinite across the axis, and an
// /Extend makes it infinite along it too) and for an extended radial one. It
// is the identity for `rect_intersect`, so a caller that intersects it with a
// clip gets the clip, which is exactly right.
shading_bounds :: proc(sh: ^Shading, ctm: Matrix) -> Rect {
    if sh == nil || !matrix_is_finite(ctm) {
        return RECT_EMPTY
    }

    r := RECT_INFINITE
    switch sh.kind {
    case .Function_Based:
        dom := Rect{sh.domain[0], sh.domain[2], sh.domain[1], sh.domain[3]}
        r = transform_rect(matrix_mul(sh.mtx, ctm), dom)
    case .Axial:
        // Unbounded across the axis whatever /Extend says.
        r = RECT_INFINITE
    case .Radial:
        if !sh.extend[0] && !sh.extend[1] {
            // Every circle between the two ends is a convex combination of a
            // point of the first disc and a point of the second, so the two
            // discs bound the whole family.
            //
            // The corners are taken with min and max rather than by unioning
            // two rectangles: a circle of radius zero -- which is how every
            // "glow" gradient starts -- is an empty rectangle, and
            // `rect_union` drops an empty operand, which would shrink the
            // bound to the other circle and clip the narrow end of the cone
            // away.
            rr0 := max(sh.coords[2], 0)
            rr1 := max(sh.coords[5], 0)
            r = transform_rect(ctm, Rect{
                min(sh.coords[0] - rr0, sh.coords[3] - rr1),
                min(sh.coords[1] - rr0, sh.coords[4] - rr1),
                max(sh.coords[0] + rr0, sh.coords[3] + rr1),
                max(sh.coords[1] + rr0, sh.coords[4] + rr1),
            })
        }
    case .Free_Triangle, .Lattice_Triangle, .Coons, .Tensor:
        // /Decode fixes the range every coordinate in the stream decodes into,
        // so the mesh is bounded without reading a byte of it.
        if len(sh.decode) >= 4 {
            r = transform_rect(ctm, rect_normalize(
                Rect{sh.decode[0], sh.decode[2], sh.decode[1], sh.decode[3]}))
        }
    }

    if sh.has_bbox {
        r = rect_intersect(r, transform_rect(ctm, sh.bbox))
    }
    if !rect_is_finite(r) {
        return RECT_INFINITE
    }
    return r
}

// ---------------------------------------------------------------------------
// Painting
// ---------------------------------------------------------------------------

@(private)
Shade_Ctx :: struct {
    doc:   ^Document,
    sh:    ^Shading,
    bmp:   ^Bitmap,
    clip:  Raster_Clip,
    alpha: f64,
    mode:  Blend_Mode,
    mask:  ^Bitmap,

    ctm: Matrix,
    // Device space back to the shading's own space.
    inv: Matrix,

    // One scratch row of straight RGBA8, indexed from clip.x0. Alpha 0 marks a
    // pixel the shading does not cover.
    row: []u8,
    // The mask's coverage for `mask_row_y`, one byte per pixel from clip.x0,
    // or nil when there is no mask.
    mask_row:   []u8,
    mask_row_y: int,

    fuel: int,
    tris: int,
}

// shading_paint composites `sh` onto `bmp` through `clip`.
//
//   ctm    maps the shading's own space to device space
//   clip   device pixels, half-open, intersected with the bitmap here
//   alpha  the constant alpha of the paint operation
//   mask   an optional soft mask aligned with `bmp` in device pixels, one
//          component read (what `mask_from_bitmap` produces). Pixels outside
//          it read as 255, so a mask smaller than the clip shows content
//          rather than erasing it.
//
// `bmp` is PREMULTIPLIED, like everything below the Device line. A singular or
// non-finite CTM paints nothing, and so does a colour space of /None.
// `context.allocator` supplies scratch only; everything it takes is released
// before this returns.
shading_paint :: proc(doc: ^Document, sh: ^Shading, bmp: ^Bitmap, ctm: Matrix,
                      clip: Raster_Clip, alpha: f64, mode: Blend_Mode, mask: ^Bitmap) -> Error {
    if sh == nil || bmp == nil || sh.cs == nil {
        return .Invalid_Argument
    }
    // A /None separation marks nothing at all on the page (8.6.6.4).
    if sh.cs.is_none {
        return .None
    }
    if !matrix_is_finite(ctm) {
        return .None
    }
    inv, invertible := matrix_invert(ctm)
    if !invertible {
        return .None
    }

    ga := blend_clamp01(alpha)
    if ga <= 0 {
        return .None
    }

    cl := Raster_Clip{
        x0 = max(clip.x0, 0),
        y0 = max(clip.y0, 0),
        x1 = min(clip.x1, bmp.width),
        y1 = min(clip.y1, bmp.height),
    }
    if cl.x1 <= cl.x0 || cl.y1 <= cl.y0 {
        return .None
    }

    // Narrowing to what the shading can possibly reach keeps a small gradient
    // on a big page from walking the whole page, and is what applies /BBox to
    // a mesh.
    b := shading_bounds(sh, ctm)
    cl = Raster_Clip{
        x0 = raster_floor_clamped(b.x0, cl.x0, cl.x1),
        y0 = raster_floor_clamped(b.y0, cl.y0, cl.y1),
        x1 = raster_ceil_clamped(b.x1, cl.x0, cl.x1),
        y1 = raster_ceil_clamped(b.y1, cl.y0, cl.y1),
    }
    if cl.x1 <= cl.x0 || cl.y1 <= cl.y0 {
        return .None
    }

    width := cl.x1 - cl.x0
    row, row_err := make([]u8, width * 4, context.allocator)
    if row_err != nil {
        return .Out_Of_Memory
    }
    defer delete(row, context.allocator)

    ctx := Shade_Ctx{
        doc        = doc,
        sh         = sh,
        bmp        = bmp,
        clip       = cl,
        alpha      = ga,
        mode       = mode,
        mask       = mask,
        ctm        = ctm,
        inv        = inv,
        row        = row,
        mask_row_y = -1,
        fuel       = SHADE_PIXEL_FUEL,
    }

    if mask != nil {
        mrow, mask_err := make([]u8, width, context.allocator)
        if mask_err != nil {
            return .Out_Of_Memory
        }
        ctx.mask_row = mrow
    }
    // Freed here rather than in the branch above: a `defer` runs at the end of
    // the block it sits in, so one inside the `if` would free the row before a
    // single span had been composited through it.
    defer if ctx.mask_row != nil { delete(ctx.mask_row, context.allocator) }

    switch sh.kind {
    case .Function_Based:
        shade_paint_function(&ctx)
    case .Axial:
        shade_paint_axial(&ctx)
    case .Radial:
        shade_paint_radial(&ctx)
    case .Free_Triangle, .Lattice_Triangle, .Coons, .Tensor:
        return shade_paint_mesh(&ctx)
    }
    return .None
}

// shade_put writes one covered pixel into the scratch row.
@(private)
shade_put :: proc(ctx: ^Shade_Ctx, x: int, rgb: [3]f64) {
    i := (x - ctx.clip.x0) * 4
    if i < 0 || i + 4 > len(ctx.row) {
        return
    }
    ctx.row[i]     = blend_quantize(rgb[0])
    ctx.row[i + 1] = blend_quantize(rgb[1])
    ctx.row[i + 2] = blend_quantize(rgb[2])
    ctx.row[i + 3] = 255
}

@(private)
shade_skip :: proc(ctx: ^Shade_Ctx, x: int) {
    i := (x - ctx.clip.x0) * 4
    if i < 0 || i + 4 > len(ctx.row) {
        return
    }
    ctx.row[i + 3] = 0
}

// shade_mask_row lifts the mask's coverage for scanline `y` into a byte per
// pixel, which is the form `blend_span_rgba8` indexes.
@(private)
shade_mask_row :: proc(ctx: ^Shade_Ctx, y: int) {
    if ctx.mask_row == nil || ctx.mask_row_y == y {
        return
    }
    ctx.mask_row_y = y
    for i in 0 ..< len(ctx.mask_row) {
        ctx.mask_row[i] = 255
    }
    if ctx.mask == nil {
        return
    }
    src := bitmap_row(ctx.mask, y)
    if src == nil {
        return
    }
    n := ctx.mask.n
    for i in 0 ..< len(ctx.mask_row) {
        x := ctx.clip.x0 + i
        if x >= ctx.mask.width {
            break
        }
        off := x * n
        if off < len(src) {
            ctx.mask_row[i] = src[off]
        }
    }
}

// shade_flush composites the covered runs of `row` between x_from and x_to.
@(private)
shade_flush :: proc(ctx: ^Shade_Ctx, y, x_from, x_to: int) {
    lo := max(x_from, ctx.clip.x0)
    hi := min(x_to, ctx.clip.x1)
    if lo >= hi {
        return
    }
    shade_mask_row(ctx, y)

    i := lo
    for i < hi {
        if ctx.row[(i - ctx.clip.x0) * 4 + 3] == 0 {
            i += 1
            continue
        }
        j := i
        for j < hi && ctx.row[(j - ctx.clip.x0) * 4 + 3] != 0 {
            j += 1
        }
        a := (i - ctx.clip.x0) * 4
        b := (j - ctx.clip.x0) * 4
        mrow: []u8
        if ctx.mask_row != nil {
            mrow = ctx.mask_row[i - ctx.clip.x0 : j - ctx.clip.x0]
        }
        blend_span_rgba8(ctx.bmp, y, i, ctx.row[a:b], ctx.alpha, ctx.mode, nil, mrow)
        i = j
    }
}

// shade_in_bbox tests a device pixel centre against /BBox, which is given in
// the shading's own space. Testing there rather than against the transformed
// box keeps a rotated CTM from painting the corners the box excludes.
@(private)
shade_in_bbox :: proc(ctx: ^Shade_Ctx, p: Point) -> bool {
    if !ctx.sh.has_bbox {
        return true
    }
    return rect_contains(ctx.sh.bbox, p)
}

// ---------------------------------------------------------------------------
// Type 1: function-based
// ---------------------------------------------------------------------------

@(private)
shade_paint_function :: proc(ctx: ^Shade_Ctx) {
    sh := ctx.sh
    // /Matrix maps the domain into the shading's space, which the CTM then
    // maps to the device; one inverse takes a pixel all the way back.
    full, ok := matrix_invert(matrix_mul(sh.mtx, ctx.ctm))
    if !ok {
        return
    }

    in_: [2]f64
    for y in ctx.clip.y0 ..< ctx.clip.y1 {
        for x in ctx.clip.x0 ..< ctx.clip.x1 {
            dev := Point{f64(x) + 0.5, f64(y) + 0.5}
            p := transform_point(full, dev)
            if !is_finite(p.x) || !is_finite(p.y) ||
               p.x < sh.domain[0] || p.x > sh.domain[1] ||
               p.y < sh.domain[2] || p.y > sh.domain[3] {
                shade_skip(ctx, x)
                continue
            }
            if sh.has_bbox && !shade_in_bbox(ctx, transform_point(ctx.inv, dev)) {
                shade_skip(ctx, x)
                continue
            }
            in_[0] = p.x
            in_[1] = p.y
            shade_put(ctx, x, shade_eval(sh, in_[:]))
        }
        shade_flush(ctx, y, ctx.clip.x0, ctx.clip.x1)
        ctx.fuel -= ctx.clip.x1 - ctx.clip.x0
        if ctx.fuel <= 0 {
            return
        }
    }
}

// ---------------------------------------------------------------------------
// Type 2: axial
// ---------------------------------------------------------------------------

@(private)
shade_paint_axial :: proc(ctx: ^Shade_Ctx) {
    sh := ctx.sh
    p0 := Point{sh.coords[0], sh.coords[1]}
    d := Point{sh.coords[2] - p0.x, sh.coords[3] - p0.y}
    len2 := d.x * d.x + d.y * d.y
    // A zero-length axis leaves the parametric variable undefined; there is no
    // colour to choose, so nothing is painted.
    if !(len2 > 0) {
        return
    }
    inv_len2 := 1 / len2

    lut: Shade_Lut
    shade_build_lut(sh, &lut)

    for y in ctx.clip.y0 ..< ctx.clip.y1 {
        for x in ctx.clip.x0 ..< ctx.clip.x1 {
            p := transform_point(ctx.inv, Point{f64(x) + 0.5, f64(y) + 0.5})
            s := ((p.x - p0.x) * d.x + (p.y - p0.y) * d.y) * inv_len2
            if !is_finite(s) || !shade_in_bbox(ctx, p) {
                shade_skip(ctx, x)
                continue
            }
            if s < 0 {
                if !sh.extend[0] {
                    shade_skip(ctx, x)
                    continue
                }
                s = 0
            } else if s > 1 {
                if !sh.extend[1] {
                    shade_skip(ctx, x)
                    continue
                }
                s = 1
            }
            shade_put(ctx, x, shade_lut_at(&lut, s))
        }
        shade_flush(ctx, y, ctx.clip.x0, ctx.clip.x1)
        ctx.fuel -= ctx.clip.x1 - ctx.clip.x0
        if ctx.fuel <= 0 {
            return
        }
    }
}

// ---------------------------------------------------------------------------
// Type 3: radial
// ---------------------------------------------------------------------------

// shade_radial_accept tests one candidate s against the two rules that decide
// whether a circle in the family exists there at all: its radius must not be
// negative (8.7.4.5.4 -- the extension stops where the circles would invert),
// and a value outside 0..1 needs the matching /Extend. The returned value is
// clamped, which is what makes an extended shading paint its end colour.
@(private)
shade_radial_accept :: proc(s, r0, dr: f64, extend: [2]bool) -> (f64, bool) {
    if !is_finite(s) {
        return 0, false
    }
    if r0 + s * dr < 0 {
        return 0, false
    }
    if s < 0 {
        return 0, extend[0]
    }
    if s > 1 {
        return 1, extend[1]
    }
    return s, true
}

// shade_radial_solve finds the largest s whose circle passes through the point,
// as the format requires: later circles paint over earlier ones, so the largest
// s is the one that shows.
//
//   a = dx^2 + dy^2 - dr^2
//   b = px*dx + py*dy + r0*dr
//   c = px^2 + py^2 - r0^2
//
// solving a*s^2 - 2*b*s + c = 0. `a` is zero exactly when the circles' radii
// grow as fast as their centres move -- the cone degenerates to a paraboloid
// and the equation is linear. A negative discriminant means no circle in the
// family reaches the point, which must paint nothing rather than an end colour.
@(private)
shade_radial_solve :: proc(a, b, c, r0, dr: f64, extend: [2]bool) -> (f64, bool) {
    if a == 0 {
        if b == 0 {
            return 0, false
        }
        return shade_radial_accept(c / (2 * b), r0, dr, extend)
    }
    disc := b * b - a * c
    if !(disc >= 0) {
        return 0, false
    }
    sq := math.sqrt(disc)
    s1 := (b + sq) / a
    s2 := (b - sq) / a
    hi := max(s1, s2)
    lo := min(s1, s2)
    if s, ok := shade_radial_accept(hi, r0, dr, extend); ok {
        return s, true
    }
    return shade_radial_accept(lo, r0, dr, extend)
}

@(private)
shade_paint_radial :: proc(ctx: ^Shade_Ctx) {
    sh := ctx.sh
    x0, y0, r0 := sh.coords[0], sh.coords[1], sh.coords[2]
    dx := sh.coords[3] - x0
    dy := sh.coords[4] - y0
    dr := sh.coords[5] - r0
    a := dx * dx + dy * dy - dr * dr

    lut: Shade_Lut
    shade_build_lut(sh, &lut)

    for y in ctx.clip.y0 ..< ctx.clip.y1 {
        for x in ctx.clip.x0 ..< ctx.clip.x1 {
            p := transform_point(ctx.inv, Point{f64(x) + 0.5, f64(y) + 0.5})
            if !is_finite(p.x) || !is_finite(p.y) || !shade_in_bbox(ctx, p) {
                shade_skip(ctx, x)
                continue
            }
            px := p.x - x0
            py := p.y - y0
            b := px * dx + py * dy + r0 * dr
            c := px * px + py * py - r0 * r0
            s, ok := shade_radial_solve(a, b, c, r0, dr, sh.extend)
            if !ok {
                shade_skip(ctx, x)
                continue
            }
            shade_put(ctx, x, shade_lut_at(&lut, s))
        }
        shade_flush(ctx, y, ctx.clip.x0, ctx.clip.x1)
        ctx.fuel -= ctx.clip.x1 - ctx.clip.x0
        if ctx.fuel <= 0 {
            return
        }
    }
}

// ---------------------------------------------------------------------------
// Triangles
// ---------------------------------------------------------------------------

@(private)
Shade_Vertex :: struct {
    // Device space: the mesh is transformed as it is read, so the whole
    // rasteriser works in one space.
    p:   Point,
    rgb: [3]f64,
}

// shade_top_left reports whether an edge owns the pixel centres that land
// exactly on it. With triangles wound to a positive area, the interior lies on
// the e >= 0 side, which makes an edge running upwards a left edge and a
// horizontal edge running right a top edge. Two triangles sharing an edge
// traverse it in opposite directions, so exactly one of them owns it.
@(private)
shade_top_left :: proc "contextless" (d: Point) -> bool {
    return d.y < 0 || (d.y == 0 && d.x > 0)
}

@(private)
shade_edge_accept :: proc "contextless" (w: f64, top_left: bool) -> bool {
    // Written so a NaN is outside rather than inside.
    return w > 0 || (w == 0 && top_left)
}

// shade_triangle scan-converts one Gouraud-shaded triangle at pixel centres.
@(private)
shade_triangle :: proc(ctx: ^Shade_Ctx, va, vb, vc: Shade_Vertex) {
    a := va
    b := vb
    c := vc
    if !is_finite(a.p.x) || !is_finite(a.p.y) ||
       !is_finite(b.p.x) || !is_finite(b.p.y) ||
       !is_finite(c.p.x) || !is_finite(c.p.y) {
        return
    }

    area := (b.p.x - a.p.x) * (c.p.y - a.p.y) - (b.p.y - a.p.y) * (c.p.x - a.p.x)
    if !is_finite(area) || area == 0 {
        return
    }
    if area < 0 {
        // Winding every triangle the same way is what makes the top-left rule
        // consistent between neighbours: an adjacent triangle on the far side
        // of a shared edge necessarily traverses that edge the other way.
        b, c = c, b
        area = -area
    }
    ctx.tris += 1

    bx0 := raster_floor_clamped(min(a.p.x, b.p.x, c.p.x), ctx.clip.x0, ctx.clip.x1)
    bx1 := raster_ceil_clamped(max(a.p.x, b.p.x, c.p.x), ctx.clip.x0, ctx.clip.x1)
    by0 := raster_floor_clamped(min(a.p.y, b.p.y, c.p.y), ctx.clip.y0, ctx.clip.y1)
    by1 := raster_ceil_clamped(max(a.p.y, b.p.y, c.p.y), ctx.clip.y0, ctx.clip.y1)
    if bx1 <= bx0 || by1 <= by0 {
        return
    }

    d0 := Point{b.p.x - a.p.x, b.p.y - a.p.y}
    d1 := Point{c.p.x - b.p.x, c.p.y - b.p.y}
    d2 := Point{a.p.x - c.p.x, a.p.y - c.p.y}
    tl0 := shade_top_left(d0)
    tl1 := shade_top_left(d1)
    tl2 := shade_top_left(d2)
    inv_area := 1 / area

    for y in by0 ..< by1 {
        py := f64(y) + 0.5
        for x in bx0 ..< bx1 {
            px := f64(x) + 0.5
            // Recomputed rather than stepped along the row on purpose: the
            // same edge evaluated from two neighbouring triangles must give
            // bit-identical results, and an incremental sum would depend on
            // where each triangle's bounding box started.
            w0 := d0.x * (py - a.p.y) - d0.y * (px - a.p.x)
            w1 := d1.x * (py - b.p.y) - d1.y * (px - b.p.x)
            w2 := d2.x * (py - c.p.y) - d2.y * (px - c.p.x)
            if !shade_edge_accept(w0, tl0) || !shade_edge_accept(w1, tl1) ||
               !shade_edge_accept(w2, tl2) {
                shade_skip(ctx, x)
                continue
            }
            if ctx.sh.has_bbox && !shade_in_bbox(ctx, transform_point(ctx.inv, Point{px, py})) {
                shade_skip(ctx, x)
                continue
            }
            // w1 is the edge opposite a, w2 opposite b, w0 opposite c.
            la := w1 * inv_area
            lb := w2 * inv_area
            lc := w0 * inv_area
            rgb := [3]f64{
                la * a.rgb[0] + lb * b.rgb[0] + lc * c.rgb[0],
                la * a.rgb[1] + lb * b.rgb[1] + lc * c.rgb[1],
                la * a.rgb[2] + lb * b.rgb[2] + lc * c.rgb[2],
            }
            shade_put(ctx, x, rgb)
        }
        shade_flush(ctx, y, bx0, bx1)
        ctx.fuel -= bx1 - bx0
        if ctx.fuel <= 0 {
            return
        }
    }
}

// ---------------------------------------------------------------------------
// Mesh streams
// ---------------------------------------------------------------------------

@(private)
Mesh_Reader :: struct {
    sh:   ^Shading,
    ctm:  Matrix,
    data: []byte,
    // Bit position, MSB first, which is how every packed field in a PDF is
    // written.
    pos:  int,
    max_coord: f64,
    max_comp:  f64,
    count: int,
}

@(private)
mesh_reader_make :: proc(sh: ^Shading, ctm: Matrix, data: []byte) -> Mesh_Reader {
    return Mesh_Reader{
        sh        = sh,
        ctm       = ctm,
        data      = data,
        max_coord = math.pow(2, f64(sh.bits_per_coord)) - 1,
        max_comp  = math.pow(2, f64(sh.bits_per_component)) - 1,
    }
}

// mesh_bits reads the next `n` bits. Returns ok = false at end of data, and
// leaves the reader at the end so every later read fails too -- a truncated
// mesh stops, it does not wrap around and paint garbage.
@(private)
mesh_bits :: proc(r: ^Mesh_Reader, n: int) -> (v: u32, ok: bool) {
    total := len(r.data) * 8
    if n <= 0 || n > 32 || r.pos + n > total {
        r.pos = total
        return 0, false
    }
    acc: u64
    left := n
    for left > 0 {
        byte_index := r.pos >> 3
        bit_off := r.pos & 7
        avail := 8 - bit_off
        take := min(avail, left)
        chunk := (u64(r.data[byte_index]) >> uint(avail - take)) & u64((1 << uint(take)) - 1)
        acc = acc << uint(take) | chunk
        r.pos += take
        left -= take
    }
    return u32(acc), true
}

// mesh_align skips to the next byte boundary. Free-form triangles and patches
// pad each record to a whole number of bytes; lattice rows do not.
@(private)
mesh_align :: proc(r: ^Mesh_Reader) {
    r.pos = ((r.pos + 7) / 8) * 8
}

@(private)
mesh_decode :: proc(v: u32, maxv, lo, hi: f64) -> f64 {
    if !(maxv > 0) {
        return lo
    }
    return lo + f64(v) * (hi - lo) / maxv
}

@(private)
mesh_point :: proc(r: ^Mesh_Reader) -> (Point, bool) {
    xv, ok0 := mesh_bits(r, r.sh.bits_per_coord)
    yv, ok1 := mesh_bits(r, r.sh.bits_per_coord)
    if !ok0 || !ok1 {
        return Point{}, false
    }
    dec := r.sh.decode
    p := Point{
        mesh_decode(xv, r.max_coord, dec[0], dec[1]),
        mesh_decode(yv, r.max_coord, dec[2], dec[3]),
    }
    return transform_point(r.ctm, p), true
}

@(private)
mesh_color :: proc(r: ^Mesh_Reader) -> ([3]f64, bool) {
    comps: [MAX_COLOR_COMPONENTS]f64
    sh := r.sh
    dec := sh.decode
    for i in 0 ..< sh.n_comp {
        cv, ok := mesh_bits(r, sh.bits_per_component)
        if !ok {
            return [3]f64{}, false
        }
        comps[i] = mesh_decode(cv, r.max_comp, dec[4 + 2 * i], dec[5 + 2 * i])
    }
    if len(sh.funcs) > 0 {
        return shade_eval(sh, comps[:1]), true
    }
    return shade_components_to_rgb(sh, comps[:sh.n_comp]), true
}

@(private)
mesh_vertex :: proc(r: ^Mesh_Reader, out: ^Shade_Vertex) -> bool {
    if r.count >= MAX_MESH_VERTICES {
        return false
    }
    p, p_ok := mesh_point(r)
    if !p_ok {
        return false
    }
    rgb, c_ok := mesh_color(r)
    if !c_ok {
        return false
    }
    out.p = p
    out.rgb = rgb
    r.count += 1
    return true
}

@(private)
shade_paint_mesh :: proc(ctx: ^Shade_Ctx) -> Error {
    sh := ctx.sh
    if sh.stream == nil || len(sh.decode) < 4 + 2 * sh.n_comp {
        return .None
    }
    data, err := stream_data(ctx.doc, sh.stream)
    if err != .None || len(data) == 0 {
        return .None
    }

    r := mesh_reader_make(sh, ctx.ctm, data)
    switch sh.kind {
    case .Free_Triangle:
        shade_mesh_free(ctx, &r)
    case .Lattice_Triangle:
        return shade_mesh_lattice(ctx, &r)
    case .Coons, .Tensor:
        shade_mesh_patches(ctx, &r)
    case .Function_Based, .Axial, .Radial:
        // Not reachable: shading_paint dispatches those elsewhere.
    }
    return .None
}

// ---------------------------------------------------------------------------
// Type 4: free-form triangles
// ---------------------------------------------------------------------------

@(private)
mesh_flagged_vertex :: proc(r: ^Mesh_Reader, out: ^Shade_Vertex) -> bool {
    // The two vertices that follow a flag-0 vertex carry a flag of their own,
    // which the format says shall be 0 and which nothing reads.
    if _, ok := mesh_bits(r, r.sh.bits_per_flag); !ok {
        return false
    }
    if !mesh_vertex(r, out) {
        return false
    }
    mesh_align(r)
    return true
}

@(private)
shade_mesh_free :: proc(ctx: ^Shade_Ctx, r: ^Mesh_Reader) {
    va, vb, vc: Shade_Vertex
    have := false

    for ctx.tris < MAX_MESH_TRIANGLES && ctx.fuel > 0 {
        flag, flag_ok := mesh_bits(r, r.sh.bits_per_flag)
        if !flag_ok {
            return
        }
        v: Shade_Vertex
        if !mesh_vertex(r, &v) {
            return
        }
        mesh_align(r)

        switch {
        case flag == 0 || !have:
            // A continuation flag with no previous triangle is malformed;
            // starting a fresh triangle is the recovery that loses least.
            v1, v2: Shade_Vertex
            if !mesh_flagged_vertex(r, &v1) || !mesh_flagged_vertex(r, &v2) {
                return
            }
            va, vb, vc = v, v1, v2
        case flag == 1:
            va, vb, vc = vb, vc, v
        case flag == 2:
            va, vb, vc = va, vc, v
        case:
            add_repair(ctx.doc, .Bad_Resource, 0)
            return
        }
        shade_triangle(ctx, va, vb, vc)
        have = true
    }
}

// ---------------------------------------------------------------------------
// Type 5: lattice-form triangles
// ---------------------------------------------------------------------------

@(private)
shade_mesh_lattice :: proc(ctx: ^Shade_Ctx, r: ^Mesh_Reader) -> Error {
    per_row := ctx.sh.vertices_per_row
    if per_row < 2 || per_row > MAX_MESH_ROW {
        return .None
    }

    rows, rows_err := make([]Shade_Vertex, per_row * 2, context.allocator)
    if rows_err != nil {
        return .Out_Of_Memory
    }
    defer delete(rows, context.allocator)
    prev := rows[:per_row]
    cur := rows[per_row:]

    have_prev := false
    for ctx.tris < MAX_MESH_TRIANGLES && ctx.fuel > 0 {
        complete := true
        for i in 0 ..< per_row {
            if !mesh_vertex(r, &cur[i]) {
                complete = false
                break
            }
        }
        // A row that ran out mid-way is dropped whole: half a row has no
        // triangles in it.
        if !complete {
            return .None
        }
        if have_prev {
            for i in 0 ..< per_row - 1 {
                shade_triangle(ctx, prev[i], prev[i + 1], cur[i])
                shade_triangle(ctx, prev[i + 1], cur[i + 1], cur[i])
            }
        }
        prev, cur = cur, prev
        have_prev = true
    }
    return .None
}

// ---------------------------------------------------------------------------
// Types 6 and 7: Coons and tensor patches
// ---------------------------------------------------------------------------

// The twelve boundary control points are written anticlockwise from one
// corner, which lays them onto the 4x4 tensor grid like this (ISO 32000-2,
// figures 46 and 47):
//
//        p1  p2  p3  p4          [0][0] [0][1] [0][2] [0][3]
//       p12          p5          [1][0]  p13    p14   [1][3]
//       p11          p6          [2][0]  p16    p15   [2][3]
//       p10  p9  p8  p7          [3][0] [3][1] [3][2] [3][3]
//
// with the corner colours c1..c4 at p1, p4, p7 and p10 in that order.
@(private)
PATCH_GRID := [12][2]int{
    {0, 0}, {0, 1}, {0, 2}, {0, 3},
    {1, 3}, {2, 3}, {3, 3},
    {3, 2}, {3, 1}, {3, 0},
    {2, 0}, {1, 0},
}

// PATCH_SHARED tells which of the previous patch's points and colours a flag
// of 1, 2 or 3 reuses as the new patch's first edge (Table 85).
@(private)
PATCH_SHARED := [4][4]int{
    {0, 0, 0, 0},      // flag 0: nothing is shared
    {3, 4, 5, 6},
    {6, 7, 8, 9},
    {9, 10, 11, 0},
}

@(private)
PATCH_SHARED_COLOR := [4][2]int{
    {0, 0},
    {1, 2},
    {2, 3},
    {3, 0},
}

@(private)
shade_mesh_patches :: proc(ctx: ^Shade_Ctx, r: ^Mesh_Reader) {
    n_points := ctx.sh.kind == .Tensor ? 16 : 12

    pts, prev_pts: [16]Point
    cols, prev_cols: [4][3]f64
    have := false

    for _ in 0 ..< MAX_MESH_PATCHES {
        if ctx.tris >= MAX_MESH_TRIANGLES || ctx.fuel <= 0 {
            return
        }

        raw_flag, flag_ok := mesh_bits(r, r.sh.bits_per_flag)
        if !flag_ok {
            return
        }
        flag := int(raw_flag)
        if flag > 3 {
            add_repair(ctx.doc, .Bad_Resource, 0)
            return
        }

        first_point := 0
        first_color := 0
        if flag != 0 && have {
            share := PATCH_SHARED[flag]
            for k in 0 ..< 4 {
                pts[k] = prev_pts[share[k]]
            }
            cshare := PATCH_SHARED_COLOR[flag]
            cols[0] = prev_cols[cshare[0]]
            cols[1] = prev_cols[cshare[1]]
            first_point = 4
            first_color = 2
        }

        for k in first_point ..< n_points {
            p, ok := mesh_point(r)
            if !ok {
                return
            }
            pts[k] = p
        }
        for k in first_color ..< 4 {
            c, ok := mesh_color(r)
            if !ok {
                return
            }
            cols[k] = c
        }
        mesh_align(r)

        shade_patch(ctx, pts, cols, ctx.sh.kind == .Tensor)
        prev_pts = pts
        prev_cols = cols
        have = true
    }
}

@(private)
bernstein3 :: proc "contextless" (t: f64) -> [4]f64 {
    mt := 1 - t
    return [4]f64{mt * mt * mt, 3 * mt * mt * t, 3 * mt * t * t, t * t * t}
}

@(private)
patch_eval :: proc(pole: ^[4][4]Point, u, v: f64) -> (p: Point) {
    bu := bernstein3(u)
    bv := bernstein3(v)
    for i in 0 ..< 4 {
        for j in 0 ..< 4 {
            w := bu[i] * bv[j]
            p.x += w * pole[i][j].x
            p.y += w * pole[i][j].y
        }
    }
    return
}

// patch_color is the bilinear interpolation of the four corner colours the
// format prescribes (8.7.4.5.7): c1 at (0,0), c2 at (0,1), c3 at (1,1) and
// c4 at (1,0), matching the corners the control points put there.
@(private)
patch_color :: proc(c: [4][3]f64, u, v: f64) -> (out: [3]f64) {
    for i in 0 ..< 3 {
        e0 := (1 - v) * c[0][i] + v * c[1][i]
        e1 := (1 - v) * c[3][i] + v * c[2][i]
        out[i] = (1 - u) * e0 + u * e1
    }
    return
}

// shade_patch subdivides one patch into triangles.
//
// The subdivision is uniform, with its step chosen from the patch's device
// size so that a step is about PATCH_STEP_PIXELS across: a patch a few pixels
// wide costs two triangles, and one covering a page is split as far as the cap
// allows. Uniform subdivision rather than adaptive flatness testing is a
// deliberate trade -- a Bezier patch's colour varies over its whole area, not
// just where it curves, so the flat parts need the same subdivision as the
// curved ones and adaptivity would buy nothing.
@(private)
shade_patch :: proc(ctx: ^Shade_Ctx, v: [16]Point, c: [4][3]f64, tensor: bool) {
    pole: [4][4]Point
    for k in 0 ..< 12 {
        g := PATCH_GRID[k]
        pole[g[0]][g[1]] = v[k]
    }
    if tensor {
        pole[1][1] = v[12]
        pole[1][2] = v[13]
        pole[2][2] = v[14]
        pole[2][1] = v[15]
    } else {
        // A Coons patch gives no interior control points; the format defines
        // them from the boundary (8.7.4.5.7).
        pole[1][1] = shade_coons_inner(pole[0][0], pole[0][1], pole[1][0], pole[0][3], pole[3][0], pole[3][1], pole[1][3], pole[3][3])
        pole[1][2] = shade_coons_inner(pole[0][3], pole[0][2], pole[1][3], pole[0][0], pole[3][3], pole[3][2], pole[1][0], pole[3][0])
        pole[2][1] = shade_coons_inner(pole[3][0], pole[3][1], pole[2][0], pole[3][3], pole[0][0], pole[0][1], pole[2][3], pole[0][3])
        pole[2][2] = shade_coons_inner(pole[3][3], pole[3][2], pole[2][3], pole[3][0], pole[0][3], pole[0][2], pole[2][0], pole[0][0])
    }

    lo := pole[0][0]
    hi := pole[0][0]
    for i in 0 ..< 4 {
        for j in 0 ..< 4 {
            p := pole[i][j]
            if !is_finite(p.x) || !is_finite(p.y) {
                return
            }
            lo.x = min(lo.x, p.x); lo.y = min(lo.y, p.y)
            hi.x = max(hi.x, p.x); hi.y = max(hi.y, p.y)
        }
    }
    extent := max(hi.x - lo.x, hi.y - lo.y)
    n := 1
    if is_finite(extent) && extent > 0 {
        steps := math.ceil(extent / PATCH_STEP_PIXELS)
        if steps >= f64(MAX_PATCH_SUBDIV) {
            n = MAX_PATCH_SUBDIV
        } else if steps > 1 {
            n = int(steps)
        }
    }
    inv_n := 1 / f64(n)

    row_a, row_b: [MAX_PATCH_SUBDIV + 1]Point
    for j in 0 ..= n {
        row_a[j] = patch_eval(&pole, 0, f64(j) * inv_n)
    }
    for i in 1 ..= n {
        u1 := f64(i) * inv_n
        u0 := f64(i - 1) * inv_n
        for j in 0 ..= n {
            row_b[j] = patch_eval(&pole, u1, f64(j) * inv_n)
        }
        for j in 0 ..< n {
            v0 := f64(j) * inv_n
            v1 := f64(j + 1) * inv_n
            p00 := Shade_Vertex{row_a[j],     patch_color(c, u0, v0)}
            p01 := Shade_Vertex{row_a[j + 1], patch_color(c, u0, v1)}
            p11 := Shade_Vertex{row_b[j + 1], patch_color(c, u1, v1)}
            p10 := Shade_Vertex{row_b[j],     patch_color(c, u1, v0)}
            shade_triangle(ctx, p00, p01, p11)
            shade_triangle(ctx, p00, p11, p10)
        }
        row_a = row_b
        if ctx.fuel <= 0 || ctx.tris >= MAX_MESH_TRIANGLES {
            return
        }
    }
}

// shade_coons_inner is the interior control point the format derives from the
// boundary of a Coons patch:
//
//   p11 = (-4*p00 + 6*(p01 + p10) - 2*(p03 + p30) + 3*(p31 + p13) - p33) / 9
//
// with the arguments named for the corner it belongs to; the other three come
// from the same expression with the grid rotated.
@(private)
shade_coons_inner :: proc "contextless" (corner, along_a, along_b, far_a, far_b, near_a, near_b, opposite: Point) -> Point {
    return Point{
        (-4 * corner.x + 6 * (along_a.x + along_b.x) - 2 * (far_a.x + far_b.x) +
            3 * (near_a.x + near_b.x) - opposite.x) / 9,
        (-4 * corner.y + 6 * (along_a.y + along_b.y) - 2 * (far_a.y + far_b.y) +
            3 * (near_a.y + near_b.y) - opposite.y) / 9,
    }
}