C# Create a mesh with a hole inside

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
    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;
    }
4 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?

I’m doing the following:

  1. loop through all the faces,
  2. take their respective center points,
  3. check if they collide with any of the inner curves’ bounding boxes,
  4. if they do, perform a more through, but also more expensive inclusion detection using the winding number
  5. if inclusion detected, remove the face from the mesh
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 + 1].Min.X && center.X < bBoxes[j + 1].Max.X && center.Y > bBoxes[j + 1].Min.Y && center.Y < bBoxes[j + 1].Max.Y)
      if (PointInPolygon(center, polylineList[j + 1].ToArray()))
        faceIndex.Add(i);
}

It all seems robust, but I’d like to avoid checking every single face against every single inner curve. For 10 000 curves it takes a significant amount of time…

https://www.cs.cmu.edu/~quake/triangle.html

Triangle.Net is good resource but C# version is buggy.

1 Like

Why you do not construct RTree for skipping checking with all bounding boxes ?

1 Like

I’m not that experienced with RTrees. How would you approach it? This method seems like a good starting point, but I could use an example implementation to understand the concept a bit better.

You will need delete unnecessary code but logic is this:

    public HashSet<long> CollisionPairs;
    public HashSet<int> CollisionPairsFE;
    public Joints plateJoints;
    public int[] CollisionTypes;


    public Dictionary<long, List<int>> componentGroups;
    public HashSet<long> componentGroupsId;
    public Joints GetNeighbours(int[] types) {

        CollisionPairs = new HashSet<long>();
        CollisionPairsFE = new HashSet<int>();
        plateJoints = new Joints();
        CollisionTypes = types;
        componentGroups = new Dictionary<long, List<int>>();
        componentGroupsId = new HashSet<long>();

     

        //Get BoundingBoxes
        Dictionary<int, BoundingBox> boundingBoxes = new Dictionary<int, BoundingBox>();

        int c = 0;
        foreach (var kv in sortedPlates) {
            var bbox = kv.Value.bbox;
            bbox.Inflate(1.1);
            boundingBoxes.Add(c, bbox);
            //Rhino.RhinoDoc.ActiveDoc.Objects.AddBrep(bbox.ToBrep());
            c++;
        }


        //Create Search
        RTree RTreeMeshBBoxes = new RTree();

        //Insert BoundingBoxes to RTree
        for (int i = 0; i < boundingBoxes.Count; i++) {
            RTreeMeshBBoxes.Insert(boundingBoxes[i], i);
        }

        //Perform search
        for (int i = 0; i < RTreeMeshBBoxes.Count - 1; i++) {
            int? index = i;//Declare nullable values
            bool ts = RTreeMeshBBoxes.Search(boundingBoxes[i], SearchCallback, index);
        }


        //Sort Closest line segments
        //Rhino.RhinoApp.WriteLine("Sort Components" + plateJoints.Count.ToString() +" "+ plateJoints[0].ToString());
        //Iterate over componets and mark not needed joints
        foreach (var group in componentGroups) {
            //Rhino.RhinoApp.WriteLine("Groups");
            double[] distances = new double[group.Value.Count];
            int[] id = Enumerable.Range(0, group.Value.Count).ToArray();
            for (int i = 0; i < distances.Length; i++)
                distances[i] = plateJoints[group.Value[i]].jointArea[0].Length;
            Array.Sort(distances, id);

            for (int i = 1; i < distances.Length; i++) {
                //Rhino.RhinoApp.WriteLine("Set To -1");
                plateJoints[group.Value[id[i]]].id = -1;
            }
        }

        //Rhino.RhinoApp.WriteLine("Before"+plateJoints.Count.ToString());
        plateJoints= plateJoints.Duplicate();
        //Rhino.RhinoApp.WriteLine("After" + plateJoints.Count.ToString());


        return plateJoints;
    }



  

    public void SearchCallback(object sender, RTreeEventArgs e) {

        //Rhino.RhinoApp.WriteLine("RTee");

        //Sort the index
        int index1 = (int)(e.Tag as int?);
        int a = index1 < e.Id ? index1 : e.Id;
        int b = index1 < e.Id ? e.Id : index1;

        #region Check indexing
        //Check1: equal elements
        if (a == b) return;
        //Rhino.RhinoApp.WriteLine("RTee not the same element " + a.ToString() + " " + b.ToString() );

        //Check2: pair is visited
        long longKey = (((long)a) << 32) + b;//long key for faceID
  
        //Rhino.RhinoApp.WriteLine(a.ToString() + " " + b.ToString() + " " + longKey.ToString() );
        bool newPair = CollisionPairs.Add(longKey);
        if (!newPair) return;
        //Rhino.RhinoApp.WriteLine("RTee collision pair exists");

        //Check3: belongs to the same component
        int[] AkeyID = this[a].KeyToIntArray();
        int[] BkeyID = this[b].KeyToIntArray();
        string Akey = string.Join(";", AkeyID);
        string Bkey = string.Join(";", BkeyID);
        if (AkeyID[0] == BkeyID[0]) return;
        #endregion

        string key = a.ToString() + ";" + b.ToString() + " " + this[a].key + " " + this[b].key;
        long longKeyGroups = (((long)AkeyID[0]) << 32) + BkeyID[0];//long key for Groups

     

        foreach (int type in CollisionTypes) {
            //Rhino.RhinoApp.WriteLine(((CollisionType)type).ToString() + " " + type.ToString());
            switch ((CollisionType)type) {
                case (CollisionType.LineLine):
                    GetPlateJointLineLine(this[Akey], this[Bkey], longKeyGroups);
                    break;

                case (CollisionType.FaceFace):
                    GetPlateJointFaceFace(this[Akey], this[Bkey]);
                    break;

                case (CollisionType.PlanePolygon):
                    
                    GetPlateJointPlanePolygon(this[Akey], this[Bkey]);
                    //Rhino.RhinoApp.WriteLine("PlanePolygon");
                    break;
            }
        }






        //Rhino.RhinoApp.WriteLine(plateJoints.Count.ToString());

    }
