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
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 = ;
Linked.P(1).PX = ;
Linked.P(1).PY = ;
%for all particles found on first frame, assign time = 0
for j = 1:1:length(PX)
Linked.P(j).T = ;
Linked.P(j).PX = [PX(j)];
Linked.P(j).PY = [PY(j)];
%Below is main loop that takes a lot of time
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
MinDistanceIndex = MinDistanceIndex(1);
%if particle in current frame is within a certain distance
%(5) from previous frame, then link positions to make a trajectory.
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;
%if any particles have not been assigned to a previous trajectory, then
%assume a new trajectory is being created (new particle has entered
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)];