QRM-Thorax Pure-Material PCCT VMI: Full-Resolution True-Scan Reference

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

using Markdown: @md_str, Markdown

using Statistics: mean, std, quantile

using Unitful: @u_str

using Random: MersenneTwister, randn!

import BasisSimulator as BS

# import CairoMakie as Mke
import WasmMakie 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

const QRM_CACHE_PATH = joinpath(@__DIR__, "data", "qrm_thorax", "qrm_thorax_1600x1100_rot_uint8.raw");

const QRM_TARGET_NX = 1600;  # full prepared cache, no extra in-plane downsample

const QRM_TARGET_NY = 1100;  # full prepared cache, no extra in-plane downsample

const QRM_TARGET_NZ = 20;    # 20 × 0.2 mm = 4 mm — short z-invariant reference phantom

const QRM_VOXEL_SIZE_CM = (0.02, 0.02, 0.02);   # (x, y, z) cm — 0.2 mm isotropic ground truth (320 × 220 × 4 mm physical extent)

const QRM_SHARED_DRIVE_DIR = "/Volumes/Molloilab/Shu Nie/water-lipid";

const QRM_SHARED_FULL_PATH = joinpath(QRM_SHARED_DRIVE_DIR, "qrm_thorax_3200x2200_rot_uint8.raw");

const QRM_SHARED_DOWN_PATH = joinpath(QRM_SHARED_DRIVE_DIR, "qrm_thorax_1600x1100_rot_uint8.raw");

mask_3d_raw = let
    isfile(QRM_CACHE_PATH) || error(
        "QRM-Thorax cache not found at $(QRM_CACHE_PATH).\n" *
            "Either:\n" *
            "  • copy the prepared cache from the lab volume:\n" *
            "      cp \"$(QRM_SHARED_DOWN_PATH)\" \"$(QRM_CACHE_PATH)\"\n" *
            "  • or run the prep notebook once to rebuild the cache from source."
    )
    cache_2d = reshape(read(QRM_CACHE_PATH), QRM_TARGET_NX, QRM_TARGET_NY)
    repeat(cache_2d; outer = (1, 1, QRM_TARGET_NZ))
end;

const ROD_HEART_CENTER_PX = (800, 654);   # ≈ (50%, 59%) of the 1600 × 1100 frame

const ROD_RADIUS_MM = 7.5;                  # mm — each rod radius (physical)

const ROD_OFFSET_MM = 25.0;                  # mm — heart-center → rod-center distance (physical)

mask_3d = let
    out = copy(mask_3d_raw)
    nx, ny, nz = size(out)

    px_mm = QRM_VOXEL_SIZE_CM[1] * 10
    cx_px, cy_px = ROD_HEART_CENTER_PX
    r_px = ROD_RADIUS_MM / px_mm
    o_px = ROD_OFFSET_MM / px_mm

    rod_specs = (
        (0.0, +o_px, UInt8(9)),    # North → basis_water
        (+o_px, 0.0, UInt8(10)),   # East  → basis_lipid
        (0.0, -o_px, UInt8(11)),   # South → basis_collagen
        (-o_px, 0.0, UInt8(12)),   # West  → gammex_472_i5_0
    )

    for (dx, dy, lab) in rod_specs
        rx, ry = cx_px + dx, cy_px + dy
        i_lo = max(1, floor(Int, rx - r_px))
        i_hi = min(nx, ceil(Int, rx + r_px))
        j_lo = max(1, floor(Int, ry - r_px))
        j_hi = min(ny, ceil(Int, ry + r_px))
        @inbounds for j in j_lo:j_hi, i in i_lo:i_hi
            ((i - rx)^2 + (j - ry)^2) <= r_px^2 || continue
            for k in 1:nz
                out[i, j, k] = lab
            end
        end
    end

    out
end;

materials_dict = Dict{Int, BS.XA.Material}(
    # Anatomy
    1 => BS.XA.Materials.air,
    2 => BS.XA.Materials.lung,
    3 => BS.XA.Materials.muscle,
    4 => BS.XA.Materials.corticalbone,
    5 => BS.XA.Materials.marrow_red,
    6 => BS.XA.Materials.water,
    7 => BS.XA.Materials.air,
    8 => BS.XA.Materials.water,
    # Pure-material rod inserts (relabeled from Ca-HA)
    9 => BS.XA.Materials.basis_water,
    10 => BS.XA.Materials.basis_lipid,
    11 => BS.XA.Materials.basis_collagen,
    12 => BS.XA.Materials.gammex_472_i5_0,
);

