# How to group and join points

It seems there are many similar questions on Mathematica.SE, but none addresses (to the best of google) exactly my problem. I have the following data:

L[1] =
{{0.0200, {-3.98, -0.902}}, {0.120, {-3.97, -0.890}}, {0.220, {-3.94, -0.860}},
{0.320, {-3.89, -0.811}}, {0.420, {-3.82, -0.744}},{0.520, {-3.71, -0.659}},
{0.620, {-3.55, -0.554}}, {0.720, {-3.32, -0.429}}, {0.820, {-2.97, -0.285}},
{0.920, {-0.119}}, {1.02, {0.0666}}, {1.12, {0.272}}, {1.22, {0.497}}, {1.32, {0.740}},
{1.42,{1.00, 3.30}}, {1.52, {1.28, 3.68}}, {1.62, {1.58, 4.04}}, {1.72, {1.89, 4.41}},
{1.82, {2.23, 4.79}}, {1.92, {2.57,5.18}}, {2.02, {2.90, 5.58}}, {2.12, {3.38,6.00}},
{2.22, {3.90, 6.43}}, {2.32, {4.44, 6.88}}, {2.42, {5.00,7.35}}, {2.52, {5.57, 7.84}},
{2.62, {6.15,8.35}}, {2.72, {6.75, 8.87}}, {2.82, {7.36, 9.42}}, {2.92, {7.99,9.98}}};
L[2] =
{{0.0200, {0.377, 0.579}}, {0.120, {0.373,0.580}}, {0.220, {0.360, 0.582}},
{0.320, {0.340, 0.585}}, {0.420, {0.334, 0.590}}, {0.520, {0.331,0.596}},
{0.620, {0.336, 0.604}}, {0.720, {0.369,0.613}}, {0.820, {0.501,0.623}},
{0.920,{0.632}}, {1.02, {0.638}}, {1.12, {0.638}},{1.22, {0.635}}, {1.32, {0.629}},
{1.42, {0.622,0.974}}, {1.52, {0.614, 0.944}}, {1.62, {0.601, 0.933}},
{1.72, {0.581, 0.929}}, {1.82, {0.545, 0.933}}, {1.92, {0.472, 0.937}},
{2.02, {0.400, 0.936}}, {2.12, {0.864, 0.934}}, {2.22, {0.889, 0.941}},
{2.32, {0.903, 0.945}}, {2.42, {0.916, 0.949}}, {2.52, {0.929, 0.958}},
{2.62, {0.934, 0.958}}, {2.72, {0.935, 0.949}}, {2.82, {0.943, 0.960}},
{2.92, {0.941, 0.955}}};
L[3] =
{{0.0200, {0, 0.0669}}, {0.120, {0, 0.0649}}, {0.220, {0,0.0601}},
{0.320, {0.00730, 0.0528}}, {0.420, {0.0350, 0.0439}},
{0.520, {0.0816, 0.0339}},{0.620, {0.155,0.0237}},
{0.720, {0.297, 0.0141}}, {0.820, {0.727, 0.00620}}, {0.920, {0.00110}},
{1.02, {0.000300}}, {1.12,{0.00460}}, {1.22, {0.0130}}, {1.32, {0.0247}},
{1.42, {0.0395,2.93}}, {1.52, {0.0578, 2.57}}, {1.62, {0.0776, 2.38}},
{1.72, {0.0977, 2.26}}, {1.82, {0.115, 2.18}}, {1.92, {0.123, 2.12}},
{2.02, {0.135, 2.06}}, {2.12, {0.543, 2.01}}, {2.22, {0.697, 1.97}},
{2.32, {0.806,1.94}}, {2.42, {0.891,1.91}}, {2.52, {0.961, 1.90}}, {2.62,
{1.01,1.86}}, {2.72, {1.05, 1.82}}, {2.82, {1.08, 1.81}}, {2.92, {1.10,1.77}}};


It is having the following structure

Row[MatrixForm /@ {L[1], L[2], L[3]}]


That is we have $3$ ragged arrays. Let us plot them:

