A Novel Method for Tracking Individuals of Fruit Fly Swarms Flying in a Laboratory Flight Arena

The growing interest in studying social behaviours of swarming fruit flies, Drosophila melanogaster, has heightened the need for developing tools that provide quantitative motion data. To achieve such a goal, multi-camera three-dimensional tracking technology is the key experimental gateway. We have developed a novel tracking system for tracking hundreds of fruit flies flying in a confined cubic flight arena. In addition to the proposed tracking algorithm, this work offers additional contributions in three aspects: body detection, orientation estimation, and data validation. To demonstrate the opportunities that the proposed system offers for generating high-throughput quantitative motion data, we conducted experiments on five experimental configurations. We also performed quantitative analysis on the kinematics and the spatial structure and the motion patterns of fruit fly swarms. We found that there exists an asymptotic distance between fruit flies in swarms as the population density increases. Further, we discovered the evidence for repulsive response when the distance between fruit flies approached the asymptotic distance. Overall, the proposed tracking system presents a powerful method for studying flight behaviours of fruit flies in a three-dimensional environment.


Introduction
Social behaviours of swarming fruit flies, Drosophila melanogaster, have recently attracted increasing scientific attention [1].Such studies may deepen our understanding of the neural mechanisms underlying social behaviours [2].The rising of such research activities is fueled by recent advances of the video tracking technology.The video tracking method can automatically measure individual's motion states using videos from different camera views [3].Previous studies have made contributions on developing methods for tracking fruit flies walking in planar arenas (2D) [4][5][6][7][8][9] or for tracking fruit flies flying in cubic or cylindrical arenas (3D) [10][11][12][13][14]. Collective behaviour [15], however, usually happens in a large group of animals, such as bird flocks [16,17], fish shoals [18][19][20], and insect swarms [21][22][23][24][25]. Acquiring individual's 3D motion state through time is still an open problem while the population size is large or the population density is high.
In order to obtain 3D quantitative motion data of swarming fruit flies, multiple synchronized and calibrated cameras are employed to capture videos.It needs methods to establish cross-view and cross-frame correspondence for methods to achieve the tracking purpose.Due to severe occlusion and mutually similar appearance, finding correspondences across multiple views is a severe challenge.And moreover, since we have to film the entire swarm in the cameras' field-of-view (FOV), each fruit fly only takes up a few pixels in the images.Such images makes resolving pixels of body and wings more challenging.From this viewpoint, previous methods [14,23,24,26,27] may have the following limitations: (I) The detected result is unstable.Due to the fast wing strokes of fruit flies (less than 4 ms per wing stroke [28]), detected objects which include pixels of the wings and body decrease the accuracy on computing their barycenters; (II) The fruit fly's orientation cannot be reconstructed.Since the pixels of body and wings cannot be resolved, there are no sufficient visual cues for estimating a fly's orientation in 3D.
We propose in this paper a multi-camera tracking system for acquiring motion data (the location and orientation) of individuals of swarms of fruit flies through time.Videos of fruit flies were captured by three synchronized and calibrated high-speed cameras, in which fruit flies were flying in a confined cubic flight arena.We propose an approach to detect pixels corresponding to a fruit fly's body and an algorithm for computing the fly's orientation.The orientation information not only helps increasing the efficiency of tracking fruit flies but also helps validating and correcting the estimated motion state.The tracking system has successfully acquired individual motion data of hundreds of fruit flies through time.We show the results of experiments on five experimental configurations to demonstrate the opportunities that our system offers for generating high-throughput quantitative motion data.We performed quantitative analysis on the kinematics and the spatial structure and the motion patterns of fruit fly swarms.We found that there exists an asymptotic distance between fruit flies in swarms as the population density increases, and we discovered that fruit flies have a strong tendency to turn away when the distance between them approached the asymptotic distance.Overall, the proposed tracking system presents a powerful method for studying flight behaviours of fruit flies in a 3D environment.

Methods
We have developed an experimental imaging system consisting of a transparent cubic flight arena, two planar lights, and three high-speed cameras (Section 2.1).The three cameras are employed to capture videos of flying fruit flies.Videos are the raw observation for the tracking system.Fig 1 shows the overview of the proposed tracking method.The first step of the system is to extract high level observations (measurements) from the raw observation (Section 2.2).Measurements are input to the proposed tracking algorithm.At the tracking step (Section 2.3), the system exploits an "one-to-one" strategy, meaning that it creates (or invokes) one tracker to track one target using measurements associated with the target.After the motion state of a target has been estimated, the system invokes a validating and correcting process (Section 2.4) to achieve results that are more accurate.