phantom_cpu = BS.create_phantom_from_mask(
    Array{Int, 3}(mask_3d),
    materials_dict,
    QRM_VOXEL_SIZE_CM,
);

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

QRM_LABEL_NAMES = Dict{Int, String}(
    1 => "1 air",
    2 => "2 lung",
    3 => "3 muscle / soft tissue",
    4 => "4 cortical bone",
    5 => "5 marrow",
    6 => "6 water rod (lung-side)",
    7 => "7 air rod",
    8 => "8 heart (water-filled)",
    9 => "9 → basis_water rod",
    10 => "10 → basis_lipid rod",
    11 => "11 → basis_collagen rod",
    12 => "12 → gammex_472_i5_0 rod",
);

scanner = let
    native_col_mm = 0.275
    native_row_mm = 0.322
    sid = 610.0
    sdd = 1113.0
    magnification = sdd / sid
    bf = 2

    pixel_col_iso = (native_col_mm * bf) / magnification
    pixel_row_iso = (native_row_mm * bf) / magnification
    n_cols = ceil(Int, 360.0 / pixel_col_iso)

    BS.Scanner(
        source_to_isocenter = sid,
        source_to_detector = sdd,

        detector_rows = 144,
        detector_cols = n_cols,
        detector_row_size = pixel_row_iso,
        detector_col_size = pixel_col_iso,
        detector_row_offset = 0.0,
        detector_col_offset = pixel_col_iso / 2,

        focal_spot_width = 0.4,
        focal_spot_length = 0.5,
        target_angle = 7.0,

        gantry_rotation_time = 0.5,
        scan_diameter = 360.0,
        gantry_aperture = 820.0,

        flat_filter_material = :aluminum,
        flat_filter_thickness = 3.0,
        # No bowtie — the physical Naeotom Alpha config (matches nb04).  The
        # decomposition uses the 1-D applied-W spectrum, so a per-pixel bowtie
        # would need the per-pixel-ŵ path to stay consistent; :none keeps the
        # setup coherent with nb04's certified chain.
        bowtie_filter = :none,

        detector_material = :cdte,
        detector_depth = 1.6,
        fill_factor_row = 0.95,
        fill_factor_col = 0.95,
        detection_gain = 1.0,
        electronic_noise = 0.0,

        detector_type = :photon_counting,
        n_energy_bins = 4,
        energy_thresholds = [20.0, 35.0, 55.0, 70.0],
        energy_resolution = 10.0,
        charge_sharing_fwhm = 0.08,
        dead_time_ns = 5.0,
        pixel_mode = :standard,

        native_dexel_col_mm = native_col_mm,
        native_dexel_row_mm = native_row_mm,
        binning_factor = bf,
    )
end;

protocol = BS.CTProtocol(
    kVp = 140,
    mA = 174.0,    # clinical PCCT dose (matches header + canonical nb04); was 5.0
    views = 1200,
    rotation_time = 0.5,
    collimation_mm = 5.0,    # nominal width; full-FOV axial cone guards are automatic
    additional_filters = [("Ti", 0.9)],
);

sim_opts = BS.SimOptions(
    fidelity = :pcct,
    seed = 1234,
    projector = :dd_fast,  # same anti-aliased DD physics, single-pass fused kernels (~47× faster poly)

    # ─── Inert for PCCT (flag exists but does nothing) ───
    use_fill_factor = false,
    use_detector_efficiency = false,
    use_optical_crosstalk = false,
    use_focal_spot = false,
    use_lag = false,
    use_heel_effect = false,

    # ─── Active for PCCT ───
    use_scatter = false,                  # EICT scatter flag — OFF (PCCT uses use_pcct_scatter)
    use_noise = true,                     # quantum noise inside simulate!()  (src :count, nr below)
    use_pcct_scatter = true,              # ← PCCT scatter injection, inside simulate!()
    use_pcct_scatter_correction = true,   # ← PCCT model-based scatter correction, inside simulate!()
    use_pcct_pileup = true,               # ← PCCT pileup forward, inside simulate!()
    use_pcct_pileup_correction = true,    # ← PCCT pileup correction (inverse S), inside simulate!()
    # DETECTOR-LEVEL CORRECTION SURROGATE — NOT a recon-level (QIR) stand-in
    # (chain is pure FBP; accuracy is independent of this knob).  Stands in
    # for the vendor's detector-side algorithms (anti-coincidence /
    # charge-sharing event reconstruction, count-rate linearization,
    # threshold compensation) whose degradations we Monte-Carlo simulate
    # but whose corrections we do not implement.
    pcct_noise_reduction = 0.7,
)

