3

I want a function that can fill small gaps between parts of a multipart polyline.

The polylines represent events on roads and are linear, not branched. The parts may be separated by small or large gaps. I want to connect any parts that are separated by a small gap less than some tolerance, while leaving the large gaps. No part is shorter than the gap tolerance.

The parts are (probably) well ordered and directed (meaning that geographically adjacent parts are adjacent in the IGeometryCollection, and parts all point in the same direction), but bonus points for an implementation that can handle misordered or flipped parts.

I'm working in ArcObjects, C#, .NET 4.0, ArcGIS 10.1SP1.

Update2: Some example test data that, when filled with a gap tolerance of 0.25, should yield 3 parts and length 7.0910109853324:

//first two parts have no gap (thus they could really be a single part, 
//but this might happen. They'll get merged in the output)
gcoll.AddGeometry(MakePath(MakePoint(0.0, 0.0), MakePoint(1.0, 0.0)));
gcoll.AddGeometry(MakePath(MakePoint(1.0, 0.0), MakePoint(2.0, 0.0)));
//third part has small gap and is flipped. 
//it is also set at an angle to check that we don't remove a vertex.
gcoll.AddGeometry(MakePath(MakePoint(3.0, 0.0), MakePoint(2.1, 0.1)));
//gcoll.AddGeometry(MakePath(MakePoint(2.1, 0.1), MakePoint(3.0, 0.0)));
//fourth part has large gap
gcoll.AddGeometry(MakePath(MakePoint(5.0, 0.0), MakePoint(6.0, 0.0)));
//fifth part has small gap
//it is also set at an angle to check that we don't remove a vertex.
gcoll.AddGeometry(MakePath(MakePoint(6.1, 0.1), MakePoint(8.0, 0.0)));
//sixth part has large gap
gcoll.AddGeometry(MakePath(MakePoint(9.0, 0.0), MakePoint(10.0, 0.0)));

2 Answers 2

2

This is an update to my old answer. The following tests ok for polylines that have gaps between their paths and are not oriented. Updated again to keep gaps that are larger than tolerance, and sort the resulting multipart polylines.

The key concept is to create a Directed Graph (Digraph) of the polyline. A hash function is used to assign node IDs points such that points near each other (based on tolerance) get the same ID. Once a Digraph is created, topological rules can be checked. This method employs two rules:

  • exactly at least two nodes have valence (or degree) of 1
  • valence for all remaining nodes is 2.

Imposing these rules make it easy to topologically sort the graph, and then chain together the coordinates for the polylines.

An edge is a connection between two nodes.

public class Edge
{
    public Edge(IPath path, int i, double tol)
    {
        this.ID = i;
        this.FnodeID = HashPoint(path.FromPoint,tol);
        this.TnodeID = HashPoint(path.ToPoint,tol);
    }

    public string FnodeID { get; set; }
    public string TnodeID { get; set; }
    public int ID { get; set; }
    public static string HashPoint(IPoint p, double tol)
    {
        int x = (int)Math.Round(p.X / tol);
        int y = (int)Math.Round(p.Y / tol);
        return string.Format("{0}_{1}", x, y);
    }
}

A DiGraph is a collection of Edges

public class DiGraph
{
    public static void TestFillGaps()
    {
        double tol = 0.1;
        IGeometryCollection gcoll = new PolylineClass();
        gcoll.AddGeometry(MakePath(MakePoint(1.0, 1.0), MakePoint(2.99, 2.99)));
        gcoll.AddGeometry(MakePath(MakePoint(3.99, 3.99), MakePoint(3.0, 3.0)));
        gcoll.AddGeometry(MakePath(MakePoint(0.0, 0.0), MakePoint(0.99, 0.99)));
        gcoll.AddGeometry(MakePath(MakePoint(10.0, 10.0), MakePoint(10.99, 10.99)));

        gcoll.AddGeometry(MakePath(MakePoint(21.0, 21.0), MakePoint(22.0, 22.0)));
        gcoll.AddGeometry(MakePath(MakePoint(22.0, 22.0), MakePoint(23.0, 23.0)));
        gcoll.AddGeometry(MakePath(MakePoint(24.0, 24.0), MakePoint(23.0, 23.0)));

        gcoll.AddGeometry(MakePath(MakePoint(11.0, 11.0), MakePoint(12.0,12.0)));
        gcoll.AddGeometry(MakePath(MakePoint(12.0, 12.0), MakePoint(13.0, 13.0)));
        gcoll.AddGeometry(MakePath(MakePoint(14.0, 14.0), MakePoint(13.0, 13.0)));


        Debug.Print("points before cleaning:");
        ListPoints((IPolyline)gcoll);
        Debug.Print("paths before cleaning:");
        ListPaths((IPolyline)gcoll);
        DiGraph d = new DiGraph((IPolyline)gcoll, tol);
        var pline = d.FillGaps();
        Debug.Print("paths after cleaning:");
        ListPaths(pline);
        Debug.Print("points after cleaning:");
        ListPoints(pline);
    }

