C# Create a mesh with a hole inside

Hey everyone,

I’m trying to create a mesh consisting of an ‘outside’ and a bunch of ‘inside’ curves, such that it would keep the holes defined by the ‘inside’ ones. Essentially a mesh counterpart of the Boundary Surfaces command. Obviously I can just create a Brep and convert it to a mesh:

Brep[] brep = Brep.CreatePlanarBreps(iCurves, 0.1);
Mesh[] mesh = Mesh.CreateFromBrep(brep[0], MeshingParameters.FastRenderMesh);

But is there a way of directly creating a mesh, which would behave in the exact same way?

I tried tessellating the mesh from a list of points, but it creates the internal faces as well:

List<Point3d> points = new List<Point3d>();
List<Polyline> polylineList = new List<Polyline>();

foreach (var crv in iCurves)
{
  var polylineCrv = crv.ToPolyline(-1, -1, 10, 1, 0, 0.1, 1, 10, true).ToPolyline();

  points.AddRange(polylineCrv);
  polylineList.Add(polylineCrv);
}

A = Mesh.CreateFromTessellation(points, polylineList, Plane.WorldXY, false);

I give it a try and didn’t manage to make it work. If you want to stick to mesh you must search for each face if it is inside or not. And delete faces according to the result.

To me this seems like a bug as documentation is clearly saying there are holes

I have developed some workaround in C# in order to be able to use Rhino3dmIO and Clipper. It is based on Winding number.

I am using this in NGon plugin for quite some time. No bugs works very well. You just need to give largest curve and holes as input.

1 Like

Awesome! Thanks @Petras_Vestartas!

Here is the code for anyone who is interested:

        iCurves = new List<Curve>();
        oCurves = new List<Curve>();
        List<Mesh> meshes = new List<Mesh>();

        Polyline pCrv = oCurves[0].ToPolyline(-1, -1, 0.25, 1, 0, 0.1, 1, 10, true).ToPolyline();
        if (!pCrv.IsClosed) pCrv.Add(pCrv[0]);

        Mesh mesh = Mesh.CreatePatch(pCrv, 0.25, null, iCurves, null, null, false, 32);
        meshes.Add(mesh);

In fact this is the code (general case: classic polygon [with holes] triangulation).

earclipper-master.zip (842.6 KB)

Hmmm,
having done some testing, it seems that the meshing part doesn’t always generate clean results. It looks like a precision issue. Depending on the distance between the internal curves the resulting mesh can either be correct or broken. See attached screenshots and .gh file.

Any ideas @Petras_Vestartas?


Patch.gh (8.2 KB)

Thanks @PeterFotiadis. Interesting approach, I haven’t read this paper before. For my particular use case, I need something blazing fast, though. There will potentially be thousands of these inclusions and here it seems, that the complexity increases exponentially with the number of holes…

Thanks @laurent_delrieu,
I’ve been thinking about this as well.

  1. Create mesh from tessellation
  2. Iterate through the faces and check whether all their vertices are contained by the inside curves
  3. Delete these faces

This looks like a lot of work, which could be avoided if these faces weren’t generated in the first place…

For who?

You mean it’s a lot of work for the person who will code that script?
Or you mean that is unnecessary task done by the CPU so you want to avoid that because is un-optimized code?

Mostly the latter. As mentioned earlier, there will potentially be thousands of these cases and I’d like the code to execute as fast as possible.

Also, traversing a native Rhino mesh is a bit cumbersome. Is there an intuitive way of selecting only the faces, which are on the inside of an edge loop? In a half-edge mesh structure like Plankton it is a trivial task. Not sure how to best approach it with vanilla Rhinocommon, though.

Looking at the code above ^ , each curve become a list of points.
Then you chain together all the lists into a big list and do the delanuay triangulation.
Edges that start or end to a point in an interior curve must have the other end that is not any other point of that same curve.
For example, an interior curve divided into 20 points, final vertex index from 200 to 219.
Vertex 210 can connect only to 209 and 211 but nothing else in the 200-219 range. If that happens, that edge is to be removed.
This should be a decently fast iterating, as it is working only with simple integers, no area/distance/intersections needed.

But yes, the holes initially are calculated anyway… not the best approach.

And here it is.
As said, this method makes a check and remove edges if they start and ends from the same curve… but! I didn’t tought about convex shapes! So, before removing the edge, its midpoint is checked to be outside of the main shape or inside an internal shape (according to edge “birthplace”).

