Topic: Cleanup self-intersecting 3d-plines

[Jeff_M]

I have a routine that creates a variable distance offset polyline from another polyline. It all works great except at angle points where, at an inside corner, the offset pline creates a loop and intersects itself. I've been beating my head against the wall on this for days now and have yet to come up with a viable solution to eliminate that loop. (There likely will be more than one of these.) Ideally, In the image below, the white crosses are my calculated offset points for the new (cyan) pline. Ideally, I'd like to not even show these temporary markers where they would be creating the loop portion of the pline.

The plines are generally not closed, and they can zig-zag creating many of these loops. I'm at a complete loss as to how to proceed.

Any suggestions would be helpful.

[kaefer]

Eliminating self-intersecting loops? That's really easy, aside from obvious definition issues of "intersection" (apparent, projected or real?): Just a) compute the self-intersection points and b) rebuild the polyline.

As for a), at least in 2d there's a viable way by converting the polyline into a CompositeCurve2d and using a  CurveCurveIntersector2d with both curves set to the former. Drawback is that I know of no way to generate meaningful curve parameters for GetSplitCurves, so we're stuck with the point, which isn't helpful to split a polyline into 3 parts.

Regarding rebuilding the polyline, using a directed acyclic graph may be algorithmic overkill, but would ensure an optimal solution.

[Bryco]

I think you need to make lines form the points with z=0, offset them, find the intersections then put the original z back in.

Code: [Select]

        [CommandMethod("rw")]
        public static void OffsetP3d()
        {
            Document doc = acadApp.DocumentManager.MdiActiveDocument;
            Editor ed = doc.Editor;
            Database db = doc.Database;
            PromptEntityOptions peo = new PromptEntityOptions("P");
            peo.SetRejectMessage("3dp");
            peo.AddAllowedClass(typeof(Polyline3d),true);
            PromptEntityResult per = ed.GetEntity(peo);
            if (per.Status != PromptStatus.OK) return;               

            using (Transaction tr = db.TransactionManager.StartTransaction())
            {
                Polyline3d p3d = tr.GetObject(per.ObjectId, OpenMode.ForWrite) as Polyline3d;
               
                BlockTableRecord space = tr.GetObject(db.CurrentSpaceId, OpenMode.ForWrite) as BlockTableRecord;
                Point3d p1=Point3d.Origin;
                int cnt = 0;

                Line l1=new Line();
                Line line=new Line();
                Polyline3d P3dOffset = new Polyline3d();
                space.AppendEntity(P3dOffset);
                tr.AddNewlyCreatedDBObject(P3dOffset, true);
                double z = 0;
                foreach(ObjectId id in p3d)
                {
                    PolylineVertex3d v = tr.GetObject(id, OpenMode.ForRead) as PolylineVertex3d;
                    if (cnt == 0)
                    {
                        p1 =new Point3d( v.Position.X,v.Position.Y,0);
                        cnt++;
                        continue;
                    }
                    Point3d p2 = new Point3d(v.Position.X, v.Position.Y, 0);
                    z = v.Position.Z;
                    Line l=new Line(p1,p2);
                    line = (Line)l.GetOffsetCurvesGivenPlaneNormal(Vector3d.ZAxis, 1)[0];
                    p1 = p2;
                    if (cnt == 1)
                    {
                        P3dOffset.AppendVertex(new PolylineVertex3d
                            (new Point3d(line.StartPoint.X, line.StartPoint.Y, z)));
                    }
                    if(cnt>1)
                    {
                        Point3dCollection pts=new Point3dCollection();
                        l1.IntersectWith(line, Intersect.ExtendBoth, pts, 1,0);
                        if (pts.Count==0) continue;
                        Point3d p = pts[0];
                        P3dOffset.AppendVertex(new PolylineVertex3d(new Point3d(p.X,p.Y,z)));
                    }
                    l1 = line;
                    cnt++;
                }
                P3dOffset.AppendVertex(new PolylineVertex3d(new Point3d(line.EndPoint.X, line.EndPoint.Y, z))); 
                tr.Commit();
            }

        }

[Jeff_M]