    private IGeometryCollection _g;
    private double _tol;
    private Dictionary<string, List<Tuple<Edge, bool>>> _dict = new Dictionary<string,List<Tuple<Edge,bool>>>(StringComparer.CurrentCultureIgnoreCase);
    public DiGraph(IPolyline p, double tol)
    {            
        _g = p as IGeometryCollection;
        _tol = Math.Abs(tol);
        for (int i = 0; i < _g.GeometryCount; i++)
        {
            var path = _g.get_Geometry(i) as IPath;
            var e = new Edge(path, i, tol);
            if (!_dict.ContainsKey(e.FnodeID))
                _dict.Add(e.FnodeID, new List<Tuple<Edge, bool>>());
            if (!_dict.ContainsKey(e.TnodeID))
                _dict.Add(e.TnodeID, new List<Tuple<Edge, bool>>());
            // bool indicates orientation of the edge
            // (true when fnode of edge is the key)
            _dict[e.FnodeID].Add(new Tuple<Edge, bool>(e, true));
            _dict[e.TnodeID].Add(new Tuple<Edge, bool>(e, false));
        }
    }

    public IPolyline FillGaps()
    {
        var outPolyline = new PolylineClass() as IPolyline;
        outPolyline.SpatialReference = ((IPolyline)_g).SpatialReference;

        List<List<string>> nodeLists = SortNodes(_dict);
        foreach (List<string> nodeList in nodeLists)
        {
            IPath path = BuildPath(nodeList, _dict, _g);
            path = Clean(path,_tol);
            ((IGeometryCollection)outPolyline).AddGeometry(path);
        }

        int neighbor1;
        int neighbor2;
        bool isFrompoint1;
        bool isFrompoint2;
        FarthestNeighbors(outPolyline,out neighbor1,out isFrompoint1, out neighbor2, out isFrompoint2);
        IPath p1 = ((IGeometryCollection)outPolyline).get_Geometry(neighbor1) as IPath;
        IPath p2 = ((IGeometryCollection)outPolyline).get_Geometry(neighbor2) as IPath;

        IPoint startPoint = isFrompoint1 ? p1.FromPoint : p1.ToPoint;
        IPoint endPoint = isFrompoint2 ? p2.FromPoint : p2.ToPoint;

        Debug.Print("start point: {0} {1} {2}",Edge.HashPoint(startPoint,_tol),neighbor1,isFrompoint1);
        Debug.Print("end point: {0} {1} {2}",Edge.HashPoint(endPoint, _tol),neighbor2,isFrompoint2);

        // flip the paths and chain them in order
        outPolyline = Orient(outPolyline, startPoint);
        return outPolyline;
    }

    public static IPolyline Orient(IPolyline p, IPoint startPoint)
    {
        // orient and order the paths in the polyline such that
        // the order is based on proximity and they are flipped so
        // the To point of a path is close to the From point of the
        // following path.
        List<int> visited = new List<int>();
        var inGc = p as IGeometryCollection;
        IGeometryCollection outGc = new PolylineClass();
        ((IGeometry)outGc).SpatialReference = p.SpatialReference;

        IPoint pnt = startPoint;
        while (visited.Count < inGc.GeometryCount)
        {
            bool isFrom;
            int idx = FindPath(p, pnt, visited, out isFrom);
            IPath outPath = inGc.get_Geometry(idx) as IPath;
            outPath = ((IClone)outPath).Clone() as IPath;
            if (!isFrom)
            {
                idx = -idx;
                outPath.ReverseOrientation();
            }
            visited.Add(idx);
            outGc.AddGeometry(outPath);
        }
        return (IPolyline)outGc;
    }

