begin
    import Pkg
    Pkg.activate(joinpath(@__DIR__, ".."))

    # a writable font cache for CairoMakie (a notebook worker may have no HOME)
    ENV["XDG_CACHE_HOME"] = get(ENV, "XDG_CACHE_HOME", joinpath(homedir(), ".cache"))
end;

import Vasculature as V

import CairoMakie as Mke

import GPUSelect

using PlutoUI: Slider, Select, TableOfContents

using Printf: @sprintf

using Serialization: serialize, deserialize

using Mmap: mmap

using Statistics: median

begin
    # licensed: these files and anything made from them stay on lab storage
    const XCAT_DIR = get(
        ENV,
        "XCAT_DIR",
        joinpath(
            homedir(),
            "Documents",
            "dev",
            "MolloiLab",
            "xcat-data",
            "vmale_50_heart_1600_1400_867",
        ),
    )
    # the full chest, arms included, at the same 0.2 mm: XCAT's heart block inside it and its surface meshes voxelized
    # around it (built once by xcat-data/work/phantom/build_chest.jl), so a scanner sees the whole body
    const XCAT_CHEST_DIR = get(
        ENV,
        "XCAT_CHEST_DIR",
        joinpath(
            homedir(),
            "Documents",
            "dev",
            "MolloiLab",
            "xcat-data",
            "work",
            "phantom",
        ),
    )
    const CACHE_DIR = joinpath(homedir(), "xcat-perfused", "cardiac_pipeline_cache")
    const PHANTOM_DIR = joinpath(homedir(), "xcat-perfused", "cardiac_pipeline")

    mkpath(CACHE_DIR)
    mkpath(PHANTOM_DIR)
end;

"Arrays to this machine's GPU (CUDA, Metal, ROCm, oneAPI), or the CPU if there is none."
const DEVICE = GPUSelect.Storage()

"Run `compute` once and keep its result in `CACHE_DIR/name.jls`; afterwards load it."
function cached(compute, name)
    path = joinpath(CACHE_DIR, name * ".jls")
    isfile(path) && return deserialize(path)

    seconds = @elapsed result = compute()
    serialize(path, result)
    @info "$name: $(round(seconds; digits = 1)) s"
    return result
end

begin
    # the chest grid: origin (mm, the centre of voxel (1, 1, 1)), voxels, voxel size (mm)
    const CHEST_GRID = readlines(joinpath(XCAT_CHEST_DIR, "chest_grid.txt"))
    const XCAT_ORIGIN = Tuple(parse.(Float64, split(CHEST_GRID[1])))
    const XCAT_DIMS = Tuple(parse.(Int, split(CHEST_GRID[2])))
    const XCAT_VOXEL_MM = parse(Float64, CHEST_GRID[3])

    xcat_labels = mmap(
        open(joinpath(XCAT_CHEST_DIR, "chest_labels_0p2.u8")),
        Array{UInt8, 3},
        XCAT_DIMS,
    )
    xcat_surfaces = V.XCAT.read_nrb(joinpath(XCAT_DIR, "vmale_50_heart.nrb"))
end;

begin
    "XCAT label groups for pictures: (name, colour, XCAT labels)."
    const LABEL_GROUPS = [
        ("myocardium", :darkorange, (15, 16)),
        ("chamber blood", :mediumpurple, (19, 20, 21, 22)),
        ("great arteries", :firebrick, (25, 28)),
        ("great veins", :royalblue, (24,)),
        ("coronary arteries", :red, (26,)),
        ("coronary veins", :deepskyblue, (27,)),
        ("lung", :gray25, (6,)),
        ("fat", :khaki, (1,)),
        ("muscle", :rosybrown, (4,)),
        ("bone", :ivory, (5, 8, 9)),
        ("other", :gray60, ()),
    ]

    "The picture group of every XCAT label (0 = air)."
    const GROUP_OF_LABEL = let group = fill(UInt8(length(LABEL_GROUPS)), 256)
        group[1] = 0
        for (k, (_, _, labels)) in enumerate(LABEL_GROUPS), label in labels
            group[label + 1] = k
        end
        group
    end

    label_group_colormap() = Mke.cgrad(
        [:black; [colour for (_, colour, _) in LABEL_GROUPS]],
        length(LABEL_GROUPS) + 1;
        categorical = true,
    )
    xcat_slice(z) = clamp(
        round(Int, (z - XCAT_ORIGIN[3]) / XCAT_VOXEL_MM) + 1,
        1,
        XCAT_DIMS[3],
    )
    xcat_x_mm = XCAT_ORIGIN[1] .+ (0:(XCAT_DIMS[1] - 1)) .* XCAT_VOXEL_MM
    xcat_y_mm = XCAT_ORIGIN[2] .+ (0:(XCAT_DIMS[2] - 1)) .* XCAT_VOXEL_MM
end;

begin
    # shades for vessels by radius, pale for the thinnest shown and deep for the largest:
    # arteries red, veins blue
    ARTERY_COLOURS = Mke.cgrad(["#f4a582", "#d6604d", "#b2182b", "#67001f"])
    VEIN_COLOURS = Mke.cgrad(["#92c5de", "#4393c3", "#2166ac", "#053061"])
end;

"""
A 3-D axis for the heart in true proportions, seen from the front: head up, the patient's left on the right. A gentle
perspective, no box or grid.
"""
function vessel_axis(position; title)
    ax = Mke.Axis3(
        position;
        title = title,
        aspect = :data,
        azimuth = 1.5π,
        elevation = 0.08π,
        perspectiveness = 0.3,
        viewmode = :fit,
        protrusions = 0,
    )
    Mke.hidedecorations!(ax)
    Mke.hidespines!(ax)
    return ax
end

growth_grid = let (points, _) = V.XCAT.tessellate(
        xcat_surfaces["dias_pericardium"];
        nu = 60,
        nv = 60,
    )
    lower = floor.(Int, vec(minimum(points; dims = 2))) .- 2
    upper = ceil.(Int, vec(maximum(points; dims = 2))) .+ 2
    V.VoxelGrid(Tuple(Float64.(lower)), 1.0, Tuple(upper .- lower .+ 1))
end

heart_masks = cached("masks") do
    pericardium = V.voxelize(
        V.XCAT.tessellate(
            xcat_surfaces["dias_pericardium"];
            nu = 160,
            nv = 160,
        )...,
        growth_grid,
    )
    V.XCAT.label_masks(
        xcat_labels,
        XCAT_ORIGIN,
        XCAT_VOXEL_MM,
        growth_grid;
        roi = pericardium,
    )
end;

begin
    myocardium = heart_masks.lvmyo .| heart_masks.rvmyo
    forbid = (heart_masks.pool .& .!myocardium) .| heart_masks.lvcav .| heart_masks.rvcav
end;

"""
The myocardium's surface: each surface voxel's centre (mm) and its number among the myocardial voxels, so per-voxel
values (such as MBF) can be painted on it.
"""
MYOCARDIUM_SURFACE = let
    neighbours = (
        CartesianIndex(1, 0, 0),
        CartesianIndex(-1, 0, 0),
        CartesianIndex(0, 1, 0),
        CartesianIndex(0, -1, 0),
        CartesianIndex(0, 0, 1),
        CartesianIndex(0, 0, -1),
    )
    on_surface(I) = any(
        !checkbounds(Bool, myocardium, I + d) || !myocardium[I + d] for d in neighbours
    )
    voxels = findall(myocardium)
    surface = [n for (n, I) in enumerate(voxels) if on_surface(I)]
    (
        points = [Mke.Point3f(Tuple(voxels[n]) .- 1 .+ growth_grid.origin) for n in surface],
        voxels = surface,
    )
