Add Dirac (#51)

* Add Dirac

* add dirac3

* refine dirac2

* fix export

* test dirc

* hide xray1

* hide phantom1

* tests without xray1 phantom1

* v0.4

Project.toml

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

src/ImagePhantoms.jl

@@ -9,13 +9,15 @@ include("object2.jl")
 include("object3.jl")
 include("jinc.jl")
 
-# shapes:
+# 2d:
+include("dirac2.jl")
 include("ellipse.jl")
 include("gauss2.jl")
 include("rect.jl")
 include("triangle.jl")
 
 # 3d:
+include("dirac3.jl")
 include("ellipsoid.jl")
 include("gauss3.jl")
 include("cuboid.jl")

src/cuboid.jl

@@ -4,10 +4,7 @@ actually a rectangular cuboid https://en.wikipedia.org/wiki/Cuboid
 =#
 
 
-#using ImagePhantoms #: Object, Object3d
-
 export Cuboid, cuboid, cube
-export phantom, radon, spectrum
 
 
 """

src/cylinder.jl

@@ -4,7 +4,6 @@ cylinder.jl
 
 
 export Cylinder, cylinder
-export phantom, radon, spectrum
 
 
 """

src/dirac2.jl

@@ -0,0 +1,67 @@
+#=
+dirac2.jl
+=#
+
+export Dirac2, dirac2
+#export Dirac
+
+
+"""
+    Dirac{D} <: AbstractShape{D}
+"""
+struct Dirac{D} <: AbstractShape{D} end
+
+const Dirac2 = Dirac{2}
+
+
+# constructor
+
+
+"""
+    dirac2(...)
+Construct `Object2d{Dirac2}` from parameters;
+here `width` is a scale factor: `δ((x - c)/w) = w * δ(x - c)`.
+"""
+dirac2(args... ; kwargs...) = Object(Dirac2(), args...; kwargs...)
+# dirac(D::Int, args... ; kwargs...) = Object(Dirac{D}(), args...; kwargs...)
+
+
+# methods
+
+
+area1(::Dirac2) = 1
+
+ℓmax1(::Dirac) = 1 # max line integral
+
+
+#=
+"""
+    phantom1(ob::Object2d{Dirac2}, (x,y))
+Evaluate Dirac at `(x,y)`,
+for unitless coordinates.
+Caution: this method is meaningless
+because the Dirac impulse is not a function!
+"""
+phantom1(ob::Object2d{Dirac2}, xy::NTuple{2,Real}) = (xy == (0,0)) ? 1 : 0
+=#
+
+
+#=
+"""
+    xray1(::Dirac2, r::Real, ϕ::RealU)
+X-ray transform (line integral) of Dirac.
+`r` should be unitless
+Caution: this method is meaningless
+because the Dirac impulse is not a function!
+"""
+xray1(::Dirac2, r::Real, ϕ::RealU) = (r == 0) ? 1 : 0
+=#
+
+
+"""
+    spectrum1(::Object2d{Dirac2}, (kx,ky))
+Spectrum of Dirac2 at `(kx,ky)`,
+for unitless spatial frequency coordinates.
+"""
+spectrum1(::Object2d{Dirac2}, kxy::NTuple{2,Real}) = 1
+#spectrum1(::Object{Dirac, D, V}, kxy::NTuple{D,Real}) where {D, V} = one(V)

src/dirac3.jl

@@ -0,0 +1,58 @@
+#=
+dirac3.jl
+=#
+
+export Dirac3, dirac3
+
+const Dirac3 = Dirac{3}
+
+
+# constructor
+
+
+"""
+    dirac3(...)
+Construct `Object3d{Dirac3}` from parameters;
+here `width` is a scale factor: `δ((x - c)/w) = w * δ(x - c)`.
+"""
+dirac3(args... ; kwargs...) = Object(Dirac3(), args...; kwargs...)
+
+
+# methods
+
+
+volume1(::Dirac3) = 1
+
+# ℓmax1(::Dirac) = 1 # max line integral
+
+
+#=
+"""
+    phantom1(ob::Object3d{Dirac3}, (x,y))
+Evaluate Dirac at `(x,y,z)`,
+for unitless coordinates.
+Caution: this method is meaningless
+because the Dirac impulse is not a function!
+"""
+phantom1(ob::Object3d{Dirac3}, xyz::NTuple{3,Real}) = (xyz == (0,0,0)) ? 1 : 0
+=#
+
+
+#=
+"""
+    xray1(::Dirac3, r::Real, ϕ::RealU)
+X-ray transform (line integral) of Dirac.
+`r` should be unitless
+Caution: this method is meaningless
+because the Dirac impulse is not a function!
+"""
+xray1(::Dirac3, u::Real, v::Real, ϕ::RealU, θ::RealU) = ((u,v) == (0,0)) ? 1 : 0
+=#
+
+
+"""
+    spectrum1(::Object3d{Dirac3}, (kx,ky,kz))
+Spectrum of Dirac3 at `(kx,ky,kz)`,
+for unitless spatial frequency coordinates.
+"""
+spectrum1(::Object3d{Dirac3}, kxyz::NTuple{3,Real}) = 1

src/ellipse.jl

@@ -2,11 +2,7 @@
 ellipse.jl
 =#
 
-
-#using ImagePhantoms #: Object, Object2d
-
 export Ellipse, ellipse, circle
-export phantom, radon, spectrum
 
 
 """

src/ellipsoid.jl

@@ -2,11 +2,7 @@
 ellipsoid.jl
 =#
 
-
-#using ImagePhantoms #: Object, Object3d
-
 export Ellipsoid, ellipsoid, sphere
-export phantom, radon, spectrum
 
 
 """

src/gauss2.jl

@@ -2,11 +2,7 @@
 gauss2.jl
 =#
 
-
-#using ImagePhantoms #: Object, Object2d
-
 export Gauss2, gauss2
-export phantom, radon, spectrum, fwhm2spread
 
 
 """

src/gauss3.jl

@@ -2,11 +2,7 @@
 gauss3.jl
 =#
 
-
-#using ImagePhantoms #: Object, Object3d
-
 export Gauss3, gauss3
-export phantom, radon, spectrum
 
 
 """

src/rect.jl

@@ -3,10 +3,7 @@ rect.jl
 =#
 
 
-#using ImagePhantoms #: Object, Object2d
-
 export Rect, rect, square
-export phantom, radon, spectrum
 
 
 """

src/triangle.jl

@@ -3,11 +3,7 @@ triangle.jl
 2D triangle shape
 =#
 
-
-#using ImagePhantoms #: Object, Object2d, AbstractShape
-
 export Triangle, triangle
-export phantom, radon, spectrum
 
 const sqrt3 = sqrt(3)