The Five-Struct API

Build a phantom, configure a scanner, run a simulation, reconstruct an image.

BasisSimulator.jl exposes its entire user surface through five structs:

This notebook walks through each one — building a Gammex Model 472 phantom and scanning it on a model of the GE Revolution Apex Elite. We'll finish by comparing a standard-dose acquisition against a low-dose one.

Notebook Setup

Activate docs/Project.toml (which has BasisSimulator from the local source tree + WasmMakie), then bring in our imports — one per cell — and finally detect a GPU backend.

Backend-agnostic device transfer: to_gpu()

BasisSimulator dispatches kernels on array type via AcceleratedKernels.jl, so the same source runs on Metal, CUDA, ROCm, oneAPI, or CPU. The only thing the user controls is which array type they hand the simulator.

The cell below probes for Metal → CUDA → AMDGPU in that order, picks the first one that's installed and functional on the host, and falls back to plain Array (CPU) otherwise. Returned as a one-liner: to_gpu(arr).

Backend detected: MtlArray

The Five Structs

01. Phantom

A Phantom is three pieces of data:

1. Labeled mask — a 3D integer array. Each integer is a region label.

2. Materials dict — maps each label to an XrayAttenuation.Material.

3. Voxel size — (dx, dy, dz) in cm.

The factory create_gammex_472(...) builds the standard Gammex Model 472 multi-energy CT phantom: a 33 cm body of solid water with calcium and iodine inserts at known concentrations (used as a clinical CT calibration standard).

The mask is UInt8 of shape (512, 512, 16), and there are 27 materials in the dictionary.

Now we move the mask to the GPU once and rebuild the Phantom struct around it. Every subsequent forward projection runs on the device.

02. Scanner

Scanner describes the hardware — properties that don't change between acquisitions: source-to-isocenter distance, detector geometry, focal spot, bowtie filter, scintillator material.

The values below model a GE Revolution Apex Elite (256-row, 0.625 mm isotropic detector pixels, curved Lumex array), verified against vendor specs and clinical data.

03. CTProtocol

CTProtocol is the acquisition struct. We define two protocols — a standard-dose scan at 200 mA and a low-dose scan at 50 mA — both at 120 kVp. We'll reconstruct each and compare noise at the end.

04. SimOptions

SimOptions is the fidelity dial: which physics effects to include in the forward model.

The :eict preset enables a clinically-realistic stack — quantum noise, beam-hardening, scatter convolution, focal-spot blur, bowtie filtration, heel effect, and detector efficiency. To study any single effect's contribution, toggle it off below.

05. ReconOptions

ReconOptions controls only the reconstructed grid: matrix size, XY FOV, and optional Z extent. The reconstruction algorithm and apodization filter belong to the reconstruction workspace constructor, where they are consumed.

Simulate & Reconstruct

01. Forward Project (Standard Dose)

Time to actually scan the phantom. Each "phase" of GPU work — forward projection, reconstruction — lives inside its own let ... end block. The shape of every cell that touches the device is the same:

output_cpu = let
    # 1. Allocate GPU workspace + run the kernel
    ws = BS.create_…_workspace(…)
    BS.simulate!(ws, …)              # or BS.reconstruct!(…)

    # 2. Copy results off the device (CPU `Array`)
    result = Array(ws.sinogram)

    # 3. Explicit GPU cleanup
    ws = nothing
    GC.gc(true)

    result
end

02. Reconstruct (Standard Dose: FBP)

Same pattern: allocate an FDK workspace, run reconstruct!, copy the volume off the device, drop GPU references, GC. The output is in linear attenuation units (cm⁻¹), which we convert to Hounsfield Units using a polychromatic effective μ_water for our 120 kVp spectrum + Gammex 472 body (33 cm diameter chord through the center voxel).

BS.compute_polychromatic_μ_water resolves the bowtie-aware source spectrum (flat filter + protocol filters + bowtie all included), pre-hardens it through water_path_cm of solid water via Beer-Lambert (mimicking the hardening a center-of-phantom ray sees), then integrates μwater(E) against the hardened spectrum. The result matches what the FBP recon actually produces for solid water — so by construction the phantom-center voxel reads ≈ 0 HU, instead of ≈ −150 HU you'd get with a monoenergetic `getreferenceμwater(70.0)`.

Both protocols share the same kVp + filtration, so we compute the reference once and reuse it across §7 and §8.

