12 · Siemens SOMATOM Definition Flash · Dual-Source Regular + VMI

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

using Markdown: @md_str, Markdown

using Statistics: mean, std, quantile

import BasisSimulator as BS

# Use CairoMakie for faithful build-time rendering. Snapshot can still isolate
# and compile independent browser-safe islands without hoisting this import.
import CairoMakie as Mke

import PlutoUI

begin
    import GPUSelect
    AT = GPUSelect.Storage()     # the backend array type, directly: MtlArray / CuArray / ROCArray
    to_gpu(x) = AT(x)
    GPU_BACKEND = (name = string(nameof(AT)),)
end;

phantom_cpu = BS.create_gammex_472(
    n_voxels = 512,
    n_slices = 16,
    fov_cm = 35.0,
    z_cm = 1.0,
);

phantom = BS.Phantom(
    to_gpu(phantom_cpu.mask),
    phantom_cpu.materials,
    phantom_cpu.voxel_size,
    phantom_cpu.origin,
    phantom_cpu.extent,
);

# Siemens SOMATOM Definition Flash: basis-spectral-denoising's `flash_scanner`, verbatim. Both
# tubes run on detector A's arc so the two channels stay co-registered ray by ray (see md above).
flash_scanner = BS.EICTScanner(
    source_to_isocenter = 595.0,
    source_to_detector = 1085.6,
    detector_rows = 64,
    detector_cols = 736,
    detector_row_size = 0.6,
    detector_col_size = 0.70473,
    detector_shape = :arc,
    focal_spot_width = 0.7,
    focal_spot_length = 0.7,
    target_angle = 7.0,
    gantry_rotation_time = 0.5,
    scan_diameter = 500.0,
    gantry_aperture = 780.0,
    flat_filter_material = :aluminum,
    flat_filter_thickness = 8.4,          # 6.8 mm Al tube assembly + 1.6 mm beam limiter
    bowtie_filter = :large_body,
    detector_material = :ufc_flash,       # the Flash's own Monte Carlo table (NOT :ufc)
    detector_depth = 1.0,
    fill_factor_row = 0.9,
    fill_factor_col = 0.9,
    electronic_noise = 1500,              # e⁻ — Stellar DAS floor; pre-2011 build ≈ 3500 (§3)
    detection_gain = 10.0,
);

begin
    VIEW_SAMPLES = 5
    flash_opts = BS.SimOptions(projector = :dd_fast, use_heel_effect = false, view_samples = VIEW_SAMPLES)
    "The same physics switches with one field replaced."
    with_option(opts, field, value) = BS.SimOptions(;
        (f => getfield(opts, f) for f in fieldnames(BS.SimOptions) if f !== field)...,
        (field => value,)...,
    )
end;

# The dual-energy acquisition: basis-spectral-denoising's `SCANNERS.flash`, verbatim. Seeds are per
# exposure: a shared seed would correlate the two channels' noise.
FLASH = (name = "Siemens SOMATOM Flash", mechanism = "dual source", kind = :eict,
    scanner = flash_scanner, opts = flash_opts, views = 1152, rotation_time = 0.5,
    collimation_mm = 4.8, slice_mm = 1.0, n_slices = 4, fov_cm = 33.0, anode_angle = 8,
    clinical_recon = "DE monoenergetic, D40f, 1.0 mm, 330 mm", ctdi_vol_mGy = 10.01,
    exposures = [
        (label = "80 kVp", kVp = 80, mA = 470.0, view_arc = 1.0, seed = 1234, filters = Tuple{String, Float64}[]),
        (label = "Sn140 kVp", kVp = 140, mA = 182.0, view_arc = 1.0, seed = 4321, filters = [("Sn", 0.4)]),
    ]);

