← Examples/02 · XCAT Anatomy + Custom Materials

XCAT Phantom + Custom Materials

Scan a digital anatomy with material attenuation you define from scratch.

Notebook 01 walked the five-struct API on the Gammex 472 calibration phantom. This notebook swaps the phantom: instead of a labeled cylinder of inserts, we load a high-resolution NCAT/XCAT voxel phantom — a digital adult torso — and assign each region (lung, blood, muscle, bone, …) its own XrayAttenuation.Material with explicit elemental composition and density.

We finish with a side-by-side FBP vs Hybrid IR reconstruction comparison on the same scan.


Phantom	XCAT adult-male 50th-percentile chest, one axial slab at full resolution (XCIST `phantoms-voxelized`, downloaded as a Julia artifact)
Materials	`XA.Materials.ncat_*` tissues, plus a custom iodinated-blood material constructed inline
Scanner	GE Revolution Apex Elite (as in notebook 01)
Reconstruction	FBP and Hybrid IR at strength 60, both after the water beam-hardening correction

Notebook Setup

As in notebook 01: the shared docs/ environment, then the device. Unitful supplies the u"eV" and u"g/cm^3" units the XA.Material constructor takes.

import BasisSimulator as BS

Choosing the device

GPUSelect.Storage() returns the array type of the GPU it finds (CuArray, MtlArray, ROCArray, oneArray), or Array on a CPU-only host; moving the phantom mask to it decides where the pipeline runs.

begin
    import GPUSelect
    AT = GPUSelect.Storage()     # the backend array type, directly: MtlArray / CuArray / ROCArray
    to_gpu(x) = AT(x)
    GPU_BACKEND = (name = string(nameof(AT)),)
end;

Backend detected: CuArray

Phantom Construction

1. Download the XCAT phantom

No local files, no manual setup. BS.load_xcat_male_chest() fetches the voxelized XCAT adult-male 50th-percentile chest phantom from the open xcist/phantoms-voxelized repository (BSD-3-Clause) into Julia's content-addressed artifact store: it downloads once, verifies the file's SHA-256 against the value pinned in the package, logs the Segars XCAT + XCIST citation, and reuses the cached copy on every later call.

Using the XCAT loaders

Four phantoms ship — :female_slab, :female_chest, :male_slab, :male_chest (see BS.xcat_phantoms()). Each load_xcat_<name>() returns the labeled pieces: mask, a materials dict (label => XA.Material), label_names, and voxel_size_cm — which you turn into a Phantom in one line (§3). Useful keywords:

downsample = n — label-preserving majority-vote shrink (physical extent preserved). This notebook uses the slab at full 0.25 mm resolution, so it passes no downsample; reach for it on the full :male_chest volume. It is lossy — a 6× round-trip agrees with full-res on ~96% of voxels — so treat it as a speed knob, not a lossless representation.
materials = Dict("ncat_blood" => …) — override the tissue → material map by XCIST name before assembly (see §2).
path = … — skip the download and read an already-extracted copy from disk.
quiet = true — silence the citation + progress logging.

2. Custom materials via `XrayAttenuation`

The loader assigns every label a sensible default tissue, but the whole point of a digital phantom is that you control the material physics. Each label maps to an XA.Material that gives the simulator its energy-dependent μ(E). Two ways to get one:

Prebuilt — XA.Materials.ncat_blood, XA.Materials.ncat_muscle, etc. XrayAttenuation ships the canonical NCAT/XCAT tissue compositions.
Constructed inline — XA.Material(name, ZA_ratio, I, density, composition) for your specific contrast bolus, calibration solution, or alloy.

a. Pull a prebuilt NCAT material

XrayAttenuation ships with the full NCAT/XCAT tissue catalog. Inspect one:

XA.Materials.ncat_blood

Field	Value
`name`	NCAT Blood
`ZA_ratio`	0.55
`I`	75.2 eV
`density`	1.06 g cm^-3
`composition`	10 elements (Z = [1, 6, 7, 8, 11, 15, 16, 17, 19, 26])

b. Construct a custom contrast material

