Skip to content

Commit bbac39c

Browse files
JeffFesslerJeffFessler
authored
Add area method and clean tests using deg2rad(1) (#42)
* area1 * area scale * |Diff| * volume1 * deg2rad * cover * clean todo
1 parent ecffa88 commit bbac39c

27 files changed

+65
-59
lines changed

‎docs/lit/examples/02-ellipse.jl

+1-1
Original file line numberDiff line numberDiff line change
@@ -114,7 +114,7 @@ p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
114114
err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
115115
@assert err < 2e-2
116116
p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact),
117-
"Difference × 10³"; xlabel, ylabel)
117+
"|Difference| × 10³"; xlabel, ylabel)
118118
jim(p1, p4, p2, p3)
119119

120120

‎docs/lit/examples/03-rect.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -112,7 +112,7 @@ p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
112112
err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
113113
@assert err < 2e-2
114114
p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact),
115-
"Difference × 10³"; xlabel, ylabel)
115+
"|Difference| × 10³"; xlabel, ylabel)
116116
jim(p1, p4, p2, p3)
117117

118118

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

‎docs/lit/examples/04-gauss.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -114,7 +114,7 @@ p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
114114
err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
115115
@assert err < 4e-4
116116
p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact),
117-
"Difference × 10³"; xlabel, ylabel)
117+
"|Difference| × 10³"; xlabel, ylabel)
118118
jim(p1, p4, p2, p3)
119119

120120

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

‎docs/lit/examples/05-triangle.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -117,7 +117,7 @@ p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
117117
err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
118118
@assert err < 2e-2
119119
p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact),
120-
"Difference × 10³"; xlabel, ylabel)
120+
"|Difference| × 10³"; xlabel, ylabel)
121121
jim(p1, p4, p2, p3)
122122

123123