Experimental imaging system
The flight arena is a Lucite cube of side length 360 mm (inner).The arena was built by gluing by five transparent acrylic planks with chloroform.

Detection
Previous studies [10][11][12][13][14]23] have demonstrated that the illumination was usually provided by front-lighting in laboratory.The front-lighting means lights and cameras are placed at the same side of the targets.It has the advantage of showing the rich texture of targets and the background.The rich texture however makes target detection and resolving body's pixels very difficult.The back-lighting is the opposite, i.e. the targets are imaged as silhouettes against plain white background.In our problem, we moved the lights to the side of targets opposite the cameras, meaning that targets were back-lit in camera views.Fig 2a shows an image in which fruit flies were back-lit in the camera view.We presented a simple and efficient detecting approach in this paper.The proposed approach makes full use of the different levels of transparency between the wings and body of a fruit fly, in which the pixels of a fly's wings have higher intensities than pixels of its body.The proposed detecting approach segments an image into blobs.Each blob is an image area with pixels' locations and intensities, and each blob represents the observation of a certain fruit fly in a camera view.To simplify notation, all definitions in this section ignore subscripts for cameras and moments.
2.2.1 Segmenting blobs.In the first step, we subtract a Gaussian background model from an acquired image.Given a pixel i which denotes the location of a certain pixel in a certain camera view, let μ(i) be its mean intensity and σ(i) be its standard deviation of intensities through time.We compute μ(i) and σ(i) according to a sequence of consecutively acquired images in the camera view.That is, μ and σ define the Gaussian background model in the camera view.
The pixel i in a certain acquired image is segmented as a foreground pixel if its intensity satisfies the constraint: where I denotes the acquired image.Connected foreground pixels are clustered as blobs b 2 {1, . .., b n }.We choose the constant value C 1 = 1.5 which guarantees that each blob includes pixels of a fly's wings and body.Fig 2b shows the segmented blobs of the image patch.2.2.2 Removing wings pixels.Fruit flies can make wing-strokes very fast (e.g. less than 4 ms [28] for one wing-stroke), meaning the wings' positions between consecutive frames are inconsistent unless the camera's frame-rate is very high, such as 2,000 frames per second.That is unusual for conventional multi-camera system.From this viewpoint, methods, such as [14,27], which included the wings pixels in the detecting results were not precise enough.
In the second step, we use local Gaussian models to remove pixels of wings from blobs.Given a blob b, let μ b denote its mean intensity and σ b denote its standard deviation of the intensities.These variables, μ b and σ b , are computed using all pixels belonging to the blob b.The local Gaussian model of the blob b is thereby defined by μ b and σ b .This second step is a refinement step and is computed as where b 2 {1, . .., b n } denotes the blobs and i 2 {i b,1 , . .., i b,n } denotes the pixels of a certain blob b.If blob(b, i) is equal to 0, we remove the pixel i from the blob b.We choose the constant value C 2 = 1.5 which guarantees most of pixels of a fly's body being preserved after refinement.Fig 2c shows the image patch overlaid with the refined blobs.This step makes full use of the different levels of transparency between a fruit fly's wings and its body (i.e. the pixels of a fly's wings have higher intensities than those of its body).This is achieved thanks to the fruit flies were back-lit in camera views.2.2.3 Fitting blob with ellipse.In order to fit a blob with an ellipse accurately, it is usually that fitting a 2D Gaussian to the locations of all pixels in the blob.Given the parameters of the best-fitting Gaussian, the parameters of the ellipse can be computed.Instead of just computing the mean and covariance of all pixels in the blob, we compute a weighted mean and covariance.Let w(i) denote the weight of pixel i, the weight w(i) is defined as the normalized distance of the pixel i's intensity to the background model which is defined by μ and σ: Therefore, the weighted mean and covariance of the pixels of the blob b are computed as: where W = ∑ i w(i), i 2 {i b,1 , . .., i b, n } defines the weight's normalization constant, and i 2 {i b,1 , . .., i b,n } denotes the pixels belonging to the blob b.The weighted mean and covariance not only improve the robustness to non uniform illumination (e.g. a single threshold for all images), but also give us the sub-pixel accuracy on fitting ellipses.The mean m 0 b defines the center of the ellipse, and the covariance S 0 b defines the axis and the direction of the ellipse.Fig 2d shows the results of fitting each blob with an ellipse.

Measurements definition.
We refer to a blob and the ellipse fitted to it as a measurement.Let χ denote a measurement, each measurement includes two components: where b denotes the blob with pixel locations and intensities, and e(b) defines the ellipse fitting to the blob.We used 3 cameras in our experiments and therefore the set of all available measurements at a certain moment was defined as {χ v,i j v 2 {1, 2, 3}, i 2 {1, . .., i v,n }}.

Tracking
Previous studies [10][11][12][13][14]23] have proposed methods that model targets as points and estimating only targets' locations in 3D.In order to encompass orientation information, we model a target as a directed line-segment in 3D (the center-axis of a target).Fig 3 shows the modeling approach.The fruit fly in our system is denoted by fx; y; z; y; ; lg where l is a constant scalar.
The constant scalar l denotes the length of center axis (In our problem, l ¼ 2:73 mm is the average body length of fruit flies).Here (x, y, z) denotes the location of the fruit fly and (θ, ϕ) denotes its orientation.The proposed tracking algorithm is based on the Bayesian inference framework.Fig 4 shows the overview of the tracking algorithm.In our tracking algorithm, the motion state of targets propagates (location (x, y, z) is predicted and orientation (θ, ϕ) is computed explicitly) in 3D space through time.The location of a target (x, y, z) at moment t is predicted by the dynamic model conditioned on its location at moment t−1.We developed a measurements association algorithm which associates 2D measurements in all views with a 3D target after its location being predicted at each moment.Measurements which have been associated with a certain target generate a matched measurement pair (MMP).The orientation of the target (θ, ϕ) is then computed using the MMP.The predicted location (x, y, z) and computed orientation (θ, ϕ) form a sampled motion state of the target (a particle).The probability of each particle is evaluated by incorporating in the MMP of the target at moment t.In brief, the proposed tracking algorithm is a particle filtering solution for a Bayesian inference problem.
2.3.1 Location prediction.We adopt the Bayesian tracking framework and construct the following eight-dimensional model to define the system state fx; y; z; x À ; y À ; z À ; y; g ð 7Þ where the " − " superscripted terms denote the location of the target at previous moment.
In order to reduce the state space, we partition the state space into two parts: the location state X(x, y, z, x − , y − , z − ) and the orientation state O(θ, ϕ).We thereby predict a target's location state using the dynamic model and the target's location state at moment t−1, and leave the orientation state being computed explicitly.
We adopt the first-order linear extrapolation (FLE) as the dynamic model.The FLE model assumes that the next location of a target is defined by the linear extrapolation of the last two locations.The transition function is defined as: where I 3 is a 3 × 3 identity matrix and v t * N(0, S) is the transition noise.
Let Z 1:t = {Z i , i = 1, . .., t} denote the set of all available observation up to moment t, the Bayesian inference can be formulated as a problem of estimating the posterior probability p (X t jZ 1:t ) [29].Under the first-order Markov assumption and the Bayes' rule, we can get the well-known equation of Bayesian filtering where p(Z t jX t ) is called the observation model.In our problem, Z t ¼ fZ v t g 3 v¼1 denotes the collection of measurements from all 3 camera views at moment t.
By using the particle filtering solution for the Bayesian inference problem, the posterior of a target's motion state at moment t, p(X t jZ 1:t ), is approximated by a set of weighted particles: fðs n t ; w n t Þjn 2 1:::ñ p g.We set ñp ¼ 200 in our experiments.Each particle s n t is associated with a weight w n t which is proportional to the likelihood of the sampled motion state given the observation at moment t, in which the sampled motion state is represented by the particle.New particles are drawn from particles in the previous step using importance sampling [29] and move independently according to the dynamic model.Given a set of particles fðs n t ; w n t Þjn 2 1:::ñ p g, the estimated motion state of the target at moment t, X t , is computed as the expectation In our problem, each particle is defined as s n t ¼ ðX n t ; O n t ; lÞ.At this point, we have the incomplete particle s n t ¼ ðX n t ; Ã; lÞ.The Ã denotes the uninitialized value.2.3.2Measurement association.At each moment, many measurements are extracted from images.We have to decide which measurement is the observation of a certain target in a certain camera view.That is, a measurement has to be assigned to a target if the target is "observed" in that camera view.This is the task of measurement association, and we propose the pixels occupancy test (POT) algorithm for solving the task.The POT algorithm simultaneously solves the association problem across all camera views.The POT algorithm is based on the idea of probability gating.A gate is a surface with constant probability [30,31].In our work it should be an ellipsoid according to the shape of a fruit fly.However, here we employ a sphere gate since there is no orientation state at this point.Let Γ define the association between the particle s n t ðX n t ; Ã; lÞ and measurements χ v,i , v = 1..3, i = 1..i v,n .It is computed as where card(Á) defines the function which counts the number of elements in a set of pixels, and χ v,i (b) denotes the blob's pixels of the measurement χ v,i .Here U 3 denotes the discrete 3D points sampled from the surface of a sphere which locates at X n t ðx; y; zÞ with the diameter equal to l; and P v is the projection matrix of camera v. Therefore, OðX n t ; w v;i Þ defines the common pixels between χ v,i (b) and pixels which are projected into camera v from U 3 ; and η is proportional to the rate of common pixels, η 2 [0, 1].
In Eq (11), η c is a scalar positive-valued percolation threshold which varies from 0 to 1 for leveraging the restriction on association.Larger η c means more restrictive association.In practice, we prefer η c = 0.5 in the first tens of frames.And then η c is set to a smaller value η c = 0.25.Particles, which cannot associate with one measurement in any views, are labeled as "missing particle" and thereby ignored, i.e. the association had to succeed across all camera views.Besides, since the perspective effect causes many 2D occlusions in dense population it is acceptable that one measurement is associated with two or more particles.Fig 5a shows the orientation O n t of the particle can be computed from the direction of the cross line between two planes F 1 and F 2 .However, the direction of the cross line may become problematic while the measurement is close to circular.If the measurements of an MMP are close to circular, the major axis of ellipse might have nothing to do with the plane in which the center-axis lies (as shown in Fig 5b).Let γ define the ratio between the major-axis to the minoraxis of an ellipse (i.e.γ measures the level of circle for a measurement), and let α define the angle between the computed orientation and the real one.We compute the error as = 1−cos (α).Fig 5c shows the error as a function of the ratio γ.It suggests the error is less than 0.01 while γ satisfies γ !1.3.In our case, we can choose two ellipses from an MMP (it includes 3 measurements from 3 camera views), which all satisfy γ !1.3.That is, the orientation O n t can be computed accurately.At this point, we have the particle s n t been completed: s n t ¼ ðX n t ; O n t ; lÞ.2.3.4Weight computation.The particle s n t ¼ ðX n t ; O n t ; lÞ represents a sampled motion state of a certain target at moment t.As aforementioned, after the location state X n t having been predicted, the orientation state O t is computed independently at moment t.That is, w n t , the weight of the particle can be computed as: where Z t denotes the observation at moment t and is defined by GðX n t Þ.The weight includes two factors: p ol , the orientation likelihood; and p al , the appearance likelihood.
Orientation Likelihood p ol defines the likelihood which is computed using the orientation's temporal consistency between consecutive moments.Let the target's orientation at moment t −1 be O t−1 (θ, ϕ), this likelihood is computed as where dðÁÞ defines the direction vector of a target's orientation.
Appearance likelihood p al defines the likelihood which is computed using the appearance's temporal consistency between consecutive moments.The target's appearance is the blobs of the MMP which is associated with the target.The intensities of all pixels in the blob jointly encode features of a target's intrinsic appearance, depth, background and illumination.These factors are approximately constant for consecutive frames while the target were filmed by high- and Camera 2(blue) respectively.The line-segment hx 1 1 ; x 1 2 i(red) denotes the ellipse's major axis in I 1 .And the line-segment hx 2 1 ; x 2 2 i(blue) denotes the ellipse's major axis in I 2 .The projection rays of those end-points define two plane Φ 1 and Φ 2 respectively, e.g. two red rays define Φ 1 .The orientation O n t can be computed from the direction of the cross line (the longer green line-segment) between two planes Φ 1 and Φ 2 .All definitions are defined in the world coordinate system.(b) The problematic orientation computation.The measurement is drawn in red, and red dashed lines are the major-axis and minor-axis of its ellipse.A generative shape according to the particle state ðX speed cameras.The likelihood is thereby computed as where ncc(Á) is the function to compute normalized cross correlation coefficient of the blobs of two measurements.Here, function appearance(v) denotes the function to get the appearance of the target in camera view v.The target's appearance is updated and stored at moment t−1.

