MarchingCubes.compute

#pragma kernel Calculate

// Define fallof function to be used. Options: GENERIC_METABALL_FUNCTION, ELECTRIC_POTENTIAL_FUNCTION
#define GENERIC_METABALL_FUNCTION

//
// structs and buffers associated with them
//

struct Values {
	float v0, v1, v2, v3, v4, v5, v6, v7;
};

struct Cubes {
	float3 centerPos;
	float3 v0, v1, v2, v3, v4, v5, v6, v7;
};

struct Vertices {
	int index;
	float3 v0, v1, v2, v3, v4, v5, v6, v7, v8, v9, v10, v11;
};

struct Blob {
	float factor;
	float3 position;
};

StructuredBuffer<Cubes> cubes;
StructuredBuffer<Blob> blobs;
StructuredBuffer<float3> vertMap;

RWStructuredBuffer<Vertices> vertices;

//
// misc variables
//

int xCubes;
int yCubes;
int numBlobs;

float threshold;

//
// lookup tables
//

static uint vertTable[256] = {
	0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
	0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
	0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
	0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
	0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
	0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
	0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
	0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
	0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
	0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
	0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
	0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
	0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
	0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
	0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
	0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
	0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
	0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
	0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
	0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
	0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
	0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
	0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
	0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
	0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
	0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
	0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
	0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
	0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
	0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
	0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
	0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0
};

static int verticesAtEndsOfVerts[24] = {
	0,  1,
	1,  2,
	2,  3,
	3,  0,
	4,  5,
	5,  6,
	6,  7,
	7,  4,
	0,  4,
	1,  5,
	2,  6,
	3,  7
};

//
// helper functions
//

int findIndex(Values vertValue) {
	int cubeIndex = 0;

	if (vertValue.v0 > threshold) {
		cubeIndex += 1;
	}
	if (vertValue.v1 > threshold) {
		cubeIndex += 2;
	}
	if (vertValue.v2 > threshold) {
		cubeIndex += 4;
	}
	if (vertValue.v3 > threshold) {
		cubeIndex += 8;
	}
	if (vertValue.v4 > threshold) {
		cubeIndex += 16;
	}
	if (vertValue.v5 > threshold) {
		cubeIndex += 32;
	}
	if (vertValue.v6 > threshold) {
		cubeIndex += 64;
	}
	if (vertValue.v7 > threshold) {
		cubeIndex += 128;
	}

	return cubeIndex;
}

float getValueAtVert(Values val, int i) {
	if (i == 0) {
		return val.v0;
	}
	else if (i == 1) {
		return val.v1;
	}
	else if (i == 2) {
		return val.v2;
	}
	else if (i == 3) {
		return val.v3;
	}
	else if (i == 4) {
		return val.v4;
	}
	else if (i == 5) {
		return val.v5;
	}
	else if (i == 6) {
		return val.v6;
	}
	else {
		return val.v7;
	}
}

Values zeroValuesStruct() {
	Values val;
	val.v0 = 0;
	val.v1 = 0;
	val.v2 = 0;
	val.v3 = 0;
	val.v4 = 0;
	val.v5 = 0;
	val.v6 = 0;
	val.v7 = 0;

	return val;
}

//
// field function evaluation
//

float metaballFalloffFunction(float factor, float3 dist) {
#ifdef GENERIC_METABALL_FUNCTION
	return factor / (dist.x * dist.x + dist.y * dist.y + dist.z * dist.z);
#endif
}

Values evaluateFieldFunction(int pos) {
	Values vertValues = zeroValuesStruct();

	for (int i = 0; i < numBlobs; i++) {
		float3 dist = cubes[pos].v0 - blobs[i].position;
		vertValues.v0 += metaballFalloffFunction(blobs[i].factor, dist);

		dist = cubes[pos].v1 - blobs[i].position;
		vertValues.v1 += metaballFalloffFunction(blobs[i].factor, dist);

		dist = cubes[pos].v2 - blobs[i].position;
		vertValues.v2 += metaballFalloffFunction(blobs[i].factor, dist);

		dist = cubes[pos].v3 - blobs[i].position;
		vertValues.v3 += metaballFalloffFunction(blobs[i].factor, dist);

		dist = cubes[pos].v4 - blobs[i].position;
		vertValues.v4 += metaballFalloffFunction(blobs[i].factor, dist);

		dist = cubes[pos].v5 - blobs[i].position;
		vertValues.v5 += metaballFalloffFunction(blobs[i].factor, dist);

		dist = cubes[pos].v6 - blobs[i].position;
		vertValues.v6 += metaballFalloffFunction(blobs[i].factor, dist);

		dist = cubes[pos].v7 - blobs[i].position;
		vertValues.v7 += metaballFalloffFunction(blobs[i].factor, dist);
	}

	return vertValues;
}

//
// computer shader kernel
//

[numthreads(8, 8, 8)]
void Calculate(uint3 id : SV_DispatchThreadID) {
	int pos = id.x + xCubes * (id.y + yCubes * id.z);

	// find scalar values of a field on this position
	Values vertValues = evaluateFieldFunction(pos);

	// partial marching cubes algorithm
	int cubeIndex = findIndex(vertValues);
	int pattern = vertTable[cubeIndex];

	float3 verticesArr[12] = {
		float3(0, 0, 0),
		float3(0, 0, 0),
		float3(0, 0, 0),
		float3(0, 0, 0),
		float3(0, 0, 0),
		float3(0, 0, 0),
		float3(0, 0, 0),
		float3(0, 0, 0),
		float3(0, 0, 0),
		float3(0, 0, 0),
		float3(0, 0, 0),
		float3(0, 0, 0)
	};

	for (int vert = 0; vert < 12; vert++) {
		if ((pattern & 1 << vert) != 0) {
			int vertIndex0 = verticesAtEndsOfVerts[(int)(vert * 2)];
			int vertIndex1 = verticesAtEndsOfVerts[(int)(vert * 2 + 1)];

			float vertValue0 = getValueAtVert(vertValues, vertIndex0);
			float vertValue1 = getValueAtVert(vertValues, vertIndex1);

			float3 vertPosition0 = cubes[pos].centerPos + vertMap[vertIndex0];
			float3 vertPosition1 = cubes[pos].centerPos + vertMap[vertIndex1];

			float delta = (threshold - vertValue0) / (vertValue1 - vertValue0);

			// linear interpolate cubes for smooth surface
			verticesArr[vert] = vertPosition0 + delta * (vertPosition1 - vertPosition0);
		}
	}

	// save processed data to vertices buffer
	vertices[pos].index = cubeIndex;
	vertices[pos].v0 = verticesArr[0];
	vertices[pos].v1 = verticesArr[1];
	vertices[pos].v2 = verticesArr[2];
	vertices[pos].v3 = verticesArr[3];
	vertices[pos].v4 = verticesArr[4];
	vertices[pos].v5 = verticesArr[5];
	vertices[pos].v6 = verticesArr[6];
	vertices[pos].v7 = verticesArr[7];
	vertices[pos].v8 = verticesArr[8];
	vertices[pos].v9 = verticesArr[9];
	vertices[pos].v10 = verticesArr[10];
	vertices[pos].v11 = verticesArr[11];
}