← Examples/02 · XCAT Phantom + 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 slab (1385 × 917 × 1, 0.25 mm in-plane, full resolution)
Materials	Mix of `XA.Materials.ncat_*` prebuilt + a custom iodinated-blood material constructed inline
Scanner	GE Revolution Apex Elite (same as notebook 01)
Recon	FBP/FDK + Hybrid IR (`:asir`-style PWLS refinement)

Notebook Setup

Same shape as notebook 01 — activate docs/Project.toml, four discrete import cells, then probe for a GPU backend. We add Unitful here for the u"eV" and u"g/cm^3" units used by the XA.Material constructor.

begin
    import Pkg
    Pkg.activate(joinpath(@__DIR__, ".."))
end

import BasisSimulator as BS

# import CairoMakie as Mke
import WasmMakie as Mke

using Statistics: std, mean

using Markdown: @md_str

using Unitful: @u_str

import Unitful: ustrip, uconvert

Backend-agnostic device transfer: `to_gpu()`

Same probe-then-load pattern as notebook 01 — checks Metal → CUDA → AMDGPU via Base.locate_package, falls back to CPU identity.

name

"MtlArray"

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: MtlArray

Phantom Construction

01. 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 standard artifact store (~/.julia/artifacts/): it downloads once (~65 MB, SHA-256-verified), logs the Segars XCAT + XCIST citation, and reuses the cached copy on every later call.

Using `load_xcat_*`

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 = "/my/local/xcat" — skip the download and read an already-extracted copy from disk.
quiet = true — silence the citation + progress logging.

02. 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⁻³
`composition`	10 elements (Z = 1678111516171926)

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.

build_iodine_blood (generic function with 1 method)

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

name

"Iodine-blood (5.0 mg/mL)"

ZA_ratio

0.54999

75.2 eV

density

1.06 g cm⁻³

composition

("16", MIME type text/plain)

0.00199057

("8", MIME type text/plain)

0.741486

("17", MIME type text/plain)

0.00298585

("1", MIME type text/plain)

0.101519

("53", MIME type text/plain)

0.00471698

("19", MIME type text/plain)

0.00199057

("6", MIME type text/plain)

0.109481

("11", MIME type text/plain)

0.000995283

("7", MIME type text/plain)

0.0328443

("15", MIME type text/plain)

0.000995283

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.

("\"ncat_spleen\"", MIME type text/plain)

name

"NCAT Spleen"

ZA_ratio

0.55024

72.5 eV

density

1.06 g cm⁻³

composition

("6", MIME type text/plain)

0.113

("7", MIME type text/plain)

0.032

("11", MIME type text/plain)

0.001

("15", MIME type text/plain)

0.003

("16", MIME type text/plain)

0.002

("8", MIME type text/plain)

0.741

("\"ncat_lung\"", MIME type text/plain)

name

"NCAT Lung"

ZA_ratio

0.55048

75.3 eV

density

0.26 g cm⁻³

composition

("6", MIME type text/plain)

0.105

("7", MIME type text/plain)

0.031

("11", MIME type text/plain)

0.002

("15", MIME type text/plain)

0.002

("16", MIME type text/plain)

0.003

("8", MIME type text/plain)

0.749

("\"ncat_dry_rib\"", MIME type text/plain)

name

"NCAT Rib"

ZA_ratio

0.519

100.0 eV

density

1.52 g cm⁻³

composition

("6", MIME type text/plain)

0.235

("7", MIME type text/plain)

0.0409

("15", MIME type text/plain)

0.0702

("16", MIME type text/plain)

0.0022

("20", MIME type text/plain)

0.1348

("12", MIME type text/plain)

0.0011

("\"ncat_adipose\"", MIME type text/plain)

name

"NCAT Fat/Adipose"

ZA_ratio

0.55579

63.2 eV

density

0.92 g cm⁻³

composition

("16", MIME type text/plain)

0.00073

("7", MIME type text/plain)

0.00797

("20", MIME type text/plain)

2.0e-5

("30", MIME type text/plain)

2.0e-5

("8", MIME type text/plain)

0.232333

("17", MIME type text/plain)

