BasisSimulator.jl
← Examples/

CatSim vs BasisSimulator: Agreement and Runtime

One scanner, one protocol, one Gammex 472 phantom, simulated and reconstructed three ways: XCIST/CatSim (GE's open-source reference simulator, Python + C), BasisSimulator on the CPU, and BasisSimulator on the GPU.

This is a cross-validation of the forward model and the reconstruction against an independent implementation. Both simulators are configured to the same beam and the same physics, both apply a water beam-hardening correction referenced to the same monoenergetic water μ, and both reconstruct with FDK and a standard kernel. The comparison therefore reads in HU: the per-rod table below puts CatSim, BasisSimulator and the XrayAttenuation theory side by side. The runtime of each pipeline comes after.

Noise is off in all three runs. The two simulators calibrate tube flux independently, so a noisy comparison would mix a flux-model difference into what should be a projector and reconstruction check.

Gammex 472 @ 128² × 8 (2.7 mm voxels)
   → GE Apex Elite geometry · 120 kVp / 200 mA / 500 views · 7 mm Al · large-body bowtie
   →┬→ CatSim (Python):  forward project → water BHC → FDK
    ├→ BasisSim CPU:     simulate! → water BHC → reconstruct!
    └→ BasisSim GPU:     simulate! → water BHC → reconstruct!   (same code, device mask)
   → images · differences · per-rod HU · runtime

Notebook setup

This is the only docs notebook that uses Python. docs/CondaPkg.toml pins gecatsim to the MolloiLab fork of XCIST (git+https://github.com/MolloiLab/main), which adds the Gammex 472 material definitions, together with numpy and pydicom. The first import PythonCall in the docs/ environment makes CondaPkg install Python and those packages into docs/.CondaPkg/ (a few minutes and about 1 GB, once). If PythonCall or gecatsim cannot be loaded, the CatSim cells skip with a notice and the BasisSimulator CPU and GPU runs still execute.

import BasisSimulator as BS
# PythonCall, loaded so that a missing Python environment degrades to "CatSim skipped"
PC = try
    Base.require(Base.PkgId(Base.UUID("6099a3de-0909-46bc-b1f4-468b9a2dfc0d"), "PythonCall"))
catch err
    @warn "PythonCall could not be loaded; the CatSim cells will skip" exception = (err, catch_backtrace())
    nothing
end;
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

const HAS_GECATSIM = PC !== nothing && try
    PC.pyimport("gecatsim")
    true
catch err
    @warn "gecatsim is not importable; the CatSim cells will skip" exception = err
    false
end;

gecatsim 1.5.4 loaded (Python 3.14.6): the CatSim cells run.

Scan and phantom set-up

One scanner, one protocol and one phantom, shared verbatim by the three pipelines, and a thin wrapper that hands them to CatSim.

1. Scanner

The GE Revolution Apex Elite geometry of notebook 01, with its large-body bowtie. BasisSimulator's :ge_revolution_large bowtie is read from the same large.txt table CatSim ships, so both simulators see the same bowtie. Electronic noise is zero, as noise is off everywhere here.

scanner = BS.EICTScanner(
    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,
    focal_spot_width = 1.0,
    focal_spot_length = 1.0,
    target_angle = 10.0,
    flat_filter_material = :aluminum,
    flat_filter_thickness = 2.5,
    bowtie_filter = :ge_revolution_large,
    detector_material = :lumex,
    detector_depth = 3.0,
    fill_factor_row = 0.9,
    fill_factor_col = 0.9,
    electronic_noise = 0,
    detection_gain = 10.0,
);

2. Protocol, physics and grid

120 kVp / 200 mA, 500 views in a 1 s rotation, 4 mm of collimation, 4.5 mm of extra aluminium. The physics is set to what both simulators model by default: polychromatic transport through the filters and bowtie, energy-dependent detector efficiency, fill factor and focal-spot blur. Noise, scatter, detector lag and the heel effect are switched off in BasisSimulator because the CatSim defaults have them off (CatSim leaves its scatter, lag and crosstalk callbacks empty). The reconstruction is 256 × 256 over 35 cm, 0.625 mm slices.

Point views on both sides

The other notebooks model a clinical scanner's detector integrating each view while the gantry turns (SimOptions(; view_samples = 5), notebook 01). This one does not: it is a like-for-like check of the projector and the reconstruction, so both simulators sample each view at a single angle. BasisSimulator keeps its default view_samples = 1, and the wrapper sets CatSim's physics.viewSampleCount (2 by default) to 1.

protocol = BS.CTProtocol(
    kVp = 120,
    mA = 200.0,
    views = 500,
    rotation_time = 1.0,
    collimation_mm = 4.0,
    additional_filters = [("Al", 4.5)],
);
sim_opts = BS.SimOptions(
    seed = 1234, projector = :dd_fast,
    use_noise = false, use_scatter = false, use_lag = false, use_heel_effect = false,
);
recon_opts = BS.ReconOptions(
    matrix_size = (256, 256, round(Int, protocol.collimation_mm / 0.625)),
    fov_cm = 35.0,
    z_cm = protocol.collimation_mm / 10,
);

3. The CatSim wrapper

A few Julia functions map EICTScanner / CTProtocol / ReconOptions / Phantom onto CatSim's configuration:

  • Geometry. BasisSimulator gives the detector pitch at the isocentre; CatSim wants it at the detector face, so it is multiplied by SDD/SID. Every column is its own module with no gap (detectorColsPerMod = 1, detectorColSkip = 0); without that CatSim braids or squashes the sinogram.

  • Beam. CatSim's 120 kVp tungsten spectrum for a 10° target (the protocol's anode angle), the scanner's flat filter plus the protocol's filters, the same large.txt bowtie, the same target angle and focal spot, and no graphite detector prefilter (BasisSimulator models none).

  • Views. One angular sample per view (physics.viewSampleCount = 1), the point views of BasisSimulator's view_samples = 1.

  • Preprocessing. CatSim's own water BHC, Prep_BHC_Accurate: a degree-5 polynomial per detector cell fitted to air scans through 1–50 cm of water, mapped to the same monoenergetic water μ that BasisSimulator's calibrate_bhc_water uses. HU conversion uses that water μ too.

function catsim_configure_scanner!(ct, scanner, protocol)
    magnification = scanner.source_to_detector / scanner.source_to_isocenter
    n_rows = protocol.collimation_mm === nothing ? scanner.detector_rows :
        round(Int, protocol.collimation_mm / scanner.detector_row_size)

    ct.scanner.sid = scanner.source_to_isocenter
    ct.scanner.sdd = scanner.source_to_detector
    ct.scanner.detectorColCount = scanner.detector_cols
    ct.scanner.detectorRowCount = n_rows
    ct.scanner.detectorColSize = scanner.detector_col_size * magnification   # isocentre → face
    ct.scanner.detectorRowSize = scanner.detector_row_size * magnification
    ct.scanner.detectorColOffset = scanner.detector_col_offset
    ct.scanner.detectorColsPerMod = 1          # every column its own module …
    ct.scanner.detectorRowsPerMod = n_rows
    ct.scanner.detectorColSkip = 0.0           # … with no inter-module gap
    ct.scanner.detectorRowSkip = 0.0

    ct.scanner.targetAngle = scanner.target_angle
    ct.scanner.focalspotWidth = scanner.focal_spot_width
    ct.scanner.focalspotLength = scanner.focal_spot_length
    ct.scanner.detectorMaterial = "Lumex"
    ct.scanner.detectorDepth = scanner.detector_depth
    ct.scanner.detectorColFillFraction = scanner.fill_factor_col
    ct.scanner.detectorRowFillFraction = scanner.fill_factor_row
    ct.scanner.detectorPrefilter = PC.pylist([])
    return ct
end;
function catsim_configure_protocol!(ct, scanner, protocol; μ_water_cm)
    ct.protocol.mA = protocol.mA
    ct.protocol.viewsPerRotation = protocol.views
    ct.protocol.viewCount = protocol.views
    ct.protocol.stopViewId = protocol.views - 1
    ct.protocol.rotationTime = protocol.rotation_time
    ct.protocol.spectrumFilename = "tungsten_tar$(protocol.anode_angle).0_$(Int(protocol.kVp))_filt.dat"
    scanner.bowtie_filter in (:ge_revolution_large, :large_body) ||
        error("this wrapper maps only the large-body bowtie")
    ct.protocol.bowtie = "large.txt"
    catsim_name(m) = Dict("aluminum" => "Al", "Al" => "Al", "copper" => "Cu", "Cu" => "Cu")[string(m)]
    filters = Any[catsim_name(scanner.flat_filter_material), scanner.flat_filter_thickness]
    for (m, t) in protocol.additional_filters
        push!(filters, catsim_name(m), t)
    end
    ct.protocol.flatFilter = PC.pylist(filters)

    ct.physics.viewSampleCount = 1      # point views, like BasisSimulator's view_samples = 1
    ct.physics.enableQuantumNoise = 0
    ct.physics.enableElectronicNoise = 0
    ct.physics.callback_post_log = "Prep_BHC_Accurate"   # CatSim's water BHC
    ct.physics.EffectiveMu = μ_water_cm                   # cm⁻¹, the mono target of the fit
    ct.physics.BHC_poly_order = 5
    ct.physics.BHC_max_length_mm = 500
    ct.physics.BHC_length_step_mm = 10
    return ct
end;
function catsim_configure_recon!(ct, recon_opts; μ_water_cm)
    xc, _, cfg = catsim_init()
    xc.source_cfg(joinpath(cfg, "Recon_Sample_2d.cfg"), ct)
    n_slices = recon_opts.matrix_size[3]
    ct.recon.fov = recon_opts.fov_cm * 10.0
    ct.recon.imageSize = recon_opts.matrix_size[1]
    ct.recon.sliceCount = n_slices
    ct.recon.sliceThickness = recon_opts.z_cm * 10.0 / n_slices
    ct.recon.reconType = "fdk_equiAngle"
    ct.recon.kernelType = "standard"
    ct.recon.unit = "HU"
    ct.recon.mu = μ_water_cm / 10.0           # cm⁻¹ → mm⁻¹
    ct.recon.huOffset = -1000
    return ct
end;

Phantom → CatSim voxelized JSON

create_gammex_472 labels the solid-water body 3, the calcium rods 10–16 and the iodine rods 20–26; the MolloiLab fork defines the matching Gammex472_* materials. CatSim's Phantom_Voxelized reads a JSON header plus one Float32 volume-fraction map per material; for a hard-segmented mask each map is the indicator mask .== label.

4. The Gammex 472 phantom

n_voxels = 128 gives 2.7 mm voxels over 35 cm: coarse enough that CatSim's voxelized projector finishes in minutes, fine enough to resolve the 28 mm rods. Eight slices cover the 4 mm beam.

phantom_cpu = BS.create_gammex_472(
    n_voxels = 128,
    n_slices = 8,
    fov_cm = 35.0,
    z_cm = protocol.collimation_mm / 10.0,
);
phantom_gpu = BS.Phantom(
    to_gpu(phantom_cpu.mask),
    phantom_cpu.materials,
    phantom_cpu.voxel_size,
    phantom_cpu.origin,
    phantom_cpu.extent,
);

5. The shared water reference

calibrate_bhc_water resolves BasisSimulator's detected spectrum per detector column and returns the per-column correction with its monoenergetic reference (the spectrum's mean energy). That reference water μ is given to CatSim as the target of its own BHC fit and used for the HU conversion of all three reconstructions.

