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

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

Answer 1

[Srivardhan Gadila]

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

[Manny Kins]

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.