Problem with spring force in growth algorithm

Dear all,

I would need a hand to help me find what is wrong with the following implementation of a sketch originally relying on the iGeo library (not compatible with Processing 3).

Simply put, the algorithm consists in iteratively adding a random number of faces to a tetrahderon. The mesh gets bigger over time and the overall shape should gradually start to look like a giant 3D blob with bumps or protrusions:

However in my case, newly added faces tend to rapidly stand out from the mesh, forming sharp triangular spikes all over the shape.

Capture%20d%E2%80%99%C3%A9cran%20(8)

see GIF

In my opinion, 2 things might cause this behavior:

  • the lack of a force continuously pressuring the vertices outward (unlike my sketch, the original script uses the push() method from the iGeo library to pressure the vertices towards their corresponding normal)

  • a difference between the physics of the VerletSprings class (my sketch is using the toxiclibs library) and the physics of the spring forces in the original script.

That second assumption seemed unlikely at first but I noticed that commenting out the push() method in the original script didn’t prevent the mesh to grow “normally” (the output looked a bit smaller in size but similar in shape, no spikes)

In order to be sure I would like to implement the snippet related to the spring force from the original script:

maxForce = 100

def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

The problem is that I’m unable to come up with a correct transcription of that code using toxiclibs:

def interact(self):
    maxForce = 100

    for l in self.links:
        dif = self.sub(cells[l])
        force = (dif.magnitude() - (self.radius + cells[l].radius)) / (self.cell1.radius + self.cell2.radius) * 300
        if force > maxForce: 
            force = maxForce
        elif force < -maxForce:
            force = -maxForce
        dif.normalizeTo(force)
        self.addForce(dif.scale(-1))
        cells[l].addForce(dif)

I believe I’m calling the correct methods:

  • iGeo dif() = toxiclibs sub()
  • iGeo len() = toxiclibs magnitude()
  • iGeo len(force) = toxiclibs normalizeTo(force) (similar to setMag(force) in Processing)
  • I’m assuming that “pushing” or “pulling” a cell (push() / pull() in iGeo) is like adding or subtracting a force to / from that cell, so I’m using toxiclibs addForce()

however the result is catastrophic, almost impossible to describe (cells are violently shaking, fleeing all over the place).

QUESTIONS:

  • Could you help me implement the interact function (spring force) correctly with toxiclibs ?
  • Do you think something else could be the cause of that strange behavior (spikes formation) ?

original sketch (Python version / Processing 2 only)
add_library('igeo')

mesh = None # global variable of mesh geometry

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    mesh = IMesh().clr(0.7,0,0)
    cell = Cell(IVec(0,0,0), 10)

