Design Process

The job of an L-System renderer is to take a description of an L-System (axiom, angle increment, and replacement rules) along with an iteration count, and produce a list of line segments to draw.

This project adapts my 2D L-system code from an earlier project to produce arrays of 3D coordinates, generating a 3D representation of these coordinates with PyCSG (a constructive solid geometry package for python), and saving these representations as STL files.

Approach
1. Expanding the L-system alphabet to 3 dimensions, and using trigonometry to create coordinates.
In 2-dimensions my code generated the L-system symbol string up to a given iteration, then interpreted the symbols as line segments and plotted them in matlab by specifying the start and end positions of each line segment.

Transitioning this process into three dimensions meant incorporating more symbols into the string, specifically to represent the angles of pitch and yaw, and doing a little extra trigonometry to produce the segment positions in 3D coordinates.


‍

‍2. Finding a CAD-compatible format to represent the shape in 3 dimensions.
After some research and experimentation with various 3d geometry packages, I settled on using the pyCSG package to generate three dimensional forms and converted these forms from vtk to stl files to create three-dimensional geometries compatible with my preferred CAD program Autodesk Fusion 360.


‍

PyCSG is a python implementation of Constructive Solid Geometry -- a modeling technique that uses a set of parameterized geometries and Boolean operations like union and intersection to combine these 3D forms. The python CSG library implements CSG operations on meshes concisely using BSP trees, and is an easily implemented version of the algorithm. 

Results

The renderer takes 3 arguments specifying the specific l-system pattern, number of iterations, and name of the new file to save the generated pattern.
‍
L-system patterns are stored in a dictionary and the rules are retrieved as needed:

    
      
     NAMED_LSYSTEMS = {
    'Test': LSystem(
        start = 'A',
        rules = dict(A='A&B', B='A^B'),
        turn_angle_deg1 = 45,
        turn_angle_deg2 = 45,),
        
    'sierpinski_arrowhead': LSystem(
        start = 'A',
        rules = dict(A='B-A-B', B='A+B+A'),
        turn_angle_deg1 = 60,
        turn_angle_deg2 = 0,),

    'dragon_curve': LSystem(
        start = 'FX',
        rules = dict(X='X+YF+', Y='-FX-Y'),
        turn_angle_deg1 = 90,
        turn_angle_deg2 = 0,),

    'barnsley_fern': LSystem(
        start = 'X',
        rules = dict(X='F+[[X]-X]-F[-FX]+X', F='FF'),
        turn_angle_deg1 = 25,
        turn_angle_deg2 = 0,)
    
			}
    
  

The core of this code is based on the two dimensional L-system renderer from my algorithmic design class; implementing a recursive approach to building these L-systems using two functions segments_recursive, and a helper function segments_recursive_helper. The former sets up some state variables and calls the latter, and the latter calls itself recursively.

    
      
    def segments_recursive(lsys, max_depth):

        cur_pos = np.array([0., 0., 0.])
        cur_angle_deg1 = 0
        cur_angle_deg2 = 0

        state_stack = [[cur_pos, cur_angle_deg1, cur_angle_deg2,]]


        segments = []
        s = lsys.start
        segments_recursive_helper(s, lsys.rules, lsys.turn_angle_deg1,lsys.turn_angle_deg2,
                                  max_depth, state_stack, segments)

        array = np.array(segments)
        return [np.array(segments) , array]
    
  

As can be seen above, the segment_recursive function acts to set up the l-systems initial conditions; it initializes the coordinates and angles, creates an empty stack to save states, and then hands off the initial axiom string and system rules to the helper function.

