What is the best way of calculating the path length of a freely moving tracked animal.

Dear Matlab Community,

I am analysing the ambulation of rats in a rectangular arena. I am using my custom tracking software, delivering pixel coordinates by manually tracking the animal due to cursor. Pixel coordinates are delivered. For linear displacements the 'norm' function gives pretty good approximation. I am wondering what would be the best if the displacement is not just linear. I would use the pdist, but I am uncertain about the distance choice?

Eucledian.......etc?

SUggestions are appreciated!

lg

[Stephen23]

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

Answer 1

[Star Strider]

What you probably want is a line integral. This integrates the hypotenuse of simultaneously varying elements along a level path.

A simple implementation of it, and an example with a circle is —

t = linspace(0, 1, 500);                                    % Time Vector
x = cos(2*pi*t);                                            % Path X-Coordinate
y = sin(2*pi*t);                                            % Path Y-Coordinate

dxdt = gradient(x, t);                                      % Numerical Derivative
dydt = gradient(y, t);                                      % Numerical Derivative

Dist = trapz(t, sqrt(dydt.^2 + dxdt.^2))                    % Calculate Distance Along The Path

Dist = 6.2830

Cdist = cumtrapz(t, sqrt(dydt.^2 + dxdt.^2));               % Calculate Cumulative Distance Along The Path

figure
plot(t, Cdist)
grid
xlabel('Time')
ylabel('Cumulative Distance')

This is the circumference of a unit circle, where the actual value is . You can also use cumtrapz for this, to get the cumulative distance as well as the total distance.

EDIT — (20 Mar 2024 at 17:36)

Adding a hypothetical animal —

t = linspace(0, 1, 500).';                                  % Time Vector
x = cumsum(randn(size(t)));                                 % Path X-Coordinate
y = cumsum(randn(size(t)));                                 % Path Y-Coordinate

x = x - min(x);
y = y - min(y);

dxdt = gradient(x, t);                                      % Numerical Derivative
dydt = gradient(y, t);                                      % Numerical Derivative

Dist = trapz(t, sqrt(dydt.^2 + dxdt.^2))                    % Calculate Distance Along The Path

Dist = 471.7434

Cdist = cumtrapz(t, sqrt(dydt.^2 + dxdt.^2));               % Calculate Cumulative Distance Along The Path

figure
plot(x, y, 'DisplayName','Path')
hold on
plot(x(1), y(1), 'sg', 'MarkerSize',10, 'MarkerFaceColor','g', 'DisplayName','Start')
plot(x(end), y(end), 'sr', 'MarkerSize',10, 'MarkerFaceColor','r', 'DisplayName','Stop')
hold off
grid
xlabel('X Coordinate')
ylabel('Y Coordinate')
title('Random Animal')
axis('equal')
legend('Location','best')

figure
plot(t, Cdist)
grid
xlabel('Time')
ylabel('Cumulative Distance')

I offset the coordinates to be purely positive. That does not make any real difference in the calculations, although it makes it more realistic with respect to an actual animal path.

EDIT — (20 Mar 2024 at 20:30)

Minor correction, and corrected typographical errors.

EDIT — (21 Mar 2024 at 10:46)

A few words on computational errors —

With respect to computational errors in this approach, it is possible to compute the relative error for the circle estimate, since we know the correct value —

t = linspace(0, 1, 500);                                    % Time Vector
x = cos(2*pi*t);                                            % Path X-Coordinate
y = sin(2*pi*t);                                            % Path Y-Coordinate

dxdt = gradient(x, t);                                      % Numerical Derivative
dydt = gradient(y, t);                                      % Numerical Derivative

format longE

CircDist = trapz(t, sqrt(dydt.^2 + dxdt.^2))                % Calculate Distance Along The Path

CircDist = 
     6.283019527771874e+00

C = 2*pi

C = 
     6.283185307179586e+00

RelErr = abs(C-CircDist)/C

RelErr = 
     2.638461220030107e-05

PctErr = RelErr * 100

PctErr = 
     2.638461220030107e-03

If this holds true for the animal path, then an estimate of that error would be (in this example) —

AnimalPathError = Dist * RelErr

AnimalPathError = 
     1.244676713605681e-02

This would be the path error, that is the estimated difference between the actual path length and the true path length, in terms of the path length units.

EDIT — (24 Mar 2024 at 00:48)

One further check on the accuracy of the distance measure is how fast the animal moves. If the measure ‘makes sense’ in that context (with time recorded as well as distance while the animal is being actively observed, and the animal’s average walking speed using the line integral method is consistent with that), then it is most likely accurate. If it significantly overstates or understates the animal’s usual walking speed, then is most likely not.

If you use cumtrapz to get the interim values for the integral (the last value is the total integral), you can estimate the animal’s instantaneous walking speed as well, to get the maximum and minimum walking speed values over the complete path.

.

[Paul]

Is trapz on the gradient the best thing to do when the underlying data that describes the curve are not smooth?

rng(100)
t = linspace(0, 1, 500).';                                  % Time Vector
x = cumsum(randn(size(t)));                                 % Path X-Coordinate
y = cumsum(randn(size(t)));                                 % Path Y-Coordinate
x = x - min(x);
y = y - min(y);
dxdt = gradient(x, t);                                      % Numerical Derivative
dydt = gradient(y, t);                                      % Numerical Derivative
Dist = trapz(t, sqrt(dydt.^2 + dxdt.^2))

Dist = 444.7534

Sum of distances between successive points is quite different, and I think would be the minimum value of the arclength for any assumption of what the assumed curve looks like between the points.

sum(vecnorm(diff([x y],1),2,2))

ans = 626.2842

For comparison, here are the results from @John D'Errico's arclength function that was linked by @Stephen23 in this comment.

arclength(x,y)

ans = 626.2842

arclength(x,y,'spline')

ans = 671.9630

arclength(x,y,'pchip')

ans = 634.3825

[Paul]

The arclength integral

applies for continously differentiable functions. Approximating the integrand with gradient for a parametric function that is not continously differentiable needs to be carefully considered. Here's an extreme example.

t = linspace(0, 1, 500).'; 
x = t;
y = 0*t;
y(2:2:end) = 1;
dxdt = gradient(x, t);                                      % Numerical Derivative
dydt = gradient(y, t);                                      % Numerical Derivative
Dist = trapz(t, sqrt(dydt.^2 + dxdt.^2))

Dist = 1.9980

sum(vecnorm(diff([x y],1),2,2))

ans = 499.0010