Analysis of the Three-Dimensional Trajectories of Organisms: Estimates of Velocity, Curvature and Torsion From Positional Information

We live in a three-dimensional world. For many organisms, this three-dimensional world is effectively two-dimensional: magnetic fields and beams of light all polarize the world, and many terrestrial organisms are constrained by gravity to a nearly planar world. Yet, for many other organisms, the third spatial dimension is as important as any other. For example, microorganisms in a puddle of water move through the water in all directions. Birds and flying insects fly over terrain and are, consequently, less constrained to a two-dimensional environment. Fish and other aquatic metazoa are similarly less constrained by gravity’s compression of life into two dimensions.

The motions of many appendages are also highly three-dimensional. The wingtip of a bird in flight, the feeding appendages of a copepod manipulating a diatom, a tree’s branches in the wind – all trace three-dimensional motions. Three-dimensional motion has, therefore, been of great interest to biologists. Nevertheless, technical difficulties in measuring three-dimensional motions and in analyzing the resulting data have limited our analyses.

The most common technique for measuring three-dimensional motion is to record positional data at discrete points in time (e.g. time, x, y and z), usually from two or more cameras having different views of the motion. Commercial systems that provide acceptable precision in the measurement of three-dimensional positional data are becoming increasingly ‘affordable’ (e.g. Motion Analysis Corporation, Santa Clara, CA, USA; Peak Performance Technologies, Inc., Englewood, CO, USA). In addition, several instruments for recording three-dimensional positional data have been developed for specific research programs (for examples, see Berg, 1978; Wehrhahn et al., 1982; Whittle, 1982; Dahmen and Zeil, 1984; Spedding et al., 1984; Spedding, 1986; Rayner and Aldridge, 1985; Strickler, 1985, 1998; Crenshaw, 1990, 1991; Tucker, 1991, 1998; Kühnel-Kratz and Häder, 1993; Britton et al., 1997; Walker and Westneat, 1997), but none has found widespread application.

The analysis of three-dimensional positional data, similarly, has no standard, a point that is addressed more fully in the Discussion. While a very small number of analyses have been published, we have not been able to find a complete treatment of the subject. The absence of convention both for collecting and, perhaps more importantly, for analyzing data has limited the study of three-dimensional motion and missed important fundamental concepts, as we will show.

This paper presents a series of algorithms for analyzing three-dimensional positional data generated by the motion of a single point in space. Employing conventions established for describing the geometry of three-dimensional curves, these algorithms provide estimates of the Frenet trihedron, speed, curvature and torsion of the trajectory – the parameters that completely describe a three-dimensional trajectory (for an introduction to these concepts, see Goetz, 1970; Gillett, 1984; Beatty, 1986). The algorithms use simple numerical differentiation using finite differences (e.g. speed is the distance between two points divided by the time interval). Curve-fitting, in the usual sense of estimating the parameters of a model function using more inputs than outputs (and, thus, permitting analysis of error in the fit), is not employed. We call this technique ‘finite helix fit’ (FHF).

The performance of FHF is evaluated using the following types of trajectories: (1) simulated trajectories for which the velocity, curvature and torsion are known; (2) simulated trajectories with noise introduced; and (3) real three-dimensional trajectories from several organisms (the flagellate Chlamydomonas reinhardtii, the ciliate Paramecium tetraurelia, spermatozoa of the sea urchin Arbacia punctulata and larvae of the ascidian Botrylloides sp.).

For the purpose of evaluating the performance of the algorithms with trajectories from true experimental conditions, the motions of several freely swimming organisms were tracked in three dimensions. Two different tracking techniques were used.

Cells were tracked using a slight modification of the technique of Crenshaw (1990, 1991). This technique uses two cameras (Panasonic WV-1854, infrared extended range Newvicon tube) with orthogonal, synchronized views to observe a three-dimensional volume. Motions of free-swimming microorganisms were videotaped (U-Matic SP videotape, NTSC video, 60 fields s⁻¹) in the center of a large (24 mm×24 mm×30 mm) observation chamber. Viscous interactions with the chamber wall were small. The only significant modification of the technique was that continuous illumination, rather than strobed light, was used.

Measurement of cell position in the video recordings was performed field-by-field using a computer (Commodore Amiga 2000 or 3000) equipped with a genlock/overlay card (Commodore Amiga 2300 Genlock) and running software developed in the laboratory of H.C.C. (copies of the software, written in AmigaBasic, are available upon request). This technique overlays the computer graphics (640×480 pixels) onto the center of the overscanned analog video image (vertical resolution 262.5 lines per field, horizontal resolution approximately 400 TV lines on the recorded image). The position of the organism in each field was the centroid of the organism as estimated by eye. Each videotaped trajectory was measured three times, and the three data sets were averaged to reduce the effect of operator-introduced errors. The single position error, averaged over the entire trajectory, was approximately 0.01 pixels; the maximum difference between repetitions at any position was 2 pixels, but this was rare. To remove motion due to convection in the chamber, a passive, neutrally dense particle was tracked simultaneously with the organism, and the motion of this passive particle was subtracted from the motion of the organism. The resulting data were the three-dimensional positions, x,y,z, of the organism recorded with a sampling frequency of 30 or 60 Hz. At the magnifications used in these experiments, the view occupied a volume of 350 μm × 350 μm (horizontal) × 500 μm (vertical), producing a spatial resolution of approximately ±2 μm in three-dimensional space.

Distortions in this system are of two types, optical distortions that distort two-dimensional measures from a single camera, and alignment distortions that distort three-dimensional measures by combining data from both cameras. Combined errors from both these distortions and from parallax were estimated by moving (rotating and translating) an object of known length in the volume and then measuring its length at multiple positions and orientations inside the volume. Variance in the measured length was always within the spatial resolution of the system, indicating that distortion and parallax errors were small.

This was similar to technique 1, but different alignment techniques and high-speed video cameras were used. The two perspectives were combined on a single split frame for analysis. The two cameras were aligned horizontally using a level and perpendicular to each other using a calibration body. The calibration body consisted of four pin heads aligned in a plane, with each pin on the corner of a square when viewed perpendicular to the plane. One camera was aligned with the pins on the first and third corners of the square and the other was aligned with the pins on the second and fourth corners.

