Document sinogram geometries (#17)

* 4space

* concise unit show

* sino plot units

* docs

* show sg

Author: JeffFessler

docs/lit/examples/02-sino-geom.jl

@@ -26,7 +26,11 @@ This page was generated from a single Julia file:
 
 # Packages needed here.
 
-using Sinograms: SinoPar, SinoMoj, SinoFanArc, SinoFanFlat, rays
+using Unitful: mm, °
+using UnitfulRecipes
+using Plots # must precede 'using Sinograms' for sino_plot_rays to work
+using Sinograms: SinoPar, SinoMoj, SinoFanArc, SinoFanFlat, SinoFan
+using Sinograms: sino_plot_rays
 using MIRTjim: jim, prompt
 using InteractiveUtils: versioninfo
 
@@ -72,25 +76,117 @@ In this package,
 is the type that describes
 a parallel-beam sinogram geometry.
 
-The built-in defaults are helpful.
+The built-in defaults provide helpful reminders about the usage.
 =#
 
 SinoPar()
 
+#=
+This package supports units via the
+[Unitful.jl](https://github.com/PainterQubits/Unitful.jl)
+and using units is recommended
+(but not required).
+
+Here is an example of how to specify
+a parallel-beam geometry.
+Everything is a named keyword argument with sensible default values.
+
+* `orbit` and `orbit_start` must both be unitless (degrees)
+  or have the same units (e.g., degrees or radians).
+* detector spacing `d` and `strip_width`
+  must both be unitless or have the same units (e.g., mm).
+* The projection angles ``ϕ`` are equally space and given by
+  `orbit_start + (0:(nb-1))/nb * orbit`.
+=#
+
+sg = SinoPar( ;
+    nb = 64, # number of radial samples ("bins")
+    na = 30, # number of angular samples
+    d = 2mm, # detector spacing
+    offset = 0.25, # quarter detector offset (unitless)
+    orbit = 180, # angular range (in degrees)
+    orbit_start = 0, # starting angle (in degrees)
+    strip_width = 2mm, # detector width
+)
+
+#=
+The struct `sg` has numerous useful properties;
+type `?SinoGeom` to see the full list.
+
+For example,
+to access the angular samples in degrees
+type `sg.ad`
+=#
 
-## todo: more details
+sg.ad
 
 
+# The following function visualizes the sampling pattern.
+
+sino_plot_rays(sg; ylims=(0,180), yticks=(0:90:180), widen=true, title="Parallel")
+
+#
+prompt()
+
+
+#=
 ## Fan-beam CT with an arc detector (3rd generation CT)
 
-SinoFanArc()
+For a fan-beam geometry,
+the arguments are the same as for `SinoPar`
+with the addition of specifying:
+* `dsd` distance from source to detector
+* `dod` distance from origin (isocenter) to detector
+* `dfs` distance from focal point of detector to source
+  (0 for a 3rd gen arc detector, and `Inf` for a flat detector)
+* `source_offset` for misaligned systems
+   where the ray from the source to the detector center
+   does not intersect the isocenter.
+   Not fully supported; submit an issue if you need this feature.
+
+Here is an example that corresponds to a GE Lightspeed CT system.
+These numbers are published in
+[IEEE T-MI Oct. 2006, p.1272-1283](http://doi.org/10.1109/TMI.2006.882141).
+=#
+
+sg = SinoFanArc( ; nb=888, na=984,
+    d=1.0239mm, offset=1.25, dsd=949.075mm, dod=408.075mm)
+
+
+# Here is a smaller example for plotting the rays.
+
+sg = SinoFanArc( ; nb=64, na=30,
+    d=20mm, offset=0.25, dsd=900mm, dod=400mm)
+sino_plot_rays(sg; ylims=(-50,400), yticks=(0:180:360), widen=true,
+    title="Fan-beam for arc detector")
 
+#
+prompt()
 
+
+#=
 ## Fan-beam CT with a flat detector
 
-SinoFanFlat()
+This geometry is the same as the arc detector
+except that `dfs=Inf`.
+=#
 
+sg = SinoFanFlat( ; nb=64, na=30,
+    d=20mm, offset=0.25, dsd=900mm, dod=400mm)
 
+
+# Here is its sampling plot
+sino_plot_rays(sg; ylims=(-50,400), yticks=(0:180:360), widen=true,
+    title="Fan-beam for flat detector")
+
+#
+prompt()
+
+
+#=
 ## Mojette sampling
+This is a specialized sampling geometry
+that is currently incompletely supported.
+=#
 
 SinoMoj()

matlab/test1.jl

@@ -4,12 +4,12 @@ using Plots
 
 # setup matlab to use MIRT
 if !@isdefined(irtdir)
-	ENV["MATLAB_ROOT"] = "/Applications/freeware/matlab"
+    ENV["MATLAB_ROOT"] = "/Applications/freeware/matlab"
 
-	irtdir = "/Users/fessler/src/matlab/alg"
-	tmp = "addpath('$irtdir')"
-	eval_string(tmp)
-	mat"setup"
+    irtdir = "/Users/fessler/src/matlab/alg"
+    tmp = "addpath('$irtdir')"
+    eval_string(tmp)
+    mat"setup"
 end
 
 fov = 500.
@@ -36,12 +36,12 @@ mat"im_mat = ellipse_im(ig, ell)"
 
 @mget im_mat
 
-plot(jim(im_julia, "julia"), jim(im_mat, "matlab"), jim(im_mat - im_julia)) 
+plot(jim(im_julia, "julia"), jim(im_mat, "matlab"), jim(im_mat - im_julia))
 
 sino_julia = ellipse_sino(sg, ell)
 
 mat"[sino_mat, pos, ang] = ellipse_sino(sg, ell)"
 mat"sizeof(sino_mat)"
 #@mget sino_mat
 
-#plot(jim(sino_julia, "julia"), jim(sino_mat, "matlab"), jim(sino_mat - sino_julia)) 
+#plot(jim(sino_julia, "julia"), jim(sino_mat, "matlab"), jim(sino_mat - sino_julia))

matlab/zwart_powell.jl

@@ -8,12 +8,12 @@ using LazyGrids: ndgrid_array
 
 # setup matlab to use MIRT
 if !@isdefined(irtdir)
-	ENV["MATLAB_ROOT"] = "/Applications/freeware/matlab"
+    ENV["MATLAB_ROOT"] = "/Applications/freeware/matlab"
 
-	irtdir = "/Users/fessler/src/matlab/alg"
-	tmp = "addpath('$irtdir')"
-	eval_string(tmp)
-	mat"setup"
+    irtdir = "/Users/fessler/src/matlab/alg"
+    tmp = "addpath('$irtdir')"
+    eval_string(tmp)
+    mat"setup"
 end
 
 

src/Sinograms.jl

@@ -34,6 +34,7 @@ If no Unitful package loaded, assume `angle` is in degrees and convert to radian
 #to_radians(angle::Real) = deg2rad(Float32(angle))
 to_radians(aa::AbstractArray{T}) where {T <: AbstractFloat} = aa * T(deg2rad(1))
 #to_radians(aa::AbstractArray{T}) where {T <: Real} = aa * deg2rad(1f0) # Float32
+_unit_precision(x::T) where {T <: Number} = T
 
 
 # support Plots iff user has loaded that package

src/fbp2/sino-geom.jl

@@ -766,7 +766,7 @@ function Base.show(io::IO, ::MIME"text/plain", sg::SinoGeom)
     println(io, "$(typeof(sg)) :")
     for f in fieldnames(typeof(sg))
         p = getproperty(sg, f)
-        t = typeof(p)
+        t = (p isa Number) ? _unit_precision(p) : typeof(p)
         println(io, " ", f, "::", t, " ", p)
     end
 end

src/fbp2/sino-plot.jl

@@ -19,16 +19,19 @@ Requires `Plots`.
 function sino_plot_rays(sg::SinoGeom; kwargs...)
     r, phi = rays(sg)
     ad = rad2deg.(phi)
-    ylims = (min(0, minimum(ad)), max(360, maximum(ad)))
-    rmax = ceil(maximum(abs, r)/10, digits=0) * 10
+    unit_a = oneunit(eltype(ad))
+    ylims = (min(zero(unit_a), minimum(ad)),
+             max(360*unit_a, maximum(ad)))
+    unit_r = oneunit(eltype(r))
+    rmax = ceil(maximum(abs, r/unit_r)/10, digits=0) * 10 * unit_r
     xlims = (-1,1) .* rmax
     scatter(
         r, ad ;
         label="",
         markersize=1, markerstrokecolor=:auto,
         markershape=:circle, linewidth=0,
         ylims, xlims, ylabel="ϕ",
-        xticks=(-1:1)*rmax, yticks=[0,360],
+        xticks=(-1:1)*rmax, yticks=[0,360] * unit_a,
         title="$(typeof(sg))",
         kwargs...
     )

src/units.jl

@@ -15,3 +15,5 @@ Convert `Unitful` quantity to radians.
 #    aa * ((one(Float32)*rad) / oneunit(eltype(aa)))
 to_radians(aa::AbstractArray{<: Unitful.Quantity{T}}) where {T <: AbstractFloat} =
     aa * ((one(T)*rad) / oneunit(eltype(aa)))
+
+_unit_precision(x::Unitful.Quantity{T}) where {T <: Number} = "Unit{$T}"

test/sys2/zwart_powell.jl

@@ -22,7 +22,7 @@ using Test: @test, @testset, @test_throws, @inferred
     sino = zwart_powell.(r, ϕ')
 
     p2 = plot()
-	for (i,θ) in enumerate(ϕ)
+    for (i,θ) in enumerate(ϕ)
         plot!(r, sino[:,i], label="ϕ=$θ")
     end
     plot(p1, p2)