curvature of a discrete function
82 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentar
Akzeptierte Antwort
Jan
am 16 Jan. 2013
Bearbeitet: Jan
am 13 Mär. 2022
Your 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) .^ 3;
curvature = num ./ denom;
curvature(denom < 0) = NaN;
Please test this, because I'm not sure if I remember the formulas correctly.
3 Kommentare
Jan
am 16 Jan. 2013
Bearbeitet: Jan
am 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
am 16 Jan. 2013
Bearbeitet: Jan
am 13 Mär. 2022
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
Weitere Antworten (2)
Roger Stafford
am 16 Jan. 2013
Bearbeitet: Bruno Luong
am 13 Mär. 2022
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.
2 Kommentare
Jan
am 16 Jan. 2013
+1, I've waited for this answer. See http://www.mathworks.com/matlabcentral/newsreader/view_thread/152405.
Moreno, M.
am 13 Mär. 2022
Bearbeitet: Jan
am 13 Mär. 2022
Siehe auch
Kategorien
Mehr zu Author Block Masks finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!