Organisms were backlit with fiberoptic illuminators and thereby appeared as silhouettes from the camera’s perspective. The cameras were equipped with macro lenses (Nikon 55 mm or Minolta 55 mm), providing a view that occupied a volume of approximately 1 cm × 1 cm (horizontal) × 1 cm (vertical) in the center of a 5 cm × 5 cm (horizontal) × 7 cm (vertical) glass tank. These dimensions were predicted to be sufficiently large to avoid viscous wall effects (Vogel, 1994).

Distortion and parallax errors were estimated as with technique 1. The reference object was 5 mm long. Variance in the measured length was always less than 1 % of the 5 mm, indicating that distortion and parallax errors were small.

Sequences of swimming were captured by high-speed video (NAC Imaging Technology, Inc., model HSV-1000) at 500 frames s⁻¹ and recorded onto S-VHS video tape (vertical resolution 262.5 lines per field, horizontal resolution approximately 350 TV lines on the recorded image). Video sequences were digitized for frame-by-frame computer analysis using a frame synchronizer (Hotronic, Inc., model AR31) and acquisition software (Scion Corp.) on a MacOS-based computer (Power Computing, model PowerTower Pro 225) with a spatial resolution of 640 horizontal pixels × 480 vertical pixels. Frames were deinterlaced prior to analysis, resulting in a final resolution of 640×240 pixels. A 2 mm long micrometer was positioned in the center of the tank to calibrate the video images. Individuals were tracked automatically by the use of a macro programmed in the NIH Image software package (version 1.61). The images were binarized using a threshold that delineated the edges of the larva’s body from the surrounding image. This produced a 1 bit image from which the centroid of the body of the larva was determined. The centroid was recorded as the position of the larva. This macro consistently tracked the same position on the body of the animals.

Strain CC-125 (wild type) was obtained from the Chlamydomonas Genetics Center (Duke University, Durham, NC, USA) and cultured in TAP medium without bubbling (see Harris, 1989). Cultures were maintained on a 12 h:12 h light:dark photoperiod. All cultures were in the logarithmic growth phase when used. Cell culture densities were typically 2×10⁶ cells ml⁻¹ when the cultures were used. Thirty minutes prior to tracking, cells were removed from the culture containers and diluted to approximately 2×10⁴ cells ml⁻¹ with fresh medium equilibrated to 100 % air. To equilibrate the cells to the new medium, the container was then covered with plastic wrap and placed in the culture room, under lights, for 30–60 min. The containers were covered with plastic wrap to prevent cooling of the medium by evaporation. A large volume of air was trapped under the plastic to minimize changes in gas composition in the container. Equilibration was performed in the light because cells swim more vigorously than when left in the dark for this period.

Three-dimensional tracking technique 1 was used for C. reinhardtii. Infrared light (created with a band-pass filter centered on 810 nm with a bandwidth of approximately 30 nm) was used for observation because C. reinhardtii, which is photosynthetic and phototactic, is known not to respond to light of this wavelength. The culture room and the room containing the three-dimensional tracking instrument were at the same temperature (24±1 °C) to reduce the formation of thermal convection currents in the observation vessel of the three-dimensional tracking instrument. The top of the container was covered with plastic wrap to reduce evaporative cooling of the medium, which also causes convection currents. All experiments were conducted within 60 min of the end of the equilibration period.

Male sea urchins were collected from Beaufort, NC, USA, and maintained in artificial sea water at 22 °C. Mild electric shock stimulated the release of ejaculate, which was collected into glass micropipettes. Micropipettes were stored on an ice bed until the sample was used. The ejaculate was diluted to a concentration of approximately 2×10⁴ cells ml⁻¹ with artificial sea water. Three-dimensional tracking technique 1 was used. Continuous illumination with white light from a projector lamp was used. All measurements were conducted in artificial sea water at 23±1 °C. A layer of mineral oil was gently placed over the top of the sea water in the vessel to reduce evaporative cooling of the medium. All samples were used within 2 h of collection.

P. tetraurelia were obtained from Judith van Houten (University of Vermont, USA) and cultured in S3 medium (Davis et al., 1998). Cultures were maintained in the dark at 24 °C. Only cultures in the logarithmic growth phase were used. A sample of culture medium was removed from a growing culture and diluted with fresh medium to a final concentration of approximately 10⁴ cells ml⁻¹. The cells were allowed to equilibrate to the fresh culture medium for 60 min before experiments.

Three-dimensional tracking technique 1 was used for P. tetraurelia under conditions identical to those described above for C. reinhardtii.

Mature colonies of Botrylloides sp. were collected from Bodega Bay, California, USA. Within 2 h of collection, colonies were placed in holding tanks located at the University of California, Berkeley, USA, with water taken from Bodega Bay and maintained at 19 °C. In the 24 h prior to filming, holding tanks were kept in the dark. Within 4 days of collection, the release of brooded larvae was stimulated by exposing colonies to bright artificial light. During filming, individuals swam in a tank with its temperature regulated at 19±1 °C. Three-dimensional tracking technique 2 was used to analyze their swimming.

All analyses were performed with Matlab (version 5.3) in Windows NT. Copies of the software are available by request from H.C.C. or can be downloaded from the following site: www.zoology.duke.edu/crenshaw.

A three-dimensional trajectory is a three-dimensional curve. The algorithms presented here are derived from standard techniques of differential geometry used to describe three-dimensional curves. For a more complete introduction to the geometry of three-dimensional curves, see Goetz (1970) or Gillett (1984). Vectors will be presented with upper case bold letters (e.g. X) and scalars with lower case italic letters (e.g. x). Derivatives with respect to time are indicated such that X′(t) is the first time derivative and X″(t) is the second time derivative, etc. A parameter that has been produced when an algorithm is applied to data, and is thus an estimate of the parameter, is designated with an asterisk; thus, s is the actual speed and s* is the estimated speed.

Three-dimensional motion permits turns in any direction. The description of three-dimensional motion is simply an extension of the description of two-dimensional motion. Now the plane in which V turns can rotate about V. The rate of rotation of this plane is torsion τ (defined below). For example, three-dimensional motion with constant, non-zero speed s, curvature κ and torsion τ is a helix.

