initial commit

.github/workflows/ci.yml

@@ -0,0 +1,64 @@
+name: CI
+on:
+  pull_request:
+    branches:
+      - master
+  push:
+    branches:
+      - master
+    tags: '*'
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.8' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
+          - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
+          # - 'nightly'
+        os:
+          - ubuntu-latest
+          - windows-latest
+        threads:
+          - '4'
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+        env:
+          JULIA_NUM_THREADS: ${{ matrix.threads }}
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info
+  docs:
+    name: Documentation
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: '1'
+      - uses: julia-actions/julia-buildpkg@latest
+      - uses: julia-actions/julia-docdeploy@latest
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}

Project.toml

@@ -3,8 +3,19 @@ uuid = "cd4f3286-b756-4d71-b99f-a4358d548556"
 authors = ["Tobias Knopp <[email protected]> and contributors"]
 version = "0.1.0"
 
+[deps]
+ImagePhantoms = "71a99df6-f52c-4da1-bd2a-69d6f37f3252"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
+Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
 [compat]
 julia = "1"
+ImagePhantoms = "0.6"
+StableRNGs = "1.0"
+ImageFiltering = "0.7"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/Shape.jl

@@ -0,0 +1,43 @@
+
+#= Based in ImagePhantoms.jl
+function makeEllipse(N::NTuple{D,Int}, radius, shift, rot) where D
+
+  ranges = ntuple(d-> 1:N[d], D)
+  ob = ellipsoid(Float32.(Tuple(shift .+ (N.÷2)) ), Float32.(Tuple(radius)), (Float32(rand(Float32)*2*pi),Float32(rand(Float32)*2*pi*0) ), 1.0f0)
+  img = phantom(ranges..., [ob], 2)
+  return img
+end
+=#
+
+function singleEllipsoid(N::NTuple{D,Int}, radius, shift, rot) where D
+
+  I = zeros(N)
+  N_ = collect(Tuple(N))
+  R = N_ .* radius
+
+  for i in CartesianIndices(N)
+    r = collect(Tuple(i)) 
+    I[i] = norm(  (rot*(r .- (N.÷2)) ./ R) .- shift ./radius)  <= 1.0
+  end
+
+  return I
+end
+
+function ellipsoidPhantom(N::NTuple{D,Int}; rng::AbstractRNG = GLOBAL_RNG,
+                                            numObjects::Int = rand(rng, 5:10)) where D
+  img = zeros(N)
+  for m=1:numObjects
+    rot = rand(rng, RotMatrix{D})
+    radius = rand(rng, D)*0.3
+    shift = rand(rng, D)*0.3 
+    value = (rand(rng, 1)).^2
+    kernelWidth = ntuple(d-> rand(rng)*N[1] / 20, D)
+    if shift > radius
+      shift = radius
+    end
+    P = singleEllipsoid(N, radius, shift, rot)
+    P = imfilter(P, Kernel.gaussian(kernelWidth))
+    img += value.*P
+  end
+  return img
+end

src/TrainingPhantoms.jl

@@ -1,5 +1,16 @@
 module TrainingPhantoms
 
-# Write your package code here.
+using ImagePhantoms
+using StableRNGs
+using Random
+using Random: GLOBAL_RNG
+using Rotations
+using LinearAlgebra
+using ImageFiltering
+
+export vesselPhantom, ellipsoidPhantom
+
+include("Vessel.jl")
+include("Shape.jl")
 
 end

src/Vessel.jl