# Same saved recon grid as notebook 07: 512 × 512 in-plane at 0.625 mm
# isotropic, 3 slices. Automatic detector guards provide full-FOV support.
recon_opts = BS.ReconOptions(
    matrix_size = (512, 512, 3),
    fov_cm = 32.0,
    z_cm = 0.1875,
);

# === SLOW CELL (~6 min) — Full measured sinogram via STANDARD simulate!() ===
# Single src call producing the 4 per-bin log-line-integral sinograms with the
# COMPLETE PCCT physics + corrections, all gated by `sim_opts` flags:
#   forward → scatter inject (use_pcct_scatter) → quantum noise (use_noise,
#   pcct_noise_reduction) → pile-up fwd (use_pcct_pileup) → pile-up correction
#   (use_pcct_pileup_correction) → scatter correction (use_pcct_scatter_correction).
# Consumed directly by the bin-combine below — the former inline forward cells
# (pile-up / scatter fwd + correct) are GONE; this is the src-proper path.
#
sim_raw = let
    @info "simulate!(): $(Int(protocol.kVp)) kVp / $(round(protocol.mA, digits = 1)) mA — full PCCT physics + corrections (nr = $(sim_opts.pcct_noise_reduction))"
    ws = BS.create_workspace(scanner, protocol, sim_opts, recon_opts, phantom)
    result = BS.simulate!(ws, phantom, protocol, sim_opts)

    bins_raw = [Array(b) for b in result.pcct_sino.bins]   # full physics + corrections (scatter+noise+pileup, per sim_opts)
    I0_bins = copy(result.I0_bins)
    geom = ws.geom
    energies = Float64.(ws.energies)
    weights = copy(ws.weights)
    # The EXACT per-bin detected spectra the forward applied (w·η·DRM with the
    # workspace's MC-LUT η + centre-pixel bowtie fold) — the decomposition
    # basis consumes THESE, so the inversion's forward model is the model
    # simulate! actually applied (matches nb04).
    W_applied = Float64.(Array(ws.W_matrix_gpu))[1:length(ws.energies), :]
    bf = scanner.binning_factor

    ws = nothing; result = nothing
    GC.gc(true)
    (
        bins_raw = bins_raw, I0_bins = I0_bins, W_applied = W_applied,
        geom = geom, energies = energies, weights = weights, bf = bf,
    )
end;

# Resample the phantom labels onto the recon grid via BS's affine
# round-trip — used downstream for ROI construction in recon coords.
phantom_in_recon = BS.resample_to_recon(
    phantom_cpu, sim_raw.geom, recon_opts.matrix_size; method = :nearest,
);

sim_lohi = let
    low_bins = [1, 2, 3]     # 123|4 partition — analytic noise optimum (matches nb04)
    high_bins = [4]

    I0_lo = Float32(sum(Float64.(sim_raw.I0_bins[low_bins])))
    I0_hi = Float32(sum(Float64.(sim_raw.I0_bins[high_bins])))

    sz = size(sim_raw.bins_raw[1])
    N_lo = zeros(Float32, sz)
    N_hi = zeros(Float32, sz)
    for b in low_bins
        I0b = Float32(sim_raw.I0_bins[b])
        @. N_lo += I0b * exp(-sim_raw.bins_raw[b])
    end
    for b in high_bins
        I0b = Float32(sim_raw.I0_bins[b])
        @. N_hi += I0b * exp(-sim_raw.bins_raw[b])
    end

    # Physical DAS floor (1 count) + first-order log-Poisson DEBIAS (matches
    # nb04): E[−log(N/I0)] = p_true + s²/(2N), s = 1 − pcct_noise_reduction.
    # Deterministic bias correction (σ untouched).
    s_nr = Float32((1.0 - sim_opts.pcct_noise_reduction)^2)
    N_lo_f = max.(N_lo, 1.0f0)
    N_hi_f = max.(N_hi, 1.0f0)
    sino_low = Float32.(.- log.(N_lo_f ./ I0_lo) .- s_nr ./ (2.0f0 .* N_lo_f))
    sino_high = Float32.(.- log.(N_hi_f ./ I0_hi) .- s_nr ./ (2.0f0 .* N_hi_f))

    (
        sino_low = sino_low, sino_high = sino_high,
        I0_lo = I0_lo, I0_hi = I0_hi,
        geom = sim_raw.geom,
    )
