Photon-Counting CT Virtual Monoenergetic Imaging

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

using Markdown: @md_str, Markdown

using Statistics: mean, std, quantile, median

import PlutoUI

import BasisSimulator as BS

# import CairoMakie as Mke
import WasmMakie as Mke

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,
);

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,

        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,
    views = 1200,
    rotation_time = 0.5,
    collimation_mm = 5.0,
    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 — all applied INSIDE simulate!() ───
    use_scatter = false,                  # EICT scatter flag — OFF (PCCT uses use_pcct_scatter)
    use_noise = true,                     # quantum noise (src :count, nr below)
    use_pcct_scatter = true,              # PCCT scatter injection
    use_pcct_scatter_correction = true,   # PCCT model-based scatter correction
    use_pcct_pileup = true,               # PCCT MC pile-up forward
    use_pcct_pileup_correction = true,    # PCCT pile-up correction (inverse S)
    # DETECTOR-LEVEL CORRECTION SURROGATE — explicitly NOT a recon-level
    # (QIR/iterative) stand-in: this chain is pure FBP end to end, and its
    # ACCURACY does not depend on this knob (pure chain, nr = 0, noise off:
    # rods within ~3 % of NIST, solid water < 3 HU).  The simulator
    # Monte-Carlo models the detector DEGRADATIONS (charge sharing,
    # fluorescence escape, pulse pileup, spectral distortion via the MC DRM)
    # but not the vendor's DETECTOR-side correction algorithms for them —
    # anti-coincidence/charge-sharing event reconstruction, count-rate
    # linearization beyond our inverse-S, threshold/spectral-distortion
    # compensation.  Those corrections recover count statistics at the
    # detector output; nr = 0.7 stands in for that recovery and nothing
    # else.
    pcct_noise_reduction = 0.7,
)

recon_opts = let
    slice_thickness_mm = 0.4
    n_recon_slices = max(1, round(Int, protocol.collimation_mm / slice_thickness_mm))
    BS.ReconOptions(
        matrix_size = (512, 512, n_recon_slices),
        fov_cm = 35.0,
        z_cm = protocol.collimation_mm / 10.0,
    )
end;

# === Forward project + full PCCT physics + corrections via simulate!() ===
# One src call: forward → scatter inject → noise → pile-up fwd → pile-up
# correction → scatter correction, all gated by the `sim_opts` flags.
sim_bins = let
    @info "Simulating: $(Int(protocol.kVp)) kVp / $(round(protocol.mA, digits = 1)) mA (PCCT 4-bin) — full physics + corrections via simulate!()"
    ws = BS.create_workspace(scanner, protocol, sim_opts, recon_opts, phantom)
    result = BS.simulate!(ws, phantom, protocol, sim_opts)

    bins = [Array(b) for b in result.pcct_sino.bins]
    I0_bins = copy(result.I0_bins)
    geom = ws.geom
    # The EXACT per-bin detected spectra the forward applied (w·η·DRM with
    # the workspace's MC-LUT η and the centre-pixel bowtie fold).  The
    # decomposition basis consumes THESE — the inversion's forward model is
    # the model simulate! actually applied, by construction.
    energies = Float64.(ws.energies)
    W_applied = Float64.(Array(ws.W_matrix_gpu))[1:length(ws.energies), :]

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

sim_bins_fbp = let
    matrix_size = recon_opts.matrix_size
    geom = sim_bins.geom

    fdk_filter = BS.CustomFilter(
        (0.0, 0.25, 0.5, 0.75, 1.0),
        (1.0, 0.75, 0.6, 0.2, 0.001),
    )

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

    [_fbp(b) for b in sim_bins.bins]
end;

bins_denoised = let
    APPLY_SVD = true   # A/B TOGGLE — false ⇒ passthrough raw bins
    if APPLY_SVD
        out = BS.apply_sino_svd_denoise_bilateral(
            sim_bins.bins;
            bilat_radius  = 3,     # search-window radius (px)
            bilat_sigma_s = 2.0,   # spatial Gaussian σ (px)
            bilat_range_k = 2.0,   # edge sensitivity (↓ protects edges harder)
        )
        (bins = out, I0_bins = sim_bins.I0_bins, geom = sim_bins.geom)
    else
        (bins = sim_bins.bins, I0_bins = sim_bins.I0_bins, geom = sim_bins.geom)
    end
end;

sino_combined = let
    low_bins = collect(1:3)
    high_bins = [4]

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

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

    N_lo_f = max.(N_lo, 1.0f0)   # physical DAS floor (1 count)
    N_hi_f = max.(N_hi, 1.0f0)
    sino_low = Float32.(.- log.(N_lo_f ./ I0_lo))
    sino_high = Float32.(.- log.(N_hi_f ./ I0_hi))

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

