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;

using PlutoUI: Slider, Select, TableOfContents

using Printf: @sprintf

using Serialization: serialize, deserialize

using Mmap: mmap

using Statistics: median, quantile

import TOML

import Vasculature as V

import CairoMakie as Mke

import GPUSelect

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", "pulmonary_pipeline_cache")
    const PHANTOM_DIR = joinpath(homedir(), "xcat-perfused", "pulmonary_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

"The six face neighbours of a voxel."
const SIX_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),
)

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"))

    # XCAT's surface meshes lie this far (mm) from its NURBS and label frame: the shift the chest volume was built with
    # (its lungs land on XCAT's own lung labels, section 2)
    const MESH_OFFSET = (131.4, 238.4, 193.8)
    const ARTERY_TUBES = ("right_pulmonary_arts", "left_pulmonary_arts")
    const VEIN_TUBES = ("right_pulmonary_veins", "left_pulmonary_veins")
    lung_meshes = V.XCAT.read_mesh_blocks(
        joinpath(XCAT_DIR, "vmale_50_heart_1_organs.raw");
        names = ("rlung", "llung"),
        offset = MESH_OFFSET,
    )
    vessel_tubes = V.XCAT.read_mesh_blocks(
        joinpath(XCAT_DIR, "vmale_50_heart_1_vessels.raw");
        names = (ARTERY_TUBES..., VEIN_TUBES...),
        offset = MESH_OFFSET,
    )
end;

begin
    "XCAT label groups for pictures: (name, colour, XCAT labels)."
    const LABEL_GROUPS = [
        ("lung", :gray40, (6,)),
        ("airways", :paleturquoise, (7, 11)),
        ("arteries", :red, (25, 28, 125)),
        ("veins", :royalblue, (24, 124)),
        ("heart blood", :mediumpurple, (19, 20, 21, 22)),
        ("myocardium", :darkorange, (15, 16, 17, 18)),
        ("fat", :khaki, (1,)),
        ("muscle", :rosybrown, (4, 13)),
        ("bone", :ivory, (5, 8, 9)),
        ("other", :gray65, ()),
    ]

    "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

    "The XCAT label at `p` (mm); 0 outside the chest."
    function xcat_label(p)
        k = round.(Int, (p .- XCAT_ORIGIN) ./ XCAT_VOXEL_MM) .+ 1
        return all(1 .<= k .<= XCAT_DIMS) ? xcat_labels[k...] : 0x00
    end
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 chest 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 = :fitzoom, # zoomed to the data, not to the box it could rotate in
        protrusions = 0,
    )
    Mke.hidedecorations!(ax)
    Mke.hidespines!(ax)
    return ax
end

growth_grid = let points = hcat(first(only(lung_meshes["rlung"])), first(only(lung_meshes["llung"])))
    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

"The lung: the 1 mm cells that are mostly XCAT's lung (label 6), the voxels the phantom paints as lung."
LUNG = cached("lung") do
    V.XCAT.label_mask(xcat_labels, XCAT_ORIGIN, XCAT_VOXEL_MM, growth_grid, (0x06,))
end;

"""
Stub spacing and length, from human pulmonary morphometry (Huang et al. 1996, `test/reference`): the segment length (between
two branchings) of a vessel of the imaged vessels' length-weighted median diameter, and of one two orders smaller, the
typical side branch, as the stub's length.
"""
STUB_LENGTHS = let
    huang = TOML.parsefile(
        joinpath(pkgdir(V), "test", "reference", "huang1996_human_pulmonary_morphometry.toml"),
    )["data"]
    # segment length (mm) at a diameter (mm) `below` orders smaller, log-linear between orders
    function segment_length(table, diameter; below = 0)
        d = log.(table["diameter_mean_mm"])
        l = log.(table["length_mean_mm"])
        x = clamp(log(diameter), d[1], d[end])
        k = clamp(searchsortedlast(d, x), 1, length(d) - 1)
        order = clamp(k + (x - d[k]) / (d[k + 1] - d[k]) - below, 1, length(d))
        j = clamp(floor(Int, order), 1, length(d) - 1)
        w = order - j
        return exp((1 - w) * l[j] + w * l[j + 1])
    end
    function median_diameter(names)
        lines = [V.XCAT.tube_centerline(points) for name in names for (points, _) in vessel_tubes[name]]
        lengths = [sum(sqrt.(sum(abs2, diff(line.X; dims = 2); dims = 1))) for line in lines]
        diameters = [2 * line.r[1] for line in lines]
        order = sortperm(diameters)
        covered = cumsum(lengths[order])
        return diameters[order][findfirst(>=(covered[end] / 2), covered)]
    end
    artery = median_diameter(ARTERY_TUBES)
    vein = median_diameter(VEIN_TUBES)
    (
        artery = (
            diameter = artery,
            spacing = segment_length(huang["arterial_segments"], artery),
            depth = segment_length(huang["arterial_segments"], artery; below = 2),
        ),
        vein = (
            diameter = vein,
            spacing = segment_length(huang["venous_segments"], vein),
            depth = segment_length(huang["venous_segments"], vein; below = 2),
        ),
    )
end

"""
The arterial inlet tree: XCAT's main pulmonary artery (its NURBS surface, from the valve), with the lung artery tubes joined
into trees and grafted onto it (`V.XCAT.graft`). Centrelines are sampled at the stub spacing.
"""
artery_base = cached("artery_base") do
    step = STUB_LENGTHS.artery.spacing
    trunk = V.XCAT.centerline(xcat_surfaces["dias_pul_art"]; step)
    hub = vec(sum(trunk.X; dims = 2)) ./ size(trunk.X, 2)
    lines = [V.XCAT.tube_centerline(points; name, step) for name in ARTERY_TUBES for (points, _) in vessel_tubes[name]]
    V.XCAT.graft(trunk, V.XCAT.tube_forest(lines; toward = hub); name = "pulmonary_arteries")
end;

