Edge division problem in Growth Algorithm

Dear all,

I’m having difficulties implementing with “toxiclibs” a sketch initially relying on the “iGeo” library.

The original script is fairly long so I tried to boil it down to its simplest form but its relative complexity requires that I give some context before going into the details of my problem.

The algorithm goes like this:

  • connect 3 cells to form a triangle
  • select one of the edges of that triangle
  • divide it into 2 and insert a new cell there
  • create a new link from that cell to its opposite vertex (forming 2 inner triangles)
  • repeat from step 2

The cells are connected by springs and their size increased at each iteration, causing the overall shape to expand and potentially allowing the endless formation of new “inner” triangles.

The output should look like this (picture taken from the iGeo website):

However in my case:

  • triangles are overlapping instead of standing side by side
  • springs seem to appear from the center of the canvas before rapidly heading back to their set location

That last point is, in my opinion, what causes the overlapping but I’m unable to understand the reason behind this behavior.

Because the algorithm is quite complex (not complicated, but fairly long) I tried to sum up on a slide the edge division mechanism:

  • links (springs) between p1, p2 and p3 form a triangle
  • a cell is selected (here p1)
  • one of its neighbor is picked (here p2)
  • a “child” cell will be placed on the edge connecting p1 and p2 (here p4)
  • its exact location is determined by the radius of its “mother” cell
  • when finally placing p4, the edge p1-p2 is removed
  • 3 new edges (or links/springs) are created instead: p1-p4, p2-p4, p3-p4

But it is at this final step that for some reason springs appear at the center of the canvas instead of at the cells’ locations.


  • What causes that behavior ?
  • Is there something else that I could be doing wrong ?
from toxi.physics3d.behaviors import GravityBehavior3D, AttractionBehavior3D
from toxi.physics3d import VerletPhysics3D, VerletParticle3D, VerletSpring3D

nCells = 240
growthInterval = 10
divisionInterval = 50
repulsiveForce = .4
springForce = .005

def setup():
    size(700, 400, P3D)
    global cam, physics, cells
    # Array containing the "Cell" objects
    cells = [Cell(0, 0, 0)]
    cam = PeasyCam(this, 80)
    physics = VerletPhysics3D()
    physics.addBehavior(AttractionBehavior3D(cells[0], cells[0].radius*2, -repulsiveForce))
def draw():
    # Updating physics engine
    for c in cells:
    # Adjusting spring length (sum of the radii of the 2 cells forming a link)
    for i in xrange(len(physics.springs)):
        physics.springs.get(i).setRestLength(physics.springs.get(i).a.radius + physics.springs.get(i).b.radius)
        # Drawing springs
        line(physics.springs.get(i).a.x(), physics.springs.get(i).a.y(), physics.springs.get(i).a.z(), physics.springs.get(i).b.x(), physics.springs.get(i).b.y(), physics.springs.get(i).b.z())
class Cell(VerletParticle3D):    
    def __init__(self, x, y, z):
        super(Cell, self).__init__(x, y, z)
        self.location = Vec3D(x, y, z)
        self.radius = 5
        self.maxRadius = 20
        self.id = 0
        self.time = 0
        self.links = []
        self.faces = []
    def grow(self):
        # Incrementing cell's lifetime
        self.time += 1
        # As long as we haven't reached the max number of cells yet --> keep growing
        if len(cells) < nCells:
            if random(1) > .5:
                if self.time > 0 and self.time%divisionInterval == 0:
                    # Point state: forming a line
                    if not self.links:
                        # downsizing mother cell when giving birth to child cell
                        self.radius *= .5
                        # location of child cell
                        newloc = self.location.add(Vec3D.randomVector().scale(self.radius))
                        child = Cell(newloc.x(), newloc.y(), self.z())
                        child.id = len(cells)
                        # creating a link between mother cell and child cell
                        self.connect(self, child, False)
                    # Line state: forming a triangle
                    elif len(self.links) == 1:
                        # downsizing mother cell when giving birth to child cell
                        self.radius *= .5
                        # location of child cell
                        newloc = self.location.add(Vec3D.randomVector().scale(self.radius))
                        child = Cell(newloc.x(), newloc.y(), self.z())
                        child.id = len(cells)
                        # creating a link between mother cell and child cell
                        # creating a link between child cell and cell on the other side of the link
                        self.connect(self, child, False)
                        self.connect(child, cells[self.links[0]], True)
                        # Storing face (triangle) information in each cell
                        for c in (self, child, cells[self.links[0]]):
                            c.faces.append({self.id, child.id, cells[self.links[0]].id})
                    # Triangle state: edge division (core of the algorithm)
                        # Picking a neighboring cell (oppcell1)
                        dlink = self.links[int(random(len(self.links)-1))]
                        oppcell1 = cells[dlink]
                        # Not selecting cells belonging to 2 faces (triangles)
                        if len(oppcell1.faces) == 2:
                            # downsizing mother cell when giving birth to child cell
                            self.radius *= .5
                            # Computing child cell location
                            dir = oppcell1.sub(self.location).normalize().scale(self.radius)
                            newloc = self.location.add(dir)
                            child = Cell(newloc.x(), newloc.y(), newloc.z())
                            child.id = len(cells)
                            # Finding the opposite cell (the other neighboring cell at the top of the triangle)
                            face = [f for f in oppcell1.faces if self.id in f and oppcell1.id in f][0]
                            oppcell2 = cells[[id for id in face if id != dlink and id != self.id][0]]
                            # Remove links (info + corresponding spring) between self and its selected neighboring cell 
                            physics.removeSpring(physics.getSpring(self, oppcell1))
                            # Creating new links for replacement
                            self.connect(self, child, False)
                            self.connect(oppcell1, child, True)
                            self.connect(oppcell2, child, True)
                            # Removing former face (triangle) information from each cell
                            for c in (self, oppcell1, oppcell2):
                            # Adding new face (inner triangle 1) information
                            for c in (self, child, oppcell2):
                                c.faces.append({self.id, child.id, oppcell2.id})
                            # Adding new face (inner triangle 2) information   
                            for c in (oppcell1, child, oppcell2):
                                c.faces.append({oppcell1.id, child.id, oppcell2.id})

                # Increasing cell size intermittently
                if self.time%growthInterval == 0:
                    if self.radius < self.maxRadius:
                        self.radius += .1
    def connect(self, c1, c2, added):
        if not added: # means that the child cell has not been included yet in the physics engine
            physics.addBehavior(AttractionBehavior3D(c2, c2.radius*2, -repulsiveForce))
        #Add spring between c1 and c2
        s = VerletSpring3D(c1, c2, c2.distanceTo(c1), springForce) 
        #Storing linkage information for each cell
    def render(self):
        point(self.x(), self.y(), self.z())
1 Like

It took me two days to find the problem, turns out I just had to replace self.location and oppcell1.location by self and oppcell1.
The devil is in the detail, as they say.