begin
    "The CTProtocol of one exposure of scanner `s`, its tube current scaled by `mA_scale` (`DOSE_SCALE`)."
    exposure_protocol(s, e; mA_scale = 1.0) = BS.CTProtocol(;
        kVp = e.kVp, mA = e.mA * mA_scale, views = s.views, rotation_time = s.rotation_time,
        collimation_mm = s.collimation_mm, additional_filters = e.filters,
        (s.anode_angle === nothing ? (;) : (; anode_angle = s.anode_angle))...,
    )
    "Reconstruction grid of scanner `s`: its clinical series' slices, field of view and 512²."
    recon_options(s) = BS.ReconOptions(
        matrix_size = (512, 512, s.n_slices),
        fov_cm = s.fov_cm, z_cm = s.n_slices * s.slice_mm / 10,
    )
end;

# One factor for both tubes: the simulated CTDIvol (32 cm body phantom, summed over the exposures)
# equals the physical scan's.
DOSE_SCALE = FLASH.ctdi_vol_mGy / sum(FLASH.exposures) do e
    p = exposure_protocol(FLASH, e)
    BS.compute_dose(BS.dose_source(FLASH.scanner, p), p).ctdi_vol_mGy
end

begin
    # dual energy: tube A (low) and tube B (high), at the physical scan's dose
    protocol_low, protocol_high = (exposure_protocol(FLASH, e; mA_scale = DOSE_SCALE) for e in FLASH.exposures)
    # regular (dual power): 120 kVp, 420 mA on each tube, on the same simulated tube
    protocol_reg = BS.CTProtocol(
        kVp = 120, mA = 420.0 * DOSE_SCALE, views = FLASH.views, rotation_time = FLASH.rotation_time,
        collimation_mm = FLASH.collimation_mm, anode_angle = FLASH.anode_angle,
    )
end;

# The clinical series' grid (`recon_options`): 4 × 1.0 mm slices over 33 cm, 512².
recon_opts = recon_options(FLASH);

# The Gammex's labels on the reconstruction grid (the phantom itself is voxelized over 35 cm):
# every ROI below is drawn on this map.
roi_labels = BS.create_gammex_472(
    n_voxels = recon_opts.matrix_size[1], n_slices = 1, fov_cm = recon_opts.fov_cm, z_cm = 0.1,
).mask[:, :, 1];

"""
    flash_detected_spectrum(protocol) -> (e, w_eta)

Display helper: tube spectrum × flat filtration × protocol filters
(IPEM, absolute flux) with the src Flash UFC MC η(E) folded in — what the
detector actually integrates (centered ray, no bowtie).
"""
function flash_detected_spectrum(protocol)
    e, w = BS.resolve_source_spectrum_without_bowtie(
        flash_opts, protocol; scanner = flash_scanner,
    )
    return e, Float64.(w) .* BS.get_ufc_flash_mc_efficiency.(e)
end;

"""
    flash_bhc_calibration(protocol, geom)

η-aware per-tube BHC: bowtie-hardened per-column spectrum × Flash UFC η(E)
→ per-column two-material polynomial.  Returns `(model, μ_water, ref_E_keV)`.
"""
function flash_bhc_calibration(protocol, geom)
    # resolve_source_spectrum_full folds bowtie AND the src Flash UFC η(E)
    # (via the same build_physics_config the forward model used).
    e, ŵ = BS.resolve_source_spectrum_full(
        flash_opts, protocol; scanner = flash_scanner, geom = geom,
    )
    e2, w_col = BS.bhc_spectrum_per_column(e, ŵ)          # [n_E, n_col]

    # Single mono-equivalent target = mean energy of the η-folded mean spectrum
    w_mean = vec(sum(w_col; dims = 2)) ./ size(w_col, 2)
    ref_E = sum(e2 .* w_mean) / sum(w_mean)

    # Parameter-free water BHC from the per-column Flash-η spectrum.
    model = BS.calibrate_bhc_water(
        e2, w_col;
        reference_energy_keV = ref_E,
    )
    return (model = model, μ_water = model.μ_water_ref, ref_E_keV = model.reference_energy_keV)
end;

