How can i calculate the length of curve?

Hi,

I have a curve which includes X (meter) and Y (meter) data. Is there any way to obtain the length of curve easily?

Thanks a lot,

X=[96.0741000000000,97.1940000000000,98.3139000000000,99.4338000000000,100.553700000000,101.673600000000,102.793500000000,103.913400000000,105.033300000000,106.153200000000,107.273100000000,108.393000000000,109.512900000000,110.632800000000,111.752700000000,112.872600000000,113.992500000000,115.112400000000,116.232300000000]

X = 1×19

96.0741 97.1940 98.3139 99.4338 100.5537 101.6736 102.7935 103.9134 105.0333 106.1532 107.2731 108.3930 109.5129 110.6328 111.7527 112.8726 113.9925 115.1124 116.2323

Y=[-4.13836296940031,-4.10455468315876,-4.05645426203322,-3.99617782198545,-3.92344322326347,-3.83385191481492,-3.73582865974161,-3.61740402741020,-3.49399064332423,-3.35231953224592,-3.20552503148528,-3.04892626846560,-2.88658570885772,-2.72091440408539,-2.55226630046971,-2.38425597793465,-2.21787687713447,-2.05656258174384,-1.89889800700337]

Y = 1×19

-4.1384 -4.1046 -4.0565 -3.9962 -3.9234 -3.8339 -3.7358 -3.6174 -3.4940 -3.3523 -3.2055 -3.0489 -2.8866 -2.7209 -2.5523 -2.3843 -2.2179 -2.0566 -1.8989

[Volcano]

Thanks for all answers...

Answer 1

[Stephen23]

A very simple approach is to download John D'Errico's excellent ARCLENGTH function:

https://www.mathworks.com/matlabcentral/fileexchange/34871-arclength

X = -1:.01:1;
Y = sqrt(1-X.^2);
L = arclength(X,Y,'spline')

L = 3.1416

L-pi

ans = 5.0768e-07

For the sample curve, this gives a more accurate solution.

[Stephen23]

Tested on the cases given here:

https://www.mathworks.com/matlabcentral/answers/1787410-how-can-i-calculate-the-length-of-curve#comment_2588390

Arbitrary spacing (L=pi/2):

X = logspace(-10,0,200);
Y = sqrt(1-X.^2);
L = arclength(X,Y,'spline')

L = 1.5709

L - pi/2

ans = 8.9196e-05

Vertical lines (L=10):

X = ones(1,200);
Y = logspace(-10,1,200);
L = arclength(X,Y,'spline')

L = 10.0000

L - 10

ans = -1.0000e-10

[VIGNESH BALAJI]

nice tool

[Alejandro Casares]

How do you find the length when you have the parametric equations of the curve?

[Torsten]

For the case of a curve in the plane:

arclength = integral_{s=a}^{s=b} sqrt(x'(s)^2+y'(s)^2) ds

Example: Half circle

syms s 
a = 0;
b = 1;
x = cos(pi*s);
y = sin(pi*s);
arclength = int(sqrt(diff(x,s)^2+diff(y,s)^2),s,a,b)

arclength = 

π

pi

ans = 3.1416

[Stephen23]

Some more comparisons with ARCLENGTH:

a = linspace(0, 2*pi, 1E+3);
X = cos(a);
Y = sin(a);

dX = gradient(X);
dY = gradient(Y);
Len = cumtrapz(hypot(dX,dY));
Len(end) - 2*pi

ans = -4.1393e-05

sum(hypot(dX,dY)) - 2*pi

ans = 0.0062

arclength(X,Y,'spline') - 2*pi

ans = -1.3272e-11

Also with unequal spacing:

X = logspace(-10,0,200);
Y = sqrt(1-X.^2);

dX = gradient(X);
dY = gradient(Y);
Len = cumtrapz(hypot(dX,dY));
Len(end) - pi/2

ans = -0.0093

sum(hypot(dX,dY)) - pi/2

ans = 0.2245

