GitHub

SPECTrecon overview

This page gives an overview of the Julia package SPECTrecon.

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

Setup

Packages needed here.

using SPECTrecon: plan_psf, psf_gauss, SPECTplan
using SPECTrecon: project, project!, backproject, backproject!
using MIRTjim: jim, prompt
using LinearAlgebra: mul!
using LinearMapsAA: LinearMapAA
using Plots: scatter, plot!, default; default(markerstrokecolor=:auto)
using Plots # @animate, gif
using InteractiveUtils: versioninfo

The following line is helpful when running this example.jl 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);

Overview

To perform SPECT image reconstruction, one must have a model for the imaging system encapsulated in a forward projector and back projector.

Mathematically, we write the forward projection process in SPECT as "y = A * x" where A is a "system matrix" that models the physics of the imaging system (including depth-dependent collimator/detector response and attenuation) and "x" is the current guess of the emission image.

However, in code we usually cannot literally store "A" as dense matrix because it is too large. A typical size in SPECT is that the image x is nx × ny × nz = 128 × 128 × 100 and the array of projection views y is nx × nz × nview = 128 × 100 × 120. So the system matrix A has 1536000 × 1638400 elements which is far to many to store, even accounting for some sparsity.

Instead, we write functions called forward projectors that calculate A * x "on the fly".

Similarly, the operation A' * y is called "back projection", where A' denotes the transpose or "adjoint" of A.

Example

To illustrate forward and back projection, it is easiest to start with a simulation example using a digital phantom. The fancy way would be to use a 3D phantom from ImagePhantoms, but instead we just use two simple cubes.

nx,ny,nz = 128,128,80
T = Float32
xtrue = zeros(T, nx,ny,nz)
xtrue[(1nx÷4):(2nx÷3), 1ny÷5:(3ny÷5), 2nz÷6:(3nz÷6)] .= 1
xtrue[(2nx÷5):(3nx÷5), 1ny÷5:(2ny÷5), 4nz÷6:(5nz÷6)] .= 2

average(x) = sum(x) / length(x)
function mid3(x::AbstractArray{T,3}) where {T}
    (nx,ny,nz) = size(x)
    xy = x[:,:,ceil(Int, nz÷2)]
    xz = x[:,ceil(Int,end/2),:]
    zy = x[ceil(Int, nx÷2),:,:]'
    return [xy xz; zy fill(average(x), nz, nz)]
end
jim(mid3(xtrue), "Middle slices of x")
Example block output

PSF

Create a synthetic depth-dependent PSF for a single view

px = 11
psf1 = psf_gauss( ; ny, px, fwhm_end = 6)
jim(psf1, "PSF for each of $ny planes")
Example block output

In general the PSF can vary from view to view due to non-circular detector orbits. For simplicity, here we illustrate the case where the PSF is the same for every view.

nview = 60
psfs = repeat(psf1, 1, 1, 1, nview)
size(psfs)
(11, 11, 128, 60)

Plan the PSF modeling (see 3-psf.jl)

plan = plan_psf( ; nx, nz, px)
1-element Vector{PlanPSF{Float32, FFTW.cFFTWPlan{ComplexF32, -1, true, 2, Tuple{Int64, Int64}}, AbstractFFTs.ScaledPlan{ComplexF32, FFTW.cFFTWPlan{ComplexF32, 1, true, 2, UnitRange{Int64}}, Float32}}} with N=128

Basic SPECT forward projection

Here is a simple illustration of a SPECT forward projection operation. (This is a memory inefficient way to do it!)

dy = 4 # transaxial pixel size in mm
mumap = zeros(T, size(xtrue)) # μ-map just zero for illustration here
views = project(xtrue, mumap, psfs, dy)
size(views)
(128, 80, 60)

Display the calculated (i.e., simulated) projection views

jim(views[:,:,1:4:end], "Every 4th of $nview projection views")
Example block output

Basic SPECT back projection

This illustrates an "unfiltered backprojection" that leads to a very blurry image (again, with a simple memory inefficient usage).

First, back-project two "rays" to illustrate the depth-dependent PSF.

tmp = zeros(T, size(views))
tmp[nx÷2, nz÷2, nview÷5] = 1
tmp[nx÷2, nz÷2, 1] = 1
tmp = backproject(tmp, mumap, psfs, dy)
jim(mid3(tmp), "Back-projection of two rays")
Example block output

