SPECTrecon

System matrix (forward and back-projector) for SPECT image reconstruction.

Rotate1D = RotateMode{:one}()

Rotate using 3-pass 1D interpolation

Rotate2D = RotateMode{:two}()

Rotate using 1-pass 2D bilinear interpolation

PlanPSF{T,Tf,Ti}( ; nx::Int, nz::Int, px::Int, pz::Int, T::Type)

Struct for storing work arrays and factors for 2D convolution for one thread. Each PSF is px × pz

• T datatype of work arrays (subtype of AbstractFloat)
• nx::Int = 128 (ny implicitly the same)
• nz::Int = nx image size is [nx,nx,nz]
• px::Int = 1
• pz::Int = px (PSF size)
• padsize::Tuple : (padup, paddown, padleft, padright)
• workmat [nx+padup+paddown, nz+padleft+padright] 2D padded image for FFT convolution
• workvecx [nx+padup+paddown,]: 1D work vector
• workvecz [nz+padleft+padright,]: 1D work vector
• img_compl [nx+padup+paddown, nz+padleft+padright]: 2D [complex] padded image for FFT
• ker_compl [nx+padup+paddown, nz+padleft+padright]: 2D [complex] padded image for FFT
• fft_plan::Tf plan for doing FFT; see plan_fft!
• ifft_plan::Ti plan for doing IFFT; see plan_ifft!

PlanRotate{T, R}

Struct for storing work arrays and factors for 2D square image rotation for one thread

• T datatype of work arrays (default Float32)
• R::RotateMode
• nx::Int image size
• padsize::Int : pad each side of image by this much
• workmat1 [nx + 2 * padsize, nx + 2 * padsize] padded work matrix
• workmat2 [nx + 2 * padsize, nx + 2 * padsize] padded work matrix
• interp::SparseInterpolator{T, 2, 1}

RotateMode

Type to control image rotation method.

• Rotate1D : 3-pass 1D interpolation
• Rotate2D : 1-pass 2D bilinear interpolation

SPECTplan

Struct for storing key factors for a SPECT system model

