# FindCurvePath for lines (rather than points)

I have need of a function to find a good ordering for a series of lines, as FindCurvePath does for points.

Sample data:

dat = {{{-2.83718,1.},{-2.83708,0.999885}},{{-2.837,0.999875},
{-2.83708,0.999885}},{{-2.83649,0.999763},{-2.83673,0.999716}},
{{-2.83673,0.999716},{-2.837,0.999875}},{{-2.83747,0.999718},{-2.83718,1.}},
{{-2.83699,0.999714},{-2.83697,0.999704}},{{-2.83696,0.999716},
{-2.8368,0.999686}},{{-2.83696,0.999716},{-2.83697,0.999704}},
{{-2.83678,0.999668},{-2.8368,0.999686}},{{-2.83702,0.999653},
{-2.83699,0.999714}},{{-2.83675,0.999644},{-2.83678,0.999668}},
{{-2.83647,0.999632},{-2.83649,0.999763}},{{-2.83647,0.999632},
{-2.8365,0.999633}},{{-2.8367,0.999603},{-2.83675,0.999644}},
{{-2.8365,0.999633},{-2.83654,0.999584}},{{-2.83666,0.99957},
{-2.8367,0.999603}},{{-2.83728,0.999697},{-2.83712,0.999592}},
{{-2.83664,0.999553},{-2.83666,0.99957}},{{-2.83654,0.999584},
{-2.83655,0.999551}},{{-2.83655,0.999551},{-2.83656,0.999549}},
{{-2.83712,0.999592},{-2.83702,0.999653}},{{-2.83656,0.999549},
{-2.83664,0.999553}}};


These lines form a single line:

ListLinePlot[dat, Frame -> True]


But they are out of order and their directions are mixed:

ListLinePlot[Join @@ dat, Frame -> True]
Graphics[Arrow @ dat, Frame -> True]


So I need not only to order the lines but to reverse some of them as well.

I also need to allow for gaps between lines. End points will not always be as close as in this example. A solution should work also on:

dat2 = dat ~Delete~ {{2}, {8}, {9}, {13}};

ListLinePlot[dat2, Frame -> True]


Additionally in practice my constituent lines are more than two points long but the end points should be sufficient for a solution. However I either need ordering and direction data that I can apply to the full lines or an algorithm that works on compound lines, not just line segments.

# Benchmarking

With multiple methods posted it is time to being benchmarking. The three methods are not entirely equivalent but I am making an effort to compare them as fairly as I can.

• Feyre's code does not return an explicit order but instead modified data
• my code is dependent on the specification of a suitable search radius
• Simon's FindShortestTour does not return an order starting with one of the ends
• I had to make a change and an addition to Simon's code to get consistent results
• I do not include application of the ordering produced by segOrder1 and segOrder2 in the benchmark but I found the overhead for that operation negligible

For random data I am using:

rdat[n_] :=
RandomSample /@ Partition[RandomReal[1, {n, 2}], 2, 1] // RandomSample


The functions as I am benchmarking them:

segOrder1[segs_, rad_: 0.0001] := (
Flatten[segs, 1]
// Nearest[# -> Automatic, #, {2, rad}] &
// Cases[{_, _}]
// Join[#, Partition[Range[2 Length@segs], 2]] &
// Graph
// FindPath[#, ## & @@ GraphPeriphery[#]] &
// First
)

segOrder2[segs_] :=
Module[{d = Flatten[segs, 1], dist},
dist[a_?OddQ, b_] /; (b == a + 1) := 0;
dist[a_, b_] := 1 + EuclideanDistance[d[[a]], d[[b]]];

FindShortestTour[Range @ Length @ d, DistanceFunction -> dist][[2]]
// If[#[[2]] === 2, Most, Rest][#] &
]