A series of discretely sampled points on a three-dimensional trajectory can be measured as follows. Two points define a line, so the points can be analyzed two points at a time, and the only parameter of motion obtained is V. Three points define a circle, so the data can be analyzed three points at a time, and the parameters obtained are V and κ. Four points define a helix, so the parameters obtained are V, κ and τ.

All three parameters, V, κ and τ, are needed to describe three-dimensional motion (see Goetz, 1970; Gillett, 1984).

Any three-dimensional curve H(t) can be defined by the motion of a reference frame, TNB (described in the next paragraph), relative to a fixed reference frame, XYZ. XYZ defines three-dimensional space. TNB is known as the Frenet trihedron, and its origin traces the curve. TNB rotates as it moves, causing T to change direction in three dimensions (see Fig. 1). This rotation is given by the curve’s curvature and torsion. First, we will define TNB, then we will define κ and τ.

The trajectory, however, is not known for experimentally measured motions of an organism. The data produced by such measurements typically are a series of discretely sampled positions in three-dimensional space gathered at constant time intervals: t_n, x_n, y_n, z_n (i.e. x_n, y_{n and}z_n are the spatial coordinates collected at time t_n).

The algorithms presented below permit estimation of TNB, V, κ and τ by a discrete analysis of the data. The general idea behind this approach was stated above: four points define a helix, so V, κ and τ are determined for the helix that fits four points on the trajectory. This helix exactly fits these four points.

The four points comprise a ‘window’. V, κ and τ are determined for that window and assigned to the portion of the trajectory between the second and third points. The window is then moved along the trajectory to estimate V, κ and τ along the entire trajectory.

The first step is to determine TNB at the second and third points in the four-point window. Consider Fig. 2, which presents a set of four points on a three-dimensional trajectory. These points can be represented as position vectors P₁–P₄. Points P₁, P₂ and P₃ lie in one plane, while P₄ lies off the plane. T, N and B are calculated as follows: T can be estimated at P_n by the vector connecting the two points P_n−1 and P_n+1:

The Frenet trihedron TNB, speed s, curvature κ and torsion τ completely describe a three-dimensional curve. While the estimates of these parameters are for the helix that fits the four points in the window, these algorithms can be applied to any discretely sampled three-dimensional trajectory. Some trajectories can create difficulties when applying the algorithms, and these are described in Table 1.

The four points used in this analysis define a ‘window’ (see Fig. 2). The distance they span along the trajectory delimits the portion of the trajectory used to estimate the parameters. However, neither the number of points nor the distance spanned is a useful measure of the size of this window when analyzing a helix. The size of this window is best expressed as a fraction of one turn of the helix that fits the four points.

Window size, therefore, varies with several parameters. For example, four consecutive points define the smallest window for any analysis; however, if four consecutive points are used, and the sampling frequency is constant, the window size still varies as the trajectory changes speed, torsion or curvature. Furthermore, doubling the sampling rate halves the window, or it is possible to skip points when selecting consecutive points for the window, which increases window size. All these parameters influence the analysis of a discretely sampled three-dimensional trajectory.

Systematic errors arise when w>0. Fig. 4 presents s*, κ* and τ* as dimensionless window size increases for a helix with the following parameters: s=68.3 μm s⁻¹, κ=0.0627 rad μm⁻¹, τ=0.0627 rad μm⁻¹, sampled at 30 Hz. (These parameters are similar to those for a free-swimming flagellate or sea urchin spermatozoan; see Figs 11 and 12.) Errors in the estimates of speed, curvature and torsion can become very large as w approaches 0.5.

Aliasing occurs when w equals or exceeds 0.5. This is equivalent to the Nyquist critical sampling frequency for one-dimensional and two-dimensional signals (see Press et al., 1988; Porat, 1997). For any oscillating signal, there must be at least two samples per cycle of the oscillation (see also Ward and Humphreys, 1981). At w=0.5, torsion changes sign (i.e. the handedness of the helix inverts). This is the three-dimensional equivalent of watching the wagon wheels go backwards in a movie. Aliasing prohibits the use of window sizes greater than 0.5.

Systematic errors in s*, κ* and τ* are also functions of the pitch angle, q of the helix. Fig. 5 presents the error in the estimates of speed, curvature and torsion for helices over a range of q spanning 0.001 (nearly 0; a line) to 1.56 (nearly p/2; a circle) and of w spanning 0.01 (nearly 0) to 0.49 (nearly 0.5). The effects of systematic error on speed and torsion are similar for all pitch angles: speed decreases and torsion increases. Systematic error of curvature changes with q such that, if q»1 rad, then there is little to no systematic error; if q>1 rad, then κ increases with increasing w; and if q<1 rad, then κ decreases with increasing w.

We have developed a routine that permits correction of systematic error. It is described in the Appendix. Fig. 6A presents a trajectory in which the parameters of motion vary as: s=120+80cost (40⩽s200), κ=0.04+0.035sin(2.333t) (0.005⩽κ ⩽0.075) and τ=0.04cos(1.666t) (−0.04⩽κ ⩽0.04), with a sampling rate of 60 Hz.

Fig. 6B–I presents s*, κ*, τ* and w* for windows having different point spacings. The actual values, the uncorrected estimates and the corrected estimates of these parameters are presented. The corrected estimates closely match the actual values.

Fig. 7 demonstrates how the algorithms perform as w approaches and exceeds 0.5. The trajectory is the same as in Fig. 6, but a 13-point window has been used. Fig. 7A–C presents the uncorrected estimates, and Fig. 7E–G presents the corrected estimates. Fig. 7H presents the estimated dimensionless window size w*, calculated using equation 29 from the uncorrected estimates of s*, κ* and τ*. Fig. 7D presents the actual dimensionless window size w, which approaches or exceeds 0.5 at three portions of the trajectory: 2 s, 8 s and 13 s. The behavior of the algorithms here is instructive. At 2 s, w is slightly greater than 0.5 (maximum value here equals 0.505) but w* is much greater than 0.5, reflecting the inaccuracy in estimating the true dimensionless window size from the estimated parameters. The errors in the corrected parameters are much larger here, demonstrating the poor behavior of the algorithms as w* approaches 0.5. Inversion of τ* does not occur because w is nearly 0.5. The corrected parameters are poor here because the correction algorithm is not applied when w*>0.49 because there is no way to correct estimates when w* equals or exceeds 0.5. At 8 s, w and w* are both greater than 0.5. Inversion of τ* is evident. This is accompanied by a reflection of the peak of w* at this point. At 13 s, w approaches, but does not exceed, 0.5 (maximum value 0.463). w* exceeds 0.5 (maximum value 0.516). No inversion of τ* occurs because w<0.5. Again, the corrected parameters are poor because the correction algorithm is not applied when w*>0.49.

