Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history 
Updates to comments, correction of long-term bug which caused some V-Systems to fail. Updated grammar descriptions to include shrinkage and expantion of vessels.
  • Loading branch information
@psweens
psweens committed Apr 14, 2023
1 parent 02c30aa commit dfa6069
Show file tree
Hide file tree
Showing 5 changed files with 221 additions and 268 deletions.
101 changes: 46 additions & 55 deletions analyseGrammar.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
DEGREES_TO_RADIANS = pi / 180

def branching_turtle_to_coords(turtle_program, d0, theta=20., phi=20.):

'''
Working with discontinuous paths i.e. tree formation
'+' : postive rotation by (deg)
Expand All @@ -22,87 +23,83 @@ def branching_turtle_to_coords(turtle_program, d0, theta=20., phi=20.):
Returns
tuple: A list of tuples which provide the coordinates of the branches
'''

# Initialize variables
saved_states = list() # Stack for saving and restoring states
stateSize = 10 # Number of elements in each state tuple
dx = 0 # x displacement
dy = 0 # y displacement
dz = 0 # z displacement
lseg = 1. # Length of each segment
rim = 400 # Rim radius for proximity calculation

# Choose the starting direction randomly
saved_states = list()
stateSize = 10
dx = 0
dy = 0
dz = 0
lseg = 1.
rim = 400

startidx = 3#random.randint(1,3)
if startidx == 1:
state = (1., 0.1, 0.1, 0, 0, d0, lseg, dx, dy, dz)
elif startidx == 2:
state = (0.1, 1., 0, 0, 0, d0, lseg, dx, dy, dz)
else:
state = (0.1, 0.1, 1., 0, 0, d0, lseg, dx, dy, dz)

yield state

index = 0 # Keeps track of the index of the command being executed
origin = (0.1, 0.1, 1.) # Origin point for proximity calculation
index = 0
origin = (0.1, 0.1, 1.)

# Loop through each command in the turtle program
for command in turtle_program:
# Get the current state of the turtle
x, y, z, alpha, beta, diam, lseg, dx, dy, dz = state

# If the command is a move forward command

if command.lower() in 'abcdefghijs': # Move forward (matches a-j and A-J)
# segment start


# Start a new segment
if command.islower():
# Get the length and diameter of the segment from the brackets
lseg, tdiam = eval_brackets(index, turtle_program)
# Calculate the displacement vector based on the current direction and segment length
dx, dy, dz = rotate(p=beta*DEGREES_TO_RADIANS,
r=alpha*DEGREES_TO_RADIANS,
v=normalise(np.array([x,y,z]),lseg))
# If the diameter is specified, update the current diameter
if tdiam > 0.0:
diam = tdiam
# Move the turtle forward by the displacement vector
dx, dy, dz = rotate(pitch_angle=beta*DEGREES_TO_RADIANS,
roll_angle=alpha*DEGREES_TO_RADIANS,
vector=normalise(np.array([x,y,z]),lseg))

if tdiam > 0.0: diam = tdiam

#dx, dy, dz, alpha = proximity(state,origin,rim)

x += dx
y += dy
z += dz

# Save the current state
state = (x, y, z, alpha, beta, diam, lseg, dx, dy, dz)

# segment end
yield state

elif command == '+': # Turn clockwise
phi, _ = eval_brackets(index, turtle_program) # Get the angle to rotate by
state = (x, y, z, alpha + phi, beta, diam, lseg, dx, dy, dz) # Update the state with new angle
elif command == '+': # Turn clockwise
phi, _ = eval_brackets(index, turtle_program)
state = (x, y, z, alpha + phi, beta, diam, lseg, dx, dy, dz)

elif command == '-': # Turn counterclockwise
phi, _ = eval_brackets(index, turtle_program) # Get the angle to rotate by
state = (x, y, z, alpha - phi, beta, diam, lseg, dx, dy, dz) # Update the state with new angle
elif command == '-': # Turn counterclockwise
phi, _ = eval_brackets(index, turtle_program)
state = (x, y, z, alpha - phi, beta, diam, lseg, dx, dy, dz)

elif command == '/': # Turn around y-axis
theta, _ = eval_brackets(index, turtle_program) # Get the angle to rotate by
state = (x, y, z, alpha, beta + theta, diam, lseg, dx, dy, dz) # Update the state with new angle
elif command == '/':
theta, _ = eval_brackets(index, turtle_program)
state = (x, y, z, alpha, beta + theta, diam, lseg, dx, dy, dz)

elif command == '[': # Remember current state
saved_states.append(state) # Push the current state onto the stack
elif command == '[': # Remember current state
saved_states.append(state)

elif command == ']': # Return to previous state
state = saved_states.pop() # Pop the previous state from the stack and update the current state
elif command == ']': # Return to previous state
state = saved_states.pop()

# Yield a tuple of NaN values to indicate the start of a new branch
nanValues = []
for i in range(stateSize): nanValues.append(float('nan'))
yield tuple(nanValues)
nanValues = []
for i in range(stateSize): nanValues.append(float('nan'))
yield tuple(nanValues)