"""
    flash_poly_recon(sino_cpu, geom, bhc) -> Array{Float32, 3}

Doctrine correction stack for one tube: knobless water sino-BHC → FDK → HU.
"""
function flash_poly_recon(sino_cpu, geom, bhc)
    matrix_size = recon_opts.matrix_size

    sino_gpu = to_gpu(sino_cpu)
    # Knobless water BHC: one sinogram-domain pass, no recon round-trip.
    sino_bhc = BS.apply_bhc_water(sino_gpu, bhc.model)
    sino_gpu = sino_bhc

    ws_fdk = BS.create_fdk_recon_workspace(sino_gpu, geom, matrix_size)
    recon_μ = BS.reconstruct!(ws_fdk, sino_gpu, geom)

    hu = Float32.(BS.to_hounsfield(Array(recon_μ); μ_water = bhc.μ_water))

    ws_fdk = nothing; sino_gpu = nothing; recon_μ = nothing
    GC.gc(true)
    return hu
end;

sim_reg = let projection = nothing, out = Dict{Symbol, Any}()
    for (tube, seed) in ((:a, 1234), (:b, 4321))
        opts = with_option(flash_opts, :seed, seed)
        ws = BS.create_workspace(flash_scanner, protocol_reg, opts, recon_opts, phantom)
        try
            t = time()
            result = BS.simulate!(ws, phantom, protocol_reg, opts; report_dose = false,
                keep_projection = projection === nothing, projection)
            projection === nothing && (projection = result.projection)   # tube B reuses tube A's
            out[tube] = (sino = Array(ws.sinogram), geom = ws.geom, seconds = time() - t)
        finally
            BS.release_backend!(ws)
        end
    end
    projection = nothing; GC.gc(true)
    (a = out[:a], b = out[:b])
end;

bhc_reg = flash_bhc_calibration(protocol_reg, sim_reg.a.geom);

hu_reg = (
    a = flash_poly_recon(sim_reg.a.sino, sim_reg.a.geom, bhc_reg),
    b = flash_poly_recon(sim_reg.b.sino, sim_reg.b.geom, bhc_reg),
);

# Dual power: the two independent same-kV acquisitions average into the
# combined image (image domain — per-ray sinogram averaging is equivalent
# to first order on the log data at these noise levels).
hu_reg_combined = 0.5f0 .* hu_reg.a .+ 0.5f0 .* hu_reg.b;

dp_stats = let
    ERODE_PX = 12.0
    sw_bool = BS.erode_mask_2d(
        roi_labels .== UInt8(BS.REGION_SOLID_WATER); erode_px = ERODE_PX,
    )
    sw_idx = findall(sw_bool)

    # Central noise ROI (12 px ≈ 7.7 mm at 0.645 mm/px) — background water.
    nx_r, ny_r, _ = size(hu_reg_combined)
    cx = nx_r ÷ 2 + 1; cy = ny_r ÷ 2 + 1
    noise_bool = falses(nx_r, ny_r)
    r² = 12.0^2
    @inbounds for j in 1:ny_r, i in 1:nx_r
        ((i - cx)^2 + (j - cy)^2) ≤ r² && (noise_bool[i, j] = true)
    end
    noise_idx = findall(noise_bool)
    n_z = size(hu_reg_combined, 3)

    function _roi(vol, idx)
        vals = Float64[Float64(vol[ci, z]) for z in 1:n_z, ci in idx]
        (mean = mean(vals), std = std(vals), n = length(vals))
    end

    water = (
        a = _roi(hu_reg.a, sw_idx),
        b = _roi(hu_reg.b, sw_idx),
        combined = _roi(hu_reg_combined, sw_idx),
    )
    noise = (
        a = _roi(hu_reg.a, noise_idx),
        b = _roi(hu_reg.b, noise_idx),
        combined = _roi(hu_reg_combined, noise_idx),
    )
    σ_single = 0.5 * (noise.a.std + noise.b.std)
    ratio = noise.combined.std / σ_single

    for (tag, s) in pairs(water)
        @info "[regular 120/120 · Flash] $(tag) SW ROI: ⟨HU⟩ = $(round(s.mean, digits = 2)), σ = $(round(s.std, digits = 2)) HU"
    end
    @info "[regular 120/120 · Flash] noise σ: A = $(round(noise.a.std, digits = 2)), " *
        "B = $(round(noise.b.std, digits = 2)), combined = $(round(noise.combined.std, digits = 2)) HU " *
        "→ ratio $(round(ratio, digits = 3)) (ideal 1/√2 = 0.707)"

    (
        water = water, noise = noise,
        noise_ratio = ratio,
        sw_mask_2d = collect(sw_bool), noise_mask_2d = noise_bool,
    )