Thanks, Bryco, but this isn't quite what I'm looking for. It misses the intersection and removes vertices that should remain....and I'm not sure why you used the offset.

Kaefer, I'm not finding this to be "really easy" which is why I've asked the question.  I have no idea what a directed acyclic graph is, even after Googling it and reading up on them...waaaaay over my old brain's comprehension level.

As for creating the intersection points, I don't necessarily need these. What I DO need is for the vertices that create the loop be omitted from the pline, which could allow the 2 vertices on either end of the crossing segments to be connected

And perhaps I didn't explain what I'm doing very well. The purpose is to draw the "daylight" line created by sloping from a known pline (which may be 2d at elevation or 3d with differing elevations) to a ground surface. In order to  provide a more accurate daylight, I allow the user to add supplemental points at a user entered interval. So where the original pline might have 5 angle points over a 600' distance, the daylight pline would have ~67 vertices, all at different offsets from the original pline. The attached image should explain it better. The 2 magenta lines are where the modified pline could connect instead of the actual intersection. It's also possible that, due to the slopes used from the original pline, that where there are 2 (or more) angle points close to each other, these could end up creating just 1 angle point along the daylight (I hope that makes sense).

Thanks again for the suggestions.

[Jeff H]

Just wondering and I wish I was smart enough to understand kaefer's solution.

I thought I had code for this, and before i try this
but would this work?
 
First find all the intersection points then iterated the polyline adding each vertex until you reach one that the segment contains a intersecting point.
Then starting with the next segment loop until you reach the next segment that contains the same intersecting point then start adding the vertices again until you find one that contains a intersecting point......
 
Are do it all at once instead breaking it up

[kaefer]

First find all the intersection points then iterated the polyline adding each vertex until you reach one that the segment contains a intersecting point.
Then starting with the next segment loop until you reach the next segment that contains the same intersecting point then start adding the vertices again until you find one that contains a intersecting point......

Yes, exactly. There would be a couple of restrictions: Loops cannot be nested, and the first and possibly the last vertex of the polyline cannot constitute a loop. That wouldn't appear to be an issue according to the other Jeff's clarification.

Maybe I'm the one who's thinking too complicated here. After having decided to deploy a CurveCurveIntersector to find the intersections it was an afterthought to use the IsOn method, to match an intersection point with the corresponding Curve2d.

I've no idea how to reason about this in an imperative setting, special cases abound and mutables are all over the place. Please excuse the choice of language.

Handling the points...

Code: [Select]

// Head of list option
let tryHead = function h::_ -> Some h | _ -> None

// Finds self-intersections and returns new array with loops removed,
// or an empty array when no loops were found
let removePoints (curve3dPts : Point3d[]) =
    let curve2ds =
        curve3dPts           
        |> Seq.pairwise
        |> Seq.map (fun (sp, ep) ->
            new LineSegment2d(
                new Point2d(sp.X, sp.Y), new Point2d(ep.X, ep.Y) )
                    :> Curve2d )
        |> Array.ofSeq
    let cc = new CompositeCurve2d(curve2ds)
    use cci2d = new CurveCurveIntersector2d(cc, cc)
    if cci2d.NumberOfIntersectionPoints = 0 then Array.empty
    else
        let intSec2dPts =
            [ for i in 0 .. cci2d.NumberOfIntersectionPoints - 1 ->
                cci2d.GetIntersectionPoint i ]
        // Map each Curve2d to its corresponding 3d start point
        // and any indices of intersections local to the Curve2d
        let pts3dIds =
            (curve3dPts.[0..curve3dPts.Length - 2], curve2ds)
            ||> Array.map2 (fun pt c ->
                let ids =
                    intSec2dPts
                    |> List.mapi (fun i ipt ->
                        if c.IsOn ipt then Some i else None)
                    |> List.choose id
                pt, ids)

        // Start point of polyline is always present in result,
        // include other points depending on their position
        // relative to a specific intersection index
        let (fstPt, fstIds) = pts3dIds.[0]
        pts3dIds.[1 .. pts3dIds.Length - 1]
        |> Array.scan (fun (res, act) (pt, ids) ->
            match act with
            | None ->                                Some pt, tryHead ids
            | Some x when List.exists ((=) x) ids -> None, None
            | Some x ->                              None, Some x )
            (Some fstPt, tryHead fstIds)
        |> Array.choose fst

