MATLAB Answers

0

Help optimising code that creates linked particle trajectories from single particle positions in each frame

Asked by Manny Kins on 7 Jul 2019
Latest activity Commented on by Manny Kins on 19 Jul 2019
The code below takes particle positions that I have found in each frame of a video and compares nearest neighbours to build particle trajectories.
It works well for a small number of particles on screen but when I have a very large number (3000 particles) it becomes very slow. Is there any part of this code I can optimise for better performance? I am not that great at vectorisation so maybe I am missing something on that front.
I have attached a file which contains the positions of the particles in the first 10 frames to test with. Currently it takes around 5 seconds per frame. I have 18000 frames for my usual videos so it takes a while!
By far the slowest part is calculating the distances between particles in the current and previous frame. The code is below.
%PositionP is n by 1 struct where n is the total number of frames being analysed.
%The struct contains PositionP.PX and PositionP.PY which are X and Y coordinates for all particles
%in each frame. So for example frame 1 contains 3462 particles, hence
%length(PositionP(1).PX) = 3462.
%Initialise an array containing all particle X and Y positions in first
%frame.
PX = PositionP(1).PX;
PY = PositionP(1).PY;
%Initialise structs for keeping particle time, X positions (PX) and Y
%positions (PY). 'Linked' will be the struct where we keep all of our
%connected (linked) particle trajectories.
Linked.P(1).T = [0];
Linked.P(1).PX = [0];
Linked.P(1).PY = [0];
%for all particles found on first frame, assign time = 0
for j = 1:1:length(PX)
Linked.P(j).T = [0];
Linked.P(j).PX = [PX(j)];
Linked.P(j).PY = [PY(j)];
end
%Below is main loop that takes a lot of time
tic
for i = 2:1:length(PositionP)
%PX and PY now become the X and Y positions of all particles in ith
%frame, starting at frame 2
PX = PositionP(i).PX;
PY = PositionP(i).PY;
%Starting from the first particle, find nearest neighbour in next frame
for j = 1:1:length(Linked.P)
%if a particle did not find a nearest neighbour and hence increase
%its trajectory last frame, then assume it is lost and do not keep
%searching in following frames,
%it does not enter loop below.
if (((max(Linked.P(j).T)) == (i-2)))
%find distance between each particle in current and previous
%frame. Linked.P(j).PX(end) represents the latest position
%of linked particle trajectory. This is the part which takes longest
Distance = sqrt((PX-Linked.P(j).PX(end)).^2 + (PY-Linked.P(j).PY(end)).^2);
%find index of smallest distance between particles in both frames
MinDistanceIndex = find(Distance==min(Distance));
%if more than zero and one min distance, just take first
if (length(MinDistanceIndex)>0)
MinDistanceIndex = MinDistanceIndex(1);
%if particle in current frame is within a certain distance
%(5) from previous frame, then link positions to make a trajectory.
if (Distance(MinDistanceIndex)<5)
Linked.P(j).T = [Linked.P(j).T (i-1)];
Linked.P(j).PX = [Linked.P(j).PX PX(MinDistanceIndex)];
Linked.P(j).PY = [Linked.P(j).PY PY(MinDistanceIndex)];
PX(MinDistanceIndex) = Inf;
PY(MinDistanceIndex) = Inf;
end
end
end
end
%if any particles have not been assigned to a previous trajectory, then
%assume a new trajectory is being created (new particle has entered
%the frame)
for k = 1:1:length(PX)
if (PX(k)<Inf && PY(k)<Inf)
j = j+1;
Linked.P(j).T = [(i-1)];
Linked.P(j).PX = [PX(k)];
Linked.P(j).PY = [PY(k)];
end
end
i
toc
end

  0 Comments

Sign in to comment.

1 Answer

Answer by Srivardhan Gadila on 19 Jul 2019
 Accepted Answer

The following changes improved the runtime
  • Compute the Distance using the hypot function
  • Use the function min for getting the minimum value and it’s index
Below is the final code after the changes
(other methods to compute distance are also mentioned in the comments)
%PositionP is n by 1 struct where n is the total number of frames being analysed.
%The struct contains PositionP.PX and PositionP.PY which are X and Y coordinates for all particles
%in each frame. So for example frame 1 contains 3462 particles, hence
%length(PositionP(1).PX) = 3462.
%Initialise an array containing all particle X and Y positions in first
%frame.
PX = PositionP(1).PX;
PY = PositionP(1).PY;
%Initialise structs for keeping particle time, X positions (PX) and Y
%positions (PY). 'Linked' will be the struct where we keep all of our
%connected (linked) particle trajectories.
Linked.P(1).T = 0;
Linked.P(1).PX = 0;
Linked.P(1).PY = 0;
%for all particles found on first frame, assign time = 0
for j = 1:1:length(PX)
Linked.P(j).T = 0;
Linked.P(j).PX = PX(j);
Linked.P(j).PY = PY(j);
end
%Below is main loop that takes a lot of time
tic
for i = 2:1:length(PositionP)
%PX and PY now become the X and Y positions of all particles in ith
%frame, starting at frame 2
PX = PositionP(i).PX;
PY = PositionP(i).PY;
%Starting from the first particle, find nearest neighbour in next frame
for j = 1:1:length(Linked.P)
%if a particle did not find a nearest neighbour and hence increase
%its trajectory last frame, then assume it is lost and do not keep
%searching in following frames,
%it does not enter loop below.
if (((max(Linked.P(j).T)) == (i-2)))
%find distance between each particle in current and previous
%frame. Linked.P(j).PX(end) represents the latest position
%of linked particle trajectory. This is the part which takes longest
% Distance = sqrt((PX-Linked.P(j).PX(end)).^2 + (PY-Linked.P(j).PY(end)).^2);
Distance = hypot(PX-Linked.P(j).PX(end), PY-Linked.P(j).PY(end));
% Distance = vecnorm([PX PY] - [Linked.P(j).PX(end) Linked.P(j).PY(end)],2,2);
% Distance = pdist2([Linked.P(j).PX(end) Linked.P(j).PY(end)], [PX PY]);
%find index of smallest distance between particles in both frames
% MinDistanceIndex = find(Distance==min(Distance));
[MinDistance, MinDistanceIndex] = min(Distance);
%if more than zero and one min distance, just take first
if (~isempty(MinDistanceIndex))
MinDistanceIndex = MinDistanceIndex(1);
%if particle in current frame is within a certain distance
%(5) from previous frame, then link positions to make a trajectory.
if (Distance(MinDistanceIndex)<5)
% if (MinDistance<5)
Linked.P(j).T = [Linked.P(j).T (i-1)];
Linked.P(j).PX = [Linked.P(j).PX PX(MinDistanceIndex)];
Linked.P(j).PY = [Linked.P(j).PY PY(MinDistanceIndex)];
PX(MinDistanceIndex) = Inf;
PY(MinDistanceIndex) = Inf;
end
end
end
end
%if any particles have not been assigned to a previous trajectory, then
%assume a new trajectory is being created (new particle has entered
%the frame)
for k = 1:1:length(PX)
if (PX(k)<Inf && PY(k)<Inf)
j = j+1;
Linked.P(j).T = (i-1);
Linked.P(j).PX = PX(k);
Linked.P(j).PY = PY(k);
end
end
i
toc
end

  1 Comment

Hi Vardhan, thanks for taking time to answer my question. The hypot fix does indeed give faster results so thank you!
I have also found that removing the second loop:
for j = 1:1:length(Linked.P)
and just making a matrix of all the end values of Linked, then comparing the matrix of end values with the new particles in PX and PY gives a great speed increase with pdist2.

Sign in to comment.