The zebrafish Danio rerio is a widely used model organism in studies of genetics, developmental biology, and recently, biomechanics. In order to quantify changes in swimming during all stages of development, we have developed a visual tracking system that estimates the posture of fish. Our current approach assumes planar motion of the fish, given image sequences taken from a top view. An accurate geometric fish model is automatically designed and fit to the images at each time frame. Our approach works across a range of fish shapes and sizes and is therefore well suited for studying the ontogeny of fish swimming, while also being robust to common environmental occlusions. Our current analysis focuses on measuring the influence of vertebra development on the swimming capabilities of zebrafish. We examine wild-type zebrafish and mutants with stiff vertebrae (stocksteif) and quantify their body kinematics as a function of their development from larvae to adult (mutants made available by the Hubrecht laboratory, The Netherlands). By tracking the fish, we are able to measure the curvature and net acceleration along the body that result from the fish's body wave. Here, we demonstrate the capabilities of the tracking system for the escape response of wild-type zebrafish and stocksteif mutant zebrafish. The response was filmed with a digital high-speed camera at 1500 frames s–1. Our approach enables biomechanists and ethologists to process much larger datasets than possible at present. Our automated tracking scheme can therefore accelerate insight in the swimming behavior of many species of (developing)fish.

Keywords:
tracking, geometric modeling, estimation, zebrafish, swimming

Capturing the detailed locomotory kinematics and behavior of biological organisms plays an important role in a wide variety of research disciplines(Fry et al., 2003; Gray, 1933; Park et al., 2001; Pfau et al., 2005). Many studies in biomechanics, neuroethology and developmental biology rely on the analysis of video sequences. Such analyses stand to benefit greatly from computer vision techniques, which can accurately and efficiently screen complex morphometric and locomotory characteristics.

Following the pioneering work of Eadweard Muybridge(Muybridge, 1887), scientists have developed a wide range of methods to record and quantify animal locomotion. Many studies have relied on tracking landmarks, such as leg joints(Full and Tu, 1991) or markers on the animal's body (Hedrick et al.,2004; Standen and Lauder,2005). However, it is not always feasible or advisable to attach markers, e.g. during field observations or on animals that have a mucous layer(e.g. many fish) or motile skin (e.g. horses). In these cases, tracking the silhouette of the animal or a body part can offer an alternative [e.g. silhouettes of ascidian larvae and swimming fish(McHenry, 2001; Tytell and Lauder, 2002)].

Tracking swimming movements by filming the swimmer's silhouette has a long tradition (Batty, 1984; Gray, 1933). Silhouettes lend themselves to automatic tracking (Tytell and Lauder, 2002). Furthermore, tracking routines can be refined by using information about the shape of the animal or limb to be tracked [e.g. insect wings (Fry et al.,2003; Willmott and Ellington,1997)]. Error reduction as a result of including shape information is particularly valuable when studying the propulsive body wave of fishes. The exact shape of a fish's body wave depends on the interaction between the solid body of the fish and the surrounding water. Exact knowledge of the wave shape can provide valuable insight into swimming performance and the interaction between undulating body and surrounding water(Cheng and Pedley, 1998; McHenry et al., 1995; Wainwright, 1983). For example, the stiffness of the fish's body affects the shape of the body wave and thereby influences the propulsive performance of the fish. An interesting case study is presented by the stocksteif mutation in zebrafish,which is characterized by an overossification of the notochord. The axial skeleton of these mutants is a stiff, bony rod, contrasting the flexible series of articulating vertebrae in their wild-type siblings.

