API Reference

The thing being scanned. A labeled voxel mask + a vector of XA.Material (one per label) + the physical voxel size in cm. The mask's array type drives the GPU/CPU choice for the rest of the pipeline.

phantom_cpu = BS.create_gammex_472(n_voxels = 512)
phantom = BS.Phantom(MtlArray(phantom_cpu.mask),
                     phantom_cpu.materials,
                     phantom_cpu.voxel_size)

see: notebooks 01, 02, 05

The hardware — geometry, source, detector, filtration. All distances in mm; detector pitch is at isocenter (not detector face — divide by magnification if you're porting from a face-pitch convention).

# GE Revolution Apex Elite (EICT, polychromatic)
scanner = BS.Scanner(
    source_to_isocenter = 625.6,
    source_to_detector  = 1100.0,
    detector_rows       = 256,
    detector_cols       = 834,
    detector_row_size   = 0.625,
    detector_col_size   = 0.6,
    detector_shape      = BS.CURVED_DETECTOR,
    bowtie_filter       = :ge_revolution_large,
    detector_material   = :lumex,
    detector_depth      = 3.0,
)

see: notebooks 01, 02, 03, 04, 05, 06

The acquisition — what tube voltage, how much current, how many views, how thick a slab. additional_filters is a list of extra material/thickness pairs layered on top of the scanner's flat filter.

protocol = BS.CTProtocol(
    kVp = 120,
    mA  = 200.0,
    views          = 500,
    rotation_time  = 1.0,
    collimation_mm = 5.0,
    additional_filters = [("Al", 4.5), ("Ti", 0.9)],   # extra Al + Ti on top of Scanner.flat_filter_*
)

Active detector rows are derived from collimation_mm via geom.n_rows = round(collimation_mm / scanner.detector_row_size).

fidelity is the master switch — :eict for energy-integrating CT (single-kVp polychromatic), :pcct for photon-counting CT (multi-bin spectral). Each fidelity preset turns on the appropriate physics; the per-effect use_* toggles override the preset for A/B testing.

sim_opts = BS.SimOptions(fidelity = :eict, seed = 1234)
# or:
sim_opts = BS.SimOptions(fidelity = :pcct, seed = 1234, pcct_noise_reduction = 0.3)

The output grid + algorithm. Recon is always centered at isocenter — there's no off-center FOV parameter (see notebook 05 for the SFOV-equivalent phantom-cropping pattern).

recon_opts = BS.ReconOptions(
    algorithm   = :fdk,
    matrix_size = (512, 512, 8),
    fov_cm      = 35.0,
    z_cm        = 0.5,
    filter      = :standard,
)

The Gammex Model 472 calibration phantom — 33 cm solid-water cylinder + 7 calcium rods (50–600 mg/mL) + 7 iodine rods (2–20 mg/mL) at 28 mm rod diameter. Returns a Phantom with semantic mask labels (10–16 = Ca, 20–26 = I).

phantom = BS.create_gammex_472(n_voxels = 256, n_slices = 8, fov_cm = 35.0, z_cm = 0.5)

Wrap any labeled UInt8 voxel grid + a label → material dict into a Phantom. Used by notebooks 02 / 05 for XCAT loading.

xcat_mask = load_xcat_bin("vmale_50.bin"; cols = 1600, rows = 1400, slices = 500)
phantom = BS.create_phantom_from_mask(xcat_mask, materials_dict, (0.04, 0.04, 0.04))

Per-voxel linear attenuation coefficient (cm⁻¹) at a given energy, derived by indexing each voxel's material's μ at energy_keV. Useful for monochromatic forward projection or theoretical reference volumes.

Linear attenuation μ (cm⁻¹) of a single material at a single energy. Backed by NIST XCOM data through XrayAttenuation.jl.

Mass attenuation μ/ρ (cm²/g). Numerically equal to compute_μ_at_energy for water (ρ ≈ 1). Used internally by the M-matrix μ_water computation.

Convenience: μ_water at a single energy. For polychromatic scans, prefer compute_polychromatic_μ_water (§5) which integrates over the spectrum + phantom hardening.

XrayAttenuation.jl is re-exported as BS.XA — you don't need to add XrayAttenuation as a separate dependency. Use it for the canonical Material constructor and the prebuilt NCAT/XCAT tissue library.

# Prebuilt tissues (full NCAT/XCAT library)
water  = BS.XA.Materials.water
muscle = BS.XA.Materials.ncat_muscle
blood  = BS.XA.Materials.ncat_blood

# Custom material from elemental mass fractions
iodine_blood = BS.XA.Material(
    "iodine_blood_5mgml",
    z_a_ratio,
    mean_excitation_energy * BS.XA.u"eV",
    1.06 * BS.XA.u"g/cm^3",
    Dict(1 => 0.105, 6 => 0.110, 7 => 0.033, 8 => 0.745, ...),  # Z → mass_frac
)

see notebook 02 for a full custom-material walkthrough

Pre-computes 3-D source/detector positions for every view angle. Built automatically inside create_*_workspace calls; you only need to construct one explicitly when you want to inspect the recon grid before running simulate! (e.g. the affine round-trip in notebook 05).

geom = BS.CTGeometry(scanner;
    n_angles       = protocol.views,
    fov_cm         = recon_opts.fov_cm,
    z_cm           = recon_opts.z_cm,
    collimation_mm = protocol.collimation_mm,
)
# geom.n_cols, geom.n_rows, geom.n_angles, geom.fov, geom.angles, ...

Detector geometry constants for Scanner.detector_shape. Curved detector ≈ all clinical CT (every column subtends the same angle from source); flat detector for cone-beam micro-CT and most academic phantom studies.

Maps 0-indexed phantom voxel (i, j, k) → world (x, y, z) in cm. Pure scale + translate; encodes phantom.voxel_size and phantom.origin.

Same shape as above but for the reconstruction grid. Origin is hard-locked to isocenter-centered (tx = -fov/2 + voxel_size/2).

Resample the phantom mask onto the recon grid using the two affines above. :nearest preserves integer labels (returns UInt8); :linear does trilinear (returns Float32). Use the linear path on a binary mask if you want true partial-volume fractions.

gt_nn  = BS.resample_to_recon(phantom_cpu, geom, recon_opts.matrix_size; method = :nearest)
gt_lin = BS.resample_to_recon(phantom_cpu, geom, recon_opts.matrix_size; method = :linear)
# Both have shape == recon volume shape — voxel-aligned for ROI extraction.

see notebook 05 for the full round-trip + custom-interpolator pattern

Raw tube spectrum at the given kVp before any flat / additional filtration. Backed by IPEM Report 78 tables. energies in keV, weights are normalized photon-count fractions per energy bin.

The post-filtration spectrum the simulator actually uses, EXCLUDING bowtie. This is what you want for BHC calibration + scatter-bin-weight derivation (both are bowtie-agnostic by construction).

Same as above but with the bowtie attenuation folded in per detector column. ŵ is 3D (n_col, n_row, n_E) when bowtie is present, 1D otherwise.

Spectrum-weighted μ_water (cm⁻¹) after Beer-Lambert hardening through water_path_cm of water. Pass the phantom diameter (full chord, not radius) — e.g. 33.0 for a Gammex 472 body. This is the value to use as to_hounsfield(...; μ_water = ...) for clean HU baselines.

μ_water = BS.compute_polychromatic_μ_water(
    sim_opts, protocol;
    scanner = scanner, geom = geom, water_path_cm = 33.0,
)

EICT (single-kVp polychromatic) workspace. Pre-allocates per-energy intensity buffers, bowtie spectral table, scatter / focal-spot / lag kernels, noise staging, and the geometry struct. Backend is inferred from phantom.mask.

PCCT (multi-bin spectral) workspace. Adds per-energy DRM tile buffers, native-dexel-resolution buffers (when binning_factor > 1), pileup matrix, charge-sharing kernel, anti-coincidence matrix. Significantly larger than the EICT workspace; sized to your protocol's view count.

ws = sim_opts.fidelity == :pcct ?
BS.create_workspace(scanner, protocol, sim_opts, recon_opts, phantom) :
BS.create_eict_workspace(scanner, protocol, sim_opts, recon_opts, phantom)

The forward-projection entry point. Dispatches on workspace type — EICTWorkspace runs the polychromatic Siddon → noise → optional corrections path; PCCTWorkspace runs spectral ray-tracing with detector-response convolution + per-bin Poisson noise.

ws = BS.create_eict_workspace(scanner, protocol, sim_opts, recon_opts, phantom)
BS.simulate!(ws, phantom, scanner, protocol, sim_opts, recon_opts)
sino = Array(ws.sino_noisy_out)   # pull off GPU

see notebooks 01–06 for the canonical call patterns

Pre-allocates the FBP-domain filter kernel + back-projection buffers. filter takes a Symbol from the preset table below or a CustomFilter.

Runs filtered back-projection in-place. Returns the μ-domain volume; convert to HU via to_hounsfield (§8).

ws_fdk  = BS.create_fdk_recon_workspace(sino_gpu, geom, matrix_size; filter = :standard)
recon_μ = Array(BS.reconstruct!(ws_fdk, sino_gpu, geom, matrix_size))
recon_HU = BS.to_hounsfield(recon_μ; μ_water = μ_water_120)

User-defined frequency-domain apodization. Two same-length tuples of (normalized frequency, gain). Linear interpolation between control points. Use it when the named presets don't quite match your scanner's vendor filter.

fdk_filter = BS.CustomFilter(
    (0.0, 0.25, 0.5,  0.75, 1.0),
    (1.0, 0.75, 0.6,  0.2,  0.001),
)
ws = BS.create_fdk_recon_workspace(sino, geom, matrix_size; filter = fdk_filter)

Zeros out voxels outside the recon's inscribed circular FOV (replaces them with sentinel_μ — typically −1024 HU = air). Run after FBP to clean up corner ringing.

Penalized weighted least-squares iterative recon (Huber prior + ordered subsets). strength ∈ {0, 1, 2, 3, 4, 5} — 0 = pure FDK pass-through, 5 = heaviest regularization. Workspace warm-starts from FDK, then runs the iterations.

Same dispatch shape as FDK. Converges in 2–5 iterations on clinical data; slowest reconstruction option — use FDK + image-domain BHC if you need speed.

ws_hir  = BS.create_hir_recon_workspace(sino, geom, matrix_size; strength = 3, filter = :standard)
recon_μ = Array(BS.reconstruct!(ws_hir, sino, geom, matrix_size))

see notebook 02 for FDK vs HIR side-by-side on XCAT

Joint SVD + bilateral denoising on a μ-domain volume pair (Clark & Badea 2023). Run before HU conversion — that's where the iodine/water noise is maximally anti-correlated. h_param controls the bilateral filter strength (lower = less smoothing).

Least-squares fit of the Ding (2012) image-domain decomposition coefficients (a₀, a₁, a₂) from a calibration table of measured rod HUs and known iodine concentrations. Returns the coefficients plus per-bin α sensitivities (HU per mg/mL) and the RMS calibration error.

Per-voxel c_iodine = a₀ + a₁·HU_low + a₂·HU_high (mg/mL). Pair with a 3-slice axial median filter (apply_median_z) to wipe single-voxel speckle without in-plane blur.

Synthesize a virtual monoenergetic image at any keV from the low-energy HU + iodine concentration map. Implements HU_E = HU_low + c_iodine·(α_E_phys − α_low_cal) — per-energy α from XrayAttenuation physics, low-bin α from the Ding fit.

Pre-allocates buffers for the Mono+ post-processing — a per-energy frequency-domain filter that anchors at the noise-quietest VMI energy and applies low-pass shaping to the others (FBP-equivalent noise behavior across the full keV sweep).

In-place Mono+ post-processing on a vector of synthesized VMIs. E_noise_opt is the anchor energy (passed through unchanged); σ_lp_px is a per-energy vector of low-pass σ in pixels.

vols_in = [vmi_dict[E] for E in energies]
σ_vec = Float64[2.0, 0.0, 2.0, 2.0]   # paired with energies
ws  = BS.create_mono_plus_workspace(vols_in[1]; n_energies = length(energies))
res = BS.apply_mono_plus!(ws, vols_in, energies; E_noise_opt = 70.0, σ_lp_px = σ_vec)

see notebooks 03 (dual-kVp) and 04 (PCCT) for full VMI pipelines

Axial-only median filter with the given z-radius (radius=1 → 3-slice window). Wipes single-voxel xy impulses without any in-plane blur — perfect for cleaning up speckle in iodine-concentration maps.

Per-voxel 1000·(μ − μ_water)/μ_water. Pass the polychromatic μ_water from §5 — using a single-energy μ_water with a polychromatic recon will land water HU off by tens of HU.

Scalar conversions for one-off calculations (e.g. theoretical HU at a given energy from a material's μ). No allocation, no array machinery.

In-place additive Gaussian — models the dose-independent DAS readout floor (≈ 28 HU on most clinical scanners). Run AFTER to_hounsfield so it lands at the right scale.

Even-polynomial radial cup correction on water-like voxels. Mostly cosmetic after BHC has done its job; useful when you want a flat HU profile across the phantom for visual inspection.

One-time fit: feed it the resolved spectrum and it returns a TwoMaterialBHC polynomial model + the calibrated polychromatic μ_water_ref. Use that μ_water_ref for the FBP-side to_hounsfield call so HU baselines stay self-consistent.

Sinogram-domain water-stage correction using the calibrated polynomial. Run BEFORE FDK; removes the bulk polychromatic bias from the line integrals.

Image-domain bone-region refinement. Smoothsteps between hu_low and hu_high to identify bone voxels, scales their residual error, and subtracts. Run AFTER FDK, BEFORE to_hounsfield.

# Calibrate once
e, w = BS.resolve_source_spectrum_without_bowtie(sim_opts, protocol; scanner = scanner)
bhc_model, μ_water_ref = BS.calibrate_bhc_two_material(e, w; order = 2)

# Apply at simulate time
BS.apply_bhc_two_material(sino, bhc_model, geom, recon_opts.matrix_size)
recon_μ = Array(BS.reconstruct!(ws_fdk, sino, geom, recon_opts.matrix_size))
BS.apply_bhc_image_domain(recon_μ, geom, recon_opts.matrix_size, μ_water_ref;
    hu_low = 50, hu_high = 150, scale_factor = 0.2)
recon_HU = BS.to_hounsfield(recon_μ; μ_water = μ_water_ref)

see notebook 02 for the full BHC pipeline on XCAT

Two efficiency modes — analytic Beer-Lambert through the scintillator depth, or MC-derived lookup tables for the GE GEMSTONE / standard GOS / CsI / CdTe scintillators. Sim-time selection via sim_opts.detector_efficiency_mode.

Pre-computed MC efficiency per (energy, depth) for the GE GEMSTONE detector. Used by MC_LUT mode.

Detective quantum efficiency at the spectrum's mean — 0 = perfect, 1 = no detection. Useful for theoretical noise-budget calculations against vendor spec sheets.

Per-energy quantum efficiency ∈ [0, 1] for a given scintillator material + depth. EICT analog of quantum_efficiency_vector (PCCT).

Construct the PCCT detector struct directly from the scanner's PCCT fields (energy_thresholds, energy_resolution, charge_sharing_fwhm, dead_time_ns). Internal helper exposed because notebook 04's μ_water + scatter-fraction derivations need it.

Monte-Carlo-derived detector response matrix. R[E, b] = probability that a photon of energy E is recorded in bin b — captures charge sharing, fluorescence escape, anti-coincidence.

Per-energy detection efficiency for the CdTe / CZT / Si crystals — Beer-Lambert through the active depth, includes K-edge and Fano weighting.

Per-bin Poisson noise applied AFTER scatter injection. Uses the detector DRM to compute the right per-bin I0 from the source spectrum.

Spectral migration matrix S[b_in, b_out] from MC simulation of dead-time pileup at the detector's per-dexel count rate. Computed once at workspace creation; used during PCCT noise application to migrate counts between bins.

Closed-form count-loss factor for the semi-non-paralyzable dead-time model. Matches the MC pileup result to within ~1% at aτ < 5. aτ is per-dexel count rate × dead-time.

Equivalent body diameter of the scanned object — drives the scatter kernel amplitude (Ohnesorge size scaling, exponent 1.5). Pass phantom.voxel_size .* 10.0 for the mm units.

Builds the Ohnesorge spatial scatter convolution kernel with geometry scaling — accounts for SID/SDD, air gap, pixel pitch, and phantom size. Pass to estimate_scatter_field!.

In-place: predicts the scatter sinogram from the recombined primary line integral. Output is the scatter field (counts), to subtract from the measurement after scaling by per-bin scatter fractions.

Per-energy scatter-to-primary spectral ratio. Used by compute_scatter_bin_weights to derive PCCT bin-specific scatter fractions.

Per-bin fraction of total scatter that falls into each PCCT energy bin (Σ_b frac_b = 1). Subtract scatter_field × I0_total × frac_b from each bin's measured counts to do exact scatter correction.

see notebook 04 §6 for the canonical PCCT scatter-correction flow

Pre-computed Gammex 472 rod HU table from a real GE Revolution Apex Elite dual-kVp Gemstone Spectral Imaging (GSI) acquisition (post-RSKR, post-cupping). Feed into fit_ding_coeffs for self-calibrated decomposition.

cal = BS.GE_REVOLUTION_APEX_ELITE_DE_CAL
# cal["I 5.0 mg/mL"] = (material = :iodine, mg_per_mL = 5.0, HU_low = ..., HU_high = ...)

Same shape but for the Siemens Naeotom Alpha PCCT scanner at 140 kVp / 174 mA / 10 mGy CTDIvol. Used in notebook 04 for self-calibrated Ding decomposition.

Same shape, 120 kVp variant of the Naeotom Alpha calibration.

Vendor-style GE Revolution Apex Elite reconstruction filter kernels indexed by vendor preset name (:standard, :soft, :bone). Bypass the Symbol path for ReconOptions.filter and pass one of these directly when you need GE-specific apodization.

ws = BS.create_fdk_recon_workspace(sino, geom, matrix_size;
filter = BS.GE_APEX_ELITE_FILTERS[:standard])

UInt8 enum constants used by create_gammex_472 to label voxels in the mask. Calcium inserts: REGION_CA_50 through REGION_CA_600 (10 → 16); iodine inserts: REGION_I_2_0 through REGION_I_20_0 (20 → 26). Match the REGION_TO_MATERIAL dict for material lookup.

# Filter for all calcium voxels
ca_mask = (UInt8(BS.REGION_CA_50)  .≤ phantom.mask) .&
          (phantom.mask .≤ UInt8(BS.REGION_CA_600))