## Abstract

Here, we present a suite of photogrammetric methods for reconstructing insect wing kinematics, to provide instantaneous topographic maps of the wing surface. We filmed tethered locusts (*Schistocerca gregaria*) and free-flying hoverflies (*Eristalis tenax*) using four high-speed digital video cameras. We digitized multiple natural features and marked points on the wings using manual and automated tracking. Epipolar geometry was used to identify additional points on the hoverfly wing outline which were anatomically indistinguishable. The cameras were calibrated using a bundle adjustment technique that provides an estimate of the error associated with each individual data point. The mean absolute three-dimensional measurement error was 0.11 mm for the locust and 0.03 mm for the hoverfly. The error in the angle of incidence was at worst 0.51° (s.d.) for the locust and 0.88° (s.d.) for the hoverfly. The results we present are of unprecedented spatio-temporal resolution, and represent the most detailed measurements of insect wing kinematics to date. Variable spanwise twist and camber are prominent in the wingbeats of both the species, and are of such complexity that they would not be adequately captured by lower resolution techniques. The role of spanwise twist and camber in insect flight has yet to be fully understood, and accurate insect wing kinematics such as we present here are required to be sure of making valid predictions about their aerodynamic effects.

## 1. Introduction

Wing deformations are a distinctive, and therefore presumably important, characteristic of insect flight. In spite of this, there are few good measurements of how the wings deform during flapping. Detailed quantitative measurements of the local three-dimensional shape of static insect wings reveal an exquisitely complex architecture (e.g. Wootton *et al.* 2000; ,Sudo *et al.* 2005), which is important in determining the aerodynamic properties of the wings under steady conditions (e.g. ,Rees 1975; ,Okamoto *et al.* 1996; ,Kesel 2000). However, it is not known how such static measurements relate to the shape of the wing in flapping flight, when dynamic structural, muscular, aerodynamic and inertial forces come into play (Combes & Daniel ,2001, ,2003*a*,,*b*). At a global level, significant spanwise twist and camber have been observed qualitatively in many insects (e.g. locusts: ,Jensen 1956, ,Baker & Cooter 1979, ,Wortmann & Zarnack 1993; hawkmoths: ,Willmott & Ellington 1997; butterflies: ,Wootton 1993; crane flies, hoverflies and bees: ,Ellington 1984), but in the absence of detailed quantitative measurements for the purposes of modelling or comparative analysis, the importance of wing deformation remains unknown and unknowable. Such measurements are clearly essential for understanding the significance of this defining characteristic of insect flight.

Current techniques for measuring insect wing deformation have significant limitations. The ‘strips’ method (Willmott & Ellington 1997) models the wing as a series of rigid strips that can each rotate independently of the others, so allowing spanwise twist to be modelled. The method first fits an outline of the flat wing to a single image. The outline is then divided into a series of strips and each one of these strips is rotated until it matches the projected outline of the wing. This technique has the advantage that it requires only one camera view, but this carries with it the considerable disadvantage that individual strips must be assumed to remain flat. In consequence, any camber of the wing will be erroneously modelled as wing twist: an error analysis of this technique indicated that a 10 per cent camber in the wing resulted in a systematic error of greater than 5° in the local angle of incidence (,Willmott & Ellington 1997).

Projected laser line, or ‘comb-fringe’, methods are the existing state of the art in measuring insect wing kinematics. These techniques involve projecting a known pattern of laser lines onto an insect wing and measuring their distortion in order to estimate the wing's deformation along these lines. These techniques have been used successfully for bumble-bees (Zeng *et al.* 2000), moths (,Sunada *et al.* 2002) and dragonflies (,Zeng *et al.* 1996; ,Song *et al.* 2001; ,Wang *et al.* 2003) and have the advantage over the strips method of being able to measure wing twist and camber independently. However, the number of lines that can be projected is constrained by the practical difficulty of creating a dense grid of suitably narrow but bright laser lines, especially as the technique also requires an insect wing with suitably dense venation for the laser to illuminate. Furthermore, as the planes of light projected by the laser are fixed in a laboratory frame of reference, rather than in a frame of reference moving with the wings, it is not possible to make measurements of specific points on the wing. The velocity and acceleration of specific parts of the wing can therefore be estimated only through interpolation. This complicates the calculations of, for example, the inertial forces acting on the wing. Finally, the error in the angle of incidence and camber is strongly dependent on the angle the wing makes with the laser, which again introduces systematic error in the angle of incidence.

A fundamentally different approach to reconstructing three-dimensional objects is to use photogrammetry. Photogrammetric methods use information from two or more cameras to estimate the coordinates of an object in three-dimensional space. Photogrammetric techniques have already been used to analyse insect wing kinematics (Zarnack 1969; ,Zanker 1990; ,Fry *et al.* 2003), but without making use of modern photogrammetric approaches to solving problems of multiple-camera geometry. In the wider biomechanical literature on flight, a more modern photogrammetric technique known as direct linear transformation (DLT) has been used to reconstruct the kinematics of birds (,Warrick & Dial 1998; Hedrick *et al.* 2002, ,2004; ,Hedrick & Biewener 2007; ,Tobalske *et al*. 2007), flying snakes (,Socha *et al.* 2005) and flying squirrels (Bishop ,2006, ,2007). An important practical limitation of DLT is that it requires knowledge of the exact three-dimensional positions of a number of known calibration points with respect to each other (,Atkinson 1996). Furthermore, DLT gives approximate photogrammetric solutions (,Atkinson 1996) in the sense that it minimizes as its cost function an algebraic distance that is not in itself a geometrically meaningful quantity (,Hartley & Zisserman 2004), although the resulting geometric error can in practice be kept small provided that a suitable set of known calibration points is used (see ,§4).