It is obvious that the estimated parameters are unreliable when w* is greater than 0.45. In fact, values of w* that exceed 0.4 should be avoided. The errors in s*, κ* and τ* are large for such large windows, so correction of these estimates is severe. As stated, our correction routines are not applied when w*>0.49.

As explained above, the window is moved along the trajectory, providing an estimate of the parameters along the trajectory. The window can be moved in one of two ways. (i) The window is moved one point at a time, such that P_n+1 becomes P_n the next time the window is applied. This maintains the finest temporal resolution for the analysis; estimates of the motion parameters are made at each point. However, estimates of the motion parameters at any point are correlated with estimates at neighboring points because a single point resides in four consecutive windows. Analyses based on correlations must take this into consideration. (ii) The window is moved four points at a time. This eliminates correlation but also decreases the resolution of the analysis; it effectively divides the sampling rate by 4.

The first of these options is preferable in our experience. As long as the results produced by the analysis are interpreted judiciously, autocorrelation is not problematic. More importantly, point-to-point variance in the estimated parameters gives an indication of the effects of noise in the data. We provide a more complete description of how to evaluate estimated parameters in the next section.

All the examples presented thus far are for data sets obtained from smooth trajectories (i.e. the data have high precision with no noise or sampling error).

The effect of noise on these algorithms can, however, be pronounced. This is expected. Noise in positional data propagates with the time derivative, so higher-order derivatives suffer larger variance when noise is present. Speed s is the magnitude of the first time derivative of the positional data (equation 10). Curvature κ is a function of the first and second time derivatives (equation 11). Torsion τ is a function of the first, second and third time derivatives (equation 12). Thus, κ* and τ* are highly sensitive to noise, with τ* being most sensitive.

The distance separating consecutive points in the window influences the impact of noise on the estimates of s*, κ* and τ*. In effect, the ratio (noise at a point:distance between points) behaves like a signal-to-noise ratio. (i) Speed s* is estimated as the distance between two points. Thus, variance in s*, caused by noise in the data, increases as the distance between the points decreases. (ii) Curvature κ* is an estimate of the rotation of the unit tangent vector T along the length of the curve. Variance in the direction of T increases as the distance separating points decreases (see equation 13). Thus, variance in κ*, caused by noise in the data, also increases as the distance between points decreases. (iii) Torsion τ* is an estimate of the rotation of the plane formed by T and N along the length of the curve. This plane is defined at P_n by the two vectors, P_n−P_n−1 and P_n+1−P_n (see Fig. 2). In effect, noise causes the plane to rotate spuriously. Thus, variance in τ*, caused by noise in the data, increases as the distance between the points decreases. The effect of noise is more pronounced when these two vectors become more nearly parallel, which occurs as curvature becomes smaller. Thus, noise-induced variance in τ* also increases as curvature becomes smaller. Note that this variance can include spurious inversions of the sense of τ*, e.g. the helix can sporadically switch from right-handed to left-handed and vice versa. This is not to be confused with inversions of handedness caused by aliasing when w exceeds 0.5.

Importantly, the distance separating consecutive points decreases as speed s decreases for constant sampling frequency, so variance in s*, κ* and τ*, caused by errors in the measurement of position, increases as s decreases.

Variance caused by noise can be reduced by one of two strategies. (i) Increase w by skipping points in the analysis, e.g. the four points in the window are actually P_n, P_n+a, P_n+2a, P_n+3a, where a is the number of points skipped. This strategy increases the distance between points, which decreases the effects of noise, as discussed above. This approach introduces the following complications. First, the spatial resolution along the length of the trajectory is decreased because the window spans a larger portion of the trajectory. Note that this increases the window size, which can inadvertently push w towards or beyond 0.5. This can also be considered to be a reduction in temporal resolution because it effectively decreases the sampling frequency; however, estimation at each point preserves temporal resolution, as discussed above. Second, there is an atypical autocorrelation. For example, if every alternate point is skipped, and the window is moved one point at a time, then estimates at consecutive points are not autocorrelated, but estimates at every alternate point are. Again, subsequent analyses must be treated with care. (ii) The positional data can be smoothed to filter out high-frequency components of the motion. This is achieved by treating each component of space as a separate signal, e.g. H(t)=[X(t),Y(t),Z(t)]. We have used two filters: a running weighted average and a Butterworth filter. They produce qualitatively similar results. Filters, however, introduce complications (see Discussion).

In our experience with noisy data, all strategies are helpful. However, increasing window size (skipping points) is the single most powerful way of reducing the effects of noise. Judicious choice of points skipped, as described in the following paragraphs, permits s*, κ* and τ* to be determined even for moderately noisy data. Bear in mind that using larger windows increases w*.

Skipping points adds a powerful analytical tool. Consider the following. If a 2-point window is used, then the data set is effectively divided into two sets of independent observations, the odd-numbered points and the even-numbered points. A 3-point window divides the data into three independent sets, etc. Variance in the results from these independent data sets arises from noise in the data; thus, variance from point to point arises either from measurement error or from the presence of high-frequency components in the motion.

To evaluate the performance of these algorithms with a noisy data set, noise was introduced to the data used in Fig. 6 as described in the Materials and Methods section. Results from analysis of this ‘noisy’ data set are presented in Fig. 8. Columns A, B and C present s*, κ* and τ*, respectively. w* is plotted below each graph for comparison. Rows i, ii and iii present the results from three different window sizes (row i, 5-point window; row ii, 10-point window; row iii, 15-point window). Fig. 8Ai,Bi,Ci demonstrates the relative sensitivities of s*, κ* and τ* to noise; the variance in τ* is the greatest. In fact, τ* does not yield usable estimates over the intervals 3–7 s and 9–12 s for this window. Over these same intervals, the estimated window size w* exhibits large variance but remains less than 0.5 except when noise is present. In fact, this large variance in w* and τ*, especially the rapid reversal of sign of τ*, is indicative of the effects of noise.

