"""

        z
        ^
        |

          4
        4-----5
       /|    /  5
      7-----6 |
      | |8  | |9
      | | 0 |10
      | 0---|-1   ->y
      |/  11|/ 1
      3-----2
     /   2
    <
   x

This Code is based on the example from
http://paulbourke.net/geometry/polygonise/
 by Paul Bourke,

 and modified by Tom Sapien to work with python within blender

 V0.1

 This code is under the terms of GPL

"""
__author__ = "Tom Sapiens"
__copyright__ = "Copyright 2011, Tom Sapiens"
__credits__ = ["Tom Sapiens", "Paul Bourke"]
__license__ = "GPL"
__version__ = "0.1"
__maintainer__ = "Tom Sapiens"
__email__ = "Tom.Sapiens@spot.colorado.edu"
__status__ = "alpha"


import time
import math
from math import sqrt,sin,cos,tan
import mathutils
import bpy
from itertools import chain
vec=mathutils.Vector
ABS=abs

def main():
    print("start calculation of isosurface")

#change this part to create your own surfaces
#####################################################
    # define a 3D scalarfield (the function which defines the shape of the isosurface)
    def scalarfield(pos):
        x,y,z=pos[0],pos[1],pos[2]
        m=2 #distance between spheres
        a= 1.0/(1+(x-m)*(x-m)+y*y+z*z)
        b= 1.0/(1+(x+m)*(x+m)+y*y+z*z)
        c= 0.5*(sin(6*x)+sin(6*z))
        csq=c**10
        return (a+b)-csq

    p0=-5,-5,-5             #first point defining the gridbox of the MC-algorithm
    p1=5,5,5                #second point defining the gridbox of the MC-algorithm
    res=100
    resolution=(res,res,res)   #resolution in x,y,z direction of the grid (10x10x10 means 1000 cubes)
    isolevel=0.3         #threshold value used for the surface within the scalarfield