Photogrammetric techniques have been used widely in clinical and veterinary gait analysis (for a review, see Allard *et al.* 1997). The use of these techniques has been facilitated by the advent of a number of turnkey systems, such as the Vicon motion measurement system (Vicon, Oxford, UK) developed for use in the entertainment industry and other motion capture applications. Unfortunately, these commercial systems are typically designed for use with large subjects wearing reflective fiducial markers under bright visible stroboscopic lighting, which makes them unsuitable for the analysis of insect wing kinematics. Our aim here is to describe the application of similarly sophisticated photogrammetric techniques to the analysis of insect wing kinematics, in order to yield results of unprecedented spatio-temporal resolution that describe the deformations undergone by the wings in flight.

In this paper, we present detailed kinematic data describing the kinematics and deformation of the wing through a single wingbeat for a sample locust and hoverfly. The full datasets comprising kinematics for a number of individuals and many wingbeats will be detailed in separate papers. There we will provide an analysis of the variation in wingbeat kinematics within and between multiple individuals, as well as placing these kinematics in the context of the dynamics of inertial and aerodynamic force production.

The nomenclature used in this article is shown in table 1.

## 2. Material and methods

### 2.1 Wind tunnel experiments with tethered locusts

Desert locusts (*Schistocerca gregaria*, Forskål) were selected for flight experiments from the gregarious breeding colony at the University of Oxford, Department of Zoology. Individuals were selected on the basis of having good wing condition and strong flight motivation. The right hindwing of each locust was marked with approximately 100 spots of 1 mm in diameter, using a fine-tipped, black permanent marker pen. The forewing was left unmarked as its natural pigmentation allowed points on its surface to be identified without the need for marking. The locusts were tethered ventrally by the plastron of the thorax using a brass mounting plate and cyanoacrylate adhesive. The mount was securely attached to a six-component strain gauge force–moment balance (I-666, FFA Aeronautical Research Institute of Sweden) in a low-speed wind tunnel, with a working section of 0.5×0.5×1 m, designed specifically for insect flight experiments. The wind speed was set to 3.3 m s^{−1} and the locust was set at a body angle of 9° (measured from the underside of the thorax relative to the free stream). These conditions have been shown to be close to the equilibrium flight speed and body angle for locusts (Taylor & Żbikowski 2005).

Two NAC Hi-DCam II cameras (NAC Image Technology, CA, USA) with Nikkor 50 mm lenses and two Photron Ultima APX cameras (Photron Ltd, Bucks, UK) with Nikkor 60 mm macro lenses were used for filming. The Photron Ultima APX camera was recorded at a resolution of 1024×1024 pixels and the Hi-DCam II recorded at 1280×512 pixels. Two cameras, one of each model, were positioned above and behind the tethered locust. The other two cameras were positioned beneath and to the right of the locust. This gave good views of the right-hand side of the locust with minimal overlap of the forewings and hindwings. For illumination, the locust was silhouetted using two ARRI 125W pocket lights pointed onto white cards on the other side of the locust. This gave even white backgrounds to the images, which contrasted well with the pigmentation and venation of the locust's wings. The lights were aligned with the centre line of the wind tunnel to reduce steering efforts by the locust. Recordings were made only when the locust was in the complete flight posture (Weis-Fogh 1956), with both hindlimbs folded underneath the abdomen, all limbs free of the tether and clear periodic forces shown on the force balance. The cameras were triggered manually and were set to record for 1000 frames at 974 fps, with synchronization provided by a timing pulse from one of the NAC Hi-DCam II cameras. A shutter speed of less than 100 μs was used in each camera, which was sufficient to eliminate motion blur.

In total, we used two male and two female locusts, analysing the first five complete wingbeats from 10 separate recordings, giving a total of 50 wingbeats for which approximately 100 points on the hindwing and approximately 15 points on the forewing could be identified. After each experiment, several synchronized images were recorded by each camera of a calibration object consisting of a flat rectilinear grid of black circular dots (22×22 dots of 1 mm in diameter spaced 2 mm apart, with a spacing accuracy of 0.0025 mm), which was held in a range of random orientations and positions in the field of view. The grid was printed on a thin sheet of glass to enable visualization from both the sides.

### 2.2 Experiments with free-flying hoverflies

Hoverflies (*Eristalis tenax* L.) were wild caught in Oxford and used for experiments on the day of capture. Similar to the forewings of the locusts, those of the hoverflies were not marked for experiments because the mass of the markers would not have been negligible in relation to the mass of the wings, and the intersections of the veins already provided a reasonable number of identifiable points for analysis. Each hoverfly was allowed to fly freely inside a 220 mm diameter flight cylinder made of 0.175 mm clear sheet polyester, and was encouraged to hover in a column of visible light generated by two light sources, placed above and below the light chamber. To prevent flicker, we used a bright AC light with high thermal inertia above the insect and a dim DC light below. The hoverflies placed inside the cylinder would fly between the two lights and often hover below the surface of the cylinder. If the hoverfly chose to land on the cylinder, then a gentle tap on the outside would cause it to resume flying.

Four Photron Ultima APX cameras were used for filming, recording at 4000 frames per second with 1024×512 pixel resolution. The cameras were positioned above the cylinder to give good views of the wings throughout the stroke cycle, with minimal obscuration by the body. The ambient lighting used to encourage hovering was insufficient to illuminate the high-speed images, and instead a 200 W infrared pulsed laser (HSI-1000; Oxford Lasers Ltd, Oxford, UK) was used. The laser was routed through four 1 mm diameter fibre optics and collimated by a series of lenses into beams approximately 50 mm wide. These were pointed directly into each of the cameras' lenses to provide extremely bright back illumination of the subject. The brightness of the laser was sufficient to enable small apertures (f16–f32) to be used so as to maximize the depth of field, while also having a short enough pulse duration (20 μs) to eliminate motion blur and overheating of the insect. The wavelength of the laser (805 nm) was beyond the range of the visible for *E. tenax* (Bishop 1974; ,Horridge *et al.* 1975), so did not interfere with the hoverflies' behaviour. Synchronization of the cameras and laser was provided by a timing pulse from one of the cameras.