class Cell(IParticle) : 
    growthDuration = 1200 #duration of growth and division
    growthInterval = 10
    divisionInterval = 50
    maxRadius = 100
    growthSpeed=0.1
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.vertex = IVertex(self.pos())
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(50) # constant force
            self.push(dif)

        self.push(self.nml().len(20)) # pressure toward normal
  
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
                if IRand.pct(50) : # random division
                    self.active = True
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(mesh.clr())
        self.vertex.pos().set(self.pos()) # update mesh vertex position
        self.vertex.nml(self.nml()) # update mesh vertex normal
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            # making a triangle loop
            child = self.createChild(IRand.dir())
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        elif len(self.links)==2 : # triangle
            # making a tetrahedron enclosure
            child = self.createChild(IRand.dir())
            f = self.faces[0]
            center = IVec.center(child.pos(),f.cell1.pos(),f.cell2.pos(),f.cell3.pos())
            if f.center().dif(center).dot(f.nml())<0 :
                f.flipNml() # adjust normal to be outward
            CellLink(f.cell1, child)
            CellLink(f.cell2, child)
            CellLink(f.cell3, child)
            f1 = CellFace(f.cell1, f.cell2, child)
            f2 = CellFace(f.cell2, f.cell3, child)
            f3 = CellFace(f.cell3, f.cell1, child)
            if f1.center().dif(center).dot(f1.nml())<0 : 
                f1.flipNml()
            if f2.center().dif(center).dot(f2.nml())<0 : 
                f2.flipNml()
            if f3.center().dif(center).dot(f3.nml())<0 : 
                f3.flipNml()
        else : # link num > 3
            # pick two links and insert two faces between
            linkIdx = self.findDividingLinks()
            linksToChild = [] # links to reconnect to child
            for i in range(linkIdx[0], linkIdx[1]+1) : 
                linksToChild.append(self.links[i%len(self.links)])
            dir = IVec()
            for l in linksToChild :  # average of links1 dir
                dir.add(l.oppositeDir(self).unit())
            dir.projectToPlane(self.nml())
            child = self.createChild(dir)
            self.insertChild(child, linksToChild)
    
    def nml(self) : # calc vertex normal from face normal
        if len(self.faces)==0 :
            return IVec(0,0,1)
        n = IVec(0,0,0)
        for f in self.faces : 
            n.add(f.nml())
        return n.unit()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        self.active = False #reset activation
        return child
    
    def sortLinks(self) : # sort links by following connected faces
        if len(self.links) <= 3 :
            return # no need to sort
        sorted = []
        currentLink = self.links[0]
        currentFace = currentLink.faces[0]
        for i in range(len(self.links)) : 
            sorted.append(currentLink)
            currentLink = currentFace.oppositeLink(currentLink.oppositeCell(self))
            currentFace = currentLink.oppositeFace(currentFace)
        self.links = sorted
    
    def findDividingLinks(self) :  # find splittable two links and return the index numbers
        linkInterval = len(self.links)//2 # interval between two links 
        
        self.sortLinks() # sort the order of links around cell
        idx = IRand.getInt(0,len(self.links)-1)
        return [ idx, idx+linkInterval ] # index can be larger than self.links.size()
    
    def insertChild(self, child, linksToChild) : # insert child cell and reconnect links
        removingFaces = []
        for f in self.faces : 
            for j in range(len(linksToChild)-1) : 
                if (f.contains(linksToChild[j]) and
                    f.contains(linksToChild[j+1])) : #a face between replacing links is removed
                    removingFaces.append(f)
                    break
        for i in range(1, len(linksToChild)-1) : 
            linksToChild[i].delete() # replacing links are once removed
            CellLink(linksToChild[i].oppositeCell(self), child) # then recreated
        cell1 = linksToChild[0].oppositeCell(self) # one of two dividing link cell
        cell2 = linksToChild[len(linksToChild)-1].oppositeCell(self) # another dividing link cell
        CellLink(child, cell1)
        CellLink(child, cell2)
        CellLink(self, child)
        for f in removingFaces : 
            f.delete() # replace face by deleting and recreating
            cells = f.oppositeCell(self)
            CellFace(child, cells[0], cells[1])
        CellFace(self, cell1, child)
        CellFace(self, child, cell2)

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.del() # stop agent
  
    def oppositeFace(self, f) : # find other face on the link
        if len(self.faces)!=2 :
            return None
        if self.faces[0] is f :
            return self.faces[1]
        return self.faces[0]
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.link1 = self.findLink(c1, c2)
        # keep order of cells in right hand screw for consistent normal
        if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
            self.cell1 = c2 
            self.cell2 = c1 
            self.cell3 = c3
        else : 
            self.cell1 = c1 
            self.cell2 = c2 
            self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
        self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
        mesh.addFace(self.face) # add triangular face to the mesh
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def nml(self) : # normal vector
        return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
    
    def flipNml(self) : # flip normal
        tmp = self.cell2
        tmpLink = self.link2
        self.cell2 = self.cell3
        self.link2 = self.link3
        self.cell3 = tmp
        self.link3 = tmpLink

    # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
    def cellOrder(self, c1, c2) : 
        if ( self.cell1 is c1 and self.cell2 is c2 or 
             self.cell2 is c1 and self.cell3 is c2 or 
             self.cell3 is c1 and self.cell1 is c2 ) : 
                return True
        return False
    
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
  
    def oppositeCell(self, c) : # find 2 cells opposite of a cell in a triangle
        if self.cell1 is c :
            return  [ self.cell2, self.cell3 ]
        if self.cell2 is c :
            return  [ self.cell3, self.cell1 ]
        if self.cell3 is c :
            return  [ self.cell1, self.cell2 ]
        return None
  
    def oppositeLink(self, c) : # find opposite link of a cell in a triangle
        if not self.link1.contains(c) :
            return self.link1
        if not self.link2.contains(c) :
            return self.link2
        if not self.link3.contains(c) :
            return self.link3
        return None
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)
        mesh.deleteFace(self.face)
