Skip to content

Commit 817eaa6

Browse files
pauljuersspauljuerss
authored
Merge pull request #13 from JuliaImageRecon/pj/ellipsoidBoundingBox
Add documentation and ellipsoid bounding box
2 parents 99ff7c4 + 1d98873 commit 817eaa6

File tree

11 files changed

+244
-46
lines changed

11 files changed

+244
-46
lines changed

‎.gitignore

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
/Manifest.toml
1+
*Manifest.toml

‎README.md

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,3 +9,20 @@
99
![Lifecycle](https://img.shields.io/badge/lifecycle-dormant-blue.svg) -->
1010
[![Build status](https://github.com/JuliaImageRecon/TrainingPhantoms.jl/workflows/CI/badge.svg)](https://github.com/JuliaImageRecon/TrainingPhantoms.jl/actions)
1111
[![codecov.io](http://codecov.io/github/JuliaImageRecon/TrainingPhantoms.jl/coverage.svg?branch=master)](http://codecov.io/github/JuliaImageRecon/TrainingPhantoms.jl?branch=master)
12+
13+
## Purpose
14+
This Julia package contains a collection of methods that generate different kinds of phantom with the aim of being used to synthesise machine learning datasets. Currently, there exist two generator methods
15+
* [Ellipsoid Phantom Generator](src/Shape.jl): ellipsoids of different size, location, orientation and value should mimic the images seen, e.g., when contrast agents are administered in the form of boluses
16+
* [Vessel Phantom Generator](src/Vessel.jl): randomly generated blood vessels
17+
18+
## Installation
19+
```julia
20+
using Pkg
21+
Pkg.add("TrainingPhantoms")
22+
```
23+
24+
## Examples
25+
To help you get started, we have provided examples for each generator. The can be found in the [examples](examples/) folder.
26+
<!--- ![ellipsoidPhantom2](https://github.com/JuliaImageRecon/TrainingPhantoms.jl/assets/115639115/863b0769-9643-4858-9201-3f94311ab2ba)
27+
![vesselPhantom](https://github.com/JuliaImageRecon/TrainingPhantoms.jl/assets/115639115/79c95f1f-b284-4562-8464-9b12ac3edf7d) --->
28+
![combined](https://github.com/JuliaImageRecon/TrainingPhantoms.jl/assets/115639115/5babfffa-c464-4b32-99b4-da8cd4e12e86)

‎examples/EllipsoidPhantoms/README.md

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
# Ellipsoid Phantom Example
2+
This example generates an ellipsoid phantom and visualizes it using a mean intensity plot and an isosurface volume rendering.
3+
4+
## Installation
5+
To use this code, you must first install `TrainingPhantoms` by running
6+
```julia
7+
using Pkg
8+
Pkg.add("TrainingPhantoms")
9+
```
10+
Afterwards, you can load the package by entering
11+
```julia
12+
using TrainingPhantoms
13+
```
14+
15+
## Execution
16+
Finally, just include the example code
17+
```julia
18+
include(joinpath(dirname(pathof(TrainingPhantoms)),"..","ellipsoidPhantoms", "example.jl"))
19+
```
20+
After the phantom was generated, it will be shown in a new window.

‎examples/EllipsoidPhantoms/example.jl

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
using Pkg
2+
Pkg.activate(joinpath(@__DIR__,".."))
3+
Pkg.instantiate()
4+
using TrainingPhantoms, Statistics
5+
include(joinpath(@__DIR__,"..","helperFunctions","meanIntensityPlotAndIsosurface.jl"))
6+
GLMakie.activate!()
7+
8+
# Generate a vessel phantom
9+
rng = StableRNGs.StableRNG(1234)
10+
N = (51,51,51)
11+
phantom = ellipsoidPhantom(N; rng=rng)
12+
13+
# Visualize the phantom
14+
meanIntensityPlotAndIsosurface(phantom)

‎examples/Project.toml

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
[deps]
2+
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
3+
GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
4+
Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
5+
StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
6+
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
7+
TrainingPhantoms = "cd4f3286-b756-4d71-b99f-a4358d548556"

‎examples/VesselPhantoms/README.md

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
# Vessel Phantom Example
2+
This example generates a vessel phantom and visualizes it using a mean intensity plot and an isosurface volume rendering.
3+
4+
## Installation
5+
To use this code, you must first install `TrainingPhantoms` by running
6+
```julia
7+
using Pkg
8+
Pkg.add("TrainingPhantoms")
9+
```
10+
Afterwards, you can load the package by entering
11+
```julia
12+
using TrainingPhantoms
13+
```
14+
15+
## Execution
16+
Finally, just include the example code
17+
```julia
18+
include(joinpath(dirname(pathof(TrainingPhantoms)),"..","VesselPhantoms", "example.jl"))
19+
```
20+
After the vessel was generated, the phantom will be shown in a new window.

‎examples/VesselPhantoms/example.jl

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
using Pkg
2+
Pkg.activate(joinpath(@__DIR__,".."))
3+
Pkg.instantiate()
4+
using TrainingPhantoms, Statistics
5+
include(joinpath(@__DIR__,"..","helperFunctions","meanIntensityPlotAndIsosurface.jl"))
6+
GLMakie.activate!()
7+
8+
# Generate a vessel phantom
9+
rng = StableRNGs.StableRNG(123)
10+
N = (51,51,51)
11+
phantom = vesselPhantom(N;
12+
start=(1, N[2]/2, N[3]/2), angles=(0.0,0.0),
13+
diameter=2.5, split_prob=0.4, change_prob=0.3,
14+
max_change=0.3, splitnr=1, max_number_splits=1, rng=rng)
15+
16+
# Visualize the phantom
17+
meanIntensityPlotAndIsosurface(phantom)

‎examples/helperFunctions/meanIntensityPlotAndIsosurface.jl

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
import CairoMakie, GLMakie, Makie
2+
import Statistics
3+
4+
function meanIntensityPlotAndIsosurface(data::AbstractArray{T,3}; func::Function=mean, fig=CairoMakie.Figure(size=(1150,300))) where T
5+
CairoMakie.empty!(fig) # ensure figure is empty
6+
imSize = size(data)
7+
xLabelVec = ["x", "x", "y"]
8+
yLabelVec = ["y", "z", "z"]
9+
dimVec = [3,2,1]
10+
sizeVec = [imSize[1:2], imSize[[1,3]], imSize[2:3]]
11+
axConfig = Dict(:aspect=>1, :ylabelrotation=>0, :xlabelsize=>20, :ylabelsize=>20)
12+
for i=1:3
13+
ax = CairoMakie.Axis(fig[1,i]; xlabel=xLabelVec[i], ylabel=yLabelVec[i], axConfig...)
14+
CairoMakie.heatmap!(ax, reshape(func(data[:,:,:], dims=dimVec[i]), sizeVec[i]), colormap=:inferno)
15+
CairoMakie.hidedecorations!(ax, label=false)
16+
end
17+
fgl = GLMakie.Figure(size=(300,300))
18+
axgl = GLMakie.Axis3(fgl[1,1], aspect=:data, azimuth=-pi/2-0.2)
19+
GLMakie.volume!(axgl, data, algorithm=:iso, isorange=0.15, isovalue=0.3, colormap=:viridis, colorrange=[0.0,0.2])
20+
GLMakie.hidedecorations!(axgl, label=false)
21+
GLMakie.hidespines!(axgl)
22+
Makie.update_state_before_display!(fgl)
23+
cb = Makie.colorbuffer(fgl.scene, backend=GLMakie, scalefactor=5)
24+
axleft, _ = CairoMakie.image(fig[1,4], cb[300:end-300,300:end-300]', axis=(yreversed=true, aspect=1))
25+
CairoMakie.hidedecorations!(axleft)
26+
CairoMakie.hidespines!(axleft)
27+
return fig
28+
end

‎src/Shape.jl

Lines changed: 86 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@ end
1313

1414
"""
1515
scaleValue(x, a, b)
16+
1617
Scale a value `x` that is assumed to lie in the interval [0,1]
1718
to the interval [a,b] for given values a and b.
1819
"""
@@ -23,33 +24,104 @@ function scaleValue(x::T, a::Real, b::Real=one(T)) where {T <: Real}
2324
end
2425

2526
"""
26-
ellipsoidPhantom(
27-
N::NTuple{D,Int};
28-
rng::AbstractRNG = GLOBAL_RNG,
29-
numObjects::Int = rand(rng, 5:10),
30-
minRadius::Real=1.0,
31-
minValue::Real=0.1,
32-
allowOcclusion::Bool=false
33-
)
27+
getRotationMatrix(D::Int, rotAngles::NTuple{Dr, <:Real}) where Dr
28+
29+
Get the rotation matrix for the given rotation angles.
30+
31+
# Arguments
32+
- `D`: dimension of the object
33+
- `rotAngles`: rotation angles
34+
"""
35+
function getRotationMatrix(D::Int, rotAngles::NTuple{Dr, <:Real}) where Dr
36+
R = zeros(Float64, (D,D))
37+
for i=1:D
38+
# readout columns of rotation matrix with unit vectors
39+
R[:,i] = [rotate_vector(ntuple(j -> (i==j ? 1.0 : 0.0), D), rotAngles)...]
40+
end
41+
return R
42+
end
43+
44+
"""
45+
getEllipsoidMatrix(radius::NTuple{Dr, <:Real}, rotAngles::NTuple{Da, <:Real}) where {Dr, Da}
46+
47+
Get the matrix that describes the ellipsoid in the rotated coordinate system.
48+
49+
# Arguments
50+
- `radius`: radii of the ellipsoid
51+
- `rotAngles`: rotation angles
52+
"""
53+
function getEllipsoidMatrix(radius::NTuple{Dr, <:Real}, rotAngles::NTuple{Da, <:Real}) where {Dr, Da}
54+
R = getRotationMatrix(Dr, rotAngles)
55+
A = Diagonal([1/radius[i]^2 for i=1:Dr])
56+
return transpose(R)*A*R
57+
end
58+
59+
"""
60+
ellipsoidBoundingBox(radius::NTuple{Dr, <:Real}, rotAngles::NTuple{Da, <:Real}) where {Dr, Da}
61+
62+
Get the bounding box of the specified ellipsoid.
63+
64+
# Arguments
65+
- `radius`: radii of the ellipsoid
66+
- `rotAngles`: rotation angles
67+
"""
68+
function ellipsoidBoundingBox(radius::NTuple{Dr, <:Real}, rotAngles::NTuple{Da, <:Real}) where {Dr, Da}
69+
B = getEllipsoidMatrix(radius, rotAngles) |> inv
70+
return sqrt.(diag(B))
71+
end
72+
73+
# rotate vector v by the given rotation angles
74+
rotate_vector(::NTuple{D, <:Real}, ::NTuple{Da, <:Real}) where {D, Da} = throw(ArgumentError("`rotate_vector` currently does not support rotations for vectors of size $D with $Da angles."))
75+
rotate_vector(v::NTuple{3, <:Real}, rotAngles::NTuple{3, <:Real}) = ImagePhantoms.Rxyz_inv(v, rotAngles...)
76+
rotate_vector(v::NTuple{2, <:Real}, rotAngles::NTuple{1, <:Real}) = ImagePhantoms.rotate2d(v, rotAngles...)
77+
78+
# Get the number of rotational degrees of freedom (according to Euclidean group)
79+
rotDOF(::NTuple{D,<:Real}) where D = Int(D*(D-1)/2)
3480

81+
"""
82+
ellipsoidPhantom(
83+
N::NTuple{D,Int};
84+
rng::AbstractRNG = GLOBAL_RNG,
85+
numObjects::Int = rand(rng, 5:10),
86+
minRadius::Real=1.0,
87+
minValue::Real=0.1,
88+
allowOcclusion::Bool=false
89+
)
90+
91+
Generate phantom composed of mulitple ellipsoids.
92+
93+
# Arguments
3594
- `N`: size of the phantom image
3695
- `rng`: random number generator
3796
- `numObjects`: number of ellipses to generate
3897
- `minRadius`: minimal radius in pixel
3998
- `minValue`: minimal value of a single ellipse
4099
- `allowOcclusion`: if `true` ellipse overshadows smaller values at its location, i.e.,
41100
new ellipses are not simply added to the exisiting object, instead the maximum is selected
101+
- `pixelMargin`: minimal distance of the object to the edge of the image
102+
103+
# Examples
104+
105+
using GLMakie, TrainingPhantoms, StableRNGs
106+
107+
im = ellipsoidPhantom((51,51); rng=StableRNG(1))
108+
f = Figure(size=(300,300))
109+
ax = Axis3(f[1,1], aspect=:data)
110+
volume!(ax, im, algorithm=:iso, isorange=0.13, isovalue=0.3, colormap=:viridis, colorrange=[0.0,0.2])
42111
"""
43112
function ellipsoidPhantom(N::NTuple{D,Int}; rng::AbstractRNG = GLOBAL_RNG,
44-
numObjects::Int = rand(rng, 5:10), minRadius::Real=1.0,
45-
minValue::Real=0.1, allowOcclusion::Bool=false) where D
113+
numObjects::Int = rand(rng, 5:10), minRadiusPixel::Real=1.0,
114+
maxRadiusPercent::Real=0.3, maxShiftPercent::Real=0.6,
115+
minValue::Real=0.1, allowOcclusion::Bool=false, pixelMargin::Real=1) where D
46116
img = zeros(N)
47117
@debug numObjects
48118
for m=1:numObjects
49-
# in 2D there is just one rotational degree of freedom
50-
rotAngles = ntuple(_ -> 2π*rand(rng), D == 2 ? 1 : D)
51-
radius = N .* scaleValue.(ntuple(_ -> 0.3*rand(rng), D), minRadius./N)
52-
shift = N .* ntuple(_ -> 0.6*(rand(rng)-0.5), D)
119+
rotAngles = ntuple(_ -> 2π*rand(rng), rotDOF(N))
120+
radius = N .* scaleValue.(ntuple(_ -> rand(rng), D), minRadiusPixel./N, maxRadiusPercent)
121+
shift = N .* ntuple(_ -> maxShiftPercent*(rand(rng)-0.5), D)
122+
boundingBox = ellipsoidBoundingBox(radius, rotAngles)
123+
# if the bounding box is too close to the edge of the image, shift the object to have `pixelMargin` distance
124+
shift = ntuple(i -> (boundingBox[i]+abs(shift[i]) < (N[i]/2 - pixelMargin)) ? shift[i] : (N[i]/2 - boundingBox[i] - pixelMargin)*sign(shift[i]), D)
53125
value = scaleValue.(rand(rng, 1), minValue)#.^2
54126
kernelWidth = ntuple(_ -> rand(rng)*N[1] / 20, D)
55127
P = singleEllipsoid(N, radius, shift, rotAngles)

‎src/Vessel.jl

Lines changed: 33 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
# This file is originally based on Matlab code written by Christine Droigk
22
"""
3-
appendRoute!(route, stepsize, angles)
3+
appendRoute!(route, stepsize, angles)
44
55
Append a new point to the route of the vessel.
66
"""
@@ -18,7 +18,7 @@ function appendRoute!(route, stepsize, angles::NTuple{1,Float64})
1818
end
1919

2020
"""
21-
changeDirection!(route, N, stepsize, angles, change_prob, max_change, rng)
21+
changeDirection!(route, N, stepsize, angles, change_prob, max_change, rng)
2222
2323
Simulate if and how the vessel changes its direction.
2424
"""
@@ -47,7 +47,7 @@ function changeDirection!(N::NTuple{D,Int}, route, stepsize, angles, change_prob
4747
end
4848

4949
"""
50-
getDiameterRoute(route, diameter, change_diameter_splitting, splitnr)
50+
getDiameterRoute(route, diameter, change_diameter_splitting, splitnr)
5151
5252
Compute the diameter anlong the route of the vessel.
5353
"""
@@ -66,33 +66,33 @@ function getDiameterRoute(route, diameter, change_diameter_splitting, splitnr)
6666
end
6767

6868
"""
69-
vesselPath(N::NTuple{D,Int}; start, angles, diameter, split_prob, change_prob,
70-
max_change, splitnr, max_number_splits, stepsize, change_diameter_splitting, split_prob_factor,
71-
change_prob_increase, steps_each_change, steps_no_change, rng)
72-
73-
74-
### Input parameters:
75-
* N: Image size, given as a D tuple
76-
* start: starting point given as a D tuple
77-
* angles: D-1 tuple of angles for the vessel's direction (xy in 2D, xy and xz in 3D)
78-
* diameter: starting diameter of vessel
79-
* split_prob: probability for a splitting of the vessel into two vessel segments. Values between 0 and 1.
80-
* change_prob: probability for directional change of the vessel route. Values between 0 and 1.
81-
* max_change: max_change * pi specifies the maximum direction-change angle.
82-
* splitnr: used for recursive call of the function. For the first call set it to 1.
83-
* max_number_splits: maximum number of splits of the vessel.
84-
* stepsize: stepsize of the vessel.
85-
* change_diameter_splitting: Indicates by how much the diameter decreases when the vessel splits
86-
* split_prob_factor: Factor by which the split probability `split_prob` is multiplied when the vessel splits
87-
* change_prob_increase: Increase of the change probability `change_prob` when the vessel splits
88-
* steps_each_change: Number of steps for the change of the vessel direction
89-
* steps_no_change: Number of steps for the case that the vessel does not change its direction
90-
* rng: Random number generator
69+
vesselPath(N::NTuple{D,Int}; start, angles, diameter, split_prob, change_prob,
70+
max_change, splitnr, max_number_splits, stepsize, change_diameter_splitting, split_prob_factor,
71+
change_prob_increase, steps_each_change, steps_no_change, rng)
72+
73+
74+
# Arguments
75+
- `N`: Image size, given as a D tuple
76+
- `start`: starting point given as a D tuple
77+
- `angles`: D-1 tuple of angles for the vessel's direction (xy in 2D, xy and xz in 3D)
78+
- `diameter`: starting diameter of vessel
79+
- `split_prob`: probability for a splitting of the vessel into two vessel segments. Values between 0 and 1.
80+
- `change_prob`: probability for directional change of the vessel route. Values between 0 and 1.
81+
- `max_change`: max_change * pi specifies the maximum direction-change angle.
82+
- `splitnr`: used for recursive call of the function. For the first call set it to 1.
83+
- `max_number_splits`: maximum number of splits of the vessel.
84+
- `stepsize`: stepsize of the vessel.
85+
- `change_diameter_splitting`: Indicates by how much the diameter decreases when the vessel splits
86+
- `split_prob_factor`: Factor by which the split probability `split_prob` is multiplied when the vessel splits
87+
- `change_prob_increase`: Increase of the change probability `change_prob` when the vessel splits
88+
- `steps_each_change`: Number of steps for the change of the vessel direction
89+
- `steps_no_change`: Number of steps for the case that the vessel does not change its direction
90+
- `rng`: Random number generator
9191
92-
Output:
93-
* route: A length N vector containing the D dimensional points of the route of the vessel. The
92+
# Output
93+
- `route`: A length N vector containing the D dimensional points of the route of the vessel. The
9494
length N depends on the random route.
95-
* diameter_route: A length N vector containing the diameter of the vessel at the positions of the route.
95+
- `diameter_route`: A length N vector containing the diameter of the vessel at the positions of the route.
9696
"""
9797
function vesselPath(N::NTuple{D,Int};
9898
start,
@@ -155,16 +155,18 @@ sphereFunction(::NTuple{2, Int}) = circle
155155
sphereFunction(::NTuple{3, Int}) = sphere
156156

157157
"""
158-
vesselPhantom(N::NTuple{D,Int}; oversampling=2, kargs...)
158+
vesselPhantom(N::NTuple{D,Int}; oversampling=2, kargs...)
159+
160+
Generate a random vessel phantom.
159161
160-
### Input parameters:
162+
# Arguments
161163
* N: Image size, given as a D tuple
162164
* oversampling: Oversampling factor for the phantom. Default is 2.
163165
* rng: Random number generator
164166
* kernelWidth: Width (standard deviation) of the Gaussian kernel (in pixel) used for smoothing the phantom. If nothing is given, a random value is chosen.
165167
* kargs...: remaining keyword arguments for `vesselPath`
166168
167-
### Example usage:
169+
# Examples
168170
169171
using GLMakie, TrainingPhantoms, StableRNGs
170172

0 commit comments

Comments
 (0)