Validation and correction
At the end of the tracking process, the motion state ðX t ; O t ; lÞ of a certain target is computed according to Eq (10) at moment t.In order to validate the motion state of a target automatically, we employ the POT algorithm (c.f.Section 2.3.2) again.The POT algorithm associates measurements with a certain target across camera views, but here the discrete 3D points U 3 in Eq (11) are now sampled from the surface of a generated shape.We implemented a parameterized generative shape U In order to invoke restrictive validation, we choose η c = 0.8 (c.f.Eq (11)).The validation is successful if the target has been associated with at least one measurement in each camera view; otherwise failures.If the target is associated with more than one measurement in a certain camera view, we choose the most likely measurement.The successful result forms only one MMP.In the second part of this step, we reconstruct a candidate 3D location by triangulation using the ellipse's center of each measurement in the MMP.If the distance between the target 3D location and the candidate 3D location is larger than l, we update the 3D location of the target using the average of these two locations.The last part of this step is to update and to store the appearance of a target.The appearance of a target is the MMP with a time stamp.The tracking system updates and stores the appearance of an active target in run-time memory.

Experiments and Results
The cameras were calibrated using the method proposed in [33,34] and synchronized by a hardware timing mechanism.All videos were captured at 100 frames per second.The captured videos were first stored in DVR Express 1 and then exported as jpeg format images (IO Industries Canada, DVR Express 1 Core Camera Link Full, monochrome, 10 × 8 bit, Full, 1TB × 3).The capturing and exporting processes were controlled by the Core View™ online console software.All software was running on a computer with Intel 1 Core™ i5-2500 CPU 3.3GHz and 8GB RAM.The proposed tracking method was implemented in the MATLAB™ environment with "Machine Vision Toolbox" [35].The evaluation of tracking performance is provided in S1 File and the raw code is provided in S2 File.The number of flies in each feeding tube was estimated at 120-150 individuals at a normal hatchability ratio after a three-day oviposition by the fresh postnatal parental generation.