"""
The venous inlet trees, one per lung: its vein tubes joined into trees, the largest taking the others, carried straight to
the nearest blood of the left atrium (XCAT's label, so the veins end where the atrium's blood begins; `V.XCAT.drain`). The
ostium's radius is the lung's two pulmonary veins (Huang 1996 venous order 15, 12.97 mm) as one cube-law equivalent.
"""
vein_bases = cached("vein_bases") do
    step = STUB_LENGTHS.vein.spacing
    points, _ = V.XCAT.tessellate(xcat_surfaces["dias_la"]; nu = 60, nv = 60)
    grid = V.grid_around(points, 1.0; margin = 10)
    atrium = V.XCAT.label_mask(xcat_labels, XCAT_ORIGIN, XCAT_VOXEL_MM, grid, (0x15,))
    wall = [c for c in findall(atrium) if any(s -> !checkbounds(Bool, atrium, c + s) || !atrium[c + s], SIX_STEPS)]
    surface = reduce(hcat, [collect(V.center(grid, Tuple(c)...)) for c in wall])
    hub = vec(sum(surface; dims = 2)) ./ size(surface, 2)
    [
        V.XCAT.drain(
            V.XCAT.tube_forest(
                [V.XCAT.tube_centerline(points; name, step) for (points, _) in vessel_tubes[name]];
                toward = hub,
            ),
            surface;
            root_radius = cbrt(2 * (12.97 / 2)^3),
            name,
        ) for name in VEIN_TUBES
    ]
end;

"""
The airway tree, from XCAT's airway labels (wall 7 and air 11) on a 0.6 mm grid: its distance from the trachea's top through
the airway, cut into 3 mm shells; each shell's connected pieces are the tree's nodes, each joined to the piece of the shell
before it that it touches. Returns the grid, the airway mask, each node's cells, parent and shell, and each cell's node.
"""
function airway_tree(; cell = 0.6, shell_mm = 3.0)
    grid = V.VoxelGrid(growth_grid.origin, cell, Tuple(ceil.(Int, growth_grid.dims ./ cell)))
    airway = V.XCAT.label_mask(xcat_labels, XCAT_ORIGIN, XCAT_VOXEL_MM, grid, (0x07, 0x0b); fraction = 1.0e-9)
    steps = [CartesianIndex(i, j, k) for i in -1:1, j in -1:1, k in -1:1 if (i, j, k) != (0, 0, 0)]
    # keep the largest connected piece
    piece = zeros(Int32, size(airway))
    sizes = Int[]
    for seed in findall(airway)
        piece[seed] == 0 || continue
        push!(sizes, 0)
        stack = [seed]
        piece[seed] = length(sizes)
        while !isempty(stack)
            c = pop!(stack)
            sizes[end] += 1
            for s in steps
                n = c + s
                checkbounds(Bool, airway, n) && airway[n] && piece[n] == 0 || continue
                piece[n] = length(sizes)
                push!(stack, n)
            end
        end
    end
    airway .= piece .== argmax(sizes)
    # distance (cells) from the trachea's top, through the airway
    top = maximum(c[3] for c in findall(airway))
    distance = fill(Int32(-1), size(airway))
    queue = [c for c in findall(airway) if c[3] >= top - 1]
    foreach(c -> distance[c] = 0, queue)
    head = 1
    while head <= length(queue)
        c = queue[head]
        head += 1
        for s in steps
            n = c + s
            checkbounds(Bool, airway, n) && airway[n] && distance[n] < 0 || continue
            distance[n] = distance[c] + 1
            push!(queue, n)
        end
    end
    shell_of(c) = floor(Int, distance[c] * cell / shell_mm)
    # nodes: the connected pieces of each shell
    node = zeros(Int32, size(airway))
    members = Vector{Vector{CartesianIndex{3}}}()
    shell = Int[]
    for seed in queue
        node[seed] == 0 || continue
        push!(members, CartesianIndex{3}[])
        push!(shell, shell_of(seed))
        stack = [seed]
        node[seed] = length(members)
        while !isempty(stack)
            c = pop!(stack)
            push!(members[end], c)
            for s in steps
                n = c + s
                checkbounds(Bool, airway, n) && airway[n] && node[n] == 0 && shell_of(n) == shell[end] || continue
                node[n] = length(members)
                push!(stack, n)
            end
        end
    end
    parent = zeros(Int, length(members))
    for (id, cells) in enumerate(members)
        for c in cells
            k = findfirst(s -> checkbounds(Bool, airway, c + s) && airway[c + s] && shell[node[c + s]] == shell[id] - 1, steps)
            k === nothing || (parent[id] = node[c + steps[k]]; break)
        end
    end
    return (; grid, airway, members, parent, shell)
end

AIRWAYS = cached("airways") do
    airway_tree()
end;

"""
The five lobar bronchi, as nodes of the airway tree. Walking down from the trachea, a split is the first node whose second
largest branch holds at least `share` of the airway below it. The carina's two branches are the main bronchi (the right at
smaller x). The right main bronchus's first split gives off the upper lobe bronchus, the branch that goes up; the bronchus
intermedius then splits into the middle lobe bronchus (the anterior one) and the lower lobe bronchus. The left main bronchus
splits into the upper and lower lobe bronchi.
"""
LOBAR_BRONCHI = let tree = AIRWAYS, n = length(tree.members)
    children = [Int[] for _ in 1:n]
    foreach(i -> tree.parent[i] > 0 && push!(children[tree.parent[i]], i), 1:n)
    # airway cells below each node, and their summed positions
    below = Float64.(length.(tree.members))
    position = zeros(3, n)
    for i in 1:n
        position[:, i] = sum(collect(Tuple(c)) for c in tree.members[i])
    end
    for i in sortperm(tree.shell; rev = true)
        p = tree.parent[i]
        p > 0 || continue
        below[p] += below[i]
        position[:, p] .+= position[:, i]
    end
    centre(i) = collect(tree.grid.origin) .+ (position[:, i] ./ below[i] .- 1) .* tree.grid.h
    function first_split(i, share)
        while true
            branches = sort(children[i]; by = c -> below[c], rev = true)
            isempty(branches) && error("no airway split below node $i")
            length(branches) >= 2 && below[branches[2]] >= share * below[i] && return branches[1:2]
            i = branches[1]
        end
        return
    end
    # a split's two branches ordered by one coordinate of their centres (1 x, 2 y front to back, 3 z up)
    ordered(pair, axis; rev = false) = pair[sortperm([centre(i)[axis] for i in pair]; rev)]
    trachea = argmax([tree.parent[i] == 0 ? below[i] : 0.0 for i in 1:n])
    right, left = ordered(first_split(trachea, 0.25), 1)
    upper, intermedius = ordered(first_split(right, 0.15), 3; rev = true)
    middle, lower = ordered(first_split(intermedius, 0.15), 2)
    left_upper, left_lower = ordered(first_split(left, 0.15), 3; rev = true)
    (
        names = ["right upper", "right middle", "right lower", "left upper", "left lower"],
        nodes = [upper, middle, lower, left_upper, left_lower],
        children = children,
    )
