Add cube (cuboid), remove MIRT dependence in docs

Project.toml

@@ -9,6 +9,6 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-LazyGrids = "0.2, 0.3"
+LazyGrids = "0.4"
 SpecialFunctions = "1.5, 2"
 julia = "1.6"

docs/Project.toml

@@ -5,7 +5,6 @@ ImageGeoms = "9ee76f2b-840d-4475-b6d6-e485c9297852"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LazyGrids = "7031d0ef-c40d-4431-b2f8-61a8d2f650db"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-MIRT = "7035ae7a-3787-11e9-139a-5545ed3dc201"
 MIRTjim = "170b2178-6dee-4cb0-8729-b3e8b57834cc"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

docs/lit/examples/10-mri-sense.jl

@@ -32,7 +32,6 @@ using ImagePhantoms: phantom, mri_smap_fit, mri_spectra
 using FFTW: fft, fftshift
 using ImageGeoms: embed
 using LazyGrids: ndgrid
-using MIRT: ir_mri_sensemap_sim
 using MIRTjim: jim, prompt
 using Random: seed!
 using Unitful: mm
@@ -83,7 +82,8 @@ dx, dy = fovs ./ (nx,ny)
 x = (-(nx÷2):(nx÷2-1)) * dx
 y = (-(ny÷2):(ny÷2-1)) * dy
 
-# Define Shepp-Logan phantom object, with random complex phases
+# Define Shepp-Logan phantom object,
+# with random complex phases
 # to make it a bit more realistic.
 
 oa = ellipse_parameters(SheppLoganBrainWeb() ; disjoint=true, fovs)
@@ -109,13 +109,39 @@ mask[[1:8;end-8:end],:] .= false
 jim(x, y, mask, "mask")
 
 
