## 3 Materials and Methods
### BeePIV
#### 3.3.6. PIV and Directional Bee Traffic  
Let \((W_{t},\ldots,W_{t+k})\) be a sequence of motion frames and let \(W_{l}\) and \(W_{l+1}\) be two consecutive motion frames in this sequence that correspond to original video frames \(F_{l}\) and \(F_{l+1}\). Let \(IA_{1}\) be a \(D\times D\) window, referred to as interrogation area or interrogation window in the PIV literature, selected from \(W_{l}\) and centered at position \((i,j)\). Another \(D^{\prime}\times D^{\prime}\) window, \(IA_{2}\), is selected in \(W_{l+1}\) so that \(D\leq D^{\prime}\). The position of \(IA_{2}\) in \(W_{l+1}\) is the function of the position of \(IA_{1}\) in \(W_{l}\) in that it changes relative to \(IA_{1}\) to find the maximum correlation peak. For each possible position of \(IA_{1}\) in \(W_{l}\), a corresponding position \(IA_{2}\) is computed in \(W_{l+1}\). \[C(r,s)=\sum\limits_{i=0}^{D-1}\sum\limits_{j=0}^{D-1}IA_{1}(i,j)IA_{2}(i+r,j+s) \tag{25}\]  
The 2D matrix correlation is computed between \(IA_{1}\) and \(IA_{2}\) with the formula in Equation (25), where \(r,s\) are integers in the interval \([-\lfloor(D+D^{\prime}-1)/2\rfloor,\ldots,\lfloor(D+D^{\prime}-1)/2\rfloor]\). In Equation (25), \(IA_{1}(i,j)\) and \(IA_{2}(i+r,j+s)\) are the pixel intensities at locations \((i,j)\) in \(W_{l}\) and \((i+r,j+s)\) in \(W_{l+1}\). For each possible position \((r,s)\) of \(IA_{1}\) inside \(IA_{2}\), the correlation value \(C(r,s)\) is computed. If the size of \(IA_{1}\) is \(M\times N\) and the size of \(IA_{2}\) is \(P\times Q\), then the size of the matrix \(C\) is \((M+P-1)\times(N+Q-1)\). The matrix \(C\) records the correlation coefficient for each possible alignment of \(IA_{1}\) with \(IA_{2}\). A faster way to calculate correlation coefficients between two image frames is to use the Fast Fourier Transform (FFT) and its inverse, as shown in Equation (26). The reason why the computation of Equation (26) is faster than the computation of Equation (25) is that \(IA_{1}\) and \(IA_{2}\) must be of the same size. \[C(r,s)=\Re[FFT^{-1}(FFT^{*}(IA_{1})\cdot FFT(IA_{2}))] \tag{26}\]  
If \(C(r_{m},s_{m})\) is the maximum value in \(C\) and \((i_{c},j_{c})\) is the center of \(IA_{1}\), the pair of \((i_{c},j_{c})\) and \((r_{m},s_{m})\) defines a displacement vector \(\vec{v}_{i_{c},j_{c},r_{m},s_{m}}\) from \((i_{c},j_{c})\) in \(W_{l}\) to \((i_{c}+r_{m},j_{c}+s_{m})\) in \(W_{l+1}\). This vector represents how particles may have moved from \(W_{l}\) to \(W_{l+1}\). The displacement vectors form a vector field used to estimate possible flow patterns. Figure 12: Motion points from motion frames \(W_{1}\) to \(W_{12}\) in Figure 11 plotted on original frames \(F_{1}\) to \(F_{12}\) in Figure 2 from which motion frames were obtained. Figure 13 shows two frames \(F_{3}\) and \(F_{4}\), two corresponding motion frames \(W_{3}\) and \(W_{4}\) obtained from them by the application of the operator in Equation (24), and the vector field \(VF_{4}\) computed from \(W_{3}\) and \(W_{4}\) by Equation (26). The two displacement vectors correspond to the motions of two bees with the left bee moving slightly left and the right bee moving down and right. In Equation (27), we define the PIV operator that applies to two consecutive motion frames \(W_{l}\) and \(W_{l+1}\) to generate a field of displacement vectors \(VF_{l+1}\). The video PIV operator \(G^{*}(V)\) applies the PIV operator \(G\) to every pair of consecutive motion frames \(W_{l}\) and \(W_{l+1}\) in a sequence of motion frames and returns a sequence of vector fields \(VF_{l+1}\), as defined in Equation (28), where \(Z=(W_{t},\ldots,W_{t+k})\) is a sequence of eroded frames. It should be noted that the number of motion frames in \(Z\) exceeds the number of the vector fields returned by the video PIV operator by exactly \(1\).