geom_inspect = BS.CTGeometry(
    scanner;
    n_angles = protocol.views,
    fov_cm = recon_opts.fov_cm,
    z_cm = recon_opts.z_cm,
    collimation_mm = protocol.collimation_mm,
);
bhc = BS.calibrate_bhc_water(sim_opts, protocol; scanner, geom = geom_inspect);

Water reference: 69.5 keV, water μ = 0.19335 cm⁻¹

Run the three pipelines

Each timing covers the whole job from the phantom to an HU volume: forward projection, preprocessing and BHC, FDK. Both BasisSimulator runs are compiled by one untimed call first, so the timings are steady-state. CatSim's includes writing and reading its intermediate files and fitting its BHC polynomials.

1. CatSim

catsim_result = !HAS_GECATSIM ? nothing : let
    work_dir = mktempdir(; prefix = "basissim_catsim_06_")
    tag = joinpath(work_dir, "gammex472")
    local recon
    elapsed = @elapsed begin
        json = export_phantom_for_catsim(phantom_cpu, work_dir, "gammex472")
        ct = catsim_create_simulation()
        catsim_configure_phantom!(ct, json)
        catsim_configure_scanner!(ct, scanner, protocol)
        catsim_configure_protocol!(ct, scanner, protocol; μ_water_cm = bhc.μ_water_ref)
        catsim_configure_recon!(ct, recon_opts; μ_water_cm = bhc.μ_water_ref)
        catsim_forward_project(ct; results_name = tag)
        recon = catsim_reconstruct_fdk(ct; results_name = tag)
    end
    rm(work_dir; recursive = true, force = true)
    (recon = recon, elapsed = elapsed)
