Fan-beam tomography: flat detector

This page describes fan-beam tomographic image reconstruction using the Julia package Sinograms.jl.

This page focuses on fan-beam with a "flat" detector, i.e., one row of a flat-panel detector as used in many cone-beam CT (CBCT) systems.

This page comes from a single Julia file: 05-fan-flat.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: 05-fan-flat.ipynb, or open it in binder here: 05-fan-flat.ipynb.

Setup

Packages needed here.

using Plots: plot, gui
using Unitful: cm
using Sinograms: SinoFanFlat, rays, plan_fbp, Window, Hamming, fbp, sino_geom_plot!
using ImageGeoms: ImageGeom, fovs, MaskCircle
using ImagePhantoms: SheppLogan, shepp_logan, 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);

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 = (0.2cm,0.2cm) )

ImageGeoms.ImageGeom{2, NTuple{2,Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}}, BitMatrix}
 dims::NTuple{2,Int64} (128, 126)
 deltas::NTuple{2,Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}} (0.2 cm, 0.2 cm)
 offsets::NTuple{2,Float32} (0.0f0, 0.0f0)
 mask: 128×126 BitMatrix {12096 of 16128}

Use SinoFanFlat to define the sinogram geometry.

rg = SinoFanFlat( ; nb = 130, d = 0.3cm, na = 100, dsd = 50cm, dod = 14cm)

SinoFanFlat{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float32} :
 nb::Int64 130
 d::Unit{Float64} 0.3 cm
 offset::Float32 0.0
 na::Int64 100
 orbit::Float32 360.0
 orbit_start::Float32 0.0
 source_offset::Unit{Float64} 0.0 cm
 dsd::Unit{Float64} 50.0 cm
 dod::Unit{Float64} 14.0 cm

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.1 / cm # 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{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 1, Float64}}:
 ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 cm, 0.0 cm), (11.776000000000002 cm, 8.694 cm), (1.5707963267948966,), 0.2 cm^-1, (1.0,), (6.123233995736766e-17,))
 ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 cm, -0.23184000000000002 cm), (11.1872 cm, 8.346240000000002 cm), (1.5707963267948966,), -0.098 cm^-1, (1.0,), (6.123233995736766e-17,))
 ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 1, Float64}((2.8160000000000003 cm, 0.0 cm), (3.968 cm, 1.3860000000000001 cm), (1.2566370614359172,), -0.002 cm^-1, (0.9510565162951535,), (0.30901699437494745,))
 ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 1, Float64}((-2.8160000000000003 cm, 0.0 cm), (5.248 cm, 2.0160000000000005 cm), (1.8849555921538759,), -0.002 cm^-1, (0.9510565162951536,), (-0.30901699437494734,))
 ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 cm, 4.41 cm), (3.2 cm, 2.6460000000000004 cm), (1.5707963267948966,), 0.001 cm^-1, (1.0,), (6.123233995736766e-17,))
 ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 cm, 1.2600000000000002 cm), (0.5888 cm, 0.5796 cm), (0.0,), 0.001 cm^-1, (0.0,), (1.0,))
 ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 cm, -1.2600000000000002 cm), (0.5888 cm, 0.5796 cm), (0.0,), 0.001 cm^-1, (0.0,), (1.0,))
 ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 1, Float64}((-1.024 cm, -7.623 cm), (0.5888 cm, 0.2898 cm), (0.0,), 0.001 cm^-1, (0.0,), (1.0,))
 ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 1, Float64}((0.0 cm, -7.623 cm), (0.2944 cm, 0.2898 cm), (0.0,), 0.001 cm^-1, (0.0,), (1.0,))
 ImagePhantoms.Object2d{ImagePhantoms.Ellipse, Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64, 1, Float64}((0.768 cm, -7.623 cm), (0.5888 cm, 0.2898 cm), (1.5707963267948966,), 0.001 cm^-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. For illustration we include Hamming window.

plan = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))
fbp_image = fbp(plan, sino)

128×126 Matrix{Unitful.Quantity{Float64, 𝐋^-1, Unitful.FreeUnits{(cm^-1,), 𝐋^-1, nothing}}}:
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1  …  0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1  …  0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
         ⋮                        ⋱                                ⋮
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1  …  0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1  …  0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1
 0.0 cm^-1  0.0 cm^-1  0.0 cm^-1     0.0 cm^-1  0.0 cm^-1  0.0 cm^-1

A narrow color window is needed to see the soft tissue structures:

clim = (0.9, 1.1) .* μ
jim(axes(ig), fbp_image, "FBP image for flat 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.