Documentation

## 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:

1. How to efficiently track 1000 boids?

2. 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 ```

Get trial now