0.00119

("\"ncat_cartilage\"", MIME type text/plain)

name

"Cartilage (ICRU-44)"

ZA_ratio

0.54595

72.6 eV

density

1.1 g cm⁻³

composition

("6", MIME type text/plain)

0.099

("7", MIME type text/plain)

0.022

("11", MIME type text/plain)

0.005

("15", MIME type text/plain)

0.022

("16", MIME type text/plain)

0.009

("8", MIME type text/plain)

0.744

("\"ncat_water\"", MIME type text/plain)

name

"Water, Liquid"

ZA_ratio

0.55508

75.0 eV

density

1.0 g cm⁻³

composition

("8", MIME type text/plain)

0.888102

("1", MIME type text/plain)

0.111898

("\"ncat_liver\"", MIME type text/plain)

name

"NCAT Liver"

ZA_ratio

0.55052

72.5 eV

density

1.06 g cm⁻³

composition

("6", MIME type text/plain)

0.139

("7", MIME type text/plain)

0.03

("11", MIME type text/plain)

0.002

("15", MIME type text/plain)

0.003

("16", MIME type text/plain)

0.003

("8", MIME type text/plain)

0.716

("\"ncat_dry_spine\"", MIME type text/plain)

name

"NCAT Spine"

ZA_ratio

0.525

95.0 eV

density

1.42 g cm⁻³

composition

("6", MIME type text/plain)

0.278

("7", MIME type text/plain)

0.0399

("15", MIME type text/plain)

0.0616

("16", MIME type text/plain)

0.0017

("20", MIME type text/plain)

0.0799

("12", MIME type text/plain)

0.0011

("\"ncat_intestine\"", MIME type text/plain)

name

"Intestine (ICRU-44)"

ZA_ratio

0.55017

72.5 eV

density

1.03 g cm⁻³

composition

("6", MIME type text/plain)

0.115

("7", MIME type text/plain)

0.022

("11", MIME type text/plain)

0.001

("15", MIME type text/plain)

0.001

("16", MIME type text/plain)

0.001

("8", MIME type text/plain)

0.751

("\"ncat_muscle\"", MIME type text/plain)

name

"NCAT Muscle"

ZA_ratio

0.54938

75.3 eV

density

1.05 g cm⁻³

composition

("6", MIME type text/plain)

0.143

… 6 more

… 2 more

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

03. 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 one 5 mm-thick axial section at full 0.25 mm in-plane resolution: real anatomy at a docs-friendly compute cost, no downsampling needed.

mask

1385×917×1 Array{UInt8, 3}: [:, :, 1] = 0x00 0x00 0x00 0x00 0x00 0x00 0x00 … 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 … 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 ⋮ ⋮ ⋱ ⋮ 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 … 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00

materials

("0", MIME type text/plain)

name

"Air, Dry (near sea level)"

ZA_ratio

0.49919

85.7 eV

density

0.001205 g cm⁻³

composition

("4", MIME type text/plain)

name

"NCAT Spine"

ZA_ratio

0.525

95.0 eV

density

1.42 g cm⁻³

composition

("5", MIME type text/plain)

name

"NCAT Rib"

ZA_ratio

0.519

100.0 eV

density

1.52 g cm⁻³

composition

("6", MIME type text/plain)

name

"NCAT Fat/Adipose"

ZA_ratio

0.55579

63.2 eV

density

0.92 g cm⁻³

composition

("2", MIME type text/plain)

name

"NCAT Muscle"

ZA_ratio

0.54938

75.3 eV

density

1.05 g cm⁻³

composition

("7", MIME type text/plain)

name

"Iodine-blood (5.0 mg/mL)"

ZA_ratio

0.54999

75.2 eV

density

1.06 g cm⁻³

composition

("9", MIME type text/plain)

name

"Cartilage (ICRU-44)"

ZA_ratio

0.54595

72.6 eV

density

1.1 g cm⁻³

composition

("8", MIME type text/plain)

name

"NCAT Heart"

ZA_ratio

0.55021

72.6 eV

density

1.05 g cm⁻³

composition

("3", MIME type text/plain)

