Need help, Python recursive subdivision

Grasshopper Developer
1

Hey,
I’ve been trying to create a python function that takes a flat list of closed triangular polylines and subdivide dem at a given parameter (in my example .5 of each curve) returning 4 new triangles for each polyline in the line. I would then like to set the number of iterations, so that 1 polyline triangle first turns into a total of 4, 16, 64, etc, forming smaller and smaller polyline triangles within the scope of the first one. I’m not a very experienced progammer and this is the first time I take on the task of creating a recursive function, so please bear with me.

I’ve managed to get the basic recursiveness of the function working, altough it the recursive part seems to only iterate over the last polyline in the list. I’ve been staring at this for hours now and would be really greatful for any help getting this right!

Here is my function:

#performs subdivision on a list of equilateral triangle polylines

def subDivEquiTriangles2(polylines, generation, recursionLevel):
    if generation > recursionLevel: return
    generation += 1

    for p in polylines:
        newPoints = []
        for i in range(len(p)-1):
            newPoints.append(Rhino.Geometry.Point3d((p[i][0] + p[i+1][0])/2, (p[i][1] + p[i+1][1])/2, (p[i][2] + p[i+1][2])/2))
        subD = [
        Rhino.Geometry.Polyline([p[0], newPoints[0], newPoints[2], p[0]]),
        Rhino.Geometry.Polyline([newPoints[2], newPoints[1], p[2], newPoints[2]]),
        Rhino.Geometry.Polyline([newPoints[0], p[1], newPoints[1], newPoints[0]]), 
        Rhino.Geometry.Polyline([newPoints[0], newPoints[1], newPoints[2], newPoints[0]])
        ]
        #REPLACE EACH VALUE IN SUBD WITH 4 NEW
        if generation <= recursionLevel:
            for i in range(len(subD)):
                subD[i:i+1] = subDivEquiTriangles2(subD, generation, recursionLevel)
    return subD

and here is a picture of what it produces:

Namnlöst-1.jpg813×610 51.1 KB

Many thanks

Jakob

2

There is a very nice post on recursion in Python on Stackoverflow.

1 Like
3

Hi Jakob

… Not sure about what’s wrong in the script,
But I would not modify a list while iterating over it (latest loop).
I think that using two lists there might be simpler (just my humble opinion :slight_smile: )

Anyway the following test Rhino script seems to work here …

HTH regards.

import Rhino
import scriptcontext

def one2four( pli ):
  p01 = ( pli[ 0 ] + pli[ 1 ] ) / 2.0
  p12 = ( pli[ 1 ] + pli[ 2 ] ) / 2.0
  p20 = ( pli[ 2 ] + pli[ 0 ] ) / 2.0
  return [
      Rhino.Geometry.Polyline( [ pli[ 0 ], p01, p20, pli[ 0 ] ] ),
      Rhino.Geometry.Polyline( [ pli[ 1 ], p12, p01, pli[ 1 ] ] ),
      Rhino.Geometry.Polyline( [ pli[ 2 ], p20, p12, pli[ 2 ] ] ),
      Rhino.Geometry.Polyline( [ p01, p12, p20, p01 ] ) ]

def dividelist( plis, times ):
  if times < 1:
    return plis
  result = []
  for pli in plis:
    result.extend( one2four( pli ) )
  if times > 1:
    result = dividelist( result, times - 1 )
  return result

def main():
  pli = Rhino.Geometry.Polyline( [
      Rhino.Geometry.Point3d( 0, 0, 0 ),
      Rhino.Geometry.Point3d( 100, 0, 0 ),
      Rhino.Geometry.Point3d( 50, 87, 0 ),
      Rhino.Geometry.Point3d( 0, 0, 0 ) ] )
  cnt = 1
  res, cnt = Rhino.Input.RhinoGet.GetInteger( 'Recursion depth ?', True, cnt )
  if res == Rhino.Commands.Result.Success:
    plis = dividelist( [ pli ], cnt )
    for pli in plis:
      scriptcontext.doc.Objects.Addpolyline( pli )
    scriptcontext.doc.Views.Redraw()

main()
1 Like
4

@MarcusStrube
Amazing tip, many thanks!

@emilio
Super helpful, thanks!
I see now that reading and writing to the same list might be a bad idea. I gave your example a go and now it works pretty well.

The code still doesn’t return a perfect icosasphere with only equilaterial triangles, but its pretty close, especially for the rougher approximations. I have no idea how hard or easy that will be to achieve, but I’ll give it a try later.