mesh_delanuay_holes_V1.gh (18.5 KB)
The script takes a list of curves as input. First curve must be the outer one.
The script doesn’t check for intersections and works on XY plane.

I’ve left inside 6 alternative of the same script: for loop, parallel loop (multithreaded), Curve or Polyline inclusion test and outputting points or not.
One of the two in the green box should be the best… probably the Curve one, for precision.

It seems to be 2x slower than original Delanuay component… not that bad, imho…
More tests needed.

Hope it helps!

Code of the Parallel.For, Curve inclusion test version:

using System.Linq;
using System.Threading.Tasks;

    private void RunScript(List<Curve> Curves, ref object Mesh)
      {
    double t = this.RhinoDocument.ModelAbsoluteTolerance;

    int cc = Curves.Count;
    Rhino.Geometry.PolylineCurve[] polylines = new Rhino.Geometry.PolylineCurve[cc];
    Rhino.Geometry.Point3d[][] pts = new Rhino.Geometry.Point3d[cc][];
    int[] pc = new int[cc]; // Point count for each polyline

    // Transforming each curve into a polyline, extracting the point array and counting them
    //for(int i = 0;i < cc;i++){  // --- FOR --- //
    Parallel.For(0, cc, (i) => { // --- PAR --- //
      polylines[i] = Curves[i].ToPolyline(-1, -1, 0.2, 1, 100, 0.1, 0.1, 1, false);
      pts[i] = polylines[i].ToPolyline().ToArray();
      //polylines[i] = ClosePolylineCurve(polylines[i]); // Closing PolylineCurve for inclusion test
      pc[i] = pts[i].Length;
      }); // --- PAR --- //
    //} // --- FOR --- //

    int total = pc.Sum(); // Total amount of points

    Rhino.Geometry.Point3d[] points = pts.SelectMany(x => x).ToArray();
    int[] cis = new int[total];
    int[] prevs = new int[total];
    int[] nexts = new int[total];
    int[] starts = new int[total];
    int[] ends = new int[total];

    int start = 0;
    for(int k = 0;k < cc;k++){
      //for(int i = 0;i < pc[k];i++){  // --- FOR --- //
      Parallel.For(0, pc[k], (i) => { // --- PAR --- //
        int w = start + i;
        int prev;
        int next;
        Neighbors(w, i, pc[k], out prev, out next);
        cis[w] = k;
        prevs[w] = prev;
        nexts[w] = next;
        starts[w] = start;
        ends[w] = start + pc[k];
        }); // --- PAR --- //
      //} // --- FOR --- //
      start += pc[k];
    }

    Rhino.Geometry.Mesh mesh = Rhino.Geometry.Mesh.CreateFromTessellation(points, null, Plane.WorldXY, false);

    Rhino.Geometry.Collections.MeshTopologyEdgeList topedges = mesh.TopologyEdges;
    int[] faces = Enumerable.Repeat(-1, topedges.Count * 2).ToArray();

    //for(int i = 0;i < topedges.Count;i++){  // --- FOR --- //
    var result = Parallel.For(0, topedges.Count, (i) => { // --- PAR --- //

      Rhino.IndexPair pair = topedges.GetTopologyVertices(i);

      if(pair.J >= starts[pair.I] & pair.J <= ends[pair.I]){
        if(!(pair.I == prevs[pair.J] || pair.I == nexts[pair.J])){

          if(cis[pair.I] == 0){
            if(Curves[0].Contains((points[pair.I] + points[pair.J]) / 2, Plane.WorldXY, t).Equals(PointContainment.Outside)){
              //if(polylines[0].Contains((points[pair.I] + points[pair.J]) / 2, Plane.WorldXY, t).Equals(PointContainment.Outside)){
              int asd = 0;
              foreach(int f in topedges.GetConnectedFaces(i)){
                faces[i * 2 + asd] = f;
                asd++;
              }
            }
          }
          else{
            if(Curves[cis[pair.I]].Contains((points[pair.I] + points[pair.J]) / 2, Plane.WorldXY, t).Equals(PointContainment.Inside)){
              //if(polylines[cis[pair.I]].Contains((points[pair.I] + points[pair.J]) / 2, Plane.WorldXY, t).Equals(PointContainment.Inside)){
              int asd = 0;
              foreach(int f in topedges.GetConnectedFaces(i)){
                faces[i * 2 + asd] = f;
                asd++;
              }
            }
          }
        }
      }
      }); // --- PAR --- //
    //} // --- FOR --- //

    faces = faces.Distinct().OrderByDescending(f => f).ToArray();

    foreach(int f in faces){
      if(f > -1){
        mesh.Faces.RemoveAt(f);
      }
    }

    Mesh = mesh;
    //  Points = points;
  }

  // <Custom additional code> 
  public void Neighbors(int val, int index, int loopsize, out int previous, out int next){
    loopsize += -1;
    if(index == 0){      previous = val + loopsize;    }
    else{      previous = val - 1;    }
    if(index == loopsize){ next = val - loopsize;}
    else{next = val + 1;}
  }
  public Rhino.Geometry.PolylineCurve ClosePolylineCurve(Rhino.Geometry.PolylineCurve p){
    Rhino.Geometry.Polyline poly = p.ToPolyline();
    poly.Add(p.PointAtStart);
    p = poly.ToPolylineCurve();
    return p;
  }