Fig. 8Aii,Bii,Cii demonstrates the effect of increasing the window size from a 5-point to a 10-point window. The variance in κ* and τ* is reduced, but aliasing is now evident at 2.5 s and at 7 s (τ* inverts). Note, however, that the periods during which aliasing is pronounced are those periods for which usable estimates of s*, κ* and τ* were obtained with the smaller window size in Fig. 8Ai,Bi,Ci. Similarly, increasing the window to a 15-point window (Fig. 8Aiii,Biii,Ciii) produces usable estimates when variance is too large for smaller windows.

Fig. 8 demonstrates that increasing the window size for the analysis decreases variance in the estimates of the parameters. However, estimates of the parameters become unreliable, including aliasing, at larger window sizes. It is possible, nevertheless, to combine the results of Fig. 8Ci,Cii,Ciii, using sections from each where variance is small and the parameters are reliable, to obtain a single graph for torsion, as in Fig. 9.

The procedure for evaluating estimates provided by a chosen window size and for combining analyses from different window sizes is as follows. Start with a small window size. Estimates are accepted if s*, κ* and τ* (especially τ*) are well-behaved (not noisy) and w is always much less than 0.5. When an n-point window (n>1) is used, then adjacent estimates can be used to ascertain the effects of noise; remember, a 2-point window separates the data into two independent sets. If two adjacent estimates give similar values, then you have confidence that the estimates are correct. Regions where s*, κ* and τ* are noisy or where w approaches 0.5 are rejected. Next, increase the window size. Repeat the above. As Fig. 8 demonstrates, increasing the window size will reduce the noise in regions that were too noisy for smaller window sizes. If, however, w* approaches 0.5 before variance in s*, κ* and τ* becomes acceptable, then there is too much noise to yield a usable result with this technique.

Fig. 10A,B presents the trajectory of a free-swimming P. tetraurelia sampled at 30 Hz. The results from analysis of this trajectory are presented in Fig. 10C–F. Results from two different window spacings are presented (8-point window and 11-point window). The behavior of τ* at approximately 4.2 s is instructive. For the 8-point window, there is a reversal of sign that is not evident for the 11-point window (Fig. 10E). w* is well below 0.5 when the reversal occurs for the 8-point window, indicating that this reversal of sign is not due to aliasing. Inspection of the trajectory for visual evidence of such changes is important, and there is a short section at approximately 4.2 s where the trajectory has a brief bend (marked by an asterisk in Fig. 10B). Whether such a brief change is real, or is an error from data collection, is a decision that must be made by the analyst.

There are several points of caution here. First, single point outliers cannot be trusted. If points are skipped in the window, then point-by-point variance is the best way to detect outliers. Second, if successive analyses are performed, each with increasing window size, and the brief change persists, then one can have confidence that this change is real. Finally, the analyst must always be aware that increasing the point spacing makes the estimate, effectively, an average over a larger portion of the helix. Transients in the trajectory will disappear as the dimensionless window size increases; they are averaged out.

Fig. 11A,B presents the trajectory from a free-swimming C. reinhardtii sampled at 60 Hz. The results from analysis of this trajectory are presented in Fig. 12C–F. Results from two different window spacings are presented (4-point window and 7-point window). Like the P. tetraurelia trajectory, there is a region where torsion inverts at approximately 3 s. However, τ* is highly variable in this region, and w* is approximately 0.5, suggesting that this inversion of τ* is due to aliasing. Increasing the point spacing further in this region does not work because this causes w* to exceed 0.5 (not shown). Noise in this portion of the trajectory, therefore, is too large for meaningful results.

Fig. 12A,B presents the trajectory of a spermatozoan of A. punctulata sampled at 30 Hz. The results from analysis of this trajectory are presented in Fig. 12C–F. The data set was filtered with a running average: P=0.25P_n−1+0.5P_n+0.25P_n+1. Results from a 4-point window are presented. This is a relatively straight helical trajectory and, as expected, the speed, torsion and curvature of this trajectory are relatively constant.

Fig. 13A,B presents the trajectory from a free-swimming Botrylloides sp. larva sampled at 500 Hz. The results from analysis of this trajectory are presented in Fig. 13C–F. The data have been filtered with the running average used in Fig. 12. Results from two different window spacings are presented (5-point window and 10-point window). The results from the 5-point window are presented to demonstrate how τ*, and thus w*, can vary widely when small windows are used to analyze noisy data. The effects of noise are reduced when the larger 10-point window is used.

Inspection of the trajectory, however, suggests that not all the high-frequency components in the trajectory are due to noise. There is a lateral displacement about the net helical trajectory (Fig. 13A,B). This arises from lateral displacements of the point tracked on the larva due to the tailbeat.

An important aspect of this analysis is that the results (s*, κ* and τ*) can be analyzed with a Fourier transform to detect periodicity in these parameters, whereas a Fourier analysis of the positional data [e.g. of X(t)] is frequently not possible because motion in one spatial dimension is strongly influenced by motions in the other two dimensions (data not shown). For example, the Botrylloides sp. larva has a tailbeat of approximately 40 Hz, as deduced from visual inspection of the video tapes. Application of a low-pass Butterworth filter to the positional data (cut-off frequency 100 Hz) prior to analysis permits higher-resolution analysis of this trajectory with a 4-point window. Fig. 14 presents a Fast Fourier Transform (FFT) of s*, κ*, τ* and w*, respectively. First, there is an approximately 17 Hz oscillation (visible as 12 maxima, especially in κ*, over the interval 0.1–0.75 s in Fig. 13D). There are six turns in the helical trajectory, suggesting that these 12 maxima, and thus the 15–20 Hz component, arise from a distortion of the trajectory. Such a distortion may have been caused by parallax or alignment errors, as discussed in the Materials and methods section, but this would be surprising in the light of the small errors measured for these effects. (In effect, parallax and camera misalignment distort a circular helix into an elliptical helix.) In addition, there is a higher-frequency component in this trajectory, evident as a 76 Hz signal in the FFT for s*, τ* and w* (Fig. 14) and as a high-frequency oscillation in these parameters (Fig. 13C,E,F). The larva’s tailbeat frequency is approximately 35–40 Hz, as determined by visual inspection of the video tapes. In effect, the tailbeat causes the point on the larva to oscillate laterally as the larva moves forward. Each cycle of lateral oscillation produces two regions of low speed in the trajectory, so the 38 Hz tailbeat produces a 76 Hz signal in the speed. The presence of this same 76 Hz signal in the torsion data suggests that the tailbeat causes the trajectory to be superhelical. In other words, the larva’s body is not only pushed side-to-side by the tailbeat but also up and down, so the trajectory of the body is a helix with small radius and an angular frequency of 38 Hz with an axis that is twisted into the major helix seen in the trajectory.

