bat365 Homepage

Tracking a Flock of Birds

This example shows how to track a large number of objects. A large flock of birds is generated and a global nearest neighbor multi-object tracker, trackerGNN, is used to estimate the motion of every bird in the flock.

Scenario Definition

The flock motion is simulated using the behavioral model proposed by Reynolds [1]. In this example, the flock is comprised of 1000 simulated birds, called boids, whose initial position and velocity was previously saved. They follow the three rules of flocking: collision avoidance, velocity matching, and flock centering. Each rule is associated with a weight and the overall behavior of the flock emerges from the relative weighting of each rule. In this case, weights are chosen that cause the flock to fly around a certain point and create a dense center. Other weight settings can cause different behaviors to emerge.

Tracking such a large and dense flock presents two challenges:

How to efficiently track 1000 boids?
How to be able to track individual boids in such a dense environment?

The following code simulates the flock behavior for 100 steps of 0.1 second. The plot on the left shows the flock as a whole and the plot on the right is zoomed in on the densest part at the flock center.

s = rng;   % Keep the current state of the random number generator
rng(2019); % Set the random number generator for repeatable results
load("initialFlock.mat","x","v");
flock = helperFlock("NumBoids",size(x,1),"CollisionAviodanceWeight",0.5,...
    "VelocityMatchingWeight",0.1,"FlockCenteringWeight",0.5,"Velocity",v,...
    "Position",x,"BoidAcceleration",1);