    private static int FindPath(IPolyline p, IPoint pnt, List<int> visited, out bool isFrom)
    {
        // return the index of the path in p that is closest to pnt
        // and is not in the visite list.
        double dist = double.MaxValue;
        var proxOp = pnt as IProximityOperator;
        int iclosest = -1;
        isFrom = true;
        var g = p as IGeometryCollection;
        for (int i = 0; i < g.GeometryCount; i++)
        {
            if (visited.Contains(i) || visited.Contains(-i))
                continue;

            var path = g.get_Geometry(i) as IPath;
            double fDist = proxOp.ReturnDistance(path.FromPoint);
            if (fDist < dist)
            {
                dist = fDist;
                iclosest = i;
                isFrom = true;
            }
            double tDist = proxOp.ReturnDistance(path.ToPoint);
            if (tDist < dist)
            {
                dist = tDist;
                iclosest = i;
                isFrom = false;
            }
        }
        return iclosest;
    }

    public static void FarthestNeighbors(IPolyline p, out int n1, out bool isFrom1, out int n2, out bool isFrom2)
    {
        // find the index of the two paths in p such that two of the paths endpoints
        // are farther from each other than any other pair of endpoints.

        // key is idx of path in polyline
        // Tuple such that: 
        //     Item2 - farthest distance from the endpoint
        //     Item3 - is true when Item2 is an From point, false otherwise.
        IGeometryCollection g = p as IGeometryCollection;
        Dictionary<int, Tuple<int, double, bool>> distDict = new Dictionary<int, Tuple<int, double, bool>>();
        for (int i = 0; i < g.GeometryCount; i++)
        {
            var path = g.get_Geometry(i) as IPath;
            var t1 = GetFarthest(path.FromPoint, g, i);
            var t2 = GetFarthest(path.ToPoint, g, i);
            if (t1.Item2 > t2.Item2)
                distDict.Add(i, t1);
            else
                distDict.Add(i, t2);
        }

        int imax = -1;            
        double distMax = double.MinValue;
        foreach (int i in distDict.Keys)
        {
            if (distDict[i].Item2 > distMax)
            {
                imax = i;
                distMax = distDict[i].Item2;
            }
        }
        n1 = imax;
        n2 = distDict[imax].Item1;
        isFrom1 = distDict[imax].Item3;
        isFrom2 = distDict[n2].Item3;
    }

    private static Tuple<int,double,bool> GetFarthest(IPoint p,IGeometryCollection g, int iskip)
    {
        bool isFrompoint = false;
        double farthestDist = double.MinValue;
        int ifarthest = -1;
        for (int i = 0; i < g.GeometryCount; i++)
        {
            if (i == iskip)
                continue;
            IPath path = g.get_Geometry(i) as IPath;
            double dist = ((IProximityOperator)p).ReturnDistance(path.FromPoint);
            if (dist > farthestDist)
            {
                farthestDist = dist;
                isFrompoint = true;
                ifarthest = i;
            }
            dist = ((IProximityOperator)p).ReturnDistance(path.ToPoint);
            if (dist > farthestDist)
            {
                farthestDist = dist;
                isFrompoint = false;
                ifarthest = i;
            }
        }
        return new Tuple<int, double, bool>(ifarthest, farthestDist, isFrompoint);
    }

    public static IPath Clean(IPath p, double tol)
    {
        var inpc = p as IPointCollection;
        var outpc = new PathClass() as IPointCollection;
        ((IGeometry)outpc).SpatialReference = p.SpatialReference;
        IPoint p1 = inpc.get_Point(0);
        for (int i = 1; i < inpc.PointCount; i++)
        {
            IPoint p2 = inpc.get_Point(i);
            if (((IProximityOperator)p2).ReturnDistance(p1) > tol)
            {
                outpc.AddPoint(p2);
                p1 = p2;
            }
        }
        if (outpc.PointCount < 2)
            return null; // throw an exception?
        return (IPath)outpc;
    }