end;

"Draw every other point of the myocardium's surface as a faint shell around the vessels."
function draw_myocardium_shell!(ax)
    return Mke.scatter!(
        ax,
        MYOCARDIUM_SURFACE.points[1:2:end];
        color = (:gray50, 0.07),
        markersize = 5,
        strokewidth = 0,
    )
end

"""
Draw a grown tree's segments of radius ≥ `min_radius` (mm): width and shade grow with radius, and the largest vessels
are drawn last so they sit on top.
"""
function draw_vessels!(ax, tree; colormap, min_radius)
    shown = filter(j -> tree.r[j] >= min_radius, sortperm(tree.r))
    segments = [
        Mke.Point3f(tree.a[:, j]) => Mke.Point3f(tree.b[:, j])
            for j in shown
    ]
    radius = tree.r[shown]
    # CairoMakie wants a colour and a width at each end of every segment
    return Mke.linesegments!(
        ax,
        segments;
        color = repeat(log10.(radius); inner = 2),
        colormap = colormap,
        colorrange = (log10(min_radius), log10(maximum(tree.r))),
        linewidth = repeat(1.5 .* sqrt.(radius ./ min_radius); inner = 2),
    )
end

"Blood (chambers and great vessels) at 0.4 mm over the growth grid's box, as `(mask, origin_mm, cell_mm)`: the blockers."
BLOOD_BLOCKERS = cached("blood") do
    first_voxel = [
        round(Int, (growth_grid.origin[d] - 0.5 - XCAT_ORIGIN[d]) / XCAT_VOXEL_MM) + 1
            for d in 1:3
    ]
    cells = [ceil(Int, (growth_grid.dims[d] + 1) / 0.4) for d in 1:3]
    is_blood = falses(256)
    for label in (19, 20, 21, 22, 24, 25, 28)
        is_blood[label + 1] = true
    end

    blood = zeros(UInt8, cells...)
    Threads.@threads for c in 1:cells[3]
        for b in 1:cells[2], a in 1:cells[1], dk in 0:1, dj in 0:1, di in 0:1
            i = first_voxel[1] + 2(a - 1) + di
            j = first_voxel[2] + 2(b - 1) + dj
            k = first_voxel[3] + 2(c - 1) + dk
            (
                1 <= i <= XCAT_DIMS[1] &&
                    1 <= j <= XCAT_DIMS[2] &&
                    1 <= k <= XCAT_DIMS[3]
            ) || continue
            is_blood[xcat_labels[i, j, k] + 1] && (blood[a, b, c] = 0x01)
        end
    end
    origin = ntuple(
        d -> XCAT_ORIGIN[d] + XCAT_VOXEL_MM * (first_voxel[d] - 1) + XCAT_VOXEL_MM / 2,
        3,
    )
    (blood, origin, 0.4)
end;

artery_inputs = cached("artery_inputs") do
    names = ["dias_lad1", "dias_lad2", "dias_lad3", "dias_lcx", "dias_rca1", "dias_rca2"]
    trees = V.XCAT.base_forest(
        [V.XCAT.centerline(xcat_surfaces[name]; step = 2.8) for name in names];
        roots = ["dias_lad1", "dias_rca1"],
    )
    stubbed = [
        V.XCAT.penetrating_stubs(
            tree,
            growth_grid,
            [heart_masks.lvmyo, heart_masks.rvmyo];
            domain = myocardium,
            forbid,
            blockers = BLOOD_BLOCKERS,
        )
            for tree in trees
    ]
    lumens = [
        V.XCAT.vessel_tube(tree, growth_grid) .|
            V.XCAT.vessel_bridges(tree, myocardium, growth_grid)
            for tree in trees
    ]
    (
        trees = first.(stubbed),
        stubs = last.(stubbed),
        pass = reduce(.|, lumens) .& .!forbid .& .!myocardium,
    )
end;

vein_inputs = cached("vein_inputs") do
    vein1 = V.XCAT.centerline(xcat_surfaces["dias_cor_vein1"]; step = 2.8)
    vein2 = V.XCAT.centerline(xcat_surfaces["dias_cor_vein2"]; step = 2.8)

    # the coronary sinus ostium: the point of cor_vein1 nearest the right atrium, joined to it
    # by a 6 mm connector
    atrium, _ = V.XCAT.tessellate(
        xcat_surfaces["dias_ra"];
        nu = 120,
        nv = 120,
    )
    nearest_atrium(p) = begin
        d = vec(sqrt.(sum(abs2, atrium .- p; dims = 1)))
        k = argmin(d)
        (d[k], atrium[:, k])
    end
    ostium = argmin([nearest_atrium(vein1.X[:, k])[1] for k in 1:size(vein1.X, 2)])
    gap, atrium_point = nearest_atrium(vein1.X[:, ostium])
    start =
        vein1.X[:, ostium] .-
        6.0 .* (vein1.X[:, ostium] .- atrium_point) ./ max(gap, 1.0e-6)

    lines = [
        V.XCAT.Centerline(
            "cs_ostium",
            hcat(start, vein1.X[:, ostium]),
            fill(vein1.r[ostium], 2),
        ),
        V.XCAT.Centerline("cor_vein1a", vein1.X[:, ostium:-1:1], vein1.r[ostium:-1:1]),
        V.XCAT.Centerline("cor_vein1b", vein1.X[:, ostium:end], vein1.r[ostium:end]),
        vein2,
    ]
    base = only(
        V.XCAT.base_forest(
            lines;
            roots = ["cs_ostium"],
            attach_tol = 4.0,
        ),
    )
    # tributaries join the imaged veins along their length, as branches join the imaged arteries (imaged vessels are fixed
    # in growth, so without stubs new vessels could join only at their tips)
    tree, stubs = V.XCAT.penetrating_stubs(
        base,
        growth_grid,
        [heart_masks.lvmyo, heart_masks.rvmyo];
        domain = myocardium,
        forbid,
        blockers = BLOOD_BLOCKERS,
    )
    lumen =
        V.XCAT.vessel_tube(tree, growth_grid) .|
        V.XCAT.vessel_bridges(tree, myocardium, growth_grid)
    (tree = tree, stubs = stubs, pass = lumen .& .!forbid .& .!myocardium)
end;

"Grow one forest: the domain and its pass-through voxels as a `.vol` file, the inlet trees as GXL files."
function grow_forest(spec, trees, pass, name)
    dir = mktempdir()
    lower = Tuple(round.(Int, growth_grid.origin))
    allow = zeros(UInt16, growth_grid.dims)
    allow[pass] .= V.NOTERM # vessels may cross these voxels, no terminal may end in them

    domain = V.write_vol(joinpath(dir, "$name.vol"), myocardium .| pass, lower)
    inlets = [
        V.XCAT.write_base_gxl(joinpath(dir, "$(name)_$k.xml"), tree)
            for (k, tree) in enumerate(trees)
    ]
    engine = V.grow(
        V.GrowthSpec(spec; device = DEVICE),
        domain,
        inlets;
        blockers = BLOOD_BLOCKERS,
        territories = allow,
    )
    return V.vtree(engine)
end