Do[
l[i] = Flatten[Table[Table[{a[[1]], b}, {b, a[[2]]}], {a, L[i]}], 1];
gr[i] =
ListPlot[
l[i], PlotRange -> {{0, 3}, All},
PlotMarkers -> {Graphics@{Black, Disk[]}, 0.02},
FrameLabel -> {"x", "\[CapitalDelta][" <>ToString[i] <> "]"},
Axes -> False, Frame -> True,
LabelStyle -> Directive[FontFamily -> "Helvetica", FontSize -> 14],
PlotStyle -> Black],
{i, 3}]
GraphicsRow[{gr[1], gr[2], gr[3]}, ImageSize -> Full]


These three panels describe 3 different properties of a physical system as a function of one parameter. On each panel there are $3$ curves.

My question is how to connect points belonging to each of $3$ curves. It is certainly possible to edit plots in e.g. adobe illustrator however I am more curios to properly extract the data corresponding to each curve so the parametric plots, e.g., $\Delta[2] (\Delta[1])$ can be created for each curve.

Just to avoid suggestions like Joined -> True or ListLinePlot[], this is not the expected result:

Reformulation

Imagine you observe 3 animals in the forrest (rabbit, wolf, fox) as a function of time $(x)$. They have 3 properties (L[1] - smell, L[2] - sound, L[3] - look). Not all animals are close enough to the detector, sometimes detectors register 1, sometimes 2 animals at a time.

The question is:

• (i) to separate signals belonging to the rabbit, wolf and fox and
• (ii) to plot the "sound" as a function of "smell" for each of them.

Clarification

As you can see, second columns of data L[1], L[2], and L[3] are lists. In fact they are ordered lists, and the ordering is the same for L[1], L[2], and L[3]. For instance, have a look at the first raw:

{-3.98, -0.902} {0.377, 0.579} {0, 0.0669}


It is known that {-3.98, 0.377, 0} belong to one group and { -0.902, 0.579, 0.0669} to another.

Expected result

There are several very nice answers below. I really appreciate! In order to clarify what exactly I am expecting I re-colored points on the plots manually.

I think the right approach would be to use information from all sets.

• @Coolwater, please read the question carefully Feb 9, 2016 at 12:12
• @Edmund, please read the question carefully Feb 9, 2016 at 12:12
• So I can see three distinct curves in L[1], and I could imagine adapting this answer to find them, but for L[2] and L[3] I do not see three distinct curves. Feb 9, 2016 at 12:44
• Let me see if I am understanding your problem: you're measuring 3 different properties of a physical system, but all your three measurements get mixed up in one bag and you want to separate them afterwards. I m really curious about how comes you can't separate them beforehand. Feb 9, 2016 at 13:06
• @Dr.belisarius, no it is not exactly how you reformulate it. Sorry for figurative language below. Imagine you observe 3 animals in the forrest (rabbit, wolf, fox) as a function of time (x). They have 3 properties (L[1] - smell, L[2] - sound, L[3] -look). Not all animals are close enough to the detector, sometimes detectors register 1, sometimes 2 animals at a time. The question is (i) to separate signals belonging to rabbit, wolf and fox and (ii) plot sound as a function of smell for each of them. Feb 9, 2016 at 13:18

This approach splits each data set in to a set of curves and then attempts to join curves whose end and start points are "close enough". The measure I have used is okay for the example data and it includes an element of rescaling with the curve data but YMMV with "real" data. It is simpler than some of the linked approaches.

Firstly we use a helper function to split the raw data. It enables us to group data points according to their position in the list for each ordinate point.