    public static List<List<string>> SortNodes(Dictionary<string, List<Tuple<Edge, bool>>> nodeDict)
    {
        List<string> dangles = new List<string>();
        foreach (string id in nodeDict.Keys)
        {
            if (nodeDict[id].Count == 1)
                dangles.Add(id);
            else if (nodeDict[id].Count > 2)
                throw new Exception("internal intersection");
        }
        if (dangles.Count < 2)
            throw new Exception("topological violation");

        // a list of paths (a path is a list of nodes)
        List<List<string>> outList = new List<List<string>>();
        var visited = new List<string>();
        foreach (string id in dangles)
        {
            if (!visited.Contains(id))
                outList.Add(SortNodes(nodeDict,id,visited));
        }
        return outList;
    }

    public static List<string> SortNodes(Dictionary<string, List<Tuple<Edge, bool>>> nodeDict, string id,List<string> visited)
    {
        // walk along the path
        List<string> outList = new List<string>();
        while (!string.IsNullOrEmpty(id))
        {
            outList.Add(id);
            visited.Add(id);
            Tuple<Edge, bool> nextT = null;
            string nextID = null;
            foreach (Tuple<Edge, bool> t in nodeDict[id])
            {
                if (t.Item2)
                    nextID = t.Item1.TnodeID;
                else
                    nextID = t.Item1.FnodeID;

                if (!visited.Contains(nextID))
                {
                    nextT = t;
                    break;
                }
            }
            if (nextT != null)
                id = nextID;
            else
                id = null;
        }
        return outList;
    }

    private static IPath BuildPath(List<string> nodeIDs, Dictionary<string, List<Tuple<Edge, bool>>> nodeDict, IGeometryCollection g)
    {
        var outPath = new PathClass() as IPointCollection;
        for (int i = 1; i < nodeIDs.Count; i++)
        {
            string fnodeID = nodeIDs[i - 1];
            string tnodeID = nodeIDs[i];
            bool reversed;
            var edge = FindEdge(fnodeID, tnodeID,nodeDict.Values.ToList(), out reversed);
            var path = g.get_Geometry(edge.ID) as IPath;
            path = ((IClone)path).Clone() as IPath;
            if (reversed)
                path.ReverseOrientation();
            outPath.AddPointCollection((IPointCollection)path);
        }
        return (IPath)outPath;
    }

    private static Edge FindEdge(string fnodeID, string tnodeID, List<List<Tuple<Edge, bool>>> edgesll, out bool reversed)
    {
        reversed = false;
        foreach (var list in edgesll)
        {
            foreach (var t in list)
            {
                if (t.Item1.FnodeID == fnodeID && t.Item1.TnodeID == tnodeID)
                    return t.Item1;
                else if (t.Item1.TnodeID == fnodeID && t.Item1.FnodeID == tnodeID)
                {
                    reversed = true;
                    return t.Item1;
                }
            }
        }
        // throw an exception?
        return null;
    }


    private static void ListPaths(IPolyline pLine)
    {
        var g = pLine as IGeometryCollection;
        for (int i = 0; i < g.GeometryCount; i++)
        {
            IPath p = g.get_Geometry(i) as IPath;
            Debug.Print("{0} : {1},{2} - {3},{4}", i, p.FromPoint.X, p.FromPoint.Y, p.ToPoint.X, p.ToPoint.Y);
        }
    }
    private static void ListPoints(IPolyline p)
    {
        var pc = p as IPointCollection;
        for (int i = 0; i < pc.PointCount; i++)
            Debug.Print("{0}: {1}, {2}", i, pc.get_Point(i).X, pc.get_Point(i).Y);
    }

