Continuous body 3-D reconstruction of limbless animals

Limbless animals such as snakes (Byrnes and Jayne, 2012; Gans, 1986; Goldman and Hu, 2010; Gray and Lissmann, 1950; Jayne, 1985, 1986; Lillywhite et al., 2000; Marvi et al., 2014; Munk, 2008; Socha, 2002)), limbless lizards (Gasc and Gans, 1990; Gans et al., 1992; Miller, 1944), worms (Karbowski et al., 2006; Park et al., 2008; Summers and O'Reilly, 1997; Dorgan, 2015) and fish (Gillis, 1998; Jayne and Lauder, 1995; Tytell, 2007; Herrel et al., 2011; Gidmark et al., 2011; Gemmell et al., 2015) can deform their long, slender bodies to move through a large diversity of environments. To move through complex environments such as branches (Jorgensen and Jayne, 2017; Byrnes and Jayne, 2012), underwater sand beds (Gidmark et al., 2011; Tatom-Naecker and Westneat, 2018), large obstacles (Gart et al., 2019), and even during gliding (Socha, 2011) and swimming (Graham and Lowell, 1987), their body must deform substantially (>10% body length) in three dimensions. In addition, how the body is oriented and rotates locally relative to the environment often strongly affects the forces that the animal generates. For example, snakes, limbless lizards and worms all have bodies with anisotropic friction properties (Hu et al., 2009; Spinner et al., 2015; Schulke et al., 2011; Shen et al., 2012). During borrowing in granular media, granular resistive stress depends sensitively on local orientation of a moving body (Li et al., 2013). In addition, pressure drag and skin friction in fluids depend on body shape and orientation (Tytell, 2007; Gemmell et al., 2015; Taylor, 1952; Jayne and Lauder, 1995; Holden et al., 2014; Danos and Lauder, 2012). Furthermore, when traversing large obstacles, where the body contacts the terrain affects a limbless animal's stability (Gart et al., 2019). Therefore, to understand limbless locomotion in these environments, it is important to accurately quantify an animal's body shape and movement in terms of both position and orientation in three dimensions.

Many previous studies have quantified limbless animal body shape and movement in three dimensions using a series of tracked points (Kwon et al., 2013; Byrnes and Jayne, 2012; Marvi et al., 2014; Socha, 2005). However, discrete points cannot accurately capture curved local body shape necessary for quantifying body–environment interactions (e.g. Movie 1) unless a very large number of markers are tracked, which is time-consuming and challenging when marker occlusion is frequent (e.g. large 3-D rotation, multiple obstacles present). Some studies used superposition of curves generated by basis functions (e.g. B-spline) (Sharpe et al., 2015; Schiebel et al., 2018; Fontaine et al., 2008; Padmanabhan et al., 2012; Yeaton et al., 2020), but the majority of these superposition methods were for planar locomotion. In addition, none of these geometric interpolation methods above capture body rotation (roll) about the longitudinal axis [although some studies measured it using computer vision techniques (Donatelli et al., 2017) or manually (Fish et al., 2007)].

A simplistic way of quantifying a long, slender body, with both position and orientation information, is to approximate it as an elastic rod. Elastic rod theories developed by Kirchhoff and Cosserat view an elastic rod as continuum on which forces are applied (for reviews, see Cao and Tucker, 2008; Dill, 1992; Zhang et al., 2019). Despite being biomechanically over-simplified – not considering musculoskeletal morphology and muscle function (Dickinson et al., 2000; but see Zhang et al., 2019) – elastic rod modeling has facilitated first-order modeling and basic understanding of the mechanics (e.g. muscle activation, resistive forces from the surrounding media) of limbless animals moving both in two (Long, 1998; Tytell et al., 2018; Ding et al., 2013; Cheng et al., 1998; Zhang et al., 2019) and three (McMillen and Holmes, 2006) dimensions.

In the field of robot motion planning, elastic rod modeling has been combined with backbone optimization to generate motion trajectories of snake-like, high degree-of-freedom (hyper-redundant) robots to achieve the desired locomotor tasks (Burdick et al., 1994; Chirikjian, 2015). The idea of backbone optimization is that, when the position and orientation of both ends of a high degree-of-freedom system are given (hereafter referred to as end constraints), the approximate shape (referred to as the backbone curve, Fig. 1A inset, dashed curve) of its actual midline in between (referred to as the midline) can be determined by minimizing a cost function. By representing a long, slender body as an elastic rod with its elastic potential energy as the cost function and applying backbone optimization over time, one can determine the evolution of body shape, or motion trajectory, that satisfies the desired end constraints (e.g. moving from point A to point B around an obstacle). Similarly, this method has enabled shape interpolation of stiff molecules such as DNA subject to end constraints (Kim and Chirikjian, 2006). It is important to note that, in backbone optimization, the use of an elastic rod is not intended to capture the physical reality of the system, but instead to simplify mathematical calculations (see Materials and Methods, Overview).

Here, we propose to use this method combining elastic rod modeling and backbone optimization to reconstruct the continuous body of limbless animals in three dimensions, using end constraints from discrete tracked markers. We emphasize that the purpose of our method here is to geometrically quantify the animal's continuous body shape and movement kinematics, rather than model its locomotion mechanics and dynamics (for doing this instead, see Boyer et al., 2010, Zhang et al., 2019). We describe our method using the example of locomotion of the variable kingsnake (Lampropeltis mexicana thayeri; Garman, 1883) with large body deformation and motion in three dimensions. We simplified the animal body as an elastic rod and applied backbone optimization to interpolate between markers placed on the body to obtain 3-D position and orientation of the continuous body. In addition, we characterized how marker separation affects positional accuracy of the interpolation, using a ground truth of the body midline positions extracted via computer vision techniques. Further, we compared our method with the widely used B-spline method (e.g. Sharpe et al., 2015; Schiebel et al., 2018; Fontaine et al., 2008; Yeaton et al., 2020) by applying both to a large dataset of kingsnakes traversing large step obstacles and comparing their position and orientation errors using tracked marker data as reference. We first interpolated the x, y, z positions together by fitting a 3-D B-spline curve to tracked marker locations, then calculated orientation from the interpolated B-spline position curve using the Bloomenthal method (Bloomenthal, 1990). Finally, we discuss why our method achieves higher interpolation accuracy than entirely geometric interpolation methods, and the benefits of such higher accuracy. We provide MATLAB codes and demonstration files for users to easily learn and apply our method to their studies (see Data availability).

The model system to describe our method is the variable kingsnake traversing large step obstacles (Gart et al., 2019). Multiple markers are attached to the animal body from head to tail to capture both 3-D position and orientation locally, and they divide the body into segments (Fig. 1A). Each elastic rod segment approximating a body segment is subject to two end constraints imposed by the 3-D position and orientation of the two markers at both its ends (Fig. 1A, insets). Interpolation is performed piecewise between each adjacent pair of markers and for both 3-D position and orientation, which is simplified using Lie group theory (with the number of equations reduced and singularity of Euler angles avoided). For an introduction of Lie group theory, see Murray et al. (1994). Below, we follow the convention in Kim and Chirikjian (2006).

Our method approximates the continuous shape of each body segment by that of a quasi-static elastic rod dominated by elastic forces, subject to end constraints from tracked markers. This is a drastic over-simplification, because in reality the deformation of an animal body segment is a result of all the forces acting on it, which include gravitational force, inertial force due to body acceleration, forces from the surrounding environment, and internal forces due to muscle activation and viscoelastic skin, tissue and skeleton (Cheng et al., 1998). The purpose of this over-simplification is to simplify derivation of the Lagrangian in Lagrange's equation (such that it is merely the elastic potential energy of the rod), which we used in backbone optimization to interpolate between end constraints. In principle, other forces can be modeled by the Lagrangian (e.g. gravitational potential energy) or as constraints using Lagrangian multipliers (e.g. non-conservative forces). However, as we demonstrate in the Results, even with such a drastic over-simplification, our method achieves higher accuracy in 3-D shape interpolation than the widely used, purely geometric B-spline method.

In our method, the elastic rod can experience twisting, shear, extension and compression aside from bending. It is unclear whether these body deformation types occur in every species – for example, twisting exists in some snakes but is limited by the vertebrae structure and varies between species (Moon, 1999; Jurestovsky et al., 2020). However, inclusion of these deformation types in the model allows more general application to different species. It also allows for better interpolation by accommodating the inevitable measurement errors of end constraints and segment lengths due to camera tracking noise and movement of the animal's skin. The reconstruction results with lateral and dorsoventral bending having the dominant contributions to elastic energy (Fig. S2) suggested that the inclusion is reasonable when applied to snakes.

For a long, slender body with very high degrees of freedom to satisfy a few discrete end constraints, many solutions exist. This is referred to as the kinematically redundant or hyper-redundant problem (Chirikjian and Burdick, 1995). A well-developed method to address this problem is backbone optimization, which minimizes a cost function to determine one or a few solutions (Chirikjian and Burdick, 1995). Such constrained optimization problems are usually solved using the Euler–Lagrange equation, which requires the use of coordinates that can result in singularity. Here, we use a coordinate-free method to avoid singularities, which consists of the Euler–Poincaré equation and a kinematic equation (Chirikjian, 2015), with the elastic energy of the rod segment as the cost function to minimize, subject to end constraints. This results in a solution of spatial velocity (spatial derivative of body frame, as viewed in the instantaneous local body frame; see Eqn 2), which can be integrated from one end constraint towards another to obtain a backbone curve.

In order for the integrated backbone curve to interpolate between the two end constraints measured by adjacent markers, we apply the method of inverse kinematics (Kim and Chirikjian, 2006) (Fig. 1B). For each body segment, we first start from one end attached to one of the measured end constraints and integrate spatial acceleration (spatial derivative of spatial velocity) to obtain a backbone curve. The spatial acceleration is calculated using an initial guess of spatial velocity at the starting end and the Euler–Poincaré equation. We then iteratively perturb the backbone curve (by perturbing the initial guess of spatial velocity) to reduce the distance of its other end to the second measured end constraint until they converge. This entire procedure is performed piecewise for all segments between adjacent markers, and the interpolated backbone curves of all segments together give the approximate continuous midline of the entire body (except the very front and rear end without markers).

To apply our method in animal experiments, the workflow is as follows (Fig. 1C). First, users need to place markers on the animal body that can provide 3-D position and orientation information, such as BEEtags (Crall et al., 2015), ArUco (Garrido-Jurado et al., 2014) and custom rigid body markers (Ravi et al., 2013). The users need to then track the markers to obtain their x–y–z coordinates and Euler angles as input. Then, they need to measure animal length, width and height (and body tapering if it is substantial) and pre-process the tracking data to obtain end constraints. After these preparations, they can run the backbone interpolation codes followed by post-processing to obtain the interpolated animal body midline for further analysis. Convergence check and twist/writhe (how much the rod twists into coils) check are performed automatically before results are saved. To validate the results, the users can project the reconstructed body midline onto videos, visually check the match, and fine-tune interpolation parameters to improve accuracy.

Below we describe our method in detail. Because our interpolation is performed piecewise on the body, the description focuses on one body segment between two adjacent markers approximated by one elastic rod segment, unless specified otherwise.

where R(s,t)=[r₁(s,t) r₂(s,t) r₃(s,t)] is a 3×3 matrix representing the marker orientation, with each column the coordinates of a unit vector along the +x,y,z axis of the body frame in the fixed world frame (Fig. 1D), p(s,t)=[x(s,t) y(s,t) z(s,t)]^T are the x,y,z coordinates of the origin of the body frame, 0^T=[0 0 0], M^T denotes the transpose of a matrix M, s is arc length from one end of the unstretched body segment (e.g. the marker closest to snout or tail tip), and t is time. Below, all lengths refer to unstretched lengths unless stated otherwise.

An elastic rod can be bent, twisted, sheared, extended or compressed under external forces and torques (Fig. 1D). Below, we derive the elastic energy of a rod segment using Lie group notations, which will be used for optimization in the next section (Kim and Chirikjian, 2006).

where K is the 6×6 stiffness matrix, which is symmetric (K=K^T) and assumed to be constant, and is the intrinsic deformation of the body frame along the rod when there is no external force. We further assumed that there is no coupling between forces along different directions and torques about different axes, i.e. K is diagonal.

For an elastic rod segment that satisfies any two end constraints, we can obtain its backbone curve using backbone optimization with elastic energy as the cost function. We omit time t in this and the next sections considering the quasi-static assumption, which works reasonably for friction-dominated snake locomotion (Hu et al., 2009). For more details of the backbone optimization method, see Kim and Chirikjian (2006).

Next, we numerically integrate Eqn 6 to obtain spatial velocity ⁠. Then, starting from one end of the segment (s=0) with a spatial velocity value chosen at this end (boundary condition), one can integrate over arc length s towards the other end (s=L) to obtain a backbone curve g(s) of the entire segment (s∈[0, L], where L is segment length). The resulting backbone curve depends on the boundary condition chosen and does not necessarily satisfy the measured end constraints from tracked markers. Next, we apply the method of inverse kinematics to converge the backbone curve to satisfy the measured end constraints.

In the first iteration (k=1), we start from one end constraint g(0) tracked by a marker (Fig. 1B) and perform the integration above to obtain a backbone curve that ends at g₁(L). Because the markers do not provide information about spatial velocity, in the first step, we use an initial guess of ⁠. As a result, the other end g₁(L) of the resulting backbone curve (Fig. 1B, bottom dashed curve) is far away from the other end constraint from the other tracked marker g(L).

Over the following steps k=2,…,n, we start with used in the previous step and iteratively perturb it to perturb the backbone curve. Using the inverse kinematics method, we gradually reduce the distance from the other end g_k(L) of the resulting backbone curve to the other measured end constraint g(L), until they converge. This results in the final backbone curve (Fig. 1B, top dashed curve) that converges to the backbone curve that exactly satisfies the two end constraints from tracked markers (Fig. 1B, solid curve).

where is the Frobenius norm of a matrix R_k(L)−R(L) that measures error in orientation, and is the Euclidean norm of a vector that measures error in position.

We set the inverse kinematics iteration to stop when d(g_n(L), g(L))<0.1. This results in a position error smaller than 0.1 mm (1% body diameter of 9 mm) and an orientation error smaller than 6 deg. If not, the iteration stops when the maximal iteration step of 65 is reached. For more details of the inverse kinematics iteration, see Kim and Chirikjian (2006).

Sometimes, an initial guess of does not produce a backbone curve that converges to the tracked marker end constraint g(L), or it results in a shape with unrealistic, large twist or writhe. To resolve this, we try different initial guesses. We try the inverse kinematics procedure at most three times for each trial before a solution is found for each body segment in each video frame, using different initial guesses, iteration step size and maximal number of iterations. We also perform a twist and writhe check (Kim and Chirikjian, 2006) using pre-defined thresholds of 0.4 (determined by visually examining body frames interpolated using different thresholds projected onto the videos) to ensure that the interpolated backbone curve is realistic.

To reduce computation time, we always first try the value from the previous successful video frame for the same segment for the following frame, because the shape change between these video frames is usually small with a sufficiently high camera frame rate. If no previous video frame succeeds or this value fails to produce realistic curves within the maximal iteration steps twice, we then try pre-defined initial guesses with extension and bending in different directions.

Several geometric parameters in the method affect the fidelity of the interpolation and should be measured as input, including the length of the body segment between markers L, body width and height, and body tapering.

Ideally, segment lengths should be directly measured piecewise between each pair of adjacent markers after marker placement. Because this was challenging for the conscious animals in our experiments, we approximated segment lengths by the largest distance observed between adjacent markers across all video frames in all trials with the same marker placement. Digital photography analysis is an alternative method (Astley et al., 2017). During reconstruction, we further fine-tuned segment lengths to improve visual match by projecting the reconstructed body onto experimental videos and evaluating whether the estimated segment lengths were too long or too short and adjusting them accordingly.

We measured the width and height of the three kingsnakes at 10 nearly equally spaced locations along the body (Fig. S1). We found that their body width is similar to body height, i.e. the body cross-section is nearly circular. In addition, there is little body tapering. Based on these measurements, for interpolation, we approximated the cross-sectional shape as a circle with a radius r=4.5 mm equal to the average of width and height measurements. For pre-processing of marker data, we calculated the body radius of each individual as a function of body length, r(s), by linear interpolation.

where κ₁ and κ₂ are lateral and dorsoventral bending stiffnesses about the x- and y-axes, κ₃ is torsional stiffness, κ₄ and κ₅ are lateral and dorsoventral shearing stiffnesses along the x- and y-axes, and κ₆ is extensive/compressive stiffness (which are assumed to be the same).

where ν is the Poisson's ratio.

With these assumptions, all elements of K scale with E. As a result, changing E simply scales elastic energy and does not change the solution.

If the cross-section is not circular, such as an ellipse, these moments of inertia need to be calculated accordingly. Our preliminary tests showed that the interpolated midline only changed by 6% body diameter in position when width:height ratio changed from 1 to 0.3, although buckling appeared when it reduced to 0.1.

We used a Young's modulus of the same order of magnitude as measured from eel (Long, 1998) and lamprey (Tytell et al., 2018) tissues, E=10⁵ Nm^–2, and a Poisson's ratio similar to that of human tissue, ν=0.3 (Choi and Zheng, 2005).

In experiments, even when care was taken during marker placement, markers can never be perfectly parallel to the local coronal plane of the animal body with the front edge perfectly perpendicular to the midline. Such an offset can affect interpolation results.

To reduce this offset, for each marker, we observed all the videos with the same marker placement and manually selected a reference video frame in which the snake body was visually straight around that marker and horizontal. To correct for yaw offset, we manually rotated the reference marker frame about its x-axis until its forward direction was visually tangent to the local snake body in this video frame. To correct for roll and pitch offsets, we calculated a correction rotation matrix ⁠, where R_ref is the rotation matrix of the marker in the reference video frame after yaw offset was corrected, and R_yaw is a rotation matrix with zero pitch and roll angles and the same yaw angle as the marker in the reference video frame after yaw offset was corrected. Then we rotated the marker frame in each video frame by applying the correction rotation matrix to its orientation. In this process, all Euler angles were calculated after rotating body frames to match the Euler angle convention: +x-axis pointing forward tangent to segment backbone curve, and +z-axis pointing upward along the axis of left–right symmetry of the cross-section.

After reducing orientation offset, we smoothed the orientation using a window average filter covering nine adjacent frames (0.09 s), by applying a mean filter MATLAB function smooth2a (https://www.mathworks.com/matlabcentral/fileexchange/23287-smooth2a) temporally to each element of the rotation matrices. Then we used a singular value decomposition (SVD)-based projection method (Belta and Kumar, 2002) to ensure that the matrix was still in the rotation matrix space.

Finally, for each snake and each marker, we translated its body frame obtained from the marker on the dorsal surface of the snake body to the body midline, along the direction normal to the marker plane by the sum of the marker thickness and the body radius r(s) at the corresponding body length to obtain end constraints of the backbone curve.

After the snake body backbone curve was reconstructed, the kinematic data were post-processed with a temporal linear interpolation to fill missing body frames resulting from missing tracking data, followed by a temporal-spatial window average smoothing of position data to reduce discontinuities produced by the piecewise numerical convergence process.

Finally, we projected the interpolated backbone curves of all body segments onto the experimental videos and visually examined the match to validate the reconstruction.

To characterize how marker separation affects positional accuracy of the interpolation, we performed experiments to extract the midline of a segment of the snake body between two markers (position only, Fig. 2).

We constructed a 51 mm (15% snout–vent length) high step using extruded T-slotted aluminum and acrylic sheets (McMaster-Carr, Elmhurst, IL, USA). The entire test surfaces were covered by black felt (BurlapFabric, Chicago, IL, USA) and lit by work lights (Designers Edge, Union City, NJ, USA) to provide high contrast of the snake body against the background. This allowed us to use computer vision techniques to extract the body midline as the ground truth.

Eight high-speed cameras (N-5A100-Gm/CXP-6-1.0, Adimec, Eindhoven, The Netherlands) recorded the test area at 100 frame s^–1 from different views to cover the entire range of 3-D body rotation (Fig. 2A). Their videos were used to track and reconstruct 3-D position and orientation of the BEEtag markers (Crall et al., 2015) on the snake body using direct linear transformation (DLT) (Hedrick, 2008). One camera was oriented in a top view with its lens axis perpendicular to the horizontal plane (Fig. 2A, white camera). Another camera was oriented in a side view with its lens axis perpendicular to the vertical plane (Fig. 2A, gray camera). Videos from these two cameras were used to extract the midline of the snake body projected into these two planes.

We used three juvenile variable kingsnakes (body length=39.6±0.4 cm, snout–vent length=34.6±0.4 cm) for midline extraction experiments. For each of the three animals, we compared our reconstructed backbone curve with the extracted midline for four different marker spacings (3, 4, 5 and 7 cm), each with five trials. This resulted in a total of 60 trials. All animal experiments were approved by and in compliance with The Johns Hopkins University Animal Care and Use Committee (protocol RE16A223).

For both the top and side view videos, we used layered image morphology operations to extract the outline of the snake body between markers. The image operations include edging, dilation, erosion and filling.

First, we converted each video frame (Fig. 2B,C) to a binary image with white edges and black background (Fig. 2D,E) by looking for local maxima of the gradient of gray intensity using the edge function in MATLAB. Then, we connected isolated small regions with white edges resulting from snakeskin texture by enlarging their boundaries (dilation), and eroded away the boundaries of the enlarged white regions to recover their original size. Next, we filled small black regions remaining inside the snake body by making them white (filling). By only keeping the largest remaining white area between the markers, we further ensured that the two outlines (left and right for top view, dorsal and ventral for side view) extracted were clean without undesired dots.

Once the two outlines of the body were found, the projection of the midline into the horizontal and vertical planes could be recovered. We evenly interpolated the same number of points along both outlines. Then for each pair of these points, we assigned a line segment with a length equal to the body radius normal to the corresponding outline, starting at each point and pointing into the body, and connected their endpoints with a line segment. The midpoint of this line segment was considered as the midline point corresponding to this pair of points. The extracted midline was then described by the sequence of these midline points (Fig. 2D,E, dark yellow curve). The extracted midline was shorter than marker separation because the markers extrude beyond the body outline and were not included in midline extraction.

Only the position data from our reconstructed backbone curve were used for comparison with the extracted midline, which lacked orientation information. For every video frame of each trial, the reconstructed backbone curve was first projected onto the video frame (Fig. 2B,C, magenta and yellow curve) using the camera model coefficients extracted by the DLT calibration. To compare with the extracted midline that was shorter, we truncated the reconstructed backbone curve using two normal vectors intersecting both ends of the extracted midline. The truncated backbone curve (Fig. 2B,C, yellow curve) was then linearly interpolated to have the same number of points as the extracted midline for calculation of pointwise errors. These errors were then averaged for each video frame and then across video frames, weighted by the extracted midline length for each trial.

To evaluate the accuracy of our method, we compared its interpolation results with those obtained using the widely used B-spline method. Most previous animal locomotion studies (Sharpe et al., 2015; Schiebel et al., 2018; Fontaine et al., 2008) only applied B-spline interpolation to position data. Here, we applied it to position data and then calculated orientation from the interpolated position curve (see ‘Implementation of B-spline’ and Yeaton et al., 2020).

For this comparison, we used both our method and the B-spline method (see ‘Definition of B-spline’) to reconstruct the continuous body of the three kingsnakes during traversal of a large step obstacle, in which the body deforms substantially in three dimensions (both position and orientation) (Gart et al., 2019). A total of 122 trials were used for statistical analysis, with a video frame rate of 100. In each trial, 10–14 BEEtag markers (Crall et al., 2015) were placed along the snake body to track local body position and orientation in 3-D. For both methods, all data were pre- and post-processed as described above.

We used each method to reconstruct a body section in the snake body, which has two markers at its two ends and another in the middle as the approximate ground truth. Then, in each interpolated curve, we selected the body frame that corresponds to the middle marker selected. We compared the interpolated with the tracked middle marker frame to obtain interpolation error for each method.

For our method, we used the front and rear of these three markers (44–73 mm apart) as end constraints to perform backbone interpolation. The middle marker in between was not used as an end constraint but to provide an approximate ground truth of local position and orientation, relative to which interpolation error was measured. We measured the length of the body segments between the front and middle markers L₁ and between the middle and rear markers L₂. We then used the interpolated body frame that lay at the same location lengthwise in the parameter space s_mid=L₁ along the interpolated backbone curve g(s), s∈[0, L₁+L₂] to provide the interpolated local position and orientation g(s_mid).

For the B-spline method, we interpolated a description of the whole body g(σ), σ∈[0, 1] using all the markers placed along the body, except the middle one used to provide the approximate ground truth mentioned above, resulting in a total number of 9–13 markers used (variable between trials). All but the middle markers were used because B-spline interpolation will produce a straight line if only two markers are used. The middle marker mentioned above lay at a location of in the σ space, where σ_kfront and σ_krear correspond to the front and rear markers (see definition of σ in ‘Definition of B-spline’ below). We then used the interpolated body frame that lay at the same location lengthwise in the σ space to provide the interpolated local position and orientation g(σ_mid).

For both methods, deviation of the interpolated local position and orientation [g(s_mid) or g(σ_mid)] from the approximate ground truth obtained from the middle marker gave the interpolation error. We used the Euclidean distance to measure position error and the Frank Park distance (Park, 1995) to measure orientation error.

For each video frame of each trial, each marker on the body except the first and the last one was used as the approximate ground truth once to obtain the errors of each method. Errors for all reference markers from all video frames and all trials were then pooled and averaged to obtain the average error for each method.

We first used B-spline (Eqn 15) to interpolate the position of the end constraints. We used the measured x,y,z values at all the selected end constraints as sample points x_k and the respective body lengths from the end constraint closest to the head to each of these end constraints as σ_k.

Next, we calculated orientation iteratively using the interpolated position curve (Bloomenthal, 1990). We first calculated the orientation of the body frame closest to the head by assuming it matched the tangent of the interpolated position curve with zero roll. We then iteratively updated the orientation of each following body frame by rotating the previous body frame to match the local tangent around an axis normal to both the previous and the local tangent. We note that rotation can be calculated independent of position by directly interpolating the measured orientation (Belta and Kumar, 2002; Pletinckx, 1989; Park and Ravani, 1997), which was widely used in temporal interpolation. However, in this case of spatial interpolation, it does not match the tangent of the interpolated position curve.

Because cubic B-spline (m=2) is widely used in animal locomotion literature (Sharpe et al., 2015; Fontaine et al., 2008; Yeaton et al., 2020), we used m=2 in our B-spline implementation for comparison with the backbone interpolation method. We used tol=0 (interpolation) because this produced smaller errors than tol>0 (smoothing). We performed additional comparisons by varying m∈{1, 2, 3} and tol∈{0, 0.2, 0.4} × body length for interpolating p (Fig. S3).

To compare our interpolation results with the extracted midline position ground truth, for each marker spacing, weighted errors from all video frames of all trials were pooled to calculate their means and standard deviations.

To compare our method with the B-spline method, for each method, interpolation errors for all reference markers from all video frames of all trials were pooled to calculate their means and standard deviations.

To test whether the position error in top and side view planes depended on marker separation, we used a linear regression for each of these measurements, with marker separation as a continuous independent factor and each position error as a continuous dependent factor.

To test whether our method and the B-spline method differed in interpolation accuracy, we used an ANOVA for both position and orientation errors, with the method as a nominal independent factor and each error as a continuous dependent factor.

All the statistical tests followed the SAS examples in McDonald (2014) and were performed using JMP Pro 13 (SAS Institute, Cary, NC, USA).

Comparison of the interpolated backbone curve with the extracted midline via computer vision techniques showed that our interpolation method achieved high positional accuracy. The average position error over the snake body segment reconstructed (3–7 cm) in both top and side views was within 1.5 mm or 17% of body diameter (9 mm) (Fig. 2F,G). In addition, as marker separation became smaller, position error decreased in both top (linear regression, R²=0.03, P<0.0001) and side (linear regression, R²=0.17, P<0.0001) view planes. Thus, users should use the smallest marker separation possible that does not significantly affect the animal's behavior.

Our interpolation method provided a good approximation of the kingsnake's continuous body with both position and orientation information while it used large 3-D body deformation to traverse a large step obstacle, which, when combined with body surface reconstruction, enabled body–terrain contact estimation (Movie 2). Compared with the widely used cubic B-spline interpolation method (with cubic B-spline basis functions and zero tolerance; Fig. 3A, blue), our method (Fig. 3A, red) achieved higher interpolation accuracy in both position and orientation, by nearly a 50% reduction in error (P<0.0001, ANOVA; Fig. 3B,C; Movie 3). In addition, varying the degree of basis functions and tolerance for the B-spline method did not reduce its error to as small as that of our method (Fig. S3B).

This high accuracy demonstrated that, when constrained by closely spaced end constraints, an over-simplified elastic rod can still achieve a higher accuracy than purely geometric methods for interpolating the continuous body of an limbless animal body with complex morphology and biomechanics.

Note that the average errors of the two methods are larger than that of our method from the extracted midline. This may be due to three reasons. First, when we compared the two methods, the tracked marker used as reference can introduce measurement error itself. Second, error is usually larger in the middle of each interpolated segment than near its ends. Third, in this dataset, besides a small, high friction step similar to the one used in the midline extraction experiment, the snake also traversed more challenging, larger steps and low friction steps. In these conditions, the body bent more substantially (Gart et al., 2019), which may lead to greater interpolation error.

The higher accuracy of our interpolation method than the purely geometric B-spline method is attributed to its physical representation of the long, slender body as an elastic rod. For example, by using stiffness values that are in the ballpark from similar limbless animals, our method does not generate unrealistically large extension or compression, which the purely geometric B-spline method may produce (Fig. 3A). Such a higher accuracy of our method is particularly useful for quantifying body–terrain contact (e.g. calculating the ground base of support; see movie 2 of Gart et al., 2019), especially as position error reduced from approximately half to a quarter body diameter.

In addition, the backbone curve description using finite elements provides more flexibility in representing a complex shape than using a limited set of basis functions. For example, we can use different stiffnesses and cross-sectional shapes and consider tapering for interpolating different body segments or even for different finite elements within a body segment.

The higher accuracy of our method came at a cost. First, markers that can provide 3-D orientation information are usually rigid and larger than non-invasive point markers, so they take more effort to attach than painting point markers and may be difficult or impossible to use on small animals. Multiple aggregated point markers or skin texture may serve as an alternative to provide local 3-D orientation information, although this adds extra effort. In addition, the computation speed of our method was two orders of magnitude lower than that of the B-spline method (2 versus 0.02 s on average to reconstruct a segment with 700 finite elements in one video frame), mainly because of the inverse kinematics iterations used. Further, it takes the experimenter more effort to measure and fine-tune segment lengths, which is not necessary for the B-spline method. Finally, there are also small discontinuities between adjacent body segments produced by the piecewise inverse kinematics iterations with a finite threshold of convergence, although this can be resolved by smoothing in post-processing. Thus, users should weigh the benefits against these costs before applying our method. If 3-D position and orientation information are desired, but a lower accuracy is acceptable to save time, users may first try the B-spline method with orientation added following our implementation.

When applied to other species, several things should be considered. First, different species can have large differences in body geometry and stiffness properties. Users should measure length, cross-sectional shape and tapering, and choose marker spacing based on body dimensions and aspect ratio. For stiffness, they may need to tune Young's modulus and Poisson's ratio, as well as use a non-circular shape to calculate second moments of inertia and cross-sectional area. In addition, maximal inverse kinematics iteration steps may be adjusted to achieve a balance between convergence success rate and reconstruction time. Furthermore, thresholds used for the writhe and twist check depend on how much the animal can deform its body.

In addition to locomotion research, our method is useful for other studies that require quantification of continuous shape of limbless biological and artificial systems in three dimensions, such as snake predation using constriction (Penning and Moon, 2017), snake thrashing (Danforth et al., 2020), root nutation (Ozkan-Aydin et al., 2019), and octopus arm (Zelman et al., 2013) and continuum robot (Laschi et al., 2012) manipulation.

We thank Henry Astley, Bob Full and Rajat Mittal for discussion; Sean Gart for providing snake step traversal experimental data for accuracy comparison with B-spline; and Henry Astley for advice on animal care and experiments.