Imaging.
The fruit flies were housed in a transparent cubic flight arena, and videos were captured by three back-lit cameras.We have collected videos according to five experiment configurations.For each session, the fruit flies were totally refreshed.At the beginning of each session, we added tubes (c.f.The fruit flies are originally in the feeding tube.) of flies to the arena through the open-top sunroof (see S1 Fig) .The course of flight usually lasted for tens of seconds.We filmed the entire course of flight on each session and discarded the first 2 seconds (200 frames) of each video to compensate for bias on flight behaviours.Table 1 shows the list of videos.The lengths of each video varied, but all were longer than 30 seconds (i.e.3000 frames).

The raw data
During experiments, we found that many fruit flies flew intermittently in the arena.The reason, we believe, is probably the absence of continuous visual stimuli or olfaction stimuli.Our system does not attempt to maintain the identity consistency on fruit flies over extended durations.In this study, if a fruit fly disappeared, for example it landed on a wall for a while, it was labeled with a different identity by the tracking system when the fly appeared again.As a result, the consecutive motion states of a fruit fly (with a unique identity) mainly continued for hundreds of frames, although we have captured rather longer video sequences (thousands of frames).The sequence of consecutive motion states are usually called trajectories.Table 1 shows the number of trajectories of each configuration.We prefer to call these results, which known as trajectories, the raw data (see S3 Fig,which shows the snapshot of the raw data of a typical configuration, T03).

