Classification of animal dive tracks via automatic landmarking, principal components analysis and clustering

. The behaviour of animals and their interactions with the environment can be inferred by tracking their movement. For this reason, biologgers are an important source of ecological data, but analysing the shape of the tracks they record is difficult. In this paper we present a technique for automatically determining landmarks that can be used to analyse the shape of animal tracks. The approach uses a parametric version of the SALSA algorithm to fit regression splines to 1-dimensional curves in N dimensions (in practice N ¼ 2 or 3). The knots of these splines are used as landmarks in a subsequent Principal Components Analysis, and the dives classified via agglomerative clustering. We demonstrate the efficacy of this algorithm on simulated 2-dimensional dive data, and apply our method to real 3-dimensional whale dive data from the Behavioral Response Study (BRS) in the Bahamas. The BRS is a series of experiments to quantify shifts in behavior due to SONAR. Our analysis of 3-dimensional track data supports an alteration in the dive behavior post-ensonification.


INTRODUCTION
Animal-borne satellite-tags or data-loggers have become an important source of ecological data on the movement of animals (Bearhop et al. 2010. From this data behaviour and habitat preference might be inferred, particularly when referenced with environmental data (Kobayashi et al. 2008. The analysis of time-series movement data includes summary statistics (Scott and Chivers 2009), spatial models (Aarts et al. 2008) and dive-profiles in 2-dimensions (Austin et al. 2006). While the data often contain genuine 3dimensional shape information, direct analysis of these shapes is challenging. However, the analysis of shape is of particular importance in conducting behavioural response experiments, such as recent behavioural response studies which investigated the effects of sonar playbacks on beaked whales (Tyack et al. 2011). An important step in the analysis of dive data is the classification of the dive shapes into a set of dive types. Typically this is done by restricting the dive data to two variables (depth and time) and then partitioning it into a number of classes.
Examples of these classifications include: Square, V, Skewed right and Skewed left (Schreer et al. 2001); u, U and W (Halsey et al. 2007); W-shaped, V-shaped, Square dives, Drift dives, U-shape and Square Root-shape (Bailleul et al. 2007); W-dives, V-dives, U-dives, u-dives and I-dives (Elliott et al. 2008). This approach has been of use when applied to relatively simple dive shapes, such as those generated by seabirds, but has some drawbacks. For example, by omitting longitude and latitude when restricting the dive data to two variables, a significant amount of information about the shape of the dive is lost. Also the determination of dive types is subjective, and potentially difficult (especially if one were to retain all the spatial coordinates in the analysis). The motivation for our work comes from the analysis of dive data to investigate the effect of ensonification on whale behaviour. These dives are relatively long in duration and tend to involve numerous changes of direction. The source of ensonification was a sound transducer deployed from a ship, and the heading of the whale relative to this ship is important when considering the shape of a dive. Thus all spatial coordinates need to be retained in the analysis of dive shape. Also it is very difficult to determine a set of dive types that succinctly and adequately classify the collected data. For these reasons we have adopted a different approach, and determine landmark points to capture the dive shapes, and use PCA and agglomerative clustering to classify the dives.
The use of spatial landmarks is now standard in modelling shape. Their popularity derives from their effectiveness in applications within computer vision, medical imaging, and anthropometry, although their use has also been considered in more closely related fields such as the study of motion trajectories (Slice 1999, Adams andCerney 2007). Ideas implemented in Active Contour Models (Kass et al. 1988) have been extended with the aid of spatial landmarks into more general deformable shape models such as Active Shape Models (Cootes and Taylor 1992, Cootes et al. 1995, van Ginneken et al. 2002, Active Appearance Models (Cootes et al. 1998, Cootes and and M-reps (Pizer et al. 1999(Pizer et al. , 2003. These deformable models have been very effective in assisting with segmentation of images, and for the analysis of object shapes when combined with Principal Component Analysis (PCA) (Pizer et al. 1999, Azouz et al. 2005, Pilgram et al. 2005 or Independent Components Analysis (ICA) (Hyvärinen andOja 2000, Ruto et al. 2006). For a survey of deformable models see (McInerney and Terzopoulos 1996).
The task of identifying suitable landmarks (or multi-scale meshes in the case of M-reps) can be difficult and time-consuming. A number of approaches have been developed, with varying success. Equally spaced landmark points were placed on the boundary of objects in Baumberg and Hogg (1994) and used as the knots in a Bspline to fit the objects' outlines. In Wang et al. (2003) surface curvature and geodesic distance were used to automatically determine a set of landmarks from a small initial set of manually determined landmarks. A combination of Procrustes analysis and thin-plate splines was used in Bookstein (1997) to identify point-to-point correspondences (''semi-landmarks'') between continuous contours of similar shape. An extension of Procrustes analysis was presented in Rangarajan et al. (1997) for identifying point-topoint correspondences for automatically derived point sets (possibly of differing cardinalities). A method for automatically identifying landmarks from the vertices of existing triangulations of 3D surfaces was given in Brett and Taylor (2000). In Hill and Taylor (1994) a two-stage process was described that incorporated a pair-wise corresponder to determine an initial set of landmarks and an iterative non-linear optimisation scheme to improve the positioning of these landmarks. Non-linear optimisation techniques were also used in Kotche and Taylor (1998) and Davies et al. (2001)-both aimed to optimise for some measure of parsimony (minimum number of parameters for the former, minimum ''description length'' for the latter) and used a genetic algorithm to find a good solution. In this paper we present a method for automatically determining landmarks that uses a parametric version of SALSA-a regression spline fitting algorithm described in Walker et al. (2011). SALSA determines the number and location of knots in a regression spline, so as to minimise some fitness criterion. The criterion used is the Bayesian Information Criterion (BIC). Like the methods described in Kotche and Taylor (1998) and Davies et al. (2001), this amounts to a non-linear optimisation approach for identifying the most parsimonious model. As with those techniques it is impractical to identify the global optimum, so instead a heuristic approach is used to determine a good solution in reasonable time. This approach includes a global exchange step inspired by Remes Exchange Algorithm for polynomial curve approximation (whereby the location of the largest residual is considered as a new knot location), a local improvement step that moves knot locations to neighbouring data points, as well as the standard forward and backward regression steps (addition and deletion of knots). A mathematical formulation of the knot location problem and further details of SALSA are given in Walker et al. (2011).
The rest of this paper is structured as follows. In Identifying Landmarks we describe our parametric SALSA algorithm for automatically determining landmarks to describe the shape of arbitrary open or closed 1-dimensional curves in N dimensions. In Principal Components Analysis for Investigation of Shape we summarise how these landmarks can be transformed and included in a PCA to analyse objects' shapes. In Simulation Example we implement our approach on a simulated set of 2-dimensional dive data from a set of a priori dive classes to illustrate its effectiveness. In Results we demonstrate how our landmarking approach can be used to analyse 3-dimensional animal tracks on real whale-dive data, and we finish with some concluding remarks in Conclusions.

SALSA for curve parameterizations
The original version of SALSA (Walker et al. 2011) fits a B-spline to a given data set so as to minimise some fitness criterion. The user specifies the fitness criterion, and the degree d of spline to be used, and SALSA determines the number and location of the knots and fits the corresponding spline. Each whale dive is a 1dimensional curve in 3-dimensional space, so it is necessary to extend the algorithm to accommodate this sort of curve. This is done by parametrizing the x, y and z co-ordinates of the dive by a fourth parameter t (which can be thought of as time in this application, although for an arbitrary closed or open curve this approach will also work). Each of these curves, x(t), y(t) and z(t), is then fitted using the SALSA algorithm, with the regression splines used for the 3 parametric curves constrained to have knots at the same values of t. To fit 3 curves at once the fitness criterion needs to be updated to incorporate a fitness measure from each curve. It is possible to use the residual sums of squares from the 3 parametric curves to compute a BIC value for the single curve in 3D, but we have found the sum of the BIC values for each parametric curve to be an effective criterion. The exchange step of SALSA is also updated to identify the value of t that yields the largest magnitude residual across all 3 parametric curves. In Fig. 1, an example of the approach for the whale dive shown in Fig. 2 is presented.
The fitted values of the degree d regression splines for x(t), y(t), and z(t) give the coordinates of a degree d spline for the original curve in 3 dimensions. Fig. 2 shows a solution for a particular whale dive. The grey curve is the original whale dive, and the black curve is the fitted regression spline, with the black dots showing the location of the knots.
The knot points of the regression spline fitted by the parametric SALSA can be used as landmarks for an analysis of the dive shape. SALSA aims to fit the model with the lowest BIC, so it is extremely unlikely that all dives being compared will generate the same number of landmarks (a requirement for the subsequent PCA). To deal with this we do an initial fit for all dives to determine the maximum number of knots K fitted across all the resulting regression splines. SALSA allows the user to define the maximum number (max) and minimum number (min) of knots to be permitted when fitting a spline. By computing a second fit for all dives with both max and min set to K we are able to produce a comparable set of landmarks across all the dives. This approach ensures we have chosen the minimum number of landmarks for all dives while still capturing sufficient shape information to satisfy the BIC criterion. For some dives we will capture more shape information than is required, but for our application this is not a problem. One could implement other approaches if desired, such as fixing the number of knots to be the minimum or average across all the fits v www.esajournals.org from the first application of SALSA.

Alignment of shapes
Before performing PCA on the landmark points of the shapes it is necessary to eliminate differences due to rotation, translation, and scale.
One approach is to use the Procrustes metric to align the shapes' landmarks (see Kendall 1984or Bookstein 1997. Typically, the landmarks for each shape are translated so their centroid lies at the origin, are scaled so that the sum of squared distances for each shape from the landmarks to the centroid is 1, and are rotated so as to minimize the sum of squares of distances between corresponding landmarks (across the Fig. 1. Parameterization of a whale dive: x, y, and z co-ordinate given in terms of time. A fit is performed using SALSA with degree 1 splines. The solution places 26 knots (shown here as dots), which can be used as landmarks for the original curve in 3D.
v www.esajournals.org shapes). It is possible that the correspondence between landmarks is determined as part of this final rotational step. In our application there is a natural correspondence between the landmarks as they are uniquely ordered by the whale dive. We align the landmarks using the first two operations described above (translation and scaling), but the final rotational step is restricted to rotations about the z-(depth) axis as all the whale dives start and end at the surface of the sea (an orientation that we do not want to change). It should be noted that the scaling step can be omitted if the goal is to include size variation in the subsequent PCA analysis.

Principal Components Analysis
The aligned landmark points (x 1,i , y 1,i , z 1,i ) ... (x n,i , y n,i , z n,i ) for each of the i ¼ 1,..., m shapes are combined as rows of a matrix where X T i (x 1,i , y 1,i , z 1,i , ... x n,i , y n,i , z n,i ). PCA is applied to X by computing the eigensystem of its covariance matrix. The resulting eigenvectors identify the directions of the major axes of the data cloud, and are ordered so that the first axis captures the maximum possible variation in the data, followed by the second, and so on. The ith eigenvalue gives the variance of the data projected onto the ith axis, and so can be used to determine the percentage of the data variation captured by each new axis (or ''component''). It is common to reduce the dimension of the problem to d by considering only the first d components. The mean shapeX can be used to approximate each shape: where U d is the matrix with the first d eigenvectors as its columns, and b is the difference in ''score'' (projection) between the mean shape and X i for the first d components, given by b ¼ U T d ðX i ÀXÞ: The choice of d is usually made based on an a priori cut-off for the percentage of data variation that should be captured by the model (for example, 90%).

SIMULATION EXAMPLE
In the simulation study presented here, three underlying dive shapes are considered: V, U and W (Schreer et al. 2001, Bailleul et al. 2007, Halsey et al. 2007, Elliott et al. 2008. We generated the V, U and W dives using cubic B-splines (via the bsfunction in R) with 1, 2 and 3 knots respectively. The function domain for all dives was [0, 1]. The coefficients used to generate the simulated dives are shown in Table 1.
The knot positions were permitted to vary normally with r ¼ 0.02 across dives, while the coefficients were also generated from a Normal distribution with mean vector shown in Table 1 and r ¼ 0.07. A random normal error term (r ¼ 0.01) was also added to each dive. These values were chosen to mimic the variability seen in observed whale dives. The set of generated dives is shown in Fig. 3.
Landmarks were assigned and PCA component scores calculated following the method described in Principal Components Analysis for Investigation of Shape. Agglomerative clustering was performed on the PCA component scores for the whale dives using the function hclust in R (R Development Core Team 2009). The dendrograms shown in Fig. 4 were obtained using average linkage to calculate the distance between clusters. The dendrogram on the left is labelled by dive type: 1 for U, 2 for V, and 3 for W. On the right is the same plot with the dive indices as labels. A natural distance at which to define clusters would be just above 0.25, resulting in: all the U dives in one cluster; three clusters containing only W dives (one containing only dive 6); two clusters containing only V dives; and 1 cluster that contains two V dives and one W dive (dive 21). Ignoring dive 21, all clusters contain only one type of dive, but more clusters than dive types result, due to the stochastic nature of the generation process. Dive 21 has a very shallow peak, which probably explains it's classification with two V dives. It appears that dives are being grouped according to the similarities of their shapes.

RESULTS
We have tested our approach on dive data for one tagged whale over one day, collected as part of the Behavioral Response Study (BRS) in the Bahamas (Tyack et al. 2011). The data was separated into 51 dives, one of which (referred to here as dive 9) occurred while the whale was ensonified with simulated Mid-Frequency Active Sonar (MFA) pulses from the study vessel. Killer whale (Orcinus orca) calls were broadcast from the study vessel during dive 19, but for all other dives no transmission was made.
The 51 dives were fitted using the parametric SALSA algorithm described in Identifying Landmarks. Each dive was fitted with 26 internal knots. A further knot was included at the start and end of each dive, resulting in 28 landmarks and 84 variables to be included in the subsequent Principal Components Analysis (PCA). The variables were transformed and an eigenanalysis conducted on their covariance matrix using the prcomp command from the base package in R (R Development Core Team 2009), as outlined in Principal Components Analysis for Investigation of Shape. A screeplot of the resulting eigenvalues is shown in Fig. 5. From a visual inspection of the plot it appears that the first six components should be used; components 7 and 8 should not be used as the corresponding eigenvalues are approximately equal and may indicate sphericity (in which the spread of the data in these 2 dimensions is approximately circular and so the associated eigenvectors are determined arbitrarily). Note: the rest of the screeplot changes smoothly, so this blip at component 8 may warrant further investigation. The first 6 compo- Fig. 4. Dendrograms showing hierarchical clustering of the PCA scores for the simulated 2-dimensional dive data. The dendrogram on the left has leaves labelled with dive types: 1, 2, and 3 represent U, V, and W dives respectively. The dendrogram on the right has leaves labelled with dive indices.
v www.esajournals.org nents account for 92.66% of the variation in the data.
Reduced space plots for all pairings of the first 6 components are given in Fig. 6. Dive 9 (the dive during which ensonification occurred) is shown as a heavy dark point. For the first 4 components this point is clearly isolated in the reduced space plots. These plots are based on data from one whale for one day, so it would be premature to draw any biological conclusions from this analysis, but these plots certainly indicate that the alteration to whale behaviour from ensonification is a possibility.
Using the function hclust in R (R Development Core Team 2009) we are able to perform agglomerative clustering on the PCA component scores for the whale dives. With the average linkage for determining the distance between clusters we obtain the dendrogram shown in Fig.  7.
The ensonified dive (dive 9) is clearly distinct from the other dives. It is one of 3 deep dives Fig. 5. Eigenvalues resulting from the PCA conducted on the transformed landmark data. The first six components account for 92.66% of the data variation.
v www.esajournals.org analysed, but this variability in scale across the dives has been removed from the data in the preprocessing step (see Identifying Landmarks). Thus the difference detected in our analysis is due to the dive shape. Closer analysis shows that unlike the other dives that are primarily headed in one direction, dive 9 includes a sharp change of direction a short time after the termination of ensonification and heads away from the transmission vessel, as shown in Fig. 8. Furthermore, the eight dives prior to ensonification (dives 1 to 8) appear in their own branch of the dendrogram, suggesting that the shape of the dives after ensonification is different to these. It is also worth mentioning here that, despite scaling to remove size effects, the dendrogram supports the classification of the other two deep dives (dives 19 and 37) into classes separate from the more shallow dives. Dive 19 is the only dive, other than dive 9, receiving audio from the transmission vessel; orca calls were broadcast during this dive. There is no obvious feature differentiating the shapes of

CONCLUSIONS
By implementing a parametric version of the SALSA algorithm we are able to automatically determine landmarks for animal tracks in N dimensions. These landmarks can then be included in a Principal Components Analysis to investigate differences in shape across the tracks. Using 2-dimensional simulated dive data, we have demonstrated the effectiveness of this approach for clustering animal tracks according to their shape. We have also applied our approach to real 3-dimensional data for whale dives, the results of which illustrate that the technique is useful for identifying changes in whale behaviour due to external stimuli. This analysis supports the possibility that ensonifica-tion affects whale behaviour. However, a more comprehensive set of whale dive data across multiple animals would be required before an inferential analysis could be conducted.

ACKNOWLEDGMENTS
The authors wish to thanks Ian Boyd, Brandon Southall and Peter Tyack for helpful comments on the manuscript. The BRS study was financially supported by the United States (U.S.) Office of Naval Research (www.onr.navy.mil) Grants N00014-07-10988, N00014-07-11023, N00014-08-10990; the U.S. Strategic Environ- Fig. 8. The animal track for dive 9. The period of ensonification is shown in black. The period prior to ensonification is shown in light grey, and the period after ensonification is shown in a darker grey. At the commencement of ensonification the animal descends. A short time after the termination of ensonification the animal heads off in a direction away from the transmission vessel.
v www.esajournals.org mental Research and Development Program (www. serdp.org) Grant SI-1539, the Environmental Readiness Division of the U.S. Navy (http://www.navy.mil/local/ n45/), the U.S. Chief of Naval Operations Submarine Warfare Division (Undersea Surveillance), the U.S. National Oceanic and Atmospheric Administration (National Marine Fisheries Service, Office of Science and Technology) (http://www.st.nmfs.noaa.gov/), U.S. National Oceanic and Atmospheric Administration Ocean Acoustics Program (http://www.nmfs.noaa.gov/ pr/acoustics/), and the Joint Industry Program on Sound and Marine Life of the International Association of Oil and Gas Producers (www.soundandmarinelife. org). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.