... and the polyline.

Code: [Select]

let removeLoops pl3dId =
    let doc = acadApp.DocumentManager.MdiActiveDocument
    let db = doc.Database
    let ed = doc.Editor
    use tr = db.TransactionManager.StartTransaction()
    let cs =
        tr.GetObject(db.CurrentSpaceId, OpenMode.ForWrite)
            :?> BlockTableRecord
    let pl3d =
        tr.GetObject(pl3dId, OpenMode.ForRead)
            :?> Polyline3d
    let pl3dPts =
        [| for v3dId in pl3d do
            let v = tr.GetObject(v3dId :?> ObjectId, OpenMode.ForRead) :?> PolylineVertex3d
            yield v.Position |]
    let arg3dPts =
        if pl3d.Closed then Array.append pl3dPts [| pl3dPts.[0] |] else pl3dPts
    match removePoints arg3dPts with
    | [||] -> ed.WriteMessage("\nPolyline not self-intersecting. ")
    | res3dPts ->
        let npl3d = new Polyline3d(Closed = pl3d.Closed)
        cs.AppendEntity npl3d |> ignore
        tr.AddNewlyCreatedDBObject(npl3d, true)
        let addVertex pt3d =
            let v3d = new PolylineVertex3d(pt3d)
            npl3d.AppendVertex v3d |> ignore
            tr.AddNewlyCreatedDBObject(v3d, true)
        for pos in res3dPts do
            addVertex pos
        // Add end point, if not implicit.
        if not pl3d.Closed then
            addVertex arg3dPts.[arg3dPts.Length - 1]
    tr.Commit()


[Jeff H]

Yes, exactly. There would be a couple of restrictions: Loops cannot be nested, and the first and possibly the last vertex of the polyline cannot constitute a loop. That wouldn't appear to be an issue according to the other Jeff's clarification.

I could not find the code I thought I had and started messing with it a little and hit a snag or two under certain conditions, but getting your conformation makes me feel better that I was not completely headed in the wrong direction.
 
When I start messing with again a little later and finish up I will post it, but I already know it is not very efficient(efficient enough for a test drawing with a handful of polylines with around 10 vertices).
I keep meaning to learn F# and cannot completely digest your code but able to pick up enough hints to see ways  simplify the approach I was taking.
 
Thanks kaefer and thanks Jeff for the fun question.

[Jeff_M]

Thanks kaefer. This will take some time for me to decipher, but I think I see where you went with it.

Just in case I don't get back here tomorrow, I will be gone for 4 days with no internet access. So it's not that I've forgotten about this or given up, I will be back (unless I crack it before I take off).

[Jeff H]

At first I was replacing all the vertices in loop with one at the intersection point and I see you just go to the next which makes it easier but how do handle this situation in the picture?
The numbers are of course the vertices.
 
Do you go from 0 to 5?
 
 
 
 
 

[kaefer]

Do you go from 0 to 5?

You've got me there. I would say I should.

The segment from Vertex 2 to Vertex 3 has two self-intersections:
0, 1, [ 1 ]
1, 2, [ ]
2, 3, [ 0; 1 ]
3, 4, [ ]
4, 5, [ 0 ]

More sensible code at the end of removePoints:

Code: [Select]

        pts3dIds.[1 .. pts3dIds.Length - 1]
        |> Array.scan (fun (res, act) (pt, ids) ->
            match act with
            | None ->                                Some pt, tryHead ids
            // Check when closing a loop for the start of another one
            | Some x when List.exists ((=) x) ids -> None, List.filter ((<>) x) ids |> tryHead
            | Some x ->                              None, Some x )
            (Some fstPt, tryHead fstIds)
        |> Array.choose fst

I find F# great for prototyping, but when too complex (or too functional) it can be hard porting it into C#.

[kaefer]

In the end it wasn't half as hard as expected, to do it in C#, with a little help of Linq. In the olden days the saying went: One can write BASIC in any language. With the advent of v3.0 it's astonishingly different.