Velocity statistics
Since each fruit fly's motion data has been obtained through time, we can thereby compute fruit flies' kinetic measurement, such as velocity v, angular velocity av and speed s which is the magnitude of velocity v at each moment.upper panel), which points to the direction of growing population size, shows the long tails of speed fluctuation.These tails are nearly exponential and grow monotonically with the population size, which suggest that fruit flies in a denser population environment perform faster manoeuvres.It probably are these long tails which cause the standard deviation of speed increases as the population size increases.On the other hand, previous studies [36,37] have demonstrated that the fruit fly exhibits a flight pattern in which straight flight sequence interspersed with rapid turn called saccades.Fig 6 shows the statistics of velocity of fruit flies in the confined arena.These results indicate that the fruit flies' kinematics are statistically similar with that of midges [23].Midges are known to form swarms under visual markers, such as stagnant water [23,24].We here, however, observed that fruit flies in the arena flew almost randomly in the arena.

Spatial structure
Even when hundreds of fruit flies were free-flying in the arena at the same time, no collision was ever observed.We speculate that the spatial organization of fruit flies may manifest the social structures.The clearest characterization of the spatial structure of fruit flies within a swarm is the probability distribution of the nearest-neighbour distance.The nearest-neighbour distance is the minimum Euclidean distance between a fruit fly and other fruit flies in a swarm.In the later text, we refer to the nearest-neighbour distance as NND.Table 1 shows the number of NNDs of each configuration.At each moment, there were many fruit flies landing on walls of the arena.Therefore, in order to avoid the bias introduced by those landing fruit flies (i.e.edge effects), we ignored a fly's motion data if its distance to any walls of the arena is less than 20 mm.Fig 7a shows the PDFs of NND distributions of all configurations.It shows that the average NND and the deviation decrease as the population size increases.Moreover, considering the population size of T05 is larger than that of T04, the similar PDFs of T04 and T05 indicate there exists an asymptotic saturation against the increasing population size.
In order to confirm whether there exists some biological forces or social forces among fruit fly swarms, we created a simulation system.The system simulated physical random particles (a.k.a Brownian motion in which social forces are absent) moving in a confined volume.The virtual volume was equal to the volume of the flight arena.The number of particles was equal to the number of fruit flies at each moment of each configuration.At the initial step, random particles were uniformly distributed in the volume.Each simulation ran on 3000 steps.We also evaluated the NND distribution of each simulation.
where x denotes the population density and y denotes the NND.The constant A determines the overall scale of the variation with x, whereas B is the decaying rate constant.The constant C gives the asymptotic value for the average NND for all configurations.The interesting thing is the constant C denotes the exclusive distance between a fruit fly and its neighbours.