end;

"""
Each 1 mm lung voxel's lobe (1-5, `LOBAR_BRONCHI.names`; 0 outside the lung): the lobe of the nearest lobar airway, by the
distance through the lung (6-connected, all lobes spreading at once from their own airways).
"""
LOBES = cached("lobes") do
    lobe = zeros(Int8, growth_grid.dims)
    for (k, node) in enumerate(LOBAR_BRONCHI.nodes)
        stack = [node]
        while !isempty(stack)
            i = pop!(stack)
            append!(stack, LOBAR_BRONCHI.children[i])
            for c in AIRWAYS.members[i]
                p = AIRWAYS.grid.origin .+ (Tuple(c) .- 1) .* AIRWAYS.grid.h
                q = round.(Int, (p .- growth_grid.origin) ./ growth_grid.h) .+ 1
                all(1 .<= q .<= growth_grid.dims) && (lobe[q...] = k)
            end
        end
    end
    queue = findall(!=(0), lobe)
    head = 1
    while head <= length(queue)
        c = queue[head]
        head += 1
        for s in SIX_STEPS
            n = c + s
            checkbounds(Bool, lobe, n) && lobe[n] == 0 && LUNG[n] || continue
            lobe[n] = lobe[c]
            push!(queue, n)
        end
    end
    lobe .* LUNG
end;

fissures(lobes, airways, grid; hilum_mm = 5.0) -> BitArray

"""
    fissures(lobes, airways, grid; hilum_mm = 5.0) -> BitArray

The fissures on `grid`: lung voxels of `lobes` beside a voxel of another lobe, except within `hilum_mm` of an airway
(`airways`, as `airway_tree` returns it), where the lobar vessels and bronchi enter the lobes.
"""
function fissures(lobes, airways, grid; hilum_mm = 5.0)
    near_airway = falses(size(lobes))
    r = ceil(Int, hilum_mm)
    for c in findall(airways.airway)
        p = airways.grid.origin .+ (Tuple(c) .- 1) .* airways.grid.h
        q = round.(Int, (p .- grid.origin) ./ grid.h) .+ 1
        for k in -r:r, j in -r:r, i in -r:r
            i^2 + j^2 + k^2 <= hilum_mm^2 || continue
            n = q .+ (i, j, k)
            all(1 .<= n .<= size(lobes)) && (near_airway[n...] = true)
        end
    end
    fissure = falses(size(lobes))
    for c in findall(!=(0), lobes)
        near_airway[c] && continue
        other(s) = checkbounds(Bool, lobes, c + s) && lobes[c + s] != 0 && lobes[c + s] != lobes[c]
        any(other, SIX_STEPS) && (fissure[c] = true)
    end
    return fissure
end

"The labels no new vessel may cross, for both trees: the airways (wall 7, air 11), the heart's chambers and the aorta."
const AIRWAYS_AND_HEART = (0x07, 0x0b, 0x13, 0x14, 0x15, 0x16, 0x1c)

"""
The blockers, at 0.4 mm, one set per tree: `AIRWAYS_AND_HEART`, the other system's imaged lumens (veins 24 and 124 for the
arteries, arteries 25 and 125 for the veins) and the fissures.
"""
BLOCKERS = cached("blockers") do
    fissure = fissures(LOBES, AIRWAYS, growth_grid)
    blocked(labels) = with_fissures(
        V.XCAT.label_blockers(xcat_labels, XCAT_ORIGIN, XCAT_VOXEL_MM, growth_grid, labels),
        fissure,
        growth_grid,
    )
    (
        artery = blocked((AIRWAYS_AND_HEART..., 0x18, 0x7c)),
        vein = blocked((AIRWAYS_AND_HEART..., 0x19, 0x7d)),
        fissure = fissure,
    )
end;

"A tree's blockers `(mask, origin, cell)` with every cell whose centre lies in a fissure voxel (on `grid`) added."
function with_fissures((mask, origin, cell), fissure, grid)
    for c in CartesianIndices(mask)
        p = origin .+ (Tuple(c) .- 1) .* cell
        k = round.(Int, (p .- grid.origin) ./ grid.h) .+ 1
        all(1 .<= k .<= size(fissure)) && fissure[k...] && (mask[c] = 0x01)
    end
    return (mask, origin, cell)
end

"The 1 mm cells no terminal or stub may end in: mostly airway, chamber or aortic blood."
FORBID = V.XCAT.label_mask(xcat_labels, XCAT_ORIGIN, XCAT_VOXEL_MM, growth_grid, AIRWAYS_AND_HEART);

"""
An inlet tree with its stubs (`V.XCAT.penetrating_stubs` at the morphometric spacing; where a vessel runs through the lung,
as almost all of them do, the stub leaves sideways) and its pass-through voxels: its lumens outside the lung, which new
vessels may cross but no terminal may end in.
"""
function inlet(tree, lengths, blockers)
    stubbed, report = V.XCAT.penetrating_stubs(
        tree,
        growth_grid,
        LUNG;
        spacing = lengths.spacing,
        depth = lengths.depth,
        reach = lengths.spacing,
        domain = LUNG,
        forbid = FORBID,
        blockers,
    )
    return (tree = stubbed, report = report, pass = V.XCAT.vessel_tube(stubbed, growth_grid) .& .!FORBID .& .!LUNG)
end