name

"NCAT Lung"

ZA_ratio

0.55048

75.3 eV

density

0.26 g cm⁻³

composition

("1", MIME type text/plain)

name

"Water, Liquid"

ZA_ratio

0.55508

75.0 eV

density

1.0 g cm⁻³

composition

label_names

("5", MIME type text/plain)

"ncat_dry_rib"

("4", MIME type text/plain)

"ncat_dry_spine"

("6", MIME type text/plain)

"ncat_adipose"

("7", MIME type text/plain)

"ncat_blood"

("2", MIME type text/plain)

"ncat_muscle"

("9", MIME type text/plain)

"ncat_cartilage"

("8", MIME type text/plain)

"ncat_heart"

("3", MIME type text/plain)

"ncat_lung"

("1", MIME type text/plain)

"ncat_water"

voxel_size_cm

0.025

5.0

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

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

0.025

5.0

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

mask

1385×917×1 Metal.MtlArray{UInt8, 3, Metal.PrivateStorage}: [:, :, 1] = 0x00 0x00 0x00 0x00 0x00 0x00 0x00 … 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 … 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 ⋮ ⋮ ⋱ ⋮ 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 … 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00

materials

name

"Air, Dry (near sea level)"

ZA_ratio

0.49919

85.7 eV

density

0.001205 g cm⁻³

composition

name

"Water, Liquid"

ZA_ratio

0.55508

75.0 eV

density

1.0 g cm⁻³

composition

name

"NCAT Muscle"

ZA_ratio

0.54938

75.3 eV

density

1.05 g cm⁻³

composition

name

"NCAT Lung"

ZA_ratio

0.55048

75.3 eV

density

0.26 g cm⁻³

composition

name

"NCAT Spine"

ZA_ratio

0.525

95.0 eV

density

1.42 g cm⁻³

composition

name

"NCAT Rib"

ZA_ratio

0.519

100.0 eV

density

1.52 g cm⁻³

composition

name

"NCAT Fat/Adipose"

ZA_ratio

0.55579

63.2 eV

density

0.92 g cm⁻³

composition

name

"Iodine-blood (5.0 mg/mL)"

ZA_ratio

0.54999

75.2 eV

density

1.06 g cm⁻³

composition

name

"NCAT Heart"

ZA_ratio

0.55021

72.6 eV

density

1.05 g cm⁻³

composition

name

"Cartilage (ICRU-44)"

ZA_ratio

0.54595

72.6 eV

density

1.1 g cm⁻³

composition

voxel_size

0.025

5.0

origin

-17.3

-11.45

0.0

extent

34.625

22.925

5.0

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

Dict{Symbol, Any}(:msg => "UndefVarError: `DataAspect` not defined in `WasmMakie`\nSuggestion: check for spelling errors or missing imports.", :stacktrace => Dict{Symbol, Any}[Dict(:call_short => "getproperty", :inlined => true, :url => nothing, :path => "./Base_compiler.jl", :source_package => nothing, :call => "getproperty", :linfo_type => "Nothing", :line => 47, :file => "Base_compiler.jl", :func => "getproperty", :parent_module => nothing, :from_c => false), Dict(:call_short => "macro expansion", :inlined => true, :url => nothing, :path => "/Users/daleblack/Documents/dev/MolloiLab/BasisSimulator.jl/docs/notebooks/02_xcat_custom_materials.jl#==#04000005-0000-4000-8000-000000000001", :source_package => nothing, :call => "macro expansion", :linfo_type => "Nothing", :line => 12, :file => "02_xcat_custom_materials.jl#==#04000005-0000-4000-8000-000000000001", :func => "macro expansion", :parent_module => nothing, :from_c => false)], :plain_error => "UndefVarError: `DataAspect` not defined in `WasmMakie`\nSuggestion: check for spelling errors or missing imports.\nStacktrace:\n [1] getproperty\n @ ./Base_compiler.jl:47 [inlined]\n [2] macro expansion\n @ ~/Documents/dev/MolloiLab/BasisSimulator.jl/docs/notebooks/02_xcat_custom_materials.jl#==#04000005-0000-4000-8000-000000000001:12 [inlined]")