Acceleration towards nearest-neighbour
For each fruit fly in a swarm, C $ 5 l À 6 l denotes the exclusive distance to other flies.If this result is significant, its effect should also be apparent in the motion statistics.We therefore measured the probability distribution of the angular direction of a fruit fly's acceleration towards its nearest neighbours.Given a reference fruit fly, we computed the angular direction of its acceleration with respect to the direction of its nearest neighbour, ñ (see Fig 8, the legend).The angular direction is defined by two angles (the azimuth and the elevation).We repeated this computation by taking each fruit fly within a configuration as reference individuals, and mapped the probability density conditioned on the NND.In

Discussion
In this paper, we detailed our tracking system which was designed for acquiring the motion data of individuals of fruit fly swarms through time.The fruit fly swarms were housed in a cubic flight arena.Three synchronized and calibrated high-speed cameras satisfy the minimum requirement of the system for resolving the ambiguity between targets.More cameras may increase the ability of resolving targets in a more dense population.The proposed tracking algorithm is easily tunable for more cameras because of the well designed measurement association algorithm.Beside, with the help of the "one-to-one" strategy, the proposed tracking algorithm is uncoupled from the population size.The computing time cost is linear to the population size.That is, the proposed tracking algorithm is suitable for tracking more targets depending the capability of computers (multi-thread).
The probability distribution of the nearest-neighbour distances shows the spatial structure of fruit flies within a swarm.This distribution is significantly different to that of random particles (non social forces simulation).The average nearest-neighbor distance of all experimental configurations show the property of asymptotically approaching saturation.Further, by evaluating the distribution of angular direction of fruit flies' acceleration towards their nearestneighbours, we found the evidence proving the existence of the asymptotical distance.Fruit flies have a strong tendency to take the repulsive response when the distances between them approached 6 l (the average body length is l ¼ 2:73 mm).This result is consistent with previous findings reported by Maimon et al. [10], who observed a strong decay in the probability of a fruit fly's approaching small post if the post subtended a visual angle of % 10°on the fly's retina.At distance between fruit flies approaching 6 l, a fly subtends a visual angles of % 10°on another fly's retina.And therefore, fruit flies exhibit a high tendency to turn away from each others.Though there exists obvious interaction between a fruit fly and other flies, swarms of fruit flies in our experiments did not show an overall polarisation.On the average, unlike bird flocks [17], here the average polarisation is % 0.02.It suggests that these fruit flies have little tendency to align their motion with their neighbours.
In conclusion, this study provided a detailed 3D tracking system for obtaining quantitative 3D motion data of individuals of the fruit fly swarms through time.The result of quantitative analysis shows that the fruit flies in a confined arena are not free particles.Behaviour rules exist in a fruit fly swarm and probably affect individual's behaviours in the swarm.Understanding the detailed origin of their behaviours will be an interesting topic for future research.
A feeding-tube matched open-top sunroof in 50 mm diameter circle was handled on escape-proof inner besieged flange.Three monochrome high-speed digital video cameras (IO Industries Canada, Flare 4M 180-CL, 2040v × 2048h pixels at 100 fps) were placed approximately 900 mm from the cubic arena against two orthogonally placed back-lit cool-running lamps.Each lamp was made by LED arrays and covered by a diffusion sheet to generate gentle and flicker-free planar illumination.The cameras were mounted with 17-35 mm lens.The supporting information (S1 Fig) shows the experimental equipment.

