Cell Division and Growth Algorithm: Adding springs

[solub]

Hi,

I’m desperately trying to connect growing cells with springs using the toxiclibs library. More specifically, I’d like to to put a spring between a “mother” cell and its “child” after division.

PROBLEMS:

• All cells are appended to an array list named collection and I don’t know how to identify which cell is a “mother” cell and which other cell is its “child”. I need to know that information in order to connect them

• I have a Cell class that inherits from VerletParticle2D but for some reason I’m unable to create springs from these verlet particles.

The final output should look like this:

add_library('verletphysics')
add_library('toxiclibscore')
from toxi.physics2d import VerletPhysics2D, VerletParticle2D

connections = []

def setup():
    global collection, physics
    size(480, 360, P2D)
    smooth(8)
    
    physics = VerletPhysics2D()
    collection = [Cell(width>>1, height>>1 , 20)]

        
def draw():
    global physics, collection
    background(255)
    
    physics.update()
        
    for i, c in enumerate(collection):
        if c.connected == False:
            physics.addParticle(c)
            c.connected = True

        c.growth()
        c.interact()
        c.display()

        for i2, c2 in enumerate(collection):
            if i != i2:
                s = VerletSpring2D(collection[i], collection[i2], 200, 1)
                physics.addSpring(s)
          
    
class Cell(VerletParticle2D):
    global physics
    growthInterval = 5
    divisionInterval = 50
    
    def __init__(self, x, y, rad):
        super(Cell, self).__init__(x, y)
        self.position = Vec2D(x, y)
        self.radius = rad
        self.time = 0
        self.connected = False
        
                                
    def growth(self):
        self.time += 1
        if len(collection) < 20:
            if self.time > 0 and self.time%Cell.divisionInterval == 0:
                if random(1) > .5:
                    dir = PVector.random2D()
                    self.radius *= .5
                    
                    child = Cell(self.position.x() + dir.x, self.position.y() + dir.y, self.radius)
                    
                    collection.append(child)
                    connections.append([self, child, False])

            
            if self.time%Cell.growthInterval == 0:
                self.radius += .5

                
    def interact(self):
        for a in collection:
            if a is not self:
                d = self.position.distanceTo(a.position)
                if d < (a.radius*.5 + self.radius*.5 + 1):
                    gap = (a.radius*.5 + self.radius*.5 + 1) - d
                    diff = Vec2D.sub(self.position, a.position).normalizeTo(gap)
                    a.position = a.position.sub(diff)                   

                                  
    def display(self):
        stroke(0)
        strokeWeight(3)
        point(self.position.x(), self.position.y())
        noFill()
        strokeWeight(1)
        ellipse(self.position.x(), self.position.y(), self.radius, self.radius)

[hamoid]

For the first problem: what if you just add a mother property to Cell and populate that property either in the constructor or by adding a setMother() function?

I don’t know enough Python to answer the second question.

[solub]

Hi @hamoid, thank you for the reply.

I eventually opted for a dictionnary to keep trace of the relations between ‘mother’ and ‘child’ cells.

Regarding the second question, trying to use the toxiclibs library was a nightmare (for this case specifically) and I decided to implement my own Spring class instead (based on this tutorial from Daniel Shiffman).

However the end result is far from what being perfect:

• springs get mixed-up
• some cells are overlapping
• others are jammed in an intertwining mess

It shouldn’t be that way obviously (cf. the example picture above) and I don’t understand the exact reason behind this chaotic behavior.

I suspect, however, the interact() function to be flawed (from line 130 to 138). Something might be off with the way I’m pushing the cells when they enter another cell’s perimeter:

   def interact(self):
        
        #Check every other cells
        for a in collection:
            if a is not self:
                
                #Distance between other cells
                d = PVector.dist(self.location, a.location)
                
                #If distance < sum of the 2 radii
                if d < (a.radius + self.radius + 2):
                    
                    #Push other cells outside its perimeter
                    gap = (a.radius + self.radius + 2) - d
                    diff = PVector.sub(self.location, a.location)
                    a.velocity.sub(diff.setMag(gap))

I would really appreciate if someone could help me find what’s going on here.

Full script with annotations

id = 0

def setup():
    global collection, d
    size(480, 360, P2D)
    smooth(8)
    
    #Array list containing cells
    collection = [Cell(width>>1, height>>1)]
    
    #Dictionnary to keep trace of the 'mother' <-> 'child' relations
    d = {}

        