end;

begin
    # Draw 0 is the noise-free expectation; draw r ≥ 1 gives each exposure the seed
    # `seed + 1000 (r − 1)`, so the two tubes stay independent within a draw, and draws of each other.
    DE_DRAWS = (reference = 0, measured = 1, calibration = 9)
    realization_seed(e, r) = e.seed + 1000 * (max(r, 1) - 1)
    opts_for(e, r) = with_option(with_option(with_option(FLASH.opts, :seed, realization_seed(e, r)),
        :use_noise, r > 0), :view_arc, e.view_arc)
end;

"""
    acquire(draws) -> Dict{Int, NamedTuple}

The dual-energy acquisition of the Gammex (both tubes of `FLASH`) for every draw in `draws`: draw 0
the noise-free expectation, draw `r ≥ 1` a noisy realisation with the seeds `realization_seed(e, r)`.
Each exposure is projected once, by the first draw, and every other draw re-uses that projection
(`keep_projection`, `projection`), only its noise drawn. Each entry holds the two log-transmission
channels, the spectral basis the decomposition inverts (shared by every draw), the geometry and the
wall time of the draw per tube.
"""
function acquire(draws)
    parts = map(FLASH.exposures) do e
        protocol = exposure_protocol(FLASH, e; mA_scale = DOSE_SCALE)
        projection = nothing
        channels = Dict{Int, Array{Float32, 3}}()
        seconds = Dict{Int, Float64}()
        model = nothing
        for r in draws
            opts = opts_for(e, r)
            t = time()
            ws = BS.create_workspace(FLASH.scanner, protocol, opts, recon_opts, phantom)
            try
                result = BS.simulate!(ws, phantom, protocol, opts; report_dose = false,
                    keep_projection = projection === nothing, projection)
                projection === nothing && (projection = result.projection)
                channels[r] = Array(ws.sinogram)
                if model === nothing
                    air = ws.bowtie_air_reference === nothing ?
                        ones(Float32, ws.geom.n_cols, ws.geom.n_rows) :
                        Float32.(Array(ws.bowtie_air_reference))
                    I0 = BS.compute_detector_I0(ws.geom, protocol, sum(ws.weights)) * Float64(ws.η_eff)
                    energies, response = BS.resolve_source_spectrum_full(
                        opts, protocol; scanner = FLASH.scanner, geom = ws.geom)
                    model = (geom = ws.geom, I0_ray = Float32.(I0 .* air),
                        energies = Float64.(energies), response = Float32.(response))
                end
            finally
                BS.release_backend!(ws)
            end
            seconds[r] = time() - t
        end
        projection = nothing
        (; channels, seconds, model...)
    end
    GC.gc(true)
    basis = BS.spectral_basis_from_acquisitions(acquisitions = [
        (energies = p.energies, response = p.response, I0_ray = p.I0_ray) for p in parts])
    Dict(r => (draw = r, channels = [p.channels[r] for p in parts], basis = basis,
               geom = first(parts).geom, seconds = [p.seconds[r] for p in parts])
         for r in draws)
end;

de_acq = acquire((DE_DRAWS.reference, DE_DRAWS.measured, DE_DRAWS.calibration));

# The measured draw: every DE result below comes from it.
de_measured = de_acq[DE_DRAWS.measured];

bhc_low = flash_bhc_calibration(protocol_low, de_measured.geom);

bhc_high = flash_bhc_calibration(protocol_high, de_measured.geom);

hu_tube = (
    low = flash_poly_recon(de_measured.channels[1], de_measured.geom, bhc_low),
    high = flash_poly_recon(de_measured.channels[2], de_measured.geom, bhc_high),
);

# Siemens linear mixed image: M = w·I_low + (1−w)·I_high (image domain, Yu 2009): the equal blend.
MIX_W_LOW = 0.5f0;

hu_mixed = MIX_W_LOW .* hu_tube.low .+ (1.0f0 - MIX_W_LOW) .* hu_tube.high;