• T datatype of work arrays
• imgsize size of image
• px,pz psf dimension
• imgr [nx, ny, nz] 3D rotated version of image
• add_img [nx, ny, nz] 3D image for adding views and backprojection
• mumap [nx,ny,nz] attenuation map, must be 3D, possibly zeros()
• mumapr [nx, ny, nz] 3D rotated mumap
• exp_mumapr [nx, nz] 2D exponential rotated mumap
• psfs [px,pz,ny,nview] point spread function, must be 4D, with px andpz` odd, and symmetric for each slice
• nview number of views, must be integer
• viewangle set of view angles, must be from 0 to 2π
• interpmeth interpolation method: :one means 1d; :two means 2d
• mode pre-allocation method: :fast means faster; :mem means use less memory
• dy voxel size in y direction (dx is the same value)
• nthread number of CPU threads used to process data; must be integer
• planrot Vector of struct PlanRotate
• planpsf Vector of struct PlanPSF

Currently code assumes the following:

• each of the nview projection views is [nx,nz]
• nx = ny
• uniform angular sampling
• psf is symmetric
• multiprocessing using # of threads specified by Threads.nthreads()

Constructor is not type stable due to the Union, but this is OK because construction is done just once before iterating.

show(io::IO, mime::MIME"text/plain", vp::Vector{<:PlanPSF})

show(io::IO, mime::MIME"text/plain", vp::Vector{<:PlanRotate})

show(io::IO, ::MIME"text/plain", plan::PlanRotate)

show(io::IO, ::MIME"text/plain", plan::PlanPSF)

show(io::IO, ::MIME"text/plain", plan::SPECTplan)

sizeof(::PlanPSF)

Show size in bytes of PlanPSF object.

sizeof(::SPECTplan)

Show size in bytes of SPECTplan object.

sizeof(::PlanRotate)

Show size in bytes of PlanRotate object.

Ablock(plan, nblocks)

Generate a vector of linear maps for OSEM. -plan: SPECTrecon plan -nblocks: Number of blocks in OSEM

backproject!(image, views, plan ; index)

Backproject multiple views into image. Array image is not initialized to zero; caller must do that.

backproject!(image, view, plan, buffer_id, viewidx)

Backproject a single view.

backproject!(image, view, plan, viewidx)

Backproject a single view.

image = backproject(views, mumap, psfs, dy; interpmeth, kwargs...)

SPECT backproject views using attenuation map mumap and PSF array psfs for pixel size dy. This method initializes the plan as a convenience. Most users should use backproject! instead after initializing those, for better efficiency.

image = backproject(views, plan ; kwargs...)

SPECT backproject views; this allocates the returned 3D array.

copy3dj!(mat2d, mat3d, j)

Non-allocating mat2d .= mat3d[:,j,:]

fft_conv!(output, image3, ker3, plans)

Mutating version of convolving a 3D image3 with a set of 2D symmetric kernels stored in 3D array ker3 using foreach.

fft_conv!(output, img, ker, plan)

Convolve 2D image img with 2D (symmetric!) kernel ker using FFT, storing the result in output.

fft_conv(img, ker; T)

Convolve 2D image img with 2D (symmetric!) kernel ker using FFT.

fft_conv_adj!(output, image3, ker3, plans)

Adjoint of fft_conv.

fft_conv_adj!(output, img, ker, plan)

Adjoint of fft_conv!.

fft_conv_adj(img, ker; T)

Adjoint of fft_conv.

fft_conv_adj2!(output, image2, ker3, plans)

Mutating version of adjoint of convolving a 2D image2 with each 2D kernel in the 3D array ker3.

fftshift2!(dst, src)

Same as fftshift in 2d, but non-allocating

imfilterz!(plan)

FFT-based convolution of plan.img_compl and kernel plan.ker_compl (not centered), storing result in plan.workmat. Over-writes plan.ker_compl.

imrotate!(output, image3, θ, plans, ntasks=nthreads)

In-place version of rotating a 3D image3 by θ in counter-clockwise direction (opposite to imrotate in Julia) using foreach with ntasks.

imrotate!(output, image3, θ, plans, :thread)

Mutating version of rotating a 3D image3 by θ in counter-clockwise direction (opposite of imrotate in Julia) writing result into output, using Threads.@threads.

imrotate!(output, img, θ, plan)

Rotate a 2D image img by angle θ ∈ [0,2π] in counter-clockwise direction (opposite to imrotate in Julia) using 3-pass 1d linear interpolations.

imrotate!(output, img, θ, plan)

Rotate a 2D image img by angle θ ∈ [0,2π] in counter-clockwise direction (opposite to imrotate in Julia) using 2d bilinear interpolation.

imrotate(img, θ; method::Symbol=:two)

Rotate a 2D image img by angle θ ∈ [0,2π] in counter-clockwise direction (opposite to imrotate in Julia) using either 2d linear interpolation (for :two) or 3-pass 1D interpolation (for :one)

imrotate_adj!(output, image3, θ, plans, ntasks=nthreads)

Adjoint of imrotate! using foreach.

imrotate_adj!(output, image3, θ, plans, :thread)

Adjoint of imrotate! using @threads.

imrotate_adj!(output, img, θ, plan)

Adjoint of imrotate! for a 2D image using 3-pass 1d linear interpolations.

imrotate_adj!(output, img, θ, plan)

Adjoint of imrotate! for a 2D image using 2D bilinear interpolation.

imrotate_adj(img, θ; method::Symbol=:two)

Adjoint of rotating a 2D image img by angle θ ∈ [0,2π] in counter-clockwise direction (opposite to imrotate in Julia) using either 2d linear interpolations or 3-pass 1D interpolation.

linearinterp!(A, s, e, len)

Assign key values in SparseInterpolator (linear) A that are calculated from s, e and len. s means start (x[1]) e means end (x[end]) len means length (length(x))

mlem!(out, x0, ynoisy, background, A; niter = 20)

Inplace version of ML-EM algorithm for emission tomography image reconstruction.

• out: Output
• x0: Initial guess
• ynoisy: (Noisy) measurements
• background: Background effects, e.g., scatters
• A: System matrix
• niter: Number of iterations

mlem(x0, ynoisy, background, A; niter = 20)

ML-EM algorithm for emission tomography image reconstruction.

• x0: Initial guess
• ynoisy: (Noisy) measurements
• background: Background effects, e.g., scatters
• A: System matrix
• niter: Number of iterations

mul3dj!(mat3d, mat2d, j)

Non-allocating mat3d[:,j,:] *= mat2d

osem!(out, x0, ynoisy, background, Ab; niter = 20)

OS-EM algorithm for SPECT reconstruction.

• out: Output
• x0: Initial estimate
• ynoisy: (Noisy) measurements
• background: Background effects, e.g., scatters
• Ab: Vector of system matrix
• niter: Number of iterations

osem(x0, ynoisy, background, Ab; niter = 20)

OS-EM algorithm for SPECT reconstruction.

• x0: Initial guess
• ynoisy: (Noisy) measurements
• background: Background effects, e.g., scatters
• Ab: Vector of system matrix
• niter: Number of iterations

pad2sizezero!(output, img, padsize)

Non-allocating version of padding: `output[paddims[1]+1 : paddims[1]+dims[1], paddims[2]+1 : paddims[2]+dims[2]]) .= img

pad_it!(X, padsize)

Zero-pad X to padsize

padrepl!(output, img, padsize)

Pad with replication from img into output

