Matlab code for finding the curvature of a curve using given data points

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.

• Curvature is not the same as variance. Also it is likely not very constant for the whole curve. There is one very easy way to check it anyway. Construct short pieces of circles and run the code on them. Oct 16 '15 at 14:23
• Yes I do that @mathreadler, and I check the matlab code, It works nicely. However I dont understand the idea used in this code to compute the curvature. Oct 16 '15 at 14:26
• Can you do least squares for a straight line fit? What difference do you find between a circle and a straight line modeling? Oct 21 '15 at 14:16

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:

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.

• Thank you for your answer, What i did really is that I wrote another matlab code that find the center and radius of the oscillating circle ( I take the number of given points as variable greater than three), I got good results that exactly mathch the results of the above code> But really I am very curios t know how this code works?? why using the mean and the variance ??!! Oct 19 '15 at 20:12
• @Nizar: Without an accompanying documentation, not being an expert in Matlab, it's virtually impossible for me to do any reverse engineering about what theory is behind that code. However, I you read the PDF document mentioned in my answer, then you will see that mean $(\mu_x,\mu_y)$ and variance $(\sigma_{xx},\sigma_{xy},\sigma_{yy})$ come with the RI method in this paper. I have no idea if that may be comparable with the theory as implemented by your Matlab program. Oct 20 '15 at 13:07
• Thak you again for your help, but sorry I have a small question in the pdf you suggested for me : What is the role of $w_i$ in the least square formula? where it come from? what did it represent?? Oct 22 '15 at 6:39
• @Nizar: The $w_i$ represent weighting coefficients; they can be conveniently taken $=1$ here. (Sorry I didn't explain them in the paper; that happens if there are no referees to check content) Oct 22 '15 at 11:20