arclength(X,Y,'spline') - pi/2

ans = 8.9196e-05

And with vertical lines:

X = ones(1,200);
Y = logspace(-10,1,200);

dX = gradient(X);
dY = gradient(Y);
Len = cumtrapz(hypot(dX,dY));
Len(end) - 10

ans = -1.0000e-10

sum(hypot(dX,dY)) - 10

ans = 0.5976

arclength(X,Y,'spline') - 10

ans = -1.0000e-10

See also: https://www.mathworks.com/matlabcentral/answers/1787410-how-can-i-calculate-the-length-of-curve#comment_2588390

https://www.mathworks.com/matlabcentral/answers/1787410-how-can-i-calculate-the-length-of-curve#comment_3040941

Answer 2

[Ankit]

X=[96.0741000000000,97.1940000000000,98.3139000000000,99.4338000000000,100.553700000000,101.673600000000,102.793500000000,103.913400000000,105.033300000000,106.153200000000,107.273100000000,108.393000000000,109.512900000000,110.632800000000,111.752700000000,112.872600000000,113.992500000000,115.112400000000,116.232300000000];
Y=[-4.13836296940031,-4.10455468315876,-4.05645426203322,-3.99617782198545,-3.92344322326347,-3.83385191481492,-3.73582865974161,-3.61740402741020,-3.49399064332423,-3.35231953224592,-3.20552503148528,-3.04892626846560,-2.88658570885772,-2.72091440408539,-2.55226630046971,-2.38425597793465,-2.21787687713447,-2.05656258174384,-1.89889800700337];

len_curve = sum(vecnorm(diff( [X(:),Y(:)] ),2,2));
% the 2-norm along the rows of a matrix: vecnorm(A,2,2) , where A is a
% vector 
% diff: Difference and approximate derivative.

[Volcano]

Thanks a lot, but there is small difference between your answer and other one.

[Guntz Romain]

le boss

Answer 3

[Star Strider]

Possibly —

X=[96.0741000000000,97.1940000000000,98.3139000000000,99.4338000000000,100.553700000000,101.673600000000,102.793500000000,103.913400000000,105.033300000000,106.153200000000,107.273100000000,108.393000000000,109.512900000000,110.632800000000,111.752700000000,112.872600000000,113.992500000000,115.112400000000,116.232300000000];
Y=[-4.13836296940031,-4.10455468315876,-4.05645426203322,-3.99617782198545,-3.92344322326347,-3.83385191481492,-3.73582865974161,-3.61740402741020,-3.49399064332423,-3.35231953224592,-3.20552503148528,-3.04892626846560,-2.88658570885772,-2.72091440408539,-2.55226630046971,-2.38425597793465,-2.21787687713447,-2.05656258174384,-1.89889800700337]

Y = 1×19

-4.1384 -4.1046 -4.0565 -3.9962 -3.9234 -3.8339 -3.7358 -3.6174 -3.4940 -3.3523 -3.2055 -3.0489 -2.8866 -2.7209 -2.5523 -2.3843 -2.2179 -2.0566 -1.8989

dX = gradient(X);                                       % Numerical Derivative
dY = gradient(Y);                                       % Numerical Derivative
Len = cumtrapz(X,hypot(dX,dY))                          % Integrate The Hypotenuse Of The Numerical Derivatives Of The Segments

Len = 1×19

0 1.2549 2.5102 3.7662 5.0231 6.2812 7.5405 8.8012 10.0634 11.3271 12.5922 13.8584 15.1256 16.3934 17.6616 18.9298 20.1976 21.4648 22.7315

figure
plot(X, Y, '.-')
hold on
plot(X, Len, '.-')
hold off
grid

.

[Volcano]

Thanks a lot, but there is small difference between your answer and other one.

[Torsten]

And what is "small" ?

[Volcano]

Ankit's and Star's answers are 20.2980 and 22.7315, respectively. What may cause this difference?

[Star Strider]