The zebrafish has long served as a convenient model to study the various aspects of fish swimming (Fuiman and Webb,1988; Thorsen et al.,2004). Automated tracking and analysis systems have been previously developed for zebrafish (Bang et al., 2002; Blaser and Gerlai,2006). However, these systems track the fish only as a point and cannot quantify body wave kinematics of swimming. Other studies of zebrafish swimming have manually tracked the fish to quantify the body posture, a method that is both time-consuming and potentially prone to subjective errors(Budick and O'Malley, 2000; McElligott and O'Malley, 2005; Müller and van Leeuwen,2004). Tytell and Lauder used a semi-automated method to estimate the fish midline by manually indentifying the snout and tail and automatically estimating the midline from the extracted silhouette(Tytell and Lauder, 2002). Other authors have relied on `skeletonizing' algorithms that dissolve a binary image representing the animal's silhouette down to its midline(Cronin et al., 2005; Geng et al., 2004; McHenry, 2001). Although this approach is automated, it will not estimate the correct midline if other objects are present in the binary image because it cannot distinguish between pixels that belong to the animal's silhouette and those that belong to a different object. As a result, these algorithms will not correctly estimate the fish's body posture when there is environmental clutter such as other fish, plants, or a hair used to initiate behavioral responses, as we did in our recordings. Automated kinematic analysis of multiple zebrafish larvae was recently demonstrated. However, this particluar analysis technique utilizes an image filter that is customized for the appearance of zebrafish larvae of a specific age (Burgess and Granato,2007). This technique does not extend nicely to zebrafish of different ages, other fish species, or when environmental clutter is present.

Here, we present a complete method for accurately and efficiently quantifying the body posture of zebrafish and other organisms with symmetric medial profiles. Our approach directly models the shape of the animal and utilizes locations of high contrast in the image to estimate its posture. The posture estimate is calculated using techniques that remain robust to clutter. The detailed swimming motion is estimated based on dorsal images of the fish recorded at sufficiently high frame rate(Harper and Blake, 1989),which enables a quantitative evaluation of the animal's kinematics and dynamics (provided the mass distribution of the animal is known). In the first section, we review the overall approach for performing model-based visual tracking. Subsequently, we develop a detailed geometric model for the zebrafish, including appropriate motion and measurement models that use information from the previous and current frame to estimate the fish's current position and posture. Finally, we demonstrate the capabilities of our tracking approach on zebrafish performing an escape response at three stages during their development (from larvae to juvenile).

Zebrafish (Danio rerio Hamilton 1822) eggs were collected after mating one stocksteif heterozygous female with one stocksteif heterozygous male. The batch of eggs contained both stocksteif mutant and wild-type embryos, but the mutant phenotype does not become apparent until 5 days post fertilization (d.p.f.). The embryos were reared at the optimal rearing temperature of 28°C. After hatching,the embryos were fed Paramecium (5 and 6 d.p.f.) and Artemia(from day 7 onwards). A fast startle response was recorded at 5, 15 and 28 d.p.f. for both wild-type and stocksteif mutant animals using a high-speed video camera (Photron, APX RS, 1500 frames s–1,1024×1024 pixels, exposure time 1/8000 s) fitted with a 105 mm Nikon lens. The startle responses were elicited by touching the animals with a horse hair. We analyzed recorded sequences from the initiation of the escape response to the moment when the fish either leaves the field of view or ceases active swimming. Our sequences therefore include stages 1 and 2 of an escape response, and usually several tail beats that are part of stage 3(Weihs, 1973) (Movie 1 in supplementary material).

Model-based image tracking and nonlinear estimation

Our goal was to quantify the body wave and swimming kinematics of a zebrafish from video sequences recorded from a top view. First, we built a geometric model of the zebrafish that is defined by the vector p, which contains the parameters that encode the body wave shape and location of the model. Then, tracking is performed by recursively estimating the fish parameters pk from measurements in the image, zk at time step k. This tracking approach is formulated within a discrete time (k), dynamic state space, framework as:
\[\ p_{k}=f(p_{k-1},{\xi}_{k-1}),z_{k}=h(p_{k},v_{k}),\]
(1)
where f is the motion model, h is the measurement model andζ k and νk are independent and identically distributed noise processes (models of the noise present in our image sequence). From a Bayesian perspective, this corresponds to estimating the probability density function P of the fish parameters given all of the image measurements up to the current time step (i.e. the posterior) via Bayes' rule:
\[\ P(p_{k}|z_{1:k})=[P(z_{k}|p_{k})P(p_{k}|z_{1:k-1})]{/}P(z_{k}|z_{1:k-1}).\]
(2)
The optimal estimate of the fish's shape and position will then correspond to the conditional mean of this probability density function, p̂k=E[pk|zk]=∫ pkP(pk|z1:k)dpk. In general, this recursive solution to the tracking problem is intractable and approximate solutions must be used instead. These approximation methods can primarily be broken into two categories: those that assume the probability density functions are normal distributions [e.g. Kalman Filters(van der Merwe and Wan, 2003)]and those that allow the probability density functions to assume arbitrary distributions [e.g. particle filters(Doucet et al., 2001)]. Although particle filters are able to solve very general estimation problems,there are many problems within visual tracking when the normal assumption holds, in which case Kalman filters provide accurate solutions and computational efficiency. The present study demonstrates that the normal assumption holds for estimating the location of zebrafish from top view images within a constrained laboratory that contains few environmental occlusions,hence we adopt Kalman filtering.

Recently, the Kalman filter has been considerably improved by statistical linearization and now performs well when applied to nonlinear motion and measurement models (Ito and Xiong,2000; van der Merwe and Wan,2003; Nörgaard et al.,2000; Sibley et al.,2006). These Sigma Point Kalman Filters (SPKF) recursively estimate the optimal fish parameters and their uncertainty by taking the weighted average of a set of regression points drawn from the probability density function that describes the fish's shape and location at the previous time step. In this study, we leverage the power of SPKFs to solve our nonlinear estimation problem. However, in order to estimate the fish parameters, we must also develop equations for the motion model f and the measurement model h that make sense for a swimming zebrafish during development. In the following sections we apply a two-step approach. First we construct a flexible geometric model of a zebrafish outline. Second,we develop appropriate motion and measurement equations for this model, which are used in conjunction with the SPKF to automatically optimize the fish parameters based on the processed image data.

Geometric modeling

The motion of zebrafish in the shallow water of a tank is largely planar,so we model the fish motion as restricted to two dimensions and assume orthographic projection in the camera model. We also assume that the zebrafish maintains a constant length and width profile during swimming. Therefore, we model the centerline of the fish as a planar, inextensible curve, which can only undergo bending. The body shape can therefore be described effectively by the bend angle
\({\Theta}\)
of the fish as a function of distance along its length, u. This function,
\({\Theta}\)
(u), is finitely parameterized using a linear combination of known basis functions where
\({\Phi}_{j}^{k}(u)\)
⁠, is the jth basis function of order k. The current implementation uses eight(N
\({\Theta}\)
=8) periodic cubic B-splines as the basis functions, which we found to be sufficient to capture the different bending modes of the swimming zebrafish and mating Caenorhabditis elegans(Fontaine et al., 2006). Other choices in the order and number of basis functions are possible. However,determining the minimum set of basis functions needed to describe all possible fish configurations is beyond the scope of this study. To increase the tracker's robustness and accuracy, our model assumes that fish have a stiff head. This assumption prevented the tracker from creating unrealistic bending deformations in the head region. This simplification is based on our experimental results (Müller and van Leeuwen, 2004), which show that the head region of freely swimming zebrafish undergoes negligible bending (zero local curvature). This corresponds to
\({\Theta}\)
(u) remaining constant in the head region, which we define as starting at the snout and extending to 20% of the fish's body length (length from snout to tail tip or caudal penduncle in the case of juvenile fish, see Fig. 5); in other words, the relative length of the stiff head region compared with the total length of the fish L is γ=0.2. We keptγ fixed for our experiments and found that γ=0.2 provided good tracking results in all three considered age groups. We implemented the head region mathematically by defining the origin of u (where uis zero) at a distance of 0.2 L behind the snout, such that the head region is described by positive u values(0≤u<0.2L) and the tail region by negative uvalues (–0.8L<u≤0). This approach simplifies the formulation of the basis functions
\({\Phi}_{j}^{k}\)
(u) and the corresponding local bend amplitudes αj along the fish's body:
\[\ {\Theta}(u)=\begin{array}{ll}{{\sum}_{j=1}^{N_{{\Theta}}}}{\alpha}_{j}{\Phi}_{j}^{k}(u)&u{\leq}0\\{{\sum}_{j=1}^{N_{{\Theta}}}}{\alpha}_{j}{\Phi}_{j}^{k}(0)&u{>}0\end{array}.\]
(3)
Fig. 1B illustrates these concepts, i.e. the definition of u as well as how the B-spline basis functions describe the body wave. The spline bases have local maxima in the tail region, but become constant in the head region.

Our complete fish model is a symmetric planar mesh with grid parameters u and v, which define a local coordinate frame in the fish body in the length and width directions, respectively. The points in this mesh[x(u,v) y(u,v)]T are continuous functions of u and v, so we can evaluate them at whatever discrete grid we choose. For our experiments, we sample 30 uniformly spaced points in the domain [–(1–γ)L,γ L] for u and three uniformly spaced points in the domain [–1, 1] for v. Because the number of sampled points in the (u,v) grid remains constant during tracking, the mesh points that make up our fish model become functions of the position and shape parameters, p. We construct the centerline of our symmetric mesh by integrating the unit tangent vector

\({\vec{e}}_{1}\)
along the grid parameter u, and the width is created by expanding the centerline in the direction of the normal vector
\({\vec{e}}_{2}\)
according to the value of R(u), which is the fish's width as a function of distance along its length. R(u) is defined as a fourth-order periodic B-spline function using 20 basis functions, and its value is calculated from the first frame of the video recording (see `Initial detection of zebrafish')and held fixed during tracking. This process is illustrated by Fig. 2 and further details are provided in the Appendix. We denote the complete fish model as H(p) where
\(p[{\vec{{\alpha}}}{\vec{T}}]^{T}\)
are the fish parameters, which include the bend angle amplitudes
\({\vec{{\alpha}}}\)
defined earlier and
\({\vec{T}}\)
⁠, the global translation vector of the entire fish.