1 Like

well… you did connect outer curve not first… anyway, yes, it had many problems…



mesh_delanuay_holes_V1.1.gh (17.3 KB)

half second for 10k curves :laughing:

Rhino.Geometry.Mesh.Faces.DeleteFaces method accepts a list of indexes… Huge difference!

But…
To have my code working, the curves must not intersect/touch each-other; fixed by uniformly flipping them before the offset step.

And unluckily your “PointInPolygon” is not reliable as much as PolylineCurve.Contains; Rounding pt coordinates seems to be not helping, that step can/should be skipped.
I tried to recreate the winding method but gives me the same result…
Let us know if you find a way to fix that!

2 Likes

Awesome @maje90!

I’ve tested it on a few examples and it works super fast! Compiled version with the winding number method is about 400 milliseconds for 16k curves. Sweet.

I can’t see any issues with the winding number approach on my machine. Were you able to identify any reproducible edge cases on your end?

There is one odd thing about your algorithm - it ignores inner curves if they are a triangle. In the following examples only the offset value was changed, so that a curve which initially was a quad turned into a triangle:


Also, could you please explain how the algorithm works? I don’t really understand what you are doing here. Would appreciate some comments in the code:

int[] cis = new int[total];
int[] AAs = new int[total];
int[] BBs = new int[total];
int[] CCs = new int[total];
int[] DDs = new int[total];

int start = 0;
for (int k = 0; k < curvesCount; k++)
{
    Parallel.For(0, pointsCount[k], (i) => {
        int w = start + i;
        cis[w] = k;
        if (i == 0)
        {
            AAs[w] = w;
            BBs[w] = w;
            CCs[w] = w + 1;
            DDs[w] = start + pointsCount[k] - 1;
        }
        else if (i == 1)
        {
            AAs[w] = w;
            BBs[w] = w;
            CCs[w] = w + 1;
            DDs[w] = start + pointsCount[k];
        }
        else if (i == pointsCount[k] - 2)
        {
            AAs[w] = start - 1;
            BBs[w] = w - 1;
            CCs[w] = w;
            DDs[w] = w;
        }
        else if (i == pointsCount[k] - 1)
        {
            AAs[w] = start;
            BBs[w] = w - 1;
            CCs[w] = w;
            DDs[w] = w;
        }
        else
        {
            AAs[w] = start - 1;
            BBs[w] = w - 1;
            CCs[w] = w + 1;
            DDs[w] = start + pointsCount[k];
        }
    });
    start += pointsCount[k];
}

mesh_delanuay_holes_V1.2.gh (40.9 KB)

Pretty much is still this ^ … but i’ll try to re-explain properly.


My method uses 4 indices to determine “delectable” edges.
Look at vertex 85.
Edges to be removed, which are connected to 85, can only have the other end being from 74 to 83, or from 87 to 96.
(Still, edges “selected” this way need to be checked with an inclusion test, the curve might have a concave shape…)
Every other edge doesn’t even need to be checked.

Particular cases need to be implemented near the start and end of the curve.
For example, from the 74 position AAs<->BBs range is null.

(“cis” is the index of the curve that generated each vertex…)

With simple triangular holes both AA-BB and CC-DD were nulled, and so no edge is “selectable”.
Also, on a triangular hole there is no edge at all inside! :sweat_smile:
Totally a flaw in my algorithm logic.


This method works on edges. Easier than on faces, but probably also slower.

Working per-face the same AA-BB-CC-DD lists can be used, but for each face 2 vertex indexes must be checked, and if not “delectable” is found, even the third vertex index need to be checked.

Also, on triangular meshes there are X faces and 1.5*X edges … if i’m not wrong… so less inclusion tests at the end.


Edit:
found how i was using wrong your PointInPolygon method, my mistake, it works fine.
mesh_delanuay_holes_V1.2.gh (15.3 KB)
It’s almost completely same as last version. 3 character fixed.
Triangles still unpredictable, face deletion instead of edge deletion still to implement.

1 Like


mesh_delanuay_holes_V1.3.gh (20.5 KB)

Yep, working with faces is faster.

More integer comparison tests but less PointInPolygon inclusion tests.
Also now it supports triangular holes.

I’ll leave to you further testing.

Anyway… 100k vertices 140k faces mesh in less than a second! :laughing:
It came out a nice tool for meshing pre-nested curves.

4 Likes