Matlab code for computing curvature equation

How can I compute with a code on Matlab the equation of the curvature k=(x' y"-y'x")/(x'^2+y'^2)^3/2 ? I have a relevant number of curves, that I divided into more parts and for each of them I measured the coordinates (x,y). I would like to calculate the curvature for each curve. I thank you very much in advance.

Answer 1

[Roger Stafford]

Your formula for curvature is that of a curve defined in terms of a parameter t in which x’, y’, x”, and y” all refer to derivatives with respect to that parameter. If you have only the x and y coordinates of points on a curve, that formula will not be very useful to you.

If you have confidence in the accuracy of your point coordinates, you can use the formula for the curvature of a circle through three adjacent points on your curve. That would be an approximation for the curvature at the middle point. That formula is: curvature equals four times the area of the triangle formed by the three points divided by the product of the triangle’s three side lengths. Suppose (x1,y1), (x2,y2), and (x3,y3) are the coordinates of three adjacent points on the curve. The computation can go as follows:

   a = sqrt((x1-x2)^2+(y1-y2)^2); % The three sides
   b = sqrt((x2-x3)^2+(y2-y3)^2);
   c = sqrt((x3-x1)^2+(y3-y1)^2);
   s = (a+b+c)/2;
   A = sqrt(s*(s-a)*(s-b)*(s-c)); % Area of triangle
   K = 4*A/(a*b*c); % Curvature of circumscribing circle

[Roger Stafford]

An alternative to that code is:

   a = sqrt((x1-x2)^2+(y1-y2)^2); % The three sides
   b = sqrt((x2-x3)^2+(y2-y3)^2);
   c = sqrt((x3-x1)^2+(y3-y1)^2);
   A = 1/2*abs((x1-x2)*(y3-y2)-(y1-y2)*(x3-x2)); % Area of triangle
   K = 4*A/(a*b*c); % Curvature of circumscribing circle

[Ipazia]

I thank you for your answer. In my data I have also the length of each curve and, since the curves are very short and approach the hypotenuse of the right triangle made by x and y, the tangent vector at the middle point of each curve is almost all the times parallel the the hypotenuse. So I can easily calculate the angle described by t with the x-axis. Is it enough to use that equation or should I use the other you wrote?

[rsnandi]

but K is coming complex in my case ....