Finite helix fit (FHF) provides estimates of the Frenet trihedron (velocity, curvature and torsion) for any arbitrary three-dimensional curve. These parameters completely define a three-dimensional curve. FHF assumes that s, κ and τ are constant over the portion of the trajectory spanned by a four-point window, which is equivalent to assuming that the first, second and third time derivatives are constant.

s, κ and τ are descriptors of the differential geometry of the curve and, therefore, not of obvious interest to most biologists. It is possible, however, to use these parameters to estimate other descriptors of motion of greater interest, such as rigid body motion.

An organism that moves through space can be considered a rigid body if the motions of appendages and deformations of the body are ignored. Thus, the body of a fly or of a ciliate is considered to be only a moving block. The motion of such a rigid body is completely described by the translational velocity V and rotational velocity W of that body. When W is described relative to the body of the organism, the three orthogonal components relative to the body of the organism are commonly referred to as ‘roll’ (rotation about the anterior/posterior axis), ‘yaw’ (rotation about the dorsal/ventral axis) and ‘pitch’ (rotation about the left/right axis). For an introduction to these concepts, see Beatty, 1986; Symon, 1971; or almost any university physics text.

The Frenet trihedron TNB can be considered a rigid body with three degrees of freedom, one of translational freedom and two of rotational freedom. The motion of a rigid body can, therefore, be understood if the orientation of TNB is known relative to the body. The following assumptions define the orientation of TNB with respect to a rigid body. (i) The point being tracked, which is thus the origin of TNB, is the center of gravity of the body. (ii) The unit tangent vector T is fixed with respect to the body of the organism. Note that this assumption restricts the organism to one degree of translational freedom; it always moves in the direction of body axis T. (iii) The unit normal N and unit binormal B vectors are fixed with respect to the body. Combined with assumption ii, this fixes the Frenet trihedron TNB to the body of the organism. The body therefore moves with two degrees of rotational freedom, one parallel Ω _‖ and the other perpendicular Ω _⊥ to the translational velocity. If the anterior/posterior axis is parallel to T, then Ω _‖ is ‘roll’ and Ω _⊥ is the resultant of ‘yaw’ and ‘pitch’.

Translational velocity V is given by equations 9 and 16 if assumption i applies.

Three-dimensional motions are performed by organisms suspended in a fluid, either air or water. Thus, the Reynolds number Re of the motion can assist in determining the validity of the assumptions listed above (for a discussion of Re, see Vogel, 1994). Assumption i is problematic in a fluid. The center of mass must include the mass of entrained fluid if the forces and moments are to be calculated (equations 32 and 37). The mass of entrained fluid can be large, especially in water and at low Re. If the centroid of the image of an organism is followed, and added masses cannot be ignored, then the kinematic descriptions of the motion (equations for the translational and rotational velocities and accelerations) are correct for the organism, but the application of equations 32 and 37 is not possible. (These equations, however, are not informative at low Re, as discussed below.)

Assumption ii requires that the organism rotate its body as the translational velocity changes direction; otherwise, the rotational velocity given by equation 33 describes the sum of the change in the orientation of the body plus the change in the direction of V with respect to the body. Assumption ii probably never applies to motions at high Re, except for the trivial case of straight-line motion. Equation 32 demonstrates that the force on the body (the sum of body forces, such as mass times gravity, and the thrust generated by the body) must change direction with respect to the body if V is to remain fixed with respect to the body. For example, a bird may fly with its beak in front, but the translational velocity V is not always beakwards; in a banking turn, where centripetal force (equation 31) is large, lift can cause V to rotate dorsally on the bird, and gravity can cause V to rotate to the bird’s ground-facing side. However, if body forces or the centripetal acceleration are small, or if the organism can change the direction of thrust production relative to its body, then the rotation of V can be small relative to the rotation of the body through the turn, and equations 31–37 will provide reasonable estimates.

Assumption ii usually applies if the motion is at low Re. At low Re, inertial forces, and thus acceleration forces, are small relative to viscous forces, so only rotation of the body causes V to change direction, assuming that the organism produces thrust in only one direction with respect to its body. Body forces on a microorganism are also usually small. Application of equations 32 and 37, however, is not informative for motions at low Re because acceleration forces are small, relative to the forces that propel the organism through the fluid, and thus constitute only a small fraction of the force on the organism.

Validation of assumption iii is similar to that of assumption ii. Assumption iii states that the unit normal vector N is fixed to the body of the organism, so N and −N emerge from opposite sides of the body, defining points U and −U on the surface of the organism. One or the other point always points towards the axis of the helix K. Thus, assumption iii is valid if one point on the organism’s body always faces the axis of the helix K (equations 25 and 26; see Crenshaw, 1989, 1990, 1993b, 1996).

The effects of noise in the positional data are of great concern for these analyses because noise creates strong variance in the estimates of curvature and, to a greater extent, of torsion. Certainly, increasing the precision of positional measurements is the best way to reduce the effects of noise. Suppression of noise can be accomplished in the analysis either by filtering the data or by increasing the point spacing of the window. Increasing the frequency of sampling will permit division of the data with a window in which points are skipped as long as the measures are independent from sample to sample. Notably, this is true even if the displacement of the object from sample to sample is less than the precision of the tracking system.