Below: a 5 mg/mL iodinated blood mixture, like a coronary CT-angiography bolus. We start from ncat_blood's composition and add elemental iodine.

Composition arithmetic

Adding iodine_mg_per_mL mg of iodine per mL of blood at density 1.06 g/cm³ gives an iodine mass fraction of iodine_mg_per_mL / (1000 · density_g_per_cm3). We rescale the other elements proportionally so the new composition still sums to 1.

function build_iodine_blood(iodine_mg_per_mL::Real)
    base = BS.XA.Materials.ncat_blood
    density = base.density                                   # ~1.06 g/cm³
    ρ_g_cm3 = ustrip(uconvert(u"g/cm^3", density))
    f_I = iodine_mg_per_mL / (1000.0 * ρ_g_cm3)          # iodine mass fraction
    scale = 1.0 - f_I

    comp = Dict{Int, Float64}()
    for (Z, frac) in base.composition
        comp[Z] = frac * scale
    end
    comp[53] = get(comp, 53, 0.0) + f_I                       # add iodine

    return BS.XA.Material(
        "Iodine-blood ($(iodine_mg_per_mL) mg/mL)",
        base.ZA_ratio,                                        # ZA shifts ~negligibly at clinical dose
        base.I,                                               # likewise for excitation energy
        density,
        comp,
    )
end;

iodine_blood = build_iodine_blood(5.0);

c. The canonical name → material map

The loader ships the canonical NCAT/XCAT mapping as BS.xcat_default_materials() — every XCIST material name → an XrayAttenuation material. This is what replaces hand-parsing a per-organ spreadsheet: the compositions already live in XrayAttenuation's NCAT catalog. Copy the dict, edit it by name, and hand it back to the loader via materials = (§d).

BS.xcat_default_materials()

XCIST name	XrayAttenuation material
`ncat_adipose`	`NCAT Fat/Adipose`
`ncat_blood`	`NCAT Blood`
`ncat_cartilage`	`Cartilage (ICRU-44)`
`ncat_dry_rib`	`NCAT Rib`
`ncat_dry_spine`	`NCAT Spine`
`ncat_heart`	`NCAT Heart`
`ncat_intestine`	`Intestine (ICRU-44)`
`ncat_liver`	`NCAT Liver`
`ncat_lung`	`NCAT Lung`
`ncat_muscle`	`NCAT Muscle`
`ncat_spleen`	`NCAT Spleen`
`ncat_water`	`Water, Liquid`

d. Override a material

Copy the default map and reassign by XCIST name — here we swap the blood pool for the 5 mg/mL iodine_blood we built in §b, so the heart chambers and great vessels light up like a contrast-enhanced CTA. The edited map is passed to the loader via materials = in §3.

materials_map = let m = BS.xcat_default_materials()
    m["ncat_blood"] = iodine_blood        # dope the blood pool with 5 mg/mL iodine contrast
    m
end;

3. Build the `Phantom`

load_xcat_male_slab returns the labeled pieces; we assemble the Phantom from them. For a GPU simulation the mask must live on the device — the workspace picks its compute backend from the mask's array type — so we to_gpu it first. The slab is a single axial section at the phantom's full in-plane resolution: real anatomy at a docs-friendly compute cost, with no downsampling.

xcat = try
    BS.load_xcat_male_slab(; materials = materials_map, quiet = true)   # full-resolution single slab
catch err
    @warn "XCAT download/load failed; the compute cells will skip" exception = err
    nothing
end;

phantom_labeled = xcat === nothing ? nothing : xcat.mask;

VOXEL_SIZE_CM = xcat === nothing ? (0.15, 0.15, 0.15) : xcat.voxel_size_cm;

phantom = xcat === nothing ? nothing :
    BS.Phantom(to_gpu(xcat.mask), xcat.materials, xcat.voxel_size_cm);

Scan Setup & Simulation

1. Scanner, protocol, sim & recon options

The GE Revolution Apex Elite of notebook 01. The protocol is a body CTA: 120 kVp / 250 mA, 5 mm collimation, 500 views in 1 s, each view integrated over the arc the gantry turns through while the detector reads it (view_samples = 5). Reconstruction: 512 × 512 over 35 cm, eight 0.625 mm slices.

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, in electrons, enters the counts before the log, as a real DAS does,
    # so it passes through the reconstruction and iterative reconstruction can act on it.
    electronic_noise = 3500.0,   # e⁻ rms
    detection_gain = 10.0,
);

