enter image description here

Random curve like the above, for example. The data source is an image file. I don't have an expression.

I want to start by curve fitting and then compute curvature. Is it possible?

  • 1
    $\begingroup$ Do you have the equation (or parametric equation) of that curve? $\endgroup$ Sep 25, 2015 at 1:22

1 Answer 1


I guess the first step would always be to find an ordered list of points along the middle of the curve. That I can help with:

First binarize and thin the image of the curve, so you get a 1-pixel wide white line:

img = Import["https://i.sstatic.net/fEf1i.jpg"];
bin = Thinning@ColorNegate[Binarize[img]]

enter image description here

Finding the white pixels in this image is easy:

pts = N[PixelValuePositions[bin, 1.]];    

Finding the "right" order of these points along the path you're looking for is a bit more tricky. I've tried FindCurvePath, with little success. The points aren't all on one line either, because Thinning produces slight artifacts at corners. So using an efficient graph algorithm is probably cleaner, anyway. For that, I convert my point set into a fully connected graph, with weights equal to the distance between each point pair:

dist = Outer[EuclideanDistance, pts, pts, 1];    
g = WeightedAdjacencyGraph[dist, DirectedEdges -> False];    

...and get the spanning tree of this graph - which is basically the shortest tree that spans this point set. For a curve like this, the tree will look something like this:

st = FindSpanningTree[g]

enter image description here

i.e. one very long line and a few short branches. The short branches are artifacts from the thinning algorithm, we're not interested in those. We're looking for the long "spine" of the tree. In other words, we're looking for the pair of leaf vertices that with the longest path in between.

Finding the leaf vertices (the endings of the tree) is simple:

endVertices = Position[VertexDegree[st], 1][[All, 1]];    
HighlightImage[bin, pts[[endVertices]]]

enter image description here

And finding the paths between all pairs of these vertices isn't hard, either:

allPaths = 
  Flatten[Outer[FindShortestPath[st, #1, #2] &, endVertices, 
    endVertices], 1];

longestPath = MaximalBy[allPaths, Length][[1]];

Graphics[{Red, Line[pts[[longestPath]]]}]

enter image description here

Reusing ideas and code from this answer of signal processing I can calculate the approximate curvature from these points directly:

curve = pts[[longestPath]];

σ = 25;

{{dx1, dy1}, {dx2, dy2}} = 
  Table[GaussianFilter[curve[[All, dim]], σ, der], {der, 
    2}, {dim, 2}];

κ = (dx1*dy2 - dx2*dy1)/((dx1^2 + dy1^2)^(3/2));
ListLinePlot[κ, PlotRange -> All, 
 PlotLabel -> "Curvature κ"]

And produce a nice-looking animation of the curvature and osculating circle:

 frames = Table[
      ListLinePlot[curve, ImageSize -> {400, 400}, 
       AspectRatio -> Automatic, 
       PlotLabel -> "Osculating circle r=1/\[Kappa]",
       Epilog -> Module[{r, pt, n}, (pt = curve[[i]];
          r = Clip[1/\[Kappa][[i]], {-1, 1}*10^4];
          n = Normalize[{dy1[[i]], -dx1[[i]]}];
          {Red, PointSize[Large], Point[pt], 
           Circle[pt - n*r, Abs[r]]})]], 
      ListLinePlot[\[Kappa], PlotRange -> All, 
       PlotLabel -> "Curvature \[Kappa]", ImageSize -> 400, 
       Epilog -> {Red, Line[{{i, 0}, {i, 10}}]}]
      }], {i, 1, Length[curve], 5}], i]

enter image description here

Instead of simple Gaussian smoothing, you could use numerical optimization to find the smoothest path that's less than one pixel from this one (see this answer).

Alternatively, you could try to approximate this point set with e.g. a NURBS curve, then calculate the curvature of that. But I don't think MMA has any built-in algorithms for this task, so you'd have to write that yourself.

  • $\begingroup$ Do you have any hints for this problem? In particular, are there any builtins I'm missing? $\endgroup$
    – Szabolcs
    Oct 2, 2015 at 12:19
  • $\begingroup$ You should try a FindShortestTour $\endgroup$
    – yode
    Jan 29, 2016 at 16:54
  • $\begingroup$ @yode: You mean instead of the spanning tree? But the shortest tour (or any tour) goes through all points, then returns to the beginning. I don't want to visit all points, I want to remove those little "branches" that Thinning produces. And I don't want my path the return to the beginning. Also, finding a the shortest tour is an NP-complete problem where only approximate solutions are realistically possible, while finding a spanning tree is cheap. $\endgroup$ Jan 31, 2016 at 9:46

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.