World's Tiniest* Marching Tetrahedron

mt.glsl

// World's Tiniest* Marching Tetrahedron
// by Leonard Ritter (leonard.ritter@duangle.com)

// this code is public domain

//      0
//      +
//    / | \
// 3 +-----+ 1
//    \ | /
//      +
//      2

const uvec4 tetbits = uvec4(0x487d210u, 0x11282844u, 0x2e71b7ecu, 0x3bde4db8u);
const uint tetbits_q = 0x16696994u;

// returns the mix factor x required to interpolate two vectors given two
// density values, such that the interpolated vector is at the position where
// a + (b-a)*x = 0
float tetlerp(float a, float b) {
    float d = (b - a);
    return (abs(d) < 1e-12)?0.0:(-a / d);
}

// takes the density values at the four endpoints (> 0 = outside / not on surface)
// and returns the indices of the triangle / quad to be generated
// returns 0 if no triangle is to be generated
// 1 if the result is a tri (interpolate edges x-y x-z x-w)
// 2 if the result is a quad (interpolate edges x-z x-w y-w y-z)
uint tetfaces(in vec4 d, out uvec4 i) {
    uint k = ((d[0] > 0.0)?2u:0u)
           | ((d[1] > 0.0)?4u:0u)
           | ((d[2] > 0.0)?8u:0u)
           | ((d[3] > 0.0)?16u:0u);
    i = (tetbits >> k) & 3u;
    return (tetbits_q >> k) & 3u;
}

const uvec4 btetbits = uvec4(0x1c008000u, 0x28404000u, 0x30c00000u, 0x14404000u);

// takes four binary values (0 = outside, 1 = surface / inside)
// returns w=1 if a triangle is to be generated,
// and the indices of the triangle in xyz
uvec4 btetfaces(in bvec4 d) {
    uint k = (d[0]?2u:0u)
           | (d[1]?4u:0u)
           | (d[2]?8u:0u)
           | (d[3]?16u:0u);
    return (btetbits >> k) & 3u;
}

/*

usage example for geometry shader:

// holds the four density values measured at the corners
vec4 d; 
// holds the coordinates of the four tetrahedral vertices
vec4 p[4];

uvec4 i;
uint c = tetfaces(d, i);
if (c == 1u) {
    // generate triangle
    gl_Position = mix(p[i.x],p[i.y],tetlerp(d[i.x],d[i.y]));
    EmitVertex();
    gl_Position = mix(p[i.x],p[i.z],tetlerp(d[i.x],d[i.z]));
    EmitVertex();
    gl_Position = mix(p[i.x],p[i.w],tetlerp(d[i.x],d[i.w]));
    EmitVertex();
    EndPrimitive();
} else if (c == 2u) {
    vec4 p0 = mix(p[i.x],p[i.z],tetlerp(d[i.x],d[i.z]));
    vec4 p2 = mix(p[i.y],p[i.w],tetlerp(d[i.y],d[i.w]));
    // generate quad
    gl_Position = p0;
    EmitVertex();
    gl_Position = mix(p[i.x],p[i.w],tetlerp(d[i.x],d[i.w]));
    EmitVertex();
    gl_Position = p2;
    EmitVertex();
    EndPrimitive();
    gl_Position = p2;
    EmitVertex();
    gl_Position = mix(p[i.y],p[i.z],tetlerp(d[i.y],d[i.z]));
    EmitVertex();
    gl_Position = p0;
    EmitVertex();
    EndPrimitive();
}

*/