sowByPosn[{x_, ys_}] := MapIndexed[Sow[{x, #1}, #2[[1]]] &, ys]


We now split the data according to the number of points for each ordinate, then generate and group data points in each sub-list according to their position:

Do[l[i] =
SplitBy[L[i], Length[#[[2]]] &]);
Print@ListLinePlot[l[i], PlotTheme -> "Detailed",
PlotLegends -> Automatic], {i, 3}]


There are five curves in each data set and three obviously correspond to a single curve.

To join the curves we compare the distances between the ends of each pair of lines to a limit value calculated from each curve.

The limit is calculated from the mean difference plus 5 times the standard deviation (somewhat arbitrarily):

stdLimit[x_] :=
Plus @@ ({1,
5} Through[{Mean, StandardDeviation}[Norm /@ Differences[x]]])


Generate a table of differences (with mean limits) for each pair of curves and select pairs that have ends closer than the limit.

Table[p[i] =
Position[Table[{Norm[Last[c1] - First[c2]],
Mean[stdLimit /@ {c1, c2}]}, {c1, l[i]}, {c2, l[i]}], {x_, l_} /;
x <= l], {i, 3}]


{{{2, 3}, {3, 4}}, {{2, 3}, {3, 4}}, {{2, 3}, {3, 4}}}

Use Mathematica's Graph functionality to group the components (probably we could do this without...) We only need to do this once with the example data -- all have the same connectivity.

components =
With[{g =
Graph[Join[# -> # & /@
Complement[Range[Length[l[1]]], Union[Flatten[p[1]]]],
Rule @@@ p[1]], VertexLabels -> "Name"]}, Print[g];
(Flatten /@ (EdgeList[Subgraph[g, #]] & /@
WeaklyConnectedComponents[g] /.
DirectedEdge[u_, v_] :> {u, v})) /. {a___, b_, b_, c___} :> {a,
b, c}]


{{2, 3, 4}, {5}, {1}}

Finally join the lines:

Table[z[i] = Join @@ Part[l[i], #] & /@ components;
Print@ListLinePlot[z[i], PlotTheme -> "Detailed",
PlotLegends -> ToString /@ components], {i, 3}];


• I think your result is the most accurate one. Thanks a lot! I need to read your algorithms in order to understand if it is applicable to a larger class of problems. Feb 9, 2016 at 18:47
• Thanks, I can see there may be problems if eg two curves end close to a beginning. In that case I think perhaps looking at the gradient may be the way to go. Another minor point is that it reorders the curves. Feb 9, 2016 at 18:58

We can use the following approach: start with an end point that belongs to a path and increment the path with nearest neighbor points that are "good candidates." A point is a "good candidate" if it is not too far away from the last point and it does not produce a sharp turn in the path. To determine "too far" we a look at the distances between the last, say, 5 point pairs.

In order this approach to work the data has to be rescaled and the paths initialized.

This answer is more-or-less complete. It can be refined further with automated rescaling and path initilization.

Below is code that implements the outlined approach.

Clear[GoodPointQ]
distNLastPoints_: 6] :=
Block[{dists, threshold, cd, v1, v2},
(*Path is too short.*)
If[Length[pointPath] < 3, Return[True]];

(*Is the distance an outlier?*)

dists = Map[Norm[#[[1]] - #[[2]]] &,
Partition[Take[Reverse[pointPath], UpTo[distNLastPoints]], 2, 1]];
(*threshold=Median[dists]+(5*Median[Abs[dists-Median[dists]]]);*)

threshold =
(*Print["norm:",Norm[
pointPath\[LeftDoubleBracket]-1\[RightDoubleBracket]-candidate]];
Print["threshold:",threshold];*)

If[Length[pointPath] > 4 &&
Norm[pointPath[[-1]] - candidate] > threshold,
Return[False]
];

(*Does it make a sharp turn?*)

v1 = pointPath[[-2]] - pointPath[[-1]];
v2 = candidate - pointPath[[-1]];
cd = Dot[v1/Norm[v1], v2/Norm[v2]];
(*Print["cd:",cd];*)
If[cd > -10^-4, Return[False]];
True
];

Clear[FindNextGoodPoint]
FindNextGoodPoint[pointPath_, nnFunc_, nCandidates_: 3,
sameDist_: 10^-9, distOutliersSDFactor_: 10, distNLastPoints_: 6] :=

Block[{cs, good},
cs = nnFunc[pointPath[[-1]], nCandidates + 1];
cs = Complement[cs, pointPath,
SameTest -> (Norm[#1 - #2] <= sameDist &)];
good =
distNLastPoints] & /@ cs;
If[Or @@ good,
{True, Append[pointPath, Pick[cs, good][[1]]]},
{False, pointPath}]
];

Clear[PathByNextGoodPoint]
PathByNextGoodPoint[startPoints : {{_?NumberQ, _?NumberQ} ..},
allPoints : {{_?NumberQ, _?NumberQ} ..}, nCandidates_: 3,

Block[{pointPath, found = True, nnFunc, k = 0},
nnFunc = Nearest[allPoints];
pointPath = startPoints;
While[found && k < 200,
k++;
{found, pointPath} =
FindNextGoodPoint[pointPath, nnFunc, nCandidates];
];
pointPath
];


Steps of applying this code follow.

Select a dataset.

data = dataOrig = Flatten[Thread /@ L[1], 1];


Normalize the data. (Standardize can be used instead of Rescale.)

Do[
data[[All, i]] =
Rescale[data[[All, i]], MinMax[data[[All, i]]], {0, 1}]
, {i, {1, 2}}]


This is useful to select paths initial points.

(*Graphics[{Point[data],MapIndexed[Text[#2,#1,{-1,1}]&,data]}]*)


Initial paths of points.

startPoints = {{data[[1]], data[[3]]}, {data[[2]],
data[[4]]}, {data[[55]], data[[53]]}};


Find the paths.

paths = PathByNextGoodPoint[#, data, 2] & /@ startPoints;


Single points can be used too.

startPoints = {{data[[1]]}, {data[[2]]}, {data[[55]]}};
paths = PathByNextGoodPoint[#, data, 2] & /@ startPoints;


Rescale to the original data.

pathsOrig = paths;
Do[
pathsOrig[[j, All, i]] =
Rescale[pathsOrig[[j, All, i]], {0, 1},
MinMax[dataOrig[[All, i]]]], {j, Length[paths]}, {i, {1, 2}}];


Plot results.

Show[{ListPlot[dataOrig,
PlotStyle -> {PointSize[0.02], GrayLevel[0.8]}],
ListPlot[{Sequence @@ pathsOrig}]}, Frame -> True]


Results with the other two sets:

The paths on last image can be somewhat improved by tweaking the parameters of the function PathByNextGoodPoint.

• I am impressed that you could get accurate results working with each set independently. Do you see a possibility to apply your method to triplets of data points? I expect it can improve the accuracy. Feb 9, 2016 at 18:52
• @yarchik I need to think more about the problem and my solution. I think the other solutions are much better. As you indicate in your comment, my code is a little too ad-hoc and low level. Feb 9, 2016 at 21:24

As David G. Stork has pointed out, Nearest Neighbours offer a good method of attack for this problem. Here I've not implemented a full NN chain approach but something a little more basic which gets most of the way there.

I'm using a 'dumb' NN but with a distance function that only allows points to be connected to another point one x-distance away, this works because your data are on a regular grid.

First I just pad the data into a full {x,y} format:

    p1 = Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@
L[1]), 1];
p2 = Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@
L[2]), 1];
p3 = Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@
L[3]), 1];


Then using the following distance function:

   distM[{u_, v_}, {x_, y_}] := With[{pen = 10},
If[Round[10 (u - x)] == -1, Abs[v - y], pen]
];
distP[{u_, v_}, {x_, y_}] := With[{pen = 10},
If[Round[10 (u - x)] == 1, Abs[v - y], pen]
];

nf1 = Nearest[p1, DistanceFunction -> distM];

Graphics[{Point@p1, Line /@ Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@Transpose[{p1, (nf1[#, {1, 1}] & /@ p1)}]), 1]}, opts]


   nf2 = Nearest[p2, DistanceFunction -> distM];

Graphics[{Point@p2, Line /@ Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@Transpose[{p2, (nf2[#, {1, 0.15}] & /@ p2)}]), 1]}, opts]


   nf3 = Nearest[p3, DistanceFunction -> distP];

Graphics[{Point@p3, Line /@ Flatten[Transpose /@ ({ConstantArray[#[[1]], Length[#[[2]]]], #[[2]]} & /@Transpose[{p3, (nf3[#, {1, 1}] & /@ p3)}]), 1]}, opts]


   opts = Sequence[Frame -> True, AspectRatio -> 1/GoldenRatio,


GridLines -> Automatic, PlotRangePadding -> None];

Imperfect but might serve as a starting point and hopefully shows a few of the possibilities.

• Very nice! (+1) If you posted it 1h earlier, I would have not spent time to program the code for my answer. Feb 9, 2016 at 18:10
• I think your code is the shortest one. It is also pretty accurate. Here I have the same question as for the approach of @Anton Antonov, can it be applied to the triples of data-points in order to improve the accuracy? Feb 9, 2016 at 18:55