2 Likes

Currently this is the method that yields the best results. I have 2 C#'s (strictly internal) that yieid results (for AEC related cases: trusses, tension membrane flat layout (for K2 relaxation etc etc) and the likes) in a couple of milliseconds (using various I9 K). But … well … I confess that I have very lttle interest about accuracy and/or following the boundary in an exact way. See what I mean (triangulation with inner plus points for K2):

Of course for non AEC matters there’s the CUDA way … far and away from what GH/R can do…

I would strongly advise to “implement” the code above in some C# of yours.

Thanks @PeterFotiadis!
Following your suggestion, I’ve tried out the EarClipper library and the results aren’t that encouraging. The algorithm sometimes struggles to produce a clean mesh with even very simple shapes, on more complex ones it does a good job, but takes relatively long to compute the results. I clocked it with a System.Diagnostics.Stopwatch to better understand what is going on. The majority of the time is spent on the triangulation process, which with the built-in Delaunay Mesh component doesn’t even trigger the Profiler (<5ms).

Mind you, this was tested with a compiled plugin just to make sure that I’m getting the most performance out of it.

Thanks @maje90!

That’s some thorough testing, you’ve done there! While tracking the EarClipper’s performance, I had an epiphany moment. We already have a list of points for both the external and internal curves, we also know, that we want our triangles to have at least one vertex attached to the outer curve. Rather than traversing an already triangulated mesh, to remove the unwanted faces, we could perform this check during the triangulation.

MeshVertexList remains the same, we just add less entries to the MeshFaceList.

I’ll try to find an open-source triangulation library and cook up a test implementation to validate this assumption.

I thought a bit about that…
But what would work (maybe) only if you always have convex shaped holes… with concave holes you’ll need again some additional check.
With simple convex shapes you can use angles to limit the triangulation algorithm, but it’s a lot of computation anyway.

Anyway, let us know if you find better solutions.
This is interesting!

You’re right, @maje90, I completely missed the concave hole issue.

Anyway, I’ve made some progress and learned 2 things:

  1. Apparently there is an issue with Point3d precision, which can lead to unexpected results while performing operations on meshes. @Krzysztof described it in more detail in the following post:
    Mesh Inclusion Error - #19 by Krzysztof
    After rounding Point3d coordinates to 4 decimal places mesh triangulation works as expected.

  2. Rather than checking individual vertices for intersection with internal edges, I’m iterating over mesh faces and perform the inclusion test on their center points. This operation is split into 2 phases:
    → coarse, where we check against a bounding box of each curve
    → detailed, more expensive Curve.Contains method is used for precision

It works quite nicely, but still, we spend the majority of the time removing unwanted faces. Comparison to the most promising algorithm from @maje90’s previous post is a bit unfair, since his approach results in a corrupt mesh. I suspect this is related to the mesh precision issue mentioned above.

Questions:

  1. Is there a smarter way of changing precision of a List<Point3d> than iterating over the entire list and changing the individual coordinate of each point?

foreach(var pt in pCrv)
{
double x = Math.Round(pt.X, 4);
double y = Math.Round(pt.Y, 4);
double z = Math.Round(pt.Z, 4);
points.Add(new Point3d(x, y, z));
}

  1. How can I iterate over a list of faces in a thread safe way? The following code converted to a Parallel.For loop breaks the mesh:
List<BoundingBox> bBoxes = new List<BoundingBox>();
foreach(var crv in iCurves)
{
  bBoxes.Add(crv.GetBoundingBox(false));
}