Creating deformable models based on medial profiles has been used in segmentation problems in medical imaging(Hamarneh and McInerney, 2001)and for tracking multiple C. elegans from microscopy images(Fontaine et al., 2006; Roussel et al., 2007). Because zebrafish are laterally symmetric about their body axis, our medial profile representation offers several advantages for the tracking framework. Each B-spline basis function is defined only over a subregion of the fish body. Therefore, the local bend amplitudes, αj, have local control over the degree of bending in the body. This property is analogous to the fish anatomy, where contractions of individual muscles affect the bending over subregions of the fish's body. In summary, our parameterization of the centerline has few degrees of freedom, requires no training data, and offers a natural and anatomically sound way to constrain the fish's length and designate certain regions as stiff.

Motion model

The motion model of the fish predicts the fish parameters at the current time step based on the parameters calculated at the previous time step. By predicting the motion, it provides a better initial estimate of the fish parameters before they are updated using the measurement model. We assume that the fish movement is a combination of undulatory motion along its body and a displacement of the whole body along the centerline axis. The undulatory motion of steady swimming in fish consists of a traveling wave of increasing amplitude from head to tail. Given our model of the zebrafish, this motion is governed by the change in local bend amplitudes, αj,from the previous time step to the current one. Given that our video sequences are filmed at a sufficiently high frame rate, we make the assumption that the current local bend amplitudes are equal to the previous ones plus a random variable drawn from a zero mean normal distribution. This is a simple way to allow in the model that we are more uncertain about the body wave shape of the fish since it has moved between frames. We further assume that the fish has constant velocity, η1, between frames during displacement along its body axis that is corrupted by acceleration noise, η2. This velocity is also estimated and thus incorporated into the set of fish parameters

\(p=[{\vec{{\alpha}}}{\vec{T}}{\eta}_{1}]^{T}\)
⁠. The complete motion model calculates the current fish parameters after the fish has undergone a total axial displacement ofη 0=η1Δt+η2(Δt2/2)from the previous frame, where Δt is the time step between subsequent frames (inverse of the camera frame rate) and η2 is modeled as a zero mean random variable drawn from a normal distribution with fixed variance. An overview of the motion model is shown in Fig. 3 with displacements exaggerated for illustration purposes.

Measurement model

The present tracker system assumes that we are able to measure the location of the fish boundary, so our measurement model consists of the boundary points, qi=[xiyi]T,along with their outward normal vectors,

\(n_{\mathrm{i}}=[n_{\mathrm{i}}^{x}n_{\mathrm{i}}^{y}]^{T}\)
along the fish's outline. In order to fit our geometric model to the image at each frame, we must take the appropriate image measurements and match them with corresponding locations in the model. To achieve this, we first segment the image using background subtraction, which produces a binary image that is used to search for edges. Next, we apply a one-dimensional edge-detector filter in the direction normal to the boundary at each of the boundary points in the fish model (Blake and Isard,1998). The distance between qi and the corresponding detected edge point, ri, is projected onto ni so that the error minimized by the SPKF is the normal displacement between the edge points and model points. This process is illustrated in Fig. 4 where the initial estimate, edge points, and final solution are overlaid on an actual image. By using the SPKF, we minimize this error to obtain an updated estimate of our fish parameters p̂k(see Appendix for details). At age 28 days, the zebrafish has fully developed pectoral and caudal fins. These fins can cause incorrect tracking because the lighting conditions can make them appear as solid as the fish's body. However,by modifying the fish model to not take edge measurements in the pectoral and caudal regions, we are able to accurately estimate the body posture of the juvenile fish (Fig. 5).

Initial detection of zebrafish

Any tracking algorithm relies on an initial estimate of the object location. To achieve this, we have developed a semi-automated initialization routine that operates on the first movie frame and extracts three important pieces of information for tracking: (1) an estimate of the background image,(2) the initial fish parameters p0, and (3) the width function R(u). The background model is estimated by selecting a region around the fish and erasing it using the built-in Matlab™ function `roifill', which smoothly interpolates inward from the pixel values on the boundary of the user-defined region. Next, the background image is used to segment the first movie frame, and the Matlab™ function `bwboundaries' calculates the fish boundary from the resulting binary image. The user is then requested to click on the snout and tail locations of the fish, which allows us to divide the boundary into the left and right discrete boundaries of the fish denoted BL(i) and BR(i),respectively.

