High-speed surface reconstruction of a flying bird using structured light

All flying animals rely to some extent on dynamic body-shape changes to propel themselves. Insects rely predominantly on passive wing morphing through aero-elastic wing deformation (Combes and Daniel, 2003; Wootton, 1992). Bats can actively change the shape of their wings through musculoskeletal control and muscle fibers in their membrane (Cheney et al., 2014). Amongst active flyers, birds can morph their wings to the greatest extent, from fully extended to completely folded in flight (Pennycuick, 2008; Williams and Biewener, 2015), but how they utilize morphing during flapping and maneuvering flight is not fully understood. Such questions have traditionally been addressed by measuring the 3D body kinematics of flying animals using semi-automated marker tracking (Hedrick, 2008; Hedrick et al., 2002; Ros et al., 2015; Tobalske et al., 2007), feature tracking (Biesel et al., 1985; Carruthers et al., 2010; Walker et al., 2009) or visual-hull based reconstruction methods (Fontaine et al., 2009; Muijres et al., 2014; Ristroph et al., 2009). However, none of these methods can directly and automatically reconstruct the wing surface at high resolution.

Structured-light-based methods record the deformation of a projected light pattern due to the animal's surface geometry for offline 3D reconstruction (Fig. 1A), generally by using one of two different pattern encoding techniques. Temporally encoded projection patterns require comparison of consecutive frames. Previous studies have shown that slowly moving human body parts and internal organs can be reconstructed using binary (Ackerman et al., 2002; McKeon and Flynn, 2010) and phase-shifted (Lohry and Zhang, 2014; Wang et al., 2013) temporal coding. During pilot experiments, we determined that this method is too slow to be automated for bird flight. Spatially encoded projection patterns can reconstruct a sequence of stand-alone frames and are hence called ‘single-shot’ (Salvi et al., 2010; Zhang, 2012), which gives the advantage of being robust to inter-frame movement. Some existing spatially encoded structured-light methods rely on binary pseudo-random dots but either have relatively low frame rate and accuracy (Saberioon and Cisar, 2016; Sarbolandi et al., 2015) or require manual digitizing of numerous points per frame (Wolf and Konrath, 2015; Zhang et al., 2008). Other existing spatial methods use grayscale patterns which cannot be projected at high frame rates (Guan et al., 2003; Lenar et al., 2013; Sagawa et al., 2012; Su and Liu, 2006). Because we found that no existing system can automatically measure dynamic shape changes at sufficiently high speeds, we developed a custom method. This new single-shot structured-light technique can automatically resolve body shape changes at high temporal and spatial resolution.