03. Repeat (Low Dose)

Same two let ... end blocks, only protocol_lowdose swapped in. Everything else (scanner, simopts, reconopts, phantom) is reused unchanged — that's the point of splitting the API into separate structs. And because each previous block already cleaned up its GPU memory, this protocol gets a fresh device with no leftover state.

04. Postprocessing: The Full Correction Pipeline

A raw FBP recon (hu_std, hu_low above) is what comes out of the FDK kernel with no clinical corrections. Real CT vendors apply a stack of corrections on top — the same stack we'll build here:

a. Calibrate the BHC model

Both protocols use 120 kVp with the same filtration, so a single BHC model covers both. The calibration is a one-time spectrum fit — input is (energies, weights) from the resolved source spectrum, output is a TwoMaterialBHC polynomial model + the calibrated μ_water_ref (which becomes our reference for to_hounsfield).

Calibrated:

• ref energy = 69.7 keV

• lac water = 0.19321 cm⁻¹

b. Standard-Dose Corrected Pipeline

Same let ... end shape as §7, but with the current correction stack inline: detected-spectrum water BHC → FDK → HU, followed by non-mutating cupping QA. Noise is already present in the simulated detector counts. GPU buffers are dropped at the end exactly like before.

c. Low-Dose Corrected Pipeline

Identical to 9b — same detected-spectrum water BHC, recon, HU conversion, and QA — only the input sinogram (sim_low.sino) and dose level differ.

Results

Compare Protocols (Before & After Correction)

The two reconstructions use identical physics, geometry, and reconstruction settings — only the tube current differs. The lower-dose scan should show visibly more noise (≈√4 = 2× higher pixel σ) but the same mean HU values inside each rod. Postprocessing flattens the radial cup and adds the system counts-domain electronic noise — visible mostly in the rim region.

Noise (σ in HU, 60×60 ROI at center of solid water):

Raw σ ratio low/std = 2.26f0× (theory: 2.0× for 4× lower mA).

Scroll Through Slices

The cone beam means edge slices are reconstructed from obliquer rays than central ones. Slide through the corrected volumes and watch the ends — this is where usable-z limits and any residual cone behaviour show up first.

Verification

Quantitative pass/fail against physics-derived expectations — no eyeballing. Per-rod theory is mono HU at the BHC reference energy from the same XrayAttenuation data that drove the simulation; tolerances cover partial volume + noise + the single-kVp limit: the knobless water BHC maps the water component exactly, and the remaining high-Z (Ca/iodine) spectral residual is physics a single-energy scan cannot resolve — dual-energy/VMI territory.

✅ NB01 VERIFICATION: PASS (5/5)

Water per slice (corrected) — central-ROI falloff = bug; peripheral-only = real cone divergence:

• full body: Float32[1.7, 1.7, 1.6, 1.7, 1.7, 1.7, 1.7, 1.7]

• central r<5 cm: Float32[6.9, 6.5, 6.6, 6.5, 6.7, 6.7, 6.9, 6.7]

Per-rod, CENTRAL slice — theory = mono HU @ 69.7 keV; raw = uncorrected recon (alignment sentinel):

Summary

This notebook walked the entire BasisSimulator.jl user surface end to end:

• Five structs — Phantom, Scanner, CTProtocol, SimOptions, ReconOptions — each owning a clearly bounded slice of the simulation contract (object, hardware, acquisition, physics, output).

• Two let ... end blocks per scan — one for simulate!, one for reconstruct! — each ending with ws = nothing; GC.gc(true) so GPU buffers actually come back. The Pluto worker stays comfortably under memory pressure even with two protocols back-to-back.

• A clinical-grade postprocessing pipeline — calibrate_bhc_water → apply_bhc_water (knobless sinogram BHC, full detected spectrum) → FDK → to_hounsfield (with the BHC model's own μ_water_ref); quantum + DAS noise arrive in the counts domain from simulate!, and cupping/DC is measured, never applied (measure_radial_cupping).

• Backend-agnostic GPU dispatch — GPUSelect.Storage() resolves at startup, so the same notebook runs on Metal, CUDA, ROCm, oneAPI, or pure CPU with no notebook-level changes.

Every reconstruction in the rest of the docs (PCCT, Hybrid IR, dual-kVp VMI, …) reuses this same pattern — workspace in a let, kernel call, copy CPU result, nothing + GC.gc(true), return.