@@ -130,7 +130,7 @@ dr = 0.2mm # radial sample spacing
130130
nr = 2^10 # radial sinogram bins
131131
r = (-nr÷2:nr÷2-1) * dr # radial samples
132132
fr = (-nr÷2:nr÷2-1) / nr / dr # corresponding spectral axis
133-
ϕ = deg2rad.(0:180) # * Unitful.rad # todo round unitful Unitful.°
133+
ϕ = deg2rad.(0:180)
134134
sino = radon(ob).(r, ϕ') # sample Radon transform of a single shape object
135135
smax = ob.value * sqrt((width[1]/2)^2 + (width[2] * sqrt(3) / 2)^2)
136136
p5 = jim(r, rad2deg.(ϕ), sino; title="sinogram",

‎docs/lit/examples/32-ellipsoid.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -167,7 +167,7 @@ vols = round.(((p -> sum(p)*prod(pg.deltas)).(proj2)..., volume) ./ 1mm^3; digi
167167

168168
# Look at a set of projections as the views orbit around the object.
169169
ϕd = 0:6:360
170-
ϕs = deg2rad.(ϕd) # * Unitful.rad # todo round unitful Unitful.°
170+
ϕs = deg2rad.(ϕd)
171171
θs = :(π/7)
172172
θ = eval(θs)
173173
projs = radon(axes(pg)..., ϕs, [θ], [ob]) # many projection views
@@ -221,7 +221,7 @@ p8 = jim(axesf(pg), sp.(proj_fft); prompt=false,
221221
err = maximum(abs, spectrum_slice - proj_fft) / maximum(abs, spectrum_slice)
222222
@assert err < 1e-3
223223
p9 = jim(axesf(pg), 1e3*abs.(proj_fft - spectrum_slice);
224-
title="Difference × 10³", xlabel, ylabel, prompt=false)
224+
title="|Difference| × 10³", xlabel, ylabel, prompt=false)
225225
jim(p6, p7, p8, p9)
226226

227227

‎docs/lit/examples/33-cuboid.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -170,7 +170,7 @@ vols = round.(((p -> sum(p)*prod(pg.deltas)).(proj3)..., volume) ./ 1mm^3; digit
170170

171171
# Look at a set of projections as the views orbit around the object.
172172
ϕd = 0:6:360
173-
ϕs = deg2rad.(ϕd) # * Unitful.rad # todo round unitful Unitful.°
173+
ϕs = deg2rad.(ϕd)
174174
θs = :(π/7)
175175
θ = eval(θs)
176176
projs = radon(axes(pg)..., ϕs, [θ], [ob]) # many projection views
@@ -224,7 +224,7 @@ p8 = jim(axesf(pg), sp.(proj_fft); prompt=false,
224224
err = maximum(abs, spectrum_slice - proj_fft) / maximum(abs, spectrum_slice)
225225
@assert err < 2e-3
226226
p9 = jim(axesf(pg), 1e6*abs.(proj_fft - spectrum_slice);
227-
title="Difference × 10⁶", xlabel, ylabel, prompt=false)
227+
title="|Difference| × 10⁶", xlabel, ylabel, prompt=false)
228228
jim(p6, p7, p8, p9)
229229

230230

‎docs/lit/examples/34-gauss3.jl

+3-3
Original file line numberDiff line numberDiff line change
@@ -128,7 +128,7 @@ p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
128128
err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
129129
@assert err < 1e-3
130130
p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact);
131-
title="Difference × 10³", xlabel, ylabel)
131+
title="|Difference| × 10³", xlabel, ylabel)
132132
jim(p1, p4, p2, p3)
133133

134134

@@ -167,7 +167,7 @@ vols = round.(((p -> sum(p)*prod(pg.deltas)).(proj2)..., volume) ./ 1mm^3; digit
167167

168168
# Look at a set of projections as the views orbit around the object.
169169
ϕd = 0:6:360
170-
ϕs = deg2rad.(ϕd) # * Unitful.rad # todo round unitful Unitful.°
170+
ϕs = deg2rad.(ϕd)
171171
θs = :(π/7)
172172
θ = eval(θs)
173173
projs = radon(axes(pg)..., ϕs, [θ], [ob]) # many projection views
@@ -221,7 +221,7 @@ p8 = jim(axesf(pg), sp.(proj_fft); prompt=false,
221221
err = maximum(abs, spectrum_slice - proj_fft) / maximum(abs, spectrum_slice)
222222
@assert err < 1e-5
223223
p9 = jim(axesf(pg), 1e6*abs.(proj_fft - spectrum_slice);
224-
title="Difference × 10⁶", xlabel, ylabel, prompt=false)
224+
title="|Difference| × 10⁶", xlabel, ylabel, prompt=false)
225225
jim(p6, p7, p8, p9)
226226

227227

‎docs/lit/examples/35-cylinder.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -170,7 +170,7 @@ vols = round.(((p -> sum(p)*prod(pg.deltas)).(proj3)..., volume) ./ 1mm^3; digit
170170

171171
# Look at a set of projections as the views orbit around the object.
172172
ϕd = 0:6:360
173-
ϕs = deg2rad.(ϕd) # * Unitful.rad # todo round unitful Unitful.°
173+
ϕs = deg2rad.(ϕd)
174174
θs = :(π/7)
175175
θ = eval(θs)
176176
projs = radon(axes(pg)..., ϕs, [θ], [ob]) # many projection views
@@ -224,7 +224,7 @@ p8 = jim(axesf(pg), sp.(proj_fft); prompt=false,
224224
err = maximum(abs, spectrum_slice - proj_fft) / maximum(abs, spectrum_slice)
225225
@assert err < 1e-3
226226
p9 = jim(axesf(pg), 1e3*abs.(proj_fft - spectrum_slice);
227-
title="Difference × 10³", xlabel, ylabel, prompt=false)
227+
title="|Difference| × 10³", xlabel, ylabel, prompt=false)
228228
jim(p6, p7, p8, p9)
229229

230230

‎src/cuboid.jl

+1-1
Original file line numberDiff line numberDiff line change
@@ -53,7 +53,7 @@ end
5353
# methods
5454

5555

56-
volume(ob::Object3d{Cuboid}) = prod(ob.width)
56+
volume1(::Cuboid) = 1 # volume of unit cube
5757

5858

5959
"""

‎src/cylinder.jl

+1-1
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ cylinder(args... ; kwargs...) = Object(Cylinder(), args...; kwargs...)
2929
# methods
3030

3131

32-
volume(ob::Object3d{Cylinder}) = π * prod(ob.width[1:2]) * ob.width[3]
32+
volume1(::Cylinder) = π
3333

3434

3535
"""

‎src/ellipse.jl

+3
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,9 @@ end
5151
# methods
5252

5353

54+
area1(::Ellipse) = π # area of unit circle
55+
56+
5457
"""
5558
phantom1(ob::Object2d{Ellipse}, (x,y))
5659
Evaluate unit circle at `(x,y)`,

‎src/ellipsoid.jl

+1-1
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@ end
5151
# methods
5252

5353

54-
volume(ob::Object3d{Ellipsoid}) = 4/3 * π * prod(ob.width)
54+
volume1(::Ellipsoid) = 4/3 * π # volume of unit sphere
5555

5656

5757
"""

‎src/gauss2.jl

+3
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,9 @@ Convert FWHM `w` to equivalent Gaussian spread `s` for ``\\exp(-π (x/s)^2)``.
5353
# methods
5454

5555

56+
area1(::Gauss2) = fwhm2spread(1)^2
57+
58+
5659
"""
5760
phantom1(ob::Object2d{Gauss2}, (x,y))
5861
Evaluate unit gauss2 at `(x,y)`,

‎src/gauss3.jl

+1-1
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ gauss3(args... ; kwargs...) = Object(Gauss3(), args...; kwargs...)
3232
# methods
3333

3434

35-
volume(ob::Object3d{Gauss3}) = fwhm2spread(1)^3 * prod(ob.width)
35+
volume1(::Gauss3) = fwhm2spread(1)^3
3636

3737

3838
"""

‎src/object2.jl

+3
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,9 @@ Utilities for 2D objects
77
export phantom, radon, spectrum
88

99

10+
area(ob::Object2d{S}) where S = area1(S()) * prod(ob.width)
11+
12+
1013
# rotate
1114

1215
"""

‎src/object3.jl

+3
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,9 @@ using LazyGrids: ndgrid
77
export phantom, radon, spectrum
88

99

10+
volume(ob::Object3d{S}) where S = volume1(S()) * prod(ob.width)
11+
12+
1013
# rotate
1114

1215
"""

‎src/rect.jl

+3
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,9 @@ end
7373
# methods
7474

7575

76+
area1(::Rect) = 1 # area of unit square
77+
78+
7679
"""
7780
phantom1(ob::Object2d{Rect}, (x,y))
7881
Evaluate unit square at `(x,y)`,

‎src/triangle.jl

+3
Original file line numberDiff line numberDiff line change
@@ -125,6 +125,9 @@ end
125125
# methods
126126

127127

128+
area1(::Triangle) = sqrt(3)/4 # area of unit base equilateral
129+
130+
128131
"""
129132
phantom1(ob::Object2d{Triangle}, (x,y))
130133
Evaluate unit triangle at `(x,y)`,

‎test/cuboid.jl

+2-1
Original file line numberDiff line numberDiff line change
@@ -158,7 +158,8 @@ end
158158
v = (-nv÷2:nv÷2-1) * dv
159159
fu = (-nu÷2:nu÷2-1) / nu / du
160160
fv = (-nv÷2:nv÷2-1) / nv / dv
161-
ϕ = deg2rad.(0:6:180) # * Unitful.rad # todo round unitful Unitful.°
161+
# ϕ = (0:6:180) * deg2rad(1)
162+
ϕ = deg2rad.(0:6:180) # helps coverage of "abs(e2) < eps(T)"
162163
θ = [π/7]
163164
sino = @inferred radon(u, v, ϕ, θ, [ob])
164165

‎test/cylinder.jl

+2-5
Original file line numberDiff line numberDiff line change
@@ -141,12 +141,9 @@ end
141141
v = (-nv÷2:nv÷2-1) * dv
142142
fu = (-nu÷2:nu÷2-1) / nu / du
143143
fv = (-nv÷2:nv÷2-1) / nv / dv
144-
ϕ = deg2rad.(0:6:180) # * Unitful.rad # todo round unitful Unitful.°
144+
ϕ = (0:30:180) * deg2rad(1)
145145
θ = [π/7]
146146
sino = @inferred radon(u, v, ϕ, θ, [ob])
147147

148-
#=
149-
todo projection slice
150-
=#
151-
148+
# todo projection slice
152149
end

‎test/ellipse.jl

+3-4
Original file line numberDiff line numberDiff line change
@@ -99,7 +99,7 @@ end
9999
for ob in obs
100100
nr, dr = 2^4, 0.02m
101101
r = (-nr÷2:nr÷2-1) * dr .+ ob.center[1]
102-
ϕ = deg2rad.(0:20:360)
102+
ϕ = (0:30:360) * deg2rad(1)
103103
@inferred IP._radon(ob, r[1], ϕ[1])
104104
sino1 = @inferred radon([r[1]], [ϕ[1]], [ob])
105105
sino = @inferred radon(r, ϕ, [ob])
@@ -118,7 +118,7 @@ end
118118
ob = shape((4m, 3m), radii, π/6, 1.0f0)
119119
img = @inferred phantom(x, y, [ob])
120120

121-
zscale = 1 / π / prod(radii) # normalize spectra by area
121+
zscale = 1 / (ob.value * IP.area(ob)) # normalize spectra by area
122122
fx = (-M÷2:M÷2-1) / M / dx
123123
fy = (-N÷2:N÷2-1) / N / dy
124124
X = myfft(img) * dx * dy * zscale
@@ -145,8 +145,7 @@ end
145145
nr = 2^10
146146
r = (-nr÷2:nr÷2-1) * dr
147147
fr = (-nr÷2:nr÷2-1) / nr / dr
148-
ϕ = deg2rad.(0:360) # * Unitful.rad # todo round unitful Unitful.°
149-
# ϕ = deg2rad.((0:360)°) # not yet due to Unitful issue
148+
ϕ = (0:30:360) * deg2rad(1)
150149
sino = @inferred radon(r, ϕ, [ob])
151150

152151
ia = argmin(abs.(ϕ .- deg2rad(55)))

‎test/ellipsoid.jl

+1-1
Original file line numberDiff line numberDiff line change
@@ -144,7 +144,7 @@ end
144144
v = (-nv÷2:nv÷2-1) * dv
145145
fu = (-nu÷2:nu÷2-1) / nu / du
146146
fv = (-nv÷2:nv÷2-1) / nv / dv
147-
ϕ = deg2rad.(0:6:180) # * Unitful.rad # todo round unitful Unitful.°
147+
ϕ = (0:30:180) * deg2rad(1)
148148
θ = [π/7]
149149
sino = @inferred radon(u, v, ϕ, θ, [ob])
150150

‎test/gauss2.jl

+2-3
Original file line numberDiff line numberDiff line change
@@ -108,7 +108,7 @@ end
108108
ob = shape((2m, 3m), width, π/6, 1.0f0)
109109
img = @inferred phantom(x, y, [ob])
110110

111-
zscale = 1 / IP.fwhm2spread(1)^2 / prod(width) # normalize spectra by area
111+
zscale = 1 / (ob.value * IP.area(ob)) # normalize spectra by area
112112
fx = (-M÷2:M÷2-1) / M / dx
113113
fy = (-N÷2:N÷2-1) / N / dy
114114
X = myfft(img) * dx * dy * zscale
@@ -135,8 +135,7 @@ end
135135
nr = 2^10
136136
r = (-nr÷2:nr÷2-1) * dr
137137
fr = (-nr÷2:nr÷2-1) / nr / dr
138-
ϕ = deg2rad.(0:360) # * Unitful.rad # todo round unitful Unitful.°
139-
# ϕ = deg2rad.((0:180)°) # not yet due to Unitful issue
138+
ϕ = (0:30:360) * deg2rad(1)
140139
sino = @NOTinferred radon(r, ϕ, [ob])
141140

142141
ia = argmin(abs.(ϕ .- deg2rad(55)))

‎test/gauss3.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -138,11 +138,11 @@ end
138138
v = (-nv÷2:nv÷2-1) * dv
139139
fu = (-nu÷2:nu÷2-1) / nu / du
140140
fv = (-nv÷2:nv÷2-1) / nv / dv
141-
ϕ = deg2rad.(0:6:180) # * Unitful.rad # todo round unitful Unitful.°
141+
ϕ = (0:30:180) * deg2rad(1)
142142
θ = [π/7]
143143
sino = @inferred radon(u, v, ϕ, θ, [ob])
144-
end
145144

146145
#=
147146
todo projection slice
148147
=#
148+
end

‎test/rect.jl

+2-3
Original file line numberDiff line numberDiff line change
@@ -98,7 +98,7 @@ end
9898
ob = shape((4m, 3m), width, π/6, 1.0f0)
9999
img = @inferred phantom(x, y, [ob])
100100

101-
zscale = 1 / prod(width) # normalize spectra by area
101+
zscale = 1 / (ob.value * IP.area(ob)) # normalize spectra by area
102102
fx = (-M÷2:M÷2-1) / M / dx
103103
fy = (-N÷2:N÷2-1) / N / dy
104104
X = myfft(img) * dx * dy * zscale
@@ -125,8 +125,7 @@ end
125125
nr = 2^10
126126
r = (-nr÷2:nr÷2-1) * dr
127127
fr = (-nr÷2:nr÷2-1) / nr / dr
128-
ϕ = deg2rad.(0:360) # * Unitful.rad # todo round unitful Unitful.°
129-
# ϕ = deg2rad.((0:360)°) # not yet due to Unitful issue
128+
ϕ = (0:30:360) * deg2rad(1)
130129
sino = @inferred radon(r, ϕ, [ob])
131130

132131
ia = argmin(abs.(ϕ .- deg2rad(55)))

0 commit comments

Comments
 (0)