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!
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.