Scan Setup & Simulation

01. Scanner, protocol, sim & recon options

Same GE Revolution Apex Elite hardware as notebook 01. The protocol is a clinical body-CTA acquisition: 120 kVp / 250 mA, 5 mm collimation, 500 views. Recon: 512 × 512, 35 cm FOV, FBP :standard filter.

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

detector_row_offset

0.0

detector_col_offset

0.25

focal_spot_width

1.0

focal_spot_length

1.0

target_angle

10.0

gantry_rotation_time

0.5

scan_diameter

500.0

gantry_aperture

700.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

detection_gain

10.0

electronic_noise

3500.0

detector_type

:energy_integrating

n_energy_bins

energy_thresholds

energy_resolution

0.0

charge_sharing_fwhm

0.0

dead_time_ns

0.0

pixel_mode

:standard

native_dexel_col_mm

1.0549872122762147

native_dexel_row_mm

1.0989450127877238

binning_factor

detector_shape

:arc

scanner = BS.Scanner(
    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,
    # DAS/electronic noise, in electrons.  This enters the COUNTS before the log
    # transform (`λ_noisy += σ_e·randn()` in simulate!), which is where a real
    # DAS contributes it — so it propagates through reconstruction and iterative
    # recon can suppress it.  Adding an equivalent floor to the HU volume *after*
    # reconstruction instead would make FBP and HIR look identical by fiat.
    electronic_noise = 3500.0,   # e⁻ — clinical GE Apex Elite DAS readout noise
    detection_gain = 10.0,
)

250.0

kVp

120.0

views

500

rotation_time

1.0

n_rotations

1.0

collimation_mm

5.0

anode_angle

additional_filters

"Al"

4.5

pitch

nothing

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

use_fill_factor

true

use_detector_efficiency

true

use_scatter

true

use_optical_crosstalk

false

use_focal_spot

true

use_noise

true

use_lag

true

use_heel_effect

true

use_pcct_pileup

false

use_pcct_pileup_correction

false

use_pcct_scatter

false

use_pcct_scatter_correction

false

pcct_noise_reduction

0.0

seed

1234

detector_efficiency_mode

:auto

projector

:dd_fast

# `projector` picks the forward ray tracer: :dd_fast (default) is distance-driven and
# anti-aliased (robust in severe beam-hardened regions), fusing the whole spectrum into
# one volume walk. :siddon point-samples one voxel per step and ALIASES in those regions —
# use it only when speed outranks accuracy. (:dd is the older, slower reference kernel:
# same physics as :dd_fast, now deprecated.)
# The BHC and Hybrid-IR cells below all read `sim_opts.projector`, so changing it
# HERE updates the whole pipeline consistently — the recon must invert the operator
# that generated the data, or the IR system matrix won't match and convergence suffers.
sim_opts = BS.SimOptions(fidelity = :eict, seed = 1234, projector = :dd_fast)

matrix_size

512

fov_cm

35.0

z_cm

0.5

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

02. Forward project

Same let ... end shape as notebook 01 — workspace, simulate!, copy off GPU, drop refs, force GC.gc(true). The XCAT phantom is ~3× the volume of the Gammex 472 demo (320×280×100 vs 512×512×16), so this step takes longer; on a CPU-only fallback it can be several minutes.

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

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

        ws = nothing
        GC.gc(true)

        result
end;

Reconstruction

01. Postprocessing (The Full Correction Pipeline)

A raw FDK or HIR recon doesn't ship as a clinical image. Real CT vendors apply a stack of corrections after the simulator's forward model:

Stage	Function	What it does
1.	`calibrate_bhc_water`	Precomputes per-column water poly→mono polynomials from the FULL detected spectrum — pure physics, zero tunables
2.	`apply_bhc_water`	Sinogram-domain water BHC — one polynomial pass before FDK/HIR
3.	`reconstruct!` (FDK or Hybrid IR)	Filtered back-projection / iterative reconstruction on the corrected sinogram
5.	`to_hounsfield` (with BHC μ_water)	Convert μ → HU using the BHC model's calibrated reference (not the 70 keV NIST value)
6.	`measure_radial_cupping`	QA metric (never applied): fitted residual cup + DC — both ≈ 0 after a correct BHC

