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.
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
VerletSpringsclass (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()= toxiclibssub() - iGeo
len()= toxiclibsmagnitude() - iGeo
len(force)= toxiclibsnormalizeTo(force)(similar tosetMag(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 toxiclibsaddForce()
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
interactfunction (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 )
