From dd335264af675d705487882086effc8a737ab51f Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Tue, 17 Jan 2023 15:19:41 -0500 Subject: [PATCH] https url (#43) * using/import * using/import * https url --- README.md | 2 +- docs/lit/examples/1-overview.jl | 196 ++++++++++++++++++-------------- docs/lit/examples/2-rotate.jl | 141 +++++++++++++---------- docs/lit/examples/3-psf.jl | 103 ++++++++++------- docs/lit/examples/4-mlem.jl | 57 ++++++---- docs/lit/examples/5-2d.jl | 128 ++++++++++++--------- docs/lit/examples/6-dl.jl | 19 ++-- docs/lit/examples/7-osem.jl | 50 +++++--- docs/src/history.md | 4 +- docs/src/index.md | 2 +- 10 files changed, 412 insertions(+), 290 deletions(-) diff --git a/README.md b/README.md index 10f49ae9..e7d33e86 100644 --- a/README.md +++ b/README.md @@ -48,7 +48,7 @@ Tested with Julia ≥ 1.6. [docs-stable-url]: https://JuliaImageRecon.github.io/SPECTrecon.jl/stable [docs-dev-img]: https://img.shields.io/badge/docs-dev-blue.svg [docs-dev-url]: https://JuliaImageRecon.github.io/SPECTrecon.jl/dev -[license-img]: http://img.shields.io/badge/license-MIT-brightgreen.svg?style=flat +[license-img]: https://img.shields.io/badge/license-MIT-brightgreen.svg [license-url]: LICENSE [colprac-img]: https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet [colprac-url]: https://github.com/SciML/ColPrac diff --git a/docs/lit/examples/1-overview.jl b/docs/lit/examples/1-overview.jl index 1d69f630..41f5922e 100644 --- a/docs/lit/examples/1-overview.jl +++ b/docs/lit/examples/1-overview.jl @@ -1,11 +1,24 @@ -#--------------------------------------------------------- -# # [SPECTrecon overview](@id 1-overview) -#--------------------------------------------------------- +#= +# [SPECTrecon overview](@id 1-overview) -# This page gives an overview of the Julia package -# [`SPECTrecon`](https://github.com/JuliaImageRecon/SPECTrecon.jl). +This page gives an overview of the Julia package +[`SPECTrecon`](https://github.com/JuliaImageRecon/SPECTrecon.jl). -# ### Setup +This page was generated from a single Julia file: +[1-overview.jl](@__REPO_ROOT_URL__/1-overview.jl). +=# + +#md # In any such Julia documentation, +#md # you can access the source code +#md # using the "Edit on GitHub" link in the top right. + +#md # The corresponding notebook can be viewed in +#md # [nbviewer](https://nbviewer.org/) here: +#md # [`1-overview.ipynb`](@__NBVIEWER_ROOT_URL__/1-overview.ipynb), +#md # and opened in [binder](https://mybinder.org/) here: +#md # [`1-overview.ipynb`](@__BINDER_ROOT_URL__/1-overview.ipynb). + +# ## Setup # Packages needed here. @@ -22,46 +35,48 @@ using InteractiveUtils: versioninfo 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](https://github.com/JuliaImageRecon/ImagePhantoms.jl), -# but instead we just use two simple cubes. +#= +## 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](https://github.com/JuliaImageRecon/ImagePhantoms.jl), +but instead we just use two simple cubes. +=# nx,ny,nz = 128,128,80 T = Float32 @@ -80,7 +95,7 @@ end jim(mid3(xtrue), "Middle slices of x") -# ### PSF +# ## PSF # Create a synthetic depth-dependent PSF for a single view px = 11 @@ -88,10 +103,12 @@ psf1 = psf_gauss( ; ny, px, fwhm_end = 6) jim(psf1, "PSF for each of $ny planes") -# 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. +#= +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) @@ -103,11 +120,13 @@ size(psfs) plan = plan_psf( ; nx, nz, px) -# ### Basic SPECT forward projection +#= +## Basic SPECT forward projection -# Here is a simple illustration -# of a SPECT forward projection operation. -# (This is a memory inefficient way to do it!) +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 @@ -120,14 +139,16 @@ size(views) jim(views[:,:,1:4:end], "Every 4th of $nview projection views") -# ### Basic SPECT back projection +#= +## Basic SPECT back projection -# This illustrates an "unfiltered backprojection" -# that leads to a very blurry image -# (again, with a simple memory inefficient usage). +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. +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 @@ -141,18 +162,20 @@ back = backproject(views, mumap, psfs, dy) jim(mid3(back), "Back-projection of ytrue") -# ### Memory efficiency +#= +## 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. +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. +Here we illustrate the memory efficient versions +that are recommended for iterative SPECT reconstruction. -# First construction the SPECT plan. +First construction the SPECT plan. +=# #viewangle = (0:(nview-1)) * 2π # default plan = SPECTplan(mumap, psfs, dy; T) @@ -168,16 +191,18 @@ project!(tmp, xtrue, plan) tmp = Array{T}(undef, nx, ny, nz) backproject!(tmp, views, plan) -@assert tmp == back +@assert tmp ≈ back -# ### Using `LinearMapAA` +#= +## 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. +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) @@ -187,7 +212,7 @@ A = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim) # Simple forward and back-projection: @assert A * xtrue == views -@assert A' * views == back +@assert A' * views ≈ back # Mutating version: tmp = Array{T}(undef, nx, nz, nview) @@ -195,19 +220,22 @@ mul!(tmp, A, xtrue) @assert tmp == views tmp = Array{T}(undef, nx, ny, nz) mul!(tmp, A', views) -@assert tmp == back +@assert tmp ≈ back + +#= +## Units -# ### Units +The pixel dimensions `deltas` can (and should!) be values with units. -# The pixel dimensions `deltas` can (and should!) be values with units. +Here is an example ... (todo) +=# -# Here is an example ... (todo) #using UnitfulRecipes #using Unitful: mm -# ### Projection view animation +# ## Projection view animation anim = @animate for i in 1:nview ymax = maximum(views) @@ -219,7 +247,7 @@ end gif(anim, "views.gif", fps = 8) -# ### Reproducibility +# ## Reproducibility # This page was generated with the following version of Julia: diff --git a/docs/lit/examples/2-rotate.jl b/docs/lit/examples/2-rotate.jl index 9d5c6be0..ac2c602e 100644 --- a/docs/lit/examples/2-rotate.jl +++ b/docs/lit/examples/2-rotate.jl @@ -1,11 +1,24 @@ -#--------------------------------------------------------- -# # [SPECTrecon rotation](@id 2-rotate) -#--------------------------------------------------------- +#= +# [SPECTrecon rotation](@id 2-rotate) -# This page explains the image rotation portion of the Julia package -# [`SPECTrecon.jl`](https://github.com/JuliaImageRecon/SPECTrecon.jl). +This page explains the image rotation portion of the Julia package +[`SPECTrecon.jl`](https://github.com/JuliaImageRecon/SPECTrecon.jl). -# ### Setup +This page was generated from a single Julia file: +[2-rotate.jl](@__REPO_ROOT_URL__/2-rotate.jl). +=# + +#md # In any such Julia documentation, +#md # you can access the source code +#md # using the "Edit on GitHub" link in the top right. + +#md # The corresponding notebook can be viewed in +#md # [nbviewer](https://nbviewer.org/) here: +#md # [`2-rotate.ipynb`](@__NBVIEWER_ROOT_URL__/2-rotate.ipynb), +#md # and opened in [binder](https://mybinder.org/) here: +#md # [`2-rotate.ipynb`](@__BINDER_ROOT_URL__/2-rotate.ipynb). + +# ## Setup # Packages needed here. @@ -19,54 +32,56 @@ default(markerstrokecolor=:auto, markersize=3) isinteractive() ? jim(:prompt, true) : prompt(:draw); -# ### Overview - -# The first step in SPECT image forward projection -# is to rotate each slice of a 3D image volume -# to the appropriate view angle. -# In principle -# one could use any of numerous candidate interpolation methods -# for this task. -# However, because emission images are nonnegative -# and maximum likelihood methods -# for SPECT image reconstruction -# exploit that nonnegativity, -# it is desirable to use interpolators -# that preserve nonnegativity. -# This constraint rules out quadratic and higher B-splines, -# including the otherwise attractive cubic B-spline methods. -# On the other hand, nearest-neighbor interpolation -# (equivalent to 0th-order B-splines) -# does not provide adequate image quality. -# This leaves 1st-order interpolation methods -# as the most viable options. - -# This package supports two 1st-order linear interpolators: -# * 2D bilinear interpolation, -# * a 3-pass rotation method based on 1D linear interpolation. - -# Because image rotation is done repeatedly -# (for every slice of both the emission image and the attenuation map, -# for both projection and back-projection, -# and for multiple iterations), -# it is important for efficiency -# to use mutating methods -# rather than to repeatedly make heap allocations. - -# Following other libraries like -# [FFTW.jl](https://github.com/JuliaMath/FFTW.jl), -# the rotation operations herein start with a `plan` -# where work arrays are preallocated -# for subsequent use. -# The `plan` is a `Vector` of `PlanRotate` objects: -# one for each thread. -# (Parallelism is across slices for a 3D image volume.) -# The number of threads defaults to `Threads.nthreads()`. - - -# ### Example - -# Start with a 3D image volume (just 2 slices here for simplicity). +#= +## Overview + +The first step in SPECT image forward projection +is to rotate each slice of a 3D image volume +to the appropriate view angle. +In principle +one could use any of numerous candidate interpolation methods +for this task. +However, because emission images are nonnegative +and maximum likelihood methods +for SPECT image reconstruction +exploit that nonnegativity, +it is desirable to use interpolators +that preserve nonnegativity. +This constraint rules out quadratic and higher B-splines, +including the otherwise attractive cubic B-spline methods. +On the other hand, nearest-neighbor interpolation +(equivalent to 0th-order B-splines) +does not provide adequate image quality. +This leaves 1st-order interpolation methods +as the most viable options. + +This package supports two 1st-order linear interpolators: +* 2D bilinear interpolation, +* a 3-pass rotation method based on 1D linear interpolation. + +Because image rotation is done repeatedly +(for every slice of both the emission image and the attenuation map, +for both projection and back-projection, +and for multiple iterations), +it is important for efficiency +to use mutating methods +rather than to repeatedly make heap allocations. + +Following other libraries like +[FFTW.jl](https://github.com/JuliaMath/FFTW.jl), +the rotation operations herein start with a `plan` +where work arrays are preallocated +for subsequent use. +The `plan` is a `Vector` of `PlanRotate` objects: +one for each thread. +(Parallelism is across slices for a 3D image volume.) +The number of threads defaults to `Threads.nthreads()`. + + +## Example + +Start with a 3D image volume (just 2 slices here for simplicity). +=# T = Float32 # work with single precision to save memory image = zeros(T, 64, 64, 2) @@ -116,10 +131,12 @@ jim(result1, "Rotated image by π/6 (3-pass 1D)") jim(result1 - result2, "Difference images") -# ### Adjoint +#= +## Adjoint -# To ensure adjoint consistency between SPECT forward- and back-projection, -# there is also an adjoint routine: +To ensure adjoint consistency between SPECT forward- and back-projection, +there is also an adjoint routine: +=# adj2 = similar(result2) imrotate_adj!(adj2, result2, π/6, plan2) @@ -135,11 +152,13 @@ jim(adj1, "Adjoint image rotation (3-pass 1D)") # so one does not expect the output here to match the original image! -# ### LinearMap +#= +## LinearMap -# One can form a linear map corresponding to image rotation using `LinearMapAA`. -# An operator like this may be useful -# as part of a motion-compensated image reconstruction method. +One can form a linear map corresponding to image rotation using `LinearMapAA`. +An operator like this may be useful +as part of a motion-compensated image reconstruction method. +=# using LinearMapsAA: LinearMapAA diff --git a/docs/lit/examples/3-psf.jl b/docs/lit/examples/3-psf.jl index e7a0dff8..c84b92f5 100644 --- a/docs/lit/examples/3-psf.jl +++ b/docs/lit/examples/3-psf.jl @@ -1,11 +1,24 @@ -#--------------------------------------------------------- -# # [SPECTrecon PSF](@id 3-psf) -#--------------------------------------------------------- +#= +# [SPECTrecon PSF](@id 3-psf) -# This page explains the PSF portion of the Julia package -# [`SPECTrecon.jl`](https://github.com/JuliaImageRecon/SPECTrecon.jl). +This page explains the PSF portion of the Julia package +[`SPECTrecon.jl`](https://github.com/JuliaImageRecon/SPECTrecon.jl). -# ### Setup +This page was generated from a single Julia file: +[3-psf.jl](@__REPO_ROOT_URL__/3-psf.jl). +=# + +#md # In any such Julia documentation, +#md # you can access the source code +#md # using the "Edit on GitHub" link in the top right. + +#md # The corresponding notebook can be viewed in +#md # [nbviewer](https://nbviewer.org/) here: +#md # [`3-psf.ipynb`](@__NBVIEWER_ROOT_URL__/3-psf.ipynb), +#md # and opened in [binder](https://mybinder.org/) here: +#md # [`3-psf.ipynb`](@__BINDER_ROOT_URL__/3-psf.ipynb). + +# ## Setup # Packages needed here. @@ -19,36 +32,38 @@ default(markerstrokecolor=:auto, markersize=3) isinteractive() ? jim(:prompt, true) : prompt(:draw); -# ### Overview +#= +## Overview -# After rotating the image and the attenuation map, -# second step in SPECT image forward projection -# is to apply depth-dependent point spread function (PSF). -# Each (rotated) image plane -# is a certain distance from the SPECT detector -# and must be convolved with the 2D PSF appropriate -# for that plane. +After rotating the image and the attenuation map, +second step in SPECT image forward projection +is to apply depth-dependent point spread function (PSF). +Each (rotated) image plane +is a certain distance from the SPECT detector +and must be convolved with the 2D PSF appropriate +for that plane. -# Because SPECT has relatively poor spatial resolution, -# the PSF is usually fairly wide, -# so convolution using FFT operations -# is typically more efficient -# than direct spatial convolution. +Because SPECT has relatively poor spatial resolution, +the PSF is usually fairly wide, +so convolution using FFT operations +is typically more efficient +than direct spatial convolution. -# Following other libraries like -# [FFTW.jl](https://github.com/JuliaMath/FFTW.jl), -# the PSF operations herein start with a `plan` -# where work arrays are preallocated -# for subsequent use. -# The `plan` is a `Vector` of `PlanPSF` objects: -# one for each thread. -# (Parallelism is across planes for a 3D image volume.) -# The number of threads defaults to `Threads.nthreads()`. +Following other libraries like +[FFTW.jl](https://github.com/JuliaMath/FFTW.jl), +the PSF operations herein start with a `plan` +where work arrays are preallocated +for subsequent use. +The `plan` is a `Vector` of `PlanPSF` objects: +one for each thread. +(Parallelism is across planes for a 3D image volume.) +The number of threads defaults to `Threads.nthreads()`. -# ### Example +## Example -# Start with a 3D image volume. +Start with a 3D image volume. +=# T = Float32 # work with single precision to save memory nx = 32 @@ -68,11 +83,13 @@ psf = repeat(psf_gauss( ; ny=nx, px), 1, 1, 1, nview) jim(psf, "PSF for each of $nx planes") -# Now plan the PSF modeling -# by specifying -# * the image size (must be square) -# * the PSF size: must be `px × pz × ny × nview` -# * the `DataType` used for the work arrays. +#= +Now plan the PSF modeling +by specifying +* the image size (must be square) +* the PSF size: must be `px × pz × ny × nview` +* the `DataType` used for the work arrays. +=# plan = plan_psf( ; nx, nz, px, T) @@ -89,10 +106,12 @@ fft_conv!(result, image, psf[:,:,:,1], plan) # mutates the first argument jim(result, "After applying PSF") -# ### Adjoint +#= +## Adjoint -# To ensure adjoint consistency between SPECT forward- and back-projection, -# there is also an adjoint routine: +To ensure adjoint consistency between SPECT forward- and back-projection, +there is also an adjoint routine: +=# adj = similar(result) fft_conv_adj!(adj, result, psf[:,:,:,1], plan) @@ -103,10 +122,12 @@ jim(adj, "Adjoint of PSF modeling") # so one does not expect the output here to match the original image! -# ### LinearMap +#= +## LinearMap -# One can form a linear map corresponding to PSF modeling using `LinearMapAA`. -# Perhaps the main purpose is simply for verifying adjoint correctness. +One can form a linear map corresponding to PSF modeling using `LinearMapAA`. +Perhaps the main purpose is simply for verifying adjoint correctness. +=# using LinearMapsAA: LinearMapAA diff --git a/docs/lit/examples/4-mlem.jl b/docs/lit/examples/4-mlem.jl index 607d5ad2..6eca9726 100644 --- a/docs/lit/examples/4-mlem.jl +++ b/docs/lit/examples/4-mlem.jl @@ -1,11 +1,24 @@ -#--------------------------------------------------------- -# # [SPECTrecon ML-EM](@id 4-mlem) -#--------------------------------------------------------- +#= +# [SPECTrecon ML-EM](@id 4-mlem) -# This page illustrates ML-EM reconstruction with the Julia package -# [`SPECTrecon`](https://github.com/JuliaImageRecon/SPECTrecon.jl). +This page illustrates ML-EM reconstruction with the Julia package +[`SPECTrecon`](https://github.com/JuliaImageRecon/SPECTrecon.jl). -# ### Setup +This page was generated from a single Julia file: +[4-mlem.jl](@__REPO_ROOT_URL__/4-mlem.jl). +=# + +#md # In any such Julia documentation, +#md # you can access the source code +#md # using the "Edit on GitHub" link in the top right. + +#md # The corresponding notebook can be viewed in +#md # [nbviewer](https://nbviewer.org/) here: +#md # [`4-mlem.ipynb`](@__NBVIEWER_ROOT_URL__/4-mlem.ipynb), +#md # and opened in [binder](https://mybinder.org/) here: +#md # [`4-mlem.ipynb`](@__BINDER_ROOT_URL__/4-mlem.ipynb). + +# ## Setup # Packages needed here. @@ -18,13 +31,15 @@ using Plots: scatter, plot!, default; default(markerstrokecolor=:auto) isinteractive() ? jim(:prompt, true) : prompt(:draw); -# ### Overview +#= +## Overview -# Maximum-likelihood expectation-maximization (ML-EM) -# is a classic algorithm for performing SPECT image reconstruction. +Maximum-likelihood expectation-maximization (ML-EM) +is a classic algorithm for performing SPECT image reconstruction. +=# -# ### Simulation data +# ## Simulation data nx,ny,nz = 64,64,50 T = Float32 @@ -43,7 +58,7 @@ end jim(mid3(xtrue), "Middle slices of xtrue") -# ### PSF +# ## PSF # Create a synthetic depth-dependent PSF for a single view px = 11 @@ -51,17 +66,19 @@ psf1 = psf_gauss( ; ny, px) jim(psf1, "PSF for each of $ny planes") -# 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. +#= +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) -# ### SPECT system model using `LinearMapAA` +# ## SPECT system model using `LinearMapAA` dy = 8 # transaxial pixel size in mm mumap = zeros(T, size(xtrue)) # zero μ-map just for illustration here @@ -76,7 +93,7 @@ odim = (nx,nz,nview) A = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim) -# ### Basic Expectation-Maximization (EM) algorithm +# ## Basic Expectation-Maximization (EM) algorithm # Noisy data using Distributions: Poisson @@ -93,21 +110,21 @@ end jim(ynoisy, "$nview noisy projection views") -# ### ML-EM algorithm - basic version +# ## ML-EM algorithm - basic version x0 = ones(T, nx, ny, nz) # initial uniform image niter = 30 if !@isdefined(xhat1) xhat1 = mlem(x0, ynoisy, background, A; niter) end +size(xhat1) -# This preferable ML-EM version preallocates the output `xhat2` +# This preferable ML-EM version preallocates the output `xhat2`: if !@isdefined(xhat2) xhat2 = copy(x0) mlem!(xhat2, x0, ynoisy, background, A; niter) end - @assert xhat1 ≈ xhat2 jim(mid3(xhat2), "ML-EM at $niter iterations") diff --git a/docs/lit/examples/5-2d.jl b/docs/lit/examples/5-2d.jl index 5229e476..49de5a0e 100644 --- a/docs/lit/examples/5-2d.jl +++ b/docs/lit/examples/5-2d.jl @@ -1,10 +1,23 @@ -#--------------------------------------------------------- -# # [SPECTrecon 2D use](@id 5-2d) -#--------------------------------------------------------- +#= +# [SPECTrecon 2D use](@id 5-2d) -# This page describes how to perform 2D SPECT forward and back-projection -# using the Julia package -# [`SPECTrecon`](https://github.com/JuliaImageRecon/SPECTrecon.jl). +This page describes how to perform 2D SPECT forward and back-projection +using the Julia package +[`SPECTrecon`](https://github.com/JuliaImageRecon/SPECTrecon.jl). + +This page was generated from a single Julia file: +[5-2d.jl](@__REPO_ROOT_URL__/5-2d.jl). +=# + +#md # In any such Julia documentation, +#md # you can access the source code +#md # using the "Edit on GitHub" link in the top right. + +#md # The corresponding notebook can be viewed in +#md # [nbviewer](https://nbviewer.org/) here: +#md # [`5-2d.ipynb`](@__NBVIEWER_ROOT_URL__/5-2d.ipynb), +#md # and opened in [binder](https://mybinder.org/) here: +#md # [`5-2d.ipynb`](@__BINDER_ROOT_URL__/5-2d.ipynb). # ### Setup @@ -24,23 +37,25 @@ using Plots: plot, default; default(markerstrokecolor=:auto) isinteractive() ? jim(:prompt, true) : prompt(:draw); -# ### Overview +#= +## Overview -# Real SPECT systems are inherently 3D imaging systems, -# but for the purpose of prototyping algorithms -# it can be useful to work with 2D simulations. +Real SPECT systems are inherently 3D imaging systems, +but for the purpose of prototyping algorithms +it can be useful to work with 2D simulations. -# Currently, "2D" here means a 3D array with `nz=1`, -# i.e., a single slice. -# The key to working with a single slice -# is that the package allows the PSFs -# to have rectangular support `px × pz` -# where `pz = 1`, i.e., no blur along the axial (z) direction. +Currently, "2D" here means a 3D array with `nz=1`, +i.e., a single slice. +The key to working with a single slice +is that the package allows the PSFs +to have rectangular support `px × pz` +where `pz = 1`, i.e., no blur along the axial (z) direction. -# ### Example +## Example -# Start with a simple 2D digital phantom. +Start with a simple 2D digital phantom. +=# T = Float32 nx,ny,nz = 128,128,1 @@ -64,21 +79,25 @@ plot(-hx:hx, tmp[:,[1:9:end-10;end]], markershape=:o, label="", prompt() -# 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. +#= +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 = 80 psfs = repeat(psf1, 1, 1, 1, nview) size(psfs) -# ### Basic SPECT forward projection +#= +## Basic SPECT forward projection -# Here is a simple illustration -# of a SPECT forward projection operation. -# (This is a memory inefficient way to do it!) +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 @@ -92,14 +111,17 @@ size(sino) jim(sino, "Sinogram") -# ### Basic SPECT back projection +#= +## Basic SPECT back projection + +This illustrates an "unfiltered backprojection" +that leads to a very blurry image +(again, with a simple memory inefficient usage). -# 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. +=# -# First, back-project two "rays" -# to illustrate the depth-dependent PSF. sino1 = zeros(T, nx, nview) sino1[nx÷2, nview÷5] = 1 sino1[nx÷2, 1] = 1 @@ -114,18 +136,20 @@ back = backproject(views, mumap, psfs, dy) jim(back, "Back-projection of ytrue") -# ### Memory efficiency +#= +## 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. +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. +Here we illustrate the memory efficient versions +that are recommended for iterative SPECT reconstruction. -# First construction the SPECT plan. +First construction the SPECT plan. +=# #viewangle = (0:(nview-1)) * 2π # default plan = SPECTplan(mumap, psfs, dy; T) @@ -141,16 +165,18 @@ project!(tmp, xtrue, plan) tmp = Array{T}(undef, nx, ny, nz) backproject!(tmp, views, plan) -@assert tmp == back +@assert tmp ≈ back -# ### Using `LinearMapsAA` +#= +## Using `LinearMapsAA` -# 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. +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) @@ -160,7 +186,7 @@ A = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim) # Simple forward and back-projection: @assert A * xtrue == views -@assert A' * views == back +@assert A' * views ≈ back # Mutating version: tmp = Array{T}(undef, nx, nz, nview) @@ -168,10 +194,10 @@ mul!(tmp, A, xtrue) @assert tmp == views tmp = Array{T}(undef, nx, ny, nz) mul!(tmp, A', views) -@assert tmp == back +@assert tmp ≈ back -# ### Gram matrix impulse response +# ## Gram matrix impulse response points = zeros(T, nx, ny, nz) points[nx÷2,ny÷2,1] = 1 diff --git a/docs/lit/examples/6-dl.jl b/docs/lit/examples/6-dl.jl index 22b35925..c60daa09 100644 --- a/docs/lit/examples/6-dl.jl +++ b/docs/lit/examples/6-dl.jl @@ -1,8 +1,6 @@ -#--------------------------------------------------------- -# # [SPECTrecon deep learning use](@id 6-dl) -#--------------------------------------------------------- - #= +# [SPECTrecon deep learning use](@id 6-dl) + This page describes how to end-to-end train unrolled deep learning algorithms using the Julia package [`SPECTrecon`](https://github.com/JuliaImageRecon/SPECTrecon.jl). @@ -16,7 +14,7 @@ This page was generated from a single Julia file: #md # using the "Edit on GitHub" link in the top right. #md # The corresponding notebook can be viewed in -#md # [nbviewer](http://nbviewer.jupyter.org/) here: +#md # [nbviewer](https://nbviewer.org/) here: #md # [`6-dl.ipynb`](@__NBVIEWER_ROOT_URL__/6-dl.ipynb), #md # and opened in [binder](https://mybinder.org/) here: #md # [`6-dl.ipynb`](@__BINDER_ROOT_URL__/6-dl.ipynb). @@ -27,20 +25,19 @@ This page was generated from a single Julia file: # Packages needed here. using LinearAlgebra: norm, mul! -using SPECTrecon +using SPECTrecon: SPECTplan, project!, backproject!, psf_gauss, mlem! using MIRTjim: jim, prompt using Plots: default; default(markerstrokecolor=:auto) -using Zygote using ZygoteRules: @adjoint using Flux: Chain, Conv, SamePad, relu, params, unsqueeze -using Flux # apparently needed for BSON @load -using NNlib +import Flux # apparently needed for BSON @load +import NNlib using LinearMapsAA: LinearMapAA using Distributions: Poisson using BSON: @load, @save import BSON # load using InteractiveUtils: versioninfo -using Downloads: download +import Downloads # download # 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. @@ -294,7 +291,7 @@ then you will see that the CNN reduces the NRMSE. A more thorough investigation would compare the CNN approach to a suitably optimized regularized approach; -see [http://doi.org/10.1109/EMBC46164.2021.9630985](http://doi.org/10.1109/EMBC46164.2021.9630985). +see [https://doi.org/10.1109/EMBC46164.2021.9630985](https://doi.org/10.1109/EMBC46164.2021.9630985). =# diff --git a/docs/lit/examples/7-osem.jl b/docs/lit/examples/7-osem.jl index b291fd1b..4dba8648 100644 --- a/docs/lit/examples/7-osem.jl +++ b/docs/lit/examples/7-osem.jl @@ -1,24 +1,36 @@ -#--------------------------------------------------------- -# # [SPECTrecon OS-EM](@id 7-osem) -#--------------------------------------------------------- - #= +# [SPECTrecon OS-EM](@id 7-osem) + This page illustrates Ordered-subset expectation-maximization (OS-EM) image reconstruction with the Julia package [`SPECTrecon`](https://github.com/JuliaImageRecon/SPECTrecon.jl). + +This page was generated from a single Julia file: +[7-osem.jl](@__REPO_ROOT_URL__/7-osem.jl). =# +#md # In any such Julia documentation, +#md # you can access the source code +#md # using the "Edit on GitHub" link in the top right. + +#md # The corresponding notebook can be viewed in +#md # [nbviewer](https://nbviewer.org/) here: +#md # [`7-osem.ipynb`](@__NBVIEWER_ROOT_URL__/7-osem.ipynb), +#md # and opened in [binder](https://mybinder.org/) here: +#md # [`7-osem.ipynb`](@__BINDER_ROOT_URL__/7-osem.ipynb). + # ### Setup # Packages needed here. -using SPECTrecon -using MIRTjim: jim, prompt -using Plots: scatter, plot!, default; default(markerstrokecolor=:auto) +using SPECTrecon: psf_gauss, SPECTplan, project!, backproject!, Ablock +using SPECTrecon: osem, osem!, mlem! using LinearMapsAA: LinearMapAA, LinearMapAO using LinearAlgebra: mul! using Distributions: Poisson +using MIRTjim: jim, prompt +using Plots: scatter, plot!, default; default(markerstrokecolor=:auto) # 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. @@ -26,17 +38,17 @@ using Distributions: Poisson isinteractive() ? jim(:prompt, true) : prompt(:draw); #= -### Overview +## Overview Ordered-subset expectation-maximization (OS-EM) is a commonly used algorithm for performing SPECT image reconstruction because of its favorable combination of image quality and speed. See -[Hudson and Larkin, 1994](http://doi.org/10.1109/42.363108). +[Hudson and Larkin, 1994](https://doi.org/10.1109/42.363108). =# -# ### Simulation data +# ## Simulation data nx,ny,nz = 64,64,50 T = Float32 @@ -55,7 +67,7 @@ end jim(mid3(xtrue), "Middle slices of xtrue") -# ### PSF +# ## PSF # Create a synthetic depth-dependent PSF for a single view px = 11 @@ -63,17 +75,19 @@ psf1 = psf_gauss( ; ny, px) jim(psf1, "PSF for each of $ny planes"; ratio=1) -# 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. +#= +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) -# ### SPECT system model using `LinearMapAA` +# ## SPECT system model using `LinearMapAA` dy = 8 # transaxial pixel size in mm mumap = zeros(T, size(xtrue)) # zero μ-map just for illustration here @@ -99,7 +113,7 @@ end jim(ynoisy, "$nview noisy projection views") -# ### OS-EM algorithm - basic version +# ## OS-EM algorithm - basic version x0 = ones(T, nx, ny, nz) # initial uniform image niter = 8 @@ -119,7 +133,7 @@ end @assert xhat1 ≈ xhat2 -# ### Compare with ML-EM +# ## Compare with ML-EM # Run 30 iterations of ML-EM algorithm. niter_mlem = 30 diff --git a/docs/src/history.md b/docs/src/history.md index 4fca9c71..36fe82eb 100644 --- a/docs/src/history.md +++ b/docs/src/history.md @@ -5,7 +5,7 @@ is based on the 1992 paper by GL Zeng & GT Gullberg "Frequency domain implementation of the three-dimensional geometric point response correction in SPECT imaging" -[(DOI)](http://doi.org/10.1109/23.173222). +[(DOI)](https://doi.org/10.1109/23.173222). Historically @@ -49,7 +49,7 @@ the work was extended to consider blob basis functions, leading to a -[2004 comparison paper](http://doi.org/10.1088/0031-9155/49/11/003). +[2004 comparison paper](https://doi.org/10.1088/0031-9155/49/11/003). Somewhere during that period the 3D SPECT projector / backprojector diff --git a/docs/src/index.md b/docs/src/index.md index 38d7aa9b..09f90c4a 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -22,7 +22,7 @@ by GL Zeng & GT Gullberg "Frequency domain implementation of the three-dimensional geometric point response correction in SPECT imaging" (IEEE Tr. on Nuclear Science, 39(5-1):1444-53, Oct 1992) -[(DOI)](http://doi.org/10.1109/23.173222). +[(DOI)](https://doi.org/10.1109/23.173222). The forward projection method works as follows. * The emission image and the attenuation map