my implementation (Python version / Processing 3)
add_library('toxiclibs')
add_library('peasycam')
from toxi.physics3d.behaviors import GravityBehavior3D, AttractionBehavior3D
from toxi.physics3d import VerletPhysics3D, VerletParticle3D, VerletSpring3D
import random as rdm

nCells = 360
growthInterval = 10
divisiointervalserval = 50
repulsiveForce = .4
nRadii = 5
springForce = .001

colors = [(random(255), random(255), random(255)) for i in xrange(nCells)]

def setup():
    perspective(60 * DEG_TO_RAD, width/float(height), 2, 6000)
    size(700, 400, P3D)
    fill(255, 40, 70)
    smooth(8)
    
    global physics, cells
    
    # Array containing the "Cell" objects
    cells = [Cell(0, 0, 0)]
    
    cam = PeasyCam(this, 80)
    cam.setWheelScale(.3)
    
    physics = VerletPhysics3D()
    physics.setDrag(.1)
    physics.addParticle(cells[0])
    physics.addBehavior(AttractionBehavior3D(cells[0], cells[0].radius*nRadii, -repulsiveForce))
    
    
def draw():
    background(255)
    
    # Updating physics engine
    physics.update()

    for c in cells:
        c.grow()
        #c.interact()
        c.render()
        
        
    # Adjusting spring length (sum of the radii of the 2 cells forming a link)
    strokeWeight(1)
    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.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%divisiointervalserval == 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.add(Vec3D.randomVector().scale(self.radius))
                        child = Cell(newloc.x(), newloc.y(), newloc.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.add(Vec3D.randomVector().normalizeTo(self.radius))
                        child = Cell(newloc.x(), newloc.y(), newloc.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(frozenset([self.id, child.id, cells[self.links[0]].id]))
                     
                   
                    # Triangle state: forming a tetrahedron  
                    elif len(self.links) == 2:
                        
                         # downsizing mother cell when giving birth to child cell
                        self.radius *= .5
                        
                        # location of child cell
                        newloc = self.computeNormal([self.id, cells[self.links[0]].id, cells[self.links[1]].id])
                        child = Cell(newloc.x(), newloc.y(), newloc.z())
                        child.id = len(cells)
                        
                        self.connect(self, child, False)
                        self.connect(cells[self.links[0]], child, True)
                        self.connect(cells[self.links[1]], child, True)
                        
                        # Adding new face (triangle 1) information
                        for c in (self, child, cells[self.links[0]]):
                            c.faces.append(frozenset([self.id, child.id, cells[self.links[0]].id]))
                            
                        # Adding new face (triangle 2) information
                        for c in (self, child, cells[self.links[1]]):
                            c.faces.append(frozenset([self.id, child.id, cells[self.links[1]].id]))
                        
                        # Adding new face (triangle 3) information
                        for c in (cells[self.links[0]], child, cells[self.links[1]]):
                            c.faces.append(frozenset([cells[self.links[0]].id, child.id, cells[self.links[1]].id]))
                        
                        
                    # Tetrahedron state: edge division (core of the algorithm)
                    else:
                        
                        # downsizing mother cell when giving birth to child cell
                        self.radius *= .5
                        
                        # Picking 2 links at random and computing the interval between (number of faces)
                        intervals, vertices = self.pickInterval()
                        print intervals
                        
                        # Computing child cell location
                        n = self.computeNormal(list(vertices)) # centered normal pointing outward
                        dir = n.sub(self)
                        newloc = self.add(dir.normalizeTo(self.radius))
                    
                        child = Cell(newloc.x(), newloc.y(), newloc.z())
                        child.id = len(cells)
                        
                        # Creating new links 
                        for i, v in enumerate(vertices):
                            if i == 0: self.connect(self, child, False)
                            else: self.connect(cells[v], child, True)
                                
                            # Adding new face (triangle) information
                            for c in (cells[vertices[i]], child, cells[vertices[(i+1)%len(vertices)]]):
                                c.faces.append(frozenset([cells[vertices[i]].id, child.id, cells[vertices[(i+1)%len(vertices)]].id]))
                          
                        
                        # Remove former links (info + corresponding spring) between self and its selected neighboring cell
                        if intervals > 1:      
                            for v in vertices[2: -1]:
                                self.links.remove(cells[v].id)
                                cells[v].links.remove(self.id)
                                physics.removeSpring(physics.getSpring(self, cells[v]))
                                
                                
                        # Removing former face (triangle) information from each cell
                        for i, v in enumerate(vertices):
                            if i > 0 and i < len(vertices) - 1:
                                for c in (self, cells[vertices[i]], cells[vertices[i+1]]):
                                    c.faces.remove(frozenset([self.id, cells[vertices[i]].id, cells[vertices[i+1]].id]))      
                          


                # Increasing cell size intermittently
                if self.time%growthInterval == 0:
                    if self.radius < self.maxRadius:
                        self.radius += .1
                    
    
    def pickInterval(self):
        nLinks = len(self.links)
        nInterval = nLinks // 2
        r1 = int(random(nLinks))
        r2 = (r1 + nInterval) % nLinks
        l1 = self.links[r1]
        l2 = self.links[r2]
        vertices = [self.id, l1]
        oppcell = l1
        visited = []
            
        while oppcell != l2:
            faces = [f for f in cells[oppcell].faces if self.id in f]
            sface = faces[0] if faces[0] not in visited else faces[1]
            visited.append(sface)
            oppcell = iter(sface - set((self.id, oppcell))).next()
            vertices.append(oppcell)
    
        interval = len(vertices) - 2

        return interval, vertices                    
    
    
    
    def connect(self, c1, c2, added):
        
        if not added: # means that the child cell has not been included yet in the physics engine
            cells.append(c2)
            physics.addParticle(c2)
            physics.addBehavior(AttractionBehavior3D(c2, c2.radius*nRadii, -repulsiveForce))
        
        #Add spring between c1 and c2
        s = VerletSpring3D(c1, c2, c2.distanceTo(c1), springForce) 
        physics.addSpring(s)  
        
        #Storing linkage information for each cell
        c1.links.append(c2.id)
        c2.links.append(c1.id)  
        
    def interact(self):
        maxForce = 100

        for l in self.links:
            dif = self.sub(cells[l])
            force = (dif.magnitude() - (self.radius + cells[l].radius)) / (self.radius + cells[l].radius) * 300
            if force > maxForce: 
                force = maxForce
            elif force < -maxForce:
                force = -maxForce
            dif.normalizeTo(force)
            self.addForce(dif.scale(-1))
            cells[l].addForce(dif)
            
        
    def computeNormal(self, indices):
        
        # Calculate normal with Newell's method
        nml = PVector()
        
        for i in xrange(len(indices)):
                    
            p0 = cells[indices[i]] #current vertex
            p1 = cells[indices[(i+1)%len(indices)]] #next vertex
            
            nml.x += (p0.y() - p1.y()) * (p0.z() + p1.z())
            nml.y += (p0.z() - p1.z()) * (p0.x() + p1.x())
            nml.z += (p0.x() - p1.x()) * (p0.y() + p1.y())
            
        # Normalizing then scaling-up to see a line
        nml = nml.normalize()#.mult(1.1)
            
        # Center of the polygon
        center = PVector()
        for i in range(len(indices)):
            center.x += cells[indices[i]].x()
            center.y += cells[indices[i]].y()
            center.z += cells[indices[i]].z()
                
        center.div(len(indices))
        
        # Centroid of tetrahedron
        centroid = PVector()
            
        # Check if normal is heading inwards or outwards
        d = PVector.dot(center - centroid, nml)
        
        if d < 0: nml.mult(-1)
             
        return Vec3D(center.x + nml.x, center.y + nml.y, center.z + nml.z)

        
    def render(self):
        strokeWeight(4)
        point(self.x(), self.y(), self.z())
    
        beginShape(TRIANGLE)
        strokeWeight(1)
        stroke(0)
        for f in self.faces:
            f = list(f)
            fill(colors[f[0]][0], colors[f[1]][1], colors[f[2]][2])
            vertex(cells[f[0]].x(), cells[f[0]].y(), cells[f[0]].z())
            vertex(cells[f[1]].x(), cells[f[1]].y(), cells[f[1]].z())
            vertex(cells[f[2]].x(), cells[f[2]].y(), cells[f[2]].z())
        endShape()

(if you want to try the interact function (line 44), line 219 and 51 need to be commented out to avoid calling the VerletSpring class at the same time )

1 Like

From visual inspection, a guess: your upper shape “grows” by splitting two joined faces into four faces by adding a center point to their shared edge, then extruding that point for all four faces.

Yet it looks like your lower image isn’t geometrically the same. It appears to have geometry like skin tags – a group of faces held on to the main body by a single shared edge. If you were using a divide algorithm like above, that should be impossible.

Hi @jeremydouglass, thank you for the reply.

Your first guess is very close to the actual algorithm, the mesh “grows” by:

  • selecting 2 edges on a dividing cell
  • inserting a new child cell
  • creating 3 new edges and 2 new faces



I have checked, double-checked and triple-checked and I am confident that the implementation of the dividing process is correct. If there’s a mistake (and there must be one!) it should be somewhere else in the code.

In the original sketch it looks like the newly added cells stay on the surface of the shape whereas in my case they immediatly stand out from the mesh.
I’m pretty sure that doesn’t have to do with the repulsive force (cells have a repulsive force to prevent overlaps), neither it has to do with the birth location.

My take is that the physics of the spring connections is critical in the shaping of the mesh but I’m still unable to implement the formula used in the original sketch.

Hi @solub,

I tried to mess around with it but I don’t know ToxicLib and never used Python before so you can imagine I did not have lots of success ^^

But, I may have some leads for you (maybe):

  • When you create a new cell and compute the new location, I added +5 in the scale factor so the new particule is not overlapping the old one.
  • I displayed each particles as sphere instead of point and deleted the triangles to have a better view of what was going on and I noticed that when the particles are at a resting place (I limited the number of cells to be 2) they are not touching each others. So there is probably something wrong with the way the springs are working.
  • Something that I don’t understand is that you add some string but you also add some behavior. Shouldn’t be the strings enough ?
  • I locked the initial particle so it will never move because when playing with the strings and behavior the particles started to fly all over the place…

Also, one thing that annoyed me is that it seems that the simulation and the display are linked together. I wanted to create a mousePressed function in which I could add a new particle anytime I want and just keep the simulation running in the draw loop but I couldn’t do it.

1 Like

@jb4x thanks for the reply, I really appreciate you took the time to help.

Indeed each particle is given a repulsive force to prevent collisions / overlapping (that’s why they’re not touching). On the other hand we don’t want the particles to get too far from each others so we connect them with springs.

Here’s a version with spheres and a mousePressed function (just replaced line 68 but it’s not ideal, a mouseReleased version should be much more convenient).

add_library('toxiclibs')
add_library('peasycam')
from toxi.physics3d.behaviors import GravityBehavior3D, AttractionBehavior3D
from toxi.physics3d import VerletPhysics3D, VerletParticle3D, VerletSpring3D

nCells = 360
growthInterval = 10
divisiointervalserval = 50
repulsiveForce = .4
nRadii = 5
springForce = .001

colors = [(random(255), random(255), random(255)) for i in xrange(nCells)]

def setup():
    perspective(60 * DEG_TO_RAD, width/float(height), 2, 6000)
    sphereDetail(9)
    size(700, 400, P3D)
    #frameRate(3)
    fill(255, 40, 70)
    smooth(8)
    
    global physics, cells
    
    # Array containing the "Cell" objects
    cells = [Cell(0, 0, 0)]
    
    cam = PeasyCam(this, 80)
    cam.setWheelScale(.1)
    
    physics = VerletPhysics3D()
    physics.setDrag(.4)
    physics.addParticle(cells[0])
    physics.addBehavior(AttractionBehavior3D(cells[0], cells[0].radius*nRadii, -repulsiveForce))
    
    
def draw():
    background(255)
    
    
    # Updating physics engine
    physics.update()
    
    for c in cells:
        c.grow()
        c.interact()
        c.render()
        

class Cell(VerletParticle3D):    
    def __init__(self, x, y, z):
        super(Cell, self).__init__(x, y, z)
        self.radius = 5
        self.maxRadius = 25
        self.id = 0
        self.time = 0
        self.links = []
        self.faces = []
        self.active = False
        
    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 keyPressed:
                if self.active:
                    
                    # 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.add(Vec3D.randomVector().scale(self.radius))
                        child = Cell(newloc.x(), newloc.y(), newloc.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.add(Vec3D.randomVector().normalizeTo(self.radius))
                        child = Cell(newloc.x(), newloc.y(), newloc.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(frozenset([self.id, child.id, cells[self.links[0]].id]))
                     
                   
                    # Triangle state: forming a tetrahedron  
                    elif len(self.links) == 2:
                        
                         # downsizing mother cell when giving birth to child cell
                        self.radius *= .5
                        
                        # location of child cell
                        newloc = self.computeNormal([self.id, cells[self.links[0]].id, cells[self.links[1]].id])
                        child = Cell(newloc.x(), newloc.y(), newloc.z())
                        child.id = len(cells)
                        
                        self.connect(self, child, False)
                        self.connect(cells[self.links[0]], child, True)
                        self.connect(cells[self.links[1]], child, True)
                        
                        # Adding new face (triangle 1) information
                        for c in (self, child, cells[self.links[0]]):
                            c.faces.append(frozenset([self.id, child.id, cells[self.links[0]].id]))
                            
                        # Adding new face (triangle 2) information
                        for c in (self, child, cells[self.links[1]]):
                            c.faces.append(frozenset([self.id, child.id, cells[self.links[1]].id]))
                        
                        # Adding new face (triangle 3) information
                        for c in (cells[self.links[0]], child, cells[self.links[1]]):
                            c.faces.append(frozenset([cells[self.links[0]].id, child.id, cells[self.links[1]].id]))
                        
                        
                    # Tetrahedron state: edge division (core of the algorithm)
                    else:
                        
                        # downsizing mother cell when giving birth to child cell
                        self.radius *= .5
                        
                        # Picking 2 links at random and computing the interval between (number of faces)
                        intervals, vertices = self.pickInterval()
                                                
                        # Computing child cell location
                        n = self.computeNormal(list(vertices)) # centered normal pointing outward
                        dir = n.normalizeTo(self.radius)
                        self.subSelf(dir)
                        newloc = self.add(dir)
                    
                        child = Cell(newloc.x(), newloc.y(), newloc.z())
                        child.id = len(cells)
                        
                        
                        # Creating new links 
                        for i, v in enumerate(vertices):
                            if i == 0: self.connect(self, child, False)
                            else: self.connect(cells[v], child, True)
                                
                            # Adding new face (triangle) information
                            for c in (cells[vertices[i]], child, cells[vertices[(i+1)%len(vertices)]]):
                                c.faces.append(frozenset([cells[vertices[i]].id, child.id, cells[vertices[(i+1)%len(vertices)]].id]))
                          
                        
                        # Remove former links (info + corresponding spring) between self and its selected neighboring cell
                        if intervals > 1:      
                            for v in vertices[2: -1]:
                                self.links.remove(cells[v].id)
                                cells[v].links.remove(self.id)
                                
                                
                        # Removing former face (triangle) information from each cell
                        for i, v in enumerate(vertices):
                            if i > 0 and i < len(vertices) - 1:
                                for c in (self, cells[vertices[i]], cells[vertices[i+1]]):
                                    c.faces.remove(frozenset([self.id, cells[vertices[i]].id, cells[vertices[i+1]].id]))      
                          
                        self.active = False
                        
                        
                if random(1) > .5:
                    self.active = True
                    
                    
                    
            # Increasing cell size intermittently
            if self.time%growthInterval == 0:
                if self.radius < self.maxRadius:
                    self.radius += .1
                                        
    
    def pickInterval(self):
        nLinks = len(self.links)
        nInterval = nLinks // 2
        r1 = int(random(nLinks))
        r2 = (r1 + nInterval) % nLinks
        l1 = self.links[r1]
        l2 = self.links[r2]
        vertices = [self.id, l1]
        oppcell = l1
        visited = []
            
        while oppcell != l2:
            faces = [f for f in cells[oppcell].faces if self.id in f]
            sface = faces[0] if faces[0] not in visited else faces[1]
            visited.append(sface)
            oppcell = iter(sface - set((self.id, oppcell))).next()
            vertices.append(oppcell)
    
        interval = len(vertices) - 2

        return interval, vertices                    
    
    
    
    def connect(self, c1, c2, added):
        
        if not added: # means that the child cell has not been included yet in the physics engine
            cells.append(c2)
            physics.addParticle(c2)
            physics.addBehavior(AttractionBehavior3D(c2, c2.radius*nRadii, -repulsiveForce))
        
        #Add spring between c1 and c2
        # s = VerletSpring3D(c1, c2, c2.distanceTo(c1), springForce) 
        # physics.addSpring(s)  
        
        #Storing linkage information for each cell
        c1.links.append(c2.id)
        c2.links.append(c1.id)  
        
    def interact(self):
        # Spring force (Hooke's law)
        for l in self.links:
            dif = cells[l].sub(self) 
            d = dif.magnitude()
            stretch = (d - (self.radius + cells[l].radius)) #/ (self.radius + cells[l].radius) * 300
            f = dif.normalize().scale(-1 * .05 * stretch)
            cells[l].addForce(f)


        
    def computeNormal(self, indices):
        
        # Calculate normal with Newell's method
        nml = PVector()
        
        for i in xrange(len(indices)):
                    
            p0 = cells[indices[i]] #current vertex
            p1 = cells[indices[(i+1)%len(indices)]] #next vertex
            
            nml.x += (p0.y() - p1.y()) * (p0.z() + p1.z())
            nml.y += (p0.z() - p1.z()) * (p0.x() + p1.x())
            nml.z += (p0.x() - p1.x()) * (p0.y() + p1.y())
            
        # Normalizing then scaling-up to see a line
        nml = nml.normalize()#.mult(1.1)
            
        # Center of the polygon
        center = PVector()
        for i in range(len(indices)):
            center.x += cells[indices[i]].x()
            center.y += cells[indices[i]].y()
            center.z += cells[indices[i]].z()
                
        center.div(len(indices))
        
        # Centroid of tetrahedron
        centroid = PVector()
            
        # Check if normal is heading inwards or outwards
        d = PVector.dot(center - centroid, nml)
        
        if d < 0: nml.mult(-1)
             
        #return Vec3D(center.x, center.y, center.z)
        return Vec3D(center.x + nml.x, center.y + nml.y, center.z + nml.z)

        
    def render(self):
        strokeWeight(4)
        point(self.x(), self.y(), self.z())
        
        pushMatrix()
        translate(self.x(), self.y(), self.z())
        strokeWeight(1)
        noStroke()
        fill(colors[self.id][0], colors[self.id][1], colors[self.id][2])
        sphere(self.radius)
        popMatrix()
    
        beginShape(TRIANGLE)
        strokeWeight(1)
        stroke(0)
        noFill()
        for f in self.faces:
            f = list(f)
            vertex(cells[f[0]].x(), cells[f[0]].y(), cells[f[0]].z())
            vertex(cells[f[1]].x(), cells[f[1]].y(), cells[f[1]].z())
            vertex(cells[f[2]].x(), cells[f[2]].y(), cells[f[2]].z())
        endShape()

Capture%20d%E2%80%99%C3%A9cran%20(20)

Also I’ve replaced toxiclibs VerletSprings by an interact function that computes the same Hooke’s law. The result is the same obviously but what I’d like to do now is implement the original sketch springs calculations (a bit different from Hooke’s law) (see the original question).

(Just realized I forgot to update the thread with the solution I had found)

The problem was towfold:

  • I had to keep the first and last edges instead of removing/rebuilding them (connecting order is important)
  • Planar force was missing
def planarForce(self):
    center, count = Vec3D(), 0
    for l in self.links:
        center.addSelf(cells[l])
        count += 1
                
    center.scaleSelf(1.0/count)
    dif = center.sub(self)
    f = dif.sub(self.getVelocity().scale(.3))
    self.addForce(f)


gif

2 Likes