Quick and dirty marching tetraheda isosurface thing

IsoSurface.cs

using UnityEngine;
using System.Collections;
using System;
using System.Collections.Generic;

[RequireComponent (typeof(MeshRenderer))]
[RequireComponent (typeof(MeshFilter))]
public class IsoSurface : MonoBehaviour 
{
	
	public Mesh inputMesh;
	public Vector3 meshOffset = new Vector3(0,0,0);
	public float meshScale = 1;
	
	public class SimpleMesh
	{
		public float[][] vertices;
		public int[][] faces;

		public int[] tris;

		public SimpleMesh()
		{

		}

		public void Prep()
		{
			var triList = new List<int>();
			for(var f = 0; f<faces.Length; f++)
			{
				var face = faces[f];

				if(face.Length == 4)
				{
					triList.Add(face[0]); triList.Add(face[1]); triList.Add(face[2]);
					triList.Add(face[0]); triList.Add(face[2]); triList.Add(face[3]);
				}
				else 
				{
					triList.Add(face[0]); triList.Add(face[1]); triList.Add(face[2]);
				}
			}

			tris = triList.ToArray();
		}
	}

	[ContextMenu ("Calculate")]
	void Calculate()
	{
		Func<float,float,float,float> sphere = (float x, float y, float z) =>
		{
			return Mathf.Sqrt(x*x+y*y+z*z) - 7.0f;
		};

		Func<float,float,float,float> cube = (float x, float y, float z) =>
		{
			float b = 5;
			if(x < -b || x > b || y < -b || y > b || z < -b || z > b) return 1;
			return -1;
		};

		var trues = 0;
		var falses = 0;

		Func<float,float,float,float> meshify = (float x, float y, float z) =>
		{
			var checkPoint = new Vector3(x,y,z)*meshScale;

			if(Physics.CheckSphere(checkPoint,10f))
			{
				trues++;
				return -1;
			}
			else
			{
				falses++;
				return 1;
			}
		};

		var func = sphere;

		if(inputMesh != null)
		{
			func = meshify;
		}

		var sourceMesh = Calculate(new int[]{32,32,32},func,new float[][]{new float[]{-1000,-1000,-1000},new float[]{1000,1000,1000}});

		sourceMesh.Prep();

		var sourceVerts = sourceMesh.vertices;
		var sourceTris = sourceMesh.tris;

		var mesh = new Mesh();

		var verts = new Vector3[sourceVerts.Length];

		for(var v = 0; v < sourceVerts.Length; v++)
		{
			verts[v] = new Vector3(sourceVerts[v][0],sourceVerts[v][1],sourceVerts[v][2]);
		}

		mesh.vertices = verts;
		mesh.triangles = sourceTris;

		GetComponent<MeshFilter>().sharedMesh = mesh;

		Debug.Log ("Calc done! Verts:" + verts.Length + " T: " + trues + " F: " + falses);
	}