poly_water_stats = let
    ERODE_PX = 12.0
    sw_bool = BS.erode_mask_2d(
        roi_labels .== UInt8(BS.REGION_SOLID_WATER); erode_px = ERODE_PX,
    )
    sw_idx = findall(sw_bool)
    n_z = size(hu_mixed, 3)

    function _stats(vol)
        vals = Float64[Float64(vol[ci, z]) for z in 1:n_z, ci in sw_idx]
        (mean = mean(vals), std = std(vals), n = length(vals))
    end
    stats = (
        low = _stats(hu_tube.low),
        high = _stats(hu_tube.high),
        mixed = _stats(hu_mixed),
    )
    for (tag, s) in pairs(stats)
        @info "[poly · Flash] $(tag) SW ROI: ⟨HU⟩ = $(round(s.mean, digits = 2)), σ = $(round(s.std, digits = 2)) HU (n = $(s.n))"
    end
    (stats..., mask_2d = collect(sw_bool))
end;

basis = de_measured.basis;

# basis-spectral-denoising's HYPR_CHAIN: ProjectionHYPR(3 × 3 Box, view_stride = 2) and ImageHYPR()
HYPR_CHAIN = BS.SpectralHYPR();

VMI_CHAIN = (method = :nchannel, controls = BS.NChannelControls(), use_tlbf = false, antialias = true);

# The Flash's fitted pair of FDK windows (basis-spectral-denoising's `RECON.flash`), fitted on
# 512² over 33 cm and rescaled onto this grid's frequency axis as its `RECON_BY` does.
PAIR_FILTER = let knots = Tuple(range(0.0, 1.0, length = 11)),
        fitted = (composite = (fc = 1.05, p = 1.5), complement = (fc = 0.8, p = 1.0)),
        g = de_measured.geom,
        r = BS.grid_bandlimit(g, recon_opts.matrix_size) / min(1.0, g.pixel_size / (33.0 / 512))
    window(q) = BS.CustomFilter(knots, Tuple(round(exp(-(x * r / q.fc)^q.p), digits = 5) for x in knots))
    (filter = BS.PairFilter(window(fitted.composite), window(fitted.complement)), kernels = fitted, scale = r)
end;

# `pair_basis` and `composite_energy` (basis-spectral-denoising's `PAIR_BASIS`), from the calibration draw
calibration = let a = de_acq[DE_DRAWS.calibration]
    decompose(denoiser) = BS.vmi_pipeline(; channels = a.channels, basis = a.basis, geom = a.geom,
        to_backend = to_gpu, matrix_size = recon_opts.matrix_size, vmi_energies = Tuple(de_vmi_energies),
        denoiser, keep_sinograms = true, VMI_CHAIN..., use_acnr = false).sinograms
    sp(d; kw...) = BS.spectral_pair(d.water, d.iodine, a.geom, recon_opts.matrix_size;
        filter = BS.SoftFilter(), antialias = VMI_CHAIN.antialias, to_backend = to_gpu, kw...)
    x = sp(decompose(nothing))                                   # the plain decomposition
    c = sp(decompose(BS.SpectralHYPR(image = nothing));          # after the projection HYPR
        basis = (Estar = x.Estar, β = x.β))
    (pair_basis = (Estar = x.Estar, β = x.β), composite_energy = c.Estar)
end;

de_vmi_energies = [50.0, 70.0, 100.0, 140.0];

de_vmi = BS.vmi_pipeline(;
    channels = de_measured.channels,
    basis,
    geom = de_measured.geom,
    to_backend = to_gpu,
    matrix_size = recon_opts.matrix_size,
    vmi_energies = Tuple(de_vmi_energies),
    denoiser = HYPR_CHAIN,
    fbp_filter = PAIR_FILTER.filter,
    pair_basis = calibration.pair_basis,
    composite_energy = calibration.composite_energy,
    keep_sinograms = true,
    VMI_CHAIN...,
);