The initial centerline of the fish will be the curve that is equidistant from the left and right boundaries. To determine this centerline, we use a modification of the integral area distance used in(Roussel et al., 2007) which iterates and converges to the fish centerline. Our method works by first finding the closest point in BL to each point in BR and vice versa. Then, we calculate the median locations between each set of corresponding points in BLand BR. These median locations are assigned as the new left and right boundaries and this process is repeated until the boundaries converge onto the true discrete centerline C(i) (see Appendix for more details). The estimated centerline curve C(u) is determined by fitting a B-spline through the points calculated from C(i), and is illustrated in Fig. 6A along with the original left and right boundaries it was derived from. The bend angle function

\({\Theta}\)
(u) is calculated from C(u), and initial bending amplitudes are calculated by projecting the bend angle function onto the basis functions from Eqn 3 and illustrated in Fig. 1B. Once the centerline of the fish has been calculated, we can estimate its optimal radius function R(u) by minimizing the normal displacement between the extracted image boundary BL,R and the model boundary dictated by the radius function in the manner described in Fig. 2. The result of this minimization is shown in Fig. 6B and demonstrates that we are able to reconstruct accurately the width profile of the fish.

In this section we present automated tracking results for wild-type and stocksteif mutant zebrafish at 5, 15 and 28 d.p.f. Fig. 7 shows the raw centerlines of the zebrafish estimated by the tracker at fixed time intervals and demonstrates the quality of the proposed method. Our method successfully tracks fast escape responses of fish larvae(Fig. 7A,D) despite occasional partial occlusions [in this case, the hair used to induce the escape response(Fig. 7E)]. Furthermore,tracking takes an average computation time of 5.5±1.7 s frame–1 on a 3.0 GHz Intel® Xeon processor, which enables us to analyze large datasets very quickly. To calculate an upper bound on the accuracy of our tracking approach, we created a synthetic movie sequence by rendering our model along a known trajectory. We then tracked this synthetic sequence and measured the average error between the known and estimated centerlines over time (Fig. 8). We were able to localize the synthetic data within 0.5% of the body length on average; our errors in estimating the real movie sequences will be slightly larger because they contain additional noise sources not present in the synthetic one.

From these centerlines, we wish to extract important kinematic parameters to gain insight into developmental influences on the propulsion mechanisms of swimming fish, and vice versa: the mechanical influences on the development of the fish. One of these parameters is body axis curvature, which provides information about the muscle strains the fish undergoes. To measure curvature, we apply spatial smoothing to the extracted centerlines and then apply temporal smoothing directly to the curvature values. Spatial smoothing is performed by fitting a cubic B-spline curve to the extracted centerlines,and the curvature is calculated directly from these smoothed curves (see Appendix for details). We calculate curvature from these smoothed centerlines instead of deriving it from the raw centerline of our geometric model because the model contains a discontinuity in the curvature at the location where it becomes stiff.

Temporal smoothing is performed by applying a low-pass filter to the curvature values across time at fixed locations along the fish's body (we use 51 uniformly spaced points). The cut-off frequency for the filter is chosen based on visual inspection of the magnitude response of the curvature's Fourier transform at each of the body locations. Fig. 9 plots the error between the filtered and unfiltered centerlines and curvature values for both wild-type and stocksteif zebrafish. With average errors around 0.1%of body length and 0.1 for the centerline and curvature, respectively, our filtering process retains all-important information and does not compromise the accuracy achieved by the tracker. The resulting curvature profiles after filtering are shown in Fig. 10. Comparing stocksteif with wild type at age 15 and 28 days (Fig. 10B,E, Fig. 10C,F) we find that the peak curvature values of stocksteif are smaller than those of wild type.

Analysis of swimming performance

To analyze the performance of our fish, we measure several additional kinematic parameters. The angular acceleration is the second temporal derivative of the fish bend angle function,

\({\ddot{{\Theta}}}\)
(u). In Fig. 11, we observe similar peak angular accelerations between wild type and stocksteif at 5 d.p.f., when the mutant phenotype has just become manifest. However, as the fish get older, large discrepancies appear in the angular accelerations. At 28 d.p.f., the peak angular accelerations of the wild type are two orders of magnitude larger than that of the stocksteif. This trend is also present in the tail beat frequency of the fish. To develop an objective method for estimating the tail beat frequency of the fish, we calculated the Fourier transform of the curvature values during continuous swimming at equally spaced locations along the fish's body. This approach does not rely on determining the tail's lateral displacement from a mean path of motion, and is therefore invariant to the spatial trajectory of the fish. The time period of continuous swimming was manually determined by inspecting the curvature profiles for regions where wavespeed remained relatively constant (see Fig. 10). Fig. 12 plots the magnitude of the frequency response along the body axis of the fish. For all fish, the body axis position at approximately 90% posterior to the snout has the largest frequency response. The tail beat frequency f is calculated by taking a weighted average of the frequencies with the maximum response at each location along the fish's body axis (see Appendix for precise definition). Again, similar tail beat frequencies are observed at 5 d.p.f. However, the 15 and 28 d.p.f. stocksteif have smaller tail beat frequencies than the wild type. We estimate the curvature wave speed by performing a linear fit to the points of zero curvature during continuous swimming (see Fig. 10), and then calculate the resulting wavelength given the tail beat frequency provided by the Fourier analysis using Eqn A7. A summary of these values is provided in Table 1.

Typical measures of escape swimming performance include displacement,speed, and acceleration of the fish center of mass (COM) when stretched straight (Domenici and Blake,1997; Walker,1998). To estimate the location of the COM, we reconstructed the fish volume from dorsal and lateral photographs assuming an elliptic cross-section (McHenry, 2001)(see supplementary material Figs S1–S3). The center of volume (COV) and center of area from the dorsal view (COA) were calculated. The COV will correspond to the COM assuming the fish has uniform density. COV calculations were only performed on wildtype zebrafish for comparison with published literature on escape responses. However, COV calculations do not lend themselves to high-throughput analysis because separate photographs of each individual fish must be acquired to account for the variation in morphometry for fish at a given age (very large variation exists for the stocksteif at a specific age). Instead, we propose to use the COA as a location for comparison between wild type and stocksteif because it is easily measured from the video sequence and has similar speed and acceleration profiles to the COV (see Fig. 13). Fig. S4 in supplementary material demonstrates that the COA location has a maximum deviation of 5–6% body length from the COV at age 5 days and the deviation decreases as the fish gets older. Fig. 13 illustrates the results of these measurements. Speed and acceleration are calculated using the MSE spline method (Walker,1998). The stocksteif consistently exhibits lower peak acceleration and speed compared with wild type at each age. Although these preliminary measurements indicate consistent discrepancies between the wild type and stocksteif due to the stiffer vertebral column present in the mutant, a larger sample must be analyzed to produce statistically significant results, such an analysis will be presented in our next paper.