I use the gradient function to calculate the derivatives., It produces a different (and in my opinion more accurate) estimate of the derivative than diff (that by definition also results in a vector that is one element shorter than the original), while the length of the gradient result is the same as the original.

[Tamas Rozsa]

@Star Strider

I believe cumtrapz() is incorrect here, one should utilize cumsum() instead:

Len = cumsum(hypot(dX,dY))

Reason: No integration is needed here (as dX and dY are not exactly derivatives either, due to the lack of spacing specification). If you consider a vertical line as example, your solution with cumtrapz() will result 0 length, which is obviously wrong. With cumsum(), it will also work, even with arbitrary spacing in the input data.

[Star Strider]

My code calculates the trapezoidal integral of the gradients (numerical derivatives) of ‘X’ and ‘Y’.

[Star Strider]

Following up —

In light of @Torsten’s Comment, using cumtrapz is correct, however including the independent variable as the first argument to it in this instance may not be —

X=[96.0741000000000,97.1940000000000,98.3139000000000,99.4338000000000,100.553700000000,101.673600000000,102.793500000000,103.913400000000,105.033300000000,106.153200000000,107.273100000000,108.393000000000,109.512900000000,110.632800000000,111.752700000000,112.872600000000,113.992500000000,115.112400000000,116.232300000000];
Y=[-4.13836296940031,-4.10455468315876,-4.05645426203322,-3.99617782198545,-3.92344322326347,-3.83385191481492,-3.73582865974161,-3.61740402741020,-3.49399064332423,-3.35231953224592,-3.20552503148528,-3.04892626846560,-2.88658570885772,-2.72091440408539,-2.55226630046971,-2.38425597793465,-2.21787687713447,-2.05656258174384,-1.89889800700337];

dX = gradient(X);                                       % Numerical Derivative
dY = gradient(Y);                                       % Numerical Derivative
Len = cumtrapz(hypot(dX,dY));                         % Integrate The Hypotenuse Of The Numerical Derivatives Of The Segments
Len_end = Len(end)

Len_end = 20.2978

figure
plot(X, Y, '.-', 'DisplayName','Data')
hold on
plot(X, Len, '.-', 'DisplayName','Length')
hold off
grid
title('Provided Data')
legend('Location','best')
axis('equal')
axis('padded')

a = linspace(0, 2*pi, 1E+3);
X = cos(a);
Y = sin(a);

dX = gradient(X);                                       % Numerical Derivative
dY = gradient(Y);                                       % Numerical Derivative
Len = cumtrapz(hypot(dX,dY));                           % Integrate The Hypotenuse Of The Numerical Derivatives Of The Segments
Len_end = Len(end)

Len_end = 6.2831

figure
plot(X, Y, '.-', 'DisplayName','Data')
hold on
plot(X, Len, '.-', 'DisplayName','Length')
hold off
grid
title('Circle (Radius = 1)')
legend('Location','best')
axis('equal')
axis('padded')

My code is unchanged, however it now includes a second data set (the circle) whose result can be checked.

.

Answer 4

[Torsten]

I'd say Ankit's solution is the more intuitive.

But Star Strider's solution should be second-order accurate while Ankit's is only first-order accurate.

X=[96.0741000000000,97.1940000000000,98.3139000000000,99.4338000000000,100.553700000000,101.673600000000,102.793500000000,103.913400000000,105.033300000000,106.153200000000,107.273100000000,108.393000000000,109.512900000000,110.632800000000,111.752700000000,112.872600000000,113.992500000000,115.112400000000,116.232300000000];
Y=[-4.13836296940031,-4.10455468315876,-4.05645426203322,-3.99617782198545,-3.92344322326347,-3.83385191481492,-3.73582865974161,-3.61740402741020,-3.49399064332423,-3.35231953224592,-3.20552503148528,-3.04892626846560,-2.88658570885772,-2.72091440408539,-2.55226630046971,-2.38425597793465,-2.21787687713447,-2.05656258174384,-1.89889800700337];
length = 0;
for i = 1:numel(X)-1
    length = length + sqrt((X(i+1)-X(i))^2 + (Y(i+1)-Y(i))^2);
