Divide vesselPath into smaller functions

src/Vessel.jl

@@ -1,9 +1,76 @@
 # This file is originally based on Matlab code written by Christine Droigk
+"""
+  appendRoute!(route, stepsize, angle_xy, angle_xz)
+
+Append a new point to the route of the vessel.
+"""
+function appendRoute!(route, stepsize, angle_xy, angle_xz)
+    push!(route, route[end] .+ stepsize .* (cos(angle_xy)*cos(angle_xz), 
+                                            sin(angle_xy)*cos(angle_xz), 
+                                            sin(angle_xz)))
+end
+
+"""
+  changeDirection!(route, N, stepsize, angle_xy, angle_xz, change_prob, max_change, rng)
+
+Simulate if and how the vessel changes its direction.
+"""
+function changeDirection!(route, N, stepsize, angle_xy, angle_xz, change_prob, max_change, rng)
+  if (rand(rng, Float64) < change_prob) && (length(route) > 2)# if directional change of the vessel
+    # randomly select new angles
+    change_angle_xy = pi*(max_change-2*max_change*rand(rng,Float64))
+    change_angle_xz = pi*(max_change-2*max_change*rand(rng,Float64))
+    # create range such that change is gradually applied
+    step_angle_xy = range(0, change_angle_xy, length=20)
+    step_angle_xz = range(0, change_angle_xz, length=20)
+
+    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
+          appendRoute!(route, stepsize, angle_xy + step_angle_xy[i], angle_xz + step_angle_xz[i])
+        elseif i>length(step_angle_xy) &&  i<=length(step_angle_xz) # only xz-angle changes
+          appendRoute!(route, stepsize, angle_xy + change_angle_xy, angle_xz + step_angle_xz[i])
+        elseif i<=length(step_angle_xy) &&  i>length(step_angle_xz) # only xy-angle changes
+          appendRoute!(route, stepsize, angle_xy + step_angle_xy[i], angle_xz + change_angle_xz)
+        end
+      end
+    end
+    # set current angles to new angles
+    angle_xy = angle_xy + change_angle_xy
+    angle_xz = angle_xz + change_angle_xz
+  else
+      # if no directional change
+      num_step = 15
+      for _=1:num_step
+        appendRoute!(route, stepsize, angle_xy, angle_xz)
+      end
+  end
+  return angle_xy, angle_xz
+end
+
+"""
+  getDiameterRoute(route, diameter, change_diameter_splitting, splitnr)
+
+Compute the diameter anlong the route of the vessel.
+"""
+function getDiameterRoute(route, diameter, change_diameter_splitting, splitnr)
+  if splitnr>1
+    # for the case that the vessel leaves the image during a split the change is set to 1
+    # otherwise `range` throws an error since start and end do not match but length is 1
+    change_diameter_splitting_ = length(route) > 1 ? change_diameter_splitting : 1
+    # gradually change diameter from old value to current one
+    diameter_route = collect(range((1/change_diameter_splitting_)*diameter, diameter, length=length(route)))
+  else
+    # no change if vessel did not already split once
+    diameter_route = diameter*ones(length(route))
+  end
+  return diameter_route
+end
 
 """
-vesselPath(N::NTuple{3,Int}; start, angle_xy, angle_xz, diameter, split_prob, change_prob, 
-  max_change, splitnr, max_number_splits, stepsize, change_diameter_splitting, split_prob_factor, 
-  change_prob_increase, rng)
+  vesselPath(N::NTuple{3,Int}; start, angle_xy, angle_xz, diameter, split_prob, change_prob, 
+    max_change, splitnr, max_number_splits, stepsize, change_diameter_splitting, split_prob_factor, 
+    change_prob_increase, rng)
 
 ### Input parameters:
 * N: Image size, given as a 3 tuple
@@ -48,92 +115,44 @@ function vesselPath(N::NTuple{3,Int};
   push!(route, start)
   while all( 1 .<= route[end] .< N) # while the route of the vessel is inside the image area
 
-    if (rand(rng, Float64) < change_prob) && (length(route) > 2)# 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 = range(0, change_angle_xy, length=20) # to obtain a piecewise change of angle
-        step_angle_xz = range(0, change_angle_xz, length=20)
-
-        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])*cos(angle_xz + step_angle_xz[i]), 
-                                 sin(angle_xy + step_angle_xy[i])*cos(angle_xz + step_angle_xz[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)*cos(angle_xz + step_angle_xz[i]), 
-                                 sin(angle_xy + change_angle_xy)*cos(angle_xz + step_angle_xz[i]), 
-                                 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])*cos(angle_xz + change_angle_xz), 
-                                 sin(angle_xy+step_angle_xy[i])*cos(angle_xz + change_angle_xz), 
-                                 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
-        num_step = 15
-        for _=1:num_step
-          push!(route, route[end] .+ stepsize .* (cos(angle_xy)*cos(angle_xz), 
-                                                sin(angle_xy)*cos(angle_xz), 
-                                                sin(angle_xz))) # use old angle
-        end
-    end
+    angle_xy, angle_xz =  changeDirection!(route, N, stepsize, angle_xy, angle_xz, change_prob, max_change, rng)
 
+    # Splitting part
     if rand(rng, Float64) < split_prob && (length(route) > 3) && (splitnr ≤ max_number_splits) # if vessel is splitting
         min_angle_diff = 0.15*pi
-        angle_diff = rand(rng, Float64)*(pi/2 - min_angle_diff) + min_angle_diff # throw dice for angle between the resulting two vessel parts
-        angle_diff_z = rand(rng, Float64)*(pi/2 - min_angle_diff) + min_angle_diff # 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
+        # Randomly select the angle between the resulting two vessel parts
+        angle_diff_xy = rand(rng, Float64)*(pi/2 - min_angle_diff) + min_angle_diff
+        angle_diff_xz = rand(rng, Float64)*(pi/2 - min_angle_diff) + min_angle_diff
+        # Randomly select the part of the angle related to the first division
+        part_a_xy = rand(rng, Float64)
+        part_a_xz = rand(rng, Float64)
 
         # 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,
+        routeA, diameterA = vesselPath(N; start=route[end], angle_xy=angle_xy-part_a_xy*angle_diff_xy, 
+            angle_xz=angle_xz-part_a_xz*angle_diff_xz, diameter=diameter*change_diameter_splitting,
             split_prob=split_prob*split_prob_factor, change_prob=change_prob+change_prob_increase, max_change, splitnr=splitnr+1, 
             max_number_splits, stepsize, change_diameter_splitting, split_prob_factor, change_prob_increase, 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, 
+        routeB, diameterB = vesselPath(N; start=route[end], angle_xy=angle_xy+(1-part_a_xy)*angle_diff_xy,  
+            angle_xz=angle_xy+(1-part_a_xz)*angle_diff_xz, diameter=diameter*change_diameter_splitting, 
             split_prob=split_prob*split_prob_factor, change_prob=change_prob+change_prob_increase, max_change, splitnr=splitnr+1, 
             max_number_splits, stepsize, change_diameter_splitting, split_prob_factor, change_prob_increase, 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
+        diameter_route = getDiameterRoute(route, diameter, change_diameter_splitting, splitnr)
         append!(route, routeA, routeB)
         append!(diameter_route, diameterA, diameterB)
         return route, diameter_route
     end
   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
-
+  diameter_route = getDiameterRoute(route, diameter, change_diameter_splitting, splitnr)
   return route, diameter_route
 end
 
 
 """
-vesselPhantom(N::NTuple{3,Int}; oversampling=2, kargs...)
+  vesselPhantom(N::NTuple{3,Int}; oversampling=2, kargs...)
 
 ### Input parameters:
 * N: Image size, given as a 3 tuple