Use helpers instead of getproperty (#47)

* Improve short CBCT demo

* Use getfield for show

* WIP to remove getproperty

* Use helpers

* update docs

* view weight

* Use parker_weight

* Coverage

* Simply fbp_ramp

* Fix ds bug

* Test _tau

* v0.4.0

* Use _ds

* Fix ramp

* Move _tau, improve coverage

* Refine DocTestSetup for import

* Setup view_weights, improve type inference, typeof

* Test _angle_weights

* Update logo code

* Update tests

* Hide jim

* fft units

* Type stable fbp_sino_filt

* Remove helper.jl

* window type inference

* fft_filter units and type inference

* fbp2 type inference

* Remove sino_filt

* Test stack of sinos

Author: JeffFessler

Project.toml

@@ -1,7 +1,7 @@
 name = "Sinograms"
 uuid = "02a14def-c6e6-4ab0-b2df-0ab64bc8cdd7"
 authors = ["fessler <[email protected]>"]
-version = "0.3.0"
+version = "0.4.0"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"

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

@@ -27,7 +27,7 @@ This page was generated from a single Julia file:
 using Unitful: mm, °
 using Plots # these 2 must precede 'using Sinograms' for sino_plot_rays to work
 using Sinograms: SinoPar, SinoMoj, SinoFanArc, SinoFanFlat, SinoFan
-using Sinograms: sino_plot_rays, rays
+using Sinograms: sino_plot_rays, rays, angles
 using MIRTjim: jim, prompt
 using InteractiveUtils: versioninfo
 
@@ -95,7 +95,7 @@ Everything is a named keyword argument with sensible default values.
   `orbit_start + (0:(nb-1))/nb * orbit`.
 =#
 
-sg = SinoPar( ;
+rg = SinoPar( ;
     nb = 64, # number of radial samples ("bins")
     na = 30, # number of angular samples
     d = 2mm, # detector spacing
@@ -105,20 +105,20 @@ sg = SinoPar( ;
 )
 
 #=
-The struct `sg` has numerous useful properties;
+The struct `rg` has numerous useful properties;
 type `?SinoGeom` to see the full list.
 
 For example,
 to access the angular samples in degrees
-type `sg.ad`
+type `angles(rg)`
 =#
 
-sg.ad
+angles(rg)
 
 
 # The following function visualizes the sampling pattern.
 
-sino_plot_rays(sg; ylims=(0,180), yticks=(0:90:180), widen=true, title="Parallel")
+sino_plot_rays(rg; ylims=(0,180), yticks=(0:90:180), widen=true, title="Parallel")
 
 #
 prompt()
@@ -144,15 +144,15 @@ 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,
+rg = 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,
+rg = 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,
+sino_plot_rays(rg; ylims=(-50,400), yticks=(0:180:360), widen=true,
     title="Fan-beam for arc detector")
 
 #
@@ -166,12 +166,12 @@ This geometry is the same as the arc detector
 except that `dfs=Inf`.
 =#
 
-sg = SinoFanFlat( ; nb=64, na=30,
+rg = 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,
+sino_plot_rays(rg; ylims=(-50,400), yticks=(0:180:360), widen=true,
     title="Fan-beam for flat detector")
 
 #
@@ -184,10 +184,10 @@ This is a specialized sampling geometry
 that is currently incompletely supported.
 =#
 
-sg = SinoMoj( ; nb=60, na=30)
+rg = SinoMoj( ; nb=60, na=30)
 
 # Here is its sampling plot
-sino_plot_rays(sg; ylims=(0,180), yticks=(0:90:180), widen=true,
+sino_plot_rays(rg; ylims=(0,180), yticks=(0:90:180), widen=true,
     title="Mojette sampling")
 
 #=
@@ -212,7 +212,7 @@ plot_ray!(r, ϕ) = plot!(
     color = :blue,
 )
 ia = 4 # pick an angle
-i = rays(sg) # iterator
+i = rays(rg) # iterator
 rs = [i[1] for i in i] # radial samples
 ϕs = [i[2] for i in i] # projection angle
 r = rs[:,ia]

docs/lit/examples/03-parallel-beam.jl

@@ -24,7 +24,7 @@ This page was generated from a single Julia file:
 
 # Packages needed here.
 
-using Sinograms: SinoPar, rays, plan_fbp, fbp
+using Sinograms: SinoPar, rays, plan_fbp, fbp, fbp_sino_filter
 using ImageGeoms: ImageGeom, fovs, MaskCircle
 using ImagePhantoms: SheppLogan, shepp_logan, radon, phantom
 using Unitful: mm
@@ -55,22 +55,22 @@ but units are optional.
 ig = ImageGeom(MaskCircle(); dims=(128,126), deltas = (2mm,2mm) )
 
 # Use `SinoPar` to define the sinogram geometry.
-sg = SinoPar( ; nb = 130, d = 2mm, na = 100)
+rg = SinoPar( ; nb = 130, d = 2mm, na = 100)
 
 # Ellipse parameters for Shepp-Logan phantom:
 μ = 0.01 / mm # typical linear attenuation coefficient
 ob = shepp_logan(SheppLogan(); fovs = fovs(ig), u = (1, 1, μ))
 
 # Radon transform of Shepp-Logan phantom:
-sino = radon(rays(sg), ob)
-jim(axes(sg), sino; title="Shepp-Logan sinogram", xlabel="r", ylabel="ϕ")
+sino = radon(rays(rg), ob)
+jim(axes(rg), sino; title="Shepp-Logan sinogram", xlabel="r", ylabel="ϕ")
 
 ## Image reconstruction via FBP
 # Here we start with a "plan",
 # which would save work if we were reconstructing many images.
 
-plan = plan_fbp(sg, ig)
-fbp_image, sino_filt = fbp(plan, sino)
+plan = plan_fbp(rg, ig)
+fbp_image = fbp(plan, sino)
 
 
 # A narrow color window is needed to see the soft tissue structures:
@@ -80,3 +80,7 @@ jim(axes(ig), fbp_image, "FBP image"; clim)
 # For comparison, here is the ideal phantom image
 true_image = phantom(axes(ig)..., ob, 2)
 jim(axes(ig)..., true_image, "True phantom image"; clim)
+
+# For fun, here is the filtered sinogram:
+sino_filt = fbp_sino_filter(sino, plan.filter)
+jim(axes(rg), sino_filt; title="Filtered sinogram", xlabel="r", ylabel="ϕ")

docs/lit/examples/04-fan-arc.jl

@@ -64,11 +64,11 @@ but units are optional.
 ig = ImageGeom(MaskCircle(); dims=(128,126), deltas = (2mm,2mm) )
 
 # Use `SinoFanArc` to define the sinogram geometry.
-sg = SinoFanArc( ; nb = 130, d = 3.2mm, na = 100, dsd = 400mm, dod = 140mm)
+rg = SinoFanArc( ; nb = 130, d = 3.2mm, na = 100, dsd = 400mm, dod = 140mm)
 
 # Examine the geometry to verify the FOV:
 jim(axes(ig), ig.mask; prompt=false)
-sino_geom_plot!(sg, ig)
+sino_geom_plot!(rg, ig)
 
 #
 prompt()
@@ -78,8 +78,8 @@ prompt()
 ob = shepp_logan(SheppLogan(); fovs = fovs(ig), u = (1, 1, μ))
 
 # Arc fan-beam sinogram for Shepp-Logan phantom:
-sino = radon(rays(sg), ob)
-jim(axes(sg), sino; title="Shepp-Logan sinogram", xlabel="r", ylabel="ϕ")
+sino = radon(rays(rg), ob)
+jim(axes(rg), sino; title="Shepp-Logan sinogram", xlabel="r", ylabel="ϕ")
 
 
 #=
@@ -88,8 +88,8 @@ Here we start with a "plan",
 which would save work if we were reconstructing many images.
 =#
 
-plan = plan_fbp(sg, ig)
-fbp_image, sino_filt = fbp(plan, sino)
+plan = plan_fbp(rg, ig)
+fbp_image = fbp(plan, sino)
 
 
 # A narrow color window is needed to see the soft tissue structures:

docs/lit/examples/05-fan-flat.jl

@@ -60,11 +60,11 @@ but units are optional.
 ig = ImageGeom(MaskCircle(); dims=(128,126), deltas = (0.2cm,0.2cm) )
 
 # Use `SinoFanFlat` to define the sinogram geometry.
-sg = SinoFanFlat( ; nb = 130, d = 0.3cm, na = 100, dsd = 50cm, dod = 14cm)
+rg = SinoFanFlat( ; nb = 130, d = 0.3cm, na = 100, dsd = 50cm, dod = 14cm)
 
 # Examine the geometry to verify the FOV:
 jim(axes(ig), ig.mask; prompt=false)
-sino_geom_plot!(sg, ig)
+sino_geom_plot!(rg, ig)
 
 #
 prompt()
@@ -74,8 +74,8 @@ prompt()
 ob = shepp_logan(SheppLogan(); fovs = fovs(ig), u = (1, 1, μ))
 
 # Arc fan-beam sinogram for Shepp-Logan phantom:
-sino = radon(rays(sg), ob)
-jim(axes(sg), sino; title="Shepp-Logan sinogram", xlabel="r", ylabel="ϕ")
+sino = radon(rays(rg), ob)
+jim(axes(rg), sino; title="Shepp-Logan sinogram", xlabel="r", ylabel="ϕ")
 
 
 #=
@@ -85,8 +85,8 @@ which would save work if we were reconstructing many images.
 For illustration we include `Hamming` window.
 =#
 
-plan = plan_fbp(sg, ig; window = Window(Hamming(), 1.0))
-fbp_image, sino_filt = fbp(plan, sino)
+plan = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))
+fbp_image = fbp(plan, sino)
 
 
 # A narrow color window is needed to see the soft tissue structures:

docs/lit/examples/06-fan-short.jl

@@ -74,14 +74,14 @@ Note that even though we specify `na = 100`
 we end up with `na = 67` views
 because of the `:short` option.
 =#
-sg = SinoFanArc( :short, ;
+rg = SinoFanArc( :short, ;
     nb = 130, d = 3.2mm, na = 100, dsd = 400mm, dod = 140mm,
 )
 
 # Examine the geometry to verify the FOV.
 # The `na=67` blue dots show the `:short` scan.
 jim(axes(ig), ig.mask; prompt=false)
-sino_geom_plot!(sg, ig)
+sino_geom_plot!(rg, ig)
 
 #
 prompt()
@@ -91,8 +91,8 @@ prompt()
 ob = shepp_logan(SheppLogan(); fovs = fovs(ig), u = (1, 1, μ))
 
 # Short arc fan-beam sinogram for Shepp-Logan phantom:
-sino = radon(rays(sg), ob)
-jim(axes(sg), sino; title="Shepp-Logan 'short' sinogram",
+sino = radon(rays(rg), ob)
+jim(axes(rg), sino; title="Shepp-Logan 'short' sinogram",
     xlabel="r", ylabel="ϕ")
 
 
@@ -102,14 +102,14 @@ Here we start with a "plan",
 which would save work if we were reconstructing many images.
 =#
 
-plan = plan_fbp(sg, ig)
+plan = plan_fbp(rg, ig)
 
 # Examine Parker weights:
-jim(axes(sg), plan.parker_weight; title = "Parker weights",
+jim(axes(rg), plan.view_weight; title = "Parker weights",
     xlabel="r", ylabel="ϕ")
 
 # Finally perform FBP:
-fbp_image, sino_filt = fbp(plan, sino);
+fbp_image = fbp(plan, sino);
 
 # A narrow color window is needed to see the soft tissue structures:
 clim = (0.9, 1.1) .* μ

docs/lit/examples/07-fdk.jl

@@ -67,25 +67,25 @@ clim = (0.95, 1.05) .* μ
 # Here is the ideal phantom image:
 oversample = 3
 true_image = phantom(axes(ig)..., ob, oversample)
-jim(axes(ig), true_image, "True 3D Shepp-Logan phantom image"; clim)
+pt = jim(axes(ig), true_image, "True 3D Shepp-Logan phantom image"; clim)
 
 # Define the system geometry
 # (for some explanation use `?CtGeom`):
 p = (ns = 130, ds = 0.3cm, nt = 80, dt = 0.4cm, na = 50, dsd = 200cm, dod = 40cm)
-cg = CtFanArc( ; p...)
+rg = CtFanArc( ; p...)
 
 # Examine the geometry to verify the FOV
 # (this is more interesting when interacting via other Plot backends):
-ct_geom_plot3(cg, ig)
+ct_geom_plot3(rg, ig)
 
 #
 prompt()
 
 
 # CBCT projections
 # using `Sinogram.rays` and `ImagePhantoms.radon`:
-proj_arc = radon(rays(cg), ob)
-jim(cg.s, cg.t, proj_arc ;
+proj_arc = radon(rays(rg), ob)
+pa = jim(axes(rg)[1:2], proj_arc ;
     title="Shepp-Logan projections (arc)", xlabel="s", ylabel="t")
 
 
@@ -103,31 +103,32 @@ which would save work if we were reconstructing many images.
 For illustration we include `Hamming` window.
 =#
 
-plan = plan_fbp(cg, ig; window = Window(Hamming(), 1.0))
+plan = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))
 fdk_arc = fdk(plan, proj_arc)
-jim(axes(ig), fdk_arc, "FDK image (arc)"; clim)
+par = jim(axes(ig), fdk_arc, "FDK image (arc)"; clim)
 
 #
 err_arc = fdk_arc - true_image
-jim(axes(ig), err_arc, "Error image (arc)"; clim = (-1,1) .* (0.05μ))
+elim = (-1,1) .* (0.02μ)
+pae = jim(axes(ig), err_arc, "Error image (arc)"; clim = elim)
 
 
 #=
 ## Repeat with flat detector geometry
 =#
 
-cg = CtFanFlat( ; p...)
-proj_flat = radon(rays(cg), ob)
-jim(cg.s, cg.t, proj_flat ;
+rg = CtFanFlat( ; p...)
+proj_flat = radon(rays(rg), ob)
+pfp = jim(axes(rg)[1:2], proj_flat ;
     title="Shepp-Logan projections (flat)", xlabel="s", ylabel="t")
 
-plan = plan_fbp(cg, ig; window = Window(Hamming(), 1.0))
+plan = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))
 fdk_flat = fdk(plan, proj_flat)
-jim(axes(ig), fdk_flat, "FDK image (flat)"; clim)
+pfr = jim(axes(ig), fdk_flat, "FDK image (flat)"; clim)
 
 #
 err_flat = fdk_flat - true_image
-jim(axes(ig), err_flat, "Error image (flat)"; clim = (-1,1) .* (0.05μ))
+pfe = jim(axes(ig), err_flat, "Error image (flat)"; clim = elim)
 
 #=
 As expected for CBCT,
@@ -138,15 +139,21 @@ the largest errors are in the end slices.
 #=
 ## Short scan
 =#
-cg = CtFanFlat(:short ; p...)
-proj_short = radon(rays(cg), ob)
-jim(cg.s, cg.t, proj_short ;
+rg = CtFanFlat(:short ; p..., na = 200)
+proj_short = radon(rays(rg), ob)
+psp = jim(axes(rg)[1:2], proj_short ;
     title="Shepp-Logan projections (flat,short)", xlabel="s", ylabel="t")
 
-plan_short = plan_fbp(cg, ig; window = Window(Hamming(), 1.0))
+#
+plan_short = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))
+psw = jim(axes(rg)[[1,3]], plan_short.view_weight[:,1,:];
+    title = "View weights (including Parker)",
+    xlabel="s", ylabel="angle")
+
+#
 fdk_short = fdk(plan_short, proj_short)
-jim(axes(ig), fdk_short, "FDK image (flat,short)"; clim)
+psr = jim(axes(ig), fdk_short, "FDK image (flat,short)"; clim)
 
 #
 err_short = fdk_short - true_image
-jim(axes(ig), err_flat, "Error image (flat,short)"; clim = (-1,1) .* (0.05μ))
+pse = jim(axes(ig), err_flat, "Error image (flat,short)"; clim = elim)

docs/make.jl

@@ -24,7 +24,7 @@ binder_root_url =
 
 
 repo = eval(:($reps))
-DocMeta.setdocmeta!(repo, :DocTestSetup, :(using $reps); recursive=true)
+DocMeta.setdocmeta!(repo, :DocTestSetup, :(using $reps; import $reps); recursive=true)
 
 for (root, _, files) in walkdir(lit), file in files
     splitext(file)[2] == ".jl" || continue # process .jl files only