Getting started
From an empty environment to a reconstructed image, a dose report and a set of virtual monoenergetic images. The listings build on one another and run as written, on a GPU or on the CPU; the sizes are small enough for a laptop.
Install
BasisSimulator.jl requires Julia 1.12 or later. From the REPL:
using Pkg
Pkg.add("BasisSimulator")
Pkg.add("GPUSelect") # picks the device array type for youAdd the package for your GPU, or none to run on the CPU. The kernels are written once against AcceleratedKernels.jl; the example notebooks run on CUDA and Metal.
Pkg.add("CUDA") # NVIDIA
Pkg.add("Metal") # Apple silicon
Pkg.add("AMDGPU") # AMD
Pkg.add("oneAPI") # IntelThe five structs
A simulation is specified by five values, each with one job:
Phantomwhat is scanned: a label mask on the device, one XrayAttenuation material per label, the voxel size
EICTScanner / PCCTScannerthe hardware: geometry, focal spot, filtration, bowtie, and the detector model of one family
CTProtocolthe acquisition: kVp, mA, views, rotation time, collimation, pitch, added filters
SimOptionsthe physics common to both families: which effects are on, the noise seed, the projector
ReconOptionsthe output grid: matrix size, field of view, z extent
Parameter sets for the clinical systems the notebooks model are on the Scanners page.
First simulation
A Gammex 472 phantom on a model of the GE Revolution Apex Elite. The workspace allocates every buffer once; simulate! writes the log line-integral sinogram into ws.sinogram and returns the acquisition's dose report.
import BasisSimulator as BS
import GPUSelect
AT = GPUSelect.Storage() # CuArray, MtlArray, ROCArray or oneArray; Array on the CPU
to_gpu(x) = AT(x)
cpu = BS.create_gammex_472(n_voxels = 256, n_slices = 4, z_cm = 1.0)
phantom = BS.Phantom(to_gpu(cpu.mask), cpu.materials, cpu.voxel_size, cpu.origin, cpu.extent)
scanner = BS.EICTScanner(source_to_isocenter = 625.6, source_to_detector = 1100.0,
detector_rows = 32, detector_cols = 834,
detector_row_size = 0.625, detector_col_size = 0.6)
protocol = BS.CTProtocol(kVp = 120, mA = 200.0, views = 500, collimation_mm = 5.0)
sim_opts = BS.SimOptions(seed = 42)
rec_opts = BS.ReconOptions(matrix_size = (512, 512, 4), fov_cm = 35.0, z_cm = 0.5)
ws = BS.create_workspace(scanner, protocol, sim_opts, rec_opts, phantom)
result = BS.simulate!(ws, phantom, protocol, sim_opts)The phantom's mask lives on the device and its materials on the host. Quantum and electronic noise are drawn in the counts domain, before the log, where a real detector produces them. A photon-counting scanner's workspace comes from the same create_workspace (below).
Reconstruct
The water beam-hardening correction is calibrated from the full detected spectrum of this acquisition (tube, filters, bowtie, heel effect, detector efficiency) and has no tunable parameters. Then FDK, and μ to Hounsfield units against the calibration's water reference:
bhc = BS.calibrate_bhc_water(sim_opts, protocol; scanner, geom = ws.geom)
sino = to_gpu(BS.apply_bhc_water(ws.sinogram, bhc))
fdk = BS.create_fdk_recon_workspace(sino, ws.geom, rec_opts.matrix_size)
hu = BS.to_hounsfield(Array(BS.reconstruct!(fdk, sino, ws.geom)); μ_water = bhc.μ_water_ref)Hybrid iterative reconstruction takes the same sinogram. Its one dial is strength, 0–100 in steps of 10: 0 is plain FBP, 60 the standard clinical setting.
hir = BS.create_hir_recon_workspace(sino, ws.geom, rec_opts.matrix_size; strength = 60)
hu_hir = BS.to_hounsfield(Array(BS.reconstruct!(hir, sino, ws.geom)); μ_water = bhc.μ_water_ref)Keep create_hir_recon_workspace(; projector) equal to SimOptions(; projector) (both default to :dd_fast) so the reconstruction inverts the operator that made the data.
Dose
Every simulate! result carries a DoseReport computed by Monte Carlo from the simulated beam, in the 32 cm body CTDI phantom by default:
result.dose # CTDIvol, DLP, CTDIw, CTDI100 centre and periphery
result.dose.ctdi_vol_mGy # 21.7 mGy for the scan above
result.dose.dlp_mGy_cm
# a head protocol, or two tubes: keywords reach compute_dose
BS.simulate!(ws, phantom, protocol, sim_opts; dose_kwargs = (; phantom = :head16))Helical and multi-kVp
One keyword, pitch, makes a protocol helical; reconstruct! then dispatches to rebinned WFBP automatically.
spiral = BS.CTProtocol(kVp = 120, mA = 200.0, views = 360,
collimation_mm = 10.0, pitch = 1.0, n_rotations = 2)The volume walk depends on the geometry and the material map, not on the spectrum, so a study of one phantom at several tube voltages walks once and reuses the path lengths, bit-identically:
paths = BS.material_paths(ws, phantom)
for kvp in (80, 100, 120, 140)
p = BS.CTProtocol(kVp = kvp, mA = 200.0, views = 500, collimation_mm = 5.0)
w = BS.create_workspace(scanner, p, sim_opts, rec_opts, phantom)
BS.simulate!(w, phantom, p, sim_opts; paths)
endPhoton counting → VMI
A PCCTScanner adds the counting detector: thresholds, energy resolution, charge sharing, dead time, and the pile-up and scatter corrections applied inside simulate!. The result holds one corrected log sinogram per bin and the air counts of every ray.
scanner_pc = BS.PCCTScanner(
source_to_isocenter = 610.0, source_to_detector = 1113.0,
detector_rows = 32, detector_cols = 1200, # 1200 × 0.3 mm covers the 33 cm phantom
detector_row_size = 0.35, detector_col_size = 0.3,
energy_thresholds = [20.0, 35.0, 55.0, 70.0], # four counting bins, keV
energy_resolution = 10.0, charge_sharing_fwhm = 0.08,
dead_time_ns = 5.0, # pile-up is modelled once there is a dead time
pileup_correction = true, scatter_correction = true,
)
ws_pc = BS.create_workspace(scanner_pc, protocol, sim_opts, rec_opts, phantom)
res_pc = BS.simulate!(ws_pc, phantom, protocol, sim_opts)A dead time switches on the Monte Carlo pile-up model, which create_workspace runs once, on the CPU, for the rates the scan will see: several minutes, after which simulate! takes seconds. dead_time_ns = 0 skips pile-up while iterating.
spectral_basis takes the detected response straight from the workspace, the model that generated the counts, so the decomposition needs no calibration scan. vmi_pipeline runs the n-channel profile-likelihood decomposition into water and iodine, reconstructs the basis pair and synthesises virtual monoenergetic images in HU:
channels = [Array(b) for b in res_pc.pcct_sino.bins]
basis = BS.spectral_basis(ws_pc; I0 = res_pc.I0_bins)
vmi = BS.vmi_pipeline(; channels, basis, geom = ws_pc.geom, to_backend = to_gpu,
matrix_size = rec_opts.matrix_size)
vmi.vmis # (512, 512, 4, 4): 40, 70, 100 and 140 keV
vmi.images # the reconstructed (water, iodine) basis pairEvery stage is a keyword: method = :cong for the Cong estimator; reduce_rows = true with use_tlbf = true for a z-invariant phantom's photon-counting configuration; recon_method = :hir to reconstruct the basis pair iteratively; vmi_energies for other energies. See ?BS.vmi_pipeline.
Dual-kVp and dual-source
The same estimator takes any number of separate acquisitions. Rapid kVp switching and two tubes both deliver a low and a high channel, each with the absolute response its own forward model applied; spectral_basis_from_acquisitions assembles them on the union of their energy grids.
function acquire(protocol)
w = BS.create_workspace(scanner, protocol, sim_opts, rec_opts, phantom)
BS.simulate!(w, phantom, protocol, sim_opts)
air = w.bowtie_air_reference === nothing ? ones(w.geom.n_cols, w.geom.n_rows) :
Array(w.bowtie_air_reference)
I0 = BS.compute_detector_I0(w.geom, protocol, sum(w.weights)) * w.η_eff
energies, response = BS.resolve_source_spectrum_full(sim_opts, protocol; scanner, geom = w.geom)
(sino = Array(w.sinogram), geom = w.geom, energies, response, I0_ray = I0 .* air)
end
# rapid kVp switching; for two tubes, add e.g. additional_filters = [("Sn", 0.4)] on the high one
low = acquire(BS.CTProtocol(kVp = 80, mA = 407 * 0.65, views = 500, collimation_mm = 5.0))
high = acquire(BS.CTProtocol(kVp = 140, mA = 405 * 0.35, views = 500, collimation_mm = 5.0))
basis_de = BS.spectral_basis_from_acquisitions(acquisitions = [
(; low.energies, low.response, low.I0_ray),
(; high.energies, high.response, high.I0_ray)])
vmi_de = BS.vmi_pipeline(; channels = [low.sino, high.sino], basis = basis_de, geom = low.geom,
to_backend = to_gpu, matrix_size = rec_opts.matrix_size)SpectralHYPR denoising
SpectralHYPR is a generalized HYPR-LR denoiser with two instances: one on the counts of each detector row before the decomposition, one on the reconstructed basis pair. Pass it as denoiser to either call above. With a PairFilter — an FDK window for the composite (the minimum-noise VMI) and one for its noise-independent complement — this is the chain of basis-spectral-denoising:
pair = BS.PairFilter(BS.SoftFilter(), BS.SoftFilter()) # the composite's window, the complement's
vmi_denoised = BS.vmi_pipeline(; channels, basis, geom = ws_pc.geom, to_backend = to_gpu,
matrix_size = rec_opts.matrix_size,
denoiser = BS.SpectralHYPR(), fbp_filter = pair)SpectralHYPR() pools each ray's count split over its 3 × 3 (column × view) neighbours in one view parity, and on the image the complement over the window the data's own risk estimate selects; ACNR (on by default) acts on the complement. The worked VMI examples use each scanner's fitted pair of windows and fix the pair's basis from a calibration draw (pair_basis, composite_energy).