end;

2. BasisSimulator, CPU and GPU

One function, two phantoms: phantom_cpu.mask is an Array, so the first run stays on the host; phantom_gpu.mask lives on the GPU, so the second runs there. Nothing else differs.

"""
    basissim_pipeline(phantom) -> HU volume

`create_workspace` → `simulate!` → water BHC → FDK → HU, on the phantom's backend.
"""
function basissim_pipeline(phantom)
    ws = BS.create_workspace(scanner, protocol, sim_opts, recon_opts, phantom)
    BS.simulate!(ws, phantom, protocol, sim_opts; report_dose = false)
    sino = BS.apply_bhc_water(ws.sinogram, bhc)
    ws_fdk = BS.create_fdk_recon_workspace(sino, ws.geom, recon_opts.matrix_size; filter = :standard)
    μ = Array(BS.reconstruct!(ws_fdk, sino, ws.geom))
    return Float32.(BS.to_hounsfield(μ; μ_water = bhc.μ_water_ref))
end;
basissim_cpu_result = let
    basissim_pipeline(phantom_cpu); GC.gc(true)                 # compile, untimed
    elapsed = @elapsed recon = basissim_pipeline(phantom_cpu)
    GC.gc(true)
    (recon = recon, elapsed = elapsed)
end;
basissim_gpu_result = let
    basissim_pipeline(phantom_gpu); GC.gc(true)                 # compile, untimed
    elapsed = @elapsed recon = basissim_pipeline(phantom_gpu)
    GC.gc(true)
    (recon = recon, elapsed = elapsed)