	SimpleMesh Calculate(int[] dims, Func<float,float,float,float> potential, float[][] bounds)
	{

		var cube_vertices = new float[][] 
		{
			new float[]{0,0,0},
			new float[]{1,0,0},
			new float[]{1,1,0},
			new float[]{0,1,0},
			new float[]{0,0,1},
			new float[]{1,0,1},
			new float[]{1,1,1},
			new float[]{0,1,1}
		};


		var tetra_list = new []
		{
			new []{0,2,3,7},
			new []{0,6,2,7},
			new []{0,4,6,7},
			new []{0,6,1,2},
			new []{0,1,6,4},
			new []{5,6,1,4} 
		};

		if(bounds == null) 
		{
			bounds = new float[][] {new float[]{0,0,0}, new float[]{dims[0],dims[1],dims[2]}};
		}

		var scale = new float[] {0,0,0};
		var shift = new float[] {0,0,0};

		for(var i = 0; i<3; i++)
		{
			scale[i] = (bounds[1][i] - bounds[0][i]) / (float)dims[i];
			shift[i] = bounds[0][i];
		}

		var vertices = new List<float[]>();
		var faces = new List<int[]>();
		var n = 0;
		var grid = new float[8];
		var edges = new int[12];
		var x = new float[] {0,0,0};

		Func<int,int,int> interp = (int i0, int i1) =>
		{
			var g0 = grid[i0];
			var g1 = grid[i1];
			var p0 = cube_vertices[i0];
			var p1 = cube_vertices[i1];
			var v = new float[] {x[0],x[1],x[2]};
			var t = g0 - g1;

			if(Math.Abs(t) > 1e-6)
			{
				t = g0 / t;
			}

			for(var i = 0; i<3; i++)
			{
				v[i] = scale[i] * (v[i] + p0[i] + t * (p1[i] - p0[i])) + shift[i];
			}

			vertices.Add(v);
			return vertices.Count-1;
		};

		for(x[2]=0; x[2]<dims[2]-1; ++x[2], n+=dims[0])
			for(x[1]=0; x[1]<dims[1]-1; ++x[1], ++n)
				for(x[0]=0; x[0]<dims[0]-1; ++x[0], ++n) 
		{
				for(var i=0; i<8; ++i) 
				{
					var cube_vert = cube_vertices[i];
					grid[i] = potential
					(
						scale[0]*(x[0]+cube_vert[0])+shift[0],
						scale[1]*(x[1]+cube_vert[1])+shift[1],
						scale[2]*(x[2]+cube_vert[2])+shift[2]
					);
				}

				for(var i=0; i<tetra_list.Length; i++) 
				{
					var T = tetra_list[i];
					var triindex = 0;
					if (grid[T[0]] < 0) triindex |= 1;
					if (grid[T[1]] < 0) triindex |= 2;
					if (grid[T[2]] < 0) triindex |= 4;
					if (grid[T[3]] < 0) triindex |= 8;

					switch (triindex) {
					case 0x00:
					case 0x0F:
						break;
					case 0x0E:
						faces.Add(new int[]{ 
						            interp(T[0], T[1])
						            , interp(T[0], T[3])
						            , interp(T[0], T[2])});
						break;
					case 0x01:
						faces.Add(new int[]{
						            interp(T[0], T[1])
						            , interp(T[0], T[2])
									, interp(T[0], T[3])});
						break;
					case 0x0D:
						faces.Add(new int[]{
						            interp(T[1], T[0])
						            , interp(T[1], T[2]) 
									, interp(T[1], T[3])});
						break;
					case 0x02:
						faces.Add(new int[]{
						            interp(T[1], T[0])
						            , interp(T[1], T[3])
									, interp(T[1], T[2])});
						break;
					case 0x0C:
						faces.Add(new int[]{
						            interp(T[1], T[2])
						            , interp(T[1], T[3])
						            , interp(T[0], T[3])
									, interp(T[0], T[2])});
						break;
					case 0x03:
						faces.Add(new int[]{
						            interp(T[1], T[2])
						            , interp(T[0], T[2])
						            , interp(T[0], T[3])
									, interp(T[1], T[3])});
						break;
					case 0x04:
						faces.Add(new int[]{
						            interp(T[2], T[0])
						            , interp(T[2], T[1])
									, interp(T[2], T[3])});
						break;
					case 0x0B:
						faces.Add(new int[]{
						            interp(T[2], T[0])
						            , interp(T[2], T[3]) 
									, interp(T[2], T[1])});
						break;
					case 0x05:
						faces.Add(new int[]{
						            interp(T[0], T[1])
						            , interp(T[1], T[2])
						            , interp(T[2], T[3])
									, interp(T[0], T[3])});
						break;
					case 0x0A:
						faces.Add(new int[]{
						            interp(T[0], T[1])
						            , interp(T[0], T[3])
						            , interp(T[2], T[3])
									, interp(T[1], T[2])});
						break;
					case 0x06:
						faces.Add(new int[]{
						            interp(T[2], T[3])
						            , interp(T[0], T[2])
						            , interp(T[0], T[1])
									, interp(T[1], T[3])});
						break;
					case 0x09:
						faces.Add(new int[]{
						            interp(T[2], T[3])
						            , interp(T[1], T[3])
						            , interp(T[0], T[1])
									, interp(T[0], T[2])});
						break;
					case 0x07:
						faces.Add(new int[]{
						            interp(T[3], T[0])
						            , interp(T[3], T[1])
									, interp(T[3], T[2])});
						break;
					case 0x08:
						faces.Add(new int[]{
						            interp(T[3], T[0])
						            , interp(T[3], T[2])
									, interp(T[3], T[1])});
						break;
					}
				}

		}

		var mesh = new SimpleMesh();
		mesh.faces = faces.ToArray();
		mesh.vertices = vertices.ToArray();

		return mesh;

	}

}