Fig 2 .
Fig 2. Detecting.(a) The image at the left was filmed by a back-lit camera and the patch marked by the green rectangle was zoomed at the right.(b) The patch was overlaid with red blobs which were detected using background subtraction.(c) The patch was overlaid with refined blobs in which the pixels of wings were removed.(d) The result of fitting each blob with an ellipse.doi:10.1371/journal.pone.0129657.g002

Fig 4 .Fig 3 .
Fig 4. Overview of the tracking algorithm.By exploiting the particle filtering solution, we use ñP ¼ 200 particles sampling the posterior distribution of the motion state of a fruit fly at moment t.Each particle includes two components: the fly's location (x, y, z) and the fly's orientation (θ, ϕ).The location component is predicted by using the dynamic model and the fly's location at moment t−1.After the location component is predicted, an MMP is associated with the particle.Then the orientation component is computed using the MMP.And then the tracking algorithm computes the weight of the particle.If there are more than one MMP (e.g.since the target is associated with two measurements in Camera 1, the algorithm thereby generates two MMPs (χ 1,1 , χ 2,2 , χ 3,3 ) and (χ 1,2 , χ 2,2 , χ 3,3 )), the algorithm choose the MMP having highest weight (e.g.(χ 1,1 , χ 2,2 , χ 3,3 )) and the orientation computed by the MMP.Finally, the motion state of the fly at moment t is the expectation computing by the weighted particles.doi:10.1371/journal.pone.0129657.g004

2 . 3 . 3
Orientation computation.Measurements defined by GðX n t Þ form an MMP (c.f.There may probably be more than one MMP.Please see Fig 4, the legend.).An MMP includes one measurement in each camera view.The ellipses of the MMP are employed to compute the target's orientation according to the pinhole camera model in Euclid geometry.The geometric interpretation of computing the orientation is shown in Fig 5a.

Fig 5 .
Fig 5. Orientation computation.(a) The geometric interpretation of the orientation computation.Here I 1 and I 2 denotes the image plane of Camera 1(red)and Camera 2(blue) respectively.The line-segment hx 1 1 ; x 1 2 i(red) denotes the ellipse's major axis in I 1 .And the line-segment hx 2 1 ; x 2 2 i(blue) denotes the ellipse's major axis in I 2 .The projection rays of those end-points define two plane Φ 1 and Φ 2 respectively, e.g. two red rays define Φ 1 .The orientation O n t can be computed from the direction of the cross line (the longer green line-segment) between two planes Φ 1 and Φ 2 .All definitions are defined in the world coordinate system.(b) The problematic orientation computation.The measurement is drawn in red, and red dashed lines are the major-axis and minor-axis of its ellipse.A generative shape according to the particle state ðX n t ; O n t ; lÞ is re-projected into Camera 2 (see S2 Fig, the generative shape model).The green pixels are re-projected pixels and the yellow pixels are identical pixels.(c) The error as a function of the ratio γ.Here γ defines the length ration between the major-axis and the minor-axis of an ellipse.
Fig 5. Orientation computation.(a) The geometric interpretation of the orientation computation.Here I 1 and I 2 denotes the image plane of Camera 1(red)and Camera 2(blue) respectively.The line-segment hx 1 1 ; x 1 2 i(red) denotes the ellipse's major axis in I 1 .And the line-segment hx 2 1 ; x 2 2 i(blue) denotes the ellipse's major axis in I 2 .The projection rays of those end-points define two plane Φ 1 and Φ 2 respectively, e.g. two red rays define Φ 1 .The orientation O n t can be computed from the direction of the cross line (the longer green line-segment) between two planes Φ 1 and Φ 2 .All definitions are defined in the world coordinate system.(b) The problematic orientation computation.The measurement is drawn in red, and red dashed lines are the major-axis and minor-axis of its ellipse.A generative shape according to the particle state ðX n t ; O n t ; lÞ is re-projected into Camera 2 (see S2 Fig, the generative shape model).The green pixels are re-projected pixels and the yellow pixels are identical pixels.(c) The error as a function of the ratio γ.Here γ defines the length ration between the major-axis and the minor-axis of an ellipse.doi:10.1371/journal.pone.0129657.g005

