Fan-beam tomography: arc detector
This page describes fan-beam tomographic image reconstruction using the Julia package Sinograms.jl.
This page focuses on fan-beam with an "arc" detector, i.e., 3rd-generation CT systems where the detector is an arc having a focal point at the X-ray source. This geometry is popular in part because it facilitates anti-scatter grids.
This page comes from a single Julia file: 04-fan-arc.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: 04-fan-arc.ipynb, or open it in binder here: 04-fan-arc.ipynb.
Setup
Packages needed here.
using Plots: plot, gui
using Unitful: mm
using Sinograms: SinoFanArc, rays, plan_fbp, fbp, sino_geom_plot!
using ImageGeoms: ImageGeom, fovs, MaskCircle
using ImagePhantoms: SheppLogan, shepp_logan, radon, phantom
using MIRTjim: jim, promptThe 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);Fan-beam sinogram of Shepp-Logan phantom
For illustration, we start by synthesizing a fan-beam sinogram of the 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=(128,126), deltas = (2mm,2mm) )ImageGeoms.ImageGeom{2, NTuple{2,Unitful.Quantity{Int64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, BitMatrix}
dims::NTuple{2,Int64} (128, 126)
deltas::NTuple{2,Unitful.Quantity{Int64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}} (2 mm, 2 mm)
offsets::NTuple{2,Float32} (0.0f0, 0.0f0)
mask: 128×126 BitMatrix {12096 of 16128}
Use SinoFanArc to define the sinogram geometry.
rg = SinoFanArc( ; nb = 130, d = 3.2mm, na = 100, dsd = 400mm, dod = 140mm)SinoFanArc{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float32} :
nb::Int64 130
d::Unit{Float64} 3.2 mm
offset::Float32 0.0
na::Int64 100
orbit::Float32 360.0
orbit_start::Float32 0.0
source_offset::Unit{Float64} 0.0 mm
dsd::Unit{Float64} 400.0 mm
dod::Unit{Float64} 140.0 mm
Examine the geometry to verify the FOV:
jim(axes(ig), ig.mask; prompt=false)
sino_geom_plot!(rg, ig)
prompt()Ellipse parameters for Shepp-Logan phantom:
μ = 0.01 / mm # typical linear attenuation coefficient
ob = shepp_logan(SheppLogan(); fovs = fovs(ig), u = (1, 1, μ))10-element Vector{ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 1, Float64}}:
ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 mm, 0.0 mm), (117.76 mm, 86.94 mm), (1.5707963267948966,), 0.02 mm^-1, (1.0,), (6.123233995736766e-17,))
ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 mm, -2.3184 mm), (111.872 mm, 83.4624 mm), (1.5707963267948966,), -0.0098 mm^-1, (1.0,), (6.123233995736766e-17,))
ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 1, Float64}((28.16 mm, 0.0 mm), (39.68 mm, 13.86 mm), (1.2566370614359172,), -0.0002 mm^-1, (0.9510565162951535,), (0.30901699437494745,))
ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 1, Float64}((-28.16 mm, 0.0 mm), (52.48 mm, 20.16 mm), (1.8849555921538759,), -0.0002 mm^-1, (0.9510565162951536,), (-0.30901699437494734,))
ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 mm, 44.099999999999994 mm), (32.0 mm, 26.459999999999997 mm), (1.5707963267948966,), 0.0001 mm^-1, (1.0,), (6.123233995736766e-17,))
ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 mm, 12.600000000000001 mm), (5.888 mm, 5.796 mm), (0.0,), 0.0001 mm^-1, (0.0,), (1.0,))
ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 mm, -12.600000000000001 mm), (5.888 mm, 5.796 mm), (0.0,), 0.0001 mm^-1, (0.0,), (1.0,))
ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 1, Float64}((-10.24 mm, -76.23 mm), (5.888 mm, 2.898 mm), (0.0,), 0.0001 mm^-1, (0.0,), (1.0,))
ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 mm, -76.23 mm), (2.944 mm, 2.898 mm), (0.0,), 0.0001 mm^-1, (0.0,), (1.0,))
ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Float64, 1, Float64}((7.68 mm, -76.23 mm), (5.888 mm, 2.898 mm), (1.5707963267948966,), 0.0001 mm^-1, (1.0,), (6.123233995736766e-17,))Arc fan-beam sinogram for Shepp-Logan phantom:
sino = radon(rays(rg), ob)
jim(axes(rg), sino; title="Shepp-Logan sinogram", xlabel="r", ylabel="ϕ")
Image reconstruction via FBP
Here we start with a "plan", which would save work if we were reconstructing many images.
plan = plan_fbp(rg, ig)
fbp_image = fbp(plan, sino)128×126 Matrix{Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(mm^-1,), 𝐋^-1, nothing}}}:
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 … 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 … 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
⋮ ⋱ ⋮
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 … 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 … 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1
0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1 0.0 mm^-1A narrow color window is needed to see the soft tissue structures:
clim = (0.9, 1.1) .* μ
jim(axes(ig), fbp_image, "FBP image for arc case"; clim)
For comparison, here is the ideal phantom image
true_image = phantom(axes(ig)..., ob, 2)
jim(axes(ig)..., true_image, "True phantom image"; clim)
This page was generated using Literate.jl.