    public static IPath MakePath(IPoint p1, IPoint p2)
    {
        IPointCollection pc = new PathClass();
        pc.AddPoint(p1);
        // include a midpoint
        pc.AddPoint(MakePoint((p2.X + p1.X) / 2.0, (p2.Y + p1.Y) / 2.0));
        pc.AddPoint(p2);
        ((IGeometry)pc).SpatialReference = p1.SpatialReference;
        return (IPath)pc;
    }

    public static IPoint MakePoint(double x, double y)
    {
        IPoint p = new PointClass();
        p.PutCoords(x, y);
        return p;
    }

}

The test uses paths that have 3 points each (midpoints added). Here are the results:

points before cleaning:
0: 1, 1
1: 1.995, 1.995
2: 2.99, 2.99
3: 3.99, 3.99
4: 3.495, 3.495
5: 3, 3
6: 0, 0
7: 0.495, 0.495
8: 0.99, 0.99
9: 10, 10
10: 10.495, 10.495
11: 10.99, 10.99
12: 21, 21
13: 21.5, 21.5
14: 22, 22
15: 22, 22
16: 22.5, 22.5
17: 23, 23
18: 24, 24
19: 23.5, 23.5
20: 23, 23
21: 11, 11
22: 11.5, 11.5
23: 12, 12
24: 12, 12
25: 12.5, 12.5
26: 13, 13
27: 14, 14
28: 13.5, 13.5
29: 13, 13
paths before cleaning:
0 : 1,1 - 2.99,2.99
1 : 3.99,3.99 - 3,3
2 : 0,0 - 0.99,0.99
3 : 10,10 - 10.99,10.99
4 : 21,21 - 22,22
5 : 22,22 - 23,23
6 : 24,24 - 23,23
7 : 11,11 - 12,12
8 : 12,12 - 13,13
9 : 14,14 - 13,13
start point: 0_0 0 False
end point: 240_240 2 False
paths after cleaning:
0 : 0,0 - 3.495,3.495
1 : 10.495,10.495 - 14,14
2 : 21.5,21.5 - 24,24
points after cleaning:
0: 0, 0
1: 0.495, 0.495
2: 1, 1
3: 1.995, 1.995
4: 3, 3
5: 3.495, 3.495
6: 10.495, 10.495
7: 10.99, 10.99
8: 11.5, 11.5
9: 12, 12
10: 12.5, 12.5
11: 13, 13
12: 13.5, 13.5
13: 14, 14
14: 21.5, 21.5
15: 22, 22
16: 22.5, 22.5
17: 23, 23
18: 23.5, 23.5
19: 24, 24
4
  • Pretty sexy K! However, if there are any gaps larger than the tolerance (which shouldn't be filled), this gives a 'not single stranded' exception. For example, add this part to the test polyline: gcoll.AddGeometry(MakePath(MakePoint(4.5, 4.5), MakePoint(6.0, 6.0)));
    – MC5
    Commented Oct 23, 2014 at 19:33
  • see update to preserve gaps larger than tolerance. Commented Oct 23, 2014 at 21:54
  • close! if a part has only a single segment, Clean() returns null. For example, comment out your 6th and 7th AddGeometry() statements in your test method. I hate to even bring this up because you have gone way above/beyond on this one!
    – MC5
    Commented Oct 23, 2014 at 23:31
  • i think that was an easy fix, replacing 'return null; // throw an exception?' with 'return p'. However, I'm not getting the right polyline length with my test data (added to question), probably something to do with forming the polyline out of points instead of segments. I'll look deeper tomorrow.
    – MC5
    Commented Oct 23, 2014 at 23:52
1

KK was quick on the trigger while I was in progress with similar.

KK's answer assumes each part is linear (consisting only of a start and end vertex). This answer is similar but does not make that assumption, but otherwise has the same limitation that it assumes that there are no misordered or flipped parts.

    /// <summary>
    /// Test PolylineGapFiller class on a polyline
    /// with ordered parts with a mix of small and large gaps.
    /// Input polyline have total length of 6.7 and 6 parts, with 4 gaps - two small, two large.
    /// Output polyline should have total length of 7.0 and 3 parts with 2 large gaps.
    /// </summary>
    [TestMethod]
    public void TestPolylineGapFiller()
    {
        IGeometryCollection inp = new PolylineClass();
        (inp as IPolyline).SpatialReference = GetSpatialReference(3071);

        //first two parts have no gap (thus they could really be a single part,
        //but this might happen. They'll get merged in the output)
        inp.AddGeometry(MakePath(MakePoint(0.0, 0.0), MakePoint(1.0, 0.0)));
        inp.AddGeometry(MakePath(MakePoint(1.0, 0.0), MakePoint(2.0, 0.0)));
        //third part has small gap
        inp.AddGeometry(MakePath(MakePoint(2.1, 0.0), MakePoint(3.0, 0.0)));
        //fourth part has large gap
        inp.AddGeometry(MakePath(MakePoint(5.0, 0.0), MakePoint(6.0, 0.0)));
        //fifth part has small gap
        inp.AddGeometry(MakePath(MakePoint(6.2, 0.0), MakePoint(8.0, 0.0)));
        //sith part has large gap
        inp.AddGeometry(MakePath(MakePoint(9.0, 0.0), MakePoint(10.0, 0.0)));
        Debug.Print("Input geometry has {0} parts and length {1}.", inp.GeometryCount, (inp as IPolyline).Length);

        PolylineGapFiller filler = new PolylineGapFiller();
        IPolyline result = filler.FillGapsOrdered((IPolyline)inp, 0.25);
        int partcount = (result as IGeometryCollection).GeometryCount;
        double length = result.Length;
        Debug.Print("Output geometry has {0} parts and length {1}.", partcount, length);
        Assert.AreEqual(7.0, length);
        Assert.AreEqual(3, partcount);
    }

    public class PolylineGapFiller
    {
    /// <summary>
    /// Fills small gaps in a polyline.
    /// This implementation assumes the parts are in order and all pointed in the same direction.
    /// </summary>
    /// <param name="p">the polyline on which to fill gaps. Make sure it has an SR</param>
    /// <param name="tol">the gap tolerance. Gaps smaller than this distance will be filled.</param>
    /// <returns>A new polyline with small gaps filled</returns>
    public IPolyline FillGapsOrdered(IPolyline inpolyline, double tol)
    {
        IGeometryCollection ingc = inpolyline as IGeometryCollection;
        if (ingc.GeometryCount < 2)
            throw new ArgumentException("You can only fill gaps on a polyline with >=2 parts.");

        IPolyline result = new PolylineClass();
        result.SpatialReference = inpolyline.SpatialReference;
        IGeometryCollection outgc = result as IGeometryCollection;
        //add the first part to the output polyline
        outgc.AddGeometry((IPath)((IClone)(ingc.get_Geometry(0) as IPath)).Clone());
        for (int i = 1; i < ingc.GeometryCount; i++)
        {
            IPath p1 = ingc.get_Geometry(i - 1) as IPath; //the prior part
            IPath p2 = ingc.get_Geometry(i) as IPath; //the current part
            double gap = ((IProximityOperator)p1.ToPoint).ReturnDistance(p2.FromPoint);
            if (gap > 0.0 && gap < tol)
            {
                //fill the gap
                outgc.AddGeometry(MakePath(p1.ToPoint, p2.FromPoint));
            }
            //add the current part to the output polyline
            outgc.AddGeometry((IPath)((IClone)(p2)).Clone());
        }
        outgc.GeometriesChanged();
        (result as ITopologicalOperator).Simplify();
        return result;
    }

    protected IPath MakePath(IPoint p1, IPoint p2)
    {
        IPointCollection pc = new PathClass();
        pc.AddPoint(p1);
        pc.AddPoint(p2);
        ((IGeometry)pc).SpatialReference = p1.SpatialReference;
        return (IPath)pc;
    }
}

Output:

Input geometry has 6 parts and length 6.7.
Output geometry has 3 parts and length 7.

Interestingly, if I hadn't set an SR on the input polyline, the ITopologicalOperator.Simplify() on the filled polyline causes vertices to get moved and the result's length gets wonky. I guess the lesson is to don't do that without an SR ;-)

Wonky Output:

Input geometry has 6 parts and length 6.7.
Output geometry has 3 parts and length 7.00000190734863.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.