artery_inputs = cached("artery_inputs") do
    inlet(artery_base, STUB_LENGTHS.artery, BLOCKERS.artery)
end;

vein_inputs = cached("vein_inputs") do
    [inlet(tree, STUB_LENGTHS.vein, BLOCKERS.vein) for tree in vein_bases]
end;

"Lobe colours: right lobes warm, left lobes cool."
const LOBE_COLOURS = ["#e41a1c", "#ff7f00", "#ffd92f", "#377eb8", "#984ea3"]

"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, blockers, 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"), LUNG .| 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, territories = allow)
    return V.vtree(engine)
end

arteries = cached("arteries") do
    # the preset is this subject's: its tissue is these lung voxels and its root the trunk's median ring
    @assert V.LUNG_GROWTH.tissue_ml ≈ count(LUNG) / 1000 && V.LUNG_GROWTH.root_caps[1] ≈ artery_base.r[2]
    grow_forest(V.LUNG_GROWTH, [artery_inputs.tree], artery_inputs.pass, BLOCKERS.artery, "arteries")
end;

veins = cached("veins") do
    @assert all(V.LUNG_VEIN_GROWTH.root_caps .≈ [tree.r[1] for tree in vein_bases])
    grow_forest(
        V.LUNG_VEIN_GROWTH,
        [v.tree for v in vein_inputs],
        reduce(.|, [v.pass for v in vein_inputs]),
        BLOCKERS.vein,
        "veins",
    )
end;

"""
The lung's surface: each surface voxel's centre (mm) and its number among the lung voxels, so per-voxel values (such as
perfusion) can be painted on it.
"""
LUNG_SURFACE = let
    function surface_of(lung, origin)
        voxels = findall(lung)
        on_surface(I) = any(s -> !checkbounds(Bool, lung, I + s) || !lung[I + s], SIX_STEPS)
        surface = [n for (n, I) in enumerate(voxels) if on_surface(I)]
        return (points = [Mke.Point3f(Tuple(voxels[n]) .- 1 .+ origin) for n in surface], voxels = surface)
    end
    surface_of(LUNG, growth_grid.origin)
end;