end;

# Bin-combine partition feeding the two Cong channels.  Must match the
# §6 combine — change here AND there together.
begin
    cong_low_bins = 1:3          # PCCT bins forming the "low"  channel (123|4, matches §6)
    cong_high_bins = 4:4          # PCCT bins forming the "high" channel
end

material_basis = let
    # Per-channel spectra = column sums of the sim's APPLIED W matrix over
    # the bin partition — exact by construction, identical to nb04.  (The
    # forward applies the MC-LUT η + centre-pixel bowtie fold; rebuilding
    # w·η·R with the analytic η here would mismatch it and reintroduce a
    # keV-dependent HU bias.  With no bowtie the per-pixel ŵ collapses to
    # this 1-D spectrum, so this is the exact effective model.)
    e = sim_raw.energies
    ΦL = Float32.(vec(sum(sim_raw.W_applied[:, collect(cong_low_bins)]; dims = 2)))
    ΦH = Float32.(vec(sum(sim_raw.W_applied[:, collect(cong_high_bins)]; dims = 2)))
    ŵ_L_f32 = ΦL ./ sum(ΦL)
    ŵ_H_f32 = ΦH ./ sum(ΦH)

    p = Float32[Float32(BS.compute_mass_μ_at_energy(BS.XA.Elements.Iodine, Float64(E))) for E in e]
    q = Float32[Float32(BS.compute_mass_μ_at_energy(BS.XA.Materials.water, Float64(E))) for E in e]

    @info "[Cong basis · applied-W] low ⟨E⟩ = $(round(sum(e .* Float64.(ŵ_L_f32)); digits = 1)) keV · " *
        "high ⟨E⟩ = $(round(sum(e .* Float64.(ŵ_H_f32)); digits = 1)) keV"

    (
        ŵ_L = ŵ_L_f32, p_L = p, q_L = q,
        ŵ_H = ŵ_H_f32, p_H = copy(p), q_H = copy(q),
    )
end;

sino_basis = let
    sino_low_gpu = to_gpu(Float32.(sim_lohi.sino_low))
    sino_high_gpu = to_gpu(Float32.(sim_lohi.sino_high))

    sino_y = similar(sino_low_gpu)
    sino_c = similar(sino_low_gpu)
    fill!(sino_y, 0.0f0); fill!(sino_c, 0.0f0)

    @info "Cong polychromatic decomposition: $(size(sim_lohi.sino_low))"
    cong_elapsed = @elapsed begin
        cong_ws = BS.create_cong_workspace(sino_low_gpu, material_basis)
        BS.apply_cong!(
            cong_ws, sino_y, sino_c, sino_low_gpu, sino_high_gpu;
            water_basis = (a = 0.0f0, c = 1.0f0),
        )
    end
    @info "Cong done in $(round(cong_elapsed, digits = 1)) s"

    result = (
        sino_iodine = Array(sino_y),
        sino_water = Array(sino_c),
        geom = sim_lohi.geom,
    )
    sino_low_gpu = nothing; sino_high_gpu = nothing
    sino_y = nothing; sino_c = nothing; cong_ws = nothing
    GC.gc(true)
    result
end;