Filtering can be accomplished by applying a filter to each spatial component (e.g. applying a running average or a Butterworth filter to X(t), Y(t) and Z(t), as with Figs 13, 14). This is problematic. Filters introduce systematic underestimates of speed and curvature and overestimates of torsion (data not shown). Underestimates of derivatives in two-dimensional signals are addressed by Harper and Blake (1989), and this same behavior is evident in filtered positional data. As an extreme example, as a stronger filter is applied to a helix, it approaches a straight line. The effect is that the angular velocity Ω is preserved, but the arc length and pitch angle grow smaller. This is especially true if a low-pass filter is used with a cut-off frequency that approaches the rotational velocity Ω of the motion. These errors are independent of those introduced by the window size. Correction of these errors would require separate correction routines for each type of filter, which was not attempted here. Nevertheless, filtering of the positional data can be helpful. Judicious application is recommended.

A preferred alternative to filtering is to increase the point spacing of the window. The systematic errors introduced by a non-zero window can be large, but they can be corrected. In addition, variance in s*, κ* and τ* as the window is moved provides insight into the performance of the analysis, especially the effects of noise in the trajectory, as discussed with the analysis of real trajectories. Increasing the point spacing of the window does decrease the spatial resolution of the analysis, as discussed above, but filtering does too, because any filter produces a filtered value of a point as a function of neighboring points, sometimes many points.

An alternative approach for determining TNB, s, κ and τ is to fit the spatial components with curves. The data are treated as H(t)=[X(t),Y(t),Z(t)], and each spatial component is then treated as a two-dimensional signal to which a curve is fitted. The fitted curve is then used with equations 6–12 to define the parameters of motion. Curve-fitting reduces the effects of noise and provides quantitative estimates of the error of the fit for each spatial component. (An example of this approach is presented below.) The model function should have continuous third derivatives if torsion is continuous (equation 12), or the fit can be reapplied to derivatives of a fitted curve to obtain higher derivatives (an example is discussed below). Fitting the trajectory by segments would be required for all but the simplest trajectories to achieve a satisfactory fit (for analyses of two-dimensional curve-fitting routines applied to biomechanical data, see Winter, 1990; Walker, 1998). Rayner and Aldridge (1985) advise against the use of polynomials and spline functions for three-dimensional analyses because they are too sensitive to noise in the data, which would lead to unreliable results when differentiated at a specific point. However, Walker (1998) obtained good results with both polynomials and splines when applied to complex two-dimensional curves. Importantly, curve-fitting and filtering use similar underlying approaches and, like filtering, curve-fitting introduces systematic errors to the estimates of s, κ and τ. An example of the use of curve-fitting is discussed below.

There is, surprisingly, only one previously published technique for analysis of the three-dimensional trajectory of a single point. Rayner and Aldridge (1985) used a two-camera system to measure the trajectories of small bats as they turned in free flight. They used curve-fitting, as described above, applying techniques described by Lanczos (1956, 1957). A parabola was fitted to each dimension by least-square minimization from a five-point window. The time derivative of the parabola was then used as the estimate of H′(t) at the middle point in the window. H″(t) was then estimated using a similar approach, but now using the velocities as the inputs, which increases the window to nine points. H′(t) was used as the translational velocity (equation 9). H″(t) was used as the tangential component of the acceleration (first term on the right-hand side of equation 31). Rayner and Aldridge (1985) also calculated the curvature of the projection of the trajectory onto the X,Y plane, where Z is vertical with respect to gravity. This two-dimensional curvature approximates the three-dimensional curvature only when the Z component of motion is small, so Rayner and Aldridge limited their analysis to trajectories with level turns. An alternative approach would be to determine the curvature of the three-dimensional trajectory directly with equation 11.

The approach of Rayner and Aldridge (1985) is a good one. Note that it can be extended. If acceleration is used as an input to their algorithm, which increases the window to 13 points, the technique can be used to estimate H″(t). Equation 12 can then be applied to determine torsion and, if assumption ii applies (which it does not for the flight of bats, as pointed out by Rayner and Aldridge, 1985), equation 33 can be used to determine rotational velocity.

We have applied their algorithm, including estimation of H″(t), to the smooth trajectory used in Fig. 6 and obtained results qualitatively similar to those obtained with FHF using a 1-point window. When applied to the noisy trajectory used in Fig. 8, it was qualitatively similar to the use of FHF with a 3-point window. However, the algorithm of Rayner and Aldridge (1985) introduces systematic errors to estimates of speed, curvature and torsion, as described above. These can be large. Fig. 15 compares the estimates of speed obtained by FHF and by Rayner and Aldridge’s algorithm. The trajectory is a straight, smooth helix with the following parameters: s=204.75 μm s⁻¹, |Ω|=2 Hz, κ=0.056 rad μm⁻¹ and τ=0.024 rad μm⁻¹. This helix was sampled at different frequencies (7.5, 15, 30 and 60 Hz) to evaluate the effects of dimensionless window size on both algorithms. Estimates of speed from the algorithm of Rayner and Aldridge (1985) are not as good as the uncorrected estimates from FHF. Corrected estimates from FHF nearly overlap the actual value (the poorest estimate being 205.1 μm s⁻¹ at w=0.267).

Buelthoff et al. (1980) (technique described in Wehrhahn et al., 1982) applied the two-dimensional analysis of Land and Collett (1974) to the three-dimensional trajectories of flies by separately analyzing the two-dimensional projections of the trajectory onto the X,Y (horizontal) and X,Z (vertical) planes. The results were then used as the horizontal and vertical components of the translational and rotational velocities. This has complications. The technique of Land and Collett (1974) correctly estimates the two-dimensional translational velocity, using the equivalent of equation 16, and the two-dimensional angular velocity by measuring the angular change in the translational velocity between points and dividing by the sample interval, invoking assumption ii. This equals curvature times speed and is equivalent to equations 19 and 35 in two dimensions (remember torsion is equal to zero for two-dimensional motion). The validity of assumption ii is questioned by Buelthoff et al. (1980) and Wehrhahn et al. (1982), but they assumed that rotation of the translational velocity vector V with respect to the body of the organism was small.