def draw():
    background(255)
        
    #Lock cell 0
    collection[0].location = PVector(width>>1, height>>1)
    
    #Grow cells 
    for i, c in enumerate(collection):
        c.update()
        c.growth()
        c.interact()
        c.display() 
       
       
    ###Connecting 'mother' and 'children' cells with springs###
    
    for e in d:
        
        #Attach a spring to every mother cells (with a length of 10)
        s = Spring(collection[e].location.x, collection[e].location.y, 10)
        
        for id in d[e]:
            
            #Connect that spring to their children
            s.connect(collection[id])
            
            #Constrain length to the sum of the 2 radii (mother radius + child radius + sum of stroke weights)
            s.constrainLength(collection[id], 10, collection[id].radius + collection[e].radius + 4)
            
            #Draw springs
            s.displayLine(collection[id])

    
    #saveFrame("movie/mitosis_####.png")  
    
   
    
class Cell(object):
    growthInterval = 5
    divisionInterval = 50
    
    def __init__(self, x, y):
        self.location = PVector(x, y)
        self.velocity = PVector()
        self.acceleration = PVector()
        self.radius = 10
        self.time = 0
        self.id = 0
        self.mass = 5
        self.damping = .1
        
        
    def applyForce(self, f):
        f = f / self.mass
        self.velocity += f
        
        
    def update(self):
        self.velocity *= self.damping
        self.location += self.velocity
        self.velocity += self.acceleration
        self.acceleration *= 0

                                
    def growth(self):
        global id
        
        #Increment cell time
        self.time += 1
        
        #While less than 40 cells
        if len(collection) < 40:
            
            #If growing time is ellapsed
            if self.time > 0 and self.time%Cell.divisionInterval == 0:
                
                #Add randomness to the division process
                if random(1) > .5:
                    
                    #Increment id
                    id += 1
                    
                    #Divide radius by 2
                    self.radius *= .5
                    
                    #Find a random vector
                    dir = PVector.random2D()
                
                    #Add this vector to new cell (child) position
                    child = Cell(self.location.x + dir.x, self.location.y + dir.y)
                    
                    #Set child's id
                    child.id = id
                    
                    #Store children cells' id for each mother cell
                    if self.id not in d:
                        d[self.id] = [child.id]
                    else:
                        d[self.id].append(child.id)
                    
                    #Store cells in array list 'collection' 
                    collection.append(child)
                    
            #Grow cell every 5 iterations 
            if self.time%Cell.growthInterval == 0:
                self.radius += .5
        
        
    def interact(self):
        
        #Check every other cells
        for a in collection:
            if a is not self:
                
                #Distance between other cells
                d = PVector.dist(self.location, a.location)
                
                #If distance < sum of the 2 radii
                if d < (a.radius + self.radius + 2):
                    
                    #Push other cells outside its perimeter
                    gap = (a.radius + self.radius + 2) - d
                    diff = PVector.sub(self.location, a.location)
                    a.velocity.sub(diff.setMag(gap))
          
                                          
    def display(self):
        stroke(0)
        strokeWeight(4)
        point(self.location.x, self.location.y)
        noFill()
        strokeWeight(1)
        ellipse(self.location.x, self.location.y, self.radius*2, self.radius*2)
   
        
                  
class Spring(object):
    def __init__(self, x, y, l):
        self.anchor = PVector(x, y)
        self.k = .2
        self.len_ = l
        
    def connect(self, b):
        force = PVector.sub(b.location, self.anchor)
        d = force.mag()
        stretch = d - self.len_
        
        force.normalize()
        force *= -1 * self.k * stretch
        b.applyForce(force)
        
    def constrainLength(self, b, minl, maxl):
        dir = PVector.sub(b.location, self.anchor)
        d = dir.mag()
        
        if d < minl:
            dir.normalize()
            dir *= minl
            b.location = PVector.add(self.anchor, dir)
            b.velocity *= 0
        if d > maxl:
            dir.normalize()
            dir *= maxl
            b.location = PVector.add(self.anchor, dir)
            b.velocity *= 0
        
    def displayLine(self, b):
        strokeWeight(1)
        stroke(255, 20, 30)
        line(b.location.x, b.location.y, self.anchor.x, self.anchor.y)
        

def keyPressed():
    noLoop()

[solub]

Trying again with the toxiclibs library.

All I’d like to know is:

• How to push surrounding cells when they enter another cell’s perimeter ? (What’s the math ?)

Context:

A large part of my problem comes from the fact that the cells’ sizes are not fixed, they change with time.
As a result, both the spring rest length and the radius of the repulsion force of each cell must vary in time.

Regarding the force, I could certainly use the AttractionBehavior class in the toxiclibs library instead of implementing my own BUT it would be then impossible to change the radius of the repulsion force dynamically (because the setRadius method doesn’t seem to be working).

As a result, cells are overlapping (fixed repulsion force radius that doesn’t fit the cells’ radius) and it quickly becomes a mess.

See animation

That’s why I need to compute the force by myself and would need your help to fix the interact() function.

 def interact(self):

        #Check every other cells
        for a in collection:
            if a is not self:

                #Distance between other cells
                d = dist(self.x(), self.y(), a.x(), a.y())

                #If distance < sum of the 2 radii
                if d < (a.radius + self.radius + 1):

                    #(NOT CORRECT) Push other cells outside its perimeter 
                    gap = (a.radius + self.radius + 1) - d
                    diff = Vec2D.sub(Vec2D(self.x(), self.y()), Vec2D(a.x(), a.y())).normalizeTo(gap)
                    a.addVelocity(diff)

(for a “Processing version” of that snippet (no use of toxicilibs classes) see previous post above)

Full script

add_library('verletphysics')
add_library('toxiclibscore')
from toxi.physics2d import VerletPhysics2D, VerletParticle2D
from toxi.physics.behaviors import GravityBehavior

id = 0

def setup():
    global collection, physics, d
    size(480, 360, P2D)
    smooth(8)
        
    physics = VerletPhysics2D()
    physics.setDrag(.4)
    physics.setWorldBounds(Rect(0, 0, width, height))
    
    
    collection = [Cell(width>>1, height>>1)]
    physics.addParticle(collection[0])
    
    d = {}

        
def draw():
    global physics, collection, connections, col
    background(255)
    
    physics.update()
    collection[0].lock()
    
    
    for c in collection:
        c.growth()
        c.display()
        c.interact()
        
    #Updating spring's rest length 

    if len(collection) > 1:
       for i in range(len(physics.springs)):
            physics.springs.get(i).setRestLength(collection[-2].radius + collection[-1].radius)
            
 
    for e in d:
        for id in d[e]:
            
            strokeWeight(1)
            line(collection[e].x(), collection[e].y(), collection[id].x(), collection[id].y())

    
class Cell(VerletParticle2D):
    growthInterval = 5
    divisionInterval = 50
    
    def __init__(self, x, y):
        super(Cell, self).__init__(x, y)
        self.position = Vec2D(x, y)
        self.radius = 5
        self.time = 0
        self.id = 0
        self.linked = False
        
                                
    def growth(self):
        global id
        self.time += 1
        if len(collection) < 20:
            if self.time > 0 and self.time%Cell.divisionInterval == 0:
                if random(1) > .5:
                    dir = Vec2D.randomVector()
                    self.radius *= .5
                    
                    id += 1
                    
                    child = Cell(self.x() + dir.x(), self.y() + dir.y())
                    child.id = id
                    
                    if self.id not in d:
                        d[self.id] = [child.id]
                    else:
                        d[self.id].append(child.id)
                    
                    collection.append(child)
                    
                    
                    #Adding child to 'physics'
                    physics.addParticle(child)
                    
                    #Adding springs between mother and child
                    s = VerletSpring2D(self, child, self.radius + child.radius + 10, .1)
                    physics.addSpring(s)
                    
            
            if self.time%Cell.growthInterval == 0:
                self.radius += .5

                
    def interact(self):
        for a in collection:
            if a is not self:
                d = dist(self.x(), self.y(), a.x(), a.y())
                if d < (a.radius + self.radius + 1):
                    gap = (a.radius + self.radius + 1) - d
                    diff = Vec2D.sub(Vec2D(self.x(), self.y()), Vec2D(a.x(), a.y())).normalizeTo(gap)
                    a.addVelocity(diff)
        
                                  
    def display(self):
        stroke(0)
        strokeWeight(3)
        point(self.x(), self.y())
        noFill()
        strokeWeight(1)
        ellipse(self.x(), self.y(), self.radius*2, self.radius*2)

[solub]

I think I’ve managed to compute the repulsion force but for some reason cells start to shake and overlap when they are attached to springs.

No springs, no overlapping -->