arteries = cached("arteries") do
    grow_forest(V.HEART_GROWTH, artery_inputs.trees, artery_inputs.pass, "arteries")
end;

veins = cached("veins") do
    grow_forest(V.HEART_VEIN_GROWTH, [vein_inputs.tree], vein_inputs.pass, "veins")
end;

"Centres (mm) of the 1 mm myocardial voxels: the tissue the flow is solved for."
tissue_centres = let voxels = findall(myocardium)
    Float32[growth_grid.origin[d] + voxels[k][d] - 1 for d in 1:3, k in eachindex(voxels)]
end

flows = cached("flow") do
    V.flow(
        V.HEART_FLOW,
        arteries,
        veins,
        tissue_centres;
        to = DEVICE,
    )
end;

"Blood flow (mL/min/g) of every 1 mm myocardial voxel in circuit `circuit`."
function voxel_mbf(circuit)
    terminal_flow = reshape(Float32[circuit.q[j] for j in circuit.terminals], 1, :)
    voxels = collect(Int32(1):Int32(size(tissue_centres, 2)))
    # mL/s per mL of tissue
    per_ml = V.OrganCT.tissue_iodine(terminal_flow, flows.pairs, flows.terr, voxels)
    return vec(per_ml) .* 60 ./ V.HEART_FLOW.rho
end

"A 1 mm map (x, y) of `values` (one per myocardial voxel) at z (mm); NaN elsewhere."
function myocardium_slice(values, z)
    k = clamp(round(Int, z - growth_grid.origin[3]) + 1, 1, growth_grid.dims[3])
    slice = fill(NaN32, growth_grid.dims[1], growth_grid.dims[2])
    for (n, voxel) in enumerate(findall(myocardium))
        voxel[3] == k && (slice[voxel[1], voxel[2]] = values[n])
    end
    return slice
end

"""
The injection: the lab's current protocol (Hadjiabdolhamid et al. 2024, Shojazadeh et al. 2026), 0.5 mL/kg of contrast at
5 mL/s, then 0.5 mL/kg of a 30:70 contrast:saline chaser at the same rate, for the body model's 70 kg: 35 mL of iohexol
350, then 35 mL at 105 mg/mL. A short bolus keeps the first pass inside the scan window: at stress, iodine starts leaving
the heart about 2.5 s after it arrives.
"""
const INJECTION = let kg = 70, conc = 350
    V.Body.Injection(; phases = [(0.5 * kg, 5, conc), (0.5 * kg, 5, 0.3 * conc)])
end

"Contrast transport through the body and the coronary circuit `circuit`, at rest or adenosine stress."
function transport_state(state, circuit)
    heart = state == "stress" ? V.Body.adenosine() : (co = 81.6f0, beat = 1.0f0)
    spec = V.TransportSpec(
        V.HEART_TRANSPORT;
        phys = circuit.phys,
        tissue_ml = Float32(size(tissue_centres, 2) / 1000),
        heart...,
    )
    return (
        state,
        circuit,
        V.transport(
            spec,
            circuit,
            INJECTION;
            times = V.HEART_PHANTOM.times,
        ),
    )
end

"""
Contrast transport at rest and stress through `circuits`, solved once and kept in `CACHE_DIR/name.jls`; returns that
file. The full results (every vessel's iodine at every frame) take several GB, so the notebook does not hold them:
`save_phantom` reads them while it writes, and the figures use `transport_curves`.
"""
function transport_cache(name, circuits)
    path = joinpath(CACHE_DIR, name * ".jls")
    isfile(path) || serialize(
        path,
        [
            transport_state("rest", circuits.rest),
            transport_state("stress", circuits.stress),
        ],
    )
    return path
end

TRANSPORT = transport_cache("transport", flows);

"Each state's LV blood and whole-myocardium iodine (mg/mL) over time, from a transport file."
function transport_curves(path)
    tissue_ml = size(tissue_centres, 2) / 1000
    return map(deserialize(path)) do (state, _, result)
        held = vec(
            sum(result.tissue; dims = 2) .+
                sum(result.artery; dims = 2) .+
                sum(result.vein; dims = 2),
        )
        (
            state = state,
            blood_times = result.body.t,
            lv = result.body.mean[:lv],
            myocardium = held ./ tissue_ml,
        )
    end
end

TRANSPORT_CURVES = transport_curves(TRANSPORT);

begin
    split_walls = true # LV and RV apart
    split_territories = true # LAD, LCx and RCA apart
    epicardial_compartments = true # the surface coronary arteries and veins as compartments of their own
end;

begin
    wall_names = split_walls ? ["LV", "RV"] : ["Myocardium"]
    territory_names = split_territories ? ["LAD", "LCx", "RCA"] : String[]
    wall_rois =
        isempty(territory_names) ? wall_names :
        ["$w · $t" for w in wall_names for t in territory_names]
    epicardial_names = epicardial_compartments ? ["$w epicardial vessels" for w in wall_names] : String[]
    ROI_NAMES = [wall_rois; "Papillary muscles"; epicardial_names]

    roi_index(name) = findfirst(==(name), ROI_NAMES)
    ROI_GROUPS = Pair{String, Vector{Int}}[]
    # a wall's own compartments: its territories, or the wall itself
    wall_ids(w) = isempty(territory_names) ? [roi_index(w)] : [roi_index("$w · $t") for t in territory_names]
    heart_ids = collect(eachindex(wall_rois))
    length(heart_ids) > 1 && push!(ROI_GROUPS, "heart myocardium" => heart_ids)
    if epicardial_compartments
        push!(ROI_GROUPS, "heart myocardium + epicardial vessels" => [heart_ids; roi_index.(epicardial_names)])
    end
    if length(wall_names) > 1
        for w in wall_names
            length(wall_ids(w)) > 1 && push!(ROI_GROUPS, "$w myocardium" => wall_ids(w))
            epicardial_compartments && push!(ROI_GROUPS, "$w myocardium + epicardial vessels" => [wall_ids(w); roi_index("$w epicardial vessels")])
        end
        if !isempty(territory_names)
            append!(ROI_GROUPS, [t => [roi_index("$w · $t") for w in wall_names] for t in territory_names])
        end
    end
    ROI_NAMES, ROI_GROUPS
end

"""
Each arterial terminal's branch (`of`: 1 LAD, 2 LCx, 3 RCA) and the LAD's first segment (`lad`). The LCA splits at its
first imaged bifurcation where both branches feed ≥ 5% of it; the more anterior one is the LAD. Terminals of the left main
before it count as LAD.
"""
coronary_branch = let tree = arteries, n = length(tree.r)
    children = [Int[] for _ in 1:n]
    for j in 1:n
        tree.parent[j] > 0 && push!(children[tree.parent[j]], j)
    end
    root_of(j) = begin
        k = j
        while tree.parent[k] != 0
            k = tree.parent[k]
        end
        k
    end
    depth(j) = begin
        k = j
        d = 0
        while tree.parent[k] != 0
            k = tree.parent[k]
            d += 1
        end
        d
    end

    # terminals below each segment, and their summed y (front to back)
    terminals_below = Float64.(length.(children) .== 0)
    y_below = terminals_below .* tree.b[2, :]
    for j in sortperm(depth.(1:n); rev = true)
        p = tree.parent[j]
        p > 0 || continue
        terminals_below[p] += terminals_below[j]
        y_below[p] += y_below[j]
    end

    roots = unique(root_of.(1:n))
    lca = roots[argmax(terminals_below[roots])]
    j = lca
    lad = lcx = 0
    while lad == 0
        big = [
            c
                for c in children[j]
                if tree.isbase[c] && terminals_below[c] > 0.05 * terminals_below[lca]
        ]
        if length(big) >= 2
            front_to_back = sortperm([y_below[c] / terminals_below[c] for c in big])
            lad, lcx = big[front_to_back[1]], big[front_to_back[2]]
        else
            imaged = [c for c in children[j] if tree.isbase[c]]
            j = imaged[argmax(terminals_below[imaged])]
        end
    end

    below(j, a) = begin
        k = j
        while k != 0
            k == a && return true
            k = tree.parent[k]
        end
        false
    end
    branch_of(j) = root_of(j) != lca ? 3 : below(j, lcx) ? 2 : 1
    (of = branch_of.(flows.rest.terminals), lad = lad)
