Hi, does anyone has any idea how to implement isosurface in grasshopper c#?
I know there is a component in Milipede but i need my own
There is not great difficulties implementing it. The difficulty is to have an optimized one. You didn’t search a lot on internet, classical reference is
http://paulbourke.net/geometry/polygonise/source1.c
private void RunScript(List<double> values, Box box, int nx, int ny, int nz, double iso, ref object M, ref object P)
{
List<Point3d> boxPoints = new List<Point3d>();
double dx = (box.X.Max - box.X.Min) / (nx - 1);
double dy = (box.Y.Max - box.Y.Min) / (ny - 1);
double dz = (box.Z.Max - box.Z.Min) / (nz - 1);
for (int z = 0; z < nz; z++)
{
for (int y = 0; y < ny; y++)
{
for (int x = 0; x < nx; x++)
{
boxPoints.Add(new Point3d(box.X.Min + dx * (double) x, box.Y.Min + dy * (double) y, box.Z.Min + dz * (double) z));
}
}
}
Rhino.Geometry.Mesh mesh = new Rhino.Geometry.Mesh();
int v0, v1, v2, v3, v4, v5, v6, v7;
for (int x = 0; x < (nx - 1); x++)
{
for (int y = 0; y < (ny - 1); y++)
{
for (int z = 0; z < (nz - 1); z++)
{
v0 = x + y * nx + z * nx * ny;
v1 = v0 + 1;
v2 = v1 + nx;
v3 = v0 + nx;
v4 = v0 + nx * ny;
v5 = v1 + nx * ny;
v6 = v2 + nx * ny;
v7 = v3 + nx * ny;
PolygoniseTri(values, iso, mesh, boxPoints, v0, v2, v3, v7);
PolygoniseTri(values, iso, mesh, boxPoints, v0, v2, v6, v7);
PolygoniseTri(values, iso, mesh, boxPoints, v0, v4, v6, v7);
PolygoniseTri(values, iso, mesh, boxPoints, v0, v6, v1, v2);
PolygoniseTri(values, iso, mesh, boxPoints, v0, v6, v1, v4);
PolygoniseTri(values, iso, mesh, boxPoints, v5, v6, v1, v4);
}
}
}
mesh.Vertices.CombineIdentical(true, true);
mesh.UnifyNormals();
mesh.Normals.ComputeNormals();
M = mesh;
P = boxPoints;
}
// <Custom additional code>
/*
Reference : Paul Bourke
Internet : http://paulbourke.net/geometry/polygonise/
Polygonise a tetrahedron given its vertices within a cube
This is an alternative algorithm to polygonisegrid.
It results in a smoother surface but more triangular facets.
+ 0
/|\
/ | \
/ | \
/ | \
/ | \
/ | \
+-------------+ 1
3 \ | /
\ | /
\ | /
\ | /
\ | /
\|/
+ 2
It's main purpose is still to polygonise a gridded dataset and
would normally be called 6 times, one for each tetrahedron making
up the grid cell.
Given the grid labelling as in PolygniseGrid one would call
PolygoniseTri(grid,iso,triangles,0,2,3,7);
PolygoniseTri(grid,iso,triangles,0,2,6,7);
PolygoniseTri(grid,iso,triangles,0,4,6,7);
PolygoniseTri(grid,iso,triangles,0,6,1,2);
PolygoniseTri(grid,iso,triangles,0,6,1,4);
PolygoniseTri(grid,iso,triangles,5,6,1,4);
*/
public void PolygoniseTri(List<double> values, double iso, Mesh mesh, List<Point3d> boxPoints, int v0, int v1, int v2, int v3)
{
/*
Determine which of the 16 cases we have given which vertices
are above or below the isosurface
*/
int vertCount = mesh.Vertices.Count;
int triindex = 0;
if (values[v0] < iso) triindex |= 1;
if (values[v1] < iso) triindex |= 2;
if (values[v2] < iso) triindex |= 4;
if (values[v3] < iso) triindex |= 8;
/* Form the vertices of the triangles for each case */
switch (triindex) {
case 0x00:
case 0x0F:
break;
case 0x0E:
case 0x01:
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v0], boxPoints[v1], values[v0], values[v1]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v0], boxPoints[v2], values[v0], values[v2]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v0], boxPoints[v3], values[v0], values[v3]));
mesh.Faces.AddFace(vertCount, vertCount + 1, vertCount + 2);
break;
case 0x0D:
case 0x02:
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v1], boxPoints[v0], values[v1], values[v0]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v1], boxPoints[v3], values[v1], values[v3]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v1], boxPoints[v2], values[v1], values[v2]));
mesh.Faces.AddFace(vertCount, vertCount + 1, vertCount + 2);
break;
case 0x0C:
case 0x03:
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v0], boxPoints[v3], values[v0], values[v3]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v0], boxPoints[v2], values[v0], values[v2]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v1], boxPoints[v3], values[v1], values[v3]));
// mesh.Faces.AddFace(vertCount, vertCount + 1, vertCount + 2);
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v1], boxPoints[v2], values[v1], values[v2]));
// mesh.Faces.AddFace(vertCount + 2, vertCount + 3, vertCount + 1);
mesh.Faces.AddFace(vertCount, vertCount + 1, vertCount + 3, vertCount + 2);
break;
case 0x0B:
case 0x04:
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v2], boxPoints[v0], values[v2], values[v0]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v2], boxPoints[v1], values[v2], values[v1]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v2], boxPoints[v3], values[v2], values[v3]));
mesh.Faces.AddFace(vertCount, vertCount + 1, vertCount + 2);
break;
case 0x0A:
case 0x05:
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v0], boxPoints[v1], values[v0], values[v1]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v2], boxPoints[v3], values[v2], values[v3]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v0], boxPoints[v3], values[v0], values[v3]));
//mesh.Faces.AddFace(vertCount, vertCount + 1, vertCount + 2);
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v1], boxPoints[v2], values[v1], values[v2]));
//mesh.Faces.AddFace(vertCount, vertCount + 3, vertCount + 1);
mesh.Faces.AddFace(vertCount, vertCount + 2, vertCount + 1, vertCount + 3);
break;
case 0x09:
case 0x06:
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v0], boxPoints[v1], values[v0], values[v1]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v1], boxPoints[v3], values[v1], values[v3]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v2], boxPoints[v3], values[v2], values[v3]));
//mesh.Faces.AddFace(vertCount, vertCount + 1, vertCount + 2);
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v0], boxPoints[v2], values[v0], values[v2]));
//mesh.Faces.AddFace(vertCount, vertCount + 3, vertCount + 2);
mesh.Faces.AddFace(vertCount, vertCount + 1, vertCount + 2, vertCount + 3);
break;
case 0x07:
case 0x08:
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v3], boxPoints[v0], values[v3], values[v0]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v3], boxPoints[v2], values[v3], values[v2]));
mesh.Vertices.Add(VertexInterp(iso, boxPoints[v3], boxPoints[v1], values[v3], values[v1]));
mesh.Faces.AddFace(vertCount, vertCount + 1, vertCount + 2);
break;
}
}
public Point3d VertexInterp(double iso, Point3d p0, Point3d p1, double v0, double v1)
{
return new Point3d(p0 + (p1 - p0) * (iso - v0) / (v1 - v0));
}
6 Likes