diff --git a/analyseGrammar.py b/analyseGrammar.py index 78e4803..f18ffd0 100644 --- a/analyseGrammar.py +++ b/analyseGrammar.py @@ -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) @@ -22,17 +23,14 @@ 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) @@ -40,69 +38,68 @@ def branching_turtle_to_coords(turtle_program, d0, theta=20., phi=20.): 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. @@ -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. @@ -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. @@ -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. @@ -244,5 +236,4 @@ def posneg(value: float) -> float: if value >= 0.: return 1. else: - return -1. - + return -1. \ No newline at end of file diff --git a/libGenerator.py b/libGenerator.py index 1c8dc9f..d548014 100755 --- a/libGenerator.py +++ b/libGenerator.py @@ -1,6 +1,5 @@ -import ast -import random import numpy as np +import random from numpy import math as npm default={"k": 3, @@ -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'] @@ -49,43 +44,37 @@ 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) @@ -93,15 +82,14 @@ def calBifurcation(d0): 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): @@ -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) diff --git a/main.py b/main.py index 4970434..6d76282 100755 --- a/main.py +++ b/main.py @@ -1,67 +1,58 @@ -# This script generates L-System networks and saves them as 3D TIFF images. - import random import numpy as np import cv2 from skimage import io - from preprocessing import resize_stacks, resize_volume from libGenerator import setProperties -from vSystem import F, I +from vSystem import F, I, R, B from analyseGrammar import branching_turtle_to_coords from visuals import plot_coords, print_coords from computeVoxel import process_network from utils import bezier_interpolation - # Default Parameters -n=20 # Number of networks to be created -d0mean = 20.0 # Mean diameter of base of tree -d0std = 5.0 # Standard deviation of base diameter -tissueVolume = (512,512,140) # Specify number of pixels in the tissue volume -outpath = "/media/sweene01/SSD/More/" # Output path +n = 5 # No. of networks to be created +d0_mean = 20.0 # Diameter of base of tree +d0_std = 5.0 # Standard deviation of base diameter +tissue_volume = (512, 512, 140) # Specify number of pixels in the tissue volume +outpath = "/home/sweene01/Documents/Test/" # Output path + +for file in range(n): + # Lindenmayer System Parameters + properties = { + "k": 3, + "epsilon": random.uniform(4, 10), # Proportion between length & diameter + "randmarg": random.uniform(2, 4), # Randomness margin between length & diameter + "sigma": 5, # Determines type deviation for Gaussian distributions + "stochparams": True, # Whether the generated parameters will also be stochastic + } + + d0 = np.random.normal(d0_mean, d0_std) # Randomly assign base diameter (no dimension) + niter = random.randint(6, 14) # Randomly assign number of V-System recursive loops + setProperties(properties) # Setting V-Sytem properties + print(f"Creating image ... {file} with {niter} iterations") + + # Run V-System grammar for n iterations + turtle_program = F(niter, d0) + + # Convert grammar into coordinates + coords = branching_turtle_to_coords(turtle_program, d0) + + # Analyse / sort coordinate data + update = bezier_interpolation(coords) + + # If you fancy, plot a 2D image of the network! + # plot_coords(newdata, array=True, bare_plot=False) # bare_plot removes the axis labels -# Lindenmayer System Parameters -properties = {"k": 3, - "epsilon": 7,#random.uniform(4,10), # Proportion between length & diameter - "randmarg": 3, # Randomness margin between length & diameter - "sigma": 5, # Determines type deviation for Gaussian distributions - "stochparams": True} # Whether the generated parameters will also be stochastic + # Run 3D voxel traversal to generate binary mask of V-System network + image = process_network(update, tVol=tissue_volume) + # Convert to 8-bit format + image = (255 * image).astype("int8") -def main(): - """Main function that generates L-System networks and saves them as 3D TIFF images.""" - for file in range(n): - # Randomly assign base diameter (no dimension) - d0 = np.random.normal(d0mean, d0std) - # Randomly assign number of L-System recursive loops - niter = random.randint(6,14) - # Setting L-Sytem properties - setProperties(properties) - - print('Creating image ... %i with %i iterations' %(file, niter)) - - # Run L-System grammar for n iterations - turtle_program = F(niter,d0) - - # Convert grammar into coordinates - coords = branching_turtle_to_coords(turtle_program,d0) - - # Analyse / sort coordinate data - update = bezier_interpolation(coords) - - # If you fancy, plot a 2D image of the network! - # plot_coords(newdata, array=True, bare_plot=False) # bare_plot removes the axis labels - - # Run 3D voxel traversal to generate binary mask of L-System network - image = process_network(update, tVol=tissueVolume) - - # Convert to 8-bit format - image = (255*image).astype('int8') - - # Save image volume - io.imsave(outpath+"Lnet_i{}_{}.tiff".format(niter,file), np.transpose(image, (2, 0, 1)), bigtiff=False) - - -if __name__ == "__main__": - main() + # Save image volume + io.imsave( + outpath + "Lnet_i{}_{}.tiff".format(niter, file), + np.transpose(image, (2, 0, 1)), + bigtiff=False, + ) diff --git a/utils.py b/utils.py index 4058d07..4d7c1b8 100644 --- a/utils.py +++ b/utils.py @@ -1,7 +1,9 @@ +import itk import numpy as np import bezier as bz import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D +from numpy import math as npm def yaw(angle_degrees): """ @@ -58,7 +60,7 @@ def roll(angle): return A -def rotate(yaw_angle=0.0, pitch_angle=0.0, roll_angle=0.0, vector=None): +def rotate(pitch_angle=0.0, roll_angle=0.0, yaw_angle=0.0, vector=None): """ Applies a sequence of yaw, pitch, and roll rotations to a vector. @@ -118,8 +120,8 @@ def magnitude(x): """ return np.sqrt(x[0]**2 + x[1]**2 + x[2]**2) - def bezier_interpolation(coords=None): + """ Interpolates a series of Bezier curves defined by input coordinates. @@ -136,14 +138,21 @@ def bezier_interpolation(coords=None): - ValueError: If the input list is empty. """ - # Error checks - if type(coords) != list: - raise TypeError("Input must be a list.") - if len(coords) == 0: - raise ValueError("Input list is empty.") - # Store coordinates - X, Y, Z, _, _, diam, lseg, _, _, _ = zip(*coords) + X = [] + Y = [] + Z = [] + diam = [] + + try: + for c in coords: + X.append(float((c[0]))) + Y.append(float((c[1]))) + Z.append(float((c[2]))) + diam.append(float((c[5]))) + except: + pass + nodes = np.vstack((X, Y, Z, diam)) # Define the Bezier curve @@ -165,7 +174,7 @@ def bezier_interpolation(coords=None): if j == 0: theChosenOne = vessel[:,j:(j+6)] else: - theChosenOne = vessel[:,j:(j+6)]wo + theChosenOne = vessel[:,j:(j+6)] points = theChosenOne # points = interpolate_this(theChosenOne) newNodes = generate_points(newNodes,points) diff --git a/vSystem.py b/vSystem.py index ff3a6a3..a2ac7f0 100755 --- a/vSystem.py +++ b/vSystem.py @@ -1,6 +1,7 @@ import random import numpy as np import random +from math import floor from libGenerator import setProperties, calBifurcation, calParam from analyseGrammar import posneg @@ -17,19 +18,12 @@ def I(n, d0, val=3): str: the resulting string representing the fractal pattern. ''' if n > 0: - # Calculate parameters for bifurcation and string formation params = calBifurcation(d0) - p1 = calParam(str.join('co/', str(int(val))), params) - - # Add a random rotation to the new line segment + p1 = calParam(str.join('co/',str(int(val))), params) rotate = np.random.uniform(22.5, 27.5) - - # Recursive call to generate the string for the next level return 'f(' + p1 + ',' + str(params['d0']) + ')' + '+(' + str(rotate) + ')' +\ - '[' + I(n-1, params['d0']) + ']' - else: - # Base case: return "I" - return 'I' + '[' + R(n-1, params['d0']) + ']' + else: return 'I' def R(n, d0): @@ -43,21 +37,15 @@ def R(n, d0): Returns: - A string description of a bifurcating tree. ''' - if n == 0: - return 'R' - - params = calBifurcation(d0) # Calculates bifurcation parameters based on parent diameter - p1 = calParam('co/' + str(int(3)), params) # Calculates the parameters for the first child branch - p2 = calParam('co/' + str(int(2)), params) # Calculates the parameters for the second child branch - - # Calls G function recursively to generate descriptions of child branches - # The descriptions are concatenated with function calls and brackets to form the final description - descrip = 'f(' + p1 + ')' + G(n, d0, val=7) + G(n-1, d0, val=7) + G(n-1, d0, val=7) \ - + '[' + B(n-1, params['d1']) + ']' \ - + 'f(' + p2 + ',' + str(params['d2']) + ')' + B(n-1, params['d2']) - - return descrip - + if n > 0: + params = calBifurcation(d0) + p1 = calParam(str.join('co/',str(int(3))), params) + p2 = calParam(str.join('co/',str(int(2))), params) + descrip = 'f(' + p1 + ')' + D(n-1, d0) +\ + D(n-1, d0) + D(n-1, d0) + '[' + B(n-1, params['d1']) +\ + ']' + 'f(' + p2 + ',' + str(params['d2']) + ')' + B(n-1, params['d2']) + return descrip + else: return 'R' def B(n, d0): @@ -74,15 +62,11 @@ def B(n, d0): Raises: - N/A """ - if n > 0: - # Recursively generate descriptions of child branches - descrip = G(n-1, d0, val=7) + G(n-1, d0, val=7) + G(n-1, d0, val=7) +\ - '/(' +str(90.0) + ')' + A(n, d0) - return descrip - else: - # Return string representation of parent branch - return 'B' + rotate = np.random.uniform(22.5, 27.5) + return D(n-1, d0) + D(n-1, d0) + D(n-1, d0) +\ + '/(' +str(rotate) + ')' + A(n-1, d0) + else: return 'B' def F(n, d0): @@ -97,26 +81,20 @@ def F(n, d0): - str: a string of characters representing the fractal generated by the given parameters. """ if n > 0: - # calculate the bifurcation parameters params = calBifurcation(d0) - - # calculate the angles for rotation and tilt + randInt = random.randint(-1,1) theta1 = params['th1'] theta2 = params['th2'] - tilt = np.random.uniform(22.5, 27.5) * random.randint(-1, 1) - - # recursively generate the fractal string using the calculated parameters - return S(n-1, d0) + '[' + '+(' + str(theta1) + ')/(' + str(tilt) + ')' + \ - F(n-1, params['d1']) + ']' + '[' + '-(' + str(theta2) + ')/(' + str(tilt) + ')' + \ - F(n-1, params['d2']) + ']' - else: - # return the base character if n <= 0 - return 'F' + tilt = np.random.uniform(22.5, 27.5)*random.randint(-1,1) + return S(n-1, d0)+'['+'+('+str(theta1)+')'+'/('+str(np.random.uniform(22.5, 27.5)*random.randint(-1,1))+')'+\ + F(n-1, params['d1'])+']'+'['+'-('+str(theta2)+')'+'/('+str(np.random.uniform(22.5, 27.5)*random.randint(-1,1))+')'+\ + F(n-1, params['d2'])+']' + else: return 'F' def S(n, d0, val=5, margin=0.5): """ - S function generates the shape of a tree's stem. + S function generates the shape of a tree branch. Args: n (int): The depth of the recursive algorithm. @@ -127,23 +105,15 @@ def S(n, d0, val=5, margin=0.5): Returns: str: A string representing the shape of a tree's stem. """ - # Check if the depth of the recursion is greater than zero - if n <= 0: - return 'S' - - # Generate a random number between 0 and 1 r = random.random() - - # Determine which sub-function to use based on the random number and margin value - if r >= 0.0 and r < margin: - return '{' + S1(n, d0, val) + '}' - if r >= margin and r < 1.0: - return '{' + S2(n, d0, val) + '}' + if r>= 0.0 and r < margin: return '{' + S1(n, d0, val) + '}' + if r >= margin and r < 1.0: return '{' + S2(n, d0, val) + '}' def S1(n, d0, val=5): """ Recursive function that generates a string representation of a tree structure. + Includes shrinkage and expansion of vessel diameters. Args: - n (int): the recursion level, i.e. the number of times the function is called recursively. @@ -154,43 +124,62 @@ def S1(n, d0, val=5): - str: string representation of a tree structure. """ if n > 0: - # Randomly rotate the branches of the tree - rotate = np.random.uniform(22.5, 27.5) * random.randint(-1, 1) - params = calBifurcation(d0) - descrip = D(n-1, params['d0'], val) + '+(' + str(rotate) + ')' + \ - D(n-1, params['d0'], val) + '-(' + str(rotate) + ')' + \ - D(n-1, params['d0'], val) + '-(' + str(rotate) + ')' + \ - D(n-1, params['d0'], val) + '+(' + str(rotate)+ ')' + \ - D(n-1, params['d0'], val) + # "Fanning" of trees + # rotate = np.random.uniform(22.5, 27.5)*random.randint(-1,1) + params = calBifurcation(d0) # calculates the branching parameters based on the starting thickness of the trunk + randInt = random.randint(-1, 1) # randomly generates an integer of either -1, 0, or 1 + theta1 = params['th1'] * randInt # calculates the angle of the first branch based on the random integer + theta2 = params['th2'] * abs(randInt) # calculates the angle of the second branch based on the absolute value of the random integer + if random.random() < 0.1: # randomly generates growth and decay parameters + growth = np.sort(np.random.uniform(1., 1.5, 3)) # generates 3 growth factors between 1 and 1.5 and sorts them in ascending order + decay = -np.sort(-np.random.uniform(1., 1.5, 2)) # generates 2 decay factors between 1 and 1.5, negates them, and sorts them in ascending order + else: + growth = np.ones(3) # sets growth factors to 1 if the random number is greater than or equal to 0.1 + decay = np.ones(2) # sets decay factors to 1 if the random number is greater than or equal to 0.1 + # recursively calls the D() function to generate descriptions for each sub-tree and concatenates them into a single string + descrip = D(n-1, growth[0]*params['d0'], val) + '+(' + str(theta1) + ')' + \ + D(n-1, growth[1]*params['d0'], val) + '-(' + str(theta2) + ')' + \ + D(n-1, growth[2]*params['d0'], val) + '-(' + str(theta1) + ')' + \ + D(n-1, decay[0]*params['d0'], val) + '+(' + str(theta2) + ')' + \ + D(n-1, decay[1]*params['d0'], val) return descrip else: - return 'S' + return 'S' # returns the string 'S' if n is 0 def S2(n, d0, val=5): """ Recursive function that generates a string representation of a tree structure. - - Parameters: - n (int): the level of recursion. Controls the depth of the tree-like structure. - d0 (float): the trunk diameter of the tree at the base. - val (int): the numerical value to assign to the final string. Default value is 5. - + Includes shrinkage and expansion of vessel diameters. + + Args: + - n (int): the recursion level, i.e. the number of times the function is called recursively. + - d0 (float): the initial diameter of the tree. + - val (int): value to be passed to function D(). Default value is 5. + Returns: - descrip (str): the final string representation of the tree-like structure. + - str: string representation of a tree structure. """ - if n > 0: + if n>0: # "Fanning" of trees - rotate = np.random.uniform(22.5, 27.5) * random.randint(-1, 1) + #rotate = np.random.uniform(22.5, 27.5)*random.randint(-1,1) params = calBifurcation(d0) - descrip = D(n-1, params['d0'], val) + '-(' + str(rotate) + ')' + \ - D(n-1, params['d0'], val) + '+(' + str(rotate) + ')' + \ - D(n-1, params['d0'], val) + '+(' + str(rotate) + ')' + \ - D(n-1, params['d0'], val) + '-(' + str(rotate) + ')' + \ - D(n-1, params['d0'], val) + randInt = random.randint(-1,1) + theta1 = params['th1']*randInt #+ np.random.uniform(-2.5, 2.5) + theta2 = params['th2']*abs(randInt) #+ np.random.uniform(-2.5, 2.5) + if random.random() < 0.1: + growth = np.sort(np.random.uniform(1., 1.5, 3)) + decay = -np.sort(-np.random.uniform(1., 1.5, 2)) + else: + growth = np.ones(3) + decay = np.ones(2) + descrip = D(n-1, growth[0]*params['d0'], val) + '-(' + str(theta1) + ')' + \ + D(n-1, growth[1]*params['d0'], val) + '+(' + str(theta2) + ')' + \ + D(n-1, growth[2]*params['d0'], val) + '+(' + str(theta1) + ')' + \ + D(n-1, decay[0]*params['d0'], val) + '-(' + str(theta2) + ')' + \ + D(n-1, decay[1]*params['d0'], val) return descrip - else: - return 'S' + else: return 'S' def D(n, d0, val=5): @@ -207,21 +196,14 @@ def D(n, d0, val=5): str: The generated L-system symbols string for the fractal tree. """ - if n > 0: - # Calculate the parameters for the bifurcation of the tree + if n>0: params = calBifurcation(d0) - - # Calculate the parameter for the function call - p1 = calParam(str.join('co/', str(int(val))), params) - - # Return the L-system string with the calculated parameters + p1 = calParam(str.join('co/',str(int(val))), params) return 'f(' + p1 + ',' + str(params['d0']) + ')' - else: - # If the number of iterations is 0, return the symbol 'D' - return 'D' + else: return 'D' -def G(n, d0, val=5, shift=18.0): +def G(n, d0, val=5): """ Generates a string representing the result of the function f with parameters p1 and d0. @@ -234,13 +216,11 @@ def G(n, d0, val=5, shift=18.0): Returns: - str: a string representing the result of the function f with parameters p1 and d0. """ - if n > 0: + if n>0: params = calBifurcation(d0) - p1 = calParam(str.join('co/', str(int(val))), params) - # Call function f with parameters p1 and d0, and add the value of shift to the result - return str(float(f(p1, params['d0'])) + shift) - else: - return 'G' + p1 = calParam(str.join('co/',str(int(val))), params) + return 'f(' + p1 + ',' + str(params['d0']) + ')' + else: return 'G' def A(n, d0): @@ -253,24 +233,18 @@ def A(n, d0): str: String description of the fractal tree """ if n > 0: - # calculate bifurcation angles and branch lengths params = calBifurcation(d0) - # recursive calls to A and S functions - return ( - S(n-1, d0) - + '[+(' + str(params['th1']) + ')' + A(n-1, params['d1']) + ']' - + '[+(' + str(params['th2']) + ')' + A(n-1, params['d2']) + ']' - ) - else: - return 'A' + return S(n-1, d0) + '[+(' + str(params['th1']) + ')' + A(n-1, params['d1']) + ']' +\ + '[+(' + str(params['th2']) + ')' + A(n-1, params['d2']) + else: return 'A' + + + -# def E(d0): -# return 'f(' + str(1) + ',' + str(d0) + ')' + 'f(' + str(1) + ',' + str(0.5*d0) + ')' + \ -# 'f(' + str(1) + ',' + str(d0) + ')' # 2D grammars - don't currently work with framework def simplest_gramma(n=None, theta1=20., theta2=20., params=None): return 'f' + '[' + '+' + F(n-1, params['d0']) + ']' + '-' + F(n-1, params['d0']) def simple_grammar(n=None, theta1=20., theta2=20., params=None): - return 'f'+'['+'+('+str(theta1)+')'+F(n-1,params['d1'])+']'+'['+'-('+str(theta2)+')'+F(n-1,params['d2'])+']' + return 'f'+'['+'+('+str(theta1)+')'+F(n-1,params['d1'])+']'+'['+'-('+str(theta2)+')'+F(n-1,params['d2'])+']' \ No newline at end of file