2D sinogram geometry¶

This page describes the 2D sinogram geometries available in the Julia package Sinograms.jl.

Setup¶

Packages needed here.

In [ ]:

using Unitful: mm, °
using Plots # these 2 must precede 'using Sinograms' for sino_plot_rays to work
using Sinograms: SinoPar, SinoMoj, SinoFanArc, SinoFanFlat, SinoFan
using Sinograms: sino_plot_rays, rays, angles
using MIRTjim: jim, prompt
using InteractiveUtils: versioninfo

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.

In [ ]:

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

2D Sinogram geometries¶

To perform 2D image reconstruction from a 2D sinogram, one must first describe how the rays in the sinogram are sampled.

Mathematically, the Radon transform $p(r,\phi)$ of $f(x,y)$ is defined for all $r$ and $ϕ$ but practical systems record $p(r,ϕ)$ only for certain finite sets of samples of those values.

Parallel-beam geometry¶

The simplest form of sinogram sampling is equally spaced samples in both $r$ and $ϕ$. This sampling is called a parallel-beam geometry. (Very few practical systems have this sampling pattern, but it is an easy place to start.) Both FBP and iterative reconstruction methods need to know which values of $r$ and $ϕ$ are sampled.

In this package, SinoPar is the type that describes a parallel-beam sinogram geometry.

The built-in defaults provide helpful reminders about the usage.

In [ ]:

SinoPar()

This package supports units via the Unitful.jl and using units is recommended (but not required).

Here is an example of how to specify a parallel-beam geometry. Everything is a named keyword argument with sensible default values.

orbit and orbit_start must both be unitless (degrees) or have the same units (e.g., degrees or radians).
detector spacing d (can be unitless or have units e.g., mm).
The projection angles $ϕ$ are equally space and given by orbit_start + (0:(nb-1))/nb * orbit.

In [ ]:

rg = SinoPar( ;
    nb = 64, # number of radial samples ("bins")
    na = 30, # number of angular samples
    d = 2mm, # detector spacing
    offset = 0.25, # quarter detector offset (unitless)
    orbit = 180, # angular range (in degrees)
    orbit_start = 0, # starting angle (in degrees)
)

The struct rg has numerous useful properties; type ?SinoGeom to see the full list.

For example, to access the angular samples in degrees type angles(rg)

In [ ]:

angles(rg)

The following function visualizes the sampling pattern.

In [ ]:

sino_plot_rays(rg; ylims=(0,180), yticks=(0:90:180), widen=true, title="Parallel")

In [ ]:

prompt()

Fan-beam CT with an arc detector (3rd generation CT)¶

For a fan-beam geometry, the arguments are the same as for SinoPar with the addition of specifying:

dsd distance from source to detector
dod distance from origin (isocenter) to detector
dfs distance from focal point of detector to source (0 for a 3rd gen arc detector, and Inf for a flat detector)
source_offset for misaligned systems where the ray from the source to the detector center does not intersect the isocenter. Not fully supported; submit an issue if you need this feature.

Here is an example that corresponds to a GE Lightspeed CT system. These numbers are published in IEEE T-MI Oct. 2006, p.1272-1283.

In [ ]:

rg = SinoFanArc( ; nb=888, na=984,
    d=1.0239mm, offset=1.25, dsd=949.075mm, dod=408.075mm)

Here is a smaller example for plotting the rays.

In [ ]:

rg = SinoFanArc( ; nb=64, na=30,
    d=20mm, offset=0.25, dsd=900mm, dod=400mm)
sino_plot_rays(rg; ylims=(-50,400), yticks=(0:180:360), widen=true,
    title="Fan-beam for arc detector")

In [ ]:

prompt()

Fan-beam CT with a flat detector¶

This geometry is the same as the arc detector except that dfs=Inf.

In [ ]:

rg = SinoFanFlat( ; nb=64, na=30,
    d=20mm, offset=0.25, dsd=900mm, dod=400mm)

Here is its sampling plot

In [ ]:

sino_plot_rays(rg; ylims=(-50,400), yticks=(0:180:360), widen=true,
    title="Fan-beam for flat detector")

In [ ]:

prompt()

Mojette sampling¶

This is a specialized sampling geometry that is currently incompletely supported.

In [ ]:

rg = SinoMoj( ; nb=60, na=30)

Here is its sampling plot

In [ ]:

sino_plot_rays(rg; ylims=(0,180), yticks=(0:90:180), widen=true,
    title="Mojette sampling")

Here is a diagram that illustrates how the radial spacing is a function of the projection view angle for the Mojette geometry. The key aspect here is that in each image row the line intersection lengths are identical.

In [ ]:

plot(xlims=(-1,1) .* 3.5, ylims = (-1,1) .* 2.5, xlabel="x", ylabel="x")
default(label="")
for y in -2:2
    plot!([-3, 3], [y, y], color=:black)
end
for x in -3:3
    plot!([x, x], [-2, 2], color=:black)
end
plot_ray!(r, ϕ) = plot!(
    r*cos(ϕ) .+ [-1, 1] * 4 * sin(ϕ),
    r*sin(ϕ) .+ [1, -1] * 4 * cos(ϕ),
    color = :blue,
)
ia = 4 # pick an angle
i = rays(rg) # iterator
rs = [i[1] for i in i] # radial samples
ϕs = [i[2] for i in i] # projection angle
r = rs[:,ia]
r = r[abs.(r) .< 3]
ϕ = ϕs[1,ia]
plot_ray!.(r, ϕ)
plot!(aspect_ratio=1, title = "Mojette line integrals")

This notebook was generated using Literate.jl.