# === Basis-sinogram covariance readout (sinogram domain, PRE-FBP) ===
# Measures the iodine↔water correlation in the basis SINOGRAMS that feed FBP — the
# quantity an ACNR-before-FBP step would act on.  Two views:
#   • GLOBAL    — cov over all rays.  With noise OFF this is the deterministic
#                 STRUCTURE correlation Cong bakes in; with noise ON it also carries it.
#   • Δcol HF   — cov of the adjacent-detector-column difference (isolates the
#                 high-frequency / noise-like component — what ACNR actually targets).
# ρ is scale-free (iodine sino tiny, water sino large), so it is the honest correlation.
# ρ < 0 ⇒ anti-correlated basis content ⇒ ACNR has something to remove.  Single-pass
# scalar accumulation — no full-sinogram temporaries (memory budget).
let
    io = sino_basis.sino_iodine
    wa = sino_basis.sino_water
    n = length(io)

    # ── GLOBAL covariance (two scalar passes) ──
    si = 0.0; sw = 0.0
    @inbounds for k in eachindex(io)
        si += io[k]; sw += wa[k]
    end
    mi = si / n; mw = sw / n
    Vi = 0.0; Vw = 0.0; Ciw = 0.0
    @inbounds for k in eachindex(io)
        a = io[k] - mi; b = wa[k] - mw
        Vi += a * a; Vw += b * b; Ciw += a * b
    end
    Vi /= n; Vw /= n; Ciw /= n
    ρ_g = Ciw / sqrt(max(Vi * Vw, 1.0e-30))

    # ── HIGH-FREQ covariance: adjacent-detector-column difference Δcol ──
    nc, nr, nv = size(io)
    sdi = 0.0; sdw = 0.0; ndf = 0
    @inbounds for v in 1:nv, r in 1:nr, c in 2:nc
        sdi += io[c, r, v] - io[c - 1, r, v]
        sdw += wa[c, r, v] - wa[c - 1, r, v]
        ndf += 1
    end
    mdi = sdi / ndf; mdw = sdw / ndf
    Vih = 0.0; Vwh = 0.0; Ciwh = 0.0
    @inbounds for v in 1:nv, r in 1:nr, c in 2:nc
        a = (io[c, r, v] - io[c - 1, r, v]) - mdi
        b = (wa[c, r, v] - wa[c - 1, r, v]) - mdw
        Vih += a * a; Vwh += b * b; Ciwh += a * b
    end
    Vih /= ndf; Vwh /= ndf; Ciwh /= ndf
    ρ_hf = Ciwh / sqrt(max(Vih * Vwh, 1.0e-30))

    @info "[basis-sino cov · GLOBAL]  σ_iod = $(round(sqrt(Vi), sigdigits = 3)) · σ_wat = $(round(sqrt(Vw), sigdigits = 3)) g/cm² · ρ(iod,wat) = $(round(ρ_g, digits = 3))"
    @info "[basis-sino cov · Δcol HF] σ_iod = $(round(sqrt(Vih), sigdigits = 3)) · σ_wat = $(round(sqrt(Vwh), sigdigits = 3)) g/cm² · ρ_hf(iod,wat) = $(round(ρ_hf, digits = 3))  ← ACNR target"
    @info "  ρ<0 ⇒ anti-correlated basis content (ACNR has something to remove).  Noise is OFF → these are the STRUCTURE baseline; turn channel noise ON to read the NOISE correlation."
    (ρ_global = ρ_g, ρ_highfreq = ρ_hf, V_iodine = Vi, V_water = Vw, C_iw = Ciw)
end

basis_volumes = let
    matrix_size = recon_opts.matrix_size
    geom = sino_basis.geom

    function _fbp(sino_cpu)
        sino_gpu = to_gpu(Float32.(sino_cpu))
        ws = BS.create_fdk_recon_workspace(
            sino_gpu, geom, matrix_size; filter = BS.SoftFilter(),
        )
        recon = Array(BS.reconstruct!(ws, sino_gpu, geom))
        ws = nothing; sino_gpu = nothing
        GC.gc(true)
        return Float32.(recon)
    end

    (
        vol_iodine_raw = _fbp(sino_basis.sino_iodine),
        vol_water_raw = _fbp(sino_basis.sino_water),
        geom = geom,
    )
end;