truLabels = string(num2str((1:flock.NumBoids)'));
bound = 20;
flockCenter = mean(x,1);
[tp1,tp2] = helperCreateDisplay(x,bound);

% Simulate 100 steps of flocking
numSteps = 100;
allx = repmat(x,[1 1 numSteps]);
dt = 0.1;
for i = 1:numSteps
    [x,v] = move(flock,dt);
    allx(:,:,i) = x;
    plotTrack(tp1.Plotters(1),x)
    inView = findInView(x,-bound+flockCenter,bound+flockCenter);
    plotTrack(tp2.Plotters(1),x(inView,:),truLabels(inView))
    drawnow
end

Tracker Definition

You define the tracker as shown in the example How to Efficiently Track Large Numbers of Objects.

You observe that the boids follow a curved path and choose a constant turn model defined by initctekf.

To limit the time required to calculate cost, you reduce the coarse cost calculation threshold in the AssignmentThreshold to a low value.

Further, you choose the more efficient Jonker-Volgenant as the assignment algorithm, instead of the default Munkres algorithm.

You want tracks to be quickly confirmed and deleted, and set the confirmation and deletion thresholds to [2 3] and [2 2], respectively.

Finally, you know that the sensor scans only a fraction of the flock at any given scan, and so you set the HasDetectableTrackIDsInput to true to be able to pass the detectable track IDs to the tracker.

The following line shows how the tracker is configured with the above properties. You can see how to generate code for a tracker in How to Generate C Code for a Tracker, and the tracker for this example is saved in the function flockTracker_kernel.m

% tracker = trackerGNN("FilterInitializationFcn",@initctekf,"MaxNumTracks",1500,...
%         "AssignmentThreshold",[50 800],"Assignment","Jonker-Volgenant",...
%         "ConfirmationThreshold",[2 3],"DeletionThreshold",[2 2],...
%         "HasDetectableTrackIDsInput",true);

Track the Flock

Next, you run the scenario and track the flock.

A simplified sensor model is simulated using the detectFlock supporting function. It simulates a sensor that scans the flock from left to right, and captures a fifth of the flock span in the x-axis in every scan. The sensor has a 0.98 probability of detection and the noise is simulated using a normal distribution with a standard deviation of 0.1 meters about each position component.

The sensor reports its currentScan bounds, which are used to provide the detectable track IDs to the tracker.

clear flockTracker_kernel
positionSelector = [1 0 0 0 0 0 0; 0 0 1 0 0 0 0; 0 0 0 0 0 1 0];
trackIDs = zeros(0,1,'uint32');
trax = zeros(0,3);
bounds = inf(3,2);
alltrax = zeros(size(allx));
allIDs = repmat({},1,numSteps);
trup2 = tp2.Plotters(1);
trap2 = tp2.Plotters(2);
trup2.HistoryDepth = 2*trap2.HistoryDepth;
clearPlotterData(tp1)
clearPlotterData(tp2)
for i = 1:numSteps
    t = i*dt;
    [detections, currentScan] = detectFlock(allx(:,:,i),t);
    bounds(1,:) = currentScan;
    tracksInScan = findInView(trax,bounds(:,1),bounds(:,2));
    [tracks,info] = flockTracker_kernel(detections,t,trackIDs(tracksInScan,1));
    trax = getTrackPositions(tracks,positionSelector);
    if ~isempty(tracks)
        trackIDs = uint32([tracks.TrackID]');
    else
        trackIDs = zeros(0,1,'uint32');
    end
    alltrax(1:size(trax,1),1:3,i) = trax;
    allIDs{i} = string(trackIDs);
    helperVisualizeDisplay(tp1,tp2,truLabels,allx,allIDs,alltrax,i)
end
rng(s); % Reset the random number generator to its previous state

Result of the Tracker in Generated Code

The following GIF shows the performance of the tracker in a mex file.

Summary

This example showed how to track a large number of objects in a realistic scenario, where a scanning sensor only reports a fraction of the objects in each scan. The example showed how to set the tracker up for large number of objects and how to use the detectable track IDs input to prevent tracks from being deleted.

References

[1] Craig W. Reynolds, "Flocks, Herds, and Schools: A Behavioral Model", Computer Graphics, Vol. 21, Number 4, July 1987.

Supporting Functions

helperCreateDisplay

The function creates the example display and returns a handle to both theater plots.

function [tp1,tp2] = helperCreateDisplay(x,bound)
f  = figure("Visible", "off");
set(f,"Position",[1 1 1425 700]);
movegui(f,"center")
h1 = uipanel(f,"FontSize",12,"Position",[.01 .01 .48 .98],"Title","Flock View");
h2 = uipanel(f,"FontSize",12,"Position",[.51 .01 .48 .98],"Title","Flock Center");
flockCenter = mean(x,1);

a1 = axes(h1,'Position',[0.05 0.05 0.9 0.9]);
grid(a1,'on')
tp1 = theaterPlot("Parent",a1);

% Flock View (Truncated)
halfspan = 250;
tp1.XLimits = 100*round([-halfspan+flockCenter(1) halfspan+flockCenter(1)]/100);
tp1.YLimits = 100*round([-halfspan+flockCenter(2) halfspan+flockCenter(2)]/100);
tp1.ZLimits = 100*round([-halfspan+flockCenter(3) halfspan+flockCenter(3)]/100);
trackPlotter(tp1,"DisplayName","Truth","HistoryDepth",0,"Marker","^","MarkerSize",4,"ConnectHistory","off");
set(findall(a1,"Type","line","Tag","tpTrackHistory_Truth"),"Color","k");
view(a1,3)
legend('Location','NorthEast')

% Flock center
a2 = axes(h2,'Position',[0.05 0.05 0.9 0.9]);
grid(a2,'on')
tp2 = theaterPlot("Parent",a2);
tp2.XLimits = 10*round([-bound+flockCenter(1) bound+flockCenter(1)]/10);
tp2.YLimits = 10*round([-bound+flockCenter(2) bound+flockCenter(2)]/10);
tp2.ZLimits = 10*round([-bound+flockCenter(3) bound+flockCenter(3)]/10);
trackPlotter(tp2,"DisplayName","Truth","HistoryDepth",0,...
    "Marker","^","MarkerSize",6,"ConnectHistory","off","FontSize",1);
set(findall(a2,"Type","line","Tag","tpTrackHistory_Truth"),"Color","k");

% Track plotters
TrackColor = [0 0.4470 0.7410]; % Blue
TrackLength = 50;
trackPlotter(tp1,"DisplayName","Tracks","HistoryDepth",TrackLength,"ConnectHistory","off",...
    "Marker",".","MarkerSize",3,"MarkerEdgeColor",TrackColor,"MarkerFaceColor",TrackColor);
set(findall(tp1.Parent,"Type","line","Tag","tpTrackHistory_Tracks"),...
    "Color",TrackColor,"MarkerSize",3,"MarkerEdgeColor",TrackColor);
trackPlotter(tp2,"DisplayName","Tracks","HistoryDepth",TrackLength,"ConnectHistory","on",...
    "Marker","s","MarkerSize",8,"MarkerEdgeColor",TrackColor,"MarkerFaceColor","none","FontSize",1);
set (findall(tp2.Parent,"Type","line","Tag","tpTrackPositions_Tracks"),"LineWidth",2);
set(findall(tp2.Parent,"Type","line","Tag","tpTrackHistory_Tracks"),"Color",TrackColor,"LineWidth",1);

view(a2,3)
legend('Location','NorthEast')
set(f,'Visible','on')
end

detectFlock

The function simulates the sensor model. It returns an array of detections and the current sensor scan limits.

function [detections,scanLimits] = detectFlock(x,t)
persistent sigma allDetections currentScan numScans
numBoids = size(x,1);
pd = 0.98;
if isempty(sigma)
    sigma = 0.1;
    oneDet = objectDetection(0,[0;0;0],"MeasurementNoise",sigma,'ObjectAttributes',struct);
    allDetections = repmat(oneDet,numBoids,1);
    currentScan = 1;
    numScans = 5;
end

% Vectorized calculation of all the detections
x = x + sigma*randn(size(x));
[allDetections.Time] = deal(t);
y = mat2cell(x',3,ones(1,size(x,1)));
[allDetections.Measurement] = deal(y{:});

% Limit the coverage area based on the number of scans
flockXSpan = [min(x(:,1),[],1),max(x(:,1),[],1)];
spanPerScan = (flockXSpan(2)-flockXSpan(1))/numScans;
scanLimits = flockXSpan(1) + spanPerScan * [(currentScan-1) currentScan];
inds = and(x(:,1)>=scanLimits(1), x(:,1)<=scanLimits(2));

% Add Pd
draw = rand(size(inds));
inds = inds & (draw<pd);
dets = allDetections(inds);
detections = num2cell(dets);

% Promote the scan count
currentScan = currentScan+1;
if currentScan > numScans
    currentScan = 1;
end
end

findInView

The function returns a logical array for positions that fall within the limits of minBound and maxBound.

function inView = findInView(x,minBound,maxBound)
inView = false(size(x,1),1);
inView(:) = (x(:,1)>minBound(1) & x(:,1)<maxBound(1)) & ...
         (x(:,2)>minBound(2) & x(:,2)<maxBound(2)) & ...
         (x(:,3)>minBound(3) & x(:,3)<maxBound(3));
end

helperVisualizeDisplay

The function displays the flock and tracks after tracking.

function helperVisualizeDisplay(tp1,tp2,truLabels,allx,allIDs,alltrax,i)
trup1 = tp1.Plotters(1);
trap1 = tp1.Plotters(2);
trup2 = tp2.Plotters(1);
trap2 = tp2.Plotters(2);
plotTrack(trup1,allx(:,:,i))
n = numel(allIDs{i});
plotTrack(trap1,alltrax(1:n,:,i))

bounds = [tp2.XLimits;tp2.YLimits;tp2.ZLimits];

inView = findInView(allx(:,:,i),bounds(:,1),bounds(:,2));
plotTrack(trup2,allx(inView,:,i),truLabels(inView))
inView = findInView(alltrax(1:n,:,i),bounds(:,1),bounds(:,2));
plotTrack(trap2,alltrax(inView,:,i),allIDs{i}(inView))
drawnow
end