Add 3D gaussian Gauss3 object

docs/lit/examples/04-gauss.jl

@@ -32,7 +32,8 @@ using MIRTjim: jim, prompt
 using FFTW: fft, fftshift
 using Unitful: mm, unit, °
 using UnitfulRecipes
-using Plots: plot, plot!, scatter!, default; default(markerstrokecolor=:auto)
+using Plots: plot, plot!, scatter!, 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.

docs/lit/examples/31-ellipsoid.jl

@@ -32,7 +32,7 @@ using FFTW: fft, fftshift
 using LazyGrids: ndgrid
 using Unitful: mm, unit, °
 using UnitfulRecipes
-using Plots: plot, plot!, scatter!, gui, default
+using Plots: plot, plot!, scatter!, default
 using Plots # gif @animate
 default(markerstrokecolor=:auto)
 
@@ -42,25 +42,27 @@ default(markerstrokecolor=:auto)
 
 isinteractive() ? jim(:prompt, true) : prompt(:draw);
 
+
 # ### Overview
 
 #=
 A basic shape used in constructing 3D digital image phantoms
-is the ellipsoid, specified by its center, radii, angle(s) and value.
+is the ellipsoid,
+specified by its center, radii, 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 an ellipsoid object,
-using parameters specified with physical units.
+using physical units.
 =#
 
 center = (20mm, 10mm, 5mm)
 radii = (25mm, 35mm, 15mm)
 ϕ0s = :(π/6) # symbol version for nice plot titles
 angles = (eval(ϕ0s), 0)
-Ellipsoid([40mm, 20mm, 10mm, 25mm, 35mm, 15mm, π/6, 0, 1.0f0]) # Vector{Number}
-Ellipsoid(40mm, 20mm, 10mm, 25mm, 35mm, 15mm, π/6, 0, 1.0f0) # 9 arguments
+Ellipsoid([20mm, 10mm, 5mm, 25mm, 35mm, 15mm, π/6, 0, 1.0f0]) # Vector{Number}
+Ellipsoid( 20mm, 10mm, 5mm, 25mm, 35mm, 15mm, π/6, 0, 1.0f0 ) # 9 arguments
 ob = Ellipsoid(center, radii, angles, 1.0f0) # tuples (recommended use)
 
 
@@ -81,7 +83,7 @@ p1 = jim(axes(ig), img;
    title="Ellipsoid, rotation ϕ=$ϕ0s", xlabel="x", ylabel="y")
 
 
-# The image integral should match the ellipsoid volume:
+# The image integral should match the object volume:
 volume = 4/3*π*prod(ob.width)
 (sum(img)*prod(ig.deltas), volume)
 
@@ -119,10 +121,10 @@ p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
 
 
 # Compare the DFT and analytical spectra to validate the code
-@assert maximum(abs, spectrum_exact - spectrum_fft) /
-        maximum(abs, spectrum_exact) < 3e-3
-p4 = jim(axesf(ig), 10^4*abs.(spectrum_fft - spectrum_exact);
-   title="Difference × 10⁴", xlabel, ylabel)
+err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
+@assert err < 3e-3
+p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact);
+   title="Difference × 10³", xlabel, ylabel)
 jim(p1, p4, p2, p3)
 
 