-# ### Sensitivity maps
+#=
+### Sensitivity maps
+
+Here we use highly idealized sensitivity maps,
+roughly corresponding to the
+[Biot-Savart law](https://en.wikipedia.org/wiki/Biot-Savart_law)
+for an infinite thin wire,
+as a crude approximation of a
+[birdcage coil](https://en.wikipedia.org/wiki/Radiofrequency_coil).
+One wire is outside the upper right corner,
+the other is outside the left border.
+=#
 
-# Here we use simulated sensitivity maps from MIRT:
+# response at (x,y) to wire at (wx,wy)
+function biot_savart_wire(x, y, wx, wy)
+    phase = cis(atan(y-wy, x-wx))
+    return oneunit(x) / sqrt(sum(abs2, (x-wx, y-wy))) * phase # 1/r falloff
+end
 
 ncoil = 2
-smap = ir_mri_sensemap_sim(dims=(nx,ny), ncoil=ncoil, orbit_start=[0])
-jim(x, y, cfun(smap), "Sensitivity maps raw")
+wire1 = (a,b) -> biot_savart_wire(a, b, maximum(x) + 8dx, maximum(y) + 8dy)
+wire2 = (a,b) -> biot_savart_wire(a, b, minimum(x) - 20dx, zero(dy))
+smap = [wire1.(x, y'), wire2.(x, y')]
+smap[1] *= cis(3π/4) # match coil phases at image center, ala "quadrature phase"
+smap = cat(dims=3, smap...)
+smap /= maximum(abs, smap)
+mag = abs.(smap)
+phase = angle.(smap)
+
+jim(
+ jim(x, y, mag, "|Sensitivity maps raw|"; color=:cividis, ncol=1, prompt=false),
+ jim(x, y, phase, "∠(Sensitivity maps raw)"; color=:hsv, ncol=1, prompt=false),
+)
 
 #=
 Typical sensitivity map estimation methods
@@ -125,7 +151,7 @@ so that the square-root of the sum of squares (SSoS) is unity:
 
 ssos = sqrt.(sum(abs.(smap).^2, dims=ndims(smap))) # SSoS
 ssos = selectdim(ssos, ndims(smap), 1)
-jim(x, y, ssos, "SSoS for ncoil=$ncoil")
+jim(x, y, ssos, "SSoS for ncoil=$ncoil"; color=:cividis, clim=(0,1))
 
 for ic=1:ncoil # normalize
     selectdim(smap, ndims(smap), ic) ./= ssos
@@ -136,29 +162,29 @@ smaps = stacker(smap) # code hereafter expects vector of maps
 jim(x, y, cfun(smaps), "Sensitivity maps (masked and normalized)")
 
 
-# ### Sensitivity map fitting using complex exponentials
-
 #=
+### Sensitivity map fitting using complex exponentials
+
 The `mri_smap_fit` function fits each `smap`
 with a linear combination of complex exponential signals.
 (These signals are not orthogonal due to the `mask`.)
-With frequencies `-7:7/N`, the maximum error is ≤ 0.2%.
+With frequencies `-9:9/N`, the maximum error is ≤ 0.4%.
 =#
 
 deltas = (dx, dy)
-kmax = 7
+kmax = 9
 fit = mri_smap_fit(smaps, embed; mask, kmax, deltas)
 jim(
  jim(x, y, cfun(smaps), "Original maps"; prompt=false, clim=(-1,1)),
  jim(x, y, cfun(fit.smaps), "Fit maps"; prompt=false, clim=(-1,1)),
  jim(x, y, cfun(100 * (fit.smaps - smaps)), "error * 100"; prompt=false),
 )
 
-# The fit coefficients are smaller near `kmax`
+# The fit coefficients are smaller near `±kmax`
 # so probably `kmax` is large enough.
 
-coefs = map(x -> reshape(x,15,15), fit.coefs)
-jim(-kmax:kmax, -kmax:kmax, cfun(coefs), "coefs", prompt=false)
+coefs = map(x -> reshape(x, 2kmax+1, 2kmax+1), fit.coefs)
+jim(-kmax:kmax, -kmax:kmax, cfun(coefs), "Coefficients")
 
 
 # ### Compare FFT with analytical spectra

docs/lit/examples/33-cuboid.jl

@@ -0,0 +1,230 @@
+#---------------------------------------------------------
+# # [Cuboid](@id 33-cuboid)
+#---------------------------------------------------------
+
+#=
+This page illustrates the `Cuboid` shape in the Julia package
+[`ImagePhantoms`](https://github.com/JuliaImageRecon/ImagePhantoms.jl).
+
+This page was generated from a single Julia file:
+[33-cuboid.jl](@__REPO_ROOT_URL__/33-cuboid.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](http://nbviewer.jupyter.org/) here:
+#md # [`33-cuboid.ipynb`](@__NBVIEWER_ROOT_URL__/33-cuboid.ipynb),
+#md # and opened in [binder](https://mybinder.org/) here:
+#md # [`33-cuboid.ipynb`](@__BINDER_ROOT_URL__/33-cuboid.ipynb).
+
+
+# ### Setup
+
+# Packages needed here.
+
+using ImagePhantoms: Cuboid, phantom, radon, spectrum
+import ImagePhantoms as IP
+using ImageGeoms: ImageGeom, axesf
+using MIRTjim: jim, prompt, mid3
+using FFTW: fft, fftshift
+using LazyGrids: ndgrid
+using Unitful: mm, unit, °
+using UnitfulRecipes
+using Plots: plot, plot!, scatter!, default
+using Plots # gif @animate
+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.
+
+isinteractive() ? jim(:prompt, true) : prompt(:draw);
+
+
+# ### Overview
+
+#=
+A basic shape used in constructing 3D digital image phantoms
+is the 3D cuboid,
+specified by its center, width, angle(s) and value.
+All of the methods in `ImagePhantoms` support physical units,
+so we use such units throughout this example.
+(Using units is recommended but not required.)
+
+Here are 3 ways to define a 3D cuboid object,
+using physical units.
+=#
+
+center = (10mm, 8mm, 6mm)
+width = (35mm, 25mm, 15mm)
+ϕ0s = :(π/6) # symbol version for nice plot titles
+ϕ0 = eval(ϕ0s)
+angles = (ϕ0, 0)
+Cuboid([10mm, 8mm, 6mm, 35mm, 25mm, 15mm, π/6, 0, 1.0f0]) # Vector{Number}
+Cuboid( 10mm, 8mm, 6mm, 35mm, 25mm, 15mm, π/6, 0, 1.0f0 ) # 9 arguments
+ob = Cuboid(center, width, angles, 1.0f0) # tuples (recommended use)
+
+
+#=
+### Phantom image using `phantom`
+
+Make a 3D digital image of it using `phantom` and display it.
+We use `ImageGeoms` to simplify the indexing.
+=#
+
+deltas = (1.0mm, 1.1mm, 1.2mm)
+dims = (2^8, 2^8+2, 48)
+offsets = (0.5, 0.5, 0.5) # for FFT spectra later
+ig = ImageGeom( ; dims, deltas, offsets)
+oversample = 3
+img = phantom(axes(ig)..., [ob], oversample)
+p1 = jim(axes(ig), img;
+   title="Cuboid, rotation ϕ=$ϕ0s", xlabel="x", ylabel="y")
+
+
+# The image integral should match the object volume:
+volume = prod(width)
+(sum(img)*prod(ig.deltas), volume)
+
+
+# Show middle slices
+jim(mid3(img), "Middle 3 planes")
+
+
+#=
+### Spectrum using `spectrum`
+
+There are two ways to examine the spectrum of this 3D image:
+* using the analytical Fourier transform of the object via `spectrum`
+* applying the DFT via FFT to the digital image.
+Because the shape has units `mm`, the spectra axes have units cycles/mm.
+Appropriate frequency axes for DFT are provided by `axesf(ig)`.
+=#
+
+vscale = 1 / volume # normalize spectra by volume
+spectrum_exact = spectrum(axesf(ig)..., [ob]) * vscale
+sp = z -> max(log10(abs(z)/oneunit(abs(z))), -6) # log-scale for display
+clim = (-6, 0) # colorbar limit for display
+(xlabel, ylabel) = ("ν₁", "ν₂")
+p2 = jim(axesf(ig), sp.(spectrum_exact), "log10|Spectrum|"; clim, xlabel, ylabel)
+
+
+# Sadly `fft` cannot handle units currently, so this function is a work-around:
+function myfft(x)
+    u = unit(eltype(x))
+    return fftshift(fft(fftshift(x) / u)) * u
+end
+
+spectrum_fft = myfft(img) * prod(ig.deltas) * vscale
+p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
+
+
+# Compare the DFT and analytical spectra to validate the code
+err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
+@assert err < 4e-2
+p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact);
+   title="|Difference| × 10³", xlabel, ylabel)
+jim(p1, p4, p2, p3)
+
+
+#=
+### Parallel-beam projections using `radon`
+
+Compute 2D projection views of the object using `radon`.
+Validate it using the projection-slice theorem aka Fourier-slice theorem.
+=#
+
+pg = ImageGeom((2^8,2^7), (0.6mm,1.0mm), (0.5,0.5)) # projection sampling
+ϕs, θs = (:(π/3), ϕ0s), (:(π/7), :(0))
+ϕ3 = ϕ0 + atan(ob.width[1], ob.width[2])
+θ3 = atan(ob.width[3], sqrt(sum(abs2, ob.width[1:2])))
+ϕ, θ = [eval.(ϕs)..., ϕ3], [eval.(θs)..., θ3]
+proj3 = [radon(axes(pg)..., ϕ[i], θ[i], [ob]) for i in 1:3] # 3 projections
+smax = ob.value * sqrt(sum(abs2, ob.width))
+p5 = jim(axes(pg)..., proj3; xlabel="u", ylabel="v", nrow = 1, title =
+    "Projections at (ϕ,θ) = ($(ϕs[1]), $(θs[1])) and ($(ϕs[2]), $(θs[2]))\n
+    and along long axis")
+
+
+#=
+Because the object has maximum chord length of
+`smax = sqrt(35^2+25^2+15^2)` ≈ 45.6mm,
+and one of the views above was along the corresponding axis,
+the maximum projection value is about that value.
+
+The integral of each projection should match the object volume:
+=#
+((p -> sum(p)*prod(pg.deltas)).(proj3)..., volume)
+
+
+# Look at a set of projections as the views orbit around the object.
+ϕd = 0:6:360
+ϕs = deg2rad.(ϕd) # * Unitful.rad # todo round unitful Unitful.°
+θs = :(π/7)
+θ = eval(θs)
+projs = radon(axes(pg)..., ϕs, [θ], [ob]) # many projection views
+
+if isinteractive()
+    jim(axes(pg)..., projs; title="projection views $(ϕd)")
+else
+    anim = @animate for ip in 1:length(ϕd)
+        jim(axes(pg), projs[:,:,ip,1]; xlabel="u", ylabel="v", prompt=false,
+            title="ϕ=$(ϕd[ip])° θ=$θs", clim = (0,1) .* smax)
+    end
+    gif(anim, "cuboid.gif", fps = 6)
+end
+
+
+#=
+The above sampling generated a parallel-beam projection,
+but one could make a cone-beam projection
+by sampling `(u, v, ϕ, θ)` appropriately.
+See Sinograms.jl.
+=#
+
+
+#=
+### Fourier-slice theorem illustration
+
+Pick one particular view and compare its FFT
+to a slice through the 3D object spectrum.
+=#
+
+ϕs, θs = :(π/3), :(π/7)
+ϕ, θ = eval.((ϕs, θs))
+proj = radon(axes(pg)..., ϕ, θ, [ob])
+p6 = jim(axes(pg), proj; xlabel="u", ylabel="v", prompt=false,
+    title = "Projection at (ϕ,θ) = ($ϕs, $θs)")
+
+e1 = (cos(ϕ), sin(ϕ), 0)
+e3 = (sin(ϕ)*sin(θ), -cos(ϕ)*sin(θ), cos(θ))
+fu, fv = ndgrid(axesf(pg)...)
+ff = vec(fu) * [e1...]' + vec(fv) * [e3...]' # fx,fy,fz for Fourier-slice theorem
+spectrum_slice = spectrum(ob).(ff[:,1], ff[:,2], ff[:,3]) * vscale
+spectrum_slice = reshape(spectrum_slice, pg.dims)
+clim = (-6, 0) # colorbar limit for display
+(xlabel, ylabel) = ("νᵤ", "νᵥ")
+p7 = jim(axesf(pg), sp.(spectrum_slice); prompt=false,
+    title = "log10|Spectrum Slice|", clim, xlabel, ylabel)
+proj_fft = myfft(proj) * prod(pg.deltas) * vscale
+p8 = jim(axesf(pg), sp.(proj_fft); prompt=false,
+     title = "log10|FFT Spectrum|", clim, xlabel, ylabel)
+
+err = maximum(abs, spectrum_slice - proj_fft) / maximum(abs, spectrum_slice)
+@assert err < 2e-3
+p9 = jim(axesf(pg), 1e6*abs.(proj_fft - spectrum_slice);
+    title="Difference × 10⁶", xlabel, ylabel, prompt=false)
+jim(p6, p7, p8, p9)
+
+
+#=
+The good agreement between
+the 2D slice through the 3D analytical spectrum
+and the FFT of the 2D projection view
+validates that `phantom`, `radon`, and `spectrum`
+are all self consistent
+for this object.
+=#

src/ImagePhantoms.jl

@@ -17,6 +17,7 @@ include("triangle.jl")
 # 3d:
 include("ellipsoid.jl")
 include("gauss3.jl")
+include("cuboid.jl")
 
 # phantoms:
 include("shepplogan.jl")