@@ -0,0 +1,138 @@
+# This file is originally based on Matlab code written by Christine Droigk
+
+"""
+vesselPath(N::NTuple{3,Int}; start, angle_xy, angle_xz, diameter, split_prob, change_prob, max_change, splitnr)
+
+Input parameters:
+* start: starting point given as a 3x1 vector
+* angle_xy: angle in radians describing the starting angle with which the
+* vessel runs. 
+* angle_xz: same, but angle from xy-plane to vessel's z
+* diameter: starting diameter of vessel
+* split_prob: probability for a splitting of the vessel into two vessel segments. Values between 0 and 1.
+* change_prob: probability for directional change of the vessel route. Values between 0 and 1.
+* max_change: max_change * pi specifies the maximum direction-change angle.
+* splitnr: used for recursive call of the function. For the first call set it to 1. 
+* N: Image size, given as a 3x1 vector
+ 
+Output:
+* route: A length N vector containing the points 3D of the route of the vessel. The
+* length N depends on the random route.
+* diameter_route: A length N vector containing the diameter of the vessel at the positions of the route.
+"""
+function vesselPath(N::NTuple{3,Int}; 
+                             start, 
+                             angle_xy, 
+                             angle_xz, 
+                             diameter, 
+                             split_prob, 
+                             change_prob, 
+                             max_change, 
+                             splitnr,
+                             rng::AbstractRNG = GLOBAL_RNG)
+
+  stepsize = 0.25
+  change_diameter_splitting = 4/5 # Indicates by how much the diameter decreases when the vessel is divided
+  route = NTuple{3,Float64}[]
+  push!(route, start)
+
+  counter = 1;
+  while all( 1 .<= route[end] .< N) # while the route of the vessel is inside the image area
+
+    if rand(rng, Float64) < change_prob # if directional change of the vessel
+        change_angle_xy = pi*(max_change-2*max_change*rand(rng,Float64)) # throw dice for the angles of directional change. For more or less variations replace (0.1-0.2*rand)
+        change_angle_xz = pi*(max_change-2*max_change*rand(rng,Float64)) # same for other angle
+        step_angle_xy = 0:(sign(change_angle_xy)*0.1*stepsize):change_angle_xy # to obtain a piecewise change of angle
+        step_angle_xz = 0:(sign(change_angle_xz)*0.1*stepsize):change_angle_xz
+
+        for i = 1:max(length(step_angle_xy),length(step_angle_xz)) # until the largest change of the two angles is reached
+          if all( 1 .<= route[end] .< N) # check whether image boundaries are reached
+            if i<=length(step_angle_xy) &&  i<=length(step_angle_xz) # both angles are still changing
+              push!(route, route[end] .+ 
+                    stepsize .* (cos(angle_xy+step_angle_xy[i]), sin(angle_xy+step_angle_xy[i]), sin(angle_xz + step_angle_xz[i])) )
+            elseif i>length(step_angle_xy) &&  i<=length(step_angle_xz) # only xz-angle changes
+              push!(route, route[end] .+ 
+                    stepsize .* (cos(angle_xy+change_angle_xy), sin(angle_xy+change_angle_xy), sin(angle_xz + step_angle_xz[i])) )
+            elseif i<=length(step_angle_xy) &&  i>length(step_angle_xz) # only xy-angle changes
+              push!(route, route[end] .+ 
+                    stepsize .* (cos(angle_xy+step_angle_xy[i]), sin(angle_xy+step_angle_xy[i]), sin(angle_xz + change_angle_xz)) )
+            end
+          end
+        end
+        # set current angles to new angles
+        angle_xy = angle_xy + change_angle_xy # set current angles to new angles
+        angle_xz = angle_xz + change_angle_xz
+    else
+        # if no directional change
+        push!(route, route[end] .+ stepsize .* (cos(angle_xy), sin(angle_xy), sin(angle_xz))) # use old angle
+    end
+
+
+    if rand(rng, Float64) < split_prob && counter > 5 # if vessel is splitting
+        angle_diff = rand(rng, Float64)*pi/2 # throw dice for angle between the resulting two vessel parts
+        angle_diff_z = rand(rng, Float64)*pi/2 # same for z-angle 
+
+        part_a = rand(rng, Float64) # Part of the angle related to the first division
+        part_a_z = rand(rng, Float64) # Same for z-angle
+
+        # Call recursively this function for each of the two vessel
+        # segments with adjusted angles and diameter. Here, an adjustment
+        # of the splitting probabilities and directional change
+        # probabilities is made.
+        routeA, diameterA = vesselPath(N; start=route[end], angle_xy=angle_xy-part_a*angle_diff, 
+            angle_xz=angle_xz-part_a_z*angle_diff_z, diameter=diameter*change_diameter_splitting,
+            split_prob=split_prob/2, change_prob=change_prob+0.01, max_change, splitnr=splitnr+1, rng)
+        routeB, diameterB = vesselPath(N; start=route[end], angle_xy=angle_xy+(1-part_a)*angle_diff,  
+            angle_xz=angle_xy+(1-part_a_z)*angle_diff_z, diameter=diameter*change_diameter_splitting, 
+            split_prob=split_prob/2, change_prob=change_prob+0.01, max_change, splitnr=splitnr+1, rng)     
+
+        if splitnr>1
+          # set diameter for the splitting parts
+          change_diameter_splitting_ = length(route) > 1 ? change_diameter_splitting : 1
+          diameter_route = collect(range((1/change_diameter_splitting_)*diameter, diameter, length=length(route)))
+        else
+          diameter_route = diameter*ones(length(route))
+        end
+        append!(route, routeA, routeB)
+        append!(diameter_route, diameterA, diameterB)
+        return route, diameter_route
+    end
+
+    counter = counter+1
+  end
+
+  if splitnr>1
+    change_diameter_splitting_ = length(route) > 1 ? change_diameter_splitting : 1
+    diameter_route = collect(range((1/change_diameter_splitting_)*diameter, diameter, length=length(route)))
+  else
+    diameter_route = diameter*ones(length(route))
+  end
+
+  return route, diameter_route
+end
+
+
+"""
+vesselPhantom(N::NTuple{3,Int}; oversampling=2, kargs...)
+
+Example usage:
+
+  using GR, TrainingPhantoms, StableRNGs
+  im = vesselPhantom((40,40,40); start=(1, 20, 20), angle_xy=0.0, angle_xz=0.0, 
+                                diameter=2, split_prob=0.5, change_prob=0.5, max_change=0.2, splitnr=1,
+                                rng = StableRNG(1));
+  isosurface(im, isovalue=0.2, rotation=110, tilt=40)
+
+"""
+function vesselPhantom(N::NTuple{3,Int}; oversampling=2, kargs...)
+  Q = zeros(N);
+
+  route, diameter_route = vesselPath(N; kargs...)
+
+  obs = [ ellipsoid( Float32.(route[i]), Float32.(ntuple(_->diameter_route[i],3)), 
+                              (0,0,0), 1.0f0) for i=1:length(route)]
+  ranges = ntuple(d-> 1:N[d], 3)
+  img = phantom(ranges..., obs, oversampling)
+  img[img .> 1] .= 1
+  return img
+end

test/runtests.jl

@@ -1,6 +1,24 @@
 using TrainingPhantoms
+using StableRNGs
 using Test
 
 @testset "TrainingPhantoms.jl" begin
-    # Write your tests here.
+  N = (40,40,40)
+
+  # Vessel Phantom
+  im = vesselPhantom(N; start=(1, 20, 20), angle_xy=0.0, angle_xz=0.0, 
+                                diameter=2, split_prob=0.5, change_prob=0.5, max_change=0.2, splitnr=1,
+                                rng = StableRNG(1));
+
+  im2 = vesselPhantom(N; start=(1, 20, 20), angle_xy=0.0, angle_xz=0.0, 
+                                diameter=2, split_prob=0.5, change_prob=0.5, max_change=0.2, splitnr=1,
+                                rng = StableRNG(1));
+  @test im ≈ im2
+
+   # Ellipsoid Phantom
+  im = ellipsoidPhantom(N; rng=StableRNG(1))
+  im2 = ellipsoidPhantom(N; rng=StableRNG(1))
+  @test im ≈ im2
+
+  #isosurface(im, isovalue=0.2, rotation=110)
 end