The experimental setup (Fig. 1A) consisted of a 3D calibrated and synchronized high-speed camera (Phantom Miro M310; Vision Research, Wayne, NJ, USA) and high-speed projector (DLP^® LightCrafter™ E4500MKII™; EKB Technologies, Bat-Yam, Israel) operating at 3200 frames s⁻¹. Calibration of the system was achieved using a modified version of the camera calibration toolbox for MATLAB (http://www.vision.caltech.edu/bouguetj/calib_doc/). All data processing was conducted in MATLAB R2015b. We analyzed the first two wingbeats after take-off of a 4-year-old female near-white Pacific parrotlet [Forpus coelestis (Lesson 1847), 27–29 g, 0.23 m wingspan], which was trained using positive reinforcement to fly between perches 0.5 m apart. All experiments were in accordance with Stanford University's Institutional Animal Care and Use Committee.

To achieve a binary single-shot light pattern for high-speed surface reconstruction, we modified a single-shot structured-light technique (Kawasaki et al., 2008). The projected pattern consists of horizontal and vertical stripes that form a grid, and in the original method, different colors were used to distinguish the horizontal and vertical stripes. We simplified this approach for use at high speed by removing the color coding, instead relying on image processing to distinguish binary stripes (Fig. 1B). Vertical stripes are equally spaced and densely packed for high spatial resolution, while horizontal stripes are unequally spaced to ensure local uniqueness. To enhance robustness, the unequally spaced horizontal stripes can be spaced such that the spacing between every set of four consecutive stripes is unique.

There are four key advantages of using this scheme. First, it is designed for full automation, which allows for high throughput of data. Second, it is single shot, which is robust for rapidly deforming objects. Third, the scheme uses binary light patterns, which allow the projectors to use their maximum frame rate. Fourth, it uses a single color that allows maximum frame rate and multiple view angles. Interference can be avoided by using different color channels, light polarization or slightly out-of-phase light pulses.

Before 3D reconstructing the surface, image processing was required to separate the vertical and horizontal stripes in the camera image, order these stripes and find their intersections. We applied the following fully automated steps (Fig. S1). The camera image was rotated to align the equally spaced stripes vertically. Next, the Laplacian of a directional Gaussian filter was applied in the horizontal and vertical directions. Adaptive thresholding was used to generate a noisy approximation of the horizontal and vertical stripes. The noise was filtered out by adding a minimum length requirement for each connected white region. Extension and closure of stripes with gaps was accomplished by choosing paths that best combine attributes of high pixel brightness and correct stripe direction. These factors were weighted, and stripes were only extended if a preset cut-off value was satisfied.

After all stripes were identified, their subpixel center was determined using the original image by quadratically fitting brightness levels perpendicular to the stripe, based on which intersection points between horizontal and vertical stripes were located. Regions near intersections produced inaccurate center lines, so these regions were interpolated and the intersections recomputed. Finally, all stripes were ordered based on connections between stripes, and discontinuous grids were ordered separately.

Brackets used in Eqns 4–7 indicate the selection of a specific [Row, Column] value in matrices.

For each potential match between known and unknown planes, the error was computed as the mean squared angle between the known and unknown planes. The correct matching sequence for the horizontal planes gave a much lower error than other possible matches due to the unequal spacing of the stripes. When the correct matching planes were found for the horizontal planes, the value of p used in Eqn 12 was then used to match the vertical planes as well.

After the unknown planes were matched with the projected planes, 3D reconstruction of the surface was straightforward. For each stripe seen on the bird, its light plane was defined. Additionally, for each point at the center of a stripe, the 3D vector along which that point must lie was specified. The intersection of the vector and plane lies on the bird's surface. We then fit a surface (average 26,907 points) to the point cloud of data (average 285 intersections and 11,405 total points) using the Gridfit toolbox in MATLAB (https://www.mathworks.com/matlabcentral/fileexchange/8998-surface-fitting-using-gridfit), which uses a modified ridge estimator with tunable smoothing. The result is shown in Fig. 2E for different wingbeat phases reconstructed with a single camera and projector.

The reconstruction accuracy achievable was estimated in a separate test with similar equipment settings. Using a sphere of known radius (22.23±0.05 mm, mean±s.d.), we found an accuracy of 0.31 mm (error) and precision of ±1.03 mm (s.d.) (see Fig. 2A,B). Errors were largest (double) in areas outside the calibration volume (the image corner regions). Additionally, occasional stripe mismatching occurred in the image processing steps, which accounts for other larger errors (Fig. S3). When processing both the sphere and bird, no manual manipulation was used and bird reconstruction was successful for 98% of the frames over four separate downstroke segments (two wingbeats of two flights).

To calculate bird-specific surface parameters, a reference frame must be defined for each frame. To accomplish this automatically, we first identified the body of the bird by identifying surface points that move less than a preset threshold distance (see Fig. S2). To compute the z-axis, the body points were fit to a plane, after which the x-axis was computed by finding the line of symmetry of the body points projected on that plane. To find a repeatable origin point, the top of the bird's head was found by fitting a paraboloid to points on the head. For the frames we analyze here, the orientation of the axis did not change significantly and was thus set constant for all frames, while the origin point was fit linearly over all relevant frames. This computed body reference frame was labeled with subscript b. Another reference frame, the wing reference frame, was used for measuring wing shape and was labeled with subscript w. It was found by rotating the body reference frame about the x_b-axis in order to best fit the (right) wing of the bird. The reference frames and the corresponding surfaces are shown in Fig. 3A,B.

Using these reference frames and the 3D data, surface metrics were computed. In Fig. 3D, the shape of the bird at its midline is tracked while in Fig. 3E,F, the shapes of dorsal airfoil slices of the wing are tracked at different spanwise positions. We determined the angles of attack spanwise along the wing (see Fig. 3G–I) based on these airfoils. The geometric angle of attack was found with a linear fit to each airfoil, while the induced angle of attack was found by computing the direction of the velocity of the wing from root to tip using the bird velocity and angular velocity of the wing. To reduce noise in these linear fits, outliers were detected using the RANSAC method (Fischler and Bolles, 1981). The effective angle of attack is the difference of the induced and geometric angles of attack. Because of the angle of the bird's wing relative to the camera and projector positions, angle of attack measurements beyond a spanwise position halfway to the wingtip are less reliable. This could be resolved in future setups by adding cameras and projectors to better image the wing under these angles.

We quantified and analyzed the 3D wing and body shape of a flying bird during a portion of four downstrokes (41% to 65% on the first downstroke and 36% to 64% on the second downstroke after take-off) using a novel high-speed, automated, 3D structured-light method. While our results are for a single projector and camera combination imaging the dorsal side of the bird, more cameras and projectors can be added to obtain full body volumetric reconstruction and provide a more complete analysis of wing morphing over an entire wingbeat cycle. In our analysis of the dorsal data, we found that the bird's tail rotates down significantly with respect to the body (1720 deg s⁻¹ in flight 1, 1550 deg s⁻¹ in flight 2) during the first downstroke after take-off (Fig. 3D) but not the second (−290 deg s⁻¹ in flight 1, 270 deg s⁻¹ in flight 2) (Fig. S3; all angular velocities are averaged over the downstroke phase analyzed). The wings rotate down at an average of 5700 deg s⁻¹ in this same portion of the downstroke. Further, the wings are tracked at different spanwise positions (Fig. 3E,F), and we see that the wing twists relative to the body through the downstroke as confirmed by the geometric angle of attack plot (Fig. 3G). Using these data, we computed the effective aerodynamic angle of attack (Fig. 3I), which remains relatively constant between 50 and 70 deg in the first downstroke and 45 and 60 deg in the second downstroke (Fig. S3C,F). These high angles during the downstroke at take-off enable the bird to support its body weight with both lift and drag, while simultaneously generating significant thrust by tilting lift forward. This is facilitated by the close to horizontal orientation of the wing surface moving predominantly downward (Fig. 3A,B) combined with the high lift and drag coefficients of a bird wing at these high angles of attack (Kruyt et al., 2014). The measurements illustrate how the system provides insight into how birds morph their wings to generate aerodynamic forces, and could give insight into maneuvering flight in future studies.

Beyond studying birds, this new structured-light method has many benefits for capturing 3D data in experimental biology and engineering. It is automated to allow for a high throughput of data, single-shot to track deforming objects, binary to allow for high-speed tracking, and uses a single color to allow for multiple view angles. While the current application of this technique is for studying bird flight, it is broadly adaptable for tracking shapes of other dynamically morphing animals, plants and objects of interest.

We thank J. W. Kruyt for initially helping us develop and test temporally encoded structured-light methods, and A. K. Stowers for reviewing the mathematics.