The target volume of the cameras formed a sphere of approximately 40 mm in diameter, centred 50 mm below the top of the flight cylinder. The cameras were triggered when the hoverfly broke a pair of crossed laser beams generated by two red laser pointers, and detected by a pair of photodiodes. The cameras save images to a continuous circular buffer that allows images captured prior to triggering to be saved. We post-triggered the cameras so that they saved images between 0.5 s either side of the trigger. In total, 172 recordings were made from 26 individual hoverflies, though only a subset of 78 recordings, from six individuals, had both wings visible in at least two cameras for at least one complete stroke cycle. From these 78 recordings, five sequences, each from a different hoverfly, were selected for the full kinematic analysis, providing 20 complete wingbeats. These sequences were classified as hovering flight under Ellington's (1984) criterion for hovering as flight with an advance ratio of less than 0.1. After each experiment, several synchronized images were recorded by each camera of a calibration object consisting of a flat rectilinear grid of black circular dots (18×18 dots 0.75 mm in diameter spaced 1.5 mm apart, with a spacing accuracy of 0.013 mm), which was held in a range of random orientations and positions in the field of view.

### 2.3 Photogrammetric model

Figure 1 shows the central perspective projection on which the photogrammetric model for each camera is based (e.g. ,Atkinson 1996; ,Hartley & Zisserman 2004). The target point *X*_{A} is projected onto the projection plane (i.e. the camera sensor) where it creates the image point *x*_{a}. The line passes through the perspective centre, *X*_{O}. The line orthogonal to the projection plane and passing through the perspective centre *X*_{O} is called the principal axis, and the point at which it intercepts the projection plane is called the principal point. The distance between the principal point and the perspective centre *X*_{O} is called the principal distance, *c*. Two coordinate systems are needed to explain the model further. The object coordinate system {*X*,*Y*,*Z*} is located arbitrarily in three-dimensional space. Within this coordinate system, the perspective centre and the target point *X*_{A} have the coordinates (*X*_{O},*Y*_{O},*Z*_{O}) and (*X*_{A},*Y*_{A},*Z*_{A}), respectively. The image coordinate system {*x*,*y*,*z*} has its origin at *X*_{O} and its *z*-axis coincides with the principal axis. The *x*- and *y*-axes of the image coordinate system are parallel to the projection plane and directed along the horizontal and vertical axes, respectively, of the camera's sensor. The image point *x*_{a} has the coordinates in the image coordinate system.

The central perspective projection model is described mathematically by the collinearity equations, so called because they express the collinearity of any point *X*_{A}, the perspective centre *X*_{O} and the corresponding image point *x*_{a} on the projection plane of a given camera,(2.1)(2.2)where *r*_{11}, *r*_{12}, etc. are the elements of the 3×3 rotation matrix ** R** mapping the orientation of the object coordinate system onto that of the image coordinate system. These are standard photogrammetric equations (for their derivation, see Atkinson 1996), and the central problem in photogrammetry is to estimate the object coordinates of multiple target points by solving the system of simultaneous equations that results when there are multiple cameras. Generalizing to the case of

*k*=1, …,

*m*cameras and

*i*=1, …,

*n*target points with position vectors

*X*_{i}=(

*X*

_{i},

*Y*

_{i},

*Z*

_{i}), we may write(2.3)(2.4)where each target point forms an image point in each camera, with the position vector . Only two camera views are necessary in order to solve these equations, but, in the presence of measurement error, the technique is more accurate if a greater number of cameras are used. A complete solution to these equations requires there to be certain known dimensional constraints upon the relationships of a subset of the target points. This is provided here by the calibration images, in which the geometry and distances between grid points are known

*a priori,*even though in each calibration image the position and orientation of the grid is not. In principle, a complete solution to the equations can then be found, although, in practice, the presence of measurement error means that this must be an approximate solution. Given the nonlinearity of the equations, this solution must be determined by numerical optimization. The optimization procedure used here is known as a bundle adjustment. This is one of the most widely used techniques in photogrammetry, principally because it makes best use of all available information and requires no prior knowledge of camera position and orientation (e.g. Atkinson 1996; ,Hartley & Zisserman 2004).

### 2.4 Camera calibration

A bundle adjustment procedure (e.g. Atkinson 1996; ,Hartley & Zisserman 2004) produces jointly optimal estimates of the camera parameters (i.e. the principal distance *c*_{k}, position vector *X*_{O,k} and rotation matrix *R*_{k} for each camera) and the object coordinates (*X*_{i}), given the known image coordinates (*x*_{i,k}) and any prior knowledge about the relationships between them (e.g. the dimensions and geometry of the calibration grid). These estimates are optimal in the sense that they minimize the difference between the image coordinates of any identified points and their reprojected image coordinates based upon the estimated object coordinates and camera parameters. This difference is known as the reprojected pixel error (RPE). The bundle adjustment used here was based on a nonlinear least-squares optimization of equations (2.3) and ,(2.4), with the addition of a set of constraint equations specifying the dimensional relationships between the grid points in each calibration image (see ,appendix A). The image coordinates of the grid points in the calibration images were digitized automatically using custom-written software in Matlab (Matlab v. 7.4, The Mathworks Inc., Natick, MA) to find the centre of each circular dot. All subsequent analyses were also performed using Matlab.

The addition of a set of constraint equations for the calibration grids reduces the number of unknowns in the collinearity equations by treating the object coordinates of the grid points not as independent points, but as a set of points specified by the rotation and translation of a grid of known geometry and dimensions. This means that, for each calibration image, only the position and orientation of the grid as a whole is solved for, rather than the object coordinates of all of the grid points independently. Treating the calibration grid as a single object of known geometry and dimensions provides a scale reference in the bundle adjustment, without the need for multiple control points of precisely known relative position in three-dimensional space, as would be required, for example, by DLT.

