XCAT → CT: Phantom Grids and the Affine Round-Trip
Scan the heart of an ultra-high-resolution XCAT chest, then lay the ground-truth labels exactly onto the reconstructed image, axial and helical.
A simulation has two voxel grids. The phantom grid is your ground truth: here an XCAT 3.0 chest at 0.2 mm, read at 0.4 mm. The reconstruction grid is what the scanner outputs: a centred stack of axial slices at the pixel size you choose. Every quantitative use of a simulation (organ ROI statistics, segmentation scoring, partial-volume analysis) needs the map between the two. This notebook shows that map and proves that it is exact:
The reconstruction grid is always centred on the isocentre;
ReconOptionshas no off-centre FOV parameter. To scan a sub-region of a large phantom, crop the phantom so the region sits at the isocentre. Only the cropped block is ray-traced.BS.phantom_to_world_affine(phantom)andBS.recon_to_world_affine(geom, matrix_size)are the two grids' 4 × 4 voxel → world (cm) matrices.BS.resample_to_recon(phantom, geom, matrix_size; method = :nearest | :linear)carries the ground truth onto the reconstruction grid.
XCAT 3.0 chest export (0.2 mm) ─▶ read at 0.4 mm
─▶ find the cardiac labels by name ─▶ bounding box + 1 cm
─▶ crop the phantom to the box (the block is centred on the isocentre)
─▶ Scan A: axial, 14 cm FOV ─▶ water BHC ─▶ FDK ─▶ HU
─▶ Scan B: helical, 4 cm of z ─▶ water BHC ─▶ WFBP ─▶ HU
─▶ resample_to_recon ─▶ labels on the HU image, organ ROI statistics, exactness auditNotebook setup
import BasisSimulator as BSbegin
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)),)
endBackend detected: CuArray
Load the XCAT phantom
1. Locate the export
The phantom is an XCAT 3.0 chest export (adult male, 50th percentile, with coronary plaque): one 8-bit little-endian activity volume, *_act_1.raw, whose dimensions are in its file name, next to a Material_Spreadsheets/ folder with one material table per contrast state. The data is licensed and not in the repository. The notebook looks in BASISSIM_XCAT_DIR, or in docs/notebooks/data/xcat/ when that variable is unset, and expects exactly one *act_1.raw there. Without it, every compute cell skips and says so.
const XCAT_DIR = get(ENV, "BASISSIM_XCAT_DIR", joinpath(@__DIR__, "data", "xcat"));const CONTRAST_STATE = "high_contrast"; # the sheets: non_contrast, low_contrast, high_contrastfind_xcat_export(dir, state) -> NamedTuple or StringThe one *act_1.raw in dir, its dimensions read off its name (…_1600x1400x867_…), its anatomy (vmale_50) and the material sheet of contrast state, or a String saying what is missing.
xcat = find_xcat_export(XCAT_DIR, CONTRAST_STATE);const HAS_XCAT = xcat isa NamedTuple;XCAT export: vmale_50_raw_1600x1400x867_8bit_little_endian_act_1.raw (1600 × 1400 × 867 voxels, 1.81 GiB) · materials: vmale_50_materials_heart_high_contrast.xlsx
2. Read the mask at 0.4 mm
The export is 0.2 mm isotropic, label 0 is air, and the file runs anterior → posterior along its second axis. Reversing that axis makes y increase towards the patient's front, so every image drawn with y up shows the sternum at the top and the patient's left on the viewer's right. A label-preserving nearest-neighbour step of 2 then reads it at 0.4 mm: still finer than the reconstruction pixels below, so the resampling has something to interpolate across.
The activity volume, y towards anterior (the file runs posterior along axis 2).
"""The activity volume, `y` towards anterior (the file runs posterior along axis 2)."""
read_xcat_mask(raw, dims) = reverse(read!(raw, Array{UInt8}(undef, dims...)); dims = 2)Nearest-neighbour 3D downsample by an integer factor: keeps labels intact.
const DOWNSAMPLE_FACTOR = 2;const VOXEL_SIZE_CM = ntuple(_ -> 0.02 * DOWNSAMPLE_FACTOR, 3); # 0.2 mm export × factorphantom_full_uhr = HAS_XCAT ?
downsample_labeled(read_xcat_mask(xcat.raw, xcat.dims), DOWNSAMPLE_FACTOR) : nothing;Phantom read: 800 × 700 × 433 voxels of 0.4 mm (231.2 MiB as UInt8), 32.0 × 28.0 × 17.3 cm, 56 labels present.
3. Materials from the XCAT sheet
Each sheet row is one organ: a name, its elemental mass fractions (one column per atomic number), its density and its label (Organ ID). Every row becomes an XrayAttenuation.Material keyed by its label. Labels the sheet does not name stay air, which is what Phantom assigns to any label without a material. The high_contrast sheet puts iodinated blood in the cardiac chambers, the aorta and the coronaries.
We need the table before cropping: the organ names are how the cardiac labels are found.
⟨Z/A⟩ of a composition (mass fractions).
Bragg additivity for the mean excitation energy of a composition.
load_xcat_materials(sheet) -> Dict{Int, XA.Material}Name | <Z> … | Density | Organ ID → one material per label.
"""
load_xcat_materials(sheet) -> Dict{Int, XA.Material}
`Name | <Z> … | Density | Organ ID` → one material per label.
"""
function load_xcat_materials(sheet)
wb = XLSX.readxlsx(sheet)
data = wb[first(XLSX.sheetnames(wb))][:]
ncol = size(data, 2)
(string(data[1, ncol - 1]) == "Density" && string(data[1, ncol]) == "Organ ID") ||
error("unexpected sheet layout: $(join(string.(data[1, :]), " | "))")
zcols = [(c, parse(Int, string(data[1, c]))) for c in 2:(ncol - 2)]
out = Dict{Int, BS.XA.Material}()
for r in 2:size(data, 1)
data[r, ncol] isa Number || continue
comp = Dict{Int, Float64}(Z => Float64(data[r, c]) for (c, Z) in zcols
if data[r, c] isa Number && data[r, c] > 0)
isapprox(sum(values(comp)), 1; atol = 1.0e-3) || error("mass fractions of $(data[r, 1]) do not sum to 1")
out[Int(data[r, ncol])] = BS.XA.Material(
string(data[r, 1]), compute_ZA_ratio(comp), compute_mean_excitation_energy(comp),
Float64(data[r, ncol - 1]) * u"g/cm^3", comp,
)
end
return out
endmaterials_full = HAS_XCAT ? load_xcat_materials(xcat.sheet) : nothing;61 materials in vmale_50_materials_heart_high_contrast.xlsx; 56 labels present in the volume; labels present but not named by the sheet (read as air): none.
4. The cardiac bounding box
The heart is found by name: myocardium (myo…) and blood pools (bldpl…) of the four chambers, the pericardium, the coronary arteries and veins with their walls, and the coronary plaque components (rca…, lcx…, lad…). The aorta is left out: it runs the length of the chest and would stretch the box to the whole export. The bounding box of those labels, padded by 1 cm on every side, is the region we scan.
A heuristic over XCAT's naming
For a phantom whose labels are named differently, replace the pattern with an explicit list of label IDs.
heart_label_ids = materials_full === nothing ? nothing : let
pattern = r"^(myo|bldpl)|pericardium|coronary|^(rca|lcx|lad)\d"i
sort(UInt8[UInt8(id) for (id, mat) in materials_full if occursin(pattern, mat.name)])
end;34 cardiac labels: 15 myoLV · 16 myoRV · 17 myoLA · 18 myoRA · 19 bldplLV · 20 bldplRV · 21 bldplLA · 22 bldplRA · 26 coronary_artery · 27 coronary_vein · 29 pericardium · 30 coronary_artery_vessel_wall · 33 coronary_wall_intima · 34 coronary_wall_adventitia · 40 rca1_lipid · 41 rca2_lipid · 44 lad2_lipid · 45 lad3_lipid · 46 rca1_calc · 47 rca2_calc · 50 lad2_calc · 51 lad3_calc · 52 rca1_fib_mat · 53 rca2_fib_mat · 56 lad2_fib_mat · 57 lad3_fib_mat · 58 rca1_adventitia · 59 rca2_adventitia · 62 lad2_adventitia · 63 lad3_adventitia · 64 rca1_intima · 65 rca2_intima · 68 lad2_intima · 69 lad3_intima
heart_bbox = (phantom_full_uhr === nothing || heart_label_ids === nothing) ? nothing : let
is_heart = falses(256)
for id in heart_label_ids
is_heart[Int(id) + 1] = true
end
nx, ny, nz = size(phantom_full_uhr)
lo, hi = [nx + 1, ny + 1, nz + 1], [0, 0, 0]
@inbounds for k in 1:nz, j in 1:ny, i in 1:nx
if is_heart[Int(phantom_full_uhr[i, j, k]) + 1]
lo .= min.(lo, (i, j, k)); hi .= max.(hi, (i, j, k))
end
end
hi[1] == 0 && error("no cardiac voxels in the phantom")
pad = round.(Int, 1.0 ./ VOXEL_SIZE_CM) # 1 cm
lo .= max.(1, lo .- pad); hi .= min.((nx, ny, nz), hi .+ pad)
(i_lo = lo[1], i_hi = hi[1], j_lo = lo[2], j_hi = hi[2], k_lo = lo[3], k_hi = hi[3])
end;Scan A: axial, zoomed cardiac FOV
1. Crop the phantom
An index operation on the mask. The cropped block is what simulate! sees, so nothing outside it is ray-traced or held in device memory. It also changes the object: the rest of the chest is no longer in the beam, so this is a scan of the cardiac block, not of a chest with a small display field. For a chest scan with a small field, keep the whole phantom and set a small ReconOptions(fov_cm = …) instead.
phantom_cropped = heart_bbox === nothing ? nothing : let b = heart_bbox
phantom_full_uhr[b.i_lo:b.i_hi, b.j_lo:b.j_hi, b.k_lo:b.k_hi]
end;Crop: 373 × 331 × 342 voxels, 14.9 × 13.2 × 13.7 cm · 17.4 % of the 242.5 M-voxel phantom
2. Build the Phantom and its world affine
Without an origin, Phantom centres the block on the isocentre (origin = -extent/2 + voxel/2), which is exactly where the centred reconstruction grid will look.
Two phantoms are built from the same crop. phantom_cpu keeps the XCAT labels, for the resampling. The simulation phantom is compact_materials(phantom_cpu) on the device: the same attenuation with the labels renumbered densely, because the single-pass :dd_fast projector handles at most 64 materials and the sheet's label IDs run to 69.
phantom_cpu = phantom_cropped === nothing ? nothing :
BS.Phantom(phantom_cropped, materials_full, VOXEL_SIZE_CM);phantom = phantom_cpu === nothing ? nothing : let c = BS.compact_materials(phantom_cpu)
BS.Phantom(to_gpu(c.mask), c.materials, c.voxel_size, c.origin, c.extent)
end;phantom_to_world_affine
The 4 × 4 matrix maps a 0-indexed phantom voxel (i, j, k) to world coordinates (x, y, z) in cm:
[ x ] [ vx 0 0 ox ] [ i ]
[ y ] = [ 0 vy 0 oy ] · [ j ]
[ z ] [ 0 0 vz oz ] [ k ]
[ 1 ] [ 0 0 0 1 ] [ 1 ](vx, vy, vz) is the voxel size and (ox, oy, oz) the world position of voxel (0, 0, 0).
A_phantom = phantom_to_world_affine(phantom_cpu) (cm)
| col 1 | col 2 | col 3 | col 4 |
|---|---|---|---|
| 0.04 | 0.0 | 0.0 | -7.44 |
| 0.0 | 0.04 | 0.0 | -6.6 |
| 0.0 | 0.0 | 0.04 | -6.82 |
| 0.0 | 0.0 | 0.0 | 1.0 |
voxel = (0.4, 0.4, 0.4) mm, extent = (14.92, 13.24, 13.68) cm
origin (voxel
(0, 0, 0)) = (-7.44, -6.6, -6.82) cmcentre of the block = (-0.0, -0.0, -0.0) cm, the isocentre
3. Scanner, protocol and a tight reconstruction grid
The GE Revolution Apex Elite of notebook 01: 256 × 0.625 mm rows, a curved Lumex detector and the large-body bowtie. The grid mapping has nothing to do with the detector type, so any scanner would do.
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 = 3500.0, # e⁻ rms, added to the counts before the log
detection_gain = 10.0,
);CTProtocol: an axial cardiac CTA
120 kVp / 250 mA, a 1 s rotation of 500 views, 5 mm of collimation. SimOptions models the detector integrating each view while the gantry turns (view_samples = 5, notebook 01).
protocol = BS.CTProtocol(
kVp = 120,
mA = 250.0,
views = 500,
rotation_time = 1.0,
collimation_mm = 5.0,
additional_filters = [("Al", 4.5)],
);sim_opts = BS.SimOptions(seed = 1234, projector = :dd_fast, view_samples = 5);ReconOptions: a 14 cm field
A 14 cm × 14 cm field on a 384 × 384 grid (0.365 mm pixels), eight 0.625 mm slices over the 5 mm beam. The field is centred on the isocentre, where the cropped block sits, and it is deliberately smaller than the block in-plane: the reconstruction grid does not have to contain the object.
recon_opts = BS.ReconOptions(
matrix_size = (384, 384, round(Int, protocol.collimation_mm / 0.625)),
fov_cm = 14.0,
z_cm = protocol.collimation_mm / 10,
);recon_to_world_affine
The reconstruction grid can be inspected before anything is simulated: CTGeometry is built from the scanner, the protocol and the grid, exactly as the workspace builds it. The affine has the same form as the phantom's, with the reconstruction voxel size (fov / matrix_size) and a centred origin (-fov/2 + voxel/2).
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,
);A_recon = recon_to_world_affine(geom, matrix_size) (cm)
| col 1 | col 2 | col 3 | col 4 |
|---|---|---|---|
| 0.0365 | 0.0 | 0.0 | -6.9818 |
| 0.0 | 0.0365 | 0.0 | -6.9818 |
| 0.0 | 0.0 | 0.0625 | -0.2188 |
| 0.0 | 0.0 | 0.0 | 1.0 |
matrix = 384 × 384 × 8, voxel = 0.365 × 0.365 × 0.625 mm, FOV = 14.0 × 14.0 × 0.5 cm
origin = -6.9818, -6.9818, -0.2188 cm: centred on the isocentre
The two grids share world coordinates but little else:
| phantom (cropped) | reconstruction | |
|---|---|---|
| voxels | 373 × 331 × 342 | 384 × 384 × 8 |
| voxel size (mm) | 0.4 × 0.4 × 0.4 | 0.365 × 0.365 × 0.625 |
| extent (cm) | 14.92 × 13.24 × 13.68 | 14.0 × 14.0 × 0.5 |
| origin (cm) | -7.44, -6.6, -6.82 | -6.982, -6.982, -0.219 |
4. Simulate and reconstruct
The chain of notebook 01: create_workspace → simulate!, then the knobless water beam-hardening correction → FDK → HU with the correction's own μ_water.
sim = phantom === nothing ? nothing : let
ws = BS.create_workspace(scanner, protocol, sim_opts, recon_opts, phantom)
t = @elapsed result = BS.simulate!(ws, phantom, protocol, sim_opts)
out = (sino = Array(ws.sinogram), geom = ws.geom, dose = result.dose, t = t)
ws = nothing
GC.gc(true)
out
end;bhc = sim === nothing ? nothing : BS.calibrate_bhc_water(sim_opts, protocol; scanner, geom = sim.geom);recon_HU = sim === nothing ? nothing : let
sino = BS.apply_bhc_water(to_gpu(sim.sino), bhc)
ws = BS.create_fdk_recon_workspace(sino, sim.geom, recon_opts.matrix_size; filter = :standard)
μ = Array(BS.reconstruct!(ws, sino, sim.geom))
ws = nothing; sino = nothing; GC.gc(true)
Float32.(BS.to_hounsfield(μ; μ_water = bhc.μ_water_ref))
end;Scan A: 834 × 10 × 500 sinogram (columns × rows × views), simulated in 50.5 s (first call, compilation included) · CTDIvol = 17.15 mGy · water BHC at 69.5 keV
5. Resample the ground truth
resample_to_recon takes each reconstruction voxel's world coordinate from A_recon, maps it to a continuous phantom index through inv(A_phantom), and samples the phantom there:
method | output | use it for |
|---|---|---|
:nearest | UInt8 labels | ROI extraction, segmentation scoring |
:linear | Float32 (trilinear) | continuous fields, and the coverage fraction of one binary mask |
Trilinear interpolation of a multi-label mask averages label IDs, which means nothing. To get a partial-volume fraction, resample a 0/1 mask of the structure instead.
gt_resampled_nn = sim === nothing ? nothing :
BS.resample_to_recon(phantom_cpu, sim.geom, recon_opts.matrix_size; method = :nearest);Fraction of each reconstruction voxel covered by the cardiac labels (trilinear on a 0/1 mask).
"""Fraction of each reconstruction voxel covered by the cardiac labels (trilinear on a 0/1 mask)."""
function cardiac_coverage(mask, geom, matrix_size)
binary = zeros(UInt8, size(mask))
for id in heart_label_ids
binary[mask .== id] .= 0x01
end
ph = BS.Phantom(binary, Dict(0 => BS.XA.Materials.air, 1 => BS.XA.Materials.water), VOXEL_SIZE_CM)
return BS.resample_to_recon(ph, geom, matrix_size; method = :linear)
endcardiac_coverage_lin = sim === nothing ? nothing :
cardiac_coverage(phantom_cropped, sim.geom, recon_opts.matrix_size);gt_resampled_nn is UInt8 384 × 384 × 8, cardiac_coverage_lin is Float32 384 × 384 × 8, recon_HU is Float32 384 × 384 × 8: index (i, j, k) is the same physical voxel in all three.
Bring your own interpolator
For a B-spline, a sinc kernel or a learned upsampler, compose the two affines yourself: M = inv(A_phantom) * A_recon takes a 0-indexed reconstruction voxel (i, j, k, 1) to a continuous phantom index, which any interpolator (Interpolations.jl, ImageTransformations.jl, your own kernel) can sample:
for k in 0:(nz - 1), j in 0:(ny - 1), i in 0:(nx - 1)
p = M * [i, j, k, 1.0] # continuous phantom index
out[i + 1, j + 1, k + 1] = my_interpolator(phantom.mask, p[1], p[2], p[3])
endM = inv(A_phantom) * A_recon
| col 1 | col 2 | col 3 | col 4 |
|---|---|---|---|
| 0.9115 | 0.0 | 0.0 | 11.4557 |
| 0.0 | 0.9115 | 0.0 | -9.5443 |
| 0.0 | 0.0 | 1.5625 | 165.0312 |
| 0.0 | 0.0 | 0.0 | 1.0 |
The centre of the reconstruction grid maps to phantom index (186.0, 165.0, 170.5), the centre of the cropped block ((186.0, 165.0, 170.5)).
The overlay
Four views of the central slice, all on the reconstruction grid:
| panel | shows |
|---|---|
| top left | the HU reconstruction |
| top right | every label, resampled with :nearest |
| bottom left | the cardiac labels (:nearest) over the HU image |
| bottom right | the cardiac coverage fraction (:linear on a 0/1 mask) over the HU image |
Organ ROI statistics from the resampled labels
With the labels on the reconstruction grid, an organ ROI is one comparison: recon_HU[gt_resampled_nn .== label]. Below, each label with enough voxels on the central slice is eroded in-plane by one voxel (to drop the partial-volume rim) and its mean HU is compared with the monoenergetic HU of its sheet material at the BHC reference energy. Expect soft tissue, fat, lung and blood within a few tens of HU. Bone reads low, the single-kVp beam-hardening residual that notebook 01 quantifies, and small or thin structures pick up partial volume from their neighbours across the 0.625 mm slice.
Scan A, slice 5: per-label ROI statistics (theory = monoenergetic HU at 69.5 keV)
| label | material | voxels | mean HU | σ (HU) | theory HU | mean − theory |
|---|---|---|---|---|---|---|
| 0 | air | 960 | -974.2 | 29.3 | -998.9 | 24.7 |
| 1 | fat_background | 3818 | -111.4 | 64.6 | -104.6 | -6.8 |
| 3 | sternum | 447 | 1411.5 | 99.1 | 1572.6 | -161.1 |
| 4 | muscle_bone_marrow | 6124 | 51.9 | 68.5 | 43.1 | 8.8 |
| 6 | lung_parenchyma | 22920 | -724.9 | 26.0 | -741.2 | 16.3 |
| 11 | airway | 458 | -937.5 | 35.5 | -998.9 | 61.5 |
| 13 | esophagus | 687 | -3.1 | 51.7 | 43.1 | -46.2 |
| 15 | myoLV | 9403 | 77.3 | 24.1 | 64.1 | 13.3 |
| 16 | myoRV | 1679 | 72.4 | 17.5 | 64.1 | 8.3 |
| 17 | myoLA | 972 | 108.1 | 52.3 | 64.1 | 44.0 |
| 18 | myoRA | 928 | 68.4 | 19.9 | 64.1 | 4.3 |
| 19 | bldplLV | 11689 | 543.7 | 27.4 | 507.4 | 36.3 |
| 20 | bldplRV | 11914 | 159.7 | 10.5 | 140.8 | 18.9 |
| 21 | bldplLA | 8821 | 548.0 | 29.4 | 507.4 | 40.6 |
| 22 | bldplRA | 5643 | 163.2 | 9.3 | 140.8 | 22.4 |
| 24 | vein | 876 | 105.5 | 77.6 | 126.7 | -21.2 |
| 28 | aorta | 555 | 529.9 | 75.3 | 507.4 | 22.5 |
| 29 | pericardium | 10159 | -100.2 | 35.9 | -104.6 | 4.3 |
Scan B: helical, extended z
The same phantom and the same three calls, now with a taller crop scanned by a helical acquisition and reconstructed into a stack of 64 axial slices to scroll through. The helical reconstruction grid has its own affine (a longer z extent than Scan A), and resample_to_recon lands the labels on every slice unchanged.
1. A taller crop
The in-plane box of Scan A, extended 4 cm beyond the cardiac extent at each end in z (where the phantom allows), so the helix has anatomy to sweep through.
heart_bbox_tall = heart_bbox === nothing ? nothing : let b = heart_bbox
extra = round(Int, 4.0 / VOXEL_SIZE_CM[3])
merge(b, (k_lo = max(1, b.k_lo - extra), k_hi = min(size(phantom_full_uhr, 3), b.k_hi + extra)))
end;phantom_cropped_tall = heart_bbox_tall === nothing ? nothing : let b = heart_bbox_tall
phantom_full_uhr[b.i_lo:b.i_hi, b.j_lo:b.j_hi, b.k_lo:b.k_hi]
end;phantom_helical_cpu = phantom_cropped_tall === nothing ? nothing :
BS.Phantom(phantom_cropped_tall, materials_full, VOXEL_SIZE_CM);phantom_helical = phantom_helical_cpu === nothing ? nothing : let c = BS.compact_materials(phantom_helical_cpu)
BS.Phantom(to_gpu(c.mask), c.materials, c.voxel_size, c.origin, c.extent)
end;recon_opts_helical = BS.ReconOptions(
matrix_size = (384, 384, 64), # 64 slices of 0.625 mm
fov_cm = 14.0,
z_cm = 4.0,
);2. The helical acquisition
pitch and n_rotations make the protocol helical: 10 mm of collimation at pitch 1.0 over 8 rotations is 8 cm of table travel, centred on the isocentre. The geometry becomes a z-ramped trajectory, :dd_fast projects it unchanged (the sub-views of the view integration follow the helix, each advanced along z by its share of the table feed), and reconstruct! recognises the helical geometry and runs rebinned weighted FBP. The reconstruction grid is still a centred stack of axial slices, so the affines apply as before.
protocol_helical = BS.CTProtocol(
kVp = 120, mA = 250.0, views = 500, rotation_time = 1.0,
collimation_mm = 10.0, additional_filters = [("Al", 4.5)],
pitch = 1.0, n_rotations = 8,
);sim_helical = phantom_helical === nothing ? nothing : let
ws = BS.create_workspace(scanner, protocol_helical, sim_opts, recon_opts_helical, phantom_helical)
t = @elapsed result = BS.simulate!(ws, phantom_helical, protocol_helical, sim_opts)
out = (sino = Array(ws.sinogram), geom = ws.geom, dose = result.dose, t = t)
ws = nothing
GC.gc(true)
out
end;recon_HU_helical = sim_helical === nothing ? nothing : let
bhc_h = BS.calibrate_bhc_water(sim_opts, protocol_helical; scanner, geom = sim_helical.geom)
sino = BS.apply_bhc_water(to_gpu(sim_helical.sino), bhc_h)
ws = BS.create_fdk_recon_workspace(sino, sim_helical.geom, recon_opts_helical.matrix_size; filter = :standard)
μ = Array(BS.reconstruct!(ws, sino, sim_helical.geom)) # helical geometry → rebinned WFBP
ws = nothing; sino = nothing; GC.gc(true)
Float32.(BS.to_hounsfield(μ; μ_water = bhc_h.μ_water_ref))
end;Scan B: 834 × 16 × 4000 sinogram, simulated in 256.4 s · CTDIvol = 17.14 mGy, DLP = 137.1 mGy·cm over 8.0 cm (pitch 1.0)
A_recon of the helical grid (cm)
| col 1 | col 2 | col 3 | col 4 |
|---|---|---|---|
| 0.0365 | 0.0 | 0.0 | -6.9818 |
| 0.0 | 0.0365 | 0.0 | -6.9818 |
| 0.0 | 0.0 | 0.0625 | -1.9688 |
| 0.0 | 0.0 | 0.0 | 1.0 |
64 slices over 4.0 cm, against Scan A's 8 over 0.5 cm; the same centred origin rule.
gt_helical_nn = sim_helical === nothing ? nothing :
BS.resample_to_recon(phantom_helical_cpu, sim_helical.geom, recon_opts_helical.matrix_size; method = :nearest);cardiac_coverage_lin_helical = sim_helical === nothing ? nothing :
cardiac_coverage(phantom_cropped_tall, sim_helical.geom, recon_opts_helical.matrix_size);3. Scroll through z
Left: the WFBP reconstruction. Middle: the cardiac labels resampled with :nearest, which snap every reconstruction voxel to the single closest phantom voxel, so boundaries stair-step onto the coarser grid. Right: the cardiac coverage from :linear on the 0/1 mask, a fully 3D partial-volume fraction that also catches voxels straddling a surface in z.
Trilinear sampling at the voxel centre equals the true covered fraction when the two grids are of similar resolution, as here (0.4 mm phantom, 0.365 mm × 0.625 mm reconstruction). For a phantom much finer than the reconstruction, a true volume fraction needs box averaging instead.
4. The mapping is exact
Three checks, for the axial and the helical geometry:
The affine is the reconstructor's grid. The FDK and WFBP back-projectors place reconstruction voxel
idx(1-indexed) at world-fov/2 + (idx − ½)·fov/n, the rulerecon_to_world_affineencodes. The table measures the largest difference.The round trip is the identity.
A⁻¹ · A · v = vto floating-point precision.The image is registered to the map. Checks 1 and 2 certify the coordinate map; this one certifies the image. Edges of the reconstruction are correlated with the label boundaries over integer in-plane shifts: the best shift should be
(0, 0). A half-voxel or rebinning offset between the forward projector and the back-projector would show up here and nowhere else.
is_helical(geom) switches the back-projection algorithm, never the grid, so the helical mapping is as exact as the axial one. Softness at a boundary is :nearest quantisation or the reconstruction's point-spread function, not the mapping.
| geometry | affine vs reconstructor grid (max, µm) | round trip A⁻¹·A·v − v (max, voxels) |
|---|---|---|
| axial | 8.9e-12 | 0.0 |
| helical | 8.9e-12 | 0.0 |
Registration of the image to the labels (central slice, shifts in voxels)
| geometry | best (dx, dy) | score at (0, 0) / best |
|---|---|---|
| axial | (0, 0) | 1.0 |
| helical | (0, 0) | 1.0 |
Summary
Crop the phantom, not the reconstruction. The reconstruction grid is always centred on the isocentre; cropping the input brings the region of interest there, and only the cropped block is ray-traced.
Two affines and one resampler.
phantom_to_world_affineandrecon_to_world_affinemap each grid's voxels to world centimetres;resample_to_reconcomposes them, with:nearestfor labels and:linearfor continuous fields or the coverage of a binary mask, andinv(A_phantom) * A_reconhands the map to any other interpolator.Exact, axial and helical. The affine matches the back-projectors' voxel rule, the round trip is the identity, and the reconstructed edges sit on the label boundaries with no shift.
What it buys you. Organ ROI statistics, segmentation scores and partial-volume fractions are one indexing expression on the reconstruction grid, as the per-label table shows.