Add cone (#62)

* Add cone

* v0.7

Project.toml

@@ -1,7 +1,7 @@
 name = "ImagePhantoms"
 uuid = "71a99df6-f52c-4da1-bd2a-69d6f37f3252"
 authors = ["Jeff Fessler <fessler@umich.edu> and contributors"]
-version = "0.6.0"
+version = "0.7.0"
 
 [deps]
 LazyGrids = "7031d0ef-c40d-4431-b2f8-61a8d2f650db"

docs/lit/examples/36-cone.jl

@@ -0,0 +1,93 @@
+#---------------------------------------------------------
+# # [Cone](@id 36-cone)
+#---------------------------------------------------------
+
+#=
+This page illustrates the `Cone` shape in the Julia package
+[`ImagePhantoms`](https://github.com/JuliaImageRecon/ImagePhantoms.jl).
+
+This page was generated from a single Julia file:
+[36-cone.jl](@__REPO_ROOT_URL__/36-cone.jl).
+=#
+
+#md # In any such Julia documentation,
+#md # you can access the source code
+#md # using the "Edit on GitHub" link in the top right.
+
+#md # The corresponding notebook can be viewed in
+#md # [nbviewer](http://nbviewer.jupyter.org/) here:
+#md # [`36-cone.ipynb`](@__NBVIEWER_ROOT_URL__/36-cone.ipynb),
+#md # and opened in [binder](https://mybinder.org/) here:
+#md # [`36-cone.ipynb`](@__BINDER_ROOT_URL__/36-cone.ipynb).
+
+
+# ### Setup
+
+# Packages needed here.
+
+using ImagePhantoms: Object, phantom, radon, spectrum
+using ImagePhantoms: Cone, cone
+import ImagePhantoms as IP
+using ImageGeoms: ImageGeom, axesf
+using MIRTjim: jim, prompt, mid3
+using FFTW: fft, fftshift, ifftshift
+using LazyGrids: ndgrid
+using Unitful: mm, unit, °
+using Plots: plot, plot!, scatter!, default
+using Plots # gif @animate
+default(markerstrokecolor=:auto)
+
+
+# The following line is helpful when running this file as a script;
+# this way it will prompt user to hit a key after each figure is displayed.
+
+isinteractive() ? jim(:prompt, true) : prompt(:draw);
+
+
+# ### Overview
+
+#=
+A basic shape used in constructing 3D digital image phantoms
+is the cone,
+specified by its center, base radii, height, angle(s) and value.
+All of the methods in `ImagePhantoms` support physical units,
+so we use such units throughout this example.
+(Using units is recommended but not required.)
+
+Here are 3 ways to define a `Object{Cone}`,
+using physical units.
+=#
+
+center = (20mm, 10mm, -15mm)
+width = (25mm, 30mm, 35mm) # x radius, y radius, height
+ϕ0s = :(π/6) # symbol version for nice plot titles
+ϕ0 = eval(ϕ0s)
+angles = (ϕ0, 0, 0)
+Object(Cone(), center, width, angles, 1.0f0) # top-level constructor
+cone( 20mm, 10mm, 5mm, 25mm, 35mm, 15mm, π/6, 0, 0, 1.0f0) # 9 arguments
+ob = cone(center, width, angles, 1.0f0) # tuples (recommended use)
+
+
+#=
+### Phantom image using `phantom`
+
+Make a 3D digital image of it using `phantom` and display it.
+We use `ImageGeoms` to simplify the indexing.
+=#
+
+deltas = (1.0mm, 1.1mm, 0.9mm)
+dims = (2^8, 2^8+2, 49) # odd
+ig = ImageGeom( ; dims, deltas, offsets=:dsp)
+oversample = 3
+img = phantom(axes(ig)..., [ob], oversample)
+p1 = jim(axes(ig), img;
+   title="Cone, rotation ϕ=$ϕ0s", xlabel="x", ylabel="y")
+
+
+# The image integral should match the object volume:
+volume = IP.volume(ob)
+(sum(img)*prod(ig.deltas), volume)
+
+
+# Show middle slices
+jim(mid3(img), "Middle 3 planes")

src/ImagePhantoms.jl

@@ -18,11 +18,12 @@ include("rect.jl")
 include("triangle.jl")
 
 # 3d:
+include("cone.jl")
+include("cuboid.jl")
+include("cylinder.jl")
 include("dirac3.jl")
 include("ellipsoid.jl")
 include("gauss3.jl")
-include("cuboid.jl")
-include("cylinder.jl")
 
 # phantoms:
 include("shepplogan.jl")

src/cone.jl

@@ -0,0 +1,89 @@
+#=
+cone.jl
+Right cone
+Base object is a cone with unit radius and unit height,
+with base at origin pointing up along z.
+=#
+
+
+export Cone, cone
+
+
+"""
+    Cone <: AbstractShape{3}
+"""
+struct Cone <: AbstractShape{3} end
+
+
+# constructor
+
+
+"""
+    cone(cx, cy, cz, wx, wy, wz, Φ, Θ, value::Number)
+    cone(center::NTuple{3,RealU}, width::NTuple{3,RealU}, angle::NTuple{3,RealU}, v)
+Construct `Object{Cone}` from parameters;
+here `width` is the base radii for `wx` and `wy` and `wz` is the height.
+"""
+cone(args... ; kwargs...) = Object(Cone(), args...; kwargs...)
+
+
+# methods
+
+
+volume1(::Cone) = π/3 # volume of unit cone
+
+ℓmax1(::Cone) = 2 # max line integral through unit cone
+
+function ℓmax(ob::Object3d{Cone})
+    rmax = maximum(ob.width[1:2])
+    return max(2 * rmax, sqrt(rmax^2 + ob.width[3]^2))
+end
+
+
+"""
+    phantom1(ob::Object3d{Cone}, (x,y,z))
+Evaluate unit cone at `(x,y,z)`,
+for unitless coordinates.
+"""
+function phantom1(ob::Object3d{Cone}, xyz::NTuple{3,Real})
+    z = xyz[3]
+    r = sqrt(sum(abs2, xyz[1:2]))
+    return (0 ≤ z ≤ 1) && r ≤ 1 - z
+end
+
+
+# radon
+
+
+# x-ray transform (line integral) of unit cone
+# `u,v` should be unitless
+#=
+function xray1(
+    ::Cone,
+    u::Ru,
+    v::Rv,
+    ϕ::RealU, # azim
+    θ::RealU; # polar
+) where {Ru <: Real, Rv <: Real}
+    T = promote_type(Ru, Rv, Float32)
+
+#   (sϕ, cϕ) = sincos(ϕ)
+#   (sθ, cθ) = sincos(θ)
+
+    throw("todo: cone radon not done")
+end
+=#
+
+
+# spectrum
+
+#=
+"""
+    spectrum1(::Object3d{Cone}, (kx,ky,kz))
+Spectrum of unit cone at `(kx,ky,kz)`,
+for unitless spatial frequency coordinates.
+"""
+function spectrum1(::Object3d{Cone}, kxyz::NTuple{3,Real})
+    throw("todo: cone spectrum not done")
+end
+=#

test/shape3.jl

@@ -4,6 +4,7 @@ using ImagePhantoms: Object3d, sphere, cube
 using ImagePhantoms: Object, AbstractShape, phantom, radon, spectrum
 using ImagePhantoms: Gauss3, gauss3
 using ImagePhantoms: Ellipsoid, ellipsoid
+using ImagePhantoms: Cone, cone
 using ImagePhantoms: Cuboid, cuboid
 using ImagePhantoms: Cylinder, cylinder
 using ImagePhantoms: Dirac3, dirac3
@@ -43,12 +44,14 @@ end
 
 # parameters for testing each shape
 # (Shape, shape, lmax, lmax1, tol1, tolk, tolp)
+# for width = (4//1, 5, 6)
 list = [
  (Dirac3, dirac3, 6, 1, Inf, Inf, Inf)
  (Ellipsoid, ellipsoid, 12, 2, 1e-2, 4e-4, 1e-3)
  (Gauss3, gauss3, IP.fwhm2spread(6), IP.fwhm2spread(1), 1e-2, 5e-4, 1e-5)
  (Cuboid, cuboid, sqrt(4^2 + 5^2 + 6^2), √3, 1e-2, 2e-2, 2e-3)
  (Cylinder, cylinder, sqrt(10^2 + 6^2), √5, 2e-2, 2e-2, 1e-3)
+ (Cone, cone, 10, 2, Inf, Inf, Inf)
 ]
 
 macro isob3(ex) # macro to streamline tests
@@ -109,12 +112,12 @@ for (Shape, shape, lmax, lmax1, tol1, tolk, tolp) in list
 
     has_xray = hasmethod(IP.xray1, (Shape, Real, Real, Real, Real))
     has_phantom = hasmethod(IP.phantom1, (typeof(shape()), NTuple{3,Real}))
+    has_spectrum = hasmethod(IP.spectrum1, (typeof(shape()), NTuple{3,Real}))
 
     @testset "construct-$shape" begin
         test3_construct(Shape, shape)
     end
 
-
     @testset "operations" begin
         test3_op(Shape, shape)
     end
@@ -161,8 +164,10 @@ for (Shape, shape, lmax, lmax1, tol1, tolk, tolp) in list
 
     volume = IP.volume(ob)
 
-    fun = @inferred spectrum([ob])
-    @test (@inferred fun(0,0,0)) ≈ ob.value * volume
+    if has_spectrum
+        fun = @inferred spectrum([ob])
+        @test (@inferred fun(0,0,0)) ≈ ob.value * volume
+    end
 
     if has_xray
         test3_radon(Shape, shape, ob)
@@ -190,10 +195,12 @@ for (Shape, shape, lmax, lmax1, tol1, tolk, tolp) in list
     volume = IP.volume(ob)
     zscale = 1 / (ob.value * volume) # normalize spectra
 
-    kspace = (@inferred spectrum(axesf(ig)..., [ob])) * zscale
-    @test spectrum(ob)(0/m, 0/m, 0/m) * zscale ≈ 1
-    @test maximum(abs, kspace) ≈ 1
-    @test kspace[L÷2+1,M÷2+1,N÷2+1] ≈ 1
+    if has_spectrum
+        kspace = (@inferred spectrum(axesf(ig)..., [ob])) * zscale
+        @test spectrum(ob)(0/m, 0/m, 0/m) * zscale ≈ 1
+        @test maximum(abs, kspace) ≈ 1
+        @test kspace[L÷2+1,M÷2+1,N÷2+1] ≈ 1
+    end
 
     if has_phantom
         oversample = 2
@@ -202,8 +209,11 @@ for (Shape, shape, lmax, lmax1, tol1, tolk, tolp) in list
         X = myfft(img) * (prod(ig.deltas) * zscale)
 
         @test abs(maximum(abs, X) - 1) < tol1
-        errk = maximum(abs, kspace - X) / maximum(abs, kspace)
-        @test errk < tolk
+
+        if has_spectrum
+            errk = maximum(abs, kspace - X) / maximum(abs, kspace)
+            @test errk < tolk
+        end
     end
   end