Springs: overlapping -->

This is all the more confusing that the springs rest lengths are updated at each iteration (mother radius + child radius).

I reallly need your help !

add_library('verletphysics')
add_library('toxiclibscore')
from toxi.physics2d import VerletPhysics2D, VerletParticle2D
from toxi.physics2d.behaviors import GravityBehavior

id = 0

def setup():
    global collection, physics, d
    size(690, 460, P2D)
    smooth(8)
        
    physics = VerletPhysics2D()
    physics.setWorldBounds(Rect(0, 0, width, height))
    
    
    collection = [Cell(width>>1, height>>1)]
    
    #Adding first particle to array list
    physics.addParticle(collection[0])
    
    #Setting a fixed repulsion force for this Cell only
    physics.addBehavior(AttractionBehavior(collection[0], collection[0].radius*5, -2.4))
    
    #Dictionnary to keep trace of 'mother' / 'child' relations
    d = {}
    

        
def draw():
    global physics, collection
    background(30)
    
    #Updating physics engine
    physics.update()
    
    #Locking first Cell
    collection[0].lock()
    
    
    #Making cells grow and repel each other
    for c in collection:
        c.growth()
        c.display()
        c.interact()
        
        
    #Updating rest lenghts
    if len(collection) > 1:
        for i in range(len(physics.springs)):
            physics.springs.get(i).setRestLength(collection[-2].radius + collection[-1].radius + 4)
            physics.springs.get(i).setStrength(.5)
   
            
 
    #Drawing connections
    for e in d:
        for id in d[e]:
            strokeWeight(2 - id*.006)
            stroke(80 + id*.6, 140 + id, 255 - id)
            line(collection[e].x(), collection[e].y(), collection[id].x(), collection[id].y())


    
class Cell(VerletParticle2D):
    growthInterval = 5
    divisionInterval = 50
    
    def __init__(self, x, y):
        super(Cell, self).__init__(x, y)
        self.position = Vec2D(x, y)
        self.radius = 5
        self.mass = self.radius
        self.time = 0
        self.id = 0        
        
                                
    def growth(self):
        global id
        self.time += 1
        if len(collection) < 40 :
            if self.time > 0 and self.time%Cell.divisionInterval == 0:
                
                if random(1) > .5:
                    
                    #Divide radius by 2
                    self.radius *= .5
                    
                    #Update id
                    id += 1
                    
                    #Create random vector
                    dir = Vec2D.randomVector()
            
                    #Adding random vector to child position
                    child = Cell(self.x() + dir.x(), self.y() + dir.y())
                    
                    #Setting child's id
                    child.id = id
                    
                    #Keep trace of 'mother' <-> 'child' relations
                    if self.id not in d:
                        d[self.id] = [child.id]
                    else:
                        d[self.id].append(child.id)
                    
                    #Append cell to array list
                    collection.append(child)
                    
                    #Adding cell to physics engine
                    physics.addParticle(child)
                    
                    #Adding spring between mother and child
                    s = VerletSpring2D(self, child, self.radius + child.radius, .5)
                    physics.addSpring(s)
                    
                
            
            if self.time%Cell.growthInterval == 0:
                self.radius += .5

                
    def interact(self):
        
        #Check every other cells
        for a in collection:
            if a is not self:
                
                #Distance between other cells
                d = dist(a.x(), a.y(), self.x(), self.y())
                if d < (a.radius + self.radius):
                    
                    ### Push other cells outside its perimeter ###
                                        
                    #Direction
                    dir = PVector.sub(PVector(self.x(), self.y()), PVector(a.x(), a.y()))
                    dis = dir.magSq()
                    dis = constrain(dis, 10, (self.radius + a.radius))
                    dir.normalize()
                    
                    #Magnitude
                    m = (self.mass * a.mass) / dis
                    
                    #Combining magnitude and direction
                    dir.mult(-m)
                    
                    #Converting to Vec2D + adding force to physics engine
                    dir = Vec2D(dir.x, dir.y)
                    a.addForce(dir)
                   
                                  
    def display(self):
        stroke(255)
        strokeWeight(3)
        point(self.x(), self.y())
        fill(300 - self.id*.4, 80 + self.id*.3, 120 + self.id, 10)
        strokeWeight(.7)
        stroke(300 - self.id*.5, 80 + self.id*.3, 120 + self.id, map(self.time, 0, 20, 0, 255))
        ellipse(self.x(), self.y(), self.radius*2, self.radius*2)