# Image-domain cov-ACNR on the FBP basis maps via the src-proper
# `BS.apply_acnr_kalender!` (denoising/acnr.jl) — no knobs, no blur.
basis_acnr = let
    APPLY_ACNR = true        # ON — image-domain edge-aware cov-ACNR

    W = copy(basis_volumes.vol_water_raw)
    I = copy(basis_volumes.vol_iodine_raw)

    if APPLY_ACNR
        info = BS.apply_acnr_kalender!(W, I)
        @info "[ACNR · Kalender-1988 true ACNR] ρ_hp(W,I)=$(round(info.ρ_hp, digits = 3))"
    else
        @info "[ACNR] OFF (passthrough)"
    end

    (vol_iodine = I, vol_water = W, geom = basis_volumes.geom)
end;

solid_water_basis = let
    WATER_ROD_LABEL = UInt8(9)
    ROI_RADIUS_PX = 8

    mask_2d = phantom_in_recon[:, :, size(phantom_in_recon, 3) ÷ 2 + 1]
    nx_r, ny_r, nz_r = size(basis_acnr.vol_water)

    rod_idx = findall(==(WATER_ROD_LABEL), mask_2d)
    isempty(rod_idx) && error(
        "solid_water_basis: no label-$(Int(WATER_ROD_LABEL)) (water rod) voxels in resampled phantom mask."
    )
    cx = sum(ci -> Float64(ci[1]), rod_idx) / length(rod_idx)
    cy = sum(ci -> Float64(ci[2]), rod_idx) / length(rod_idx)

    sw_bool = falses(nx_r, ny_r)
    r² = Float64(ROI_RADIUS_PX)^2
    @inbounds for j in 1:ny_r, i in 1:nx_r
        ((i - cx)^2 + (j - cy)^2) ≤ r² && (sw_bool[i, j] = true)
    end

    n_voxels = count(sw_bool)
    @info "solid_water_basis: water-rod core ROI center = " *
        "($(round(cx, digits = 1)), $(round(cy, digits = 1))) recon voxels, " *
        "$(n_voxels) voxels in core"

    sw_idx = findall(sw_bool)
    function _mean(vol)
        s = 0.0; n = 0
        for z in 1:nz_r, ci in sw_idx
            s += vol[ci, z]; n += 1
        end
        return s / n
    end

    c_w = Float64(_mean(basis_acnr.vol_water))
    c_i = Float64(_mean(basis_acnr.vol_iodine))
    @info "solid_water_basis: ⟨c_water⟩_water-rod = $(round(c_w, digits = 4)) g/cm³, " *
        "⟨c_iodine⟩_water-rod = $(round(c_i, digits = 6)) g/cm³"

    (
        c_water = c_w, c_iodine = c_i, n_voxels = length(sw_idx) * nz_r,
        mask_2d = collect(sw_bool),
    )
end;

# Per-rod basis-decomp sanity check.  For each of the 4 rods, measure
# the mean (c_water, c_iodine) inside an 8-px-radius core ROI at the
# rod centroid and compare to the rod material's density.
#
# Expected:
#   rod 9  (basis_water)    → c_water ≈ 1.00 g/cm³, c_iodine ≈ 0
#   rod 10 (basis_lipid)    → c_water ≈ 0.92,        c_iodine ≈ 0
#   rod 11 (basis_collagen) → c_water ≈ 1.26,        c_iodine ≈ 0
#   rod 12 (gammex_472_i5_0)→ c_water ≈ 1.0,         c_iodine ≈ 0.005 (5 mg/mL)
let
    mask_2d = phantom_in_recon[:, :, size(phantom_in_recon, 3) ÷ 2 + 1]
    nx, ny = size(mask_2d)
    ROI_R_PX = 8

    n_z = size(basis_acnr.vol_water, 3)

    function rod_centroid(label::UInt8)
        idx = findall(==(label), mask_2d)
        cx = sum(ci -> Float64(ci[1]), idx) / length(idx)
        cy = sum(ci -> Float64(ci[2]), idx) / length(idx)
        return (cx, cy)
    end

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

    function _mean(vol, roi)
        s = 0.0; n = 0
        for z in 1:n_z, ci in roi
            s += vol[ci, z]; n += 1
        end
        return s / n
    end

    @info "Per-rod basis-decomp diagnostic (8-px core ROI in RECON coords, mean over z):"
    for (lab, name, mat) in zip(ROD_LABELS, ROD_NAMES, ROD_MATERIALS)
        cx, cy = rod_centroid(lab)
        roi = _disc_idx(cx, cy)
        c_w = _mean(basis_acnr.vol_water, roi)
        c_i = _mean(basis_acnr.vol_iodine, roi)
        ρ = round(BS.XA.val(mat.density), digits = 3)
        @info "  rod $(lab) ($(rpad(name, 9))): c_water = $(round(c_w, digits = 4)) g/cm³, " *
            "c_iodine = $(round(c_i, digits = 6)) g/cm³  (truth ρ = $(ρ))"
    end
