I have looked up various ways to achieve mesh generation for formulae that describe 3D surfaces, and landed on Marching Cubes (Implementation pictured below)
I had considered generating a point cloud earlier and using the convex hull algorithm to generate a surface enclosing all points, however not every 3D surface plot has a single closed surface, so that idea had to be scrapped, and why I ended up with Marching Cubes.
Since Implicit surfaces also need to be supported, I settled on sampling points at fixed intervals between specific bounds, and then performing Marching Cubes on that data.
However, the mesh I generate:
- Roughly follows the contour of the surface (artifact of the Marching Cubes algorithm), and
- Has plenty of holes, which I suspect is due to precision issues.
I’m kind of stuck on how to optimize this.
I have heard of the “Marching Tetrahedra” algorithm, which could help get smoother surfaces, however I still have no idea how to manage accuracy of the surface w.r.t. the equation of the surface.
The only one way I could think of was performing gradient descent on every point after generation, which seems too computationally expensive.
Any suggestions would be highly appreciated.
Here’s the script I use to generate the surfaces (Mind you; I am new to C# , so please go easy on me :,) ) :
using System;
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
public delegate float OnSurfaceDelegate(float x, float y, float z);
public class MarchingCubes {
static readonly int[][] TRIANGULATIONS = {
new int[]{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
/* 254 other triangulation "edge" indices */
new int[]{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
};
static readonly (int, int, int)[] POINTS = {
(0, 0, 0),
(0, 0, 1),
(1, 0, 1),
(1, 0, 0),
(0, 1, 0),
(0, 1, 1),
(1, 1, 1),
(1, 1, 0),
};
static readonly (int, int)[] EDGES = {
(0, 1),
(1, 2),
(2, 3),
(3, 0),
(4, 5),
(5, 6),
(6, 7),
(7, 4),
(0, 4),
(1, 5),
(2, 6),
(3, 7),
};
OnSurfaceDelegate OnSurface;
float bounds;
float resolution;
List<List<List<float>>> samples;
public MarchingCubes( OnSurfaceDelegate OnSurface, float bounds = 5.0f, float resolution = 0.5f) {
this.OnSurface = OnSurface;
this.bounds = bounds;
this.resolution = resolution;
}
void SamplePoints() {
int samplesPerAxis = (int) (bounds * 2 / resolution) + 1;
samples = new List<List<List<float>>>();
for (float i = -bounds; i <= bounds; i += resolution) {
List<List<float>> temp2D = new ();
for (float j = -bounds; j <= bounds; j += resolution) {
List<float> temp1D = new ();
for (float k = -bounds; k <= bounds; k += resolution) {
temp1D.Add(OnSurface(i, j, k));
}
temp2D.Add(temp1D);
}
samples.Add(temp2D);
}
}
int TriangulationIndex(Vector3 coord) {
int idx = 0b00000000;
for (int i = 0; i < POINTS.Length; i++) {
var (xOffset, yOffset, zOffset) = POINTS[i];
idx |= (
samples[(int) ((coord.x + bounds) / resolution + xOffset)]
[(int) ((coord.y + bounds) / resolution + yOffset)]
[(int) ((coord.z + bounds) / resolution + zOffset)] <= 0.01f ? 1 : 0
) << i;
}
return idx;
}
public (Vector3[], int[]) GetMesh() {
List<Vector3> points = new ();
List<int> tris = new ();
SamplePoints();
for (int i = 0; i + 1 < samples.Count; i++) {
for (int j = 0; j + 1 < samples.Count; j++) {
for (int k = 0; k + 1 < samples.Count; k++) {
Vector3 coords = new Vector3(i, j, k) * resolution - new Vector3(bounds, bounds, bounds);
int triangleIndex = TriangulationIndex(coords);
for (int edgeIndex = 0; edgeIndex < TRIANGULATIONS[0].Length; edgeIndex += 3) {
if (TRIANGULATIONS[triangleIndex][edgeIndex] == -1) { continue; }
List<Vector3> singleTriangle = new ();
for (int offset = 2; offset >= 0; offset--) {
int edge = TRIANGULATIONS[triangleIndex][edgeIndex + offset];
(int left, int right) = EDGES[edge];
(int leftX, int leftY, int leftZ) = POINTS[left];
(int rightX, int rightY, int rightZ) = POINTS[right];
singleTriangle.Add((
(new Vector3(leftX, leftY, leftZ) + new Vector3(rightX, rightY, rightZ)) / 2 - new Vector3(0.5f, 0.5f, 0.5f)
) * resolution + coords);
}
foreach (Vector3 pt in singleTriangle) {
points.Add(pt);
tris.Add(points.Count - 1);
}
}
}
}
}
samples.Clear();
return (points.ToArray(), tris.ToArray());
}
}
[RequireComponent(typeof(MeshFilter))]
public class MeshGenerator : MonoBehaviour {
Mesh mesh;
Vector3[] verts;
int[] tris;
int r = 1;
int R = 2;
void Awake() {
mesh = GetComponent<MeshFilter>().mesh;
}
float OnSurface(float x, float y, float z) {
return z *z + (float) Math.Pow(R - Math.Sqrt(x * x + y * y), 2) - r * r;
// return z * z + y * y + x * x - R * R;
// return x * x + y * y + z * z + (float) (Math.Sin(4 * x) + Math.Sin(4 * y) + Math.Sin(4 * z)) - R;
}
void Start() {
// GenerateMeshData();
// CreateMesh();
MarchingCubes generator = new MarchingCubes(OnSurface, 6f, 0.05f);
mesh.Clear();
mesh.indexFormat = UnityEngine.Rendering.IndexFormat.UInt32;
(verts, tris) = generator.GetMesh();
mesh.vertices = verts;
mesh.triangles = tris;
mesh.RecalculateNormals();
}
void GenerateMeshData() {
verts = new Vector3[]{new Vector3(0, 0, 0), new Vector3(0, 0, 1), new Vector3(1, 0, 0)};
tris = new int[]{0, 1, 2};
}
void CreateMesh() {
mesh.Clear();
mesh.vertices = verts;
mesh.triangles = tris;
mesh.RecalculateNormals();
}
}