[hamoid]

To make it run I had to add 2D to the second import and add a third one:

from toxi.physics2d.behaviors import AttractionBehavior2D

the replace AttractionBehavior with AttractionBehavior2D.

Question: why the following -2 and -1?

        physics.springs.get(i).setRestLength(collection[-2].radius + collection[-1].radius + 4)

Wouldn’t that reduce the rest length of all springs to that of the smallest shapes?

[solub]

Hi hamoid,

@hamoid

To make it run I had to add 2D to the second import and add a third one:

from toxi.physics2d.behaviors import AttractionBehavior2D

Thanks for spotting this issue. On my computer AttractionBehavior2D doesn’t exist so I have to stick with AttractionBehavior. Anyway that import is not really important since I’m using it for setting the behavior of the very first cell only.

You’re right , collection[-1] and collection[-2] were refering to the wrong cell objects. I’m now iterating through another array list that contains every mother-child couples to fit the order of physics.springs. (There must be a clever way to do this though).

Yet, I still end-up with overlapping cells.

I truely wonder how Satoru Sugihara (author of the iGeo library) did to avoid any overlapping in his sketch

[hamoid]

It looks like now consecutive cells no longer overlap, right? The overlapping now happens across branches maybe? Is it comparing every cell with every other cell? or only with the next/previous one?

[solub]

@hamoid

It looks like now consecutive cells no longer overlap, right?

True !

@hamoid

The overlapping now happens across branches maybe?

Yes

@hamoid

Is it comparing every cell with every other cell? or only with the next/previous one?

It IS comparing every cell with every other cell:

def interact(self):
        
        #Check every other cells
        for a in collection:
            if a is not self:
                
                #Distance between other cells
                d = dist(a.x(), a.y(), self.x(), self.y())
                if d < (a.radius + self.radius):
                    
                    ### Push other cells outside its perimeter ###
                                        
                    #Direction
                    dir = PVector.sub(PVector(self.x(), self.y()), PVector(a.x(), a.y()))
                    dis = dir.magSq()
                    dis = constrain(dis, 10, (self.radius + a.radius))
                    dir.normalize()
                    
                    #Magnitude
                    m = (self.mass * a.mass) / dis
                    
                    #Combining magnitude and direction
                    dir.mult(-m)
                    
                    #Converting to Vec2D + adding force to physics engine
                    dir = Vec2D(dir.x, dir.y)
                    a.addForce(dir)

[hamoid]

Mmm… looking at the image… where should those poor cells go to? I haven’t studied the code, but maybe it’s strictly maintaining the separations with those springs. And unless the separations can be extended further, there’s no place where the cells could go to to avoid overlap just by rotating branches. Maybe that’s the issue?

As a test, you could lower the forking probability and see if that solves it (for now).

[jeremydouglass]

A guess: perhaps you are not replicating what IParticle.push() does correctly.

You are looping through cells and allowing them to push other cells outside their perimeter. However, that act allows them to push cells into another cell perimeter. If your loop order is stable, that will create at least one winner but potentially many losers in the get-out-of-my-perimeter game – the cells that get checked last in their neighborhood win, everyone else (potentially) loses.

I haven’t examined your reference, but brainstorming a few approaches to handle that might include (?):

1. walk from the root node (which pushes anything) down to level level 2 – these can’t push the root, but they can push anything else. Then level 3 – these can’t push anything from level 2 or the root. Since only peers and leaf nodes are pushable, overlaps should propagate out. That makes sense assuming you have a root node heirarchy, of course – if you have a mandala it doesn’t work.

2. Rather than setting something outside the perimeter as a response to the overlap, move it proportionate to the overlap – e.g. move it 0.2 of the overlap distance. This means that strong overlaps will produce more displacement force than weak overlaps. You might want to combine that kind of approach with:

3. shuffle the cells array before checking for interactions. If there is no stable voting order, there aren’t permanently weak and strong cells that can get stuck in stable bad configurations, so solutions are more likely. Keep in mind, if your springs are too strong then some configurations do not have non-overlapping solutions. To avoid those, you might need to:

4. only create a new node in place if there is available non-overlapping space for it – which might involve displacing all nodes to the left or right by x in order to make space at the insertion point, then allowing the springs to snap them back. In essence, this is an already-valid-compactor, rather than a fix-the-invalid expander. On no frame will there be overlap.