end;

"Each 1 mm voxel's territory (1–3), carried 25 mm out from the myocardium so the lumens on its surface get one too."
territory_map = let weights = zeros(Float32, 3, size(tissue_centres, 2)),
        map = zeros(Int8, growth_grid.dims)
    for (voxel, terminal, w) in flows.pairs
        weights[coronary_branch.of[terminal], voxel] += w
    end
    voxels = findall(myocardium)
    for (n, voxel) in enumerate(voxels)
        map[voxel] = argmax(@view weights[:, n])
    end

    steps = (
        CartesianIndex(1, 0, 0),
        CartesianIndex(-1, 0, 0),
        CartesianIndex(0, 1, 0),
        CartesianIndex(0, -1, 0),
        CartesianIndex(0, 0, 1),
        CartesianIndex(0, 0, -1),
    )
    front = voxels
    for _ in 1:25
        next = CartesianIndex{3}[]
        for voxel in front, step in steps
            neighbour = voxel + step
            checkbounds(Bool, map, neighbour) && map[neighbour] == 0 || continue
            map[neighbour] = map[voxel]
            push!(next, neighbour)
        end
        front = next
    end
    map
end;

"Each 1 mm voxel's wall (1 LV, 2 RV), carried 25 mm out from the myocardium: the wall a surface vessel lies beside."
wall_map = let map = zeros(Int8, growth_grid.dims)
    map[heart_masks.lvmyo] .= 1
    map[heart_masks.rvmyo .& .!heart_masks.lvmyo] .= 2
    steps = (
        CartesianIndex(1, 0, 0),
        CartesianIndex(-1, 0, 0),
        CartesianIndex(0, 1, 0),
        CartesianIndex(0, -1, 0),
        CartesianIndex(0, 0, 1),
        CartesianIndex(0, 0, -1),
    )
    front = findall(!=(0), map)
    for _ in 1:25
        next = CartesianIndex{3}[]
        for voxel in front, step in steps
            neighbour = voxel + step
            checkbounds(Bool, map, neighbour) && map[neighbour] == 0 || continue
            map[neighbour] = map[voxel]
            push!(next, neighbour)
        end
        front = next
    end
    map
end;

"""
The coronary vessels a wall wraps, as one segmenting a patient's myocardium would include them: from each vessel voxel
(lumen or vessel wall) of `box`, 26 rays run through the vessel to the first voxel that is not vessel. A vessel's cross-section
(its connected voxels in one axial slice) joins the wall whose myocardium more than half of its rays end in, on average;
surface vessels, wrapped on one side only, stay out. Returns 0 (neither), 1 (LV) or 2 (RV) for every voxel of `box`.
"""
function enclosed_vessels(box; reach_mm = 8.0)
    vessel(label) = label == 0x1a || label == 0x1b || label == 0x1e
    directions = [
        (i, j, k) ./ sqrt(i^2 + j^2 + k^2)
            for i in -1:1, j in -1:1, k in -1:1 if (i, j, k) != (0, 0, 0)
    ]
    steps = round(Int, reach_mm / XCAT_VOXEL_MM)
    sides = (
        CartesianIndex(1, 0),
        CartesianIndex(-1, 0),
        CartesianIndex(0, 1),
        CartesianIndex(0, -1),
    )
    wall = zeros(UInt8, size(box))
    Threads.@threads for k in axes(box, 3)
        # the share of each vessel voxel's rays that end in each wall
        lv_share = zeros(Float32, size(box, 1), size(box, 2))
        rv_share = zeros(Float32, size(box, 1), size(box, 2))
        for j in axes(box, 2), i in axes(box, 1)
            vessel(box[i, j, k]) || continue
            for d in directions
                for s in 1:steps
                    q = round.(Int, (i, j, k) .+ s .* d)
                    checkbounds(Bool, box, q...) || break
                    label = box[q...]
                    vessel(label) && continue
                    label == 0x0f && (lv_share[i, j] += 1 / length(directions))
                    label == 0x10 && (rv_share[i, j] += 1 / length(directions))
                    break
                end
            end
        end
        # one decision per cross-section: flood each connected patch of vessel in this slice
        seen = falses(size(box, 1), size(box, 2))
        patch = CartesianIndex{2}[]
        for j in axes(box, 2), i in axes(box, 1)
            (vessel(box[i, j, k]) && !seen[i, j]) || continue
            empty!(patch)
            stack = [CartesianIndex(i, j)]
            seen[i, j] = true
            while !isempty(stack)
                c = pop!(stack)
                push!(patch, c)
                for side in sides
                    n = c + side
                    checkbounds(Bool, seen, n) || continue
                    (seen[n] || !vessel(box[n, k])) && continue
                    seen[n] = true
                    push!(stack, n)
                end
            end
            lv = sum(lv_share[c] for c in patch) / length(patch)
            rv = sum(rv_share[c] for c in patch) / length(patch)
            lv + rv > 0.5 || continue
            for c in patch
                wall[c, k] = lv >= rv ? 0x01 : 0x02
            end
        end
    end
    return wall
end

"""
Where the LV and RV walls would be, at 0.2 mm: each wall closed by 3 mm, so a vessel lumen in or overlapping it is in it,
and the vessels each wall wraps (`enclosed_vessels`). And the papillary muscles: the LV myocardium that a 10 mm closing of
the LV cavity fills (thin muscle surrounded by blood).
"""
WALL_ENVELOPES = cached("wall_envelopes_enclosed") do
    walls = findall(label -> label == 0x0f || label == 0x10, xcat_labels)
    margin = ceil(Int, 12 / XCAT_VOXEL_MM)
    lower = max.(Tuple(minimum(walls)) .- margin, 1)
    upper = min.(Tuple(maximum(walls)) .+ margin, XCAT_DIMS)
    box = xcat_labels[CartesianIndices(UnitRange.(lower, upper))]
    (
        # the box's first voxel centre in mm, so the masks stay right whatever grid the labels come on
        origin_mm = XCAT_ORIGIN .+ (lower .- 1) .* XCAT_VOXEL_MM,
        lv = V.XCAT.closing(box .== 0x0f, 3 / XCAT_VOXEL_MM),
        rv = V.XCAT.closing(box .== 0x10, 3 / XCAT_VOXEL_MM),
        papillary = V.XCAT.closing(box .== 0x13, 10 / XCAT_VOXEL_MM) .& (box .== 0x0f),
        enclosed = enclosed_vessels(box),
    )
end;

make_roi(labels, envelopes, territories, walls, grid; papillary, epicardial, split_walls, split_territories) -> roi

