How do I calculate radius of curvature from discrete samples?

Question

I've got a series of (x,y) samples from a plane curve, taken from real measurements, so presumably a bit noisy and not evenly spaced in time.

x = -2.51509  -2.38485  -1.88485  -1.38485  -0.88485  -0.38485   0.11515 0.61515   1.11515   1.61515 ...

y =  -48.902  -48.917  -48.955  -48.981  -49.001  -49.014  -49.015  -49.010 -49.001  -48.974 ...

If I plot the whole series, it looks like a nice oval, but if I look closely, the line looks a bit wiggly, which is presumably the noise.

How would I go about extracting an estimate of the radius of curvature of the underlying oval?

Any programming language would be fine!

Answer 1

Roger Stafford gave some MATLAB code here:

http://www.mathworks.com/matlabcentral/newsreader/view_thread/152405

Which I verbosed a bit to make this function:

# given a load of points, with x,y coordinates, we can estimate the radius
# of curvature by fitting a circle to them using least squares.
function [r,a,b]=radiusofcurv(x,y)
  # translate the points to the centre of mass coordinates
  mx = mean(x);
  my = mean(y);
  X = x - mx; Y = y - my; 

  dx2 = mean(X.^2);
  dy2 = mean(Y.^2);

  # Set up linear equation for derivative and solve
  RHS=(X.^2-dx2+Y.^2-dy2)/2; 
  M=[X,Y];
  t = M\RHS;

  # t is the centre of the circle [a0;b0]
  a0 = t(1); b0 = t(2);

  # from which we can get the radius
  r = sqrt(dx2+dy2+a0^2+b0^2); 

  # return to given coordinate system
  a = a0 + mx;
  b = b0 + my; 

endfunction

It seems to work quite well for my purposes, although it gives very strange answers for e.g. collinear points. But if they're from a nice curve with a bit of noise added, job pretty much done.

It should adapt pretty easily to other languages, but note that \ is MATLAB/Octave's 'solve using pseudo-inverse' function, so you'll need a linear algebra library that can calculate the pseudo inverse to replicate it.

Answer 2

I think you can calculate the ellipse which satisfies the least-square of the difference between your data and ellipse equation. Then use the major and minor axis of the ellipse. Here is a paper I found after a Google search.