Convergence of iterative nonlinear optimization procedures is improved by specifying reasonable initial estimates of the unknowns. These initial estimates were obtained by solving the constraint equations for each camera individually, so as to give initial estimates of the principal distance *c*_{k} and the position and orientation of the calibration grid in each image relative to the image coordinate system of the camera viewing it. These were then used to provide initial estimates of the position vector *X*_{O,k} and the rotation matrix *R*_{k} of cameras *k*=2, …, 4 relative to camera *k*=1. A full bundle adjustment was then run for all four cameras simultaneously, using these initial estimates as starting values for the unknown camera parameters. This provides the estimates of the camera parameters that are optimal in a least-squares sense and amounts to a complete calibration of the cameras. Given this calibration, it is possible to take a new set of image coordinates (i.e. the points on the insects' wings) and to estimate their object coordinates as a nonlinear least-squares solution of the collinearity equations with the camera parameters fixed. Fixing the camera parameters in this way means that the error in the image coordinates of one target point will not affect the accuracy of the object coordinates estimated for another target point. This would not be the case if the camera parameters were left free to vary.

The photogrammetric model presented so far is the simplest possible model and essentially assumes a pinhole camera. It is possible to introduce parameters accounting for lens distortion and offset of the camera sensor from the principal axis, but the effect of including such parameters was expected to be minimal in this instance, because the sensor sizes of the cameras were small compared with the focal length of the lenses used, and the lenses were stopped down to small apertures. Introducing distortion parameters also increases the order of the optimization problem, making convergence slower and increasing the risk of convergence upon a solution that is locally rather than globally optimal. By way of validation, we repeated the bundle adjustment for a subset of the data while allowing for lens offset, second-order tangential distortion and sixth-order radial lens distortion (Brown 1966; ,Atkinson 1996). The inclusion of these additional camera parameters made negligible difference to the final estimates of the object coordinates and basic camera parameters, and we therefore excluded them from the main analysis.

### 2.5 Tracking

The data for the four locusts were analysed using custom-written semi-automatic tracking software. Image coordinates of marked points on the hindwings and identifiable natural pigmentation on the forewings were first digitized manually for one wingbeat for each of the locusts (figure 2*a*). Some of the points became invisible at times (e.g. as the hindwing folded into the body on the upstroke), and were treated as missing data points. The three-dimensional object coordinates of the manually digitized points were then estimated using the methods described in §2.4. The resulting kinematics were used to provide a reference wingbeat for each locust to enable efficient automatic tracking of the same points on other wingbeats. The purpose of identifying a reference wingbeat was to solve the data association problem of automated multi-target tracking, in which it is difficult to be certain of tracking the same individual target, out of a set of targets, from one frame to the next. We used the reference wingbeat to provide a fixed template image of the appearance of the point at different stages of the wingbeat, and to predict the area of the frame in which a given point would fall on other wingbeats. Automatic pattern recognition software was then used to identify the exact image coordinates of the tracked point.

In order to predict the area of the frame in which a given target point would fall, the image coordinates of the wing tip were digitized manually for every frame. This was the only manual input that the semi-automatic tracking software required. The resulting three-dimensional object coordinates were used to identify the two closest matching frames of the reference wingbeat, and the expected three-dimensional positions of the various target points were linearly interpolated between frames and projected back into the image coordinates using equations (2.3) and ,(2.4). For each point, a template image of that point in the reference wingbeat was created from the closest matching frame. This template image was compared with a larger search area around the expected position of the point. The exact image coordinates of each target point were estimated by finding the portion of the search area with the lowest mean square difference in greyscale value from the template image.

Tracking of the hoverfly data could not be readily automated, because the wing kinematics used by the hoverflies were too variable to enable the use of a reference wingbeat. Instead, we manually digitized the image coordinates of 12 anatomically identifiable points in the venation of each wing (white circles in figure 2*b*). The image coordinates of a further 10 points around the wing outline were identified using a technique known as epipolar geometry, which is described in §2.6. Finally, the image coordinates of the base of the antennae, the tip of the labrum and the tip of the abdomen were digitized manually to give measurements of the body kinematics. The three-dimensional object coordinates of the manually digitized points were then estimated using the methods described in ,§2.4. A limitation of using only natural markers in the case of the hoverfly is that the wing outline has few anatomically identifiable features. We overcame this problem by using a technique known as epipolar geometry.

### 2.6 Epipolar geometry

Epipolar geometry refers to the intrinsic projective geometry between two camera views (Hartley & Zisserman 2004) and can be used to simplify the problem of identifying the image coordinates of an object in one calibrated camera view if its image coordinates have already been identified in another camera view (,figure 3). For example, having taken the image coordinates of an arbitrary point in one camera view, it is possible to identify the epipolar line along which the same point must lie in a second camera view. In the case that the arbitrary point is known to lie on the outline of the wing, the intersection of the epipolar line with the wing outline precisely defines the image coordinates of that point. We used this method to determine the image coordinates, and hence object coordinates, of multiple anatomically indistinguishable points along the outline of the hoverflies' wings. We first digitized 10 arbitrary points along the wing outline (shown as white lines intersecting the wing outline in ,figure 2*b*) in one camera view. We then used epipolar geometry to draw the epipolar lines corresponding to these points in the other three camera views. This allowed us to identify the image coordinates of the same object points in the other three camera views as the intersections of the epipolar lines with the wing outline. More generally, this technique can be used to determine the three-dimensional position of points on an object with linear or curvilinear features where it is not otherwise possible to match point features between camera views; for example, a similar technique was recently used to identify points along the midline of an octopus's arm by Yekutieli *et al.* (2007).

### 2.7 Reconstruction of wing topography

Once the three-dimensional coordinates of the tracked points had been determined, they were forwards–backwards filtered using the coefficients of a third-order low-pass Butterworth filter with a −3 dB cut-off frequency of 150 Hz for the locusts and 1000 Hz for the hoverflies. We determined which cut-off frequency to use by examining the autocorrelation function of the residual noise removed by filtering, and selecting the lowest cut-off frequency at which the autocorrelation function had the properties of white noise. This is the lowest frequency at which random noise is filtered from the data without removing any of the underlying signal, and, in this sense, is the optimal cut-off frequency to use. Frequency-based filtering is appropriate in this case because it preserves all of the magnitude and phase characteristics of the signal up to about the seventh or eighth harmonic of wingbeat frequency, while eliminating most of the error. This error is expected to be uncorrelated from one frame to the next, by removing any frequency content down to two or three times below the Nyquist frequency.

We then used three-dimensional cubic splines to fit the wing outline, and also to fit the main wing veins in the locust. A 100×100 interpolated point mesh was then fitted to these splines to provide a surface map of the wing. The mesh was spaced evenly along the wing and the morphological angle of incidence at each spanwise station was measured, defined as the angle between the horizontal plane and the chordline joining the leading and trailing edge. This mesh accurately described the topography of the locust wing surface and was of sufficient resolution to visualize the corrugation between the veins in the case of the hindwing fan. In the case of the hoverfly, fewer points were tracked on each wing, and having fitted a 100×50 point mesh in the same way as for the locust, we then fitted a more accurate surface mesh by modelling the wing as a deformable surface without stretch. This was done by generating a flat, two-dimensional spline model of the wing veins and the outline from a scanned image of the wing of each hoverfly. This spline model was then distorted in three dimensions to fit the original 100×50 point surface mesh, while preserving the total surface area of the wing. A new 100×100 point mesh was then fitted to the veins and the wing outline to create an accurate topographical mesh of the wing. Figure 2*c* shows the veins and the wing outline projected back into two dimensions using equations (2.3) and ,(2.4). As can be seen, there is a close match between the veins of the model and the veins of the image, indicating that the model provides an accurate fit to the real hoverfly wing.

## 3. Results

### 3.1 Topographic wing kinematics

The full datasets, totalling 50 wingbeats from five locusts and 20 wingbeats from four hoverflies, will be detailed in separate papers. In this section, we provide a quantitative summary of the changes in spanwise twist and camber that occur through a sample wingbeat from one of the locusts and one of the hoverflies. These wingbeats were selected because they are representative of the data as a whole, and allow us to compare our results with those obtained using other methods for measuring changes in camber and/or spanwise twist in §4. ,Figures 4 and ,5 plot a surface map of the wing at 10 stages of the sample wing strokes. The surface maps in ,figures 4 and ,5 are generated using 100×100 point meshes, so can be thought of as plotting the morphological angle of incidence of a set of 100 cambered blade elements for each wing. The blade elements are coloured according to the morphological angle of incidence of the local chordline with respect to the horizontal, so as to visualize the pattern of spanwise twist. In the case of the locust, the chordlines were taken parallel to the sagittal plane of the insect, while for the hoverfly, the chordlines were taken perpendicular to the long axis of the wing. This meant that the chordlines were approximately aligned with the relative flow, whether this was due primarily to the freestream, as for the locust, or primarily due to the sweeping of the wings, as for the hoverfly. The morphological angle of incidence was measured from the anatomical undersurface of the wing, so that when the wing was inverted on the upstroke for the hoverfly, the morphological angle of incidence was greater than 90°. This definition simplifies the description of morphological twisting of the wing, but note that, for aerodynamic purposes, the angle of incidence would be measured with respect to the functional, rather than anatomical, undersurface. ,Figures 6 and ,7 plot only those blade elements falling at 10 per cent intervals along the length of the wing, coloured according to the morphological angle of incidence of the local mesh elements, so as to visualize the local aerofoil sections and camber distribution.

#### 3.1.1 Tethered locusts

The locust forewing displays prominent spanwise twist on the upstroke, when the twist is in a positive sense from root to tip (i.e. the angle of incidence increases from root to tip), but is relatively untwisted on the downstroke, with an approximately constant morphological angle of incidence across the span (figure 4). At its maximum, the forewing is twisted by approximately 30° from root to tip. This occurs at , where is the proportion of the way through the wingbeat from the start of the hindwing upstroke. By contrast, the hindwing shows prominent spanwise twist on both the upstroke and the downstroke, but whereas it is twisted in a positive sense from root to tip through most of the upstroke, it is twisted in a negative sense throughout the entire downstroke (figure 4). At its maximum, the hindwing is twisted by approximately 20° on the upstroke, and its peak twist approximately coincides with that of the forewing at . On the downstroke, the hindwing is maximally twisted by approximately −20°, which occurs at .

The locust forewing is positively cambered across its entire span for almost all of the wingbeat (figure 6). Consistent with ,Jensen's (1956) classic result, a z-shaped profile was observed only during the latter stages of the upstroke (=0.3–0.5), and then only on the most proximal sections of the wing. This ‘z-profile’ is due to the bending of the medial and vannal folds, marked ‘v’ and ‘m’, respectively, in figures 4 and ,6. The positive camber present during the rest of the stroke was maximal close to the wing root, and peaked at the middle of the downstroke () at approximately 8 per cent, measured as the maximum deviation of the aerofoil perpendicular to the chordline expressed as a percentage of the local chord length. The degree of camber decreases towards the tip, as the medial and vannal folds of the forewing, which contribute to its camber, become reduced in size. The locust hindwing is positively cambered throughout the entire stroke, except close to the tips, where the camber briefly becomes negative part way through the upstroke (figure 6). As with the forewing, the degree of camber decreases towards the tip. Unlike the forewing, which is most strongly cambered at the middle of the downstroke, the hindwing only reaches its maximum camber at the end of the downstroke (), when it reaches approximately 10 per cent near the root.

#### 3.1.2 Free-flying hoverflies

The hoverfly wings display prominent spanwise twist throughout the entire wingbeat (figure 5), which is strongest at stroke reversal when a torsional wave passes along the wing, reversing the sign of the twist. Morphologically, the twist is positive from root to tip on the upstroke and negative on the downstroke, but in aerodynamic terms, the twist is effectively negative throughout most of the wingbeat. This is because the inversion of the wing on the upstroke means that an increase in the morphological angle of incidence (measured from the anatomical undersurface to the horizontal) amounts to a decrease in the aerodynamic angle of incidence (measured from the functional undersurface to the horizontal). Thus, the torsional wave at stroke reversal maintains an aerodynamically similar twist distribution on both the upstroke and the downstroke. At its maximum, the morphological twist exceeds 50° from root to tip. The maximum morphological twist occurs following the reversal from the downstroke to the upstroke at . Even during the translatory phases of the wingbeat, the morphological twist of the wing is of the order of ±20°(at =0.33 and 0.78).

The hoverfly wings are strongly cambered throughout the entire wingbeat (figure 7), except in the most distal regions, which remain approximately flat throughout the stroke. The sign of the morphological camber rapidly switches from negative to positive between the upstroke and the downstroke, which in aerodynamic terms means that the wing is positively cambered throughout the wingbeat, given that the wing is inverted on the upstroke. Maximum camber occurs at the reversal from the upstroke to the downstroke () at approximately a quarter of the way along the wing, where it reaches approximately 12 per cent, measured as the maximum deviation of the aerofoil perpendicular to the chordline expressed as a percentage of the local chord length. During the translatory phases of the wingbeat, the highest camber is found nearer the middle of the wing, and reaches a maximum absolute value of approximately 8 per cent (). The ‘alula’ (marked ‘a’ in figures 2*c*, 5 and ,7) is a simple hinged flap at the base of the wing whose movements do not appear to have been described, even qualitatively, elsewhere. Although the alula is most often parallel to the plane of the wing, it occasionally flips to become almost perpendicular to the wing (as in ,figures 5 and ,7).

### 3.2 Error analysis

A key advantage of photogrammetric methods over the strips and projected laser line methods is that as well as being able to provide an estimate of the typical three-dimensional error, it is also possible to estimate the error associated with every single data point. This can be done by making use of the difference between the actual and reprojected image points, which is known as the RPE. We use this here to provide a detailed error analysis of the method we have described. There are three principal sources of error when using photogrammetric techniques to determine the three-dimensional coordinates of multiple targets: systematic error arising from the calibration if the camera parameters are incorrectly estimated; random error arising from incorrect data association if targets are not correctly tracked; and random error arising from image pixelation and errors in finding the exact centre of the same target in different camera views. All of these will result in error in the three-dimensional coordinate measurements, and so their effects must be quantified in order to determine the accuracy of the method. We will deal with each in turn.

#### 3.2.1 Calibration

Images of the calibration grid that had not been used in estimating the camera parameters were used as validation data. The camera calibration parameters were used to determine the three-dimensional coordinates of the grid points for these validation data, treating each grid point as an independent data point. We used the measured mean dot spacing to estimate the systematic error in the scaling of the calibration. The measured mean dot spacing was found to be 2.0005 and 1.50003 mm for the grid used for the locusts and the hoverflies, respectively. Both of these estimates are well within the stated accuracy of the spacing of the calibration grids (2 mm±0.0025 and 1.5 mm±0.013 mm, respectively). Furthermore, the measured local dot spacing was only weakly correlated with the three-dimensional location in space of the dots (for the locusts, the product moment correlation coefficients of the dot spacing with the *x*-, *y*- and *z*-coordinates of the dots were −0.004, 0.068 and −0.033, respectively; for the hoverflies, the correlation coefficients were −0.030, −0.017 and 0.057, respectively), indicating that there is negligible bias in the scaling across the object coordinate space, as would be caused, for example, by failing to account adequately for lens distortion.

#### 3.2.2 Automated tracking

Any error arising from incorrect data association between cameras results in an apparent inconsistency between camera views. In this case, the least-squares optimization still produces a single estimate of the three-dimensional coordinates of the object, but when projected back into the view of each camera as reprojected image coordinates, there will be a large difference between the reprojected image coordinates and the actual image coordinates of the identified points. In the automated tracking software used for the locusts, any points with a RPE of above 4 pixels were assumed to be errors. Overall, the tracking software had a success rate of approximately 96 and 92 per cent for the hindwing and the forewing, respectively, the lower value for the forewing being due to the increased difficulty in tracking natural features compared with high-contrast, near-circular marked spots. Screened outliers were individually examined by a human operator, and either manually digitized or treated as missing data points if they were not sufficiently visible to be accurately digitized. Since the outlier rejection threshold was chosen to be approximately the diameter of the marked spots on the hindwings, we can be confident that there are no incorrectly associated data points remaining in the data after screening.

#### 3.2.3 Image coordinate errors

The RPE provides a direct estimate of the image coordinate error of a given point in a given camera view, but it does not give a direct estimate of the error in object coordinate space, which is more useful for assessing the accuracy of the three-dimensional data. In order to determine the relationship between the RPE and the three-dimensional error, we first reprojected the image coordinates for all of the measured three-dimensional data points. We then added Gaussian noise to these image coordinates in both the *x*- and *y*-directions, and estimated the new three-dimensional object coordinates and new RPE. We then used the distance between the new and old three-dimensional object position as a measure of the three-dimensional error resulting from the known added noise. This was done as a Monte Carlo simulation, in which we adjusted the variance of the added noise. The three-dimensional error varies slightly among the different axes of the three-dimensional object coordinate system due to their differing alignments with respect to the cameras' two-dimensional projection planes in which the image coordinates are measured. Nevertheless, the arrangement of the cameras is such that the mean three-dimensional error provides a good summary of the measurement error.

Figure 8 plots the mean three-dimensional error against the mean RPE for all of the data points combined, at each level of variance of the added noise (note that both the three-dimensional error and the RPE are unsigned, since they measure distances). The slope of the linear regression is the best estimate of the relationship between the RPE and the three-dimensional measurement error across the object coordinate space, and can be used to estimate the three-dimensional measurement error of a given data point based upon its own RPE. ,Figure 9*a*,*c* plots the histograms of the RPE for the locust and hoverfly data points, respectively. Using the regressions from figure 8, we can convert from units of pixels (in image coordinates) to millimetres (in object coordinates), as shown in ,figure 9*b*,*d*. The mean three-dimensional error (table 2) for all of the locust data points was 0.11 mm (approx. 0.24% of the wing length, corresponding to a RPE of 0.88 pixels), while the mean three-dimensional error for the hoverfly was 0.030 mm (approx. 0.20% of the wing length, corresponding to a RPE of 0.46 pixels). The lower absolute error for the hoverfly data is due to the smaller working volume, but because the wings are an order of magnitude smaller, the percentage error is similar.

We tested the effect of these three-dimensional measurement errors on the measured angles of incidence by generating a topographic map of the wing surfaces and calculating the angles of incidence for the run of the Monte Carlo simulation in which we had added Gaussian noise corresponding most closely to the RPE in the actual dataset. We then compared these angles with those that had been measured for the actual data and used the difference between them as an estimate of the angle error. The angle error is an inverse function of the local chord length, but, at worst, had a standard deviation of 0.40°, 0.51° and 0.88° for the locust hindwing, the locust forewing and the hoverfly, respectively (table 2).

## 4. Discussion

### 4.1 Importance of measuring topographic wing kinematics

The wing kinematics we have presented here are of unprecedented spatio-temporal resolution, and demonstrate the prominence of variable spanwise twist and camber in two morphologically and phylogenetically distant insects. In aeronautical engineering, twist is an important aspect of wing design. In the case of rotating wings, such as propeller and helicopter blades, a washout distribution (i.e. the decreasing angle of incidence from root to tip) is used to counteract the increasing angle of attack from root to tip that would otherwise result in forward flight from the increasing ratio of blade self-motion to freestream velocity. A higher degree of washout can also increase the efficiency of a helicopter in hover, albeit at a cost to forward flight performance (Leishman 2000). In a fixed-wing aircraft, a washout distribution reduces bending moments by concentrating lift at the root of the wing, increases lateral stability by reducing the likelihood of tip stall and can provide longitudinal stability if combined with backward sweep of the tips (,Mises 1959). Unfortunately, little is known about the aerodynamic effects of wing twist under unsteady, separated flow conditions at low Reynolds numbers. Nevertheless, our locust and hoverfly wing kinematics display such prominent spanwise twist through the wingbeat that it surely cannot safely be neglected in future aerodynamic analyses (,figures 4 and ,5). ,Figure 10 plots the local angle of incidence against spanwise station at mid-downstroke, mid-upstroke and the stroke reversals for both the insects. The graphs show that spanwise twist is both pronounced and complex, with evidence of significant high-order nonlinearities. This makes it vital that future investigations of the aerodynamic effects of wing twist use kinematic data with sufficiently high spatio-temporal resolution to capture these complexities.

The effect of camber is complex. Suitable positive camber can increase the lift-to-drag ratio at a given angle of attack and reduce the risk of flow separation by generating a more favourable pressure distribution. Camber can also enhance the stiffness of a wing, and can provide longitudinal stability if the trailing edge is reflexed upwards. Hence, nearly all modern aircraft use cambered wings. In insect flight, although wing camber is obviously important in determining the aerodynamic forces on insect wings in steady flows at low angle of incidence (Rees 1975; ,Okamoto *et al.* 1996; ,Kesel 2000), its effects in unsteady flows at high angle of incidence have proven more controversial (,Dickinson & Götz 1993; ,Sunada *et al.* 1993; ,Usherwood & Ellington 2002; ,Wang *et al.* 2003). ,Figures 6 and ,7 indicate that camber is pronounced with a complex spatial distribution. Although this does not in itself imply that camber is important in unsteady flows at high angle of incidence, it indicates that most aerodynamic models of flapping insect flight, which are currently based on flat-plate kinematics (mechanical models: ,Dickinson & Götz 1993, ,Dickinson *et al.* 1999, Sane & Dickinson ,2001, ,2002, ,Lehmann *et al.* 2005, ,Lehmann & Pick 2007, ,Lu & Shen 2008; analytical models: ,Ansari *et al.* 2006*a*,,*b*; numerical models: ,Liu *et al.* 1998, ,Ramamurti & Sandberg 2002, ,Wu & Sun 2004, ,Wang & Sun 2005, ,Ramamurti & Sandberg 2007, but see ,Wang *et al.* (2003) for a notable exception), only roughly approximate the true kinematics. Insect wings are complex inhomogeneous structures built from anisotropic materials (,Smith *et al.* 2000), which deform owing to a combination of muscular, elastic, inertial and aerodynamic forcing (,Miyan & Ewing 1985; ,Ennos 1988; ,Wootton *et al.* 2000). For this reason, although mechanical models with membrane wings are becoming more common in the flapping flight literature, none has yet come close to simulating the intricacies of a variable twist and camber distribution in a real insect wing. However, there is no reason at all why current numerical methods should not simulate the effects of variable camber and/or twist (e.g. insects: ,Wang *et al.* 2003; fishes: ,Bozkurttas *et al.* 2006, ,Mittal *et al.* 2006).

Besides providing quantitative measurements of gross surface topography, the wing kinematics we describe here are of sufficiently high spatial resolution to investigate the motion of specific smart structures in the wings. For example, in the case of the locust, it is possible to measure how the surface area of the hindwing changes as a result of the ‘umbrella effect’ (Wootton 1995; ,Wootton *et al.* 2000). It is also possible to measure changes in the angles of the medial and vannal folds of the forewing, thereby quantifying changes to the z-profile which the forewing adopts during the upstroke (,Jensen 1956). In the case of the hoverfly, the most prominent smart structure is the alula, which is a simple hinged flap at the base of the wing. The movements of this structure do not appear to have been described even qualitatively to date, and although the alula is most often parallel to the plane of the wing, it displays marked kinematic variations, flipping to become almost perpendicular to the wing in some sequences (,figures 5 and ,7). We hypothesize that the alula functions as a flow control device, and will test this hypothesis separately in future work. A full discussion of the kinematic results of this work will be provided elsewhere in separate papers linking the kinematics of wing deformation to the dynamics of force production for the complete datasets for hoverflies and locusts.

### 4.2 Accuracy of photogrammetric and other methods

The only other kinds of technique that have been able to measure both the camber and twist of a flapping insect wing to date are projected laser line methods (Zeng *et al.* 2000; ,Song *et al.* 2001; ,Sunada *et al.* 2002; ,Wang *et al.* 2003). Error analyses of projected laser line methods indicate comparable percentage errors with those given here for real insect wings (,Zeng *et al.* 1996; ,Song *et al.* 2001; ,Wang *et al.* 2003), although the errors vary systematically with wing orientation, which is not the case using photogrammetric techniques. However, despite their potential accuracy, projected laser line methods have a number of technical limitations that make them less well suited to reconstructing insect wing kinematics than the photogrammetric method we present here. Projected laser line methods are limited in their spatial resolution by the number of laser lines that can be projected onto the wing, typically giving camber and angle of incidence at only five spanwise stations. Furthermore, if the wing long axis is close to being parallel or perpendicular to the projected laser fringe, then a second laser source must be used to enable measurements of wing profile to be made (,Zeng *et al.* 2000). This is clearly a problem for insects such as hoverflies with large stroke amplitudes. Multiple-camera photogrammetry avoids this problem because clear views of the wings are not needed in all camera views simultaneously. Finally, laser projection methods provide instantaneous measurements of the camber and the angle of incidence of wing sections in planes that are fixed in space, and they cannot therefore be used to track the kinematics of points or sections fixed with respect to the wing surface.

Photogrammetric methods are not in themselves new to biomechanics. However, most recent photogrammetric methods used in biomechanics have adopted a simpler method than the bundle adjustment technique we present here, known as DLT (e.g. Warrick & Dial 1998; ,Hedrick *et al.* 2002; ,Hsieh 2003; ,Socha *et al.* 2005; Bishop ,2006, ,2007; ,Tobalske *et al.* 2007). This involves linearizing equations ,(2.1) and ,(2.2) in 11 parameters that can be solved using classical linear least squares. DLT has the advantage that it is simple and efficient to implement, but the linear least-squares solution minimizes an algebraic distance as its cost function, so does not necessarily result in solutions that minimize a meaningful and statistically well-behaved geometric distance, such as the RPE (,Hartley & Zisserman 2004). DLT is therefore frequently used to provide initial estimates of the camera parameters before bundle adjustment is performed to ‘fine-tune’ the calibration (,Hartley & Zisserman 2004; ,Reyes & Bayro-Corrochano 2006). A more serious practical limitation of DLT is that it requires precise knowledge of the three-dimensional position of a number of calibration points whose coordinates in the object coordinate system are assumed to be known relative to each other without error, but, in practice, this assumption is never strictly fulfilled. Furthermore, it can often be experimentally impractical to provide a three-dimensional calibration object covering the entire viewable space, especially when working in the field or in constrained laboratory set-ups. By contrast, the bundle adjustment technique we use here simply requires several images of a two-dimensional calibration grid of known internal dimensions held in arbitrary positions and orientations throughout the object space, which greatly reduces the effort required to calibrate the cameras.

Although the bundle adjustment techniques we have applied to insect flight are standard in photogrammetry, we have also presented a suite of other tools, including the use of a template wingbeat to solve problems of data association in three-dimensional tracking. Similar tracking tools could be used in any problem involving approximately periodic kinematics, such as often arise in biomechanics. The epipolar geometry technique that we have used to enable the matching of points along curvilinear lines without any otherwise identifiable features could also be used widely in biomechanics to measure the kinematics of moving edges of other kinds besides the outline of an insect's wing (e.g. striped or patterned surfaces of animals in motion). We hope that such techniques will now become more widely adopted by biomechanists beyond the clinical and veterinary gait analysis community (see §1), in order to provide the accurate high-resolution three-dimensional kinematics that are so essential to solving many of the outstanding problems in the field. To this end, Matlab-executable software implementing the bundle adjustment used in this study is available for downloading on request.

## Acknowledgments

This research was funded by the EPSRC and MoD under grant GR/S23049/01 to A.L.R.T. and G.K.T. G.K.T. is a Royal Society University Research Fellow and RCUK Academic Fellow. We thank Rafał Żbikowski, Andrew Moore, Kevin Knowles, Nicholas Lawson, Brian White, Iain Wallace and Salman Ansari for their contributions to collaborative aspects of this research project. We especially thank Iain Wallace for his assistance with some of the experimental work. We are also grateful to Charles Bibby for his advice on photogrammetric aspects of the work. We thank the EPSRC equipment loan pool and Nicholas Lawson for lending cameras and lighting.

## Footnotes

- Received June 6, 2008.
- Accepted July 1, 2008.

- © 2008 The Royal Society