"""
    make_roi(labels, envelopes, territories, walls, grid; papillary, epicardial, split_walls, split_territories) -> roi

The compartment rule `roi(p)` for the 0.2 mm voxel at `p` (0 = none): the papillary muscles (compartment `papillary`); else
its wall (XCAT's label, or for a vessel lumen or vessel wall the wall envelope it lies in) and its 1 mm voxel's territory.
A coronary artery or vein lumen outside both, on the heart's surface, is the epicardial compartment of the wall beside it
(`walls`): `epicardial` holds their ids (LV, RV; 0 for none).
"""
function make_roi(
        labels,
        envelopes,
        territories,
        walls,
        grid;
        papillary,
        epicardial,
        split_walls,
        split_territories,
    )
    lv = envelopes.lv
    rv = envelopes.rv
    papillary_mask = envelopes.papillary
    return function roi(p)
        voxel = round.(Int, (p .- XCAT_ORIGIN) ./ XCAT_VOXEL_MM) .+ 1
        all(1 .<= voxel .<= XCAT_DIMS) || return 0
        label = labels[voxel...]
        k = round.(Int, (p .- envelopes.origin_mm) ./ XCAT_VOXEL_MM) .+ 1
        inside = all(1 .<= k .<= size(lv))
        inside && papillary_mask[k...] && return papillary

        wall = label == 0x0f ? 1 : label == 0x10 ? 2 : 0
        # coronary artery, vein or vessel wall: in a wall's envelope, or wrapped by a wall
        if wall == 0 && inside && label in (0x1a, 0x1b, 0x1e)
            wall = lv[k...] ? 1 : rv[k...] ? 2 : Int(envelopes.enclosed[k...])
        end
        g = round.(Int, p .- grid.origin) .+ 1
        # a surface coronary artery or vein lumen: the epicardial compartment of the wall beside it
        if wall == 0 && (label == 0x1a || label == 0x1b) && all(1 .<= g .<= size(walls)) && walls[g...] > 0
            return epicardial[split_walls ? Int(walls[g...]) : 1]
        end
        wall == 0 && return 0

        territory = all(1 .<= g .<= size(territories)) ? Int(territories[g...]) : 0
        territory == 0 && return 0

        wall = split_walls ? wall : 1
        return split_territories ? 3 * (wall - 1) + territory : wall
    end
end

roi = make_roi(
    xcat_labels,
    WALL_ENVELOPES,
    territory_map,
    wall_map,
    growth_grid;
    papillary = findfirst(==("Papillary muscles"), ROI_NAMES),
    epicardial = Tuple(
        something(findfirst(==("$w epicardial vessels"), ROI_NAMES), 0) for w in (split_walls ? ("LV", "RV") : ("Myocardium", "Myocardium"))
    ),
    split_walls,
    split_territories,
)

"""
Categorical colours of the compartments: hue for the territory (LAD red, LCx green, RCA blue), dark for the LV and light
for the RV, gold for the papillary muscles, violet for the epicardial vessels (dark beside the LV, light beside the RV).
"""
function roi_colormap()
    lv = ["#c62828", "#2e7d32", "#1565c0"]
    rv = ["#ef9a9a", "#a5d6a7", "#90caf9"]
    walls =
        split_territories ? (split_walls ? [lv; rv] : lv) :
        (split_walls ? ["#c62828", "#90caf9"] : ["#c62828"])
    epicardial = ["#6a1b9a", "#ce93d8"][1:count(endswith(" epicardial vessels"), ROI_NAMES)]
    return Mke.cgrad(
        [:black; walls; "#ffc400"; epicardial],
        length(ROI_NAMES) + 1;
        categorical = true,
    )
end

begin
    "1 mm myocardial voxel number of each growth-grid voxel (0: not myocardium)."
    const TISSUE_VOXEL_NUMBER = V.XCAT.tissue_numbers(myocardium)

    # the imaged segments of the grown trees, for naming XCAT's drawn coronary lumens
    const IMAGED_ARTERY_SEGMENTS = V.XCAT.imaged_segments(arteries)
    const IMAGED_VEIN_SEGMENTS = V.XCAT.imaged_segments(veins)
end;

"""
Which great vessel each vessel voxel belongs to, and where it lies along it, on a 1 mm grid: its tree and the blood (mL)
between it and the heart. The trees grow by connection through XCAT's vessel labels from where they meet the heart: the
aorta from the LV, the pulmonary arteries from the RV, the venae cavae and brachiocephalic veins from the RA, the pulmonary
veins from the LA (`V.XCAT.great_vessel_transit`). Each state's transit is volume over flow (`V.pool_delay`).

The aorta is timed by the body model's volume, not XCAT's. XCAT vmale_50's aorta is about three standard deviations
narrower than adult men's (21–22 mm ascending and 15 mm descending, against 33 ± 4 and 24 ± 3 mm; Wolak et al. 2008), so its
lumen holds 26.7 mL from the valve to the arch's top. The whole-body model that makes the arterial curves holds 60 mL there
(aortic root 25 mL, arch 35 mL). Each aortic voxel keeps its place along XCAT's aorta, and the blood upstream of it is
scaled by that ratio (`V.XCAT.scale_tree`), so the arterial curves and their timing come from one body.
"""
GREAT_VESSEL_TRANSIT = let
    xcat = cached("great_vessel_transit") do
        V.XCAT.great_vessel_transit(xcat_labels, XCAT_ORIGIN, XCAT_VOXEL_MM)
    end
    body = V.circulation(V.HEART_TRANSPORT)
    volume(segment) = body.segs[findfirst(s -> s.name === segment, body.segs)].v
    V.XCAT.scale_tree(xcat, 0x01, (volume(:aortic_root) + volume(:arch)) / xcat.apex_ml)
end;

"The heart's phantom spec: the body pools, with the descending aorta timed from this subject's arch."
PHANTOM_SPEC = let pools = copy(V.HEART_PHANTOM.pools)
    i = findfirst(p -> p.name == "descending aorta", pools)
    d = pools[i]
    pools[i] = V.Pool(
        d.name,
        d.material,
        d.segment,
        d.what;
        through = d.through,
        lead_ml = GREAT_VESSEL_TRANSIT.apex_ml,
        lead_through = d.lead_through,
    )
    V.PhantomSpec(V.HEART_PHANTOM; pools)
end;

make_classifier(labels, grid, numbers, artery_segments, vein_segments, transit) -> classify

