Determine the location of max curvature for a set of data

In order to find the max curvature, I am using polyfit on my data with which I can calculate the derivatives and therefor the curvature. Unfortunately it is not precise enough and so the max curvature is not right there, where it is supposed to be. I tried higher order, diff, and other fits, however I couldn't find a satisfying solution. Which bugs me quite a bit, because it seems to me like a simple task (at least i thought so). You can easily determine the point of max curvature by just looking at the plot. At the Picture you can see, that the max curvature is at around 7.5, but it should be at 11.

Maybe somebody can Point me in the right direction.

Answer 1

[John D'Errico]

If I had a nickel for every time someone asked a similar question, I won't say I'd be rich, but at least I'd get tired of rolling up those nickels and carrying them to the bank.

A problem is that it is easy to see something in your mind. You say, thats what I want to see. But doing the computation can be sometimes less easy. Computing the derivative(s) of a noisy function, especially when the data is not equally spaced is what is called an ill-posed problem. It amplifies any noise in the data. And that amplification can be significant. Finding the location of a second derivative max for noisy data can be nasty as hell to do.

So what would I suggest?

Don't use polyfit!!!!!!! Using a high order polynomial here is absolutely insane. Raising the order of the polynomial is insanity raised to a power. Run as fast as you can, away from polyfit!

I don't have your data, so I can't give you much of an example. In general, a far better choice will be a spline model. If you attach your data in a .mat file to a comment, I can give you a decent example of how I would approach the problem.

x = linspace(-5,5,100);
y = exp(-x.^2) + sin(x/2);
plot(x,y,'o'),grid on

This curve has NO noise in it at all. So I can use a simple interpolating spline to fit it.

pp = spline(x,y);

Now, find the location of the max and min of the second derivative of the curve. I'll use my slmpar tool, as found in my SLM toolbox, on the File Exchange.

[fppmax,fppmaxloc] = slmpar(pp,'maxfpp')
fppmax =
       1.0405
fppmaxloc =
      -1.2626

[fppmin,fppminloc] = slmpar(pp,'minfpp')
fppmin =
      -2.0011
fppminloc =
     0.050505

In fact, I'd probably not have picked that spot for the maximum of the second derivative by eye, but that location does make sense.

Dp you want to see the second derivative function plotted? This will give you a hint as to whether it was successful.

If my data was noisy, but equally spaced in x (as it is here, but with noise added on) then I would use either a smoothing spline of some sort (SLM will do this nicely) or I would use a Savitsky-Golay filter.

If the data is unequally spaced AND noisy, then a smoothing or least squares spline of some sort is your only viable choice.

[Kaindl]

The smoothing spline + equal spacing does the trick for me.

 % Equal spacing
xEqualSpaced = linspace(x(1), x(end), length(x))';
y = interp1(x, y, xEqualSpaced);
% Smoothing spline
[pp, ~] = csaps( xEqualSpaced, y ,0.9 );

The 2nd derivative with fnder

p_der2 = fnder(pp,2);

Thanks for the support.

[Bruce Liu]

Thanks for sharing,

How to plot the second figure?

[Cong Liu]

Thank you for sharing.

Why do you use fnder to find the derivative? Can you use ODE45? what is the differences between these two?

Answer 2

[Image Analyst]

For what it's worth, take a look at Dirk's file Exchange:

http://www.mathworks.com/matlabcentral/fileexchange/32696-2d-line-curvature-and-normals

Answer 3

[Mads]

Hello

I can't seem to get neither the SLM toolbox or Dirks stuff to locate the corner of this L-curve. Can anyone help?

Here are the data points

eta =
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2960e-003
     1.2959e-003
     1.2957e-003
     1.2954e-003
     1.2947e-003
     1.2932e-003
     1.2900e-003
     1.2833e-003
     1.2693e-003
     1.2409e-003
     1.1848e-003
     1.0829e-003
   920.7605e-006
   712.8970e-006
   511.8961e-006
   366.8875e-006
   280.2828e-006
   227.0583e-006
   189.2376e-006
   160.2095e-006
   136.8679e-006
   116.7839e-006
    99.4034e-006
    84.8592e-006
    72.8964e-006
    63.6574e-006
    57.0217e-006
    52.1719e-006
    48.3709e-006
    44.9284e-006
    40.8770e-006
    35.4289e-006

rho =
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1455e-003
     3.1456e-003
     3.1457e-003
     3.1462e-003
     3.1476e-003
     3.1506e-003
     3.1552e-003
     3.1602e-003
     3.1648e-003
     3.1696e-003
     3.1756e-003
     3.1838e-003
     3.1957e-003
     3.2144e-003
     3.2435e-003
     3.2874e-003
     3.3530e-003
     3.4455e-003
     3.5701e-003
     3.7457e-003
     4.0160e-003
     4.5048e-003
     5.6426e-003
     8.5212e-003

[Image Analyst]

How about looking at distances of the curve from its asymptotes? Try this:

format long g;
format compact;
fontSize = 15;

eta =[...
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2960e-003
  1.2959e-003
  1.2957e-003
  1.2954e-003
  1.2947e-003
  1.2932e-003
  1.2900e-003
  1.2833e-003
  1.2693e-003
  1.2409e-003
  1.1848e-003
  1.0829e-003
  920.7605e-006
  712.8970e-006
  511.8961e-006
  366.8875e-006
  280.2828e-006
  227.0583e-006
  189.2376e-006
  160.2095e-006
  136.8679e-006
  116.7839e-006
  99.4034e-006
  84.8592e-006
  72.8964e-006
  63.6574e-006
  57.0217e-006
  52.1719e-006
  48.3709e-006
  44.9284e-006
  40.8770e-006
  35.4289e-006]

rho =[...
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1455e-003
  3.1456e-003
  3.1457e-003
  3.1462e-003
  3.1476e-003
  3.1506e-003
  3.1552e-003
  3.1602e-003
  3.1648e-003
  3.1696e-003
  3.1756e-003
  3.1838e-003
  3.1957e-003
  3.2144e-003
  3.2435e-003
  3.2874e-003
  3.3530e-003
  3.4455e-003
  3.5701e-003
  3.7457e-003
  4.0160e-003
  4.5048e-003
  5.6426e-003
  8.5212e-003]
plot(eta, rho, 'b-', 'LineWidth', 2);
grid on;
axis equal % Make axes equal, otherwise "corner" depends on scaling.
xlabel('eta', 'FontSize', fontSize);
ylabel('rho', 'FontSize', fontSize);

cornerX = min(eta);
cornerY = min(rho);
distancesFromCorner = sqrt((eta - cornerX) .^ 2 + (rho - cornerY) .^ 2);
% Find min distance
[minDistance, indexOfMin] = min(distancesFromCorner);
minEta = eta(indexOfMin)
minRho = rho(indexOfMin)
% Start x axis at 0
xlim([0, max(eta)]);
% Start y axis at 3
ylim([0.003, max(rho)]);
% Get location of axes:
xl = xlim;
yl = ylim;
% Plot a vertical line there.
line([minEta, minEta], [yl(1), minRho], 'Color', 'r', 'LineWidth', 2);
% Plot a horzontal line there.
line([xl(1), minEta], [minRho, minRho], 'Color', 'r', 'LineWidth', 2);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);

[Mads]

Awesome, beautiful. Thanks. If I could set this as "Accepted" I would.