sim_lohi = let
    SVD2_SIGMA = 0.5f0   # subtle 2-channel Gaussian σ (px) on the combined pair
    out = BS.apply_sino_svd_denoise(
        [sino_combined.sino_low, sino_combined.sino_high]; σ_px = SVD2_SIGMA,
    )
    (
        sino_low = out[1], sino_high = out[2],
        I0_lo = sino_combined.I0_lo, I0_hi = sino_combined.I0_hi,
        geom = sino_combined.geom,
    )
end;

# Bin-combine partition feeding the two Cong channels.  Must match the
# `_combine` calls in §7 — change here AND there together.
begin
    # Partition from analytic noise optimization (iodine σ_y and water σ_c
    # over candidate 2-channel splits behind 33 cm water): moving bin 3
    # (55–70 keV — spectrally muddy in the high channel) into LOW raises the
    # separation determinant 0.392 → 0.523, cutting σ_y ≈ 16 % and σ_c ≈ 15 %
    # while keeping every photon.
    low_bins = 1:3          # PCCT bins forming the "low"  channel
    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.  (The old manual w·η·R
    # rebuild used the analytic η while the sim applies the MC-LUT η and
    # the bowtie-centre fold; together with the pre-air-cal normalization
    # this produced the historical +28…47 HU solid-water bias.)
    e = sim_bins.energies
    ΦL = Float32.(vec(sum(sim_bins.W_applied[:, collect(low_bins)]; dims = 2)))
    ΦH = Float32.(vec(sum(sim_bins.W_applied[:, collect(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)   # iodine basis line integrals (g/cm²)
    sino_c = similar(sino_low_gpu)   # water  basis line integrals (g/cm²)
    fill!(sino_y, 0.0f0); fill!(sino_c, 0.0f0)

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

    sino_iodine_cpu = Array(sino_y)
    sino_water_cpu = Array(sino_c)
    @info "[Cong decomp] ⟨∫ρ_I·dr⟩ = $(round(mean(sino_iodine_cpu), sigdigits = 4)) g/cm²   " *
        "⟨∫ρ_W·dr⟩ = $(round(mean(sino_water_cpu), sigdigits = 4)) g/cm²"

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

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

    # PER-BASIS apodization (NOT per-keV).  Applied ONCE to each basis sinogram;
    # every VMI is VMI_E = W + α(E)·I.  A SOFT iodine filter crushes the low-keV
    # noise (α(50) large); a SOFT water filter holds down the σ_W floor that
    # dominates the high-keV VMIs.  Selectivity is emergent from α(E).
    #
    # ── TUNE iodine softness here: lower mid/high values ⇒ softer ⇒ less
    #    low-keV noise (but softer iodine detail).
    iodine_filter = BS.CustomFilter(
        (0.0, 0.25, 0.5, 0.75, 1.0),
        (1.0, 0.40, 0.12, 0.03, 0.001),
    )
    water_filter = BS.SoftFilter()

    function _fbp(sino_cpu, filt)
        sino_gpu = to_gpu(Float32.(sino_cpu))
        ws = BS.create_fdk_recon_workspace(
            sino_gpu, geom, matrix_size; filter = filt,
        )
        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, iodine_filter),
        vol_water_raw = _fbp(sino_basis.sino_water, water_filter),
        geom = geom,
    )
end;

# Kalender-1988 true ACNR on the FBP basis maps via `BS.apply_acnr_kalender!`.
basis_acnr = let
    APPLY_ACNR = true          # Kalender-1988 true ACNR; false ⇒ passthrough
    ACNR_PASSES = 4            # geometric shrink of residual anti-correlation per pass
    ACNR_BETA_MAX = 20.0       # ← THE monotonicity lever: caps the per-pixel regression
    #   coefficient.  The default (8) UNDER-corrects, leaving Cov≈−6.3 when
    #   monotonicity needs |Cov| ≤ 5σ_I² ≈ 2.7 (α★ sat at ~100 keV → 140 keV tail
    #   flipped up).  Raising it lets strongly anti-correlated pixels actually be
    #   regressed out, pulling Cov under the bar in few passes.

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

    if APPLY_ACNR
        # TRUE ACNR (Kalender 1988): per-pixel local regression between the
        # two maps' high-frequency channels — anti-correlated (noise) content
        # subtracted exactly, structure pixels clamped to zero correction and
        # bit-untouched.  Zero blur by construction.
        info = BS.apply_acnr_kalender!(W, I; passes = ACNR_PASSES, beta_max = ACNR_BETA_MAX)
        @info "[ACNR · Kalender-1988 true ACNR] ρ_hp(W,I)=$(round(info.ρ_hp, digits = 3)) · σ_hp(W)=$(round(info.σ_hW, sigdigits = 3)) σ_hp(I)=$(round(info.σ_hI, sigdigits = 3))"
    else
        @info "[ACNR] OFF (passthrough)"
    end

    (vol_iodine_raw = I, vol_water_raw = W, geom = basis_volumes.geom)
end;

solid_water_basis = let
    # Mid-slice SW mask, broadcast across all recon z (Gammex 472 SW
    # background is z-invariant).  Then DEEPLY erode (σ = 12 px ≈ 8 mm
    # at 0.683 mm/px) so we sample only deep-interior SW voxels — well
    # clear of rod edges where partial-volume mixing with iodine / Ca
    # contaminates the basis means.
    ERODE_PX = 12.0

    mask_2d_raw = phantom_cpu.mask[:, :, size(phantom_cpu.mask, 3) ÷ 2]
    sw_bool_raw = (mask_2d_raw .== 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."
    )
    @info "solid_water_basis: SW mid-slice voxel count $(n_raw) → $(n_eroded) " *
        "after $(ERODE_PX)-px erosion"

    sw_idx = findall(sw_bool)
    n_z = size(basis_acnr.vol_water_raw, 3)
    function _mean(vol)
        s = 0.0; n = 0
        for z in 1:n_z, ci in sw_idx
            s += vol[ci, z]; n += 1
        end
        return s / n
    end

    c_w = Float64(_mean(basis_acnr.vol_water_raw))
    c_i = Float64(_mean(basis_acnr.vol_iodine_raw))
    @info "solid_water_basis: ⟨c_water⟩_SW = $(round(c_w, digits = 4)) g/cm³, " *
        "⟨c_iodine⟩_SW = $(round(c_i, digits = 6)) g/cm³"

    (
        c_water = c_w, c_iodine = c_i, n_voxels = length(sw_idx) * n_z,
        mask_2d = collect(sw_bool),
    )   # for downstream viz
end;

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

vmi_HU_by_keV = let
    # `BS.synth_vmi_2basis` expects c_iodine in mg/mL; our basis maps
    # are in g/cm³ (= g/mL).  Multiply by 1000 to convert.
    c_iodine_mg_per_mL = basis_acnr.vol_iodine_raw .* 1000.0f0

    out = Dict{Float64, Array{Float32, 3}}()
    for E in pcct_vmi_energies
        # Diagnostic-only — the SW-ROI synth-evaluated μ_water vs the
        # textbook mono divisor.  Drift = residual basis-decomp bias.
        μρ_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_raw, c_iodine_mg_per_mL;
            energy_keV = E,
        )
    end
    out
end;

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 = phantom_cpu.mask[:, :, size(phantom_cpu.mask, 3) ÷ 2]
    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 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

    out = Dict{Symbol, NamedTuple}()
    for group in (:Ca, :I)
        labels = ROD_LABELS[group]
        n_rods = length(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 pairs(labels)
            mat = materials[Int(lab) + 1]   # mask_value + 1
            for (j, E) in pairs(pcct_vmi_energies)
                meas[i, j] = measured_hu(vmi_HU_by_keV[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;

const WATER_NOISE_ROI_RADIUS_PX = 12;   # ≈8.2 mm at 0.683 mm/px (FOV 35 cm / 512)

# Central circular noise ROI in the solid-water background (image center =
# isocenter = phantom center for the centered Gammex 472).
water_noise_roi = let
    nx_r, ny_r, nz_r = size(basis_acnr.vol_water_raw)
    cx = nx_r ÷ 2 + 1
    cy = ny_r ÷ 2 + 1

    roi_bool = falses(nx_r, ny_r)
    r² = Float64(WATER_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 "water_noise_roi: center = ($(cx), $(cy)), radius = $(WATER_NOISE_ROI_RADIUS_PX) px, " *
        "$(n_vox) vx × $(nz_r) z = $(n_vox * nz_r) total"

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

# Per-keV HU noise (σ) + mean over the LARGE ERODED solid-water region
# (`solid_water_basis.mask_2d`) — the canonical water ROI for ALL noise
# measurements (not the tiny central circle).
vmi_noise_by_keV = let
    roi_idx = findall(solid_water_basis.mask_2d)
    nz_r = size(vmi_HU_by_keV[70.0], 3)

    out = Dict{Float64, NamedTuple}()
    for E in pcct_vmi_energies
        vol = vmi_HU_by_keV[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 "water-region noise @ $(Int(E)) keV: ⟨HU⟩ = $(round(μ, digits = 2)),  σ = $(round(σ, digits = 2)) HU  (n = $(length(vals)))"
    end
    out
end;