C# isosurface

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));
  }
3 Likes