Forward, backward and central differences

Hi, I am a biology student and new to MATLAB. I am working on a project to study the behavior of small aquatic fishes. I am have captured a video of 2000 frames and the time difference between frames is 1/500 sec. I have a software which gives me the positions of the fishes. I have a MAT file which contains the data of the centroid positions of the fishes. The structure array looks like below.

Fish(1).frames=[10, 11,12,13,14]; 
Fish(1).position1=[20.5, 20.9, 21.3, 21.7,22.1];
Fish(1).position2=[4.1, 4.27, 4.53, 4.79,4.88];
Fish(2).frames=[1,2,3,4,5,6,7,8,9];
Fish(2).position1=[10.4,10.6,10.7,10.8,11.1, 11.47,11.82,12.31,12.44];
Fish(2).position2=[4.2, 4.7, 4.8, 4.9,5.2,5.28,5.63,5.89,6.01];

The size of the structure is 30,000 X 1. I want to calculate the velocity of the fishes based on position1 (velocity1) and position2(velocity2). I want to apply forward and backward differences to the end data points and central difference to the inside data. I don't know how to do this. Please, help me.

Answer 1

[Bruno Luong]

dt=1/500;
for k=1:length(Fish)
    Fish(k).velocity1 = gradient(Fish(k).position1,dt);
    Fish(k).velocity2 = gradient(Fish(k).position2,dt);
end

[Tomer]

Thank you.

[Bruno Luong]

In the gradient doc it states:

Algorithms

gradient calculates the central difference for interior data points. For example, consider a matrix with unit-spaced data, A, that has horizontal gradient G = gradient(A). The interior gradient values, G(:,j), are

G(:,j) = 0.5*(A(:,j+1) - A(:,j-1));

The subscript j varies between 2 and N-1, with N = size(A,2).

gradient calculates values along the edges of the matrix with single-sided differences:

G(:,1) = A(:,2) - A(:,1);

G(:,N) = A(:,N) - A(:,N-1);

If you specify the point spacing, then gradient scales the differences appropriately. If you specify two or more outputs, then the function also calculates differences along other dimensions in a similar manner. Unlike the diff function, gradient returns an array with the same number of elements as the input.

[Tomer]

Thanks.

Answer 2

[Adam Danz]

Assuming the camera is looking down on the fish, since you do not have position 3 (vertical position within the water) you can only compute the x and y velocity components. This could result in an underestimation of velocity. Imagine a fish is swiming quickly in an upward and forward direction. As the angle from the ground increases toward 90 degrees, the forward component becomes smaller relative to the vertical component even though the fish is accelerating. Some systems include a 2nd camera to capture the vertical movement as well in order to get around that problem.

Nevertheless, the solution to your task can be broken down into a few manageable parts.

First, given the (x,y) coordinates and assuming no change in z, you'll need to compute the distance moved between coordinates. So if you have n coordinates, you'll have n-1 distances. Hint: distance formula.

Second, you'll need to compute the time between each position datapoint. That should be easy because you have a fixed frame rate and you have frame counts. Again, if you have n timestamps, you'll end up with n-1 durations.

Now you have a vector of distances and a vector of times. Lastly, velocity is just distance divided by time.

That should get you started. Feel free to followup if you get stuck.

[Tomer]

The data is after making the corrections you have mentioned. I want to use the above mentioned schemes to have the same dimension of the velocity field as the positions. I don't know how this can be done in a loop for the entire structure array.

[Adam Danz]

To compute velocity, you need distance. You only have position. Distance is the change in position so if you have n positions, you'll end up with n-1 distances. The only way around that is if you have additional informaiton such as the velocity at time 0 or a starting coordinate that comes before position #1. Another way around that is to interpolate the velocity curve from n-1 to n number of data points.

You could loop through each fish but you certainly do not have to loop through each position. Use the vectors of (x,y) positions to compute a vector of distances in 1 line of code. Same with computing duration from the time stamps.

To compute velocity for 1 fish, you'll need ~3 lines of code.

[Tomer]

Each Fish(i) represents the path of one single fish. Fish(1).position1(1) is the starting x-coordinate of that particular fish. Also, you have the time 1/500 sec for movement between Fish(1).position1(1) to Fish(1).position1(2).

The same is for other Fish(2) so on ,they are data of different fishes.

I am a biology student and doesn't know much about numerical techniques and MATLAB.

Can you help me with writing a loop for this. Forward and backward differences at the end points of the data and central differences to inner data.

[Adam Danz]

Assuming Fish(1).position1 are x coordinates and Fish(1).position2 are y coordinates, this loop will add a new field to each structure, named "velocity". "Velocity" will be a vector of length n-1 where n is the number of position and frame elements.

nFish = numel(Fish); 
for i = 1:nFish
    frameRate = 1/500; % sec
    distance = sqrt((Fish(i).position1(2:end) - Fish(i).position1(1:end-1)).^2 + ...
        (Fish(i).position2(2:end) - Fish(i).position2(1:end-1)).^2);
    time = (Fish(i).frames(2:end) - Fish(i).frames(1:end-1)).* frameRate;
    Fish(i).velocity = distance ./ time;
end

[Tomer]

Thank you. Does this gives the x-velocity and y-velocity? I want to calculate x-velocity based on x-positions and y-velocity based on y-positions.

[Adam Danz]

No, the above gives total velocity. If you want the x and y components of velocity, you'll need to adapt the code above. The distance calculation is simplified to

distance_x = Fish(i).position1(2:end) - Fish(i).position1(1:end-1) ; 
distance_y = Fish(i).position2(2:end) - Fish(i).position2(1:end-1) ; 
% time =  same thing as above
Fish(i).velocity_x = distance_x ./ time;
% same for velocity_y....

Again, velocity is just distance divided by time. That's why velocity is reported in units such as miles per hour (miles is distance, divided by time, hours). If you want x and y velocities, you need x and y distances which are just the different between each x and y position value. Taking some time to understand that will be beneficial to the rest of your analysis.

[Tomer]

Thanks.

[Adam Danz]

No problem.

Feel free to come back and check in if you get stuck.

When you get a chance, you might want to accept the answers that were helpful to you. Here are a list of your questions:

https://www.mathworks.com/matlabcentral/answers/?term=asked_by_id%3A295611