"""
    make_classifier(labels, grid, numbers, artery_segments, vein_segments, transit) -> classify

The classifier `classify(p)`: what is at the 0.2 mm voxel centre `p`. Myocardium is its 1 mm voxel, a coronary lumen its
imaged segment, a great-vessel voxel its vessel's pool and its place along the vessel (`transit`); every other XCAT
label is decided by the label alone (a body pool, `V.XCAT_POOL_OF_LABEL`, or a static material). Everything it reads is
captured with its type, so it compiles to fast code.
"""
function make_classifier(
        labels,
        grid,
        numbers,
        artery_segments,
        vein_segments,
        transit,
    )
    # by the label alone: 0 needs the voxel's geometry
    by_label = fill((V.VOX_POOL, Int32(24), 0.0f0), 256)
    needs_geometry = falses(256)
    for label in 0x00:0xff
        if label in (0x0f, 0x10, 0x1a, 0x1b, 0x18, 0x19, 0x1c, 0x7c, 0x7d)
            needs_geometry[label + 1] = true
        elseif haskey(V.XCAT_POOL_OF_LABEL, label)
            by_label[label + 1] = (V.VOX_POOL, Int32(V.XCAT_POOL_OF_LABEL[label]), 0.0f0)
        elseif haskey(V.XCAT_STATIC, Int(label))
            by_label[label + 1] = (V.VOX_STATIC, Int32(label), 0.0f0)
        end
    end
    vessel_pool(p, t) = V.XCAT.great_vessel_pool(transit, p, t)

    return function classify(p)
        voxel = round.(Int, (p .- XCAT_ORIGIN) ./ XCAT_VOXEL_MM) .+ 1
        all(1 .<= voxel .<= XCAT_DIMS) || return (V.VOX_STATIC, Int32(0), 0.0f0)
        label = labels[voxel...]
        needs_geometry[label + 1] || return by_label[label + 1]

        if label == 0x0f || label == 0x10
            number = V.XCAT.nearest_tissue_voxel(numbers, round.(Int, p .- grid.origin) .+ 1)
            return number > 0 ?
                (V.VOX_TISSUE, number, 0.0f0) :
                (V.VOX_POOL, Int32(26), 0.0f0)
        end
        label == 0x1a && return (
            V.VOX_VESSEL,
            Int32(first(V.XCAT.nearest_segment(artery_segments, p))),
            0.0f0,
        )
        label == 0x1b && return (
            V.VOX_VESSEL,
            -Int32(first(V.XCAT.nearest_segment(vein_segments, p))),
            0.0f0,
        )
        label == 0x1c && return vessel_pool(p, 0x01)
        label == 0x19 && return vessel_pool(p, V.XCAT.vessel_tree_at(transit, p, (0x01, 0x02)))
        label == 0x7d && return vessel_pool(p, 0x02)
        label == 0x7c && return vessel_pool(p, 0x05)
        return vessel_pool(p, V.XCAT.vessel_tree_at(transit, p, (0x03, 0x04, 0x05, 0x06, 0x07)))
    end
end

classify = make_classifier(
    xcat_labels,
    growth_grid,
    TISSUE_VOXEL_NUMBER,
    IMAGED_ARTERY_SEGMENTS,
    IMAGED_VEIN_SEGMENTS,
    GREAT_VESSEL_TRANSIT,
)

"""
The bolus-tracking points: the descending aorta at the myocardium's mid level, and the ascending aorta 15 mm above where it
leaves the heart (clear of the root's sinuses), the centre of its cross-section there.
"""
TRACKING_POINTS = let voxels = findall(myocardium)
    z_mid =
        (minimum(v[3] for v in voxels) + maximum(v[3] for v in voxels)) / 2 +
        growth_grid.origin[3] - 1
    k = xcat_slice(z_mid)
    aorta = [
        (i, j)
            for j in 1:XCAT_DIMS[2], i in 1:XCAT_DIMS[1]
            if xcat_labels[i, j, k] == 0x1c
    ]
    back = maximum(last, aorta)
    descending = [
        (i, j)
            for (i, j) in aorta
            if j > back - 15 / XCAT_VOXEL_MM # within 15 mm of its back
    ]
    descending_point = (
        XCAT_ORIGIN[1] + (sum(first, descending) / length(descending) - 1) * XCAT_VOXEL_MM,
        XCAT_ORIGIN[2] + (sum(last, descending) / length(descending) - 1) * XCAT_VOXEL_MM,
        XCAT_ORIGIN[3] + (k - 1) * XCAT_VOXEL_MM,
    )

    # the ascending aorta: the aortic tree's cells before the arch's top, 15 mm above the lowest of them
    transit = GREAT_VESSEL_TRANSIT
    step = round(Int, transit.h / XCAT_VOXEL_MM)
    ascending = [
        c
            for c in CartesianIndices(transit.tree)
            if transit.tree[c] == 0x01 && !transit.descending[c] && !isnan(transit.volume[c]) &&
            transit.volume[c] < transit.apex_ml && xcat_labels[((Tuple(c) .- 1) .* step .+ 1)...] == 0x1c
    ]
    z_low = minimum(c[3] for c in ascending)
    level = [c for c in ascending if c[3] == z_low + round(Int, 15 / transit.h)]
    ascending_point = Tuple(
        transit.origin_mm[d] + (sum(c[d] for c in level) / length(level) - 1) * transit.h for d in 1:3
    )
    ["descending aorta" => descending_point, "ascending aorta" => ascending_point]
end

"""
Write one scenario's phantom into `PHANTOM_DIR/name`, once, from its transport file (read here, and let go once the
phantom is written).
"""
function save_phantom(
        name,
        transport_file;
        spec = PHANTOM_SPEC,
        classify = classify,
    )
    dir = joinpath(PHANTOM_DIR, name)
    isfile(joinpath(dir, "phantom.toml")) && return dir

    scans = deserialize(transport_file)
    seconds = @elapsed V.phantom(
        spec,
        scans,
        classify,
        V.VoxelGrid(XCAT_ORIGIN, XCAT_VOXEL_MM, XCAT_DIMS);
        pairs = flows.pairs,
        terr = flows.terr,
        nvox = size(tissue_centres, 2),
        out = dir,
        roi = roi,
        roi_names = ROI_NAMES,
        roi_groups = ROI_GROUPS,
        tracking = TRACKING_POINTS,
        facts = ["organ" => "heart (XCAT vmale_50, 0.2 mm)", "model" => "Vasculature.jl"],
        verbose = true,
    )
    @info "$name: $(round(seconds; digits = 1)) s"
    scans = nothing
    GC.gc()
    return dir
end

HEALTHY_DIR = save_phantom("healthy", TRANSPORT);

"""
The heart's box of a saved phantom, the only part read back: the growth grid's (the pericardium's) box and the
descending aorta's tracking point, with 15 mm to spare. Returns the bundle (labels and
compartments) and its first z (mm). Ground-truth MBF is each compartment's own (its beds' flow over the mass of
everything in it, tissue and vessels), from the bundle's tables, so the `mbf_<state>.f32.zst` volumes are not read back.
"""
function inspect(dir)
    tracking = V.DynamicPhantom.load(dir; volume = false).meta["tracking"]
    aorta_mm = tracking["descending aorta"]["point_mm"]
    heart_lower_mm = collect(growth_grid.origin)
    heart_upper_mm = heart_lower_mm .+ collect(growth_grid.dims) .- 1
    lower_mm = [min.(heart_lower_mm[1:2], aorta_mm[1:2]); heart_lower_mm[3]] .- 15
    upper_mm = [max.(heart_upper_mm[1:2], aorta_mm[1:2]); heart_upper_mm[3]] .+ 15
    function voxel(mm, d)
        number = round(Int, (mm - XCAT_ORIGIN[d]) / XCAT_VOXEL_MM) + 1
        return clamp(number, 1, XCAT_DIMS[d])
    end
    bundle = V.DynamicPhantom.load(
        dir;
        rois = true,
        zrange = voxel(lower_mm[3], 3):voxel(upper_mm[3], 3),
        xrange = voxel(lower_mm[1], 1):voxel(upper_mm[1], 1),
        yrange = voxel(lower_mm[2], 2):voxel(upper_mm[2], 2),
    )
    return (bundle = bundle, z0 = bundle.origin_mm[3])
end

healthy = inspect(HEALTHY_DIR);