Now back-project all the views of the phantom.

back = backproject(views, mumap, psfs, dy)
jim(mid3(back), "Back-projection of ytrue")
Example block output

Memory efficiency

For iterative reconstruction, one must do forward and back-projection repeatedly. It is more efficient to pre-allocate work arrays for those operations, instead of repeatedly making system calls.

Here we illustrate the memory efficient versions that are recommended for iterative SPECT reconstruction.

First construction the SPECT plan.

#viewangle = (0:(nview-1)) * 2π # default
plan = SPECTplan(mumap, psfs, dy; T)
SPECTplan{Float32}
 imgsize::Tuple{Int64, Int64, Int64} (128, 128, 80)
 px::Int64 11
 pz::Int64 11
 nview::Int64 60
 viewangle::StepRangeLen{Float32, Float64, Float64, Int64} 0.0f0:0.10471976f0:6.178466f0
 interpmeth::Symbol two
 mode::Symbol fast
 dy::Float32 4.0
 nthread::Int64 1
 mumap: 128×128×80 Array{Float32, 3}
 (8960151 bytes)

Mutating version of forward projection:

tmp = Array{T}(undef, nx, nz, nview)
project!(tmp, xtrue, plan)
@assert tmp == views

Mutating version of back-projection:

tmp = Array{T}(undef, nx, ny, nz)
backproject!(tmp, views, plan)
@assert tmp ≈ back

Using LinearMapAA

Calling project! and backproject! repeatedly leads to application-specific code. More general code uses the fact that SPECT projection and back-projection are linear operations, so we use LinearMapAA to define a "system matrix" for these operations.

forw! = (y,x) -> project!(y, x, plan)
back! = (x,y) -> backproject!(x, y, plan)
idim = (nx,ny,nz)
odim = (nx,nz,nview)
A = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim)
LinearMapAO: 614400 × 1310720 odim=(128, 80, 60) idim=(128, 128, 80) T=Float32 Do=3 Di=3
map = 614400×1310720 LinearMaps.FunctionMap{Float32,true}(#5, #7; issymmetric=false, ishermitian=false, isposdef=false)

Simple forward and back-projection:

@assert A * xtrue == views
@assert A' * views ≈ back

Mutating version:

tmp = Array{T}(undef, nx, nz, nview)
mul!(tmp, A, xtrue)
@assert tmp == views
tmp = Array{T}(undef, nx, ny, nz)
mul!(tmp, A', views)
@assert tmp ≈ back

Units

The pixel dimensions deltas can (and should!) be values with units.

Here is an example with units ... (todo)

using UnitfulRecipes using Unitful: mm

Projection view animation

anim = @animate for i in 1:nview
    ymax = maximum(views)
    jim(views[:,:,i],
        "SPECT projection view $i of $nview",
        clim = (0, ymax),
    )
end
gif(anim, "views.gif", fps = 8)
Example block output

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')
12-element Vector{SubString{String}}:
 "Julia Version 1.10.3"
 "Commit 0b4590a5507 (2024-04-30 10:59 UTC)"
 "Build Info:"
 "  Official https://julialang.org/ release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LIBM: libopenlibm"
 "  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)"
 "Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()
Status `~/work/SPECTrecon.jl/SPECTrecon.jl/docs/Project.toml`
  [fbb218c0] BSON v0.3.9
  [31c24e10] Distributions v0.25.108
  [e30172f5] Documenter v1.4.1
  [587475ba] Flux v0.14.15
  [71a99df6] ImagePhantoms v0.8.1
  [599c1a8e] LinearMapsAA v0.12.0
  [98b081ad] Literate v2.18.0
  [170b2178] MIRTjim v0.24.0
  [872c559c] NNlib v0.9.17
  [91a5bcdd] Plots v1.40.4
  [ab1be465] SPECTrecon v0.2.0 `~/work/SPECTrecon.jl/SPECTrecon.jl`
  [1986cc42] Unitful v1.20.0
  [42071c24] UnitfulRecipes v1.6.1
  [e88e6eb3] Zygote v0.6.70
  [700de1a5] ZygoteRules v0.2.5
  [b77e0a4c] InteractiveUtils
  [37e2e46d] LinearAlgebra

This page was generated using Literate.jl.