Non-Self-Intersecting closed PolyLine from planar set of points

Hi everybody

I have not found in rhinocommon method to create polyline from an arbitrary set of planar points in a way that polyline is not self intersecting.
There are also some other requirements while creating polyline, such as find shortest-length polyline or longest-lenght polyline, etc.
Any suggestions?

Thanks

Sample file with points: PolyLine01.3dm (186.0 KB)

Hi @RadovanG,

i don’t think there is a RhinoCommon function to get there in one step. Not sure if this helps, but if you find a circle radius by increasing it gradually, and then boolean union all circles using all points as centers you might find a curve which encloses all points. (Blue in my example).

If you then would iterate over all points and exclude the ones which are not within the distance of the circle radius you could find the points to connect. The connection order could be derived from a curve closest point operation and the parameter on the curve. The process (boolean union of all circles) might be computational intensive but it would lead to the green curve in my example file.

Outline.3dm (328.7 KB)

_
c.

Hi Clement

Thanks for suggestion.
I have already found another fast algorithm to create closed polyline from set of planar points.
Actualy I remebered long time ago on one mathematical competition there was problem to be solved that sounded like this:
If we have 3000 points in plane and none 3 of them are colinear then prove that there is way to construct 1000 triangles which do not intersect and do not share any point.
Solution is simple, just construct line which is not parallel to any of lines created by any of two points from available set of 3000 poinnts. Then start to move such line from one"edge" of set to the other and as you move you will pick up point by point. Each time you pick i*3rd point you cosntruct i-th triangle.
So I use similar aproach to construct closed polyline. Here is the code with explanation inside comments:

    public static List<Point3d> CreateClosedPolyLine(List<Point3d> points, double tolerance)
    {
        List<Point3d> pointsCulled = Point3d.CullDuplicates(points, tolerance).ToList();
        //we sort points by defaukt comparer, see Point3d.GreaterThan Operator 
        pointsCulled.Sort();
        if (pointsCulled.Count<3)
        {
            return null;
        }
        //calculate average-point, lets call it centroidPoint
        var centroidPoint = pointsCulled.Aggregate((pointSum, point) => pointSum += point);
        centroidPoint /= pointsCulled.Count;
        //now we  split points in two sets, upper and lower, base on condition if point is above or below centroid point ( comapring Y coord.)
        var upperPoints = pointsCulled.Where(p => p.Y > centroidPoint.Y).ToList();
        var lowerPoints = pointsCulled.Where(p => p.Y <= centroidPoint.Y).ToList();
        // Now we can connect points in lowerPoints set form left to right 
        // and points in upperPoints from right to left 
        // and merge these two result in final closed polyline => 
        //    we will conect most-left point form upperPoints set with most-left point from lowerPoints set 
        //    and most-right point form lowerPoints set with most-right point from upperPoints set.
        // !!! But we have to handle situation when there is more than one most-left and most-right points in boths sets 
        //       (more points with same X but different Y ).
        // So, we have to connect most-left point form upperPoints set (it has lowest Y) with most-left point from lowerPoints set that has highest Y
        //   which means that we have to reverse order of most-left points in lowerPoints set 
        //   and reverse order of most-right point in upperPoints set 

        Point3d lowerSet_MostLeftPoint = lowerPoints.Min();
        var lowerMostLeftPoints = lowerPoints.Where(p=> p.X == lowerSet_MostLeftPoint.X ).Reverse().ToList();
        var lowerTheRestOfPoints = lowerPoints.Where(p => p.X > lowerSet_MostLeftPoint.X).ToList();

        Point3d upperSet_MostRightPoint = upperPoints.Max();
        var upperMostRightPoints = upperPoints.Where(p => p.X == upperSet_MostRightPoint.X).Reverse();
        var upperTheRestOfPoints = upperPoints.Where(p => p.X < upperSet_MostRightPoint.X);

        //Now we merge everything into one list
        List<Point3d> polyLinePoints = new List<Point3d>(lowerMostLeftPoints);
        polyLinePoints.AddRange(lowerTheRestOfPoints);
        polyLinePoints.AddRange(upperMostRightPoints.Reverse());
        polyLinePoints.AddRange(upperTheRestOfPoints.Reverse());
        //and to make closed polyline we have to add one more point that is equal to first point  
        polyLinePoints.Add(lowerMostLeftPoints[0]);

        return polyLinePoints;
    }