#end of isosurface definition
###############################################


    isosurface(p0,p1,resolution,isolevel,scalarfield)

    print("end test"#   Determine the index into the edge table which
    #   tells us which vertices are inside of the surface
    cubeindex = 0
    if cornervalues[0] < isolevel: cubeindex = cubeindex | 1
    if cornervalues[1] < isolevel: cubeindex = cubeindex | 2
    if cornervalues[2] < isolevel: cubeindex = cubeindex | 4
    if cornervalues[3] < isolevel: cubeindex = cubeindex | 8
    if cornervalues[4] < isolevel: cubeindex = cubeindex | 16
    if cornervalues[5] < isolevel: cubeindex = cubeindex | 32
    if cornervalues[6] < isolevel: cubeindex = cubeindex | 64
    if cornervalues[7] < isolevel: cubeindex = cubeindex | 128

    # Cube is entirely in/out of the surface
    if edgetable[cubeindex] == 0: return []

    vertlist=[[]]*12
    # Find the vertices where the surface intersects the cube
    if (edgetable[cubeindex] & 1):    vertlist[0]  = vertexinterp(isolevel,cornerpos[0],cornerpos[1],cornervalues[0],cornervalues[1])
    if (edgetable[cubeindex] & 2):    vertlist[1]  = vertexinterp(isolevel,cornerpos[1],cornerpos[2],cornervalues[1],cornervalues[2])
    if (edgetable[cubeindex] & 4):    vertlist[2]  = vertexinterp(isolevel,cornerpos[2],cornerpos[3],cornervalues[2],cornervalues[3])
    if (edgetable[cubeindex] & 8):    vertlist[3]  = vertexinterp(isolevel,cornerpos[3],cornerpos[0],cornervalues[3],cornervalues[0])
    if (edgetable[cubeindex] & 16):   vertlist[4]  = vertexinterp(isolevel,cornerpos[4],cornerpos[5],cornervalues[4],cornervalues[5])
    if (edgetable[cubeindex] & 32):   vertlist[5]  = vertexinterp(isolevel,cornerpos[5],cornerpos[6],cornervalues[5],cornervalues[6])
    if (edgetable[cubeindex] & 64):   vertlist[6]  = vertexinterp(isolevel,cornerpos[6],cornerpos[7],cornervalues[6],cornervalues[7])
    if (edgetable[cubeindex] & 128):  vertlist[7]  = vertexinterp(isolevel,cornerpos[7],cornerpos[4],cornervalues[7],cornervalues[4])
    if (edgetable[cubeindex] & 256):  vertlist[8]  = vertexinterp(isolevel,cornerpos[0],cornerpos[4],cornervalues[0],cornervalues[4])
    if (edgetable[cubeindex] & 512):  vertlist[9]  = vertexinterp(isolevel,cornerpos[1],cornerpos[5],cornervalues[1],cornervalues[5])
    if (edgetable[cubeindex] & 1024): vertlist[10] = vertexinterp(isolevel,cornerpos[2],cornerpos[6],cornervalues[2],cornervalues[6])
    if (edgetable[cubeindex] & 2048): vertlist[11] = vertexinterp(isolevel,cornerpos[3],cornerpos[7],cornervalues[3],cornervalues[7])

    #Create the triangle
    triangles = []
    #for (i=0;triTable[cubeindex][i]!=-1;i+=3) {
    i=0
    while tritable[cubeindex][i] != -1:
        triangles.append([vertlist[tritable[cubeindex][i  ]],
                                   vertlist[tritable[cubeindex][i+1]],
                                   vertlist[tritable[cubeindex][i+2]]])
        i+=3

    return triangles

def vertexinterp(isolevel,p1,p2,valp1,valp2):
   if (ABS(isolevel-valp1) < 0.00001):
      return p1
   if (ABS(isolevel-valp2) < 0.00001):
      return p2
   if (ABS(valp1-valp2) < 0.00001):
      return p1
   mu = (isolevel - valp1) / (valp2 - valp1);
   p=vec([0,0,0])
   p.x = p1.x + mu * (p2.x - p1.x);
   p.y = p1.y + mu * (p2.y - p1.y);
   p.z = p1.z + mu * (p2.z - p1.z);

   return p

def genobject(objname,verts=[],faces=[],edges=[]):
    me = bpy.data.meshes.new(objname)  # create a new mesh
    me.from_pydata(verts,edges,faces)
    me.update()      # update the mesh with the new data
    ob = bpy.data.objects.new(objname,me) # create a new object
    ob.data = me          # link the mesh data to the object
    scene = bpy.context.scene           # get the current scene
    scene.objects.link(ob)                      # link the object into the scene
    return ob

def creategeometry(verts):
    faces=[]
    faceoffset=0
    for ver in verts:
        if len(ver)==4:
            faces.append((faceoffset+0,faceoffset+1,faceoffset+2,faceoffset+3))
            faceoffset+=4
        elif len(ver)==3:
            faces.append((faceoffset+0,faceoffset+1,faceoffset+2))
            faceoffset+=3
    return list(chain.from_iterable(verts)),faces

def genobjandremovedoubles(verts):
    verts,faces=creategeometry(verts)
    block=genobject("block",verts,faces)
    selectobj(block)
    bpy.ops.object.editmode_toggle()
    bpy.ops.mesh.select_all(action='SELECT')
    bpy.ops.mesh.remove_doubles(mergedist=0.01)
    bpy.ops.object.editmode_toggle()
    return block

def selectobj(obj):
    bpy.ops.object.select_all(action="DESELECT")
    obj.select=True
    bpy.context.scene.objects.active=obj

#a threshold function:
borders=[5,5,5]

def arange(start, stop, step):
     r = start
     while r < stop:
        yield r
        r += step

def cellloop(p0,p1,r):
    for z in arange(p0[2],p1[2],r[2]):
     for y in arange(p0[1],p1[1],r[1]):
      for x in arange(p0[0],p1[0],r[0]):
        yield x,y,z

def cornerloop(x,y,z):
    for cz in (0,z):
     for cy,cx in zip((0,y,y,0),(0,0,x,x)):
       yield cx,cy,cz

def isosurface(p0,p1,resolution,isolevel,isofunc):
    r=[(x1-x0)/sw for x0,x1,sw in zip(p0,p1,resolution)]

    triangles=[]
    for x,y,z in cellloop(p0,p1,r):
       cornervalues=[]
       cornerpos=[]
       for cx,cy,cz in cornerloop(r[0],r[1],r[2]):
          pos=vec((x+cx,y+cy,z+cz))
          cornerpos.append(pos)
          cornervalues.append(isofunc(pos))
       triangles.extend(polygonise(cornervalues,cornerpos,isolevel))

    genobjandremovedoubles(triangles)

if __name__=="__main__":
    main()