Where the DAS noise goes

Quantum (Poisson) and DAS/electronic noise are both injected by simulate! in the counts domain, before the log transform — electronic noise via scanner.electronic_noise (electrons). That is where a real DAS contributes it, and it means reconstruction sees it and can suppress it.

Do not substitute add_system_noise_floor! here. That helper adds white Gaussian noise to the HU volume after reconstruction, so it is identical in FBP and HIR by construction and makes iterative recon look useless — an earlier version of this notebook set electronic_noise = 0 and added a 28 HU floor post-recon, which reported HIR at ~1.6 % noise reduction instead of the ~26 % it actually delivers. Reserve that helper for effects that genuinely survive reconstruction (calibration drift, ring-correction residuals).

The same pipeline applies to both FBP and HIR — only the recon workspace inside the let block differs (create_fdk_recon_workspace vs create_hir_recon_workspace). See notebook 01 §9 for the same pattern on the Gammex 472 phantom.

a. Calibrate the BHC model

One-time spectrum fit at 120 kVp + 7 mm Al filtration. Returns a TwoMaterialBHC polynomial model + the calibrated μ_water_ref used by both FBP and HIR pipelines below for to_hounsfield.

bhc_calibration = sim === nothing ? nothing : let
        # KNOBLESS water BHC — pure physics, zero tunables, zero 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.7 keV,
lac water = 0.19321 cm⁻¹.

b. Polychromatic `μ_water` from the XCAT body — informational

The BHC's μ_water_ref above (the monoenergetic μwater at the spectrum-mean energy) is what §8 + §9 use as the HU divisor — that's the right reference post-BHC, where the recon reads as approximately monochromatic at `refE_keV`.

For comparison, here's what the polychromatic-effective μwater would be — the value an uncorrected polychromatic FBP would land on for solid water at this body chord. Useful as a sanity check; not plumbed into `tohounsfield`.

BS.compute_polychromatic_μ_water does the spectrum + Beer- Lambert hardening analytically; BS.estimate_phantom_diameter_cm reads the body's chord length straight off the XCAT mask, so the calibration scales naturally if you swap phantoms (vmale50 → vfemale50, etc.) without hardcoding any cm.

μ_water_poly_uncorrected = phantom === nothing ? nothing : let
        voxel_size_mm = VOXEL_SIZE_CM .* 10.0
        body_diameter_cm = BS.estimate_phantom_diameter_cm(
            phantom_labeled, voxel_size_mm
        )
        BS.compute_polychromatic_μ_water(
            sim_opts, protocol;
            scanner = scanner,
            geom = sim.geom,
            water_path_cm = body_diameter_cm,
        )
end;

Analytic poly μ_water (XCAT body chord, no-BHC reference for comparison):

μ_water_poly_uncorrected = $(μ_water_poly_uncorrected === nothing ? "—" : "$(round(μwaterpoly_uncorrected, digits = 5)) cm⁻¹")
bhc_calibration.μ_water = $(bhc_calibration === nothing ? "—" : "$(round(bhccalibration.μwater, digits = 5)) cm⁻¹") (used downstream)

02. FBP with the full correction pipeline