We have presented a method for estimating the body posture of zebrafish within laboratory environments using flexible geometric models and nonlinear estimation. Given the generalized mathematical framework used to model the fish's appearance, this method should track any fish species with a symmetric medial profile and that swims by undulating the body. The discrete time,dynamic state, space model also provides a general framework for performing statistical inference that is robust to outliers and enables tracking during partial occlusions. This is a strong improvement over previously developed methods that are either manual (Budick and O'Malley, 2000; McElligott and O'Malley, 2005; Müller and van Leeuwen, 2004), require a perfectly segmented image with no environmental clutter (Cronin et al.,2005; Geng et al.,2004; McHenry,2001), or are customized for fish of a specific size and appearance (Burgess and Granato,2007).

The assumption of planar motion is a limitation in the proposed method,which arises directly from the recorded material – top-view shots with a single camera. For behaviors that contain large out-of-plane motions, the geometric model will not accurately represent the appearance of the fish. For such behaviors, a 3D version of our model-based tracker would be required; our principle is sound, but should be extended. The tracker may also fail if a significant portion of the fish becomes completely occluded by environmental clutter (e.g. the entire head). Nevertheless, this could be improved by extending the algorithm to use a more advanced nonlinear estimator, a more advanced motion model, or a direct model of the occlusions. However, in laboratory settings, where many environmental parameters are controlled, this technique represents an accurate and fully automatic approach to quantify behavior and will facilitate studies requiring the analysis of many and long image sequences.

In addition to traditional swimming performance indicators we explored two new ways of analyzing the kinematics data: by plotting the angular acceleration as a function of time and the frequency response along the body. We also used more objective mathematical definitions and corresponding algorithms to quantify standard variables such as tail beat frequency, wave speed and length. The center of area (COA) of the fish dorsal view was proposed as a valid location for comparison between wild-type and stocksteif fish because its ease of measurement lends itself to high-throughput analysis. Our preliminary comparison between wild-type and stocksteif swimming performance indicators already suggest significant differences. Hence, this preliminary data analysis of swimming fish illustrates the capabilities of the automatic fish tracker and bodes well for gaining a complete understanding of how stiffness of the vertebral column affects swimming performance.