"""
One phantom slice in six panels, each with its own scale underneath: labels by kind, compartments and ground-truth MBF in
`state`; the idealized 70 keV image and iodine at time `t`; ground-truth MBF in the other state.
"""
function slice_figure(view, state, t, z)
    bundle = view.bundle
    k = clamp(round(Int, (z - view.z0) / XCAT_VOXEL_MM) + 1, 1, size(bundle.labels, 3))
    is, js = axes(bundle.labels, 1), axes(bundle.labels, 2)
    xs = bundle.origin_mm[1] .+ (is .- 1) .* XCAT_VOXEL_MM
    ys = bundle.origin_mm[2] .+ (js .- 1) .* XCAT_VOXEL_MM
    labels = bundle.labels[is, js, k]
    other = state == "rest" ? "stress" : "rest"

    # per-label lookups: kind, idealized HU and iodine at time t
    kind_colours = [:gray30, :gray65, :mediumpurple, :darkorange, :red, :deepskyblue]
    kind_names = [
        "static (air, cartilage)",
        "body tissue",
        "blood",
        "myocardium",
        "artery",
        "vein",
    ]
    function kind_code(row)
        row.kind == "static" && return 1
        row.kind == "pool" && return occursin("blood", row.base) ? 3 : 2
        return Dict("tissue" => 4, "artery" => 5, "vein" => 6)[row.kind]
    end
    kind_of = zeros(UInt8, 65536)
    for (label, row) in bundle.table
        kind_of[label + 1] = label == 0 ? 0 : kind_code(row)
    end
    # each compartment's ground truth over all of it, tissue and vessels: its beds' flow over the mass of everything in it
    mbf = Dict(st => zeros(Float32, length(ROI_NAMES) + 1) for st in ("rest", "stress"))
    for st in ("rest", "stress"), (r, name) in enumerate(ROI_NAMES)
        mbf[st][r + 1] = bundle.roi_table[(st, name)].mbf
    end
    rois = bundle.rois[is, js, k]
    hu = fill(-1000.0f0, 65536)
    for (label, value) in V.DynamicPhantom.hu(bundle, state, t)
        hu[label + 1] = value
    end
    iodine = zeros(Float32, 65536)
    for (label, value) in V.DynamicPhantom.iodine(bundle, state, t)
        iodine[label + 1] = value
    end

    fig = Mke.Figure(; size = (1500, 1120))
    function panel(
            row,
            column,
            title,
            image;
            colormap,
            colorrange,
            label = nothing,
            legend = nothing,
        )
        cell = fig[row, column] = Mke.GridLayout()
        ax = Mke.Axis(
            cell[1, 1];
            title,
            aspect = Mke.DataAspect(),
            yreversed = true,
        )
        map = Mke.heatmap!(
            ax,
            xs,
            ys,
            image;
            colormap,
            colorrange,
        )
        Mke.hidedecorations!(ax)
        if legend === nothing
            Mke.Colorbar(
                cell[2, 1],
                map;
                label,
                vertical = false,
                flipaxis = false,
                width = Mke.Relative(0.9),
                tellwidth = false,
            )
        else
            colours, names = legend
            Mke.Legend(
                cell[2, 1],
                [Mke.PolyElement(; color = c) for c in colours],
                names;
                orientation = :horizontal,
                nbanks = 3,
                framevisible = false,
                tellwidth = false,
            )
        end
        return ax
    end

    panel(
        1,
        1,
        "Labels by kind",
        [kind_of[l + 1] for l in labels];
        colormap = Mke.cgrad([:black; kind_colours], 7; categorical = true),
        colorrange = (-0.5, 6.5),
        legend = (kind_colours, kind_names),
    )
    panel(
        1,
        2,
        "Compartments",
        rois;
        colormap = roi_colormap(),
        colorrange = (-0.5, length(ROI_NAMES) + 0.5),
        legend = (roi_colormap().colors[2:end], ROI_NAMES),
    )
    panel(
        1,
        3,
        "Ground-truth MBF by compartment, $state",
        [mbf[state][r + 1] for r in rois];
        colormap = :turbo,
        colorrange = (0, 6),
        label = "mL/min/g",
    )
    panel(
        2,
        1,
        "Idealized 70 keV, $state, t = $t s",
        [hu[l + 1] for l in labels];
        colormap = :grays,
        colorrange = (-200, 600),
        label = "HU",
    )
    panel(
        2,
        2,
        "Iodine, $state, t = $t s",
        [iodine[l + 1] for l in labels];
        colormap = :turbo,
        colorrange = (0, 8),
        label = "mg/mL",
    )
    panel(
        2,
        3,
        "Ground-truth MBF by compartment, $other",
        [mbf[other][r + 1] for r in rois];
        colormap = :turbo,
        colorrange = (0, 6),
        label = "mL/min/g",
    )
    Mke.Label(
        fig[0, :],
        "z = $(round(z; digits = 1)) mm";
        fontsize = 20,
        font = :bold,
    )
    # three equal columns: no legend or colour bar widens its own
    for column in 1:3
        Mke.colsize!(fig.layout, column, Mke.Relative(1 / 3))
    end
    return fig
end

"Each compartment's iodine, with the tracking arrival and peak, and the heart's iodine ledger at stress."
function curves_figure(view)
    bundle = view.bundle
    fig = Mke.Figure(size = (1500, 540))
    for (column, state) in enumerate(("rest", "stress"))
        ax = Mke.Axis(
            fig[1, column];
            title = "Compartments, $state",
            xlabel = "time after injection (s)",
            ylabel = "iodine (mg)",
        )
        for name in ROI_NAMES
            row = bundle.roi_table[(state, name)]
            Mke.lines!(
                ax,
                bundle.times,
                row.mass;
                linewidth = 2,
                label = @sprintf("%s · %.2f mL/min/g", name, row.mbf),
            )
        end
        timing = V.DynamicPhantom.scan_times(bundle; state, roi = "roi:descending aorta")
        Mke.vlines!(
            ax,
            [timing.arrival, timing.peak_time];
            color = :gray40,
            linestyle = :dash,
        )
        Mke.xlims!(ax, 0, 40)
        column == 1 && Mke.axislegend(
            ax;
            position = :lt,
            labelsize = 11,
        )
    end

    ledger = bundle.ledger["stress"]
    ax = Mke.Axis(
        fig[1, 3];
        title = "Iodine ledger, stress",
        xlabel = "time after injection (s)",
        ylabel = "iodine (mg)",
    )
    entries = [
        ("delivered_mg", :black, "delivered", :solid),
        ("held_mg", :darkorange, "held", :solid),
        ("outflow_mg", :royalblue, "out through the veins", :solid),
        ("painted_mg", :red, "painted", :dash),
    ]
    for (column, colour, name, style) in entries
        Mke.lines!(
            ax,
            ledger["t_s"],
            ledger[column];
            color = colour,
            linewidth = 2,
            linestyle = style,
            label = name,
        )
    end
    Mke.xlims!(ax, 0, 40)
    Mke.axislegend(ax; position = :lt)
    return fig
end

