CBCT FDK

This page describes image reconstruction for cone-beam computed tomography (CBCT) with both arc and flat detectors using the Julia package Sinograms.jl.

This page comes from a single Julia file: 07-fdk.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: 07-fdk.ipynb, or open it in binder here: 07-fdk.ipynb.

Setup

Packages needed here.

using Plots: plot, gui
using Unitful: cm
using Sinograms: CtFanArc, CtFanFlat # CtPar
using Sinograms: rays, plan_fbp, Window, Hamming, fdk, ct_geom_plot3
using ImageGeoms: ImageGeom, MaskCircle, fovs
using ImagePhantoms: ellipsoid_parameters, ellipsoid
using ImagePhantoms: radon, phantom
using MIRTjim: jim, prompt

The following line is helpful when running this file as a script; this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

CBCT projections of 3D Shepp-Logan phantom

For illustration, we start by synthesizing CBCT projections of the 3D Shepp-Logan phantom.

For completeness, we use units (from Unitful), but units are optional.

Use ImageGeom to define the image geometry.

ig = ImageGeom(MaskCircle(); dims=(64,62,30), deltas = (1,1,2) .* 0.4cm)

ImageGeoms.ImageGeom{3, NTuple{3,Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}}, BitArray{3}}
 dims::NTuple{3,Int64} (64, 62, 30)
 deltas::NTuple{3,Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}} (0.4 cm, 0.4 cm, 0.8 cm)
 offsets::NTuple{3,Float32} (0.0f0, 0.0f0, 0.0f0)
 mask: 64×62×30 BitArray{3} {84840 of 119040}

Ellipsoid parameters for 3D Shepp-Logan phantom:

μ = 0.1 / cm # typical linear attenuation coefficient
params = ellipsoid_parameters( ; fovs = fovs(ig), u = (1, 1, μ))
ob = ellipsoid(params)

12-element Vector{ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}}:
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((0.0 cm, 0.0 cm, 0.0 cm), (8.831999999999999 cm, 11.408000000000001 cm, 10.8 cm), (0.0, 0.0, 0.0), 0.2 cm^-1, (0.0, 0.0, 0.0), (1.0, 1.0, 1.0))
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((0.0 cm, -0.22816 cm, 0.0 cm), (8.478720000000001 cm, 10.8376 cm, 10.56 cm), (0.0, 0.0, 0.0), -0.098 cm^-1, (0.0, 0.0, 0.0), (1.0, 1.0, 1.0))
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((-2.8160000000000003 cm, 0.0 cm, -3.0 cm), (5.248 cm, 1.9840000000000002 cm, 2.52 cm), (-1.2566370614359172, 0.0, 0.0), -0.002 cm^-1, (-0.9510565162951535, 0.0, 0.0), (0.30901699437494745, 1.0, 1.0))
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((2.8160000000000003 cm, 0.0 cm, -3.0 cm), (3.968 cm, 1.364 cm, 2.64 cm), (1.2566370614359172, 0.0, 0.0), -0.002 cm^-1, (0.9510565162951535, 0.0, 0.0), (0.30901699437494745, 1.0, 1.0))
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((0.0 cm, 4.34 cm, -3.0 cm), (2.688 cm, 3.1 cm, 4.199999999999999 cm), (0.0, 0.0, 0.0), 0.001 cm^-1, (0.0, 0.0, 0.0), (1.0, 1.0, 1.0))
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((0.0 cm, 1.2400000000000002 cm, -3.0 cm), (0.5888 cm, 0.5704 cm, 0.552 cm), (0.0, 0.0, 0.0), 0.001 cm^-1, (0.0, 0.0, 0.0), (1.0, 1.0, 1.0))
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((-1.024 cm, -7.502 cm, -3.0 cm), (0.5888 cm, 0.2852 cm, 0.24 cm), (0.0, 0.0, 0.0), 0.001 cm^-1, (0.0, 0.0, 0.0), (1.0, 1.0, 1.0))
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((0.0 cm, -1.2400000000000002 cm, -3.0 cm), (0.5888 cm, 0.5704 cm, 0.552 cm), (0.0, 0.0, 0.0), 0.001 cm^-1, (0.0, 0.0, 0.0), (1.0, 1.0, 1.0))
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((0.0 cm, -7.502 cm, -3.0 cm), (0.2944 cm, 0.2852 cm, 0.276 cm), (0.0, 0.0, 0.0), 0.001 cm^-1, (0.0, 0.0, 0.0), (1.0, 1.0, 1.0))
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((0.768 cm, -7.502 cm, -3.0 cm), (0.5888 cm, 0.2852 cm, 0.24 cm), (-1.5707963267948966, 0.0, 0.0), 0.001 cm^-1, (-1.0, 0.0, 0.0), (6.123233995736766e-17, 1.0, 1.0))
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((0.768 cm, -1.302 cm, 0.75 cm), (0.7168000000000001 cm, 0.49600000000000005 cm, 1.2000000000000002 cm), (-1.5707963267948966, 0.0, 0.0), 0.002 cm^-1, (-1.0, 0.0, 0.0), (6.123233995736766e-17, 1.0, 1.0))
 ImagePhantoms.Object3d{ImagePhantoms.Ellipsoid, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 3, Float64}((0.0 cm, 1.2400000000000002 cm, 7.5 cm), (0.7168000000000001 cm, 0.6944 cm, 1.2000000000000002 cm), (0.0, 0.0, 0.0), -0.002 cm^-1, (0.0, 0.0, 0.0), (1.0, 1.0, 1.0))