end;

Results

Image agreement

The central slice of each reconstruction in one HU window, and the difference of each BasisSimulator image from CatSim's. CatSim writes its image with its own axis convention; it is brought onto BasisSimulator's orientation by the flip or transpose (of the eight in-plane symmetries) that best matches the BasisSimulator image, which the figure reports. The thin rings along the body and rod edges are where the two images differ in edge sharpness; the flat interiors are what the per-rod table below measures.

Per-rod HU

Mean HU inside each rod (the phantom labels resampled onto the reconstruction grid, eroded by one voxel) on the central slice, next to the rod's theoretical monoenergetic HU at the water reference energy. A water-only BHC leaves dense calcium and iodine rods below their theory in both simulators (the rod's own beam hardening beyond the water curve; notebook 01 discusses it). What this table tests is whether the two simulators agree with each other.

regionvoxelstheory (HU)CatSim (HU)BasisSim (HU)BasisSim − CatSim÷ theory
solid water body389560.03.80.3-3.5
Gammex 472 (50.0 mg/ml Calcium)252175.0170.7175.24.52.6 %
Gammex 472 (100.0 mg/ml Calcium)248315.0291.0303.312.23.9 %
Gammex 472 (200.0 mg/ml Calcium)252617.0550.5580.630.14.9 %
Gammex 472 (300.0 mg/ml Calcium)248907.0792.7836.343.64.8 %
Gammex 472 (400.0 mg/ml Calcium)2481197.01028.21081.953.64.5 %
Gammex 472 (500.0 mg/ml Calcium)2521487.01258.71322.263.54.3 %
Gammex 472 (600.0 mg/ml Calcium)2481792.01506.91581.174.24.1 %
Gammex 472 (2.0 mg/ml Iodine)24856.050.753.83.15.6 %
Gammex 472 (2.5 mg/ml Iodine)25669.059.863.53.75.4 %
Gammex 472 (5.0 mg/ml Iodine)256132.0111.4121.39.97.5 %
Gammex 472 (7.5 mg/ml Iodine)244195.0156.8173.516.78.6 %
Gammex 472 (10.0 mg/ml Iodine)256257.0205.4229.924.59.5 %
Gammex 472 (15.0 mg/ml Iodine)256385.0301.7338.336.69.5 %
Gammex 472 (20.0 mg/ml Iodine)248526.0411.9461.549.79.4 %

