From 39bd244dec9165909638a1f2ca4fab8cd47515e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Paul=20J=C3=BCr=C3=9F?= Date: Mon, 26 Feb 2024 18:49:06 +0100 Subject: [PATCH] Update vesselPath to make parameter selection more stable --- src/Vessel.jl | 63 +++++++++++++++++++++++++++++++-------------------- 1 file changed, 38 insertions(+), 25 deletions(-) diff --git a/src/Vessel.jl b/src/Vessel.jl index 0443e5c..54825d4 100644 --- a/src/Vessel.jl +++ b/src/Vessel.jl @@ -1,9 +1,12 @@ # 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) +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 * start: starting point given as a 3x1 vector * angle_xy: angle in radians describing the starting angle with which the * vessel runs. @@ -13,7 +16,12 @@ Input parameters: * 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 +* max_number_splits: maximum number of splits of the vessel. +* stepsize: stepsize of the vessel. +* change_diameter_splitting: Indicates by how much the diameter decreases when the vessel splits +* split_prob_factor: Factor by which the split probability `split_prob` is multiplied when the vessel splits +* change_prob_increase: Increase of the change probability `change_prob` when the vessel splits +* rng: Random number generator Output: * route: A length N vector containing the points 3D of the route of the vessel. The @@ -29,21 +37,22 @@ function vesselPath(N::NTuple{3,Int}; change_prob, max_change, splitnr, + max_number_splits = Inf, + stepsize = max(N...)/200, + change_diameter_splitting = 0.6, + split_prob_factor = 0.5, + change_prob_increase = 0.01, 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 + 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 = 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 + 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 @@ -70,15 +79,18 @@ function vesselPath(N::NTuple{3,Int}; angle_xz = angle_xz + change_angle_xz else # if no directional change - push!(route, route[end] .+ stepsize .* (cos(angle_xy)*cos(angle_xz), + 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 - - 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 + 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 @@ -89,10 +101,12 @@ function vesselPath(N::NTuple{3,Int}; # 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) + 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, - split_prob=split_prob/2, change_prob=change_prob+0.01, max_change, splitnr=splitnr+1, rng) + 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 @@ -105,8 +119,6 @@ function vesselPath(N::NTuple{3,Int}; append!(diameter_route, diameterA, diameterB) return route, diameter_route end - - counter = counter+1 end if splitnr>1 @@ -132,15 +144,16 @@ Example usage: 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) ] +function vesselPhantom(N::NTuple{3,Int}; oversampling=2, rng = GLOBAL_RNG, kargs...) + route, diameter_route = vesselPath(N; rng, kargs...) + # add small sphere for every entry in the route + obs = [ sphere( Float32.(route[i]), Float32(diameter_route[i]), 1.0f0) for i=eachindex(route) ] ranges = ntuple(d-> 1:N[d], 3) img = phantom(ranges..., obs, oversampling) + # filter for smoothing, offset to ensure minimal filter width + filterWidth = (1.0-0.2)*rand(rng) + 0.2 + kernelWidth = ntuple(_ -> filterWidth*N[1] / 20, 3) + img = imfilter(img, Kernel.gaussian(kernelWidth)) img[img .> 1] .= 1 return img end