Here, we present some of the detailed equations related to the geometric,motion and measurement models used in our approach for tracking zebrafish. The tangent and normal vectors of the model's centerline are:
\[\ {\vec{e}}_{1}(u)=\left[\begin{array}{c}\mathrm{cos}({\Theta}(u))\\\mathrm{sin}({\Theta}(u))\end{array}\right],{\vec{e}}_{2}(u)=\left[\begin{array}{c}-\mathrm{sin}({\Theta}(u))\\\mathrm{cos}({\Theta}(u))\end{array}\right],\]
(A1)
and the complete fish model is given by:
\[\ H(p)={\beta}\left({{\int}_{0}^{u}}{\vec{e}}_{1}({\hat{u}})\mathrm{d}{\hat{u}}+vR(u){\vec{e}}_{2}(u)\right)+{\vec{T}},\]
(A2)
where
\({\Theta}(u)\)
(u) is defined in Eqn 3, H(p) is the complete fish model, and β is a constant scaling term (pixels mm–1) to scale the model to the appropriate size depending on the camera magnification. For our experiments, the scaling term β is calculated from the test image of a calibration grid.
The equations for our motion model predict the current fish parameters given the previous ones. The noise vector
\({\xi}=[{\Delta}{\vec{{\alpha}}}{\Delta}{\vec{T}}{\eta}_{2}]^{T}\)
is included in the motion model to account for the uncertainty in the fish parameters, and is drawn from a zero mean multivariate normal distribution with fixed diagonal covariance:
\[\ p_{k}=f(p_{k-1},{\xi}_{k-1})\left[\begin{array}{c}{\vec{{\alpha}}}_{k}\\T_{x_{k}}\\T_{y_{k}}\\{\eta}_{1_{k}}\end{array}\right]=\left[\begin{array}{c}{\Lambda}(u)^{-1}{\Lambda}(u+{\eta}_{0}){\vec{{\alpha}}}_{k-1}\\T_{x_{k-1}}+{{\int}_{0}^{{\eta}_{0}}}\mathrm{cos}({\Phi}({\hat{u}}){\vec{{\alpha}}}_{k-1})\mathrm{d}{\hat{u}}\\T_{y_{k-1}}+{{\int}_{0}^{{\eta}_{0}}}\mathrm{sin}({\Phi}({\hat{u}}){\vec{{\alpha}}}_{k-1})\mathrm{d}{\hat{u}}\\{\eta}_{1_{k-1}}\end{array}\right]+\left[\begin{array}{c}{\Delta}{\alpha}_{k-1}\\{\Delta}T_{x_{k-1}}\\{\Delta}T_{y_{k-1}}\\{\eta}_{2_{k-1}}{\Delta}t\end{array}\right].\]
(A3)
Here, Φ(u) is a N
\({\Theta}\)
dimensional row vector of B-spline bases and
\[{\Lambda}(u)=\left[\begin{array}{c}{\Phi}(u_{1})\\{\cdot}\\{\cdot}\\{\cdot}\\{\Phi}(u_{s})\end{array}\right]\]
is a s by N
\(N_{{\Theta}}\)
matrix of B-spline bases evaluated at s sampled grid points in u (we set s=30). The predicted local bend amplitudes are calculated by projecting onto this basis in the domain of displacement u+η0, calculating the new
\({\Theta}\)
(u+η0), and projecting back onto the original basis, Φ(u). Likewise, the predicted translation is calculated by integrating the tangent vector over the axial displacement.
The measurement model minimizes the normal displacement between the model boundary points q and the edge feature points r detected in the current image. We use the SPKF to minimize the error E=½||nT(q–r)||2(i.e. the sum of squares normal distance), and thus we obtain an updated estimate of our model location p̂k by iterating the equation:
\[\ P_{i+1}={\hat{p}}_{\overline{k}}+P_{pz}P_{zz}^{-1}[z-h(p_{i})-P_{pz}^{T}P_{pp}^{-1}({\hat{p}}_{\overline{k}}-p_{i})],\]
(A4)
where Pab=E[(â–a)(b̂–b)T]is the covariance matrix associated with the random variables a and b and
\({\hat{p}}_{\overline{k}}\)
is the predicted state estimate from our motion model. At the first iteration, Pi+1=
\({\hat{p}}_{\overline{k}}\)
and after convergence (typically 5 iterations for our system), the updated state estimate is set to the current value, p̂k=pi+1.[For more detail on how to calculate the appropriate covariance matrices using the statistical linearization approach, see elsewhere(van der Merwe and Wan, 2003; Sibley et al., 2006).]
The initialization process used to determine the initial shape of the fish is described by the following pseudo-code:
\[\mathrm{while}{\parallel}B_{\mathrm{L}}^{j}-B_{\mathrm{R}}^{j}{\parallel}{>}{\varepsilon}{\ }C_{\mathrm{L}}(i)=\left(B_{\mathrm{L}}^{j}(i)+{\lambda}_{B_{\mathrm{R}}^{j}}(B_{\mathrm{L}}^{j}(i))\right){/}2{\ }C_{\mathrm{R}}(i)=\left(B_{\mathrm{R}}^{i}(i)+{\lambda}_{B_{\mathrm{L}}^{j}}(B_{\mathrm{R}}^{j}(i))\right){/}2{\ }B_{\mathrm{L}}^{j+1}=C_{\mathrm{L}}{\ }B_{\mathrm{R}}^{j+1}=C_{\mathrm{R}}{\ }j=j+1{\ }\mathrm{end}\]
We first calculate the initial left and right boundaries
\(B_{\mathrm{L}}^{0}\)
and
\(B_{\mathrm{R}}^{0}\)
in a semi-automated fashion, then create a correspondence between them using the nearest neighbor functionλ, where λA(B) is the element in A that is closest to B. Then, the computed centerlines become the new boundaries (i.e.
\(B_{\mathrm{L}}^{j+1}=C_{\mathrm{L}},B_{\mathrm{R}}^{j+1}=C_{\mathrm{R}}\)
⁠),and the process is repeated until the average RMS error between the boundaries is less than a threshold (i.e. they have converged on top of each other). Once we have found the centerline through this iterative process, we must determine the width profile of the fish R(u)=ΛR(u)S,where S is a 20×1 vector of control points, andΛ R is the matrix of B-spline bases associated with the radius function. The optimal radius function is found by minimizing the squared normal displacement between the extracted image boundary BL,R and the model boundary ML,Rdictated by the radius function while constraining the width profile to positive values:
\[\ \begin{array}{c}\\\mathrm{argmin}\\s\end{array}{{\sum}_{i}}\left(N_{i}^{T}\left(M_{\mathrm{L}}(u_{i})-{\lambda}_{B_{\mathrm{L}}}\left(M_{\mathrm{L}}(u_{i})\right)\right)\right)^{2}+{{\sum}_{i}}\left(-N_{i}^{T}\left(M_{\mathrm{R}}(u_{i})-{\lambda}_{B_{\mathrm{R}}}\left(M_{\mathrm{R}}(u_{i})\right)\right)\right)^{2},\]
(A5)
\[\mathrm{subject\ to}{\ }R(u){\geq}0\]
where ML,R(ui)=C(ui)±N(ui)ΛR(ui)Sand Ni is the vector normal to C(ui) from the local Frenet frame. The result of this minimization is shown in Fig. 6B and shows that we are able to reconstruct the shape accurately.
We use the spatially smoothed centerlines of the fish to calculate the curvature values directly. A planar B-spline curve has the formγ(u)=Λ(u)X where Λ(u)denotes the s×m matrix of B-spline bases evaluated at s sampled grid points in u, and X denotes the m×2 matrix of control points. Let t(u) and n(u) denote the unit tangent and normal vectors to the curveγ(u), then t′(u)=κ(u).n(u)from the geometry of planar curves. Therefore, the curvature,κ(u), can be derived directly from the B-spline bases as follows:
\[\ {\gamma}^{{^\prime}}(u)={\Lambda}^{{^\prime}}(u)X,\frac{{\partial}}{{\partial}u}\left({\parallel}{\gamma}^{{^\prime}}(u){\parallel}t(u)\right)=\frac{{\partial}}{{\partial}u}({\Lambda}^{{^\prime}}(u)X),{\parallel}{\gamma}^{{^\prime}}(u){\parallel}{\kappa}(u)n(u)={\Lambda}^{{^{\prime\prime}}}(u)X,{\kappa}(u)=\frac{{\langle}{\Lambda}^{{^{\prime\prime}}}(u)X,n(u){\rangle}}{{\parallel}{\gamma}^{{^\prime}}(u)}.\]
(A6)
The tail beat frequency of the fish is calculated as follows. Let F[.] denote the Fourier transform of a function,and κ(u,t) the curvature as a function of position along body axis and time. Then Ki(ω)=|F[κ(u,t)|u=ui]|is the magnitude of the Fourier transform at the ith location along the fish's body, and
\({\omega}_{\mathrm{max}}^{j}=\begin{array}{c}\\\mathrm{argmax}\\{\omega}\end{array}\mathrm{K}_{i}({\omega})\)
is the frequency that maximizes the magnitude response at the ith location. The tail beat frequency f is calculated as:
\[\ f=\frac{{{\sum}_{i=1}^{N}}K_{i}({\omega}_{\mathrm{max}}^{i}){\cdot}{\omega}_{\mathrm{max}}^{i}}{{{\sum}_{i=1}^{N}}K_{i}({\omega}_{\mathrm{max}}^{i})},\]
(A7)
using N locations sampled along the fish's body. This is an average frequency over the length of the body that is weighted by the magnitude response.

