4
$\begingroup$

TL;DR

How to implement in mathematica a tool such as WebPlotDigitizer by Ankit Rohatgi which would allow us to trace automatically curves from either (randomly ordered) data or imported images.

enter image description here

Context

As a Follow up of this question and both nice answers given there, I am after an algorithm which could (as automatically as possible) trace the different sets of curves in plots such as these (calorific curves of distribution of intermediate black holes in a Galactic centre).

enter image description here

The challenge is that these points are found by some complex optimisation routine and are over-numerous in places. So I am interested in resampling the different curves. Note importantly that the curves cross (e.g. near (-0.015,2)). This can be seen as a reverse engineering problem to extract a compact re-parametrisation of the different curves (hence the link to WebPlotDigitizer).

Question

How can one achieve the automatic curve matching/ sampling procedure?


Attempt

Let me define a toy problem as follows: let me produce two curves:

data = Flatten[{Table[{x, Sin[x^2/15]}, {x, 0, 20, 0.01}],
    Table[{x, 1 + Cos[x]}, {x, 0, 20, 0.1}]}, 1];

and draw random points from those.

pts = RandomChoice[data, Length[data]];

ListPlot[pts]

enter image description here

Note that on purpose the sampling is not the same for the two curves. One could imaging one of two situations: either have access to pts or ListPlot[pts]//Image (but obviously not data).

Following this (excellent) answer I can get a good resampling as follows:

tour = Last[FindShortestTour[pts]];
tourPts = Extract[pts, List /@ tour];
peaks = Ordering[EuclideanDistance @@@ Partition[tourPts, 2, 1], -2];
{firstCurve, secondCurve} = 
  TakeDrop[RotateLeft[tourPts, peaks[[1]]], Abs[Subtract @@ peaks]];


np[f_][u_, dt_] := u + dt/Norm[f'[u]]
equallySpacedPts[pts_, dt_] := 
 With[{bsf = BSplineFunction[pts]}, 
  bsf /@ Most[NestWhileList[np[bsf][#, dt] &, 0, # < 1 &]]]

equallySpacedPts[#,0.25]&/@{firstCurve,secondCurve}//ListLinePlot

enter image description here

BUT It is clear that the two curves are not (always) properly matched, e.g. near x=3 or x=15.

I understand that this is not a trivial matter, but it should be IMHO of general interest for mathematica to be able to stand up to this challenge with minimum manual input (?).


Comment:

It might be possible to use a R package via digitizeR and R integration in Mathematica, but obviously a standalone implementation would be preferable.

$\endgroup$
8
  • 1
    $\begingroup$ Take a look at this answer for (maybe) some helpful clues. mathematica.stackexchange.com/questions/149244/… $\endgroup$
    – MikeY
    Jul 12, 2021 at 16:20
  • 1
    $\begingroup$ Also this one, which uses the previous one's result. mathematica.stackexchange.com/questions/146193/… $\endgroup$
    – MikeY
    Jul 12, 2021 at 16:22
  • $\begingroup$ @MikeY sounds promising indeed. Would you care writing an answer? $\endgroup$
    – chris
    Jul 12, 2021 at 16:28
  • $\begingroup$ Is the big plot a single curve that spirals around, or two curves spiraling, or a set of individual ellipses? Checking... $\endgroup$
    – MikeY
    Jul 12, 2021 at 18:15
  • $\begingroup$ @MikeY the initial plot is only illustrative: it is a 3D set of tilted circles $\endgroup$
    – chris
    Jul 12, 2021 at 18:20

1 Answer 1

1
$\begingroup$

Here's one approach that seems to work reasonably well:

cleaned = DeleteDuplicates@SortBy[First]@pts;
curves = {};
(pt \[Function] Module[
     {best},
     best = MinimalBy[#[[2]] &]@MapIndexed[{#2[[1]], VectorAngle[(#[[-1]] - pt), Subtract @@ #[[-2 ;;]]], EuclideanDistance[pt, #[[-1]]]} &]@curves;
     If[best === {} || best[[1, 3]] > 1,
      AppendTo[curves, {pt - {0.1, 0}, pt}],
      AppendTo[curves[[best[[1, 1]]]], pt]
      ]
     ]) /@ cleaned;
curves = Rest /@ curves;

ListPlot@curves

enter image description here

We first sort the points by their x-coordinate and remove any duplicates. We now initialize a list of collected curves (initially empty). Then, we go through the points one by one and find the curve where adding the new point would cause the smallest change in direction (where "change in direction" is computed using VectorAngle). The idea here is that we'd like to avoid sharp corners in the lines if possible. If curves was empty, or the existing lines are all too far away (this threshold might need to be adjusted depending on the curves at hand), we create a new entry in curves. We add two points to make sure the angle computation doesn't fail (this point is later removed). Otherwise, we add the point to the best match.

The biggest issue currently is the detection of new lines: The current method based on the distance to existing curves might not always work as expected if the points are spaced too far apart on some of the curves.

$\endgroup$
2
  • $\begingroup$ Nice! Thank you for this elegant and short answer. I will let others a chance to give alternative answers and hope to motivate WRI to integrate a generic solution in upcoming release of Mathematica which would possibly integrate some of the code provided by WebPlotDigitizer (?). $\endgroup$
    – chris
    Jul 15, 2021 at 9:49
  • $\begingroup$ In other words, I am personally interested in your answer but I think being able to do this from images should prove generally useful. $\endgroup$
    – chris
    Jul 15, 2021 at 9:59

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.