– @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()
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).