Water: BasisSim − CatSim = -3.5 HU. Over the 14 rods the difference is 2.6 to 9.5 % of the rod's theoretical HU. Both simulators read every dense rod below its theory, and they differ mainly in how far: the two use independent source-spectrum models (IPEM tables in BasisSimulator, CatSim's own tungsten spectra) and independently fitted water corrections, and a spectrum difference would show exactly there, in the beam-hardening residual that grows with rod density. BasisSimulator CPU vs GPU: largest voxel difference 0.00922 HU over the whole volume (the same code on two backends; floating-point summation order is all that differs).

Runtime

Wall-clock of each end-to-end pipeline, on a log scale so the GPU bar stays visible. The CPU bar is the like-for-like comparison with CatSim (both run on the host); the GPU bar is what BasisSimulator is built for. This page was rendered on an NVIDIA RTX PRO 6000 (CUDA).

pipelinewall-clockspeed-up vs CatSim
CatSim (voxelized projector, Python + C)48.57 s1.0× (reference)
BasisSimulator.jl, CPU (16 Julia threads)7.57 s6.4×
BasisSimulator.jl, CuArray0.76 s64.3×

Why BasisSimulator is faster

  • No disk round trip. CatSim writes .air, .offset, .scan and .prep files and reads the .prep back for FDK, then writes the image volume. BasisSimulator keeps the sinogram on the device from forward projection through FBP.

  • One volume walk for the whole spectrum. The default :dd_fast projector accumulates the path length through each material once and weights all energies from it, instead of tracing every energy separately.

  • The same kernels on any backend. CatSim's voxelized projector is C on the host. BasisSimulator's kernels are written once with AcceleratedKernels.jl and run on CUDA, Metal, ROCm, oneAPI or the CPU, as the CPU and GPU runs above show.

The phantom here is deliberately small so that CatSim finishes in minutes. For a larger phantom or a clinical-size detector, scale n_voxels and the scanner; the wrapper forwards any EICTScanner / CTProtocol pair, with the bowtie restriction noted in catsim_configure_protocol!.