protocol = BS.CTProtocol(
    kVp = 120,
    mA = 250.0,
    views = 500,
    rotation_time = 1.0,
    collimation_mm = 5.0,
    additional_filters = [("Al", 4.5)],
);

# `projector` picks the forward ray tracer: :dd_fast (default) is distance-driven and
# anti-aliased, and walks the volume once for the whole spectrum; :siddon point-samples
# the volume and can alias in strongly beam-hardened regions. The Hybrid-IR cell reads
# `sim_opts.projector`, so its system matrix always matches the operator that made the data.
# `view_samples = 5`: the detector integrates each view while the gantry turns (notebook 01).
sim_opts = BS.SimOptions(seed = 1234, projector = :dd_fast, view_samples = 5);

recon_opts = BS.ReconOptions(
    matrix_size = (512, 512, 8),
    fov_cm = 35.0,
    z_cm = 0.5,
);

2. Forward project

The notebook 01 pattern: create_workspace, simulate!, copy the sinogram off the device, release the device buffers. simulate! returns the acquisition's dose report.

sim = phantom === nothing ? nothing : let
        @info "Simulating XCAT body CTA: 120 kVp / 250 mA…"
        ws = BS.create_workspace(scanner, protocol, sim_opts, recon_opts, phantom)
        out = BS.simulate!(ws, phantom, protocol, sim_opts)

        result = (sino = Array(ws.sinogram), geom = ws.geom, dose = out.dose)

        ws = nothing
        GC.gc(true)

        result
end;

DoseReport (body32, 120 kVp, 250.0 mAs/rotation, axial, increment 5.0 mm) CTDIvol = 17.15 mGy DLP = 8.58 mGy·cm over 0.5 cm CTDIw = 6.861 mGy/100 mAs (centre 4.151, periphery 8.217; ±0.48 % statistical) free-in-air kerma at isocentre = 16.16 mGy/100 mAs N·T = 5.0 mm, beam 5.0 mm, filters [("aluminum", 2.5), ("Al", 4.5)], bowtie :ge_revolution_large

sim === nothing ? md"" : sim.dose

Reconstruction

1. The reconstruction chain

The chain of notebook 01, applied to both algorithms:

step	function	what it does
1	`calibrate_bhc_water`	per detector column, the polynomial mapping polychromatic water line integrals to monochromatic ones, from the full detected spectrum; no tunable parameters
2	`apply_bhc_water`	applies it to the sinogram, before reconstruction
3	`reconstruct!`	FBP (`create_fdk_recon_workspace`) or Hybrid IR (`create_hir_recon_workspace`)
4	`to_hounsfield`	μ → HU with the correction's own water μ at its reference energy

Noise is already in the counts

Quantum noise and the electronic noise of scanner.electronic_noise are both added by simulate! to the detector counts, before the log, where a real detector adds them. The reconstruction sees all of it, so iterative reconstruction can reduce all of it. Adding noise to the HU volume after reconstruction would be identical for FBP and HIR by construction.

a. Calibrate the water BHC

One calibration for the 120 kVp beam (2.5 mm flat filter + 4.5 mm Al + bowtie). It returns a WaterBHC with one polynomial per detector column and its reference μ_water_ref, used by both reconstructions below.

bhc_calibration = sim === nothing ? nothing : let
        # Parameter-free water BHC — physics only, no tunables, no thresholds:
        # per-column poly→mono polynomials from the FULL detected spectrum
        # (tube × filters × bowtie × heel × η(E) — whatever sim_opts enabled).
        model = BS.calibrate_bhc_water(
            sim_opts, protocol;
            scanner = scanner, geom = sim.geom,
        )

        (
            model = model,
            μ_water = model.μ_water_ref,
            ref_E_keV = model.reference_energy_keV,
        )
end;

Calibrated:

ref energy = 69.5 keV,
lac water = 0.19335 cm⁻¹.