x, y, z, alpha, beta, diam, lseg, dx, dy, dz = state # Update the current state with previous state
yield state # Yield the current state
x, y, z, alpha, beta, diam, lseg, dx, dy, dz = state
yield state

index += 1 # Increment the command index
# Note: silently ignore unknown commands
index += 1

def eval_brackets(index, turtle):

"""
Extracts values within brackets in the turtle program and evaluates them to return tuple of values.
Expand Down Expand Up @@ -145,7 +142,6 @@ def eval_brackets(index, turtle):
if neg1 == 1: aa *= -1 # if neg1 flag is set, negate first value
return aa, 0.0


def randomposneg():
"""
Return either 1 or -1 with a 50/50 probability.
Expand All @@ -155,8 +151,6 @@ def randomposneg():
"""
return 1 if random.random() < 0.5 else -1



def raddist(origin, location, shell=80, core=False):
"""
Calculate the distance between two points and check if it falls within a specified range.
Expand All @@ -181,8 +175,6 @@ def raddist(origin, location, shell=80, core=False):
# check if location is outside the shell
return distance > shell



def proximity(state: tuple, origin: np.ndarray, rim: float) -> tuple:
"""
Calculates the proximity of a point to a given origin within a specified rim.
Expand Down Expand Up @@ -244,5 +236,4 @@ def posneg(value: float) -> float:
if value >= 0.:
return 1.
else:
return -1.

return -1.
66 changes: 27 additions & 39 deletions libGenerator.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
import ast
import random
import numpy as np
import random
from numpy import math as npm

default={"k": 3,
Expand All @@ -19,16 +18,12 @@ def setProperties(properties):
Returns:
None
Raises:
None
"""

# Check if properties is None and assign default value
if properties is None:
if properties == None:
properties = default

# Define global variables based on the input dictionary
global k, epsilon, randmarg, sigma, stochparams

k = properties['k']
epsilon = properties['epsilon']
randmarg = properties['randmarg']
Expand All @@ -49,59 +44,52 @@ def calParam(text, params):
'''

txt = text[:]
for i in params:
txt = txt.replace(i, str(params[i]))
try:
result = ast.literal_eval(txt)
return str(params['co'] / result)
except:
return 'Error: Invalid expression'
for i in params: txt = txt.replace(i, str(params[i]))

return str(params['co'] / eval(txt))

def calBifurcation(d0):

'''
'''
Calculates the diameters and angles of bifurcation given an input diameter
Parameters:
Args:
d0 (float): input diameter
Returns:
resp (dict): a dictionary containing the calculated values for d1, d2, d0, th1, th2, and co
'''

resp = {}

# Calculate optimal diameter and adjust if stochastic parameter is set to True

dOpti = d0 / 2 ** (1.0 / k)
if stochparams:
d1 = abs(np.random.normal(dOpti, dOpti / sigma))
else:
d1 = dOpti # Optimal diameter

# Ensure that d1 is not greater than d0
if d1 >= d0:
d1 = dOpti

# Calculate second diameter using the first diameter and the rate of symmetry between daughters
if stochparams: d1 = abs(np.random.normal(dOpti, dOpti / sigma))
else: d1 = dOpti # Optimal diameter

if d1 >= d0: d1 = dOpti # Elimate possibility of d1 being greater than d0

d2 = (d0 ** k - d1 ** k) ** (1.0 / k) # Calculate second diameter
# alpha = abs(np.random.uniform(1., 0.25)) * (d2 / d1) # Rate of symmetry of daughters (=1 symmetrical ?)
alpha = d2 / d1
d2 = (d0 ** k - d1 ** k) ** (1.0 / k)

# Equations which mimic bifurcation angles in the human body (Liu et al. (2010) and Zamir et al. (1988))

'''
Equations which mimic bifurcation angles in the human body
Liu et al. (2010) and Zamir et al. (1988)
'''
xtmp = (1 + alpha * alpha * alpha) ** (4.0 / 3) + 1 - alpha ** 4
xtmpb = 2 * ((1 + alpha * alpha * alpha ) ** (2.0 / 3))
a1 = npm.acos(xtmp / xtmpb)

xtmp = (1 + alpha * alpha * alpha) ** (4.0 / 3) + (alpha ** 4) - 1
xtmpb = 2 * alpha * alpha * ((1 + alpha * alpha * alpha) ** (2.0/3))
a2 = npm.acos(xtmp / xtmpb)

# Store calculated values in a dictionary and return

resp["d1"] = d1
resp["d2"] = d2
resp["d0"] = d0
resp["th1"] = a1 * 180 / npm.pi
resp["th2"] = a2 * 180 / npm.pi
resp["co"] = getLength(d0)

return resp

def getLength(d0):
Expand All @@ -114,6 +102,6 @@ def getLength(d0):
Returns:
float: The length of the branch.
"""
c0 = d0 * epsilon # calculate c0 based on the parent branch diameter and epsilon
return np.random.uniform(c0, randmarg * 2)
c0 = d0 * epsilon
# abs(np.random.normal(50,10))
return np.random.uniform(c0 - randmarg, c0 + randmarg)
Loading

0 comments on commit dfa6069

Please sign in to comment.