In the meantime, here is the full python script if you feel like giving it a go :wink:
It takes an integer value for the level of recursion, and returns an icosasphere subdivided by that number:

import math
import Rhino

#Define function to shift a list in place
def shiftInPlace(l, n):
    n = n % len(l)
    head =l[:n]
    l[:n] = []
    l.extend(head)
    return l

#Define a function to move points to the unit sphere
#Remove this bit if subdivision is beeing performed on another set of triangles than a icosasphere
def point2unitsphere(Point):
    l = math.sqrt(Point[0] * Point[0] + Point[1] * Point[1] + Point[2] * Point[2])
    return Rhino.Geometry.Point3d(Point[0]/l, Point[1]/l, Point[2]/l)

#Define function to draw the points of an icosahedron
def icosahedron():
    t = ( 1 + math.sqrt(5)) / 2
    points = []
    polylines = []
    for k in reversed(range(-2, 1)):
        c = []
        for i in reversed(range(-1, 2, 2)):
            for j in range(-1, 2, 2):
                c.append([j, t * i, 0])
        for j in c:
            points.append(point2unitsphere(shiftInPlace(j, k)))

    polylines.append(Rhino.Geometry.Polyline([points[0], points[11], points[5], points[0]]))
    polylines.append(Rhino.Geometry.Polyline([points[0], points[5], points[1], points[0]]))
    polylines.append(Rhino.Geometry.Polyline([points[0], points[1], points[7], points[0]]))
    polylines.append(Rhino.Geometry.Polyline([points[0], points[7], points[10], points[0]]))
    polylines.append(Rhino.Geometry.Polyline([points[0], points[10], points[11], points[0]]))

    polylines.append(Rhino.Geometry.Polyline([points[1], points[5], points[9], points[1]]))
    polylines.append(Rhino.Geometry.Polyline([points[5], points[11], points[4], points[5]]))
    polylines.append(Rhino.Geometry.Polyline([points[11], points[10], points[2], points[11]]))
    polylines.append(Rhino.Geometry.Polyline([points[10], points[7], points[6], points[10]]))
    polylines.append(Rhino.Geometry.Polyline([points[7], points[1], points[8], points[7]]))

    polylines.append(Rhino.Geometry.Polyline([points[3], points[9], points[4], points[3]]))
    polylines.append(Rhino.Geometry.Polyline([points[3], points[4], points[2], points[3]]))
    polylines.append(Rhino.Geometry.Polyline([points[3], points[2], points[6], points[3]]))
    polylines.append(Rhino.Geometry.Polyline([points[3], points[6], points[8], points[3]]))
    polylines.append(Rhino.Geometry.Polyline([points[3], points[8], points[9], points[3]]))

    polylines.append(Rhino.Geometry.Polyline([points[4], points[9], points[5], points[4]]))
    polylines.append(Rhino.Geometry.Polyline([points[2], points[4], points[11], points[2]]))
    polylines.append(Rhino.Geometry.Polyline([points[6], points[2], points[10], points[6]]))
    polylines.append(Rhino.Geometry.Polyline([points[8], points[6], points[7], points[8]]))
    polylines.append(Rhino.Geometry.Polyline([points[9], points[8], points[1], points[9]]))

    return polylines

#subdivides a single triangular polyline into 4 smaller
def subdivide_triangle(polylines, recursionLevel):

    #Define function to subdive a single closed triangular polyline to 4 smaller ones
    def one2four(polyline):
        p01 = point2unitsphere((polyline[0] + polyline[1]) / 2)
        p12 = point2unitsphere((polyline[1] + polyline[2]) / 2)
        p20 = point2unitsphere((polyline[2] + polyline[0]) / 2)

        return [
        Rhino.Geometry.Polyline([polyline[0], p01, p20, polyline[0]]),
        Rhino.Geometry.Polyline([p20, p12, polyline[2], p20]),
        Rhino.Geometry.Polyline([p01, polyline[1], p12, p01]), 
        Rhino.Geometry.Polyline([p01, p12, p20, p01])
        ]

    def dividelist(polylines, recursionLevel):
        if recursionLevel < 1: return polylines
        result = []
        for polyline in polylines:
            result.extend(one2four(polyline))
        if recursionLevel > 1:
            result = dividelist(result, recursionLevel - 1)
        return result

    polylines = dividelist(polylines, recursionLevel)
    return polylines

Icosahedron = icosahedron()
a = subdivide_triangle(Icosahedron, recursionLevel)
5

I think only for 4, 8 or 20 equilateral triangles, you keep the vertices on a sphere.
Otherwise the triangles will be somehow deformed.

Thanks for the script !

Regards