end
length

length = 20.2980

Answer 5

[Tamas Rozsa]

Based on @Star Strider's answer, but with correct result:

dX = gradient(X);
dY = gradient(Y);
% option 1
Len = cumsum(hypot(dX,dY))  % if sublengths of all segments also needed
% option 2
Len = sum(hypot(dX,dY))     % if only total length needed

As @Star Strider also highlighted in comment, gradient() may be substituted with diff(), but gradient() may give more satisfactory (i.e., smoother) result in most cases. [UPDATE: in some cases, and depending on the actual use-case]

Unlike @Star Strider's original answer, this code gives correct result even in case of arbitrary spacing in the input data as well as in case of vertical line segments.

[Paul]

Why not just compare the methods to a known result? For example, find the arclength of the top half of the unit cirle. Answer should be pi.

format short e

% the data
X = -1:.01:1;
Y = sqrt(1-X.^2);

% Tamas Rozsa solution
dX = gradient(X);
dY = gradient(Y);
Len = sum(hypot(dX,dY));
[Len , abs(Len - pi)]

ans = 1×2

3.2824e+00 1.4077e-01

% Star Strider solution
dX = gradient(X);                                       % Numerical Derivative
dY = gradient(Y);                                       % Numerical Derivative
Len = cumtrapz(X,hypot(dX,dY));                         % Integrate The Hypotenuse Of The Numerical Derivatives Of The Segments
[Len(end) , abs(Len(end) - pi)]

ans = 1×2

3.1409e-02 3.1102e+00

% Ankit solution
len_curve = sum(vecnorm(diff( [X(:),Y(:)] ),2,2));
[len_curve , abs(len_curve - pi)]

ans = 1×2

3.1413e+00 2.9409e-04

% trapezoidal integration of the arclength integrand, function with
% continous derivative
dX = gradient(X);
dY = gradient(Y);              
len = trapz(X,sqrt(1 + (dY./dX).^2));
[len , abs(len - pi)]

ans = 1×2

3.1409e+00 6.4720e-04

[Tamas Rozsa]

@Paul

Because it also depends on the input data.

My main message was not really about the accuracy, but to point out that @Star Strider's solution is buggy and shall never be used by anyone in that way.

It accidentally gave kinda fair result for the OP input data, but fails heavily e.g. in the two scenarios I've mentioned:

• arbitrary spacing, e.g.:

    X = logspace(-10,0,200);
    Y = sqrt(1-X.^2);

(result is 0.11 instead of ca. 1.57 (pi/2))

• vertical line segments, e.g.:

    X = ones(1,200);
    Y = logspace(-10,1,200);

(result is 0 instead of ca. 10).

The accuracy is another topic; diff() is indeed better here than gradient().

[Paul]

I think the example I showed reinforces your concerns.

Answer 6

[Walter Roberson]

No, you cannot really get the length of a curve defined by a finite list of points. A finite list of 2D points does not define a curve: a finite list of 2D points defines a polygon at best (possibly a self-interesecting one.)

In order to get a curve length, you either have to be given a curve equation, or else you have to be willing to approximate the true curve length by using a model. The model might over-estimate or under-estimate the true curve length.

Mathematically it is impossible to be given a finite set of points that are finitely expressed, and use them to come up with "the" defining curve. Mathematically given any finite set of points that are finitely expressed, there are an uncountable infinity of curves that go through the given points. (Allowing for the possibility that there is noise or round-off or truncation error in the list of coordinates does not increase the number of possible curves, as uncountable infinity is the largest infinity already until you get into abstractions such as Aleph-One )

The answers posted by the other participants are either finding total segment lengths (treating the point list as a polygon whose perimeter is to be found), or else are using different models of how to interpolate the points into a curve. They produce different results because they use different interpolation methods. That does not make any of them "wrong", just different. Unless you know the form of the original function, you just have to accept that the problem is under-specified.