LIST OF ABBREVIATIONS AND SYMBOLS

     
  • BL,R

    initial left and right boundaries of fish image (in pixels)

    •  
  • C(i)

    initial centerline of fish image (in pixels)

    •  
  • C(u)

    initial centerline of fish image modeled as B-spline

    •  
  • COA

    center of area

    •  
  • COM

    center of mass

    •  
  • COV

    center of volume

    •  
  • \({\vec{e}}_{1}\)

    tangent vector to fish model centerline

    •  
  • \({\vec{e}}_{2}\)

    normal vector to fish model centreline

    •  
  • f(.,.)

    motion model of zebrafish

    •  
  • h(.,.)

    measurement model of zebrafish

    •  
  • H(p)

    complete geometric model of zebrafish

    •  
  • L

    total length

    •  
  • ni

    outward normal vectors of model boundary

    •  
  • \(N_{{\Theta}}\)

    number of B-spline basis functions used to model centerline bend angle

    •  
  • p0

    initial fish parameters

    •  
  • pk

    fish parameters at the kth time step

    •  
  • Pab

    covariance matrix associated with the random variables a and b

    •  
  • qi

    model boundary point

    •  
  • ri

    edge feature point in image

    •  
  • R(u)

    width profile of zebrafish model

    •  
  • \({\vec{T}}\)

    global translation vector of fish model

    •  
  • (u,v)

    grid parameterization of model mesh

    •  
  • zk

    image measurements at the kth time step

    •  
  • αj

    local bend amplitudes

    •  
  • \({\vec{{\alpha}}}\)

    bend angle amplitude

    •  
  • β

    model scaling term (pixels mm–1)

    •  
  • γ

    ratio of stiff region of zebrafish model to entire length

    •  
  • Δt

    time between frames

    •  
  • η0

    axial displacement of fish between frames

    •  
  • η1

    axial velocity of fish between frames

    •  
  • η2

    axial acceleration fish between frames, sampled from normal distribution

    •  
  • \({\Theta}\)
    (u)

    bend angle of centerline in zebrafish model

    •  
  • λA(B)

    nearest-neighbor function; the element in A that is closest to B

    •  
  • Λ(u)

    s by

    \(N_{{\Theta}}\)
    matrix of B-spline basis functions evaluated at s sampled grid points in u

    •  
  • \({\Phi}_{j}^{k}\)
    (u)

    B-spline basis function of order k

The authors would like to thank B. Walderich, H. M. Maischein and J. Odenthal of the Hubrecht laboratory for graciously allowing us to use video of the stocksteif mutant for demonstrating the ability of our tracker system. We also thank Ansa Wasim for recording the movie sequences and Henk Schipper for capturing the photographs used in the COV calculations. E.I.F.,A.H.B. and J.W.B. thank the Beckman Institute of Caltech for providing financial support of this project. D.L. and J.L.v.L. are funded by NOW-ALW grant 817.02.012.

