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