end

pcct_vmi_energies = [40.0, 70.0, 100.0, 140.0];

vmi_HU_final = let
    c_iodine_mg_per_mL = basis_acnr.vol_iodine .* 1000.0f0

    out = Dict{Float64, Array{Float32, 3}}()
    for E in pcct_vmi_energies
        μρ_w = BS.compute_mass_μ_at_energy(BS.XA.Materials.water, E)
        μρ_I = BS.compute_mass_μ_at_energy(BS.XA.Elements.Iodine, E)
        μ_water_anchor = solid_water_basis.c_water * μρ_w +
            solid_water_basis.c_iodine * μρ_I
        Δ_pct = 100.0 * (μ_water_anchor - μρ_w) / μρ_w
        @info "VMI synth @ $(Int(E)) keV: divisor = $(round(μρ_w, digits = 5)) cm⁻¹ " *
            "(mono μρ_water);  SW-ROI anchor = $(round(μ_water_anchor, digits = 5)) " *
            "→ Δ = $(round(Δ_pct, digits = 2))%"

        out[E] = BS.synth_vmi_2basis(
            basis_acnr.vol_water, c_iodine_mg_per_mL; energy_keV = E,
        )
    end
    out
end;

const HEART_NOISE_ROI_OFFSET_PX = (0, -40);   # (dx, dy) recon vx — default: below water-rod centroid

const HEART_NOISE_ROI_RADIUS_PX = 12;   # ≈7.5 mm at 0.625 mm/voxel

heart_noise_roi = let
    WATER_ROD_LABEL = UInt8(9)

    mask_2d = phantom_in_recon[:, :, size(phantom_in_recon, 3) ÷ 2 + 1]
    nx_r, ny_r, nz_r = size(basis_acnr.vol_water)

    rod_idx = findall(==(WATER_ROD_LABEL), mask_2d)
    isempty(rod_idx) && error(
        "heart_noise_roi: no label-$(Int(WATER_ROD_LABEL)) (water rod) voxels in resampled phantom mask."
    )
    rod_cx = sum(ci -> Float64(ci[1]), rod_idx) / length(rod_idx)
    rod_cy = sum(ci -> Float64(ci[2]), rod_idx) / length(rod_idx)

    dx, dy = HEART_NOISE_ROI_OFFSET_PX
    cx = rod_cx + dx
    cy = rod_cy + dy

    roi_bool = falses(nx_r, ny_r)
    r² = Float64(HEART_NOISE_ROI_RADIUS_PX)^2
    @inbounds for j in 1:ny_r, i in 1:nx_r
        ((i - cx)^2 + (j - cy)^2) ≤ r² && (roi_bool[i, j] = true)
    end

    n_vox = count(roi_bool)
    @info "heart_noise_roi: water-rod centroid = ($(round(rod_cx, digits = 1)), $(round(rod_cy, digits = 1))), " *
        "offset = $(HEART_NOISE_ROI_OFFSET_PX), ROI center = ($(round(cx, digits = 1)), $(round(cy, digits = 1))), " *
        "$(n_vox) vx × $(nz_r) z = $(n_vox * nz_r) total"

    (
        rod_center_xy = (rod_cx, rod_cy),
        center_xy = (cx, cy),
        mask_2d = roi_bool,
        n_voxels = n_vox,
        n_total = n_vox * nz_r,
    )
end;

