The Five-Struct API
Describe a scan with five structs, simulate it, reconstruct it, and read off its dose.
Every simulation in BasisSimulator.jl is specified by five plain Julia structs, one for each thing a CT scan needs:
| Struct | What it describes | Built here with |
|---|---|---|
Phantom | the object: a labeled voxel mask, one XrayAttenuation.Material per label, the voxel size | create_gammex_472 |
EICTScanner / PCCTScanner | the hardware: geometry, detector, focal spot, filtration, bowtie | EICTScanner(; ...) |
CTProtocol | the acquisition: kVp, mA, views, rotation time, collimation, pitch | CTProtocol(; ...) |
SimOptions | the physics: which effects the forward model includes, the noise seed, the projector | SimOptions(; ...) |
ReconOptions | the output grid: matrix size, in-plane FOV, z extent | ReconOptions(; ...) |
Three calls turn them into an image:
Phantom ─┐
Scanner ─┤
CTProtocol ─┼─▶ create_workspace ─▶ simulate! ─▶ reconstruct! ─▶ to_hounsfield ─▶ HU
SimOptions ─┤ │
ReconOptions ─┘ └─▶ result.dose: CTDIvol, DLPThis notebook scans a Gammex 472 calibration phantom on a model of the GE Revolution Apex Elite twice, at 200 mA and at 50 mA. It reports the dose of each scan, reconstructs both with filtered back-projection (FBP) and with Hybrid IR, and checks the result against physics-derived expectations: water at 0 HU, noise scaling with dose, and every calcium and iodine rod against its theoretical HU.
Notebook setup
The notebook activates the shared docs/ environment (which pulls BasisSimulator from this source tree) and loads CairoMakie for the figures.
import BasisSimulator as BSChoosing the device
BasisSimulator writes its kernels once, with AcceleratedKernels.jl, and runs them on whatever array type it is handed. GPUSelect.Storage() returns the array type of the GPU it finds (CuArray, MtlArray, ROCArray or oneArray), or Array when there is none. Moving the phantom mask to that type is the only thing that decides where the whole pipeline runs.
begin
import GPUSelect
AT = GPUSelect.Storage() # CuArray / MtlArray / ROCArray / oneArray, or Array on a CPU-only host
to_gpu(x) = AT(x)
GPU_BACKEND = (name = string(nameof(AT)),)
end;Backend detected: CuArray
The five structs
1. Phantom
A Phantom holds three things:
a labeled mask: a 3D array of unsigned integers, one region label per voxel;
a material per label: an
XrayAttenuation.Material, which gives the simulator μ(E) at every energy of the spectrum;the voxel size
(dx, dy, dz)in cm. The phantom is centred on the isocentre unless you pass an origin.
create_gammex_472 builds the Gammex Model 472 multi-energy phantom: a 33 cm solid-water body with seven calcium rods (50–600 mg/mL) on an inner ring and seven iodine rods (2–20 mg/mL) on an outer ring, each 28 mm across.
phantom_cpu = BS.create_gammex_472(
n_voxels = 512, # in-plane voxels across the 35 cm field
n_slices = 16, # axial slices
fov_cm = 35.0,
z_cm = 1.0,
);The mask is UInt8 of size 512 × 512 × 16, with voxels of 0.684 × 0.684 × 0.625 mm and 27 entries in its material table (entry L + 1 is the material of label L; labels the phantom doesn't use are air).
Phantom is parameterised on its mask's array type, so moving the mask to the device once is enough: every forward projection of this phantom then runs there.
phantom = BS.Phantom(
to_gpu(phantom_cpu.mask),
phantom_cpu.materials,
phantom_cpu.voxel_size,
phantom_cpu.origin,
phantom_cpu.extent,
);
2. The scanner: EICTScanner
The scanner is the hardware: everything that stays the same from one scan to the next. There are two families, one per detector type, and they share a ScannerGeometry (source and detector distances, detector array, focal spot, gantry, flat filter, bowtie):
EICTScanner: an energy-integrating scintillator detector (material, depth, fill factors, detection gain, electronic noise);PCCTScanner: a photon-counting detector (sensor, energy thresholds, energy resolution, charge sharing, dead time and pile-up, pixel binning).
Each constructor rejects the other family's keywords, so a photon-counting parameter can't end up on a scintillator by accident. Geometry fields read the same on both: scanner.detector_rows.
The values below model the GE Revolution Apex Elite: 256 rows of 0.625 mm, 834 columns, a curved (arc) Lumex detector and GE's large-body bowtie. Distances are in mm; the detector pitch is given at the isocentre.
scanner = BS.EICTScanner(
# Geometry (mm)
source_to_isocenter = 625.6,
source_to_detector = 1100.0,
# Detector array (pitch at the isocentre, mm)
detector_rows = 256,
detector_cols = 834,
detector_row_size = 0.625,
detector_col_size = 0.6,
# Focal spot
focal_spot_width = 1.0,
focal_spot_length = 1.0,
target_angle = 10.0,
# Filtration
flat_filter_material = :aluminum,
flat_filter_thickness = 2.5,
bowtie_filter = :ge_revolution_large,
# Scintillator
detector_material = :lumex,
detector_depth = 3.0,
fill_factor_row = 0.9,
fill_factor_col = 0.9,
# Data-acquisition electronics: added to the counts before the log, as a real DAS does
electronic_noise = 3500.0, # e⁻ rms
detection_gain = 10.0, # e⁻ / keV
);Scanner: 256 rows × 834 columns, detector shape :arc, bowtie :ge_revolution_large, source–isocentre 625.6 mm.
Hardware here, acquisition in the protocol
Anything that changes from scan to scan (kVp, mA, views, rotation time, collimation, pitch) belongs to CTProtocol. That split is what lets one scanner run a sweep of protocols, as the two scans below do.
3. CTProtocol
The acquisition. Two protocols share everything except the tube current: a standard-dose scan at 200 mA and a low-dose scan at 50 mA, both 120 kVp, 500 views in a 1 s rotation, 5 mm of collimation and 4.5 mm of extra aluminium filtration.
protocol_standard = BS.CTProtocol(
kVp = 120,
mA = 200.0,
views = 500,
rotation_time = 1.0,
collimation_mm = 5.0, # nominal beam width N·T
additional_filters = [("Al", 4.5)], # on top of the scanner's 2.5 mm flat filter
);protocol_lowdose = BS.CTProtocol(
kVp = 120,
mA = 50.0, # a quarter of the dose
views = 500,
rotation_time = 1.0,
collimation_mm = 5.0,
additional_filters = [("Al", 4.5)],
);Collimation selects the detector rows
The detector is 256 × 0.625 mm = 160 mm wide; collimation_mm = 5.0 uses only the central rows. The simulator adds symmetric guard rows automatically where the cone beam needs them to cover the whole reconstruction cylinder at the edge of the field. Leave collimation_mm unset to use the full detector; add pitch and n_rotations for a helical scan (notebook 11).
4. SimOptions
The physics of the forward model. Every effect is on by default except optical crosstalk: fill factor, energy-dependent detector efficiency, scatter, focal-spot blur, quantum and electronic noise, detector lag and the anode heel effect. The polychromatic source spectrum, the bowtie and beam hardening are always part of the model. Turn an effect off to measure its contribution.
projector = :dd_fast (the default) is the anti-aliased distance-driven projector that walks the volume once for the whole spectrum; :siddon is the point-sampled ray tracer, kept for comparison.
The detector integrates while the gantry turns. A clinical detector reads each view for the whole view period, during which the gantry sweeps the view spacing Δθ = 360° / views. Each reading is therefore the transmitted intensity averaged over that arc: the object is blurred along the direction of rotation, more the farther it sits from the isocentre, while the noise, counted once per view, stays independent from view to view. view_samples samples that average at sub-views spread across the arc (the midpoint rule, in the intensity domain, before scatter, noise and the detector effects); view_arc is the fraction of Δθ the detector integrates over, 1.0 for a detector that reads for the whole view period and the duty cycle of one energy for rapid kVp switching. The default, view_samples = 1, is an instantaneous point view. Every clinical scanner in these notebooks is modelled with view_samples = 5, which costs five forward projections per view.
sim_opts = BS.SimOptions(
seed = 1234, # fixed noise realisation: same seed + same inputs = same sinogram
projector = :dd_fast,
view_samples = 5, # the detector integrates over the arc of each view (5 sub-views)
# Switch individual effects off, e.g.:
# use_scatter = false,
# use_heel_effect = false,
# use_focal_spot = false,
# use_lag = false,
# use_noise = false, # the noiseless expected sinogram
);Detector-specific physics lives on the scanner
SimOptions carries only the effects both detector families share. Pile-up, its correction, photon-counting scatter correction and the count-noise blend are fields of PCCTScanner; detector lag applies to the scintillator only.
5. ReconOptions
The output grid: matrix size, in-plane field of view and z extent, centred on the isocentre. Here the 5 mm beam is reconstructed as eight 0.625 mm slices on a 512 × 512 grid over 35 cm (0.68 mm pixels). The algorithm and its filter are chosen where they are used, when the reconstruction workspace is created.
recon_opts = let
slice_thickness_mm = 0.625
n_slices = round(Int, protocol_standard.collimation_mm / slice_thickness_mm)
BS.ReconOptions(
matrix_size = (512, 512, n_slices),
fov_cm = 35.0,
z_cm = protocol_standard.collimation_mm / 10,
)
end;Simulate
Workspace, then simulate!
create_workspace does all the set-up once: it builds the scan geometry, resolves the source spectrum through the filters and bowtie, and allocates every device buffer on the phantom's backend. simulate! then runs the forward model into ws.sinogram, the log-transformed line integrals -log(I/I₀) of size (columns, rows, views), noise included. The same constructor serves both detector families: for a PCCTScanner it builds the photon-counting workspace, and simulate! returns one sinogram per energy bin.
simulate! returns the dose of the acquisition it simulated. The CTDI is computed the way IEC 60601-2-44 defines it: the simulated beam is transported by Monte Carlo through the 32 cm PMMA body phantom, and the air kerma is integrated over a 100 mm pencil chamber at the centre and at the periphery. It therefore follows this scanner's spectrum, filtration and bowtie rather than a lookup table.
Each scan runs inside a let block that keeps only CPU copies of what it needs, so the device buffers are released before the next scan starts.
sim_std = let
ws = BS.create_workspace(scanner, protocol_standard, sim_opts, recon_opts, phantom)
result = BS.simulate!(ws, phantom, protocol_standard, sim_opts)
out = (sino = Array(ws.sinogram), geom = ws.geom, dose = result.dose)
ws = nothing
GC.gc(true) # hand the device memory back now
out
end;sim_std.dose
The low-dose scan
The same two calls with protocol_lowdose. The phantom, scanner, physics and grid are reused unchanged, which is the point of keeping them in separate structs.
sim_low = let
ws = BS.create_workspace(scanner, protocol_lowdose, sim_opts, recon_opts, phantom)
result = BS.simulate!(ws, phantom, protocol_lowdose, sim_opts)
out = (sino = Array(ws.sinogram), geom = ws.geom, dose = result.dose)
ws = nothing
GC.gc(true)
out
end;sim_low.dose| protocol | mAs / rotation | CTDIvol (mGy) | DLP (mGy·cm) | scan length (cm) |
|---|---|---|---|---|
| standard dose | 200.0 | 13.72 | 6.86 | 0.5 |
| low dose | 50.0 | 3.43 | 1.72 | 0.5 |
The beam is identical, so CTDIvol scales with the tube current alone: the low-dose scan delivers 0.25 of the standard dose. CTDIw of this beam is 6.86 mGy per 100 mAs. The Monte Carlo's statistical uncertainty is ±0.48 %. Pass dose_kwargs = (; phantom = :head16) to simulate! for the 16 cm head phantom.
Repeated noise draws: keep_projection and projection
Everything before the noise (the view-integrated forward projection and the physics that precedes the noise, scatter included) is the same for every noise draw of one phantom, protocol and workspace. simulate!(...; keep_projection = true) returns it as result.projection, and simulate!(...; projection) draws a new realization from it, with the seed of the SimOptions it is given, without projecting again. A noise study over many seeds, or a noise-free reference (use_noise = false) that shares the noisy scan's projection, pays for the projection once.
noise_draws = let
ws = BS.create_workspace(scanner, protocol_standard, sim_opts, recon_opts, phantom)
# one full simulation that keeps everything before the noise …
t_full = @elapsed (kept = BS.simulate!(ws, phantom, protocol_standard, sim_opts;
report_dose = false, keep_projection = true))
# … and a second noise realization drawn from it with another seed
opts_2 = BS.SimOptions(; seed = 5678, projector = sim_opts.projector,
view_samples = sim_opts.view_samples)
BS.simulate!(ws, phantom, protocol_standard, opts_2; report_dose = false, projection = kept.projection)
sino_reused = Array(ws.sinogram)
# the same seed simulated from scratch, for comparison
BS.simulate!(ws, phantom, protocol_standard, opts_2; report_dose = false)
identical = Array(ws.sinogram) == sino_reused
t_reuse = @elapsed BS.simulate!(ws, phantom, protocol_standard, opts_2;
report_dose = false, projection = kept.projection)
ws = nothing; kept = nothing
GC.gc(true)
(t_full = t_full, t_reuse = t_reuse, identical = identical)
end;A full simulation of the standard scan took 10.53 s; a noise draw from its kept projection took 0.02 s, 432.4× less. The draw from the projection is bit-identical to simulating the same seed from scratch.
Reconstruct
Water beam-hardening correction
The simulated sinogram is polychromatic: the spectrum hardens as it crosses the patient, so the line integrals grow more slowly than the path length, and an uncorrected reconstruction reads low in the middle of the body (cupping). calibrate_bhc_water removes that for water. For every detector column it resolves the full detected spectrum (tube, filters, bowtie, heel effect and detector efficiency: exactly the factors sim_opts switched on) and fits the polynomial that maps polychromatic water line integrals to monochromatic ones at the spectrum's mean energy. It has no tunable parameters, and nothing in it is fitted to the data being corrected. apply_bhc_water applies it to a sinogram.
The reference energy and the water attenuation μ at that energy, which come with the correction, are what to_hounsfield needs to put water at 0 HU. Both protocols share a beam, so one calibration serves both.
bhc = BS.calibrate_bhc_water(sim_opts, protocol_standard; scanner, geom = sim_std.geom);Water BHC: 834 per-column polynomials · reference energy 69.5 keV · water μ = 0.19335 cm⁻¹
FBP and Hybrid IR through one reconstruct!
reconstruct! dispatches on the reconstruction workspace, so the two algorithms differ only in the workspace you create:
create_fdk_recon_workspace: filtered back-projection (Feldkamp–Davis–Kress for an axial scan, weighted FBP for a helical one) with the:standardkernel by default;create_hir_recon_workspace(...; strength): Hybrid IR, an FBP start followed by ordered-subsets penalised weighted least squares with an edge-preserving Huber prior.strengthis a percentage in steps of 10, read like the vendor dials (GE ASIR-V, Siemens SAFIRE):0is plain FBP,60the standard clinical setting. HIR's system matrix uses the same projector as the simulation (projector = sim_opts.projector), so it inverts the operator that made the data.
The function below is the whole chain, water BHC → reconstruction → HU, and returns a CPU volume.
"""
reconstruct_hu(sim, bhc; algorithm = :fbp, strength = 60) -> Array{Float32, 3}
Water BHC, then FBP (`algorithm = :fbp`) or Hybrid IR (`algorithm = :hir`), then HU.
"""
function reconstruct_hu(sim, bhc; algorithm::Symbol = :fbp, strength::Integer = 60)
sino = BS.apply_bhc_water(to_gpu(sim.sino), bhc)
ws = if algorithm === :fbp
BS.create_fdk_recon_workspace(sino, sim.geom, recon_opts.matrix_size)
elseif algorithm === :hir
BS.create_hir_recon_workspace(sino, sim.geom, recon_opts.matrix_size;
strength, projector = sim_opts.projector)
else
throw(ArgumentError("algorithm must be :fbp or :hir"))
end
μ = BS.reconstruct!(ws, sino, sim.geom)
hu = Float32.(BS.to_hounsfield(Array(μ); μ_water = bhc.μ_water_ref))
ws = nothing; sino = nothing; μ = nothing
GC.gc(true)
return hu
end;hu_fbp_std = reconstruct_hu(sim_std, bhc; algorithm = :fbp);hu_hir_std = reconstruct_hu(sim_std, bhc; algorithm = :hir, strength = 60);hu_fbp_low = reconstruct_hu(sim_low, bhc; algorithm = :fbp);hu_hir_low = reconstruct_hu(sim_low, bhc; algorithm = :hir, strength = 60);Results
Where the noise is measured
Noise is only meaningful inside one material. The phantom labels are carried onto the reconstruction grid with resample_to_recon, and the solid-water label is eroded so that only voxels surrounded by water remain: the rods and their partial-volume edges are excluded. All statistics below use the central slice; the edge slices of a cone-beam reconstruction see fewer complete rays.
FBP vs Hybrid IR at two doses
The four reconstructions share physics, geometry and grid. Down the columns, only the tube current changes; across the rows, only the reconstruction. Quartering the dose should double the FBP noise (noise ∝ 1/√mAs) and leave every mean HU where it was. Hybrid IR at strength 60 should lower the noise at both doses without moving the means.
Water, central slice, 151452 voxels deep inside the solid-water body:
| protocol | FBP mean (HU) | FBP σ (HU) | HIR mean (HU) | HIR σ (HU) | HIR noise reduction |
|---|---|---|---|---|---|
| 200 mA | 1.4 | 54.4 | 1.5 | 34.8 | 36.1 % |
| 50 mA | 1.7 | 121.1 | 1.7 | 76.4 | 36.9 % |
FBP noise ratio, 50 mA / 200 mA = 2.22×. Quantum noise alone predicts 2×; the fixed electronic noise of the detector weighs more at low dose and pushes the ratio a little above that. The package's calibration band for HIR at strength 60 is 25–35 % noise reduction.
Dose and noise at other tube currents
The two scans anchor a simple rule: CTDIvol is proportional to mAs (same beam, same phantom), and quantum-limited FBP noise scales as 1/√mAs. The table applies it to this protocol at other tube currents. The detector's electronic noise makes the real noise somewhat higher than this at low mA, as the measured 50 mA scan shows.
| tube current (mA) | CTDIvol (mGy) | expected FBP water noise (HU) |
|---|---|---|
| 25 | 1.72 | 153.9 |
| 50 | 3.43 | 108.8 |
| 100 | 6.86 | 77.0 |
| 200 | 13.72 | 54.4 |
| 300 | 20.58 | 44.4 |
| 400 | 27.45 | 38.5 |
Measured at 50 mA: 121.1 HU.
@bind tube_mA PlutoUI.Slider(10:10:500; default = 200, show_value = true)At 200 mA: CTDIvol is 13.72 mGy and the expected quantum-limited FBP water noise is 54.4 HU, 1.0 times that of the 200 mA scan.
Scroll through the slices
The reconstructed slab is 5 mm thick. The edge slices are reconstructed from more oblique rays than the central ones, so this is where any z-dependence would show first.
Verification
Pass/fail against physics-derived expectations. The theory for each rod is its monoenergetic HU at the BHC reference energy, computed from the same XrayAttenuation data that drove the simulation. The tolerance is max(15 HU, 15 %): a water-only BHC maps water exactly, but a dense calcium or iodine rod hardens the beam beyond the water curve. That residual grows with density (−12.4 % at Ca 600 mg/mL and −14.5 % at I 20 mg/mL in the table below) and a single-energy scan cannot resolve it; quantitative high-Z imaging is what the dual-energy and photon-counting VMI notebooks (03, 04) are for. Each rod must also read the same at both doses and in both reconstructions.
✅ VERIFICATION PASS (8/8)
| check | value | expected | pass |
|---|---|---|---|
| water mean, FBP 200 mA (HU) | 1.44 | [-6.0, 6.0] | ✅ |
| water mean, HIR 200 mA (HU) | 1.5 | [-6.0, 6.0] | ✅ |
| FBP noise ratio 50 / 200 mA | 2.22 | [1.7, 2.4] | ✅ |
| HIR noise < FBP noise, 200 mA (ratio) | 0.64 | [0.0, 0.95] | ✅ |
| CTDIvol ratio 50 / 200 mA | 0.25 | [0.2499, 0.2501] | ✅ |
| radial cupping QA, FBP (HU, worst slice) | 12.02 | [0.0, 15.0] | ✅ |
| DC offset QA, FBP (absolute HU, worst slice) | 5.88 | [0.0, 6.0] | ✅ |
| rods passing (theory, dose- and algorithm-invariant) | 14.0 | [14, 14] | ✅ |
Per rod, central slice. Theory is the monoenergetic HU at 69.5 keV.
| rod | theory | FBP 200 mA | FBP − theory | FBP 50 mA | HIR 200 mA | pass |
|---|---|---|---|---|---|---|
| Gammex 472 (50.0 mg/ml Calcium) | 175 | 174 | -2.0 (-0.9 %) | 176 | 174 | ✅ |
| Gammex 472 (100.0 mg/ml Calcium) | 315 | 298 | -17.0 (-5.5 %) | 297 | 299 | ✅ |
| Gammex 472 (200.0 mg/ml Calcium) | 617 | 577 | -39.0 (-6.4 %) | 586 | 577 | ✅ |
| Gammex 472 (300.0 mg/ml Calcium) | 907 | 828 | -79.0 (-8.7 %) | 835 | 828 | ✅ |
| Gammex 472 (400.0 mg/ml Calcium) | 1197 | 1071 | -126.0 (-10.6 %) | 1080 | 1071 | ✅ |
| Gammex 472 (500.0 mg/ml Calcium) | 1487 | 1313 | -175.0 (-11.8 %) | 1329 | 1313 | ✅ |
| Gammex 472 (600.0 mg/ml Calcium) | 1792 | 1570 | -222.0 (-12.4 %) | 1594 | 1571 | ✅ |
| Gammex 472 (2.0 mg/ml Iodine) | 56 | 53 | -3.0 (-4.9 %) | 54 | 53 | ✅ |
| Gammex 472 (2.5 mg/ml Iodine) | 69 | 63 | -6.0 (-8.8 %) | 63 | 63 | ✅ |
| Gammex 472 (5.0 mg/ml Iodine) | 132 | 119 | -12.0 (-9.5 %) | 120 | 119 | ✅ |
| Gammex 472 (7.5 mg/ml Iodine) | 195 | 170 | -25.0 (-12.9 %) | 172 | 169 | ✅ |
| Gammex 472 (10.0 mg/ml Iodine) | 257 | 225 | -32.0 (-12.4 %) | 228 | 225 | ✅ |
| Gammex 472 (15.0 mg/ml Iodine) | 385 | 331 | -54.0 (-14.0 %) | 335 | 331 | ✅ |
| Gammex 472 (20.0 mg/ml Iodine) | 526 | 450 | -76.0 (-14.5 %) | 452 | 450 | ✅ |
Summary
Five structs describe a scan.
Phantom(the object),EICTScanner/PCCTScanner(the hardware, sharing aScannerGeometry),CTProtocol(the acquisition),SimOptions(the physics, including the view integration of a turning gantry,view_samples = 5) andReconOptions(the output grid). Changing the dose meant changing oneCTProtocoland nothing else.Three calls produce an image.
create_workspace(for either detector family) sets everything up once,simulate!runs the forward model and returns the dose report, andreconstruct!runs FBP or Hybrid IR, depending on the workspace it is given.A noise draw can reuse the projection.
keep_projection = truekeeps everything before the noise, andprojectiondraws another seed from it, bit-identically to a full simulation.Every scan reports its dose. CTDIvol and DLP come from a Monte Carlo of the simulated beam in the IEC body phantom, so they follow the scanner's spectrum, filtration and bowtie.
The standard reconstruction chain is
calibrate_bhc_water→apply_bhc_water→reconstruct!→to_hounsfieldwith the correction's own water μ. Quantum and electronic noise are already in the simulated counts, andmeasure_radial_cuppingis a QA measurement, not a correction.Where it runs is decided by the phantom mask's array type; this page was rendered on an NVIDIA RTX PRO 6000 (CUDA), and the same notebook runs on Metal, ROCm, oneAPI or the CPU.
Every other notebook reuses this pattern: build the structs, create a workspace, simulate!, reconstruct!.