Virtual monoenergetic imaging from any set of spectral channels (dual-kVp, dual-source or photon-counting bins). vmi_pipeline is the whole chain: an optional SpectralHYPR denoiser, channel merging, K-channel (:nchannel) or Cong (:cong) projection-domain decomposition, basis reconstruction, and synthesis of the VMI stack. Each stage is also public on its own.
src/reconstruction/vmi/nchannel.jl
vmi_pipeline(; channels, basis, geom, to_backend = identity, kwargs...)
Measured channels to virtual monoenergetic images:
channels ─ hypr_lr ─┬─ merge_channels ─┬─ reduce_detector_rows ─┬─ decompose_nchannel ─┐
└──────────────────┘ └─ decompose_cong ─────┤
└─ tlbf_denoise ─┬─ reconstruct_basis_slice ─ apply_acnr_kalender! ────────────────┬─ synthesize_vmi_stack
└─ spectral_pair ─ acnr_complement! ─ image_hypr ─────────────────┘
Every stage's settings are keywords, so the function makes no decision the caller cannot see and override, which is what makes it usable as the inner call of an ablation sweep.
Required: channels (vector of K corrected log-transmission sinograms), basis (spectral_basis), geom (the acquisition geometry), matrix_size.
Stages, all optional:
denoiser = nothing or a SpectralHYPR — generalized HYPR-LR: its projection-domain instance (hypr_lr) on the counts of each detector row before anything else, with the dispersion measured from the acquisition's air rays by default, and its image-domain instance (image_hypr) on the reconstructed pair, after ACNR. Either can be nothing.
method = :nchannel or :cong; controls = NChannelControls().
merge_groups = nothing — channel groups summed before decomposing, e.g. [1:2, 3:4].
reduce_rows = false, rows = : — sum detector rows in counts, for a z-invariant object whose slices are all the same. Without it every row is kept and the volume is reconstructed slice by slice onto matrix_size.
use_tlbf = false, tlbf_alpha1, tlbf_alpha2, tlbf_radius — photon-counting only; filters each detector row in its own (column, view) plane.
use_acnr = true, acnr_passes = 4, acnr_beta_max = 20, acnr_hp_sigma_px = 1.5, acnr_window = 4. On a pair reconstructed by spectral_pair — a PairFilter, or the image-domain instance — ACNR acts on the complement only (acnr_complement!), before the image-domain instance; on a plain FDK pair it is basis-vmi's apply_acnr_kalender!.
matrix_size (required) — the reconstruction grid, (nx, ny, nz); this function never sees a workspace's ReconOptions, so the caller states the grid. One slice when the rows were reduced. fbp_filter = SoftFilter() — one window, or a PairFilter, a window for the composite and one for its complement; pair_basis = nothing, or (Estar = …, β = …) to fix the pair spectral_pair otherwise measures (a noise-free acquisition, reconstructed as its noisy counterpart); composite_energy = nothing, or the energy of the composite ACNR and the image-domain instance act on, fixed for a scanner and protocol (spectral_pair); antialias = true, recon_rows = geom.n_rows.
recon_method = :fbp (the published chain) or :hir, with hir_strength = 60, recon_projector = :dd_fast and hir_reference_kev = 70. :hir reconstructs the basis pair with the penalized iterative reconstructor and nothing else changes: T-LBF, ACNR and the synthesis run as before. Each material is reconstructed in μ-equivalent units at the reference energy, weighted by its own inverse variance from the K-channel Fisher information (uniformly for :cong, which has none), through a one-row geometry when the rows were reduced — see reconstruct_basis_slice for why each of those is needed. Run once each way for the analytic and iterative stacks of one acquisition.
vmi_energies = (40, 70, 100, 140).
keep_sinograms adds the decomposed pair; keep_diagnostics adds the estimator's per-ray maps (Fisher information, quality flags) for :nchannel.
The published photon-counting configuration of basis-vmi is vmi_pipeline(; channels, basis, geom, to_backend, reduce_rows = true, use_tlbf = true). The generic chain of basis-spectral-denoising, the same for photon-counting, rapid kVp-switching and dual-source acquisitions and with every detector row kept, is vmi_pipeline(; channels, basis, geom, to_backend, matrix_size, denoiser = SpectralHYPR(), fbp\_filter = PairFilter(composite, complement)).
Returns (vmis, energies, images = (water, iodine), quality, elapsed_s, settings), with vmis (nx, ny, nz, n_energies) in HU.
spectral_basis(; energies, response, I0, tolerance = 5e-5) -> NamedTuple
Assemble the estimator's spectral input from an ABSOLUTE per-energy response.
energies: length nE grid, keV.
response: (nc, nr, nE, K); nc or nr may be 1, meaning one response serves every ray.
I0: (nc, nr, K) absolute air counts, matching the ray axes of response.
tolerance: the largest relative disagreement allowed between sum(response; dims = 3) and I0. The likelihood needs absolute responses, not independently normalised spectra, so a disagreement is an error rather than something to renormalise away.
Returns (E, Φ, μρ_I, μρ_W, I0, μI_eff, μW_eff, normal_II, normal_IW, normal_WW, I0_relerr, n\_channels, ray\_resolved). The effective-energy quantities initialise the solver only.
spectral_basis(ws::PCCTWorkspace; I0 = ws.I0, tolerance = 5e-5)
Basis straight from a photon-counting workspace: the response the simulation applied (ws.W_matrix_gpu on ws.energies, through the per-ray source transmission ws.bowtie_spectral when the scanner has a bowtie) and the air response of every ray. Because this is the model that generated the data, the decomposition inverts exactly what was simulated and needs no calibration scan.
src/reconstruction/vmi/basis.jl
p_photoelectric(E_keV)
Photoelectric cross-section per electron at energy E_keV (Cong Eq 3c).
q_compton(E_keV)
Klein-Nishina Compton cross-section per electron at energy E_keV (Cong Eq 3d; f_kn from Eq 3e).
src/reconstruction/vmi/cong.jl
CongWorkspace{T, AŴ_L, AŴ_H, A1}
Pre-staged spectral basis on the sinogram's backend. Cong is per-pixel parallel and writes directly to the user's outputs, so the workspace has no sino-shape scratch — only the six basis tables.
Type parameters: • T — element type (Float32 enforced). • AŴ_L — type of staged ŵ_L (1D Vector{Float32} for centered spectrum, 3D [n_col,n_row,n_E] for per-ray bowtie). • AŴ_H — same for ŵ_H (must share dimensionality with AŴ_L). • A1 — type of staged p_L/q_L/p_H/q_H (always 1D Vector).
Construct with create_cong_workspace.
create_cong_workspace(sino_template, basis) -> CongWorkspace
Stage the Cong basis (ŵ_L, p_L, q_L, ŵ_H, p_H, q_H) onto the same backend as sino_template.
Cong is per-pixel parallel with no per-tile scratch, so this workspace's sole purpose is amortizing the basis-staging cost across repeated apply_cong! calls (e.g., hyperparameter sweeps, multiple scans sharing the same basis).
Arguments
sino_template : witness for the target backend.
basis : NamedTuple with ŵ_L, p_L, q_L, ŵ_H, p_H, q_H, all Float32 — built by your basis-construction cell (see notebook 06 / 07 for the canonical setup).
ŵ_L / ŵ_H may be 1D [n_E] (centered spectrum) or 3D [n_col, n_row, n_E] (per-ray bowtie). Both must share dimensionality.
apply_cong!(ws, sino_y, sino_c, sino_low, sino_high;
water_basis,
newton_max_iter = 12,
newton_tol = eps(Float32),
y_max_factor = 0.99,
y_max_cap = 1f7) -> (sino_y, sino_c)
Per-ray Cong 2022 decomposition, writing results into sino_y and sino_c. All four sino arrays share shape (nx, nv, nr) and Float32.
Arguments
ws : CongWorkspace from create_cong_workspace(...). REQUIRED — no internal allocation path.
sino_y, sino_c : output photoelectric and Compton line integrals, mutated in place.
sino_low, sino_high : measured low- / high-kVp sinograms (log line integrals).
Keyword arguments
water_basis : Float32 NamedTuple (a, c), water's coefficients in the (p, q) basis the workspace was built with — (a = 0f0, c = 1f0) for the material basis (p = iodine, q = water) that decompose_cong passes.
newton_max_iter : inner Newton iterations on the Eq 8 quintic.
newton_tol : |Δ| break tolerance for the Newton inner loop.
y_max_factor : safety factor on the outer Brent upper bound y_max = factor · p_L_meas / min(p_L).
y_max_cap : hard ceiling on y_max.
src/reconstruction/vmi/image_domain_decomp.jl
synth_vmi_2basis!(HU_E, c_water, c_iodine; energy_keV,
water_material = XA.Materials.water,
iodine_material = XA.Elements.Iodine)
-> HU_E
Per-energy VMI synthesis from a two-basis (water + iodine) image-domain decomposition, written into HU_E:
μ(E) = c_water · μρ_water(E) + c_iodine · 1e-3 · μρ_iodine(E)
HU(E) = 1000·(μ(E) − μρ_water(E)) / μρ_water(E)
= 1000·(c_water − 1) + c_iodine · α(E), α(E) = μρ_iodine(E) / μρ_water(E)
μρ are the mass attenuation coefficients of water_material and iodine_material at energy_keV. c_water is in g/mL (≈ 1 for water) and c_iodine in mg/mL. HU_E, c_water and c_iodine are 3-D Float32 arrays of the same size (an error is thrown otherwise); any array backend works.
synth_vmi_2basis(c_water, c_iodine; energy_keV,
water_material = XA.Materials.water,
iodine_material = XA.Elements.Iodine) -> AbstractArray{Float32, 3}
Allocating form of synth_vmi_2basis!: a new Float32 volume like c_water (same size and array backend) with the virtual monoenergetic HU at energy_keV, HU(E) = 1000·(c_water − 1) + c_iodine · μρ_iodine(E)/μρ_water(E), c_water in g/mL and c_iodine in mg/mL (3-D Float32 arrays of the same size).
src/reconstruction/vmi/mono_plus.jl
MonoPlusWorkspace{T, A3}
Pre-allocated scratch + output volumes for apply_mono_plus!.
Layout: • 3 vol-shape Float32 scratch: lp_buf, lp_opt, hp_opt. • out_vols: Vector{Array{Float32, 3}} of length n_energies, each shape matching the input volumes — reused across calls. • kernel_cache: lazy Dict{Float64, Matrix{Float64}} keyed by σ.
Construct with create_mono_plus_workspace.
create_mono_plus_workspace(vol_template; n_energies, mem_budget_GB = nothing)
-> MonoPlusWorkspace
Allocate the Mono+ workspace on the same backend as vol_template.
Arguments
vol_template : an AbstractArray{<:AbstractFloat, 3} whose shape + backend the workspace should match (one example VMI volume).
Keyword arguments
n_energies : number of target energies — sizes out_vols.
mem_budget_GB : explicit budget override. nothing = auto-probe via Sys.free_memory() · 0.6. Even though Mono+ runs on CPU (FFTW), the output volumes can be large enough on big inputs to OOM the host — the constructor will retry-with-shrink-on-failure.
Mono+ does not auto-tile; if even allocation at tile_size = 1 fails, the total working set (3 scratch + n_energies output volumes) just doesn't fit. Downsize the volume or reduce n_energies.
apply_mono_plus!(ws, volumes, energies;
E_noise_opt = 70.0,
σ_lp_px = 2.0,
verbose = true) -> NamedTuple
Compute Mono+(E) for each energy in energies using Grant 2014's frequency-split rule:
Mono+(E) = LP_σ(VMI_E) + VMI_opt − LP_σ(VMI_opt)
Mono+(E_opt) = VMI_opt (identity at the noise-optimal reference)
Mutates ws.out_vols in place and returns a NamedTuple with energies, volumes (= references into ws.out_vols), σ_lp_px, E_noise_opt. Copy if you need to retain across subsequent calls on the same ws.
Arguments
ws : MonoPlusWorkspace from create_mono_plus_workspace. REQUIRED — no internal allocation path.
volumes : Vector{Array{Float32, 3}} aligned to energies.
energies : Vector of target keV values.
Keyword arguments
E_noise_opt : noise-optimal reference energy (must appear in energies).
σ_lp_px : Gaussian LP σ in pixels. Either a scalar (same σ at every energy) or a per-energy vector (length == energies). σ = 0 is a valid identity case (Mono+(E) = VMI_E exactly, no FFT).
verbose : log timing + identity-sanity message.
phantom_mask: optional Bool array same shape as the volumes. When provided, voxels where the mask is false are reverted to the raw input VMI_E after the Mono+ formula. Eliminates the bright ring at the phantom-air boundary that LP-based frequency splitting injects via HP_opt = VMI_opt − LP_σ(VMI_opt). Pass an eroded phantom mask (≥ ~3·σ pixels of erosion — see BS.erode_mask_3d) so the residual edge spike falls outside the masked region. nothing (default) → no masking, identical to legacy behavior.
apply_mono_plus_regression!(ws, volumes, energies;
E_noise_opt = 70.0,
σ_lp_px = 1.5,
window = 4,
beta_max = 4.0,
verbose = true,
phantom_mask = nothing) -> NamedTuple
Structure-preserving Mono+. Same LP/HP skeleton as [apply_mono_plus!] (Grant 2014) — the quantitative low-frequency band is kept from the target energy E, so rod HU is untouched — but the hard high-frequency swap is replaced by a per-pixel local regression of the target HF onto the low-noise reference HF, the exact analog of apply_acnr_kalender!:
Mono+_res(E) = LP_σ(VMI_E) + β(x)·HP_σ(VMI_opt)
= VMI_E − ( hE − β(x)·hOpt )
β(x) = clamp( Σ_win hE·hOpt / Σ_win hOpt² , 0, β_max )
with hE = VMI_E − LP_σ(VMI_E), hOpt = VMI_opt − LP_σ(VMI_opt).
Why this reduces noise with zero resolution cost (by construction, not by an edge heuristic):
Real edge (present in both energies) → hE, hOpt strongly correlated → β ≈ the true local contrast ratio CE/Copt. β·hOpt rebuilds the target energy's TRUE edge amplitude (fixing plain Mono+'s contrast loss for fine iodine detail) but carried on the low-noise reference → denoised. Edge geometry = hOpt's geometry exactly (β only scales amplitude) → no blur, no misregistration.
Flat / noise-only region → the two energies' HF are independent noise realizations → uncorrelated → β ≈ 0 → HF ≈ 0 → HF noise removed (strictly better than plain Mono+, which injects hOpt's residual HF noise everywhere).
E == E_opt → hE ≡ hOpt → β ≡ 1 → output = VMI_opt identically (Grant's reference-identity property retained).
β_max is additionally variance-limited to the global HF-energy ratio λ = √(Σ hE² / Σ hOpt²) (Kalender's original bound), so a noisy local β̂ can never over-amplify the reference HF.
Keyword arguments
E_noise_opt : noise-optimal reference energy (must appear in energies).
σ_lp_px : Gaussian LP σ (px) for the HF split. Scalar or per-energy vector. σ = 0 ⇒ identity at that energy (Mono+res(E) = VMI\E).
window : half-width of the (2w+1)² regression window (px). Sizes only the ESTIMATE of β — the correction stays a per-pixel linear combination, so no signal smoothing occurs regardless of window.
beta_max : hard ceiling on the HF amplification β (further capped by the per-energy global HF ratio λ).
phantom_mask: optional eroded Bool mask; voxels where false revert to the raw VMI_E (kills the LP boundary ring — same semantics as apply_mono_plus!).
Mutates ws.out_vols in place; copy the returned volumes to retain them across subsequent calls on the same ws. Returns a NamedTuple with energies, volumes, σ_lp_px, E_noise_opt, and per-energy diagnostics β_mean / β_frac_active (fraction of pixels with β > 0).
apply_mono_plus_structure!(ws, volumes, energies;
E_noise_opt = 70.0,
σ_lp_px = 6.0,
window = 5,
beta_max = 6.0,
verbose = true,
phantom_mask = nothing) -> NamedTuple
The powerful resolution-preserving Mono+. Same Grant-2014 frequency-split skeleton and the same per-pixel contrast-matched transplant as [apply_mono_plus_regression!], but run in the regime that actually denoises:
Mono+(E) = LP_σ(VMI_E) + β(x)·HP_σ(VMI_opt), HP_σ = I − LP_σ
Two things make it much stronger than either prior variant:
WIDE transplant band (low cutoff / large σ). Post-FBP CT noise is high-frequency-weighted (ramp filter), so most of the noise variance lives above even a low cutoff. The generic Gaussian Mono+ and the narrow-σ regression Mono+ used a HIGH cutoff (σ≈1.5) → they only transplanted the top sliver of the spectrum → ~all the noise stayed in the LP band (impotent, esp. at low keV where σVMI is dominated by the broadband α(E)²σ²I term). Here σ is LARGE: the LP band keeps only the coarse quantitative base (per-region HU ⇒ rod accuracy preserved) and the entire broadband detail is transplanted from the low-noise 70 keV anchor. In a flat region this collapses the output to a heavily-smoothed base ⇒ σ plummets (≈ 15–20× for ramp-weighted noise).
Contrast-matched transplant (local β) — the anatomical-identity prior, quantified. β(x) = clamp(Σwin hE·hOpt / Σwin hOpt², 0, βmax) is the LOCAL contrast ratio: ≈0 in flat regions (⇒ pure smoothed base, noise killed), ≈ α(E)/α(Eopt) at real iodine/anatomy edges (⇒ the anchor's sharp low-noise edge, re-scaled to the TARGET energy's amplitude). This is why it does not blur: edges come from the full-bandwidth anchor at the target's own contrast. βmax is variance-limited to the global HF ratio λ = √(Σhᴇ²/Σhopt²) so a noisy β̂ cannot over-amplify.
Net: flat → smoothed base (huge low-keV denoise); edges → low-noise anchor structure at correct target contrast (full resolution). E == E_opt ⇒ β ≡ 1 ⇒ output = VMI_opt identically.
Only the default σ (6.0 vs 1.5) and window (5 vs 4) differ from apply_mono_plus_regression! in mechanism; σ_lp_px is the master dial — larger = wider transplant = more denoising (until the LP base gets too coarse to hold per-rod HU). All arguments as in apply_mono_plus_regression!.
src/reconstruction/vmi/nchannel.jl
NChannelControls(; kwargs...)
Numerical controls of the K-channel estimator. The defaults are the published ones.
iodine_bounds = (-0.10, 0.40) g/cm². Scaled to per-ray noise, not to geometry, so negative values are admitted: clipping at zero would bias the mean of a noisy zero upward.
water_bounds = (-2, 50) g/cm². The ceiling must exceed the longest water-equivalent path through the object; 50 leaves about 1.5× margin for a 33 cm body.
outer_iterations = 16, inner_iterations = 12, bisection_iterations = 28
max_iodine_step = 0.05, max_water_step = 5 (per Newton step, g/cm²)
parameter_tolerance = 5e-5 (relative), fisher_condition_limit = 1e8
air_gate = 0 — rays whose every |h_k| is below this are set to zero without solving. Zero disables the gate; it is off in every published configuration.
tile_views = 8 — views sent to the backend per kernel launch.
spectral_basis_from_bins(; energies, W_applied, I0, transmission = nothing, tolerance = 5e-5)
Basis for the K energy windows of one photon-counting acquisition. W_applied[e, k] is the detected response the forward model applied to the unattenuated central beam and I0[col, row, k] the air count of every ray and bin. With transmission[col, row, e] — the source transmission (bowtie × heel) the simulation applied per ray — the response is ray-resolved, W_applied[e, k]·transmission[col, row, e], which is exactly what each ray saw; without it one response serves every ray and I0 must then be constant across the fan.
spectral_basis_from_acquisitions(; acquisitions, tolerance = 5e-5)
Basis for K separate acquisitions: rapid kVp switching, or two tubes. Each entry of acquisitions is a NamedTuple with
energies: that acquisition's own grid, keV;
response: (nc, nr, nE) relative detected spectrum per ray;
I0_ray: (nc, nr) absolute air counts for that acquisition.
Grids are merged onto their sorted union, so 80 and 140 kVp coexist without resampling either. Each response is normalised over energy and scaled by its own I0_ray, which is where mA, duty cycle and rotation time enter: they change how many photons an acquisition delivers, not the shape of its spectrum.
reduce_detector_rows(; channels, basis, rows = :, floor_counts = 1e-12)
Combine detector rows into one by summing counts, carrying Φ and I0 through the same sum.
For an object invariant over the active longitudinal extent the rows are repeated measurements of one in-plane ray, so Y_k = Σ_r I0_{k,r} · exp(-h_{k,r}) is again Poisson with the summed rate and the likelihood is unchanged except that Φ_k and I0_k are the summed ones. Summing transmission and re-logging gives the same answer only when I0 does not vary across rows; it does whenever a bowtie is modelled, so the sum is taken in counts unconditionally.
Returns (channels, basis, n_rows, selected_rows).
merge_channels(; channels, basis, groups, floor_counts = 1e-12)
Sum channels into groups, in counts, carrying Φ and I0 with them: counts add, transmission does not, so merged Φ and I0 are summed over the group and never averaged. groups = [1:2, 3:4] collapses four photon-counting windows into a low and a high channel, which is how a two-measurement comparator is given the same acquisition.
The groups must partition the channels, each appearing exactly once: a channel in two groups would enter the likelihood twice as though it were an independent measurement, and one in none would be discarded silently.
Returns (channels, basis, groups).
prepare_channels(; channels, basis, merge_groups = nothing, reduce_rows = false, rows = :)
Everything the decomposition needs done to its inputs first: sum channel groups, then sum detector rows. Both reductions add counts and rebuild the response so that sum(Φ[:, k]) still equals that channel's I0. Returns (channels, basis, n_rows).
nchannel_profile_tile!(sino_I, sino_W, fisher_AA, fisher_AC, fisher_CC, quality_flag,
score_norm, outer_count, inner_count, hs::NTuple{K},
Φ, μρ_I, μρ_W, I0, μI_eff, μW_eff,
normal_II, normal_IW, normal_WW, controls::NChannelControls)
One tile of the K-channel profile-likelihood estimator; every array lives on one backend. hs holds the K corrected log-transmission tiles (n_col, n_row, n_view_tile).
Quality flag bits: 1 iodine at a bound, 2 water at a bound, 4 not converged, 8 ill-conditioned Fisher matrix or invalid initialiser, 16 the summed counts are not attainable anywhere inside the bounds, 32 a non-finite result.
Counts may be fractional after detector corrections, so nothing here assumes integers.
decompose_nchannel(; channels, basis, controls = NChannelControls(), to_backend = identity,
tile_views = controls.tile_views, keep_diagnostics = false)
The K-channel estimator over a full sinogram. channels is a vector of K corrected log-transmission arrays (n_col, n_row, n_view); basis comes from spectral_basis or one of its builders; to_backend uploads an array to the compute backend (CuArray, MtlArray, …; the default stays on the CPU).
Views are tiled so the backend never holds a whole sinogram of workspace, and every tile's buffers are released without a full garbage collection.
Returns (sino_iodine, sino_water, quality, elapsed_s, method = :nchannel, n_channels) with both sinograms in g/cm². keep_diagnostics = true adds the per-ray maps.
cong_material_basis(basis) -> NamedTuple
Translate a two-channel spectral_basis into the (ŵ, p, q) form apply_cong! wants. Φ carries absolute per-energy air counts; Cong wants the shape only, normalised so Σ_E ŵ = 1, because the absolute scale enters through the measured line integrals. p and q are the iodine and water mass attenuations: the same material-direct pair the K-channel estimator uses, so both invert the identical physical basis. A response with singleton ray axes collapses to a 1-D spectrum; a ray-resolved one stays 3-D, which is how a bowtie reaches the solver.
decompose_cong(; channels, basis, to_backend = identity, newton_max_iter = 12,
newton_tol = eps(Float32), y_max_factor = 0.99, y_max_cap = 1f7)
The Cong 2022 two-measurement estimator (apply_cong!) behind the same signature and return shape as decompose_nchannel, so one configuration can be sent through either. Exactly two channels are required; pass merge_groups to prepare_channels to sum a photon-counting acquisition into two.
angular_antialias_response(; n_views, matrix_nx) -> Vector{Float64}
Deterministic angular response: unit gain through mode ⌈π·matrix_nx/4⌉, then a raised-cosine roll-off across the oversampling margin only. It is all ones when the acquisition does not oversample the target grid, in which case the filter is an identity and is skipped. This is not a tunable denoiser: a 1200-view acquisition oversamples what a 512-pixel grid can represent, and the roll-off discards only the angular frequencies that grid cannot carry without aliasing.
reconstruct_basis_slice(sino, geom, matrix_size; to_backend = identity,
filter = SoftFilter(), n_rows = geom.n_rows, antialias = true,
method = :fbp, hir_strength = 60, projector = :dd_fast,
hir_weights = :uniform, scale = 1)
Reconstruct one sinogram, whether it carries a measured channel or an estimated basis material.
method = :fbp is FDK as basis-vmi does it: a single-row sinogram (what reduce_detector_rows leaves) is repeated to n_rows and backprojected onto the requested grid with filter.
method = :hir is the penalized iterative reconstructor at hir_strength (0-100, see get_hir_params). Three things make it sound on a basis sinogram, none of which FDK needs:
The weights. HIR's default data weighting is exp(-y), the Poisson heuristic for a log-transmission. A basis sinogram is g/cm², not a transmission — through a body it reads 20 to 30 for water and a fraction of one for iodine — so that weighting would switch the data term off inside the object for one material and leave it uniform for the other. hir_weights is what is passed to reconstruct! instead: :uniform, or the per-ray inverse-variance map from the decomposition's Fisher information (see vmi_pipeline).
The units. HIR's Huber threshold and regularisation strength were tuned on attenuation images in cm⁻¹. scale multiplies the sinogram before reconstruction and divides the image after, so with scale = μ/ρ of the material at a reference energy the reconstruction happens in μ-equivalent units and the strength dial means what its table says. FDK is linear and needs no such thing.
The geometry. A single-row sinogram is NOT repeated: it is reconstructed through a one-row geometry with a volume one row thick. Repeating one row to n_rows asserts that every row saw the same thing, which a one-slice forward model cannot reproduce at the outer rows; FDK ignores the contradiction, PWLS would fit it.
hir_weights and scale are ignored by :fbp.
synthesize_vmi_stack(water, iodine, energies) -> Array{Float32,4}
Virtual monoenergetic images in HU from a reconstructed basis pair: (nx, ny, nz, n_energies), every slice of the pair at every energy. water is a density in g/mL and iodine a mass density in g/cm³; the latter is converted to mg/mL for synth_vmi_2basis.
src/reconstruction/vmi/phantom_mask.jl
erode_mask_2d(mask2d::AbstractMatrix{Bool}; erode_px::Real) -> BitMatrix
"Soft" erosion of a 2D Bool mask by approximately erode_px voxels. Implementation: FFT-Gaussian-blur the float-cast mask with σ = erode_px pixels, then threshold at 0.999 — interior voxels (whose neighborhood is fully inside the mask) survive; edge voxels (whose neighborhood is partially outside) get rejected.
erode_px ≤ 0 returns a copy of the input mask (no-op).
erode_mask_3d(mask3d::AbstractArray{Bool, 3}; erode_px::Real,
broadcast_in_z::Bool = true) -> BitArray{3}
Per-slice 2D erosion broadcast across all z slices.
broadcast_in_z = true (default): erode the FIRST slice only and replicate (correct for z-invariant masks). Faster.
broadcast_in_z = false: erode each slice independently. Use when the mask varies in z.
src/reconstruction/vmi/roots_kernels.jl
brent_solve(f, a, b; xabstol=eps(T), xreltol=eps(T), maxiters=100)
Find a root of f in the bracket [a, b] using Brent's method.
Returns (root, converged::Bool). On failure (non-finite value, or maxiters hit without nextfloat(a) == b), converged = false and root is the best of (a, b) by |f|.
Identical algorithm to Roots.find_zero(f, (a, b), Roots.Brent()) for Float64 / Float32 inputs. Tolerances default to eps(T) (matches Roots.jl default_tolerances at convergence.jl:55-56).
Caller contract
f(a) and f(b) must have opposite signs (bracket must be valid).
f must be pure (no allocations, no I/O) for GPU execution.