3. 1 . 1
Animals.Wild-type Canton fruit flies were fed basic medium (corn flour broth with brown sugar) 3-5 days after eclosion.The flies were housed at 25°C room temperature and 60% relative humidity under 12 hour light-dark circadian clock without sexual discrimination.
Fig 6a shows the statistics of speed.It shows that the measured mean speed μ s is greater than 400 mm/s (see the lower panel of Fig 6a), while the measured standard deviation σ s increases from 200 mm/s to 400 mm/s as the population size increases.Considering the z-score (a.k.a the standard score) of speed denotes the fluctuation of speed, we computed the z-scores of all configurations and reported the PDFs of z-scores of each configuration in Fig 6a (the upper panel).The z-score of speed z s is computed as z s = (s −μ s )/σ s .Fig 6a shows the fluctuation of fruit flies' speed usually narrows down (i.e. the amplitude of fluctuation decreases) as the population size increases.That is, the flight of fruit flies in a denser population environment is less volatile than those in a less dense environment.It suggests fruit flies have probably awareness of the environment.The black arrow in Fig 6a (the Fig 6d shows the statistic of angular velocity.The measured mean angular velocity μ av is greater than 400 degree/s (see the lower panel of Fig 6d), while Fig 6d (the upper panel)shows that the angular velocity av is less than the mean μ av in most of moments.That is, the angular velocity is usually less than 400°per seconds.But the nearly exponential tails (indicated by dashed-line in grey) suggest that fruit flies have also often taken the rapid turn.Considering the velocity components, Fig6b and 6cshow the statistics of velocity components: vy (the horizontal velocity component, the other horizontal velocity component vx is statistically the same) and vz (the vertical velocity component).Anisotropy of the motion is evident.The fruit flies move more actively in horizontal than in vertical.The PDFs of z-scores of the horizontal velocity z vy = (vy−μ vy )/σ vy and of the vertical velocity z vz = (vz−μ vz )/σ vz have similar Gaussian shapes (but not Gaussian, see the upper panel of Fig6b and 6c).Even though these PDFs deviate from Gaussian values at the tails, they suggest the mean velocity μ vy and μ vz of all configurations over time are approximately zero (see the lower panel of Fig 6band 6c).That is, fruit flies in the arena do not show an overall polarisation, as shown in the inset of Fig 6a.Here the polarisation is defined as F ¼ j P N i¼1 ṽi =v i j=N, where N is the number of fruit flies and ṽi is the velocity of fruit fly i.The polarisation measures the degree of alignment of the

Fig 6 .
Fig 6.Statistics of fruit flies' kinematics.The measured data of different configurations are color coded.The upper row shows the PDFs of measured zscores z and the lower panel shows the measured mean μ and standard deviation σ.(a) The PDFs of z-scores of speed, z s , of all configurations.The black arrow shows the direction of growing population size.The inset shows the polarisation value in six seconds (more data is not present).(b, c) The PDFs of zscores of velocity components: (b) the horizontal component, z vy , and (c) the vertical direction, z vz , of all configurations.The distributions are nearly Gaussian (an empirical Gaussian curve is shown in grey), with small deviation in the tails.(d) The PDFs of z-scores of angular velocity, z av , of all configurations.The dashed-line shows the nearly exponential long tails on the high angular velocity.doi:10.1371/journal.pone.0129657.g006 Fig 7b shows the comparison between the PDF of NND distribution of T03 and that of the correspondent simulation.The NNDs should follow a Poisson distribution for random particles.But the Poisson distribution fitting the distribution of the simulation compares poorly with that of fruit flies.The distribution of fruit

Fig 7 .
Fig 7. The NND distributions.The PDFs are represented by normalized histograms.(a) The PDFs ofNND distributions of all configurations.(b) The comparison between the PDFs of NND distributions of T03 and the correspondent simulation (random particles).(c) The relationship between the average NND and the population densities at each moment of all configurations.(d) The average NND y as a power function of the population density x. doi:10.1371/journal.pone.0129657.g007 Fig 7d shows the best fitting result.The best fitting C is C = 15.11± 0.89.It means that fruit flies have a high tendency to take the repulsive response when the distance between them approached C = 15.11± 0.89 mm.Considering the average body length of the fruit flies ( l ¼ 2:73 mm), C can also be expressed as C $ 5 l À 6 l.
Fig 8a and 8b the NND condition (between 6 l and 8 l) is omitted for clarity.Fig 8a shows the distribution in which all data were computed conditioned on the NND is < 6 l; in this way, that of Fig 8b were computed conditioned on the NND is > 8 l.Both distributions are anisotropic, but the former has three obvious clusters.The left and right clusters (see Fig 8a) indicate a strong tendency towards a repulsive response.On the other hand, Fig 8b does not exhibit special manoeuvres.Fruit flies' acceleration scattered when they met their nearest neighbours at distance > 8 l.Fig 8a shows another data cluster which is the cluster centering on h0, 0i.It indicates that there were fruit flies accelerated directly towards their nearest neighbours.This manoeuvre can be thought as representing "chasing" behaviour.

Table 1 .
Video sets and raw results.