# === Image-domain basis covariance (heart ROI) — the VMI-noise theorem inputs ===
# The quantity that ACTUALLY governs the U:
#   σ_HU(E)² = 1e6·V_w + α(E)²·V_i + 2e3·α(E)·C_iw ,   min at  α* = −1e3·C_iw/V_i
#   monotonic-DECREASING ⟺ α* ≤ α(140) ⟺ C_iw not too negative.
# Measured over the heart ROI on the ACNR basis maps VMI synth consumes
# (c_iodine in mg/mL).  ρ_basis<0 = anti-correlated basis noise = the U + the ACNR
# target.  Small ROI → tiny vectors, no memory concern.
let
    roi = findall(heart_noise_roi.mask_2d)
    nz = size(basis_acnr.vol_water, 3)
    cw = Float64[Float64(basis_acnr.vol_water[p, z])           for z in 1:nz, p in roi]
    ci = Float64[Float64(basis_acnr.vol_iodine[p, z]) * 1000.0 for z in 1:nz, p in roi]  # mg/mL
    mw = mean(cw); mi = mean(ci)
    Vw = mean((cw .- mw) .^ 2); Vi = mean((ci .- mi) .^ 2)
    Ciw = mean((cw .- mw) .* (ci .- mi))
    ρ_b = Ciw / sqrt(max(Vw * Vi, 1.0e-30))
    αf(E) = Float64(BS.compute_mass_μ_at_energy(BS.XA.Elements.Iodine, E)) /
        Float64(BS.compute_mass_μ_at_energy(BS.XA.Materials.water, E))
    α_star = -1000.0 * Ciw / max(Vi, 1.0e-30)
    verdict = α_star ≤ αf(140.0) ? "MONOTONIC-decreasing predicted ✓" :
        α_star ≥ αf(40.0) ? "monotonic-INCREASING" :
        "U-shape (min near α=α*)"
    @info "[basis cov · heart ROI, image]  σ_water = $(round(sqrt(Vw), sigdigits = 3)) g/mL · " *
        "σ_iod = $(round(sqrt(Vi), sigdigits = 3)) mg/mL · ρ_basis = $(round(ρ_b, digits = 3)) · " *
        "water floor = $(round(1000 * sqrt(Vw), digits = 1)) HU"
    @info "  α* = $(round(α_star, digits = 1))  vs  α(40)=$(round(αf(40.0), digits = 1)), " *
        "α(70)=$(round(αf(70.0), digits = 1)), α(140)=$(round(αf(140.0), digits = 1))  →  $(verdict)"
    (V_water = Vw, V_iodine = Vi, C_iw = Ciw, ρ_basis = ρ_b, α_star = α_star)
end

vmi_noise_by_keV = let
    roi_idx = findall(heart_noise_roi.mask_2d)
    nz_r = size(basis_acnr.vol_water, 3)

    out = Dict{Float64, NamedTuple}()
    for E in pcct_vmi_energies
        vol = vmi_HU_final[E]
        vals = Float64[Float64(vol[ci, z]) for z in 1:nz_r, ci in roi_idx]
        μ = mean(vals)
        σ = std(vals)
        out[E] = (mean = μ, std = σ, n = length(vals))
        @info "heart noise @ $(Int(E)) keV: ⟨HU⟩ = $(round(μ, digits = 2)),  σ = $(round(σ, digits = 2)) HU  (n = $(length(vals)))"
    end
    out
end;

const ROD_LABELS = (UInt8(9), UInt8(10), UInt8(11), UInt8(12));

const ROD_NAMES = ("water", "lipid", "collagen", "iodine_5");

const ROD_MATERIALS = (
    BS.XA.Materials.basis_water,
    BS.XA.Materials.basis_lipid,
    BS.XA.Materials.basis_collagen,
    BS.XA.Materials.gammex_472_i5_0,
);

rod_data = let
    mask_2d = phantom_in_recon[:, :, size(phantom_in_recon, 3) ÷ 2 + 1]
    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 ROD_LABELS)

    μ_water_E = Dict(
        E => BS.compute_μ_at_energy(BS.XA.Materials.water, E)
            for E in pcct_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

    n_rods = length(ROD_LABELS)
    n_E = length(pcct_vmi_energies)
    meas = zeros(Float64, n_rods, n_E)
    theo = zeros(Float64, n_rods, n_E)
    for (i, lab) in enumerate(ROD_LABELS)
        mat = ROD_MATERIALS[i]
        for (j, E) in enumerate(pcct_vmi_energies)
            meas[i, j] = measured_hu(vmi_HU_final[E], lab)
            theo[i, j] = theoretical_hu(mat, E)
        end
    end
    (
        labels = ROD_LABELS, names = ROD_NAMES,
        measured = meas, theoretical = theo,
    )
end;