Extension of the technique of Land and Collett (1974) to three dimensions, as performed by Buelthoff et al. (1980), provides correct estimates of translational velocity but misinterprets the rotational velocity. This technique measures only the curvature of the trajectory. Thus, the rotational, or angular, velocity from their analysis is similar to the perpendicular component of the rotational velocity (equation 35) and neglects the contribution from torsion (equation 34). Their technique is similar, but not identical, to equation 35 because the use of equation 19 in three dimensions introduces severe systematic errors, as discussed after equation 19. These errors vary strongly with both dimensionless window size and pitch angle (data not shown), and both these parameters appear to vary in the trajectories published by Buelthoff et al. (1980) and by Wehrhahn et al. (1982).

Finite helix fit yields estimates of the speed, curvature and torsion of a trajectory created by the motion of a single point. These parameters can be used to calculate other descriptors of the kinematics of motion, such as the translational and rotational velocities of a rigid body and of the kinetics of motion (e.g. force or moments acting on the body). Extension to these other descriptors of motion, however, requires large assumptions about the underlying motion. Consequently, results from tracking a single point on a body should be treated with caution.

When more than one point is sampled on an organism or appendage at each point in time, then other analyses are available. A complete description of the three-dimensional motion of a rigid body requires that three points on the object be sampled per point in time. For such data, the Euler angles for the rotation can be determined or a ‘finite helical axis’ analysis can be performed. Finite helical axis analysis is a technique that has recently become popular for analyzing joint motion, and it is a technique that shares many ideas with the algorithms presented here (see Woltring et al., 1994; Bull and Amis, 1998).

To correct systematic errors in the estimates of s*, κ* and τ*, we needed correction surfaces formed by having evenly spaced values either of q* or of w*. It was possible to begin a priori with evenly spaced q or w, but not for the estimates of these parameters. We constructed a surface with evenly spaced values of w* as follows. We calculated the error in the estimates of speed, curvature and torsion for a series of 32 helices varying in pitch angle from 0.01 to 1.56 (nearly zero to nearly p/2). Each helix in the series was composed of 10⁵ points per rotation of the helix. The errors were determined for the smallest window size (four consecutive points formed the window such that w* is nearly zero). The window was increased by skipping one point (a 2-point window), and the errors were determined again. This process was repeated, increasing w by one point at each iteration, until w=0.5 (not w*). This produced an array of errors for finely, but irregularly, spaced w* for each helix in the series of varying q. From this array, the w* closest to a preselected interval of w* was then written to a new array. (The difference between the actual w* and the preselected w* was never more than 0.0001.) This produced arrays of errors for fixed intervals of w* for each helix. Note, however, that these arrays had irregularly spaced q*. These arrays of regularly spaced w* but unevenly spaced q* are referred to below as the ‘correction arrays’. (The resulting two-dimensional surfaces are similar to those presented in Fig. 5; data not shown.)

The correction arrays were then used to correct errors in the estimates of s*, κ* and τ* as follows: (i) s*, κ* and τ* were determined for the window; (ii) w* and q* were determined from these parameters; (iii) the correction arrays were then used to estimate the error associated with s*, κ* and τ* as follows (described only for s*, but the process was the same for κ* and τ*). (a) The correction arrays have w* at fixed intervals. The two values of w* that bracket the measured w* were selected, providing two sets of data for which the error of s* (given as s*/s) varies with q*. (b) Each of these two sets of data was fitted with a fourth-order polynomial. (c) Using the measured q*, the error was then determined for each polynomial, providing an estimate of error at the two bracketing values of w*. (d) The error associated with the measured w* was then found by linear interpolation between the bracketing values of w*. (e) The error in s* was then removed by dividing the estimate by the error (corrected s*=s*/error).

 
• ′

  derivative with respect to time

 
• *

  superscripted; an estimated parameter

 
• <2>

  subscripted; an estimate for the region between the second and third point in the window

 
• A

  translational acceleration

 
• a

  the number of points skipped plus one between points used in the analysis window

 
• B

  unit binormal vector of the Frenet trihedron

 
• C

  equal to P_n−P_n−1

 
• D

  equal to P_n+1−P_n

 
• d

  scaling factor

 
• F

  force

 
• f

  sample frequency

 
• G

  number generated by random number generation (range −1 to +1)

 
• H(t)

  a three-dimensional curve

 
• I

  moment of inertia

 
• K

  axis of a helix

 
• l

  arc length

 
• m

  mass

 
• N

  unit normal vector of the Frenet trihedron

 
• P_n

  the nth point in a data set

 
• P_n^noise

  the nth point in a data set with noise introduced

 
• Q_<2>*

  a vector that is perpendicular to N₂* and N₃* and is, therefore, parallel or antiparallel to K

 
• Re

  Reynolds number

 
• r_i

  noise

 
• s

  speed

 
• T

  unit tangent vector of the Frenet trihedron

 
• T_n*, N_n*,

  values of T, N and B, respectively, at P_nB_n*

 
• t

  time

 
• U, −U

  point on body where N emerges

 
• V

  translational velocity

 
• w

  dimensionless window size

 
• X, Y, Z

  vectors describing three-dimensional space

 
• x, y, z

  three-dimensional coordinates in space

 
• x_n, y_n, z_n

  spatial coordinates collected at time t_n

 
• α

  rotational acceleration

 
• Γ

  moment acting on a body

 
• κ

  curvature

 
• y_<2>*

  the angle between T_<2>* and Q_<2>*

 
• θ

  pitch angle of a helix, the angle between T and K

 
• τ

  torsion

 
• Ω

  angular velocity of a helical trajectory

 
• Ω_||, Ω_⊥

  values of Ω parallel and perpendicular to T

We thank Sonja Burbano (Duke University) for assistance in collecting the data for C. reinhardtii and Gabriella Schmajuk for assistance in collecting the data for P. tetraurelia. Two anonymous reviewers provided helpful comments. G. Schmajuk was supported by a Howard Hughes Summer Research Fellowship. Stephen Katz and Mathieu Kemp provided many helpful comments. C.N.C. was supported by a predoctoral training grant from the US National Institutes of Health. M.J.M. was supported by a US National Science Foundation Predoctoral Fellowship. This work has been supported with funds to H.C.C. from the US National Science Foundation (DCB-8819271) and from the US Defense Advanced Research Projects Agency (US Army contract DAAN02-98-C-4030).