But this is just one solution.
I need to find shortest path, shortest nonintersecting closed polyline.
I gueass that I will have to iterate over all posible solutions, and find shortest polygon…

!!! There is bug.
Actauly I can not use imaginary horizontal line sto split sets as describe in previous code, but I have to create line determined with most left point and most right point and then split points in upper and lower set by this line.
Updated code is here:

    public static List<Point3d> CreateClosedPolyLine(List<Point3d> points, double tolerance)
    {
        List<Point3d> pointsCulled = Point3d.CullDuplicates(points, tolerance).ToList();
        //we sort points by default comparer, see Point3d.GreaterThan Operator 
        pointsCulled.Sort();
        if (pointsCulled.Count < 3)
        {
            return null;
        }
        //we determine most left and most right point
        Point3d A = pointsCulled.First(); // most left
        Point3d B = pointsCulled.Last(); // most right
        if (A.X == B.X)
        {
            //points are laying on vertical line
            return pointsCulled;
        }
        // Now we calclulate equation of line determined with these two points
        // Equation is Y= k*X + q, if we now two points then 
        // k = (y2-y1)/(x2-x1)
        // q = ( x2*y1 -x1*y2 )/ (x2-x1)
        double k = (B.Y - A.Y) / (B.X - A.X);
        double q = (B.X * A.Y - A.X * B.Y) / (B.X - A.X);
        //now we  split points in two sets, upper and lower, base on condition if point is above or below line ( comapring Y coord.)
        var upperPoints = pointsCulled.Where(p => p.Y > k * p.X + q).ToList();
        var lowerPoints = pointsCulled.Where(p => p.Y <= k * p.X + q).ToList();
        // Now we can connect points in lowerPoints set form left to right 
        List<Point3d> polyLinePoints = new List<Point3d>();
        polyLinePoints.AddRange(lowerPoints);
        //    and points in upperPoints from right to left 
        upperPoints.Reverse();
        polyLinePoints.AddRange(upperPoints);
        //    and add one more time start point at the end of list to have closed polyline 
        polyLinePoints.Add(lowerPoints[0]);

        return polyLinePoints;
    }

Hi lement

I studied your suggestion and applied basic idea of creating outer polygon.
Actually if I

  1. create outer CONVEX polygon (outer = all points are inside polygon or on corners of polygon) p000
  2. again create outer CONVEX polygon p001 of remaining points
  3. “merge” p000 and p001 into p002
  4. again create outer CONVEX polygon p003 of remaining points
  5. then “merge” p002 and p003 into p004
  6. and repeat repeat step 4) and 5) untill no more points left…
    Resulting polygon will be closed nonintersecting and I guess short in length comparing to all posible solutions.
    PolyLine02.3dm (306.0 KB)

Hi @RadovanG,

the last file posted above above looks identical to the result you would get when using grasshopper ConvexHull component:

I wonder if you could share the code so i could try to convert it into python ?

_
c.

Hi Clement