A narrow grayscale window is needed to see the soft tissue structures

clim = (0.95, 1.05) .* μ

(0.095 cm^-1, 0.10500000000000001 cm^-1)

Here is the ideal phantom image:

oversample = 3
true_image = phantom(axes(ig)..., ob, oversample)
pt = jim(axes(ig), true_image, "True 3D Shepp-Logan phantom image"; clim)

Define the system geometry (for some explanation use ?CtGeom):

p = (ns = 130, ds = 0.3cm, nt = 80, dt = 0.4cm, na = 50, dsd = 200cm, dod = 40cm)
rg = CtFanArc( ; p...)

CtFanArc{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float32, CtSourceCircle} :
 ns::Int64 130
 nt::Int64 80
 ds::Unit{Float64} 0.3 cm
 dt::Unit{Float64} 0.4 cm
 offset_s::Float32 0.0
 offset_t::Float32 0.0
 na::Int64 50
 orbit::Float32 360.0
 orbit_start::Float32 0.0
 source_offset::Unit{Float64} 0.0 cm
 dsd::Unit{Float64} 200.0 cm
 dod::Unit{Float64} 40.0 cm
 src::CtSourceCircle CtSourceCircle()

Examine the geometry to verify the FOV (this is more interesting when interacting via other Plot backends):

ct_geom_plot3(rg, ig)

prompt()

CBCT projections using Sinogram.rays and ImagePhantoms.radon:

proj_arc = radon(rays(rg), ob)
pa = jim(axes(rg)[1:2], proj_arc ;
    title="Shepp-Logan projections (arc)", xlabel="s", ylabel="t")

There is no "inverse crime" here because we compute the projection views using the analytical phantom geometry, but then reconstruct on a discrete grid.

Image reconstruction via FBP / FDK

We start with a "plan", which would save work if we were reconstructing many images. For illustration we include Hamming window.

plan = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))
fdk_arc = fdk(plan, proj_arc)
par = jim(axes(ig), fdk_arc, "FDK image (arc)"; clim)

err_arc = fdk_arc - true_image
elim = (-1,1) .* (0.02μ)
pae = jim(axes(ig), err_arc, "Error image (arc)"; clim = elim)

Repeat with flat detector geometry

rg = CtFanFlat( ; p...)
proj_flat = radon(rays(rg), ob)
pfp = jim(axes(rg)[1:2], proj_flat ;
    title="Shepp-Logan projections (flat)", xlabel="s", ylabel="t")

plan = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))
fdk_flat = fdk(plan, proj_flat)
pfr = jim(axes(ig), fdk_flat, "FDK image (flat)"; clim)

err_flat = fdk_flat - true_image
pfe = jim(axes(ig), err_flat, "Error image (flat)"; clim = elim)

As expected for CBCT, the largest errors are in the end slices.

Short scan

rg = CtFanFlat(:short ; p..., na = 200)
proj_short = radon(rays(rg), ob)
psp = jim(axes(rg)[1:2], proj_short ;
    title="Shepp-Logan projections (flat,short)", xlabel="s", ylabel="t")

plan_short = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))
psw = jim(axes(rg)[[1,3]], plan_short.view_weight[:,1,:];
    title = "View weights (including Parker)",
    xlabel="s", ylabel="angle")

fdk_short = fdk(plan_short, proj_short)
psr = jim(axes(ig), fdk_short, "FDK image (flat,short)"; clim)

err_short = fdk_short - true_image
pse = jim(axes(ig), err_flat, "Error image (flat,short)"; clim = elim)

This page was generated using Literate.jl.