"Draw every fourth point of the lungs' surface as a faint shell around the vessels."
function draw_lung_shell!(ax)
    return Mke.scatter!(ax, LUNG_SURFACE.points[1:4:end]; color = (:gray50, 0.04), markersize = 4, 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

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

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

"Blood flow (mL/min per mL of lung) of every 1 mm lung voxel in circuit `circuit`."
function voxel_perfusion(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
    return vec(V.OrganCT.tissue_iodine(terminal_flow, flows.pairs, flows.terr, voxels)) .* 60
end

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

"""
The injection: the lab's current protocol (Hadjiabdolhamid et al. 2024, Shojazadeh et al. 2026), the same as the heart's:
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.
"""
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 lung circuit `circuit`, in its one (passive) state."
function transport_state(circuit)
    spec = V.TransportSpec(V.LUNG_TRANSPORT; phys = circuit.phys, tissue_ml = Float32(size(tissue_centres, 2) / 1000))
    return ("rest", circuit, V.transport(spec, circuit, INJECTION; times = V.LUNG_PHANTOM.times))
end

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

TRANSPORT = transport_cache("transport", flows.rest);

"The RV and LA blood and the whole lung's iodine (mg/mL) over time, and the cardiac output, from a transport file."
function transport_curves(path)
    tissue_ml = size(tissue_centres, 2) / 1000
    _, circuit, result = only(deserialize(path))
    held = vec(sum(result.tissue; dims = 2) .+ sum(result.artery; dims = 2) .+ sum(result.vein; dims = 2))
    return (
        blood_times = result.body.t,
        rv = result.body.mean[:rv],
        la = result.body.mean[:la],
        lung = held ./ tissue_ml,
        co_l_min = 60 * result.body.circ.co / 1000,
    )
end

TRANSPORT_CURVES = transport_curves(TRANSPORT);

split_lobes = true # the five lobes apart; false: one compartment per lung

begin
    ROI_NAMES = [split_lobes ? LOBAR_BRONCHI.names : ["right lung", "left lung"]; "right hilar vessels"; "left hilar vessels"]
    ROI_GROUPS = let n = length(ROI_NAMES) - 2, right = split_lobes ? [1, 2, 3] : [1], left = split_lobes ? [4, 5] : [2]
        [
            "lungs" => collect(1:n),
            "lungs + hilar vessels" => collect(1:(n + 2)),
            (split_lobes ? ["right lung" => right, "left lung" => left] : Pair{String, Vector{Int}}[])...,
            "right lung + hilar vessels" => [right; n + 1],
            "left lung + hilar vessels" => [left; n + 2],
        ]
    end
    ROI_NAMES, ROI_GROUPS
end

"""
Where the lungs would be with their vessels: XCAT's lung closed by 5 mm on the 1 mm grid, so a vessel lumen inside or
overlapping the lung is in it (the hilar vessels outside the lung are not); and the lobes carried 15 mm out from the lung, so
those lumens get the lobe they lie in.
"""
LUNG_ENVELOPE = cached("lung_envelope") do
    envelope = V.XCAT.closing(V.XCAT.label_mask(xcat_labels, XCAT_ORIGIN, XCAT_VOXEL_MM, growth_grid, (0x06,); fraction = 1.0e-9), 5)
    (inside = envelope, lobe = carried(LOBES, 15))
end;

"`labels` with each nonzero label carried `steps` voxels out into the zeros around it, nearest first."
function carried(labels, steps)
    out = copy(labels)
    front = findall(!=(0), out)
    for _ in 1:steps
        next = CartesianIndex{3}[]
        for c in front, s in SIX_STEPS
            n = c + s
            checkbounds(Bool, out, n) && out[n] == 0 || continue
            out[n] = out[c]
            push!(next, n)
        end
        front = next
    end
    return out
end

make_roi(labels, envelope, grid, classify, hilar; split_lobes) -> roi

"""
    make_roi(labels, envelope, grid, classify, hilar; split_lobes) -> roi

The compartment rule `roi(p)` for the 0.2 mm voxel at `p` (0 = none): lung tissue takes its lobe; a vessel lumen takes the
lobe it lies in when it is inside the lungs' envelope (both on `grid`); outside it, an imaged pulmonary segment's lumen
(`classify`) takes its lung's hilar vessels (`hilar`, 1 right, 2 left, 0 the main pulmonary artery). (Only tissue and
vessel voxels take a compartment.) Everything it reads is captured with its type, so it compiles to fast code.
"""
function make_roi(labels, envelope, grid, classify, hilar; split_lobes)
    inside = envelope.inside
    lobe = envelope.lobe
    origin = grid.origin
    n_lungs = split_lobes ? 5 : 2
    return function roi(p)
        voxel = round.(Int, (p .- XCAT_ORIGIN) ./ XCAT_VOXEL_MM) .+ 1
        all(1 .<= voxel .<= XCAT_DIMS) || return 0
        k = round.(Int, p .- origin) .+ 1
        all(1 .<= k .<= size(lobe)) || return 0
        if labels[voxel...] == 0x06 || inside[k...]
            l = Int(lobe[k...])
            l == 0 && return 0
            return split_lobes ? l : (l <= 3 ? 1 : 2)
        end
        kind, j, _ = classify(p)
        kind == V.VOX_VESSEL || return 0
        side = j > 0 ? hilar.artery[j] : hilar.vein[-j]
        return side == 0 ? 0 : n_lungs + side
    end
end

"""
Each imaged segment's lung (1 right, 2 left; 0 the main pulmonary artery and every grown segment), by the inlet tree's
vessel nearest its middle: the right and left pulmonary artery tubes, and each lung's venous tree.
"""
HILAR_SIDE = let
    function sides(tree, bases)
        side = zeros(Int8, length(tree.r))
        for j in findall(tree.isbase)
            middle = (tree.a[:, j] .+ tree.b[:, j]) ./ 2
            side[j] = argmin(((d, s),) -> d, ((minimum(sum(abs2, X .- middle; dims = 1)), s) for (X, s) in bases))[2]
        end
        return side
    end
    trunk = artery_base.vessel .== "dias_pul_art"
    right = startswith.(artery_base.vessel, "right")
    left = startswith.(artery_base.vessel, "left")
    (
        artery = sides(arteries, [(artery_base.X[:, trunk], 0), (artery_base.X[:, right], 1), (artery_base.X[:, left], 2)]),
        vein = sides(veins, [(vein_bases[1].X, 1), (vein_bases[2].X, 2)]),
    )
end;

roi = make_roi(xcat_labels, LUNG_ENVELOPE, growth_grid, classify, HILAR_SIDE; split_lobes)

"Categorical colours of the compartments (`LOBE_COLOURS`, or one per lung; teal for the hilar vessels, dark on the right)."
roi_colormap() = Mke.cgrad(
    [:black; split_lobes ? LOBE_COLOURS : ["#e41a1c", "#377eb8"]; "#00838f"; "#80deea"],
    length(ROI_NAMES) + 1;
    categorical = true,
)

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

    # the imaged segments of the grown trees, indexed for naming XCAT's drawn pulmonary lumens near them
    const NEAR_MM = 2.0 # a drawn lumen is an imaged segment's when it lies within this of the segment's surface
    # (the arteries' reaches 10 mm: inside the main pulmonary artery's surface every voxel is the trunk's, wherever its
    # grown radius puts the segment's own surface)
    const ARTERY_INDEX = V.XCAT.segment_index(V.XCAT.imaged_segments(arteries); reach = 10.0)
    const VEIN_INDEX = V.XCAT.segment_index(V.XCAT.imaged_segments(veins); reach = NEAR_MM)
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 (`V.XCAT.great_vessel_transit`). Each state's transit is volume over flow (`V.pool_delay`).

XCAT vmale_50's aorta is about three standard deviations narrow (21–22 mm ascending and 15 mm descending, against
33 ± 4 and 24 ± 3 mm in adult men; Wolak et al. 2008): it holds 27 mL from the valve to the arch's top, where the body
model the bolus is carried through holds 60 mL. The aortic tree's volumes are scaled to the body model's
(`V.XCAT.scale_tree`), as for the heart, so the aorta's timing comes from the same body as everything else.
"""
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.LUNG_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 main pulmonary artery's lumen, at 0.4 mm: inside XCAT's `dias_pul_art` surface. XCAT's label volume paints this lumen
as the right ventricle's blood (label 20), so its surface, not its label, tells the trunk from the ventricle; its volume
(mL) is the blood between the valve and the pulmonary arteries the label volume does show.
"""
PULMONARY_TRUNK = cached("pulmonary_trunk") do
    points, triangles = V.XCAT.tessellate(xcat_surfaces["dias_pul_art"]; nu = 160, nv = 160)
    grid = V.grid_around(points, 0.4; margin = 2)
    inside = V.voxelize(points, triangles, grid)
    (grid = grid, inside = inside, ml = count(inside) * 0.4^3 / 1000)
end;

"""
The lungs' phantom spec: the body pools, with the descending aorta timed from this subject's arch and the pulmonary
arteries the label volume shows timed from the valve (the trunk's blood before them).
"""
PHANTOM_SPEC = let pools = copy(V.LUNG_PHANTOM.pools)
    function lead(name, ml)
        i = findfirst(p -> p.name == name, pools)
        d = pools[i]
        return pools[i] = V.Pool(d.name, d.material, d.segment, d.what; through = d.through, lead_ml = ml, lead_through = d.lead_through)
    end
    lead("descending aorta", GREAT_VESSEL_TRANSIT.apex_ml)
    lead("pulmonary arteries", Float32(PULMONARY_TRUNK.ml))
    V.PhantomSpec(V.LUNG_PHANTOM; pools)
end;

make_classifier(labels, grid, numbers, artery_index, vein_index, transit, trunk) -> classify

"""
    make_classifier(labels, grid, numbers, artery_index, vein_index, transit, trunk) -> classify

The classifier `classify(p)`: what is at the 0.2 mm voxel centre `p`. Lung tissue is its 1 mm voxel; a pulmonary vessel
lumen (XCAT's labels 25 and 125 on the pulmonary arterial tree, 24 and 124 on the pulmonary venous tree, and the ventricle's
label inside the main pulmonary artery's surface, `trunk`) is the imaged segment it lies in or beside, else its vessel's pool
and its place along the vessel; every other great-vessel voxel is its
pool; every other XCAT label is decided by the label alone (a body pool, `V.XCAT_POOL_OF_LABEL`, the heart lumped, 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_index, vein_index, transit, trunk)
    # 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 (0x06, 0x14, 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)
    ventricle = (V.VOX_POOL, Int32(V.XCAT_POOL_OF_LABEL[0x14]), 0.0f0)
    # a pulmonary lumen: the imaged segment it lies in or beside, else its pool along the vessel
    function imaged_or_pool(p, index, sign, t)
        j, distance = V.XCAT.nearest_segment(index, p)
        return j > 0 && distance <= NEAR_MM ? (V.VOX_VESSEL, sign * Int32(j), 0.0f0) : vessel_pool(p, t)
    end

    in_trunk(p) = (
        k = round.(Int, (p .- trunk.grid.origin) ./ trunk.grid.h) .+ 1;
        all(1 .<= k .<= size(trunk.inside)) && trunk.inside[k...]
    )

    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]

        # the main pulmonary artery, painted as the ventricle's blood: its imaged trunk segment
        if label == 0x14
            in_trunk(p) || return ventricle
            j, _ = V.XCAT.nearest_segment(artery_index, p)
            return j > 0 ? (V.VOX_VESSEL, Int32(j), 0.0f0) : ventricle
        end
        if label == 0x06
            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 == 0x7d && return imaged_or_pool(p, artery_index, Int32(1), 0x02)
        label == 0x7c && return imaged_or_pool(p, vein_index, Int32(-1), 0x05)
        label == 0x1c && return vessel_pool(p, 0x01)
        if label == 0x19
            t = V.XCAT.vessel_tree_at(transit, p, (0x01, 0x02))
            return t == 0x02 ? imaged_or_pool(p, artery_index, Int32(1), t) : vessel_pool(p, t)
        end
        t = V.XCAT.vessel_tree_at(transit, p, (0x03, 0x04, 0x05, 0x06, 0x07))
        return t == 0x05 ? imaged_or_pool(p, vein_index, Int32(-1), t) : vessel_pool(p, t)
    end
end

classify = make_classifier(
    xcat_labels,
    growth_grid,
    TISSUE_VOXEL_NUMBER,
    ARTERY_INDEX,
    VEIN_INDEX,
    GREAT_VESSEL_TRANSIT,
    PULMONARY_TRUNK,
)

"""
The bolus-tracking points: the main pulmonary artery at the middle of its imaged trunk (the centre of its lumen there), and
the right ventricle at its blood's deepest point (farthest from its wall, on a 1 mm grid).
"""
TRACKING_POINTS = let
    trunk = findall(==("dias_pul_art"), artery_base.vessel)
    artery = Tuple(artery_base.X[:, trunk[(length(trunk) + 1) ÷ 2]])
    ["main pulmonary artery" => artery, "right ventricle" => deepest_point(xcat_labels, 0x14)]
end

"The point (mm) of label `label`'s voxels deepest inside them (farthest from their surface), on a 1 mm grid."
function deepest_point(labels, label)
    step = round(Int, 1 / XCAT_VOXEL_MM)
    coarse = labels[1:step:end, 1:step:end, 1:step:end]
    blood = coarse .== label
    depth = fill(-1, size(blood))
    queue = [c for c in findall(blood) if any(s -> !checkbounds(Bool, blood, c + s) || !blood[c + s], SIX_STEPS)]
    foreach(c -> depth[c] = 0, queue)
    head = 1
    while head <= length(queue)
        c = queue[head]
        head += 1
        for s in SIX_STEPS
            n = c + s
            checkbounds(Bool, blood, n) && blood[n] && depth[n] < 0 || continue
            depth[n] = depth[c] + 1
            push!(queue, n)
        end
    end
    return Tuple(XCAT_ORIGIN .+ (Tuple(argmax(depth)) .- 1) .* (step * XCAT_VOXEL_MM))
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" => "lungs (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 lungs' box of a saved phantom, 80 mm of it about the main pulmonary artery's tracking point (z), the only part read back.
Returns the bundle (labels and compartments) and its first z (mm). Ground truth is each compartment's own, from the bundle's
tables.
"""
function inspect(dir; half_mm = 40)
    point = V.DynamicPhantom.load(dir; volume = false).meta["tracking"]["main pulmonary artery"]["point_mm"]
    lower = collect(growth_grid.origin) .- 5
    upper = collect(growth_grid.origin) .+ collect(growth_grid.dims) .+ 5
    voxel(mm, d) = clamp(round(Int, (mm - XCAT_ORIGIN[d]) / XCAT_VOXEL_MM) + 1, 1, XCAT_DIMS[d])
    bundle = V.DynamicPhantom.load(
        dir;
        rois = true,
        zrange = voxel(point[3] - half_mm, 3):voxel(point[3] + half_mm, 3),
        xrange = voxel(lower[1], 1):voxel(upper[1], 1),
        yrange = voxel(lower[2], 2):voxel(upper[2], 2),
    )
    return (bundle = bundle, z0 = bundle.origin_mm[3])
end

healthy = inspect(HEALTHY_DIR);

"Per-label lookups of a bundle at time `t`: kind (for pictures), idealized 70 keV HU and iodine (mg/mL)."
function label_lookups(bundle, t)
    kind_of = zeros(UInt8, 65536)
    for (label, row) in bundle.table
        label == 0 && continue
        kind_of[label + 1] =
            row.kind == "static" ? 1 :
            row.kind == "pool" ? (occursin("blood", row.base) ? 3 : 2) :
            Dict("tissue" => 4, "artery" => 5, "vein" => 6)[row.kind]
    end
    hu = fill(-1000.0f0, 65536)
    for (label, value) in V.DynamicPhantom.hu(bundle, "rest", t)
        hu[label + 1] = value
    end
    iodine = zeros(Float32, 65536)
    for (label, value) in V.DynamicPhantom.iodine(bundle, "rest", t)
        iodine[label + 1] = value
    end
    return (; kind_of, hu, iodine)
end

# the kinds of label, for pictures
const KIND_COLOURS, KIND_NAMES = (
    [:gray30, :gray65, :mediumpurple, :gray45, :red, :deepskyblue],
    ["static (air, cartilage)", "body tissue", "blood", "lung", "artery", "vein"],
)

"""
One heatmap panel with its own scale underneath (a colour bar, or a legend for categories), at `fig[row, column]`; `xs`
and `ys` (mm) span `image`.
"""
function panel!(fig, row, column, title, xs, ys, image; colormap, colorrange, label = nothing, legend = nothing, yreversed = true)
    cell = fig[row, column] = Mke.GridLayout()
    ax = Mke.Axis(cell[1, 1]; title, aspect = Mke.DataAspect(), yreversed)
    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

"""
One axial phantom slice in six panels: labels by kind, compartments and ground-truth perfusion; the idealized 70 keV image
at time `t` in a mediastinal and a lung window, and iodine at `t`.
"""
function slice_figure(view, 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]
    rois = bundle.rois[is, js, k]
    look = label_lookups(bundle, t)
    # each compartment's ground truth over all of it: its beds' flow over the mass of everything in it
    truth = [0.0f0; [Float32(bundle.roi_table[("rest", name)].mbf) for name in ROI_NAMES]]

    fig = Mke.Figure(; size = (1500, 1000))
    kinds = [look.kind_of[l + 1] for l in labels]
    panel!(
        fig, 1, 1, "Labels by kind", xs, ys, kinds;
        colormap = Mke.cgrad([:black; KIND_COLOURS], 7; categorical = true), colorrange = (-0.5, 6.5), legend = (KIND_COLOURS, KIND_NAMES)
    )
    panel!(
        fig, 1, 2, "Compartments", xs, ys, rois;
        colormap = roi_colormap(), colorrange = (-0.5, length(ROI_NAMES) + 0.5), legend = (roi_colormap().colors[2:end], ROI_NAMES)
    )
    panel!(
        fig, 1, 3, "Ground-truth perfusion by compartment", xs, ys, [truth[r + 1] for r in rois];
        colormap = :turbo, colorrange = (0, 12), label = "mL/min/g"
    )
    hu = [look.hu[l + 1] for l in labels]
    panel!(fig, 2, 1, "Idealized 70 keV, t = $t s", xs, ys, hu; colormap = :grays, colorrange = (-200, 600), label = "HU")
    panel!(fig, 2, 2, "Idealized 70 keV, lung window, t = $t s", xs, ys, hu; colormap = :grays, colorrange = (-1000, 200), label = "HU")
    panel!(
        fig, 2, 3, "Iodine, t = $t s", xs, ys, [look.iodine[l + 1] for l in labels];
        colormap = :turbo, colorrange = (0, 8), label = "mg/mL"
    )
    Mke.Label(fig[0, :], "z = $(round(z; digits = 1)) mm"; fontsize = 20, font = :bold)
    for column in 1:3
        Mke.colsize!(fig.layout, column, Mke.Relative(1 / 3))
    end
    return fig
end

slice_figure(healthy, healthy_time, healthy_z)

"""
One coronal plane of a saved phantom, through the main pulmonary artery's tracking point, at time `t`: compartments,
ground-truth perfusion by compartment, the idealized 70 keV image in a lung window, and iodine. Reads the whole phantom
once for this one plane.
"""
function coronal_figure(dir, t)
    point = V.DynamicPhantom.load(dir; volume = false).meta["tracking"]["main pulmonary artery"]["point_mm"]
    j = clamp(round(Int, (point[2] - XCAT_ORIGIN[2]) / XCAT_VOXEL_MM) + 1, 1, XCAT_DIMS[2])
    voxel(mm, d) = clamp(round(Int, (mm - XCAT_ORIGIN[d]) / XCAT_VOXEL_MM) + 1, 1, XCAT_DIMS[d])
    bundle = V.DynamicPhantom.load(
        dir;
        rois = true,
        yrange = j:j,
        xrange = voxel(growth_grid.origin[1] - 5, 1):voxel(growth_grid.origin[1] + growth_grid.dims[1] + 5, 1),
        zrange = voxel(growth_grid.origin[3] - 5, 3):voxel(growth_grid.origin[3] + growth_grid.dims[3] + 5, 3),
    )
    xs = bundle.origin_mm[1] .+ (axes(bundle.labels, 1) .- 1) .* XCAT_VOXEL_MM
    zs = bundle.origin_mm[3] .+ (axes(bundle.labels, 3) .- 1) .* XCAT_VOXEL_MM
    labels = bundle.labels[:, 1, :]
    rois = bundle.rois[:, 1, :]
    look = label_lookups(bundle, t)
    truth = [0.0f0; [Float32(bundle.roi_table[("rest", name)].mbf) for name in ROI_NAMES]]

    fig = Mke.Figure(; size = (1600, 560))
    panel!(
        fig, 1, 1, "Compartments", xs, zs, rois;
        colormap = roi_colormap(), colorrange = (-0.5, length(ROI_NAMES) + 0.5), legend = (roi_colormap().colors[2:end], ROI_NAMES), yreversed = false
    )
    panel!(
        fig, 1, 2, "Ground-truth perfusion by compartment", xs, zs, [truth[r + 1] for r in rois];
        colormap = :turbo, colorrange = (0, 12), label = "mL/min/g", yreversed = false
    )
    panel!(
        fig, 1, 3, "Idealized 70 keV, lung window, t = $t s", xs, zs, [look.hu[l + 1] for l in labels];
        colormap = :grays, colorrange = (-1000, 200), label = "HU", yreversed = false
    )
    panel!(
        fig, 1, 4, "Iodine, t = $t s", xs, zs, [look.iodine[l + 1] for l in labels];
        colormap = :turbo, colorrange = (0, 8), label = "mg/mL", yreversed = false
    )
    Mke.Label(fig[0, :], "coronal y = $(round(point[2]; digits = 1)) mm"; fontsize = 20, font = :bold)
    for column in 1:4
        Mke.colsize!(fig.layout, column, Mke.Relative(1 / 4))
    end
    return fig
end

coronal_figure(HEALTHY_DIR, 14)

"Each compartment's iodine, with the tracking arrival and peak, and the lungs' iodine ledger."
function curves_figure(view)
    bundle = view.bundle
    fig = Mke.Figure(size = (1400, 520))
    ax = Mke.Axis(fig[1, 1]; title = "Compartments", xlabel = "time after injection (s)", ylabel = "iodine (mg)")
    for (k, name) in enumerate(ROI_NAMES)
        row = bundle.roi_table[("rest", name)]
        Mke.lines!(
            ax,
            bundle.times,
            row.mass;
            linewidth = 2,
            color = roi_colormap().colors[k + 1],
            label = @sprintf("%s · %.2f mL/min/g", name, row.mbf),
        )
    end
    timing = V.DynamicPhantom.scan_times(bundle; state = "rest", roi = "roi:main pulmonary artery")
    Mke.vlines!(ax, [timing.arrival, timing.peak_time]; color = :gray40, linestyle = :dash)
    Mke.xlims!(ax, 0, 40)
    Mke.axislegend(ax; position = :rt, labelsize = 11)

    ledger = bundle.ledger["rest"]
    ax = Mke.Axis(fig[1, 2]; title = "Iodine ledger", 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 pulmonary 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 lungs through the pulmonary veins, 1, 2 and 3 s after arrival
and at the arterial peak (bolus tracking in the main pulmonary artery).
"""
function outflow_check(view)
    bundle = view.bundle
    ledger = bundle.ledger["rest"]
    function share_out(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
    timing = V.DynamicPhantom.scan_times(bundle; state = "rest", roi = "roi:main pulmonary artery")
    shares = [@sprintf("%.0f%%", 100 * share_out(t)) for t in [timing.arrival .+ [1, 2, 3]; timing.peak_time]]
    text = """
    !!! note "Iodine already out of the lungs"
        | | arrival + 1 s | arrival + 2 s | arrival + 3 s | arterial peak |
        |---|---|---|---|---|
        | arrival at $(round(timing.arrival; digits = 1)) s, peak at $(round(timing.peak_time; digits = 1)) s | $(join(shares, " | ")) |

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

curves_figure(healthy)

outflow_check(healthy)

"""
The occluded artery: of the imaged arterial segments whose terminals lie ≥ 95% in the left lower lobe and that feed at most
a third of it, the one feeding the most; its lumen closed to `area` of its own over its length (at most 10 mm), as `circuit`
takes a narrowing.
"""
OCCLUSION = occluded_artery(arteries, LOBES, growth_grid; lobe = 5, area = 1.0e-3)

"""
The occluded artery of `tree`: the imaged segment that feeds the most of lobe `lobe` (`lobes` on `grid`) among those whose
terminals lie ≥ 95% in it and that feed at most a third of it; with `circuit`'s narrowing to `area` of its lumen.
"""
function occluded_artery(tree, lobes, grid; lobe, area)
    n = length(tree.r)
    lobe_at(p) = (k = round.(Int, (p .- grid.origin) ./ grid.h) .+ 1; all(1 .<= k .<= grid.dims) ? Int(lobes[k...]) : 0)
    children = [Int[] for _ in 1:n]
    foreach(j -> tree.parent[j] > 0 && push!(children[tree.parent[j]], j), 1:n)
    depth = [
        (
            k = j; d = 0; while tree.parent[k] != 0
                k = tree.parent[k]; d += 1
            end; d
        ) for j in 1:n
    ]
    counts = zeros(Int, 5, n)
    for j in 1:n
        isempty(children[j]) && (l = lobe_at(tree.b[:, j]); l > 0 && (counts[l, j] += 1))
    end
    for j in sortperm(depth; rev = true)
        p = tree.parent[j]
        p > 0 && (counts[:, p] .+= counts[:, j])
    end
    total = sum(counts[lobe, j] for j in 1:n if isempty(children[j]))
    candidates = [
        j for j in 1:n if tree.isbase[j] && sum(counts[:, j]) > 0 &&
            counts[lobe, j] >= 0.95 * sum(counts[:, j]) && counts[lobe, j] <= total / 3
    ]
    j = candidates[argmax([counts[lobe, c] for c in candidates])]
    return (segment = j, share = counts[lobe, j] / total, stenosis = [(j, area, 10.0)])
end

"""
The occluded lungs' circuit. The right ventricle keeps the cardiac output: the inlet pressure is raised (by the secant rule)
until the lungs carry the healthy flow again.
"""
flows_occluded = cached("flow_occluded") do
    target = V.organ_flow(flows.rest)
    solve(p_in) = V.Flow.circuit(
        arteries,
        flows.terr;
        state = :rest,
        veins,
        vleaf = flows.vleaf,
        reference = flows.rest,
        phys = V.FlowSpec(flows.rest.phys; p_in = Float32(p_in)),
        stenosis = OCCLUSION.stenosis,
        to = DEVICE,
    )
    p0 = Float64(flows.rest.phys.p_in)
    c0 = solve(p0)
    p1 = p0 + 1
    c1 = solve(p1)
    for _ in 1:12
        f0 = V.organ_flow(c0) - target
        f1 = V.organ_flow(c1) - target
        abs(f1) <= 1.0e-4 * target && break
        p0, c0, p1 = p1, c1, p1 - f1 * (p1 - p0) / (f1 - f0)
        c1 = solve(p1)
    end
    (rest = c1, stress = c1, p_in = p1)
end;

TRANSPORT_OCCLUDED = transport_cache("transport_occluded", flows_occluded.rest);

OCCLUDED_DIR = save_phantom("lll_occlusion", TRANSPORT_OCCLUDED);

occluded = inspect(OCCLUDED_DIR);

slice_figure(occluded, occluded_time, occluded_z)

coronal_figure(OCCLUDED_DIR, 14)

curves_figure(occluded)

outflow_check(occluded)