Page MenuHome
Files F7873

implicit_marchingtetra.py
Public
Actions

implicit_marchingtetra.py
View Options

"""
This script is based on J M Soler's polygonize.py script, available at
file http://jmsoler/free.fr/util/blenderfile//py/polygonize04.py
webpage http://jmsoler.free.fr/didacticiel/blender/tutor/cpl_surfaceimplicite.htm
which in turn was adapted from a Kevin Loney's macro for Povray,
found in his public files directory:
http://www.geocities.com/qsquared_1999/pubfiles/pubfiles.html
This is basically a class that creates a Marching Tetrahedra object, to which you add
implicit objects for polygonisation.
The resulting surface is a mesh in my vefm format which can be converted to a Blender mesh.
"""
from implicit_vefm import vertex,face
class marchingtetra:
def __init__(self,bounds,resolution,threshold,mesh,interp,accuracy):
self.xmin=bounds[0]
self.xmax=bounds[1]
self.ymin=bounds[2]
self.ymax=bounds[3]
self.zmin=bounds[4]
self.zmax=bounds[5]
xres=resolution[0]
yres=resolution[1]
zres=resolution[2]
self.xres=xres
self.xstep=(self.xmax-self.xmin)/xres
self.ystep=(self.ymax-self.ymin)/yres
self.zstep=(self.zmax-self.zmin)/zres
self.corner=8*[[0.0,0.0,0.0]] ## Vertex coordinates of each corner of the current sample cube
self.fieldvalue=8*[[0.0,0.0,0.0]] ## Function values at each corner
self.threshold=threshold ## Global blob threshold value
self.mesh=mesh
self.implist=[] ## List of implicit objects to be polygonised
if interp==0:
self.interp=self.approx ## Interpolate values for fast,loose approximation of surface
else:
self.interp=self.exact ## Exact evaluation of surface.
self.accuracy=accuracy
def samplegrid(self):
## Call this method to evaluate and polygonise the implicit objects in the implist.
X=self.xmin
xcount=0
while X<self.xmax:
Y=self.ymin
while Y<self.ymax:
Z = self.zmin
while Z<self.zmax:
self.corner[0] = [X, Y, Z]
self.corner[1] = [X+self.xstep, Y, Z]
self.corner[2] = [X+self.xstep, Y, Z+self.zstep]
self.corner[3] = [X, Y, Z+self.zstep]
self.corner[4] = [X, Y+self.ystep, Z]
self.corner[5] = [X+self.xstep, Y+self.ystep, Z]
self.corner[6] = [X+self.xstep, Y+self.ystep, Z+self.zstep]
self.corner[7] = [X, Y+self.ystep, Z+self.zstep]
self.samplecube()
Z=Z+self.zstep
Y=Y+self.ystep
X=X+self.xstep
xcount+=1
print "Done a slice...",xcount,"/",self.xres
def samplecube(self):
## Does this sample cube contain any of the surface we're looking for?
lesser=0
greater=0
i = 0
while(i < 8):
self.fieldvalue[i] = self.calculatefield(self.corner[i][0], self.corner[i][1], self.corner[i][2]) #-self.threshold
if self.fieldvalue[i]<0.0:
lesser=lesser+1
if self.fieldvalue[i]>0.0:
greater=greater+1
i=i+1
## Found some surface, polygonise
if(lesser&greater):
self.dotetrahedrons()
def calculatefield(self,x,y,z):
## Evaluate the implicit objects in the list and return the sum of their functions
result=0
for candidate in self.implist:
field=candidate.calculate(vertex(x,y,z))
result=result+field
return result-self.threshold
def dotetrahedrons(self):
## Split the sample cube into six tetrahedrons, sample each to find surface
self.sampletetra(0, 2, 3, 7)
self.sampletetra(0, 2, 6, 7)
self.sampletetra(0, 4, 6, 7)
self.sampletetra(0, 6, 1, 2)
self.sampletetra(0, 6, 1, 4)
self.sampletetra(5, 6, 1, 4)
def sampletetra(self,v0, v1, v2, v3):
##
tetratype = 0;
newvert = 3*[0.0]
if(self.fieldvalue[v0] < 0):
tetratype = tetratype + 1
if(self.fieldvalue[v1] < 0):
tetratype = tetratype + 2
if(self.fieldvalue[v2] < 0):
tetratype = tetratype + 4
if(self.fieldvalue[v3] < 0):
tetratype = tetratype + 8
if tetratype==1 or tetratype==14:
newvert[0] = self.interp(self.fieldvalue[v0], self.fieldvalue[v1], self.corner[v0], self.corner[v1])
newvert[1] = self.interp(self.fieldvalue[v0], self.fieldvalue[v2], self.corner[v0], self.corner[v2])
newvert[2] = self.interp(self.fieldvalue[v0], self.fieldvalue[v3], self.corner[v0], self.corner[v3])
self.makeface(newvert[0], newvert[1], newvert[2])
elif tetratype==2 or tetratype==13:
newvert[0] = self.interp(self.fieldvalue[v1], self.fieldvalue[v0], self.corner[v1], self.corner[v0])
newvert[1] = self.interp(self.fieldvalue[v1], self.fieldvalue[v3], self.corner[v1], self.corner[v3])
newvert[2] = self.interp(self.fieldvalue[v1], self.fieldvalue[v2], self.corner[v1], self.corner[v2])
self.makeface(newvert[0], newvert[1], newvert[2])
elif tetratype==3 or tetratype==12:
newvert[0] = self.interp(self.fieldvalue[v0], self.fieldvalue[v3], self.corner[v0], self.corner[v3])
newvert[1] = self.interp(self.fieldvalue[v0], self.fieldvalue[v2], self.corner[v0], self.corner[v2])
newvert[2] = self.interp(self.fieldvalue[v1], self.fieldvalue[v3], self.corner[v1], self.corner[v3])
self.makeface(newvert[0], newvert[1], newvert[2])
newvert[0] = newvert[2]
newvert[2] = newvert[1]
newvert[1] = self.interp(self.fieldvalue[v1], self.fieldvalue[v2], self.corner[v1], self.corner[v2])
self.makeface(newvert[0], newvert[1], newvert[2])
elif tetratype==4 or tetratype==11:
newvert[0] = self.interp(self.fieldvalue[v2], self.fieldvalue[v0], self.corner[v2], self.corner[v0])
newvert[1] = self.interp(self.fieldvalue[v2], self.fieldvalue[v1], self.corner[v2], self.corner[v1])
newvert[2] = self.interp(self.fieldvalue[v2], self.fieldvalue[v3], self.corner[v2], self.corner[v3])
self.makeface(newvert[0], newvert[1], newvert[2])
elif tetratype==5 or tetratype==10:
newvert[0] = self.interp(self.fieldvalue[v0], self.fieldvalue[v1], self.corner[v0], self.corner[v1])
newvert[1] = self.interp(self.fieldvalue[v2], self.fieldvalue[v3], self.corner[v2], self.corner[v3])
newvert[2] = self.interp(self.fieldvalue[v0], self.fieldvalue[v3], self.corner[v0], self.corner[v3])
self.makeface(newvert[0], newvert[1], newvert[2])
newvert[2] = newvert[1]
newvert[1] = self.interp(self.fieldvalue[v1], self.fieldvalue[v2], self.corner[v1], self.corner[v2])
self.makeface(newvert[0], newvert[1], newvert[2])
elif tetratype==6 or tetratype==9:
newvert[0] = self.interp(self.fieldvalue[v0], self.fieldvalue[v1], self.corner[v0], self.corner[v1])
newvert[1] = self.interp(self.fieldvalue[v1], self.fieldvalue[v3], self.corner[v1], self.corner[v3])
newvert[2] = self.interp(self.fieldvalue[v2], self.fieldvalue[v3], self.corner[v2], self.corner[v3])
self.makeface(newvert[0], newvert[1], newvert[2])
newvert[1] = self.interp(self.fieldvalue[v0], self.fieldvalue[v2], self.corner[v0], self.corner[v2])
self.makeface(newvert[0], newvert[1], newvert[2])
elif tetratype==7 or tetratype==8:
newvert[0] = self.interp(self.fieldvalue[v3], self.fieldvalue[v2], self.corner[v3], self.corner[v2])
newvert[1] = self.interp(self.fieldvalue[v3], self.fieldvalue[v0], self.corner[v3], self.corner[v0])
newvert[2] = self.interp(self.fieldvalue[v3], self.fieldvalue[v1], self.corner[v3], self.corner[v1])
self.makeface(newvert[0], newvert[1], newvert[2])
def makeface(self,a, b, c):
## Create a face from the verts
v1=vertex(a[0],a[1],a[2])
v2=vertex(b[0],b[1],b[2])
v3=vertex(c[0],c[1],c[2])
self.mesh.verts.append(v1)
self.mesh.verts.append(v2)
self.mesh.verts.append(v3)
triangle=[v1,v2,v3]
newface=face(triangle)
self.mesh.faces.append(newface)
def approx(self,a, b, A, B):
approximate=[0.0,0.0,0.0]
for coord in range(3):
approximate[coord]=B[coord]-A[coord]
approximate[coord]=approximate[coord]*((-a)/(b-a))+A[coord]
return approximate
def exact(self,a, b, A, B):
if a<b:
lo=A
hi=B
else:
lo=B
hi=A
# print "lo",lo
# print "hi",hi
exactfunc=self.settle(lo,hi)
return exactfunc
def settle(self,low,high):
s1=(low[0]+high[0])/2.0
s2=(low[1]+high[1])/2.0
s3=(low[2]+high[2])/2.0
mid=[s1,s2,s3]
# print "hi",high
# print "low",low
# print "mid",mid
test=self.calculatefield(s1,s2,s3)
# print "test",test
if test>-self.accuracy and test<self.accuracy:
# print "final mid",mid
return mid
elif test>0.0:
return self.settle(low,mid)
elif test<0.0:
return self.settle(mid,high)

File Metadata

Mime Type
text/plain
Storage Engine
local-disk
Storage Format
Raw Data
Storage Handle
c3/cd/d5819ea2316f0ed5666c6404be5f

Event Timeline