# The same draw and reconstruction without the denoising (no HYPR, no ACNR), and the noise-free
# reference (draw 0, the same projection), reconstructed alike.
de_none, de_ref = (BS.vmi_pipeline(;
        channels = de_acq[r].channels, basis, geom = de_measured.geom, to_backend = to_gpu,
        matrix_size = recon_opts.matrix_size, vmi_energies = Tuple(de_vmi_energies),
        denoiser = nothing, use_acnr = false, fbp_filter = PAIR_FILTER.filter,
        pair_basis = calibration.pair_basis, composite_energy = calibration.composite_energy,
        VMI_CHAIN...,
    ) for r in (DE_DRAWS.measured, DE_DRAWS.reference));

solid_water_basis = let
    ERODE_PX = 12.0

    sw_bool_raw = (roi_labels .== UInt8(BS.REGION_SOLID_WATER))
    sw_bool = BS.erode_mask_2d(sw_bool_raw; erode_px = ERODE_PX)

    n_raw = count(sw_bool_raw); n_eroded = count(sw_bool)
    n_eroded == 0 && error(
        "solid_water_basis: deep erosion (σ = $(ERODE_PX) px) wiped out the SW " *
            "ROI (raw count = $(n_raw)).  Reduce erode_px or check phantom mask."
    )

    sw_idx = findall(sw_bool)
    _mean(vol) = mean(vol[ci, z] for z in axes(vol, 3) for ci in sw_idx)

    c_w = Float64(_mean(de_vmi.images.water))
    c_i = Float64(_mean(de_vmi.images.iodine))
    @info "solid_water_basis: ⟨c_water⟩_SW = $(round(c_w, digits = 4)) g/cm³, " *
        "⟨c_iodine⟩_SW = $(round(1000c_i, digits = 3)) mg/mL"

    (
        c_water = c_w, c_iodine = c_i, n_voxels = length(sw_idx) * size(de_vmi.images.water, 3),
        mask_2d = collect(sw_bool),
    )
end;

# The VMI stack as one volume per energy (keV → (nx, ny, nz) HU).
vmi_HU_final = Dict(
    Float64(E) => de_vmi.vmis[:, :, :, k] for (k, E) in pairs(de_vmi.energies)
);

ROD_LABELS = (
    Ca = (UInt8(10), UInt8(11), UInt8(12), UInt8(13), UInt8(14), UInt8(15), UInt8(16)),
    I = (UInt8(20), UInt8(21), UInt8(22), UInt8(23), UInt8(24), UInt8(25), UInt8(26)),
);

ROD_NAMES = (
    Ca = ("50 mg/mL", "100 mg/mL", "200 mg/mL", "300 mg/mL", "400 mg/mL", "500 mg/mL", "600 mg/mL"),
    I = ("2.0 mg/mL", "2.5 mg/mL", "5.0 mg/mL", "7.5 mg/mL", "10.0 mg/mL", "15.0 mg/mL", "20.0 mg/mL"),
);