Point3d center;
List<int> faceIndex = new List<int>();
for (int i = 0; i < mesh.Faces.Count; i++)
{
  center = mesh.Faces.GetFaceCenter(i);
  for (int j = 0; j < iCurves.Count; j++)
  {
    if (center.X > bBoxes[j].Min.X && center.X < bBoxes[j].Max.X && center.Y > bBoxes[j].Min.Y && center.Y < bBoxes[j].Max.Y)
    {
      if (iCurves[j].Contains(center, Plane.WorldXY, 0.01) == PointContainment.Inside)
        faceIndex.Add(i);
    }
  }
}

mesh.Faces.DeleteFaces(faceIndex);[meshing.gh|attachment]

meshing.gh (14.2 KB)

Maybe you can use clipper for point inclusion it will be probably faster than parallel Curve contains method from RhinoCommon. It is not so difficult to set it up for C#.

You can check point inclusion for polylines with holes:

public static int PointInPolygon(IntPoint pt, Path path)
      {
        //returns 0 if false, +1 if true, -1 if pt ON polygon boundary
        //See "The Point in Polygon Problem for Arbitrary Polygons" by Hormann & Agathos
        //http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.88.5498&rep=rep1&type=pdf
        int result = 0, cnt = path.Count;
        if (cnt < 3) return 0;
        IntPoint ip = path[0];
        for (int i = 1; i <= cnt; ++i)
        {
          IntPoint ipNext = (i == cnt ? path[0] : path[i]);
          if (ipNext.Y == pt.Y)
          {
            if ((ipNext.X == pt.X) || (ip.Y == pt.Y &&
              ((ipNext.X > pt.X) == (ip.X < pt.X)))) return -1;
          }
          if ((ip.Y < pt.Y) != (ipNext.Y < pt.Y))
          {
            if (ip.X >= pt.X)
            {
              if (ipNext.X > pt.X) result = 1 - result;
              else
              {
                double d = (double)(ip.X - pt.X) * (ipNext.Y - pt.Y) -
                  (double)(ipNext.X - pt.X) * (ip.Y - pt.Y);
                if (d == 0) return -1;
                else if ((d > 0) == (ipNext.Y > ip.Y)) result = 1 - result;
              }
            }
            else
            {
              if (ipNext.X > pt.X)
              {
                double d = (double)(ip.X - pt.X) * (ipNext.Y - pt.Y) -
                  (double)(ipNext.X - pt.X) * (ip.Y - pt.Y);
                if (d == 0) return -1;
                else if ((d > 0) == (ipNext.Y > ip.Y)) result = 1 - result;
              }
            }
          }
          ip = ipNext;
        }
        return result;
      }
1 Like

Thanks @Petras_Vestartas!

In the meantime I have implemented the winding number approach from this paper as suggested by @laurent_delrieu. It is almost 10 times faster than the Curve.Contains method.

    public static bool PointInPolygon(Point3d p, Point3d[] polylineArray)
    {
        int n = polylineArray.Length - 1;
        int windingNumber = 0;    // the winding number counter

        // loop through all edges of the polygon
        for (int i = 0; i < n; i++)
        {   // edge from V[i] to V[i+1]
            if (polylineArray[i].Y <= p.Y)
            {         // start y <= P.y
                if (polylineArray[i + 1].Y > p.Y)      // an upward crossing
                    if (IsLeft(polylineArray[i], polylineArray[i + 1], p) > 0)  // P left of edge
                        ++windingNumber;            // have a valid up intersect
            }
            else
            {                       // start y > P.y (no test needed)
                if (polylineArray[i + 1].Y <= p.Y)     // a downward crossing
                    if (IsLeft(polylineArray[i], polylineArray[i + 1], p) < 0)  // P right of edge
                        --windingNumber;            // have a valid down intersect
            }
        }
        if (windingNumber != 0)
            return true;
        else
            return false;

    }

    private static int IsLeft(Point3d p0, Point3d p1, Point3d p2)
    {
        double calc = ((p1.X - p0.X) * (p2.Y - p0.Y)
          - (p2.X - p0.X) * (p1.Y - p0.Y));
        if (calc > 0)
            return 1;
        else if (calc < 0)
            return -1;
        else
            return 0;
    }
5 Likes

So you need 1st to loop through points and check if they are inside bigger curve and then 2nd time to loop if point is not inside inner curve?