C# Create a mesh with a hole inside

[mrhe (Mariusz Hermansdorfer)]

@Petras_Vestartas

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…

[Antonio_Cersosimo]

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

[Petras_Vestartas]

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

[Petras_Vestartas]

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

[mrhe (Mariusz Hermansdorfer)]

@Petras_Vestartas

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

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.

https://developer.rhino3d.com/api/RhinoCommon/html/M_Rhino_Geometry_RTree_CreateMeshFaceTree.htm

[Petras_Vestartas]

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());

    }

[maje90 (Riccardo Majewski)]

@mrhe

his approach results in a corrupt mesh

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

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!

[mrhe (Mariusz Hermansdorfer)]

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)

[maje90 (Riccardo Majewski)]

@maje90

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.

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

[maje90 (Riccardo Majewski)]

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!
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.

[maje90 (Riccardo Majewski)]

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!
It came out a nice tool for meshing pre-nested curves.

[mrhe (Mariusz Hermansdorfer)]

Thanks a million for your awesome contribution @maje90!

Didn’t have too much time to do proper testing, but on first impression it works very well. To make things run even faster, I’m using Bowerbird’s Offset component. Under the hood, it leverages the Clipper library to clean up geometry, unify curve directions and create reliable offsets. All this much faster than the built-in Offset Curve component.

Clipper for Grasshopper does the same thing:

arendvw/clipper

oberbichler/Bowerbird