segReorder[dat_] :=
Module[{newdat, z, k, temp, it},
newdat = {dat[[1]]};
z = 1; k = 1;
While[k < Length@dat,
temp = Select[dat, FreeQ[Join[Reverse /@ newdat, newdat], #] &];
it = Table[
RegionDistance[Line@newdat[[k]], temp[[i, j]]], {i, Length[dat] - k}, {j, 2}];
z = Position[it, Min@it][[1, 1]];
If[it[[z, 1]] > it[[z, 2]], AppendTo[newdat, Reverse@temp[[z]]],
AppendTo[newdat, temp[[z]]]]; k++;];
newdat
]


Confirmation that they are working on this data:

SeedRandom[1]
dat = rdat[20];
o1 = segOrder1[dat];
o2 = segOrder2[dat];
newdat = segReorder[dat];

Partition[Flatten[dat, 1][[#]], 2] & /@ {o1, o2};
Append[%, newdat];
Graphics[Arrow@#, ImageSize -> 200] & /@ % // Row


Benchmark Plot

Needs["GeneralUtilities"]

BenchmarkPlot[{segOrder1, segOrder2, segReorder}, rdat, 5]


• You can use a loop with minimizing RegionDistance but my attempts got rather messy. Easy if you know begin and end lines though. Jan 25, 2017 at 19:49
• @Feyre If you feel like posting something don't worry about it being messy. I am curious to see multiple approaches. Jan 25, 2017 at 19:53
• Duplicate here maybe.
– yode
Feb 16, 2017 at 10:39
• @yode That would appear to be a more advanced problem than mine. I only need a single line start to finish, with no revisited points. I find that problem very interesting as well however and you just got my vote on it. I am sorry to see that there has not been more interest generally, but perhaps it was just posted at a poor time. Feb 16, 2017 at 11:38
• Related: (111460), (118132), and to a lesser degree (102618), (222252) Feb 16, 2017 at 11:39

Using FindShortestTour with a custom distance function:

d = Flatten[dat, 1];

dist[a_?OddQ, b_] /; (b == a + 1) := 0.0001 EuclideanDistance[d[[a]], d[[b]]]

dist[a_, b_] := EuclideanDistance[d[[a]], d[[b]]]

o = Most@FindShortestTour[Range[Length@d], DistanceFunction -> dist][[2]]
(* {1, 2, 4, 3, 8, 7, 6, 5, 24, 23, 25, 26, 29, 30, 37, 38, 39, \
40, 43, 44, 35, 36, 31, 32, 27, 28, 21, 22, 17, 18, 14, 13, 15, 16, \
12, 11, 20, 19, 42, 41, 34, 33, 9, 10} *)

Graphics[Arrow /@ Partition[d[[o]], 2]]


Update

A revised version which addresses Mr.Wizard's observations. Performance is still poor though.

segOrder2[segs_] :=
Module[{d = Flatten[segs, 1], dist, o},
dist[a_?OddQ, b_] /; (b == a + 1) := 0;
dist[a_, b_] := 1 + EuclideanDistance[d[[a]], d[[b]]];
o = FindShortestTour[Range[Length@d], DistanceFunction -> dist][[2]] //
If[#[[2]] === 2, Rest, Most][#] &;
RotateLeft[o, 2 Ordering[dist @@@ Partition[o, 2], -1] - 1]]

• Great! I had a feeling there was a more direct approach but I could not see my way to it. I knew I needed a way to "contract" the distance between points joined by existing lines but I struggled to do it. Treating odd indices seems obvious now but it was not at the time, like so many great ideas. Jan 26, 2017 at 5:23
• I'm probably having another moment of obtuseness but why can one not use dist[a_?OddQ, b_] /; (b == a + 1) := 0? Jan 26, 2017 at 6:16
• There is one issue with this as written: the order always starts with {1} whereas in practice it should start with one of the ends, i.e. {9} or {33}. Jan 26, 2017 at 7:29
• To handle adjacent line segments that start and end on the same point I needed to add a constant to the distance, e.g. dist[a_, b_] := 1 + EuclideanDistance[d[[a]], d[[b]]] Jan 26, 2017 at 7:54
• @Mr.Wizard, you can use zero distance. An earlier attempt did not work with zero as it allowed lines to be traversed multiple times in both directions. I hadn't realised that my workaround was no longer necessary. The start point is a problem. FindShortestTour identifies a complete circuit so it could perhaps be post-processed to split at the largest gap. Jan 26, 2017 at 19:46

This approach generates the data into newdat.

newdat = {dat[[1]]};
z = 1; k = 1;
While[k < Length@dat,
temp = Select[dat, FreeQ[Join[Reverse /@ newdat, newdat], #] &];
it = Table[
RegionDistance[Line@newdat[[k]], temp[[i, j]]], {i,
Length[dat] - k}, {j, 2}];
z = Position[it, Min@it][[1, 1]];
If[it[[z, 1]] > it[[z, 2]], AppendTo[newdat, Reverse@temp[[z]]],
AppendTo[newdat, temp[[z]]]]; k++;]


And the results:

ListLinePlot[Join @@ newdat, Frame -> True]
Graphics[Arrow@newdat, Frame -> True]


For the reduced data one arrow stays reversed.

With the missing piece from How do I "read out" the vertex names on this graph? I can self-answer using Nearest and Graph. Please don't let this post discourage answering as I am eager to see other approaches.

Now as a function with at least a little reusability. The second parameter is the search radius.

segOrder[segs_, rad_: 0.0001] := (
Flatten[segs, 1]
// Nearest[# -> Automatic, #, {2, rad}] &
// Cases[{_, _}]
// Join[#, Partition[Range[2 Length@segs], 2]] &
// Graph
// FindPath[#, ## & @@ GraphPeriphery[#]] &
// First
)

ListLinePlot[Part[Join @@ dat, segOrder[dat]], Frame -> True]


It works on the set with gaps given a sufficient radius:

ListLinePlot[Part[Join @@ dat2, segOrder[dat2, 0.0001]], Frame -> True]


### Extension

Here is my application of this ordering to the sorting (and joining) of longer lines.

lineSort[lines_, r_: 0.0001] :=
lines[[All, {1, -1}]] ~segOrder~ r ~Partition~ 2 //
Cases[ {a_, b_} :> lines[[⌈a/2⌉, b - a ;; a - b ;; b - a]] ]


Now I can do things like this:

geo = Import["http://www.rr4w.com/kml/9.kml"];

Cases[geo, Line[x_] :> x, {-4}] // lineSort // Catenate;

Graphics[{
Thickness[1/150],
Line[%, VertexColors -> Array[ColorData["Rainbow"], Length@%, {0, 1}]]
}]


• How about a larger test set to enable benchmarking? I'm very interested to learn about the speed of different approaches (very relevant for CNC path optimization...) Jan 25, 2017 at 20:29
• @YvesKlett Do you have a set to propose? If uniformly distributed pseudorandom data useful? Jan 25, 2017 at 20:38
• ... will take a look if I find a wicked set tomorrow (sorry, very busy at the office recently) . However, pseudorandom should work as well. Jan 25, 2017 at 20:41
• @YvesKlett A basic benchmark has been added to the Question. I look forward to a real world set of a larger size than my original sample. Jan 26, 2017 at 9:19

For the sake of general,such as the line include not only $$2$$ points,the every line maybe is a mess order(include first one),I make some adjust like following.Maybe too long,but powerful.

### If you have a gap data like dat2

We can convert it into non-gap data first by ConnectComponentPoints

ConnectComponentPoints[p_] :=
Module[{f, var1, var2, nearePoint, graph}, f = Nearest /@ Most[p];
var2 = Drop[p, #] & /@ Range[Length[p] - 1];
var1 = MapThread[Catenate /@ # /@ #2 &, {f, var2}];
nearePoint =
Catenate[
Map[First[MinimalBy[#, EuclideanDistance @@ # &]] &,
Flatten[{var1, var2}, List /@ {2, 3, 4, 1, 5}], {2}]];
graph =
CompleteGraph[Length[p],
EdgeWeight -> EuclideanDistance @@@ nearePoint];
Join[p,
nearePoint[[EdgeIndex[graph, #] & /@
EdgeList[FindSpanningTree[graph]]]]]]

ListLinePlot[ConnectComponentPoints[dat2]]


http://o8aucf9ny.bkt.clouddn.com/2017-02-18-20-30-50.png

### If you have a non-gap data like dat:

ConnectLines[dat_] :=
Module[{g =
SimpleGraph[RelationGraph[IntersectingQ, dat],
VertexLabels -> "Index"], path},
path = FindShortestPath[g, Sequence @@ GraphPeriphery[g]];
Append[First /@ #, #[[-1, -1]]] &[
FoldPairList[(Transpose[
DeleteDuplicates[
SortBy[Tuples[{##}], N[EuclideanDistance @@ #] &],
ContainsAny]] // {Reverse[#], #2} & @@ # &) &, First[path],
Rest[path], Identity]]]

newdat = ConnectLines[dat];
ListLinePlot[Join @@ newdat, Frame -> True]
Graphics[Arrow@newdat, Frame -> True]


http://o8aucf9ny.bkt.clouddn.com/2017-02-18-20-28-53.png

• It seems I do not have RelationGraph in Mathematica 10.1. Based on the appearance of IntersectingQ I am guessing that this code does not handle the case with gaps, dat2, correct? Feb 16, 2017 at 1:38
• @Mr.Wizard Yes,It cannot.
– yode
Feb 16, 2017 at 2:01
• @Mr.Wizard But your this post is about no gap case or I misunderstand?
– yode
Feb 16, 2017 at 2:57
• You mean my self-answer? It does work on the case with gaps as I showed. As noted in the question I also need to allow for gaps between lines. End points will not always be as close as in this example. Nevertheless I will vote for your answer for showing me RelationGraph`. Feb 16, 2017 at 6:24
• @Mr.Wizard In your gap case,do you wanna leave the maximal gap as your self-answer or whatever it is?
– yode
Feb 16, 2017 at 8:19