rod_data = let
    materials = phantom_cpu.materials
    mask_2d = roi_labels
    nx, ny = size(mask_2d)
    ROI_RADIUS_PX = 8

    function rod_centroid(label::UInt8)
        idx = findall(==(label), mask_2d)
        isempty(idx) && error("rod_centroid: no voxels with label $label")
        cx = sum(ci -> Float64(ci[1]), idx) / length(idx)
        cy = sum(ci -> Float64(ci[2]), idx) / length(idx)
        return (cx, cy)
    end

    function rod_roi_mask(label::UInt8)
        cx, cy = rod_centroid(label)
        i_lo = max(1, floor(Int, cx - ROI_RADIUS_PX))
        i_hi = min(nx, ceil(Int, cx + ROI_RADIUS_PX))
        j_lo = max(1, floor(Int, cy - ROI_RADIUS_PX))
        j_hi = min(ny, ceil(Int, cy + ROI_RADIUS_PX))
        roi = CartesianIndex{2}[]
        r² = Float64(ROI_RADIUS_PX)^2
        for j in j_lo:j_hi, i in i_lo:i_hi
            ((i - cx)^2 + (j - cy)^2) ≤ r² && push!(roi, CartesianIndex(i, j))
        end
        return roi
    end

    rod_rois = Dict(
        lab => rod_roi_mask(lab)
            for lab in vcat(collect(ROD_LABELS.Ca), collect(ROD_LABELS.I))
    )

    μ_water_E = Dict(
        E => BS.compute_μ_at_energy(BS.XA.Materials.water, E)
            for E in de_vmi_energies
    )

    function theoretical_hu(material, E::Float64)
        μ = BS.compute_μ_at_energy(material, E)
        return 1000.0 * (μ - μ_water_E[E]) / μ_water_E[E]
    end

    function measured_hu(vmi_vol, label::UInt8)
        roi = rod_rois[label]
        s = 0.0; n = 0
        for z in 1:size(vmi_vol, 3), ci in roi
            s += vmi_vol[ci, z]; n += 1
        end
        return s / n
    end

    out = Dict{Symbol, NamedTuple}()
    for group in (:Ca, :I)
        labels = ROD_LABELS[group]
        n_rods = length(labels)
        n_E = length(de_vmi_energies)
        meas = zeros(Float64, n_rods, n_E)
        theo = zeros(Float64, n_rods, n_E)
        for (i, lab) in pairs(labels)
            mat = materials[Int(lab) + 1]   # mask_value + 1
            for (j, E) in pairs(de_vmi_energies)
                meas[i, j] = measured_hu(vmi_HU_final[E], lab)
                theo[i, j] = theoretical_hu(mat, E)
            end
        end
        out[group] = (
            labels = labels, names = ROD_NAMES[group],
            measured = meas, theoretical = theo,
        )
    end
    out
end;

vmi_noise_by_keV = let
    roi_idx = findall(solid_water_basis.mask_2d)
    nz_r = size(de_vmi.vmis, 3)
    vals(run, k) = Float64[Float64(run.vmis[ci, z, k]) for z in 1:nz_r, ci in roi_idx]

    out = Dict{Float64, NamedTuple}()
    for (k, E) in enumerate(de_vmi_energies)
        p, n, r = vals(de_vmi, k), vals(de_none, k), vals(de_ref, k)
        out[E] = (mean = mean(p), reference = mean(r), std = std(p .- r), std_none = std(n .- r),
            n = length(p))
        @info "water region @ $(Int(E)) keV: ⟨HU⟩ = $(round(mean(p), digits = 2)) " *
            "(reference $(round(mean(r), digits = 2))), σ = $(round(std(p .- r), digits = 2)) HU " *
            "(no denoising $(round(std(n .- r), digits = 2)) HU, n = $(length(p)))"
    end
    out
end;