2. FBP

Water BHC → FDK with the :standard kernel → HU, releasing the device buffers at the end.

hu_fbp = sim === nothing ? nothing : let
    sino_gpu = BS.apply_bhc_water(to_gpu(sim.sino), bhc_calibration.model)
    ws_fdk = BS.create_fdk_recon_workspace(sino_gpu, sim.geom, recon_opts.matrix_size)
    recon_μ = BS.reconstruct!(ws_fdk, sino_gpu, sim.geom)
    hu = Float32.(BS.to_hounsfield(Array(recon_μ); μ_water = bhc_calibration.μ_water))
    ws_fdk = nothing; sino_gpu = nothing; recon_μ = nothing
    GC.gc(true)
    hu
end;

3. Hybrid IR

The same chain with create_hir_recon_workspace(...; strength = 60) in place of the FDK workspace.

What Hybrid IR does

An FBP start, then ordered-subsets penalised weighted least squares with an edge-preserving Huber prior: the open-literature (Fessler) counterpart of the vendor hybrid IR algorithms (GE ASIR-V, Siemens SAFIRE, Philips iDose⁴, Canon AIDR 3D), whose noise-reduction range its strength table targets.

One dial: strength

A percentage in steps of 10, read like the GE ASIR-V dial: 0 is pure FBP (the iterative loop is skipped), 60 the standard clinical setting used here, 100 the maximum noise reduction. It moves the regularisation weight, the Huber threshold, the relaxation and the number of epochs together along one calibrated trajectory.

Projector consistency

Hybrid IR's data term uses a forward projector, which must be the one that generated the sinogram. The cell passes projector = sim_opts.projector, so changing the simulation's projector keeps the reconstruction consistent. FBP has no forward projector and is unaffected.

hu_hir = sim === nothing ? nothing : let
    sino_gpu = BS.apply_bhc_water(to_gpu(sim.sino), bhc_calibration.model)
    ws_hir = BS.create_hir_recon_workspace(
        sino_gpu, sim.geom, recon_opts.matrix_size; strength = 60, projector = sim_opts.projector,
    )
    recon_μ = BS.reconstruct!(ws_hir, sino_gpu, sim.geom)
    hu = Float32.(BS.to_hounsfield(Array(recon_μ); μ_water = bhc_calibration.μ_water))
    ws_hir = nothing; sino_gpu = nothing; recon_μ = nothing
    GC.gc(true)
    hu
end;

Results

Compare FBP vs Hybrid IR

Both reconstructions share the whole chain except the algorithm. The soft-tissue window shows the iodine-enhanced blood pools (the cardiac chambers) and the difference in noise texture between the two.

Measure noise inside one material, not inside a box

A box near the image centre spans lung to iodinated blood, so its standard deviation is anatomy, not noise, and would hide what HIR does. The cell below carries the phantom labels onto the reconstruction grid with BS.resample_to_recon and erodes the myocardium label, so only voxels surrounded by myocardium are measured.

Noise inside the myocardium (12336 voxels, all slices)

algorithm	mean (HU)	σ (HU)
FBP	47.8	18.4
Hybrid IR, strength 60	47.9	13.3

Hybrid IR lowers the noise by 27.4 % and leaves the mean where it was. The package's calibration band for strength 60 is 25–35 %.

Summary

Three ideas on top of the five-struct API of notebook 01:

Voxel phantoms from an artifact. BS.load_xcat_male_slab() (and the other load_xcat_* loaders) download a published XCAT phantom once, verify it, and return the labeled mask, a material per label and the voxel size; BS.Phantom(to_gpu(mask), materials, voxel_size) makes it simulatable. Any other labeled mask works the same way.
Materials you define. XA.Materials.ncat_* for the standard tissues, or XA.Material(name, ZA, I, density, composition) for a contrast bolus, a calibration solution or an alloy; the composition is a Dict of atomic number → mass fraction, and the loader's materials = keyword swaps it in by XCIST name.
Hybrid IR. create_hir_recon_workspace(sino, geom, matrix; strength) in place of the FDK workspace: the same reconstruct! call, lower noise at the same mean HU.