The helper function calls itself recursively when string replacement is required, producing an array of 3d coordinates.

    
      
    def segments_recursive_helper(s, rules, turn_angle_deg1,turn_angle_deg2,
                              remaining_steps,
                              state_stack,
                              segments):

        # for each symbol in input
        for symbol in s:

            # see if we can run a replacement rule for this symbol
            if remaining_steps > 0 and symbol in rules:

                # get the replacement
                replacement = rules[symbol]

                # recursively call this function with fewer remaining steps
                segments_recursive_helper(replacement, rules, 
                                          turn_angle_deg1, turn_angle_deg2,
                                          remaining_steps-1,
                                          state_stack, segments)

            else: # execute symbol directly

                cur_state = state_stack[-1]

                if symbol.isalpha(): 

                    # gives us the new x, y, z coords to add to current position:
                    cur_horiz_theta = cur_state[HORIZ_ANGLE] * np.pi / 180
                    cur_vert_theta = cur_state[VERT_ANGLE] * np.pi / 180
                    h = np.cos(cur_vert_theta)
                    offset = np.array([np.cos(cur_horiz_theta)*h, np.sin(cur_horiz_theta)*h, np.sin(cur_vert_theta) ])

                    #keep appending the old and new position to the segment list in pairs 
                    #-- each segment starts where the last one ended
                    new_pos = cur_state[POS_IDX] + offset
                    segments.append([cur_state[POS_IDX], new_pos])
                    cur_state[POS_IDX] = new_pos

                elif symbol == '+':

                    cur_state[HORIZ_ANGLE] += turn_angle_deg1

                #Additonal specifications for a 3d system -------------------

                elif symbol == '&': # pitch down by specified angle

                    cur_state[VERT_ANGLE] += turn_angle_deg2

                elif symbol == '^': # pitch up by specified angle

                    cur_state[VERT_ANGLE] += turn_angle_deg2


                elif symbol == '-':

                    cur_state[HORIZ_ANGLE] -= turn_angle_deg1

                elif symbol == '[':

                    state_stack.append([cur_state[0], cur_state[1], cur_state[2]])

                elif symbol == ']':

                    state_stack.pop()

                else:

                    raise RuntimeError('invalid symbol: ' + symbol)
    
  

    
      
    #Combine cylinders
    def combine_cyl(b, a):

        new_picture1 = a.union(b).toPolygons();
        #creates one combined polygon in the end 
        #all future operations have to be done on previous conglomerate, 
        #so the system keeps growing as one combined shape
        return new_picture1

    # Generating this whole thing in 3D 
    def generate_in_3D(array, segments):    
        print("original array:", array)

        arg1_a = array[0][0]
        arg1 = arg1_a.tolist()


        print ('arg1:', arg1)
        arg2_a = array[0][1]
        arg2 = arg2_a.tolist()

        print ('arg2:', arg2)
        radius = 1.0
        args = 'start =%s , end =%s , radius = %s' % (arg1, arg2, radius)

        cylinder_previous = CSG.cylinder(args)
        for i in (len(segments)-1):
            arg1_a = array[i+1][0]
            arg1 = arg1_a.tolist()
            arg2_a = array[i+1][1]
            arg2 = arg2_a.tolist()
            args = 'start =%s , end =%s , radius = %s' % (arg1, arg2, radius)

            cylinder_current = CSG.cylinder(args)



            #combine cylinders
            new_picture = combine_cyl(cylinder_current, cylinder_previous)


            cylinder_previous = cylinder_current

        return new_picture
    
  

    
      
def coordinates_to_vtk(new_array_3D, name):
     new_array_3D.saveVTK('%s.vtk'%(name))
     print('file %s .vtk written'%(name))
    
def vtk_to_stl(new_file_vtk, name):
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(new_file_vtk)

    surface_filter = vtk.vtkDataSetSurfaceFilter()
    surface_filter.SetInputConnection(reader.GetOutputPort())

    triangle_filter = vtk.vtkTriangleFilter()
    triangle_filter.SetInputConnection(surface_filter.GetOutputPort())

    writer = vtk.vtkSTLWriter()
   
    writer.SetFileName(name)
    writer.SetInputConnection(triangle_filter.GetOutputPort())
    writer.SetFileTypeToBinary()
    writer.Write()
    print ('file %s .stl written'%(name))
    
  

When no replacement is found or the desired number of iterations has been reached, the array is returned and the task of producing the Constructive Solid Geometry representation of the system is delegated to the generate_in_3D function. Each coordinate in the L-system array is used to dictate the start or end of a 3d form.

In this case I chose to use cylinders to represent each segment of the L-system. The pyCSG package requires three arguments: the coordinates of the start and end of the cylinder, and its radius, which I kept constant.




The final 3D shape was built iteratively, with separate function combine_cyl which used the .union command to add each new segment to a growing conglomerate shape. Each time the function is called it adds the latest cylinder to the previously saved shape, and returns the new combined polygon in the end.

In the early testing phases of this project I found that using the Vision Toolkit (VTK) to view output 3D forms was a convenient intermediate step for reviewing the produced l-system forms. Once my code consistently produced 3D L-systems, I then implemented a method to convert these VTK files to STLs. For that reason my code produces a VTK file as well as an STL file which is the commonly used CAD-compatible file format.

EXAMPLE RENDERING: Modified-Fern

        -     Axiom: X
        -      Angle increment: 25°
        -     Rules:
                            -     X → X+[[F]-F]^X[-XF]&F
                            -     F → FF           
‍
        Rendering of first 3 iterations: