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

*No products are associated with this question.*

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.

Jan Simon
on 16 Jan 2013

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

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.

Jan Simon
on 16 Jan 2013

Opportunities for recent engineering grads.

## 1 Comment

## Matt J (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/58964#comment_122940

Your function looks 1D to me.