"""
The share of the delivered iodine that has already left the heart through the coronary sinus, 1, 2 and 3 s after
arrival, 1 s before the arterial peak and at the peak.
"""
function outflow_check(view)
    bundle = view.bundle
    function share_out(ledger, t)
        ts = ledger["t_s"]
        k = clamp(searchsortedlast(ts, t), 1, length(ts) - 1)
        w = (t - ts[k]) / (ts[k + 1] - ts[k])
        out = (1 - w) * ledger["outflow_mg"][k] + w * ledger["outflow_mg"][k + 1]
        delivered = (1 - w) * ledger["delivered_mg"][k] + w * ledger["delivered_mg"][k + 1]
        return out / max(delivered, eps())
    end
    rows = map(("rest", "stress")) do state
        ledger = bundle.ledger[state]
        timing = V.DynamicPhantom.scan_times(bundle; state, roi = "roi:descending aorta")
        times = [timing.arrival .+ [1, 2, 3]; timing.peak_time - 1; timing.peak_time]
        shares = [@sprintf("%.0f%%", 100 * share_out(ledger, t)) for t in times]
        @sprintf(
            "| **%s** (arrival at %.1f s, peak at %.1f s) | %s |",
            state,
            timing.arrival,
            timing.peak_time,
            join(shares, " | "),
        )
    end
    text = """
    !!! note "Iodine already out of the heart"
        | | arrival + 1 s | arrival + 2 s | arrival + 3 s | 1 s before the peak | arterial peak |
        |---|---|---|---|---|---|
    $(join("    " .* rows, "\n"))

        A measurement that counts the iodine the heart has taken up misses this share.
    """
    return Markdown.parse(text)
end

const DIAMETER_STENOSIS = 0.9

"The lesion along the LAD's first 10 mm of imaged path: its segments, their imaged radii and lengths, and `circuit`'s list."
LESION = let tree = arteries, n = length(tree.r)
    children = [Int[] for _ in 1:n]
    for j in 1:n
        tree.parent[j] > 0 && push!(children[tree.parent[j]], j)
    end
    segment_length(j) = sqrt(sum(abs2, tree.b[:, j] .- tree.a[:, j]))

    path = Int[]
    covered = 0.0f0
    j = coronary_branch.lad
    while covered < 10.0f0
        push!(path, j)
        covered += segment_length(j)
        imaged = [c for c in children[j] if tree.isbase[c]]
        isempty(imaged) && break
        j = imaged[argmax([segment_length(c) * tree.r[c]^2 for c in imaged])]
    end

    # the imaged lumen radius at each segment: XCAT's LAD centreline radius nearest to it
    lad = V.XCAT.centerline(xcat_surfaces["dias_lad1"]; step = 0.5)
    imaged_radius(j) = lad.r[
        argmin(vec(sum(abs2, lad.X .- (tree.a[:, j] .+ tree.b[:, j]) ./ 2; dims = 1))),
    ]
    radii = Float32[imaged_radius(j) for j in path]
    lengths = Float32[segment_length(j) for j in path]
    area_fraction(j, r) = min(1.0f0, (r * (1 - DIAMETER_STENOSIS) / tree.r[j])^2)
    (
        path = path,
        radii = radii,
        stenosis = [(j, area_fraction(j, r), l) for (j, r, l) in zip(path, radii, lengths)],
    )
end

flows_lad90 = cached("flow_lad90") do
    stress = V.Flow.circuit(
        arteries,
        flows.terr;
        state = :stress,
        veins,
        vleaf = flows.vleaf,
        reference = flows.stress,
        kappa = flows.stress.resist.kappa,
        stenosis = LESION.stenosis,
        phys = flows.stress.phys,
        to = DEVICE,
    )
    rest = V.Flow.circuit(
        arteries,
        flows.terr;
        state = :rest,
        veins,
        vleaf = flows.vleaf,
        reference = flows.rest,
        kappa = flows.stress.resist.kappa,
        stenosis = LESION.stenosis,
        phys = flows.rest.phys,
        demand = flows.demand,
        to = DEVICE,
    )
    (rest = rest, stress = stress)
end;

TRANSPORT_LAD90 = transport_cache("transport_lad90", flows_lad90);

begin
    "Plaque materials, as static labels: fibrous tissue (collagen), a lipid core, calcification (hydroxyapatite)."
    const PLAQUE = Dict(40 => "basis_collagen", 41 => "basis_lipid", 42 => "hydroxyapatite")
    const STENOSIS_SPEC = V.PhantomSpec(
        PHANTOM_SPEC;
        materials = merge(V.XCAT_STATIC, PLAQUE),
    )

    "Two unit vectors perpendicular to `axis` and to each other."
    function perpendiculars(axis)
        a = axis ./ sqrt(sum(abs2, axis))
        e = abs(a[3]) < 0.9 ? [0.0, 0.0, 1.0] : [1.0, 0.0, 0.0]
        u = e .- sum(e .* a) .* a
        u ./= sqrt(sum(abs2, u))
        v = [
            a[2] * u[3] - a[3] * u[2],
            a[3] * u[1] - a[1] * u[3],
            a[1] * u[2] - a[2] * u[1],
        ]
        return u, v
    end

    """
    The plaque at `p` inside XCAT's LAD lumen along the lesion: calcium in a 60° arc at the outer rim, a lipid core behind the
    cap on the opposite side, fibrous tissue elsewhere; `nothing` where `p` is lumen or outside the lesion.
    """
    function plaque(p)
        for (j, r) in zip(LESION.path, LESION.radii)
            a = arteries.a[:, j]
            b = arteries.b[:, j]
            axis = b .- a
            s = clamp(sum((p .- a) .* axis) / sum(abs2, axis), 0, 1)
            offset = p .- (a .+ s .* axis)
            distance = sqrt(sum(abs2, offset))
            (0 < s < 1 && distance <= r + 0.2) || continue

            lumen = r * (1 - DIAMETER_STENOSIS)
            distance <= lumen && return nothing
            u, v = perpendiculars(axis)
            angle = atan(sum(offset .* v), sum(offset .* u))
            # 0 at the lumen, 1 at the vessel wall
            depth = (distance - lumen) / max(r - lumen, 1.0e-6)
            abs(angle) < π / 6 && depth > 0.6 && return 42
            abs(abs(angle) - π) < π / 3 && 0.25 < depth < 0.75 && return 41
            return 40
        end
        return nothing
    end

    "The box (mm) around the lesion's segments and their imaged radius: plaque is only looked for inside it."
    const LESION_BOX = let
        ends = reduce(hcat, [hcat(arteries.a[:, j], arteries.b[:, j]) for j in LESION.path])
        reach = maximum(LESION.radii) + 0.2
        (vec(minimum(ends; dims = 2)) .- reach, vec(maximum(ends; dims = 2)) .+ reach)
    end

    """
    The healthy classifier `healthy`, with plaque where the lesion narrows XCAT's LAD lumen. A factory, like the healthy
    one: the returned closure holds `healthy` and the lesion's box as concrete values, so the billions of voxel calls
    are compiled for them rather than looked up as untyped notebook globals each time.
    """
    function make_plaque_classifier(healthy, box)
        lower, upper = box
        return function (p)
            what = healthy(p)
            if what[1] == V.VOX_VESSEL && what[2] > 0 && all(lower .<= p .<= upper)
                material = plaque(collect(Float64, p))
                material !== nothing && return (V.VOX_STATIC, Int32(material), 0.0f0)
            end
            return what
        end
    end

    classify_with_plaque = make_plaque_classifier(classify, LESION_BOX)
end;

STENOSIS_DIR = save_phantom(
    "lad_stenosis_90",
    TRANSPORT_LAD90;
    spec = STENOSIS_SPEC,
    classify = classify_with_plaque,
);

stenosis = inspect(STENOSIS_DIR);