verification = let
    checks = Tuple{String, Bool, String}[]

    # 1. Regular (dual-power) water accuracy
    for (tag, s) in (("tube A", dp_stats.water.a), ("tube B", dp_stats.water.b),
        ("combined", dp_stats.water.combined))
        push!(checks, (
            "regular water $(tag)", abs(s.mean) ≤ 5.0,
            "⟨HU⟩ = $(round(s.mean, digits = 2)) (gate ±5)",
        ))
    end

    # 2. Dual-power √2 noise reduction
    push!(checks, (
        "dual-power σ ratio", 0.62 ≤ dp_stats.noise_ratio ≤ 0.80,
        "`σ_combined / σ_single` = $(round(dp_stats.noise_ratio, digits = 3)) (ideal 0.707, gate [0.62, 0.80])",
    ))

    # 3. DE poly water accuracy
    for (tag, s) in (("80 kVp", poly_water_stats.low), ("Sn140 kVp", poly_water_stats.high),
        ("mixed", poly_water_stats.mixed))
        push!(checks, (
            "DE poly water $(tag)", abs(s.mean) ≤ 5.0,
            "⟨HU⟩ = $(round(s.mean, digits = 2)) (gate ±5)",
        ))
    end

    # 4. VMI water accuracy per keV
    for E in de_vmi_energies
        m = vmi_noise_by_keV[E].mean
        push!(checks, (
            "VMI water @ $(Int(E)) keV", abs(m) ≤ 10.0,
            "⟨HU⟩ = $(round(m, digits = 2)) (gate ±10)",
        ))
    end

    # 5. Bias against the noise-free reference (same projection, same reconstruction)
    bias = [vmi_noise_by_keV[E].mean - vmi_noise_by_keV[E].reference for E in de_vmi_energies]
    push!(checks, (
        "VMI bias vs noise-free reference", maximum(abs, bias) ≤ 2.0,
        "⟨HU⟩ − ⟨HU⟩_ref = " * join([string(round(b, digits = 2)) for b in bias], ", ") * " (gate ±2)",
    ))

    # 6. The noise is flat in energy: the chain carries the composite's noise to every VMI
    σs = [vmi_noise_by_keV[E].std for E in de_vmi_energies]
    σn_all = [vmi_noise_by_keV[E].std_none for E in de_vmi_energies]
    flat = maximum(σs) / minimum(σs)
    push!(checks, (
        "VMI noise flat in energy", flat ≤ 1.10,
        "σ = " * join([string(round(σ, digits = 1)) for σ in σs], ", ") *
            " HU, max / min = $(round(flat, digits = 3)) (gate ≤ 1.10; no denoising " *
            "$(round(maximum(σn_all) / minimum(σn_all), digits = 2)))",
    ))

    # 7. The denoising lowers the noise at every keV
    σn = [vmi_noise_by_keV[E].std_none for E in de_vmi_energies]
    push!(checks, (
        "denoising lowers σ", all(σs .< σn),
        "σ / σ_none = " * join([string(round(a / b, digits = 2)) for (a, b) in zip(σs, σn)], ", "),
    ))

    # 8. The chain that ran
    chain = de_vmi.settings
    push!(checks, (
        "VMI chain",
        chain.n_channels == 2 && chain.denoiser !== nothing &&
            chain.denoiser.projection !== nothing && chain.denoiser.image !== nothing &&
            chain.acnr !== nothing && chain.acnr.on === :complement &&
            chain.recon.filter isa BS.PairFilter &&
            chain.pair.basis == calibration.pair_basis &&
            chain.pair.Estar == calibration.composite_energy &&
            FLASH.opts.view_samples == 5 && all(e -> e.view_arc == 1.0, FLASH.exposures),
        "K = $(chain.n_channels), projection + image HYPR, ACNR on the `$(chain.acnr === nothing ? "none" : chain.acnr.on)`, " *
            "`PairFilter`, calibration pair E* = $(round(chain.pair.basis.Estar, digits = 1)) keV, " *
            "composite $(round(chain.pair.Estar, digits = 1)) keV, `view_samples` = $(FLASH.opts.view_samples)",
    ))

    # 9. Per-rod regression gates
    function _fit(x, y)
        x̄ = mean(x); ȳ = mean(y)
        β = sum((x .- x̄) .* (y .- ȳ)) / sum((x .- x̄) .^ 2)
        α = ȳ - β * x̄
        ŷ = α .+ β .* x
        r² = 1 - sum((y .- ŷ) .^ 2) / sum((y .- ȳ) .^ 2)
        (slope = β, r² = r²)
    end
    for group in (:Ca, :I), (j, E) in pairs(de_vmi_energies)
        d = rod_data[group]
        f = _fit(Vector{Float64}(vec(d.theoretical[:, j])), Vector{Float64}(vec(d.measured[:, j])))
        push!(checks, (
            "$(group) regression @ $(Int(E)) keV",
            (0.85 ≤ f.slope ≤ 1.15) && (f.r² ≥ 0.99),
            "slope = $(round(f.slope, digits = 3)), R² = $(round(f.r², digits = 4)) (gates slope [0.85, 1.15], R² ≥ 0.99)",
        ))
    end

    n_pass = count(c -> c[2], checks)
    for (name, ok, detail) in checks
        @info "$(ok ? "✅ PASS" : "❌ FAIL") · $(name) — $(detail)"
    end
    overall = n_pass == length(checks)
    @info (overall ? "🎉 VERIFICATION PASS $(n_pass)/$(length(checks))" :
        "🚨 VERIFICATION FAIL $(n_pass)/$(length(checks))")

    (checks = checks, n_pass = n_pass, n_total = length(checks), pass = overall)
end;