Parker weighting for CBCT (#46)

* Test f32

* Construct :short cases

* Move parker shape to plan

* 2D parker weights

* Parker weighting for 3d

* Test Parker weighting

* Add short scan example

Author: JeffFessler

docs/lit/examples/07-fdk.jl

@@ -133,3 +133,20 @@ jim(axes(ig), err_flat, "Error image (flat)"; clim = (-1,1) .* (0.05μ))
 As expected for CBCT,
 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 ;
+    title="Shepp-Logan projections (flat,short)", xlabel="s", ylabel="t")
+
+plan_short = plan_fbp(cg, ig; window = Window(Hamming(), 1.0))
+fdk_short = fdk(plan_short, proj_short)
+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μ))

src/fbp2/fbp.jl

@@ -3,7 +3,7 @@
 export fbp
 
 using ImageGeoms: ImageGeom
-# using Sinograms: SinoGeom, SinoPar, SinoFan, SinoMoj, parker_weight, _reale, FBPplan
+# using Sinograms: SinoGeom, SinoPar, SinoFan, SinoMoj, _reale, FBPplan
 
 
 """
@@ -21,6 +21,8 @@ out
 - `sino_filt::Matrix{<:Number}`   filtered sinogram(s)
 
 """
+fbp
+
 function fbp(
     plan::FBPNormalPlan{<:SinoPar},
     sino::AbstractMatrix{<:Number},
@@ -62,14 +64,11 @@ function fbp(
     ig::ImageGeom,
     sino::AbstractMatrix{<:Number},
     filter::AbstractVector{<:Number},
-    parker_weight::AbstractArray{<:Number} = ones(size(sino,2)),
+    parker_weight::AbstractMatrix{<:Number} = ones(1,1),
 )
     dims(sg) == size(sino) || error("bad sino size")
-    parker = length(parker_weight) == sg.na ? parker_weight' :
-        size(parker_weight) == dims(sg) ? parker_weight :
-        error("bad parker size $(size(parker_weight))")
 
-    sino_filt = sino .* parker
+    sino_filt = sino .* parker_weight
     sino_filt = fbp_sino_filter(sino_filt, filter)
 
     image = fbp_back(sg, ig, sino_filt)

src/fbp2/parker.jl

@@ -6,28 +6,41 @@ export parker_weight
 using LazyGrids: ndgrid
 
 
-# parallel-beam case: returns a Vector
+"""
+    parker_weight(sg::SinoGeom ; T = Float32)
+Compute Parker weighting for non-360° orbits.
+See http://doi.org/10.1118/1.595078.
+Returns `Matrix{T}` of size:
+- (1,1) for `SinoPar` with typical 180 or 360 orbit
+- (1,na) for `SinoPar` with atypical orbit
+- (1,1) for `SinoFan` with typical 360 orbit
+- (ns,na) for `SinoFan` with typical 360 orbit
+"""
+parker_weight
+
+
+# parallel-beam case
 function parker_weight_par(
     orbit::RealU,
     ad::AbstractVector{<:RealU}, # angles in degrees: sg.ad - sg.orbit_start
     ;
     T::Type{<:Real} = Float32,
 )
 
-    wt = ones(T, length(ad))
-
     if (orbit ÷ 180) * 180 == orbit
-        return wt # no weighting needed
+        return ones(T, 1, 1) # no weighting needed
     end
 
-    if orbit < 180
+    if abs(orbit) < 180
         @warn("orbit $orbit < 180")
-        return wt # nonuniform weighting would not help
+        return ones(T, 1, 1) # nonuniform weighting would not help
     end
 
+    orbit = abs(orbit)
     orbit > 360 && error("only 180 ≤ |orbit| ≤ 360 supported for Parker weighting")
     extra = orbit - 180 # extra beyond 180
 
+    wt = ones(T, 1, length(ad)) # (1,na)
     ad = abs.(ad)
     ii = ad .< extra
     @. wt[ii] = abs2(sin(ad[ii] / extra * π/2))
@@ -38,17 +51,11 @@ function parker_weight_par(
 end
 
 
-"""
-    parker_weight(sg::SinoGeom ; T = Float32)
-Compute Parker weighting for non-360° orbits.
-Returns `Vector{T}` of length `sg.na` for `SinoPar`.
-Returns `Matrix{T}` of size `dims(sg)` for `SinoFan`.
-"""
-parker_weight(sg::SinoPar ; T::Type{<:Real} = Float32)::Vector{T} =
+parker_weight(sg::SinoPar ; T::Type{<:Real} = Float32)::Matrix{T} =
     parker_weight_par(sg.orbit, sg.ad .- sg.orbit_start ; T)
 
 
-function parker_weight_fan(
+function parker_weight_fan_short(
     nb::Int,
     na::Int,
     orbit::RealU,
@@ -65,7 +72,7 @@ function parker_weight_fan(
         @warn("orbit $orbit is less than a short scan $orbit_short")
 
     orbit > orbit_short + rad2deg(ar[2] - ar[1]) &&
-        @warn("orbit $orbit exeeds short scan $orbit_short by %g views")
+        @warn("orbit $orbit exceeds short scan $orbit_short by %g views")
     #   (orbit - orbit_short) / rad2deg(ar[2] - ar[1]))
 
     bet = ar .- ar[1] # trick: force 0 start, so this ignores orbit_start!
@@ -96,22 +103,47 @@ end
 
 
 function parker_weight(sg::SinoFan; T::Type{<:Real} = Float32)::Matrix{T}
-    wt = ones(T, dims(sg))
     if (sg.orbit ÷ 360) * 360 == sg.orbit
-        return wt
+        return ones(T, 1, 1) # no weighting needed
     end
-    return parker_weight_fan(
+    return parker_weight_fan_short(
         sg.nb, sg.na, sg.orbit, sg.orbit_short,
         sg.ar, sg.gamma, sg.gamma_max; T
     )
 end
 
 
-function parker_weight(sg::SinoMoj ; T::Type{<:Real} = Float32)::Vector{T}
+function parker_weight(sg::SinoMoj ; T::Type{<:Real} = Float32)
     orbit = abs(sg.orbit)
     na = sg.na
     ((sg.orbit ÷ 180) * 180 == orbit) ||
         throw("No Parker weighting for Mojette geometry with orbit=$orbit")
-    wt = ones(T, na)
+    wt = ones(T, 1, 1)
     return wt
 end
+
+
+function parker_weight_fan_short(cg::CtFan; kwargs...)
+    weight = parker_weight_fan_short(
+        cg.ns, cg.na, cg.orbit, cg.orbit_short,
+        cg.ar, cg.gamma, cg.gamma_max; kwargs...,
+    )
+    weight .*= 360 / cg.orbit_short # trick due to scaling in cbct-back
+    return weight
+end
+
+
+"""
+    parker_weight(cg::CtFan; T::Type{<:Real} = Float32, kwargs...)
+For 3D case, return `Array{T,3}` where size is
+- `(1,1,1)` typical fan case with 360° orbit
+- `(ns,1,na)` atypical fan case including short scan
+"""
+function parker_weight(cg::CtFan; T::Type{<:Real} = Float32, kwargs...)
+    if (cg.orbit ÷ 360) * 360 == cg.orbit
+        return ones(T, 1, 1, 1) # (1,1,1) for type stability
+    end
+    weight = parker_weight_fan_short(cg; kwargs...)
+    weight = reshape(weight, dims(cg)[1], 1, dims(cg)[3]) # (ns,1,na)
+    return weight
+end

src/fbp2/plan2.jl

@@ -35,7 +35,7 @@ struct FBPNormalPlan{
     sg::S
     ig::I
     filter::H # frequency response Hk of apodized ramp filter, length npad
-    parker_weight::P # typically a 1D or 2D array of nonnegative reals
+    parker_weight::P # typically a Matrix of nonnegative reals
 #   moj::Moj # todo
 
 #=

src/fbp3/cbct-back.jl

@@ -35,7 +35,7 @@ function cbct_back(
         (oneunit(Td) * oneunit(To) / oneunit(Td) + oneunit(Toffset)))
 
     return cbct_back_fan(proj,
-        cg.ar,
+        cg.ar, # "betas"
         cg.dsd, cg.dso,
 #       cg.offset_source::RealU,
         cg.ds, cg.dt,
@@ -129,7 +129,7 @@ function cbct_back_fan!(
     cosβ = cos.(betas)
 
     # scale projections by dβ for Riemann-like integration
-    dβ = diff(betas)
+    dβ = diff(betas) # typically 2π/na for 360° orbit
     dβ = [dβ[1]; dβ] # todo
     proj = proj .* reshape(dβ, 1, 1, :) / 2
 

src/fbp3/fdk.jl

@@ -14,8 +14,10 @@ function fdk_weight_cyl_arc(
     t::RealU,
     dsd::RealU,
     dso::RealU,
+#   T::Type{<:AbstractFloat} = Float32,
 )
-    return (dso/dsd) * cos(s / dsd) / sqrt(abs2(t / dsd) + 1)
+    T = Float32
+    return T((dso/dsd) * cos(s / dsd) / sqrt(abs2(t / dsd) + 1))
 end
 
 function fdk_weight_cyl_flat(
@@ -24,15 +26,18 @@ function fdk_weight_cyl_flat(
     dsd::RealU,
     dso::RealU,
 )
-    return dso / sqrt(abs2(s) + abs2(t) + abs2(dsd))
+    T = Float32
+    return T(dso / sqrt(abs2(s) + abs2(t) + abs2(dsd)))
 end
 
 # (ns,nt) matrix
 fdk_weight_cyl(cg::CtFanArc) =
-    fdk_weight_cyl_arc.(cg.s, cg.t', cg.dsd, cg.dso)
+    fdk_weight_cyl_arc.(cg.s, cg.t', cg.dsd, cg.dso)::Matrix{Float32}
+#   Iterators.map((st) -> fdk_weight_cyl_arc(st..., cg.dsd, cg.dso),
+#       Iterators.product(ct_geom_s(cg), ct_geom_t(cg))) # eltype = Any !?
 
 fdk_weight_cyl(cg::CtFanFlat) =
-    fdk_weight_cyl_flat.(cg.s, cg.t', cg.dsd, cg.dso)
+    fdk_weight_cyl_flat.(cg.s, cg.t', cg.dsd, cg.dso)::Matrix{Float32}
 
 
 """
@@ -55,8 +60,12 @@ References: Feldkamp, Davis, Kress, JOSA-A, 1(6):612-9, June 1984.
 function fdk(plan::FDKplan, proj::AbstractArray{<:Number,3})
     cg = plan.cg
 
-    # step 1: apply weights
+    size(proj) == dims(plan.cg) ||
+        error("size mismatch $(size(proj)) $(dims(plan.cg))")
+
+    # step 1: apply cone-beam weights
     proj = proj .* fdk_weight_cyl(cg) # todo: precompute with plan!
+    proj .*= plan.parker_weight # todo: before or after filtering?
 
     # step 2: filter each projection view
     proj = fbp_sino_filter(proj, plan.filter)

src/fbp3/plan3.jl

@@ -35,14 +35,6 @@ struct FDKplan{
 end
 
 
-function parker_weight(cg::CtFan; kwargs...)
-    return parker_weight_fan(
-        cg.ns, cg.na, cg.orbit, cg.orbit_short,
-        cg.ar, cg.gamma, cg.gamma_max; kwargs...,
-    )
-end
-
-
 """
     plan = plan_fbp(cg, ig; window=Window(), ...)
 

src/geom/common.jl

@@ -12,7 +12,7 @@ export angles
 Base.show(io::IO, ::MIME"text/plain", st::RayGeom) = _show(io,  MIME("text/plain"), st)
 
 
-function _getproperty(st::RayGeom, s::Symbol, arg)
+function _getproperty(st::RayGeom, s::Symbol, arg) # not type stable :(
     d = Dict(arg)
     return haskey(d, s) ? d[s](st) : getfield(st, s)
 end