Here is code for method like ConvexHull.

    public static List<Point3d> CreateConvexHull(List<Point3d> points, double tolerance)
    {
        List<Point3d> pointsCulled = Point3d.CullDuplicates(points, tolerance).ToList();
        //we sort points by defaukt comparer, see Point3d.GreaterThan Operator 
        pointsCulled.Sort();             
        if (pointsCulled.Count < 3)
        {
            return null;
        }
        if (!Point3d.ArePointsCoplanar(pointsCulled, tolerance))
        {
            return null;
        }
        //we determine most left and most right point
        Point3d A = pointsCulled.First(); // most left
        Point3d B = pointsCulled.Last(); // most right
        if (A.X == B.X)
        {
            //points are laying on vertical line
            return pointsCulled;
        }
        // Now we calclulate equation of line determined with points A and B
        // Equation is Y= k*X + q, if we now two points then 
        // k = (y2-y1)/(x2-x1)
        // q = ( x2*y1 -x1*y2 )/ (x2-x1)
        double k = (B.Y - A.Y) / (B.X - A.X);
        double q = (B.X * A.Y - A.X * B.Y) / (B.X - A.X);

        //Create resulting list
        List<Point3d> polyLinePoints = new List<Point3d>();
        //add A point to result
        polyLinePoints.Add(A);
        //and process all points below A-B line
        Vector3d previousVector = new Vector3d(0, -1, 0); ;
        for (int i = 1; i < pointsCulled.Count ; i++)
        {
            double minAngle = Math.PI;
            int foundPointIndex = -1;
            Point3d lastPoint = polyLinePoints[polyLinePoints.Count - 1];
            for (int j = i; j < pointsCulled.Count; j++)
            {
                var px = pointsCulled[j].X;
                var py = pointsCulled[j].Y;                    
                //check that point is below A-B line,
                if (py <= k * px + q || j== pointsCulled.Count-1)
                {
                    Vector3d currVector = pointsCulled[j] - lastPoint;
                    double angle = Vector3d.VectorAngle(previousVector, currVector, Plane.WorldXY);
                    if (angle < minAngle)
                    {
                        minAngle = angle;
                        foundPointIndex = j;
                    }
                }
            }
            //process found whose vector has minimum angle 
            if (foundPointIndex> -1)
            {
                polyLinePoints.Add(pointsCulled[foundPointIndex]);
                i = foundPointIndex;
                previousVector = polyLinePoints[polyLinePoints.Count - 1] - polyLinePoints[polyLinePoints.Count - 2];
            }
        }
        //and process all points above A-B line  (A and B should be already in resulting list polyLinePoints)
        previousVector = polyLinePoints[polyLinePoints.Count - 1] - polyLinePoints[polyLinePoints.Count - 2];
        //we process from right to left
        for (int i = pointsCulled.Count-2; i >=0; i--)
        {
            double minAngle = Math.PI;
            int foundPointIndex = -1;
            Point3d lastPoint = polyLinePoints[polyLinePoints.Count - 1];                
            for (int j = i; j >=0 ; j--)
            {
                var px = pointsCulled[j].X;
                var py = pointsCulled[j].Y;
                //check if point is above A-B line
                if (py > k * px + q || j == 0)
                {
                    Vector3d currVector = pointsCulled[j] - lastPoint;
                    double angle = Vector3d.VectorAngle(previousVector, currVector, Plane.WorldXY);
                    if (angle < minAngle)
                    {
                        minAngle = angle;
                        foundPointIndex = j;
                    }
                }
            }
            //process found point whose vector has minimum angle 
            if (foundPointIndex > -1)
            {
                polyLinePoints.Add(pointsCulled[foundPointIndex]);
                i = foundPointIndex;
                previousVector = polyLinePoints[polyLinePoints.Count - 1] - polyLinePoints[polyLinePoints.Count - 2];
            }
        }
        return polyLinePoints;
    }

Code can be shorten and probably more optimised but for easy readibility I post it like this.
Algorithm is ismple. We find most left point, A, and most right point, B, which will be definitely part of ConvexHull polyline. First we process points below A-B line and later points above A-B line.
Point A is first point, we find next point P by creating vector A-P and checking for vector that has minimum angle between prevoius vector (in case of first point A previous vector can be {0,-1,0} which is unit vector in -Y axis direction). And so on…
Points should be in XY plane.

Radovan

1 Like

@RadovanG, thanks. I’ll try to convert it to python.

_
c.