Automated Reconstruction of Three-Dimensional Fish Motion, Forces, and Torques

Fish can move freely through the water column and make complex three-dimensional motions to explore their environment, escape or feed. Nevertheless, the majority of swimming studies is currently limited to two-dimensional analyses. Accurate experimental quantification of changes in body shape, position and orientation (swimming kinematics) in three dimensions is therefore essential to advance biomechanical research of fish swimming. Here, we present a validated method that automatically tracks a swimming fish in three dimensions from multi-camera high-speed video. We use an optimisation procedure to fit a parameterised, morphology-based fish model to each set of video images. This results in a time sequence of position, orientation and body curvature. We post-process this data to derive additional kinematic parameters (e.g. velocities, accelerations) and propose an inverse-dynamics method to compute the resultant forces and torques during swimming. The presented method for quantifying 3D fish motion paves the way for future analyses of swimming biomechanics.

with rŷ and rẑ the radii in respectivelyŷ-andẑ-direction and ψ = [0, 2π] where the sign in the equation forẑ depends on the considered quadrant: + in the positiveẑ-quadrant, − in the negativeẑ-quadrant. For the fin fold and body cross-sections in the tail, the exponent n = 2, making it a regular ellipse. The elliptic exponent for the eyes is 2.3. The slices of the body with a bulging yolk sack have their egg-like cross-sectional shape defined by a cubic spline fit (using MATLAB's spline) through 4 points: (x 0 , rŷ,ẑ 0 + h shift ), (x 0 , 0,ẑ 0 + rẑ), (x 0 , −rŷ,ẑ 0 + h shift ), (x 0 , 0,ẑ 0 − rẑ). To prevent edge effects of the spline, we perform the fit on three repetitions of the points, and only use the middle section as a cross-section.
We build up the body with N lon cross-sections and N circ points along the circumference of each cross-section. We can connect each set of two cross-sections with (N lon − 1)N circ quadrilateral faces, or 2(N lon − 1)N circ triangular faces. To close the surface, we cap the front-and aftmost cross-sections with N circ triangular faces each.

Parameterisation
The body curvature κ is prescribed in a small number (< N lon , 7 for the zebrafish larva) of control points along the length of the fish. We define a normalised parameter along the fish length s = [0, 1] and use it to create a cubic spline (MATLAB's spline) that interpolates the curvature control points. We analytically integrate each piecewise polynomial of the curvature spline to obtain the local deformation angle: 1/8 with the body length. We model the fish's centreline as a sequence of straight segments i of length ∆ = /(N lon − 1), and apply the local deformation angle to each: in coordinates relative to the fish's head. Every cross-section i consists of N circ circumferential points to which we assign the index j. We rotate every point i, j on the cross-section with the rotation of the centreline, such that it remains perpendicular: We now have a deformed fish in the (x,ŷ,ẑ) coordinate system, which we place in world coordinates with the position x snout and rotation ϕ at the snout. The angles ϕ roll , ϕ pitch , ϕ yaw (abbreviated as ϕ r , ϕ p , ϕ y for reasons of clarity) specifying the rotation are converted to a rotation matrix: cos ϕ p cos ϕ y − cos ϕ p sin ϕ y sin ϕ p sin ϕ r sin ϕ p cos ϕ y + cos ϕ r sin ϕ y − sin ϕ r sin ϕ p sin ϕ y + cos ϕ r cos ϕ y − sin ϕ r cos ϕ p − cos ϕ r sin ϕ p cos ϕ y + sin ϕ r sin ϕ y cos ϕ r sin ϕ p sin ϕ y + sin ϕ r cos ϕ y cos ϕ r cos ϕ p We compute the final position of each point on the fish with: Finally, all points are reassembled into a three-dimensional surface of quadrilaterals during tracking and of triangles during post-processing.

Fitting the body model
Objective function The objective function that we aim to minimise in order to achieve an optimal fit is defined as follows: Here, Ω is the list of optimised parameters, f GoF (Ω) is a term indicating how accurately the fish model approaches the high-speed video images, and f reg (Ω) is a term that penalises "unsmoothness". The two terms are discussed in detail below.

Goodness of fit
To compare the virtual images to the segmented high-speed video images, we project the three-dimensional shape of the fish onto virtual cameras. We define a left-handed camera coordinate system (ξ, η, ϑ), respectively rows, columns, and a coordinate perpendicular to the image plane in the viewing direction of the camera. A transformation matrix from camera to world coordinate system is given by: Since we are considering a camera system with parallel light, the distance of an object from the camera has no influence on its projected size. The light is collimated to within 10 • divergence, so perspective effects are negligible at the millimetre scales of zebrafish motion. Projecting objects onto the camera is therefore only a coordinate transformation to camera coordinates (i.e., a rotation, translation and rescaling): where x object are the object coordinates, x cam the camera coordinate system origin (i.e., the top left corner of the image), d scale a scaling factor in unit distance per pixel. The 1 is added because by convention, the top left corner is denoted as having row-column (i.e., ξ, η) coordinates (1, 1). Because the transformation matrix T cam is orthogonal, we can invert the matrix by taking its transpose. The image plane coordinates of the projected point are given by (ξ, η); the third coordinate ϑ can be ignored. When considering cases with significant magnification effects (e.g., when no collimated light setup is used), we use a simple linear pinhole model, where image magnification is linear with distance from the camera. In that case, we use the ϑ coordinate to determine the scaling factor for the current projected point and thus compute its projection on the image plane.
To compute the goodness of fit, we project the three-dimensional surface of the fish onto the camera plane using the procedure outlined above. We then overlay the fish silhouette in image coordinates onto the segmented high-speed video frames. All pixels in the combined image are counted, and the pixels that overlap between the model fish and the actual fish are subtracted from this number. The final result is the number of mismatching pixels -an expression for the goodness of fit.

Penalty function
The regularising term in the objective function is of the form where w(s) is a weighting function that allows stronger suppression of curvature near the head than near the tail. For the zebrafish, we prescribe w(s) as a cubic spline fit (MATLAB's spline) in logarithmic space of control points along the body. In our case of the larval zebrafish, we prescribe two points at s = {0, 1}. Since a cubic spline fit of two points is linear, the result is a function of the form w(s) = c 0 e −c1s .

Initialisation
The initial set of model parameters for the first frame ω 0 is computed from two user-clicked points per camera: the snout and the tail. Rewriting Eqn. 9 gives, for every camera, Manually indicating a point in an image will give us the first two components of ξ, so we can use the first two equations of this system for each camera. This results in an overdetermined linear system (for N cam > 2): we have 3 unknowns (x, y, z) and 2N cam equations. We solve this system with least-squares (MATLAB's \ operator) to obtain the snout and tail positions in 3D. The snout position can be directly used to initialise the position parameter for the optimisation. We set the initial roll angle to 0, the pitch angle to and the yaw angle to The arctangents are computed using MATLAB's atan2 function, that takes into account the sign of the nominator and denominator to compute the angle in the correct quadrant.
We approximate the fish length as follows: This approximation is used as an initial condition for an optimisation in the first frame of the sequence, where we allow the fish length to vary. The final, optimised fish length is then kept fixed in the rest of the video sequence.

Solution extrapolation
We extrapolate the solution from previous time steps to initialise the next time step. We use the following: where i t ∈ N denotes the time step, Ω init,it is the set of parameters used to initialise the optimisation at time step i t and Ω init,it is the final, optimised set of parameters at time step i t .

Post-processing and inverse dynamics Smoothing
In general, the solution from the optimisation procedure is insufficiently smooth in time to accurately compute (second) derivatives. We perform post hoc smoothing to remove spurious high-frequency components, but retain (most of) the physical signal. We use a regularised data fit similar to [1], where smoothing is presented as a minimisation problem: where φ sm are the smoothed data points, φ data are the measured data points, λ sm is a parameter controlling the amount of regularisation and n is the order of the derivative used for the penalty function. It has been shown [1] that the solution for a discrete number of points is given by where D is a matrix that approximates the required n th -order derivative. We assume here that our fitted points are located at the same time instants as the data points. For the D matrix, we use 4 subsequent one-sided differences (using MATLAB's diff) to approximate the fourth order derivative, ensuring that our second derivatives are smooth.
To combat edge effects of the smoothing, we separate two cases: starts, where the animal accelerates from a stationary position, and "continuous" swimming, where it enters and exits the field of view while swimming. In the case of starts, we know that the animal is stationary before the time series starts. We therefore add 25 time points before t = 0 which all have the value at t = 0, thus forcing the solution to have a realistic zero gradient at the start. In the case of "continuous" swimming, we can simply cut off the first few (in our case 5) points, where the edge effects are strongest. For both cases, we cut off the same number of points off the end of the time series.

Derivatives
Differentiation is required to compute velocity and acceleration from the centre of mass position, and torque from the angular momentum vector. We use the same procedure for any time series, for first derivatives: for second derivatives:

Resultant forces and torques
The resultant force F net on the fish can be computed from Newton's second law: where m is fish mass, ρ is fish density, and a CoM denotes the acceleration of the centre of mass; the volume V and position of the centre of mass r CoM can be computed using the method by [2] from the triangulated surface of the body model. Every triangular face is made into a tetrahedron by taking the coordinate system origin as a fourth vertex, see Fig 1A. Calculations can be performed separately on each of these tetrahedra, and subsequently summed to yield the result for the complete fish. The method by [2] only computes the centre of mass and the moment of inertia matrix; we extend their method to compute angular momentum for an arbitrary triangulated polyhedron of constant density. Angular momentum is defined as: where the asterisks denotes quantities measured with respect to the centre of mass. The integrand is linear, so we can evaluate Eqn. 21 for each tetrahedron separately and then sum the results. To simplify this computation, we transform every tetrahedron into a unit tetrahedron in the origin of a new coordinate system (x , y , z ), see Fig 1B. This is a linear transformation, with the following transformation matrix from unit to world tetrahedron: The transformed integral for every tetrahedron becomes Transformed to a unit tetrahedron, linear interpolation of quantities known at its vertices reduces to a single matrix multiplication. The interpolating matrices for the position and velocity are given by: With interpolation, Eqn. 23 can be written as: Evaluating this integral (using Maple 18, Maplesoft, Waterloo, Ontario, Canada) gives the following expressions for the x-, y-and z-components of the angular momentum vector: The total angular momentum is then given by the sum over all tetrahedra k: With the angular momentum known in each time step, we can compute the net torque by taking its time derivative: τ net = dL dt . (31)

Fish reference system
To make interpretation of the forces and torques more intuitive, we transform them to a local fish coordinate system (x fish , y fish , z fish ) that rotates with the changing orientation of the fish. To create this fish axis system, we first rotate the world coordinate system so it is aligned with the fish's head, using the rotation matrix in Eqn. 5. This transforms the fish to the (x,ŷ,ẑ) coordinate system defined in section Creating the parameterised fish model, where thex,ŷ-plane is aligned with the deformation plane of the fish, andẑ is perpendicular to it. Thex-coordinate is aligned with the head, which rotates with a large amplitude and makes interpretation of the time-variation of the forces and torques difficult. In order to define an angle that is related to the entire body, we use a "body angle" [3], defined as: where I CoM,i (s) is the moment of inertia about the centre of mass per unit length of the cross-section at s and α(s) is the angle between a line projected from the centre of mass to the centre of segment and horizontal. The body angle is effectively a moment of inertia-weighted average of the local angles of each segment. We can express the local moment of inertia per unit length as the contribution of the cross-section about its own vertical centre axis I cross,i and its contribution with respect to the centre of mass where µ is a mass per unit length and r 2 CoM,i the distance between the centre of mass in the x fish , y fish -plane and the cross-section centre. We compute the moment of inertia per unit length and area, and thus mass, assuming constant density, of each cross-section with the expressions for arbitrary polygons [4]. The local segment angle α is computed by projecting a line from the centre of mass to the segment centre in the deformation plane (i.e., inx,ŷ,ẑ-coordinates), and computing the angle it makes with respect to the head: making sure we get the correct sign depending on the quadrant we are in by using MATLAB's atan2 function.
We compute the body angle by evaluating the integral from Eqn. 32 using MATLAB's trapz function. This results in a series of body angles over time that may still have discontinuities where the centre of mass crosses the centreline; in these points, the angle "reflects" around an unknown line since all the projected lines flip direction. To remove these discontinuities, we average a one-sided derivative left and right of the discontinuity. We then extrapolate the solution from the point left of the discontinuity: where i is the index of the point just right of the discontinuity. The entire sequence to the right of the discontinuity is shifted, retaining gradient information. We iterate over all discontinuities in the sequence, yielding a body angle that is continuous from the first angle onwards. Since we start with a straight fish in thex,ŷ,ẑ-coordinates, we define this first angle as 0.
We rotate the head-aligned coordinate system along its vertical axis by the computed body angle. We represent the rotation as an axis and angle: respectively ω z head (abbreviated in Eqn. 36 as ω), the unit orientation vector of the z head -axis and ϕ body (abbreviated in Eqn. 36 as ϕ).This results in a rotation matrix of the following form: x (1 − cos ϕ) ω x ω y (1 − cos ϕ) −ω z sin ϕ ω x ω z (1 − cos ϕ) +ω y sin ϕ ω x ω y (1 − cos ϕ) +ω z sin ϕ cos ϕ +ω 2 y (1 − cos ϕ) ω y ω z (1 − cos ϕ) −ω x sin ϕ ω x ω z (1 − cos ϕ) −ω y sin ϕ ω y ω z (1 − cos ϕ) +ω x sin ϕ cos ϕ +ω 2 z (1 − cos ϕ) 7/8 Applying this rotation matrix to the axes of the head coordinate system yields our final fish coordinate system (x fish , y fish , z fish ), in which we express our forces and torques. We transform the force vector by projecting it onto each unit orientation vector of the axes, e.g., We define an additional axis in the direction of the velocity vector: where v CoM is the instantaneous velocity vector of the centre of mass. We project the force to this axis to have an indication of the effective force propelling the fish.