@@ -144,9 +146,10 @@ p5 = jim(axes(pg)..., proj2; xlabel="u", ylabel="v", title =
 #=
 Because the ellipsoid has major axis of length 70mm
 and one of the two views above was along that axis,
-the maximum projection value is about 70mm.
+the maximum projection value is about
+70mm.
 
-The integral of each projection should match the ellipsoid volume:
+The integral of each projection should match the object volume:
 =#
 ((p -> sum(p)*prod(pg.deltas)).(proj2)..., volume)
 
@@ -205,10 +208,10 @@ proj_fft = myfft(proj) * prod(pg.deltas) * vscale
 p8 = jim(axesf(pg), sp.(proj_fft); prompt=false,
      title = "log10|FFT Spectrum|", clim, xlabel, ylabel)
 
-@assert maximum(abs, spectrum_slice - proj_fft) /
-        maximum(abs, spectrum_slice) < 1e-3
-p9 = jim(axesf(pg), 10^4*abs.(proj_fft - spectrum_slice);
-    title="Difference × 10⁴", xlabel, ylabel, prompt=false)
+err = maximum(abs, spectrum_slice - proj_fft) / maximum(abs, spectrum_slice)
+@assert err < 1e-3
+p9 = jim(axesf(pg), 1e3*abs.(proj_fft - spectrum_slice);
+    title="Difference × 10³", xlabel, ylabel, prompt=false)
 jim(p6, p7, p8, p9)
 
 
@@ -217,5 +220,6 @@ 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.
+are all self consistent
+for this object.
 =#

docs/lit/examples/34-gauss3.jl

@@ -0,0 +1,226 @@
+#---------------------------------------------------------
+# # [Gauss3](@id 34-gauss3)
+#---------------------------------------------------------
+
+#=
+This page illustrates the `Gauss3` shape in the Julia package
+[`ImagePhantoms`](https://github.com/JuliaImageRecon/ImagePhantoms.jl).
+
+This page was generated from a single Julia file:
+[34-gauss3.jl](@__REPO_ROOT_URL__/34-gauss3.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 # [`34-gauss3.ipynb`](@__NBVIEWER_ROOT_URL__/34-gauss3.ipynb),
+#md # and opened in [binder](https://mybinder.org/) here:
+#md # [`34-gauss3.ipynb`](@__BINDER_ROOT_URL__/34-gauss3.ipynb).
+
+
+# ### Setup
+
+# Packages needed here.
+
+using ImagePhantoms: Gauss3, 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 gaussian,
+specified by its center, fwhm, 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 gaussian object,
+using physical units.
+=#
+
+center = (20mm, 10mm, 5mm)
+width = (25mm, 35mm, 15mm)
+ϕ0s = :(π/6) # symbol version for nice plot titles
+angles = (eval(ϕ0s), 0)
+Gauss3([40mm, 20mm, 2mm, 25mm, 35mm, 12mm, π/6, 0, 1.0f0]) # Vector{Number}
+Gauss3(20mm, 20mm, 2mm, 25mm, 35mm, 12mm, π/6, 0, 1.0f0) # 9 arguments
+ob = Gauss3(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 = 2
+img = phantom(axes(ig)..., [ob], oversample)
+p1 = jim(axes(ig), img;
+   title="Gauss3, rotation ϕ=$ϕ0s", xlabel="x", ylabel="y")
+
+
+# The image integral should match the object volume:
+volume = IP.fwhm2spread(1)^3 * 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 < 1e-3
+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 = (:(π/2), ϕ0s), (:(π/7), :(0))
+ϕ, θ = [eval.(ϕs)...], [eval.(θs)...]
+proj2 = [radon(axes(pg)..., ϕ[i], θ[i], [ob]) for i in 1:2] # 2 projections
+smax = ob.value * IP.fwhm2spread(maximum(ob.width))
+p5 = jim(axes(pg)..., proj2; xlabel="u", ylabel="v", title =
+    "Projections at (ϕ,θ) = ($(ϕs[1]), $(θs[1])) and ($(ϕs[2]), $(θs[2]))")
+
+
+#=
+Because the object has maximum FWHM of 35mm,
+and one of the two views above was along the corresponding axis,
+the maximum projection value is about
+`fwhm2spread(35mm) = 35mm * sqrt(π / log(16))` ≈ 37.25mm.
+
+The integral of each projection should match the object volume:
+=#
+((p -> sum(p)*prod(pg.deltas)).(proj2)..., 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, "gauss3.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 < 1e-5
+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

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

src/gauss2.jl

@@ -82,7 +82,7 @@ Returns function of `(x,y)` for making image.
 """
 function phantom(ob::Object2d{Gauss2})
     return (x,y) -> ob.value * # trick due to fwhm to spread scaling:
-        exp(-π * sum(abs2.(coords(ob, x, y))) / fwhm2spread(1)^2)
+        exp(-π * sum(abs2, coords(ob, x, y)) / fwhm2spread(1)^2)
 end