Same let ... end shape as notebook 01 — detected-spectrum water BHC → FDK → HU, with explicit GPU cleanup and non-mutating cupping QA at the end. HU conversion uses bhc_calibration.μ_water (the BHC's calibrated μ_water_ref at the spectrum-mean energy), since the post-BHC recon reads as approximately monochromatic at that reference.

hu_fbp = sim === nothing ? nothing : let
        matrix_size = recon_opts.matrix_size

        # 1. Sinogram-domain BHC (returns a CPU array)
        sino_gpu = to_gpu(sim.sino)
        sino_bhc = BS.apply_bhc_water(sino_gpu, bhc_calibration.model)
        sino_gpu = to_gpu(sino_bhc)

        # 2. FDK on the corrected sinogram
        ws_fdk = BS.create_fdk_recon_workspace(sino_gpu, sim.geom, matrix_size)
        recon_μ = BS.reconstruct!(ws_fdk, sino_gpu, sim.geom)

    # (image-domain BHC removed — audit: scaled self-subtraction that
    #  deflated dense-material HU; the knobless sinogram water BHC is the
    #  whole correction)

        # 4. μ → HU using BHC's calibrated μ_water_ref
        hu = Float32.(BS.to_hounsfield(Array(recon_μ); μ_water = bhc_calibration.μ_water))

        # 5. DAS/electronic noise is injected in the counts domain by simulate!
        # (scanner.electronic_noise), not bolted onto the HU volume here.

        # 6. Residual radial cupping is measured as QA, never applied
        # (cupping is QA-only now — measure_radial_cupping; never applied)

        # GPU cleanup
        ws_fdk = nothing
        sino_gpu = nothing
        recon_μ = nothing
        GC.gc(true)

        hu
end;

03. Hybrid IR with the full correction pipeline

Identical pipeline to §8 — sino BHC → recon → image BHC → HU → noise floor → cupping — with create_hir_recon_workspace(...; strength = 60) swapped in for the FDK workspace.

What HIR adds over FBP

Hybrid IR is the clinical workhorse: vendors call it ASIR, AIDR3D, iDose⁴, SAFIRE depending on brand. All share the same architecture — FDK init + iterative PWLS refinement with a Huber prior. It costs a few FDK passes but returns ~30 % lower pixel σ at matched resolution.

One dial: `strength`

strength is a percentage in 10 % steps, read exactly like the GE ASIR-V dial: 0 is pure FBP (the PWLS loop is skipped entirely), 60 is the standard clinical setting used here, 100 is maximum noise reduction. It moves λ, the Huber edge threshold δ, the relaxation and the epoch count together along one tuned trajectory — there is no second knob to keep in sync.

Projector consistency (IR only)

Hybrid IR's data-fidelity term uses a forward projector A, so it must use the same projector that generated the sinogram. We pass projector = sim_opts.projector to create_hir_recon_workspace, so flipping sim_opts.projector (e.g. to :siddon) automatically keeps the recon consistent. FBP (§8) is immune — it has no forward A.

hu_hir = sim === nothing ? nothing : let
        matrix_size = recon_opts.matrix_size

        # 1. Sinogram-domain BHC
        sino_gpu = to_gpu(sim.sino)
        sino_bhc = BS.apply_bhc_water(sino_gpu, bhc_calibration.model)
        sino_gpu = to_gpu(sino_bhc)

        # 2. Hybrid IR (FBP init + PWLS refinement, strength = 60 %)
        # projector MUST match the forward sim so the IR system matrix A inverts
        # the operator that generated the data (read from sim_opts, single source).
        ws_hir = BS.create_hir_recon_workspace(
            sino_gpu, sim.geom, matrix_size; strength = 60, projector = sim_opts.projector,
        )
        recon_μ = BS.reconstruct!(ws_hir, sino_gpu, sim.geom)

        # FBP's `reconstruct!` masks voxels outside the inscribed scan FOV
        # (`apply_fov_mask!`); HIR doesn't, so the iterative refinement
        # leaves garbage in the corners.  Apply the same mask for parity.
        BS.apply_fov_mask!(recon_μ, sim.geom)
    
        # 3. μ → HU using BHC's calibrated μ_water_ref
        hu = Float32.(BS.to_hounsfield(Array(recon_μ); μ_water = bhc_calibration.μ_water))

        # GPU cleanup
        ws_hir = nothing
        sino_gpu = nothing
        recon_μ = nothing
        GC.gc(true)

        hu
end;

Results

Compare FBP vs Hybrid IR

Both reconstructions go through the identical correction pipeline (water sino-BHC + recon + HU; counts-domain noise; cupping QA). The soft-tissue window shows the iodine-enhanced blood pools (XCAT labels 19–22 = bldplLV/RV/LA/RA) and the texture / noise contrast between the two algorithms.

Measure noise inside a material, not inside a box

A fixed box near the image center spans −794 to +353 HU here — lung through iodinated blood — so its σ (≈160 HU) is anatomy, not noise, and HIR appears to do nothing. The cell below instead uses BS.resample_to_recon to carry the phantom labels onto the recon grid, then erodes the myocardium label so only voxels surrounded by myocardium survive. That is the only place where σ means noise.

Dict{Symbol, Any}(:msg => "UndefVarError: `DataAspect` not defined in `WasmMakie`\nSuggestion: check for spelling errors or missing imports.", :stacktrace => Dict{Symbol, Any}[Dict(:call_short => "getproperty", :inlined => true, :url => nothing, :path => "./Base_compiler.jl", :source_package => nothing, :call => "getproperty", :linfo_type => "Nothing", :line => 47, :file => "Base_compiler.jl", :func => "getproperty", :parent_module => nothing, :from_c => false), Dict(:call_short => "macro expansion", :inlined => true, :url => nothing, :path => "/Users/daleblack/Documents/dev/MolloiLab/BasisSimulator.jl/docs/notebooks/02_xcat_custom_materials.jl#==#11000002-0000-4000-8000-000000000001", :source_package => nothing, :call => "macro expansion", :linfo_type => "Nothing", :line => 15, :file => "02_xcat_custom_materials.jl#==#11000002-0000-4000-8000-000000000001", :func => "macro expansion", :parent_module => nothing, :from_c => false)], :plain_error => "UndefVarError: `DataAspect` not defined in `WasmMakie`\nSuggestion: check for spelling errors or missing imports.\nStacktrace:\n [1] getproperty\n @ ./Base_compiler.jl:47 [inlined]\n [2] macro expansion\n @ ~/Documents/dev/MolloiLab/BasisSimulator.jl/docs/notebooks/02_xcat_custom_materials.jl#==#11000002-0000-4000-8000-000000000001:15 [inlined]")

Noise (σ in HU, 12336 voxels deep inside the myocardium, mean 47.8 HU):

Algorithm	σ (HU)
FBP	21.8
Hybrid IR	19.3
Reduction	11.3%

Compare against the strength = 60 target band (25–35 %). Both quantum and DAS/electronic noise enter in the counts domain, so HIR acts on all of it — no post-hoc correction needed to read this number.

Summary

This notebook layered three new ideas on top of the five-struct API walked in notebook 01:

External voxel phantoms — load any labeled .bin mask (XCAT, ICRP, custom), then hand it to BS.Phantom(...) along with a materials Dict. No internal preprocessing; whatever labeling you give, you get back.
Custom materials via XrayAttenuation — XA.Materials.ncat_* for canonical tissues, plus XA.Material(name, ZA, I, density, composition) built inline whenever you need a non-standard mixture (contrast bolus, calibration solution, alloy). Composition is just a Dict{Int,Float64} of atomic-number → mass-fraction.
Hybrid IR — BS.create_hir_recon_workspace(sino, geom, matrix; strength) swaps in for the FDK workspace and runs PWLS refinement after the FBP init. Same reconstruct! call site, lower pixel noise.

Every other piece — the let ... end GPU pattern, μ → HU conversion, the correction pipeline from §9 of notebook 01 — carries over unchanged.

XCAT Phantom + Custom Materials

Notebook Setup

Backend-agnostic device transfer: to_gpu()

Phantom Construction

01. Download the XCAT phantom

02. Custom materials via XrayAttenuation

a. Pull a prebuilt NCAT material

b. Construct a custom contrast material

c. The canonical name → material map

d. Override a material

03. Build the Phantom

Scan Setup & Simulation

01. Scanner, protocol, sim & recon options

02. Forward project

Reconstruction

01. Postprocessing (The Full Correction Pipeline)

a. Calibrate the BHC model

b. Polychromatic μ_water from the XCAT body — informational

02. FBP with the full correction pipeline

03. Hybrid IR with the full correction pipeline

Results

Compare FBP vs Hybrid IR

Summary

Backend-agnostic device transfer: `to_gpu()`

02. Custom materials via `XrayAttenuation`

03. Build the `Phantom`

b. Polychromatic `μ_water` from the XCAT body — informational