Bang, P. I., Yelick, P. C., Malicki, J. J. and Sewell, W. F.(
2002
). High-throughput behavioral screening method for detecting auditory response defects in zebrafish.
J. Neurosci. Methods
118
,
177
-187.
Batty, R. (
1984
). Development of swimming movements and musculature of larval herring (Clupea harengus).
J. Exp. Biol.
110
,
217
-229.
Blake, A. and Isard, M. (
1998
).
Active Contours
. London: Springer.
Blaser, R. and Gerlai, R. (
2006
). Behavioral phenotyping in zebrafish: comparison of three behavioral quantification methods.
Behav. Res. Methods
38
,
456
-469.
Budick, S. and O'Malley, D. (
2000
). Locomotor repertoire of the larval zebrafish: swimming, turning and prey capture.
J. Exp. Biol.
203
,
2565
-2579.
Burgess, H. A. and Granato, M. (
2007
). Modulation of locomotor activity in larval zebrafish during light adaptation.
J. Exp. Biol.
210
,
2526
-2539.
Cheng, J. and Pedley, T. (
1998
). A continuous dynamic beam model for swimming fish.
Philos. Trans. R. Soc. Lond. B Biol. Sci.
353
,
981
-997.
Cronin, C. J., Mendel, J. E., Muhktar, S., Kim, Y. M., Stirbl,R. C., Bruck, J. and Sternberg, P. W. (
2005
). An automated system for measuring parameters of nematode sinusoidal movement.
BMC Genet.
6
,
5
.
Domenici, P. and Blake, R. (
1997
). The kinematics and performance of fish fast-start swimming.
J. Exp. Biol.
200
,
1165
-1178.
Doucet, A., de Freitas, N. and Gordon, N.(
2001
).
Sequential Monte-Carlo Methods in Practice
. New York: Springer-Verlag.
Fontaine, E., Burdick, J. W. and Barr, A.(
2006
). Automated tracking of multiple C. elegans. In
Engineering in Medicine and Biology Society, 28th Annual International Conference of the IEEE
, pp.
3716
-3719. doi:
https://doi.org/10.1109/IEMBS.2006.260657
.
Fry, S., Sayaman, R. and Dickinson, M. (
2003
). The aerodynamics of free-flight maneuvers in Drosophila.
Science
300
,
495
-498.
Fuiman, L. and Webb, P. (
1988
). Ontogeny of routine swimming activity and performance in Zebra danios (Teleostei,Cyprinidae).
Anim. Behav.
36
,
250
-261.
Full, R. and Tu, M. (
1991
). Mechanics of rapid running insect-2-legged, 4-legged and 6-legged locomotion.
J. Exp. Biol.
156
,
215
-231.
Geng, W., Cosman, P., Berry, C. C., Feng, Z. and Schafer, W. R. (
2004
). Automatic tracking, feature extraction and classification of C. elegans phenotypes.
IEEE Trans. Biomed. Eng.
51
,
1811
-1820.
Gray, J. (
1933
). Studies in animal locomotion. I. The movement of fish with special reference to the eel.
J. Exp. Biol.
10
,
88
-104.
Hamarneh, G. and McInerney, T. (
2001
). Controlled shape deformations via medial profiles. In
Vision Interface
, pp.
252
-258.
Harper, D. G. and Blake, R. W. (
1989
). On the error involved in high-speed film when used to evaluate maximum accelerations of fish.
Can. J. Zool.
67
,
1929
-1936.
Hedrick, T., Usherwood, J. and Biewener, A.(
2004
). Wing inertia and whole-body acceleration: an analysis of instantaneous aerodynamic force production in cockatiels (Nymphicus hollandicus).
J. Exp. Biol.
207
,
1689
-1702.
Ito, K. and Xiong, K. (
2000
). Gaussian filters for nonlinear filtering problems.
IEEE Trans. Automat. Contr.
45
,
910
-927.
McElligott, M. B. and O'Malley, D. M. (
2005
). Prey tracking by larval zebrafish: axial kinematics and visual control.
Brain Behav. Evol.
66
,
177
-196.
McHenry, M. (
2001
). Mechanisms of helical swimming: asymmetries in the morphology, movement and mechanics of larvae of the ascidian Distaplia occidentalis.
J. Exp. Biol.
204
,
2959
-2973.
McHenry, M. J., Pell, C. A. and Long, J. H., Jr(
1995
). Mechanical control of swimming speed: stiffness and axial wave form in undulating fish models.
J. Exp. Biol.
198
,
2293
-2305.
Müller, U. K. and van Leeuwen, J. L.(
2004
). Swimming of larval zebrafish: ontogeny of body waves and implications for locomotory development.
J. Exp. Biol.
207
,
853
-868.
Muybridge, E. (
1887
).
Animal Locomotion
. Philadelphia: University of Pennsylvania.
Nörgaard, M., Poulsen, N. K. and Ravn, O.(
2000
). New developments in state estimation for nonlinear systems.
Automatica
36
,
1627
-1638.
Park, K. J., Rosen, M. and Hedenstrom, A.(
2001
). Flight kinematics of the barn swallow (Hirundo rustica) over a wide range of speeds in a wind tunnel.
J. Exp. Biol.
204
,
2741
-2750.
Pfau, T., Witte, T. H. and Wilson, A. M.(
2005
). A method for deriving displacement data during cyclical movement using an inertial sensor.
J. Exp. Biol.
208
,
2503
-2514.
Roussel, N., Morton, C., Finger, F. and Roysam, B.(
2007
). A computational model for C. Elegans locomotory behavior: application to multiworm tracking.
IEEE Trans. Biomed. Eng.
54
,
1786
-1797.
Sibley, G., Sukhatme, G. and Matthies, L.(
2006
). The iterated sigma point Kalman filter with applications to long range stereo.
Robotics Science and Systems
, http://www.roboticsproceedings.org/rss02/p34.pdf.
Standen, E. and Lauder, G. (
2005
). Dorsal and anal fin function in bluegill sunfish Lepomis macrochirus:three-dimensional kinematics during propulsion and maneuvering.
J. Exp. Biol.
208
,
2753
-2763.
Thorsen, D. H., Cassidy, J. J. and Hale, M. E.(
2004
). Swimming of larval zebrafish: fin-axis coordination and implications for function and neural control.
J. Exp. Biol.
207
,
4175
-4183.
Tytell, E. D. and Lauder, G. V. (
2002
). The C-start escape response of Polypterus senegalus: bilateral muscle activity and variation during stage 1 and 2.
J. Exp. Biol.
205
,
2591
-2603.
van der Merwe, R. and Wan, E. (
2003
). Sigma-point Kalman filters for probabilistic inference in dynamic state-space models.
Proceedings of the Workshop on Advances in Machine Learning
, Montreal, Canada.
Wainwright, S. (
1983
). To bend a fish. In
Fish Biomechanics
(ed. P. Webb and D. Weihs), pp.
68
-91. New York: Praeger.
Walker, J. A. (
1998
). Estimating velocities and accelerations of animal locomotion: a simulation experiment comparing numerical differentiation algorithms.
J. Exp. Biol.
201
,
981
-995.
Weihs, D. (
1973
). The mechanism of rapid starting of slender fish.
Biorheology
10
,
343
-350.
Willmott, A. and Ellington, C. (
1997
). Measuring the angle of attack of beating insect wings: robust three-dimensional reconstruction from two-dimensional images.
J. Exp. Biol.
200
,
2693
-2704.
© The Company of Biologists Limited 2008
2008

Supplementary information

Supplemental Figure S1 (Adobe PDF)

Fig. S1. Example of center of volume (COV) and center of area (COA) calculation for zebrafish at age 5 days. Volume assumes an elliptical cross-section of the fish. COV will coincide with the COM is a uniform density is assumed.

- pdf file
Supplemental Figure S2 (Adobe PDF)

Fig. S2. Example of center of volume (COV) and center of area (COA) calculation for zebrafish at age 15 days. Volume assumes an elliptical cross-section of the fish. COV will coincide with the COM is a uniform density is assumed.

- pdf file
Supplemental Figure S3 (Adobe PDF)

Fig. S3. Example of center of volume (COV) and center of area (COA) calculation for zebrafish at age 28 days. Volume assumes an elliptical cross-section of the fish. COV will coincide with the COM is a uniform density is assumed.

- pdf file
Supplemental Figure S4 (Adobe PDF)

Fig. S4. Location of center of volume (COV) and center of area (COA) measured posterior from the snout of the fish. For ages 5, 15 and 28 days, we measured data from 3, 5 and 4 fish, respectively. Horizontal bars indicate the mean value. As the fish get older, the COV and COA become closer to the same location.

- pdf file
Movie 1

Movie 1. 15 d.p.f. wild-type zebrafish. Tracking results for 15 d.p.f. wild-type zebrafish. Note that accurate tracking is achieved despite the partial occlusion by the horse hair.

- mov file