I have data points $(x,y)$ for a plane curve, and I would like to find its curvature. While I was googling to check how could I start, I found this matlab code:
mx = mean(x); my = mean(y)
X = x - mx; Y = y - my; % Get differences from means
dx2 = mean(X.^2); dy2 = mean(Y.^2); % Get variances
t = [X,Y]\(X.^2-dx2+Y.^2-dy2)/2; % Solve least mean squares problem
a0 = t(1); b0 = t(2); % t is the 2 x 1 solution array [a0;b0]
r = sqrt(dx2+dy2+a0^2+b0^2); % Calculate the radius
a = a0 + mx; b = b0 + my; % Locate the circle's center
curv = 1/r; % Get the curvature
No specific explanation was added to this code except:
"The circle defined by center $(a, b)$ and radius $r$ will yield the least mean square value for the expression $(x-a)^2 + (y-b)^2 - r^2$ among all possible parameters, $a$, $b$, and $r$, over the points in vectors $x$ and $y$. The circle's curvature will be $1/r$."
The code worked very well, I checked this by taking as an example a circle.
I also tried to follow what they did step by step to understand the concept used behind this code; however I was stuck on the use of the mean. I know least square is used to solve systems, but why use mean least square? Why is the mean and variance included?
Would someone help me please. Any help will be highly appreciated.
The Matlab code may be explained by the fact that the curvature is one divided by the radius of an Osculating circle . Since all what you have is a bunch of data points, it has to be decided which data points may be relevant for the curvature at some place at the curve. At least three of them are needed to determine a circle. If more data points are taken, then a least squares procedure will be needed in order to determine a Best Fit Circle to them. For a better understanding of the latter, one might read the following document:
-
Best Fit Circles made Easy
There is no unambiguous recipe for how many data points to take for an osculating circle at some place at the curve. It all depends upon your data and all you can have is a decent approximation. Now I don't know how Matlab works, but it seems to me that these data points form the arrays $x$ and $y$ in the code. However, speaking for myself, I'd rather ignore this particular implementation and start from scratch, with my own understanding of the matter.Update. To demonstrate that Robert Israel's method (as explained in the reference) does indeed work:
The code itself is to large to fit into the margin. Serious: it really has no sense to publish all those (Delphi Pascal) programming details: I'm pretty sure that wouldn't add much to understanding.
But if someone asks me for the source code, then he or she can have it for free.