Improve type inference (#44)

* infer radon

* infer ellipse(Matrix)

* cannot infer ellipse(Matrix)

* infer phantom-over

* infer disk phantom

* infer shepp logan

* ellipse vector tuple

* fix usage

* infer xray1

* unify 2d tests

* unify 2d tests further

* unified shape2 tests

* purge old 2d tests

* radon gateway for 3d infer

* try to infer radon for ellipsoid

* toward 3d test unify

* radon_type

* type annotation

* infer tests

* clean 3d tests

* final tests

* rm todo

* rm todo

* v0.3.0

* type infer shepp_logan at last!

* fix M,N bug

* doc

* cover

* type infer

* infer!

* test/

Author: JeffFessler

Project.toml

@@ -1,7 +1,7 @@
 name = "ImagePhantoms"
 uuid = "71a99df6-f52c-4da1-bd2a-69d6f37f3252"
 authors = ["Jeff Fessler <[email protected]> and contributors"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 LazyGrids = "7031d0ef-c40d-4431-b2f8-61a8d2f650db"

README.md

@@ -59,8 +59,8 @@ function disk_phantom(title::String)
     x = (-M÷2:M÷2-1) * dx
     y = (-N÷2:N÷2-1) * dy
     params = disk_phantom_params( ; rhead = () -> rand(100:105))
-    objects = Ellipse(params) # vector of Ellipse objects
-    img = phantom(x, y, objects)
+    objects = ellipse(params) # vector of Object2d{Ellipse}
+    img = phantom(x, y, objects) # sampled at all (x,y) pairs
     jim(x, y, img; title, clim=(0,1300))
 end
 anim = @animate for i in 1:8

docs/lit/examples/02-ellipse.jl

@@ -127,7 +127,7 @@ dr = 0.2mm # radial sample spacing
 nr = 2^10 # radial sinogram bins
 r = (-nr÷2:nr÷2-1) * dr # radial samples
 fr = (-nr÷2:nr÷2-1) / nr / dr # corresponding spectral axis
-ϕ = deg2rad.(0:180) # * Unitful.rad # todo round unitful Unitful.°
+ϕ = deg2rad.(0:180)
 sino = radon(ob).(r, ϕ') # sample Radon transform of a single shape object
 smax = ob.value * 2 * maximum(radii)
 p5 = jim(r, rad2deg.(ϕ), sino; title="sinogram",

docs/lit/examples/07-shepp.jl

@@ -159,9 +159,9 @@ Sometimes it can be more convenient
 to have the middle ellipses be non-overlapping:
 =#
 
-ob = ellipse_parameters(SheppLoganBrainWeb(), disjoint=true)
-ob[:,end] = 1:10
-ob = ellipse(ob)
+params = ellipse_parameters(SheppLoganBrainWeb(), disjoint=true)
+params = [(p[1:5]..., i) for (i, p) in enumerate(params)]
+ob = ellipse(params)
 x = LinRange(-0.4, 0.4, 206)
 y = LinRange(-0.5, 0.5, 256)
 oversample = 3

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

@@ -88,10 +88,11 @@ with random complex phases
 to make it a bit more realistic.
 =#
 
-oa = ellipse_parameters(SheppLoganBrainWeb() ; disjoint=true, fovs)
+params = ellipse_parameters(SheppLoganBrainWeb() ; disjoint=true, fovs)
 seed!(0)
-oa[:,end] = [1; randn(ComplexF32, 9)] # random phases
-oa = ellipse(oa)
+phases = [1; rand(ComplexF32,9)] # random phases
+params = [(p[1:5]..., phases[i]) for (i, p) in enumerate(params)]
+oa = ellipse(params)
 oversample = 3
 image0 = phantom(x, y, oa, oversample)
 cfun = z -> cat(dims = ndims(z)+1, real(z), imag(z))

src/cuboid.jl

@@ -84,12 +84,12 @@ end
 # `u,v` should be unitless
 function xray1(
     ::Cuboid,
-    u::Real,
-    v::Real,
+    u::Ru,
+    v::Rv,
     ϕ::RealU, # azim
-    θ::RealU, # polar
-)
-    T = promote_type(eltype(u), eltype(v), Float32)
+    θ::RealU; # polar
+) where {Ru <: Real, Rv <: Real}
+    T = promote_type(Ru, Rv, Float32)
 
     (sϕ, cϕ) = sincos(ϕ)
     (sθ, cθ) = sincos(θ)
@@ -119,7 +119,7 @@ function xray1(
     if abs(e3) < eps(T)
         ℓ *= (-1/2 ≤ v < 1/2)
     end
-    return ℓ
+    return T(ℓ)::T
 end
 
 

src/disk-phantom.jl

@@ -12,7 +12,7 @@ using Random: seed!
 """
     params = disk_phantom_params( ; ...)
 
-Generate `ndisk × 6` ellipse phantom parameters
+Generate `ndisk` ellipse phantom parameter 6-tuples
 for a head-sized disk plus many disks within it,
 designed so that the disks have some minimum separation `minsep`
 to avoid overlap and to simplify patch-based model fitting.
@@ -31,7 +31,15 @@ to avoid overlap and to simplify patch-based model fitting.
 
 The function options can be replaced with rand() for other Distributions.
 """
-function disk_phantom_params( ;
+function disk_phantom_params( ; kwargs...)
+    params = _disk_phantom_params( ; kwargs...)
+
+    out = ellipse_parameters_uscale(params, (1,1), 1, 1, 1)
+    return out
+end
+
+
+function _disk_phantom_params( ;
     fov::RealU = 240,
 #   rhead::RealU = 100, # background radius for "head" [mm]
     rhead::Function = () -> 100, # background radius for "head" [mm]
@@ -63,9 +71,9 @@ function disk_phantom_params( ;
 
     (seed != 0) && seed!(seed)
 
-    for id=1:ndisk
+    for id in 1:ndisk
         trial = 0
-        for ii=1:maxtry
+        for ii in 1:maxtry
 #           rc = randu(0, rhead - rmin - minsep, f = x->x^0.5)
             rc = cdisk()
             phi = rand() * 2π

src/ellipse.jl

@@ -59,8 +59,8 @@ phantom1(ob::Object2d{Ellipse}, xy::NTuple{2,Real}) = (sum(abs2, xy) ≤ 1)
 
 # x-ray transform (line integral) of unit circle
 # `r` should be unitless
-function xray1(::Ellipse, r::Real, ϕ::RealU)
-    T = promote_type(eltype(r), Float32)
+function xray1(::Ellipse, r::R, ϕ::RealU) where {R <: Real}
+    T = promote_type(R, Float32)
     r2 = r^2
     return r2 < 1 ? 2 * sqrt(one(T) - r2) : zero(T)
 end

src/ellipsoid.jl

@@ -61,14 +61,14 @@ phantom1(ob::Object3d{Ellipsoid}, xyz::NTuple{3,Real}) = (sum(abs2, xyz) ≤ 1)
 # `u,v` should be unitless
 function xray1(
     ::Ellipsoid,
-    u::Real,
-    v::Real,
+    u::Ru,
+    v::Rv,
     ϕ::RealU, # irrelevant
     θ::RealU, # irrelevant
-)
-    T = promote_type(eltype(u), eltype(v), Float32)
+) where {Ru <: Real, Rv <: Real}
+    T = promote_type(Ru, Rv, Float32)
     r2 = u^2 + v^2
-    return r2 < 1 ? 2 * sqrt(one(T) - r2) : zero(T)
+    return (r2 < 1 ? 2 * sqrt(one(T) - r2) : zero(T))::T
 end
 
 

src/gauss2.jl

@@ -68,9 +68,10 @@ phantom1(ob::Object2d{Gauss2}, xy::NTuple{2,Real}) =
 
 # x-ray transform (line integral) of unit gauss2
 # `r` should be unitless
-function xray1(::Gauss2, r::Real, ϕ::RealU)
+function xray1(::Gauss2, r::R, ϕ::RealU) where {R <: Real}
+    T = promote_type(R, Float32)
     s = fwhm2spread(1)
-    return s * exp(-π * abs2(r / s))
+    return T(s * exp(-π * abs2(r / s)))
 end
 
 

src/mri-sense.jl

@@ -84,7 +84,7 @@ function mri_smap_fit(
 
     smaps_fit = [embed(B * coef, mask) for coef in coefs]
     ss(v) = sum(x -> sum(abs2, x), v)
-    nrmse = sqrt( ss(smaps_fit - smaps) / ss(smaps) )
+    nrmse = Float64( sqrt( ss(smaps_fit - smaps) / ss(smaps) ))::Float64
     return (; B, ν, coefs, nrmse, smaps = smaps_fit)
 end
 

src/object.jl

@@ -16,7 +16,8 @@ abstract type AbstractObject end
 
 
 """
-    Object{S, D, V, ...}(center, width, angle, value) <: AbstractObject
+    Object{S, D, V, C, A, Da}(center, width, angle, value) <: AbstractObject
+        where {S <: AbstractShape, V <: Number, C,A <: RealU}
 General container for 2D and 3D objects for defining image phantoms.
 
 * `center::NTuple{D,C}` coordinates of "center" of this object
@@ -83,16 +84,16 @@ end
 
 
 """
-    Object2d = Object{S,2} where S <: AbstractShape
+    Object2d{S,V,C} = Object{S,2,V,C} where {S <: AbstractShape, V,C <: Number}
 For 2D objects
 """
-const Object2d = Object{S,2} where S
+const Object2d{S,V,C} = Object{S,2,V,C}
 
 """
-    Object3d = Object{S,3} where S <: AbstractShape
+    Object3d{S,V,C} = Object{S,3,V,C} where {S <: AbstractShape, V,C <: Number}
 For 3D objects
 """
-const Object3d = Object{S,3} where S
+const Object3d{S,V,C} = Object{S,3,V,C}
 
 
 # constructors
@@ -229,3 +230,14 @@ translate(ob::Object{S,D}, shift::NTuple{D,RealU}) where {S, D} =
     Object(S(), ob.center .+ shift, ob.width, ob.angle, ob.value)
 translate(ob::Object{S,2}, x::RealU, y::RealU) where {S} = translate(ob, (x,y))
 translate(ob::Object{S,3}, x::RealU, y::RealU, z::RealU) where {S} = translate(ob, (x,y,z))
+
+
+"""
+    radon_type(::Object)
+Determine the element type of the Radon transform of an object
+(including units if applicable).
+Ensures that its precision is at least `Float32`.
+"""
+function radon_type(::Object{S, D, V, C, A}) where {S, D, V <: Number, C <: RealU, A <: RealU}
+    return eltype(oneunit(C) * oneunit(V) * one(A) * one(Float32))
+end

src/object2.jl

@@ -84,7 +84,7 @@ function phantom(
     x::AbstractVector,
     y::AbstractVector,
     oa::Array{<:Object2d},
-    oversample::Int;
+    oversample::Int ;
     T::DataType = promote_type(eltype.(oa)..., Float32),
 )
 
@@ -201,7 +201,7 @@ function radon(
     ϕ::AbstractVector,
     oa::Array{<:Object2d},
 )
-    return sum(ob -> radon(ob).(r,ϕ'), oa)
+    return radon(oa).(r,ϕ')
 end
 
 

src/object3.jl

@@ -88,10 +88,10 @@ function phantom(
     x::AbstractVector,
     y::AbstractVector,
     z::AbstractVector,
-    oa::Array{<:Object3d},
+    oa::Array{<:Object3d{S,V}},
     oversample::Int;
-    T::DataType = promote_type(eltype.(oa)..., Float32),
-)
+    T::DataType = promote_type(V, Float32),
+) where {S, V <: Number}
 
     oversample < 1 && throw(ArgumentError("oversample $oversample"))
     dx = x[2] - x[1]
@@ -104,7 +104,8 @@ function phantom(
     o3 = oversample^3
     ophantom = ob -> (x,y,z) ->
         T(sum(phantom(ob).(ndgrid(x .+ dx*tmp, y .+ dy*tmp, z .+ dz*tmp)...)) / o3)
-    return sum(ob -> ophantom(ob).(ndgrid(x,y,z)...), oa)
+    out = sum(ob -> ophantom(ob).(ndgrid(x,y,z)...), oa)
+    return out::Array{T, 3}
 end
 
 
@@ -173,6 +174,13 @@ function _xray(
 end
 
 
+# this gateway seems to help type inference
+function _radon(ob::Object3d{S}, u::RealU, v::RealU, ϕ::RealU, θ::RealU) where S
+    T = radon_type(ob)
+    return T(ob.value * _xray(S(), ob.center, ob.width, ob.angle, u, v, ϕ, θ))::T
+end
+
+
 """
     radon(ob::Object3d)::Function
 Return function of `(u,v,ϕ,θ)` that user can sample
@@ -185,8 +193,7 @@ for an object ``f(x,y,z)``.
 Then as `ϕ` increases, the line integrals rotate counter-clockwise.
 """
 function radon(ob::Object3d{S}) where S
-    return (u,v,ϕ,θ) -> ob.value *
-        _xray(S(), ob.center, ob.width, ob.angle, u, v, ϕ, θ)
+    return (u::RealU, v::RealU, ϕ::RealU, θ::RealU) -> _radon(ob, u, v, ϕ, θ)
 end
 
 
@@ -198,7 +205,7 @@ to make projection views
 for an array of 3D objects.
 """
 function radon(oa::Array{<:Object3d})
-    return (u,v,ϕ,θ) -> sum(ob -> radon(ob)(u,v,ϕ,θ), oa)
+    return (u::RealU, v::RealU, ϕ::RealU, θ::RealU) -> sum(ob -> radon(ob)(u,v,ϕ,θ), oa)
 end