Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

To resolve issues starting MATLAB on Mac OS X 10.10 (Yosemite) visit: http://www.mathworks.com/matlabcentral/answers/159016

curvature of a discrete function

Asked by David Kusnirak on 16 Jan 2013

Hello,

I need to compute a curvature of a simple 2D discrete function like this one:

x=1:0.5:20;
y=exp(x);

can anybody help how to do that? thanks

1 Comment

Matt J on 16 Jan 2013

Your function looks 1D to me.

David Kusnirak

Tags

Products

No products are associated with this question.

3 Answers

Answer by Jan Simon on 16 Jan 2013
Edited by Jan Simon on 16 Jan 2013
Accepted answer

You function seems to be a 1D function.

Are you looking for the 2nd derivative? While diff calculates the one-sided differential quotient, gradient uses the two-sided inside the interval:

gradient(gradient(y))

If you mean the curvature as reciprocal radius of the local fitting circle:

dx  = gradient(x);
ddx = gradient(dx);
dy  = gradient(y);
ddy = gradient(dy);
num   = dx .* ddy - ddx .* dy;
denom = dx .* dx + dy .* dy;
denom = sqrt(denom);
denom = denom * denom * denom;
curvatur = num ./ denom;
curvature(denom < 0) = NaN;

Please test this, because I'm not sure if I remember the formulas correctly.

3 Comments

Matt J on 16 Jan 2013

CONV might be better then. GRADIENT is kind of slow

conv(x,[.25 0 -.5 0 .25],'same')
Jan Simon on 16 Jan 2013

Therefore I'm using an efficient C-Mex function: FEX: DGradient, which is 10 to 20 times faster and handles unevenly spaced data more accurate.

Jan Simon on 16 Jan 2013
x = rand(1, 1e6);
tic; ddx = gradient(gradient(x)); toc
tic; ddx = DGradient(DGradient(x)); toc
tic; ddx = conv(x,[.25 0 -.5 0 .25],'same')); toc
 Elapsed time is 0.251547 seconds.
 Elapsed time is 0.025728 seconds.
 Elapsed time is 0.028208 seconds

Matlab 2009a/64, Core2Duo, Win7

Jan Simon
Answer by Roger Stafford on 16 Jan 2013

Let (x1,y1), (x2,y2), and (x3,y3) be three successive points on your curve. The curvature of a circle drawn through them is simply four times the area of the triangle formed by the three points divided by the product of its three sides. Using the coordinates of the points this is given by:

 K = 2*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1)) ./ ...
  sqrt(((x2-x1).^2+(y2-y1).^2)*((x3-x1).^2+(y3-y1).^2)*((x3-x2).^2+(y3-y2).^2));

You can consider this as an approximation to the curve's curvature at the middle point of the three points.

Roger Stafford
Answer by Matt J on 16 Jan 2013
Edited by Matt J on 16 Jan 2013
diff(y,2)./0.5^2

0 Comments

Matt J

Contact us