Isosurface: Marching Tetrahedrons. Statt Marching Cubes.
Statt riesiger Case-Liste gibts bei Tetraedern effektiv nur 3 Fälle: alle Punkte auf derselben Seite (nichts rendern), 1 Punkt auf der anderen Seite (Dreieck rendern), 2 Punkte auf der anderen Seite (Quad rendern). Jeweils gedreht und negiert: 14 Renderfälle, keine doppeldeutigen.
Dafür entstehen 6 Tetras je Cube.
Die Fälle habe ich mir per Skizze selbst hergeleitet, weil selbermachen und so. Das ist im Gegensatz zu Marching Cubes (MC) tatsächlich überschaubar.
Bei Paul Bourke musste ich nur mal nachschauen, wie denn der Space überhaupt in Tetras aufgeteilt wird. Was ich mir auch nochmal für meine Zwecke skizziert habe.
(die oberen 80% befassen sich mit MC, ganz runterscrollen zu "Polygonising a Scalar Field Using Tetrahedrons")
Edit: Hier noch mein Testcode.
Bei mir ist vieles (z.B. SDF) eigentlich in komplexeren Klassen und -Systemen gekapselt; damit trotzdem was geht, habe ich hier alles Notwendige einfach prozedural durchgereicht und parametrisiert. Wenig optimiert, hoffentlich noch lesbar.
Wer noch weiter mit SDF-Funktionen spielen möchte: Inigo Quilez ist die Anlaufstelle dafür.
Vector4 und TriMesh sollten selbsterklärend sein.
Code: Alles auswählen
class MarchingTetraDemo
{
private List<float> blobs = new List<float>();
public TriMesh mesh = new TriMesh();
public MarchingTetraDemo()
{
float smooth = 1.0f;
float iso = 0.2f;
demoBlobs();
scan(5f, 0.5f, smooth, iso);
smoothNormals(smooth, iso);
}
//generate some random blobs
private void demoBlobs()
{
Random rnd = new Random();
float spread = 2.5f;
for (int i = 0; i < 40; i++)
{
blobs.Add((float)rnd.NextDouble() * spread * 2 - spread);
blobs.Add((float)rnd.NextDouble() * spread * 2 - spread);
blobs.Add((float)rnd.NextDouble() * spread * 2 - spread);
blobs.Add((float)rnd.NextDouble() * 0.4f + 0.3f);
}
}
//overwrite face normals
private void smoothNormals(float smooth, float iso)
{
mesh.normalList.Clear();
for (int i = 0; i < mesh.pointList.Count; i++)
{
mesh.normalList.Add(-getGradient(mesh.pointList[i], smooth, iso));
}
}
//--------------------------------- marching tetrahedrons algorithm
//scan isosurface
private void scan(float scanWidth, float cellSize, float smooth, float iso)
{
for (float x = -scanWidth; x < scanWidth; x += cellSize)
{
for (float y = -scanWidth; y < scanWidth; y += cellSize)
{
for (float z = -scanWidth; z < scanWidth; z += cellSize)
{
splitCube(new Vector4(x, y, z, 1), cellSize, smooth, iso);
}
}
}
}
//split cube into 6 tetrahedrons
private void splitCube(Vector4 p, float scanSize, float smooth, float iso)
{
float r = scanSize * 0.5f;
Vector4 v0 = p + new Vector4(-r, r, -r, 1);
Vector4 v1 = p + new Vector4(r, r, -r, 1);
Vector4 v2 = p + new Vector4(r, -r, -r, 1);
Vector4 v3 = p + new Vector4(-r, -r, -r, 1);
Vector4 v4 = p + new Vector4(-r, r, r, 1);
Vector4 v5 = p + new Vector4(r, r, r, 1);
Vector4 v6 = p + new Vector4(r, -r, r, 1);
Vector4 v7 = p + new Vector4(-r, -r, r, 1);
float d0 = unionDist(v0, smooth, iso);
float d1 = unionDist(v1, smooth, iso);
float d2 = unionDist(v2, smooth, iso);
float d3 = unionDist(v3, smooth, iso);
float d4 = unionDist(v4, smooth, iso);
float d5 = unionDist(v5, smooth, iso);
float d6 = unionDist(v6, smooth, iso);
float d7 = unionDist(v7, smooth, iso);
tetra(v3, v0, v1, v7, d3, d0, d1, d7);
tetra(v1, v2, v3, v7, d1, d2, d3, d7);
tetra(v6, v2, v1, v7, d6, d2, d1, d7);
tetra(v7, v6, v5, v1, d7, d6, d5, d1);
tetra(v5, v4, v7, v1, d5, d4, d7, d1);
tetra(v7, v4, v0, v1, d7, d4, d0, d1);
}
//generate iso polygon for a tetra cell
private void tetra(Vector4 v0, Vector4 v1, Vector4 v2, Vector4 v3, float d0, float d1, float d2, float d3)
{
int mask = 0;
if (d0 < 0) mask |= 1;
if (d1 < 0) mask |= 2;
if (d2 < 0) mask |= 4;
if (d3 < 0) mask |= 8;
//major case A: 0 vs. 4 (all inside or outside) -> do nothing
if (mask == 0 || mask == 15) return;
switch (mask)
{
//major case B: 1 vs. 3 -> triangle
case 1: addTri(v0, v1, v2, v3, d0, d1, d2, d3, 1); break;
case 2: addTri(v1, v0, v3, v2, d1, d0, d3, d2, 1); break;
case 4: addTri(v2, v0, v1, v3, d2, d0, d1, d3, 1); break;
case 8: addTri(v3, v0, v2, v1, d3, d0, d2, d1, 1); break;
//tri cases inverted
case 14: addTri(v0, v3, v2, v1, d0, d3, d2, d1, -1); break;
case 13: addTri(v1, v2, v3, v0, d1, d2, d3, d0, -1); break;
case 11: addTri(v2, v3, v1, v0, d2, d3, d1, d0, -1); break;
case 7: addTri(v3, v1, v2, v0, d3, d1, d2, d0, -1); break;
//major case C: 2 vs. 2 -> quad
case 3: addQuad(v0, v1, v2, v3, d0, d1, d2, d3); break;
case 5: addQuad(v0, v2, v3, v1, d0, d2, d3, d1); break;
case 6: addQuad(v1, v2, v0, v3, d1, d2, d0, d3); break;
//quad cases inverted
case 12: addQuad(v2, v3, v0, v1, d2, d3, d0, d1); break;
case 10: addQuad(v1, v3, v2, v0, d1, d3, d2, d0); break;
case 9: addQuad(v3, v0, v2, v1, d3, d0, d2, d1); break;
}
}
private void addTri(Vector4 v0, Vector4 v1, Vector4 v2, Vector4 v3, float d0, float d1, float d2, float d3, float dir)
{
if (dir < 0) mesh.addTriangle(surface(v0, v3, d0, d3), surface(v0, v2, d0, d2), surface(v0, v1, d0, d1));
else mesh.addTriangle(surface(v3, v0, d3, d0), surface(v2, v0, d2, d0), surface(v1, v0, d1, d0));
}
private void addQuad(Vector4 v0, Vector4 v1, Vector4 v2, Vector4 v3, float d0, float d1, float d2, float d3)
{
mesh.addQuad(surface(v2, v1, d2, d1), surface(v3, v1, d3, d1), surface(v3, v0, d3, d0), surface(v2, v0, d2, d0));
}
//interpolate edge
private Vector4 surface(Vector4 v0, Vector4 v1, float d0, float d1)
{
float f = -d0 / (d1 - d0);
return v0 + (v1 - v0) * f;
}
//--------------------------------- some sdf stuff
//classic sphere
private float sphereDist(Vector4 samplePos, Vector4 center, float rad, float iso = 0)
{
return (samplePos.Xyz - center.Xyz).Length - rad - iso;
}
//smooth union of spheres
private float unionDist(Vector4 samplePos, float smooth, float iso)
{
float d = float.MaxValue;
for (int i = 0; i < blobs.Count; i += 4)//data: x,y,z,rad
{
float d2 = sphereDist(samplePos, new Vector4(blobs[i], blobs[i + 1], blobs[i + 2], 1), blobs[i + 3]);
if (smooth > 0) d = smoothUnion(d, d2, smooth);
else d = Math.Min(d, d2);
}
return d;
}
//smoothing func
private float smoothUnion(float d1, float d2, float k)
{
float h = 0.5f + 0.5f * (d2 - d1) / k;
if (h < 0) h = 0;
if (h > 1) h = 1;
return d2 * (1 - h) + d1 * h - k * h * (1.0f - h);
}
//field gradient
private static float e = 0.001f;
private static Vector4 exp = new Vector4(e, 0, 0, 0);
private static Vector4 exn = new Vector4(-e, 0, 0, 0);
private static Vector4 eyp = new Vector4(0, e, 0, 0);
private static Vector4 eyn = new Vector4(0, -e, 0, 0);
private static Vector4 ezp = new Vector4(0, 0, e, 0);
private static Vector4 ezn = new Vector4(0, 0, -e, 0);
private Vector3 getGradient(Vector4 samplePos, float smooth, float iso)
{
float x0 = unionDist(samplePos + exp, smooth, iso);
float x1 = unionDist(samplePos + exn, smooth, iso);
float y0 = unionDist(samplePos + eyp, smooth, iso);
float y1 = unionDist(samplePos + eyn, smooth, iso);
float z0 = unionDist(samplePos + ezp, smooth, iso);
float z1 = unionDist(samplePos + ezn, smooth, iso);
Vector3 dir = new Vector3(x1 - x0, y1 - y0, z1 - z0).Normalized();
return dir;
}
}