Loop detection could still use a bit polishing:

Code: [Select]

        // Omit points between pairs of equal intersection indices
        static IEnumerable<Point3d> ScanForLoops(IEnumerable<Pt3dIds> pis, Pt3dIds pi0)
        {
            Func<int[],int?> tryFirst = ids =>
            {
                return ids.Length > 0 ? new Nullable<int>(ids[0]) : null;
            };
            int? x = tryFirst(pi0.ids);
            yield return pi0.pt;
            foreach(Pt3dIds pi in pis)
            {
                if(!x.HasValue)
                {
                    x = tryFirst(pi.ids);
                    yield return pi.pt;
                }
                else if(pi.ids.Contains(x.Value))
                {
                    x = tryFirst(pi.ids.Where(y => y != x.Value).ToArray());
                }
            }
        }

Complete code file attached.

[Jeff_M]

Thank you kaefer! I won't be able to look at this until sometime over the weekend, but I most definitely will!

Jeff

[Jeff H]

I never knew about the CurveCurveIntersector2d & CurveCurveIntersector3d class until I saw kaefer's solution and his solution looks like the way to go.
 
This one is stripped down, and uses IntersectWith  which kaefer avoids(look in docs at CurveCurveIntersector2d & CurveCurveIntersector3d),
and I am not smart enough to keep it all straight in my head so I had to create some ADT's -- PolyLineSegment & HighestIntersectingVertex

 
You pass in the polyline and one of the NonSelfIntersectingPolylineType values.

Code: [Select]

    public enum NonSelfIntersectingPolylineType
    {
        RemoveLoops = 0,
        ReplaceLoopsWithIntersectingPoints
    }

In the pic below:
If RemoveLoops is used the result is the cyan polylines
If ReplaceLoopsWithIntersectingPoints is used the result is the red polylines

Part of code and the rest is attached

Code: [Select]

        public Polyline CreateNonIntersectingPolyline(Polyline poly, NonSelfIntersectingPolylineType type)
        {
            createPolyLineSegments(poly);
            Polyline pl = new Polyline();           
            for (int i = 0; i < PolyLineSegments.Count; i  )
            {
                PolyLineSegment pls = PolyLineSegments[i];
                pl.AddVertexAt(pl.NumberOfVertices, new Point2d(pls.VertexPoint.X, pls.VertexPoint.Y), 0, 0, 0);
               
                if (pls.HasHigherIntersectingVertex())
                {
                   
                    while (pls.HasHigherIntersectingVertex())
                    {
                        if (type == NonSelfIntersectingPolylineType.ReplaceLoopsWithIntersectingPoints)
                        {
                            pl.AddVertexAt(pl.NumberOfVertices, new Point2d(pls.HighestIntersection.IntersectionPoint.X, pls.HighestIntersection.IntersectionPoint.Y), 0, 0, 0);
                        }
                        i = pls.HighestIntersection.VertexIndex;
                        pls = PolyLineSegments[i];
                    }                 
                }
            }
            return pl;
        }

 

[Jeff_M]

OK, finally got a chance to do some testing. Kaefer's c# code is working in most cases. The one case I've found it to fall down is when the first/last segments cross. This can occur if the user is modeling a pond and daylighting to the inside of a closed polygon. If I open that polygon and keep the start/end points separate by a bit to eliminate that crossing then the other corners are trimmed correctly. But keeping the closed polygon results in a new 3dpolyline having just the start & end vertices.

I will do my best to learn from what you've provided (I'm seeing things in your code I've never seen/used before) and try to fix this.

Jeff, off to test your code now. This looks similar to what I have been attempting. Will be back with results....

Again, I can't thank you both enough for helping me solve this riddle!

[Jeff_M]

OK, testing Jeff's and I'm seeing the same problem when used on an offset of a closed polygon. It also only works with Polylines and not the Polyline3d, but this is pretty easy to workaround.

I will be messing around with these tonight to try to solve this last little problem.

Here's a picture of the results I get with offset of a close polygon. Note the white X's that should be connected are not, but the 'overlaps' at the start/end are included when they shouldn't be.