padzero!(output, img, pad_x, pad_y)

Mutating version of padding a 2D image by filling zeros. Output has size (size(img, 1) + padsize[1] + padsize[2], size(img, 2) + padsize[3] + padsize[4]).

plan_psf( ; nx::Int, nz::Int, px::Int, pz::Int, nthread::Int, T::Type)

Make Vector of structs for storing work arrays and factors for 2D convolution with SPECT depth-dependent PSF model, threaded across planes parallel to detector. Option

• nx::Int = 128
• nz::Int = nx
• px::Int = 1
• pz::Int = px PSF size is px × pz
• T : datatype of work arrays, defaults to Float32
• nthread::Int # of threads, defaults to Threads.nthreads()

plan_rotate(nx::Int; T::Type{<:AbstractFloat}, method::Symbol)

Make Vector of PlanRotate structs for storing work arrays and factors for threaded rotation of a stack of 2D square images.

Input

• nx::Int must equal to ny (square images only)

Option

• T : datatype of work arrays, defaults to Float32
• method::Symbol : default is :two for 2D interpolation; use :one for 3-pass rotation with 1D interpolation
• nthread::Int : default to Threads.nthreads() The other useful option is 1 when rotating just one image.

plus1di!(mat2d, mat1d)

Non-allocating mat2d[i, :] += mat1d

plus1dj!(mat2d, mat1d)

Non-allocating mat2d[:, j] += mat1d

plus2di!(mat1d, mat2d, i)

Non-allocating mat1d += mat2d[i,:]

plus2dj!(mat1d, mat2d, j)

Non-allocating mat1d += mat2d[:,j]

plus3di!(mat2d, mat3d, i)

Non-allocating mat2d += mat3d[i,:,:]

plus3dj!(mat2d, mat3d, j)

Non-allocating mat2d += mat3d[:,j,:]

plus3dk!(mat2d, mat3d, k)

Non-allocating mat2d += mat3d[:,:,k]

project!(views, image, plan; index)

Project image into multiple views with indexes index (defaults to 1:nview). The 3D views array must be pre-allocated, but need not be initialized.

project!(view, plan, image, buffer_id, viewidx)

SPECT projection of image into a single view with index viewidx. The view must be pre-allocated but need not be initialized to zero.

project!(view, plan, image, viewidx)

SPECT projection of image into a single view with index viewidx. The view must be pre-allocated but need not be initialized to zero.

views = project(image, mumap, psfs, dy; interpmeth, kwargs...)

Convenience method for SPECT forward projector that does all allocation including initializing plan.

In

• image : 3D array (nx,ny,nz)
• mumap : (nx,ny,nz) 3D attenuation map, possibly zeros()
• psfs : 4D PSF array
• dy::RealU : pixel size

Option

• interpmeth : :one or :two

views = project(image, plan ; kwargs...)

Convenience method for SPECT forward projector that allocates and returns views.

psf_gauss( ; ny, px, pz, fwhm_start, fwhm_end, fwhm, fwhm_x, fwhm_z, T)

Create depth-dependent Gaussian PSFs having specified full-width half-maximum (FHWM) values.

Options

• 'ny::Int = 128'
• 'px::Int = 11' (should be odd)
• 'pz::Int = px' (should be odd)
• 'fwhm_start::Real = 1'
• 'fwhm_end::Real = 4'
• 'fwhm::AbstractVector{<:Real} = range(fwhmstart, fwhmend, ny)'
• 'fwhm_x::AbstractVector{<:Real} = fwhm,
• 'fwhmz::AbstractVector{<:Real} = fwhmx'
• 'T::Type == Float32'

Returned psf is [px, pz, ny] where each PSF sums to 1.

rot180!(B::AbstractMatrix, A::AbstractMatrix)

In place version of rot180, returning rotation of A in B.

rot_f90!(output, img, m)

In-place version of rotating an image by 90/180/270 degrees

rot_f90_adj!(output, img, m)
The adjoint of rotating an image by 90/180/270 degrees

rotate_x!(output, img, tan_θ, interp)

rotate_x_adj!(output, img, tan_θ, interp)

rotate_y!(output, img, sin_θ, interp)

rotate_y_adj!(output, img, sin_θ, interp)

rotl90!(B::AbstractMatrix, A::AbstractMatrix)

In place version of rotl90, returning rotation of A in B.

rotr90!(B::AbstractMatrix, A::AbstractMatrix)

In place version of rotr90, returning rotation of A in B.

scale3dj!(mat2d, mat3d, j, s)

Non-allocating mat2d = s * mat3d[:,j,:]

spawner(fun!, nthread::Int, ntask::Int)

Apply fun!(buffer_id, task) for task in 1:ntask where buffer_id is in {1, …, nthread}

A single-thread version of this would simply be fun!.(1, 1:ntask)