A comparison of Fourier methods for the description of wing shape in mosquitoes (Diptera: Culicidae)

A Comparison of Fourier Methods for the Description of Wing Shape in Mosquitoes (Diptera: Culicidae) ® F. James Rohlf; James W. Archie Systematic Zoology, Vol. 33, No. 3 (Sep., 1984), 302-317. Stable URL: http://links.jstor.org/sici?sici=0039-7989%28198409%2933%3A3%3C302%3AACOFMF%3E2.0.CO%3B2-T Systematic Zoology is currently published by Society of Systematic Biologists. Your use of the JSTOR archive indicates your acceptance of JSTOR's Terms and Conditions of Use, available at http://www.jstor.org/about/terms.html. JSTOR's Terms and Conditions of Use provides, in part, that unless you have obtained prior permission, you may not download an entire issue of a journal or multiple copies of articles, and you may use content in the JSTOR archive only for your personal, non-commercial use. Please contact the publisher regarding any further use of this work. Publisher contact information may be obtained at http://www.jstor.org/journals/ssbiol.html. Each copy of any part of a JSTOR transmission must contain the same copyright notice that appears on the screen or printed page of such transmission. JSTOR is an independent not-for-profit organization dedicated to creating and preserving a digital archive of scholarly journals. For more information regarding JSTOR, please contact support@jstor.org. http://www.j stor.org/ FriJul 30 16:15:31 2004 Syst. ZooL, 33(3):302-317/ 1984 A COMPARISON OF FOURIER METHODS FOR THE DESCRIPTION OF WING SHAPE IN MOSQUITOES (DIPTERA: CULICIDAE) F. James Rohlf1 and James W. Archie2 department of Ecology and Evolution, State University of New York, Stony Brook, New York 11794; and ^Department of Zoology, University of Hawaii, Honolulu, Hawaii 96822 Abstract.—Outlines of wings from 127 species of North American mosquitoes were digitized. Comparisons were made among several different methods of reducing the information in the resulting coordinates to a series of descriptors that could be used in multivariate analysis. Methods included Fourier analysis of both radii and tangent angle change functions. In addi- tion, the relatively new method of elliptic Fourier analysis was tried. Cluster and ordination analyses based on the various sets of descriptors summarized well the pattern of similarities and differences in wing shapes, but clusters of similar wings do not agree well with traditional taxonomic groupings. The use of elliptic Fourier descriptors appears to be especially promising for future work. [Fourier analysis; morphometries; mosquitoes; image analysis; feature extrac- tion.] The present study is concerned with methods for quantitative analysis of simi- larities and differences in wing shape among 127 species of mosquitoes from America north of Mexico (see Fig. 1 for examples of the shapes of wings found among these species). There are many methods for quantifying outline shapes (a feature extraction problem in image anal- ysis; Rohlf and Ferson, 1983) so that their variation can be analyzed using tech- niques of multivariate analysis. The con- ventional morphometric approach would be to define a series of characters—such as length to width ratios, and more complex ad hoc indices—to describe the ways in which wings differ in shape. This is less likely to be successful in cases like the present one where the diversity of shapes is not large and, hence, the morphometric trends are apt to be rather subtle. We are interested here in methods that allow one to capture the entire shape outline in a systematic manner and with a desired de- gree of precision. We have investigated the relative usefulness of five different feature extraction methods based on Fourier anal- yses of various functions of the coordi- nates of points around the outline of each wing. Kaesler and Waters' (1972) analysis was one of the first applications of Fourier de- scriptors to study morphological shapes in systematics. Since that time there have been many such applications (some of them are listed in Rohlf and Ferson, 1983). The present study represents an attempt to define and evaluate new types of taxo- nomic characters that can be recorded au- tomatically using optical scanners on dig- ital computers and, thus, make larger quantitative studies more feasible. Meth- ods based on Fourier coefficients are not the only ones that can be used. Another approach is to use moment invariants (Hu, 1962). In this method an image is treated as if it were a bivariate density function (image darkness treated as proportional to probability density) and then the x and y moments are computed. Functions of these moments are computed that are invariant to differences in image location, size, ro- tation, reflection, and contrast. One limi- tation of this approach is that one cannot easily reconstruct an image from the de- scriptors. Empirical evaluations of their usefulness in biological image analysis will be given by Ferson et al. (in prep.) and Rohlf and Ferson (in prep.). Rohlf and So- kal (1967) described another method in which the presence or absence of a visible feature at a series of randomly or system- 1984 FOURIER METHODS AND WING SHAPES 303 atically selected locations were used as characters. This method requires careful alignment of images and does not permit an accurate reconstruction of the image. While this method seems to "work," it is not very satisfying intellectually. There have been few numerical taxo- nomic studies of mosquitoes. Rohlf (1963) studied 48 species of Aedes based upon characters from both the adults and larvae, while Steward (1968) reported on a similar analysis using 82 characters for 42 Cana- dian species of this genus. Hendrickson and Sokal (1968) studied 29 species from the genera Psorophora and Aedes based on 158 adult morphological characters. Rohlf (1967) reported on a classification of 45 species in several genera using only pupal characters. The images of these same pu- pae were also studied in Rohlf and Sokal (1967). Rohlf (1977) analyzed data on 63 species of Aedes using 14 adult thoracic se- tal characters. Moss et al. (1979) studied 25 species of the genus Toxorhynchites from the Orient using 79 adult and 26 pupal characters. A follow-up study (which in- cluded a numerical cladistic analysis) by Simon et al. (1982) included 32 species of Toxorhynchites based on 100 characters from the immature stages. While general agree- ment has been found with previous clas- sifications, there have also been a number of important differences (Crovello, 1969; Nielsen, 1969). Much further work is clearly needed both at the higher and low- er taxonomic levels. MATERIALS AND METHODS Since this study is concerned with the possible usefulness of the wing outline to determine broad relationships within the family as a whole (rather than to distin- guish between similar species), the data were obtained from plates 1 to 127 in Car- penter and LaCasse (1955). These corre- spond to the 127 species of mosquitoes from America north of Mexico included in that volume. The use of these published illustrations made it possible for this study to be comprehensive without the compli- cations of mixing images from different authors or from microscope slides pre- AN3 T012 WY14 PS42 PS43 AE46 AE8 1 CX1 16 Fig. 1. Examples of mosquito wing shapes based on the RP data set. The central "dot" corresponds to the centroid of each wing and the "dot" to its left corresponds to the point at which the fifth longitu- dinal vein branches. Codes explained in Table 1. pared using different procedures. (Studies concerned with the problem of distin- guishing between pairs of similar species are underway using actual wings from many replicate specimens digitized auto- matically from TV images. These will be reported separately.) In various diagrams in this paper, the species are identified by a two-letter designation for each genus (see Table 1) followed by a number that cor- responds to the plate number in Carpenter and LaCasse (1955). The original data were obtained using a coordinate digitizer connected to an HP9830 desk-top computer that was con- trolled by a program written in BASIC. The digitizer is accurate to 0.01 inches with a 304 SYSTEMATIC ZOOLOGY VOL. 33 Table 1. List of genera with codes and species code numbers. Numbers Anopheles AN 1-11 Toxorhynchites TO 12 Wyeomyia WY 13-16 Uranotaenia UR 17-20 Culiseta CA 21-27 Orthopodomyia OR 28 Mansonia MA 29-31 Psorophora PS 32-43 Aedes AE 44-100 Culex CX 101-125 Deinocerites DE 126-127 field of 17 by 17 inches (effectively a 1,700 by 1,700 pixel grid). The outline of a right wing of each species was traced using a cursor and the x- and y-coordinates were read "continuously" by the computer (i.e., the digitizer furnished coordinates as rap- idly as the computer could request them). The wings were oriented more or less hor- izontally and the point at which the fifth longitudinal vein branches was taken as the origin (this landmark was used, since it is centrally located on most wings). These initial coordinates were converted to polar-coordinate form as the points were digitized. Angles were measured in a clockwise direction starting from the point at which the leading edge of the wing joins the thorax. In order to eliminate the effect of unequal sampling due to the variation in speed at which different parts of the wing were traced, coordinates were com- puted for the end points of 100 equally- spaced radii (trigonometric interpolation was used to give equal angles). Thus, the basic data set was obtained in radius-func- tion form (it consisted of 100 distances from the branching point of the fifth hor- izontal vein to the edge of the wing, every 3.6°) for each of 127 species of mosquitoes (see Fig. 2A). This was repeated five times for each wing (to reduce operator error in tracing) and the averaged lengths of the radii were transmitted to the campus' UNIVAC 11/82 computer and the depart- ment's MODCOMP IV/25 and IV/35 min- icomputers for further processing. These radii (see Fig. 2B) were analyzed by Fou- rier analysis. The coefficients for the first 15 harmonics were computed for each of the 127 wings and used as descriptors. Since the "zeroth" harmonic is just a func- tion of the size of the wing, all coefficients were scaled so that the zeroth harmonic was equal to 1.0. The resulting matrix of descriptors for harmonics 1 through 15 is referred to below as the "RP" (raw polar) data set. A second data set was constructed math- ematically from the first by translating the origin of the system of polar coordinates to the centroid (center of gravity) of each wing (see Fig. 2C). The location of the cen- troid was determined by a numerical in- tegration. Anstey and Delmet (1973) in- dicated that this choice of an origin was preferable to that of a morphological land- mark since its use increases the relative amount of information in the harmonics above the first. Actually, the use of the centroid simply ensures that the results are just a function of the wing shape and not of changes in the location of the landmark within a wing (this simple adjustment seems to have been overlooked by Roberts et al., 1983:379). This proved to be impor- tant since species were found that had very similar wing outlines but differed in the location of the branching point of the fifth horizontal vein. A new set of points were computed by interpolation so that the ra- dii would be equally-spaced (every 3.6° around the new origin). A matrix of Fourier coefficients for har- monics 0 through 15 was calculated from these data and is denoted below as the "CP" (centered polar) data set. One of the major advantages of these various types of Fourier coefficients as descriptors is that one can reconstruct the outlines of the wings both from the various data matrices and from cluster's centroids and points along principal component axes. This is useful as a check on the numerical results (and as an aid to interpretation). Figure 3A shows the separate contribution of har- monics 0 through 8 to the average mos- quito wing (the wing one obtains if one uses the average Fourier coefficients based 1984 FOURIER METHODS AND WING SHAPES 305 Fig. 2. Representations of the coordinate systems used for the data sets in the present study. (A) Polar coordinates relative to a landmark. (B) The radius from A shown as a function of the angle. (C) Polar coor- dinates relative to the centroid of the wing. (D) Angle of a tangent vector as a function of distance along the perimeter of the wing. The elliptic Fourier coefficients are based simply on the x- and y-coordinates of points along the wing outline and are not shown. on all 127 species), and Figure 3B provides the successive reconstructions of the wing outline based on harmonics 0 through 8 and 15. Since the zeroth harmonic is again just a function of the size of the wing, all coefficients were scaled by dividing them by the square root of the area of the wing. While it is not important when consid- ering wing outlines (which are simple convex images), the two methods de- scribed above are limited to images in which the radius functions are single-val- ued. These methods cannot be used for shapes that have many concave curves (such as maple leaves). Zahn and Roskies (1972) suggested a more general way to describe the shape of an object based on the cumulative change in angle of a vector tangent to the outline of the object as a function of distance around the periphery of the wing (see Fig. 2D). The wing is first scaled so that its perimeter is 2ir. Then the function 0(0 = d(t) - 6(0) - t (1) is computed for 100 equally-spaced values of t. In this function, t is the distance along the periphery of the wing from the start- ing point (ranges from 0 to 2w radians), 0(0) is the angle of a tangent at the starting point, and 6(t) is the angle of a tangent vector at a distance t from the starting point. Bookstein et al. (1982) referred to this as an intrinsic representation. The third, "TA" (tangent angle), data set con- sists of the coefficients for harmonics 0 through 15 based on a Fourier analysis of the function 0(f) described above. Figure 4 shows the contributions of the first eight harmonics to the description of the aver- 306 SYSTEMATIC ZOOLOGY VOL. 33 A B Fig. 3. Fourier harmonics for the average mos- quito wing based on polar coordinates relative to the centroid of the wing. (A) Contributions of the first eight harmonics. The relative sizes of the figures re- flect the magnitudes of their individual contribu- tions. (B) Reconstructions of the wing outline based on the cumulative contributions of the first eight har- monics and for the first 15 harmonics. age wing (defined as above) and its recon- struction using the sums of the contribu- tions of various harmonics. The Fourier analyses applied above can be thought of simply as multiple regres- sion analyses in which the best fitting val- ues (in a least squares sense) of the coef- ficients at and b, are found in the following trigonometric regression equation, A B 2 O CT3 3 O <=> Fig. 4. Fourier harmonics for the average mos- quito wing based on angle of a tangent vector as a function of distance along the perimeter of the wing. (A) Contributions of the first eight harmonics. Note that the contribution for the first harmonic does not result in a closed curve. (B) Reconstructions of the wing outline based on the cumulative contributions of the first eight harmonics and for the first 15 har- monics. The curves do not close due to the effects of the first harmonic. k y = 0o + 2 (0,cos[#] z=l + fe,sin[if]). (2) In this equation the angle t varies from 0 to 27r. The zeroth harmonic describes the contribution of a centered circle, the first harmonic an offset circle, the second a 1984 FOURIER METHODS AND WING SHAPES 307 "figure 8," the third a trefoil, the fourth a quatrefoil, etc. (see Fig. 3A). The contri- butions are less obvious for the TA data set (see Fig. 4A), since the function being fitted, 0(f), is not just a simple radius. The coefficients can be computed by a variety of computational algorithms. The algo- rithm given in Ralston (1965) was used for data sets RP, CP, and TA described above (this algorithm takes advantage of the fact that the data points were equally-spaced). Since there are relatively few data points, there is no great advantage in using the fast Fourier transformation type algo- rithms for any of the data sets we used. The fourth, fifth, and sixth data sets were constructed using the x- and y-coordinates (rectilinear) of each data point in the orig- inal basic data matrix. The 'TP" (elliptic Fourier) data matrix consists of coefficients of the elliptic Fourier series of the gener- alized chain-encoded outline of a figure as suggested by Kuhl and Giardina (1982). Their algorithm does not require the points to be spaced equally and can fit an arbi- trary closed contour—given enough har- monics. The coefficients of the nth har- monic in a Fourier series expansion for the x-projections of an outline are A = and 2n2ir2 ~ Atp •(cos[2rniVT] - cos[2mr£p - 1/T]) 2^ (3) 2«27T2 ~2 Atp ■(sin[2nirtp/T] - sin[2tt7rfp - 1/T]), (4) where k is the number of steps in the trace around the outline (indexed by p), and Axp is the displacement along the x-axis of the contour between steps p — 1 and p. Atp is the length of the linear segment between these steps, tp is the accumulated length of such segments, and T = tk is the total length of the contour as approximated by the trace polygon. The coefficients for the y-coordinates, C„ and D„, are found in the same way using the incremental changes in the y-direction. The resulting coeffi- cients are not dependent on the x or y dis- placement of the wing, and the original wings were aligned so as to have the same orientation and starting point for the trace of their outline. A related method called "dual-axis Fou- rier shape analysis," DAFSA, has been proposed by Moellering and Rayner (1981, 1982). As presented by them, this tech- nique requires that the outline must be a closed plane curve that is sampled at uni- form arc-length intervals. The x- and y-co- ordinates are treated as pairs of complex numbers (x + iy) expressed as a function of the arc-length distance of each point from an arbitrary reference point. This function is then fitted by Fourier analysis. Kincaid and Schneider (1983) seem to have used the same method. The "EPO" data matrix is the "EP" matrix with the terms for the zeroth harmonics deleted to remove the effects of these terms (the starting position for the digitization trace along the outline). Figure 5 shows the individual and the cumulative contri- butions of harmonics one through eight to the reconstruction of the average mosqui- to wing. The coefficients for the "EP" data set can be mathematically adjusted (i.e., normal- ized) to be invariant to size, rotation, and starting position of the outline trace. The use of the following matrix transforma- tion yields the "NEP" (normalized elliptic Fourier) data set. c„ dt _W cos \p sin \p E*lsin \p cos \p Cn D, cos n6 —sin n6\ sin n6 cos n6 r (5) where n is the harmonic order of the four coefficients {a, b, c, and d), and E* = {a*2 + c*2)05 (6) 308 SYSTEMATIC ZOOLOGY VOL. 33 is the magnitude of the semimajor axis of the best-fitting ellipse. The angle of rota- tion (between 0 and 2ir radians) of this el- lipse is i^ = arc tan(c*/fl*), (7) and the angle of rotation (from 0 to 7r ra- dians) of the starting point from the new starting point at the end of the ellipse (standardized arbitrarily to one end for a given set of images) is 0 = 0.5 arctan^xBi + QDJ/ [A,2 + Q2 - B> - Dffj. (8) The values of a* and c* are given by the matrix equation (?)-(£ £)(£::} After this transformation, three coeffi- cients are constant (ax = 1, bx = 0, cx = 0) and can be ignored. The last coefficient for the first harmonic, du gives the eccentric- ity of the ellipse. Each of the matrices of Fourier coeffi- cients were then used as input to several methods of multivariate analysis. The ef- fect of standardization of the characters (coefficients) was investigated by using both standardized and unstandardized data (the variances among species for each har- monic are shown in Fig. 6). The use of unstandardized data is of interest here since the original coordinates were scaled so that the wings all had the same area and the coefficients of the various harmonics were all in comparable units. An un- weighted pair group (UPGMA) cluster analysis based on a matrix of average taxo- nomic distances among the 127 species of mosquitoes was performed and cophenet- ic correlation coefficients were computed as measures of goodness of fit (for a gen- eral account, see Sneath and Sokal, 1973). Next a principal components analysis (PCA) was performed on the variance— covariance matrix for the characters and the 127 species were then projected onto A B Fig. 5. Elliptic Fourier harmonics for the average mosquito wing based on x- and y-coordinates of the wing. (A) Contributions of harmonics one through eight. (B) Reconstructions of the wing outline based on the cumulative contributions of the first 8 har- monics and for the first 20 harmonics. the resulting axes. The configuration of points obtained from the data sets were compared visually and by Gower's (1971) M2 criterion which measures the sum of the squared distances between corre- sponding points after one PCA matrix has been rotated to superimpose as well as possible on another matrix. Smaller M2 values correspond to better agreement be- tween two PCA solutions. RESULTS The effect of standardization was ex- amined first. The cophenetic correlation 1984 FOURIER METHODS AND WING SHAPES 309 0.0006 o. o o 04-; 0.0 002 v, B P71 FY! r 10 12 0.00 15- 0.0010 0.0005 0.0006 0.0004 0.0002 1 1 T"1 T"1 —1-1-1-1—1 2 4 6 8 10 12 14 1 Fig. 6. Plot of the variance among wings of the 127 species of mosquitoes for the harmonics used in this study. The variances for each coefficient have been summed to give a total contribution for each harmonic. Variances for the: (A) RP data set; (B) CP data set; (C) TA data set; (D) EF data set; and (E) NEF data set. coefficients from the UPGMA cluster anal- yses were similar whether or not stan- dardization of the variables was per- formed. The average correlation for unstandardized data was 0.705 versus 0.691 for standardized data. The average matrix correlation between the distance matrices and distances based upon three-dimen- sional principal component analyses (see Rohlf, 1972) was higher for unstandard- ized data than for standardized data (0.991 versus 0.909). The percentage of variance explained by the first three principal com- ponent axes was also higher for unstan- dardized data than it was for standardized data. Standardization increases the rela- 310 SYSTEMATIC ZOOLOGY VOL. 33 Table 2. Correlations among taxonomic distance matrices for data sets based upon unstandardized data. Data set3 1 RP 2 CP 3 TA 4 EF 5 EFO 2 CP 0.592 3 TA 0.180 0.310 4 EF 0.561 0.798 0.297 5 EFO 0.560 0.798 0.289 1.000 6 NEF 0.469 0.827 0.296 0.907 0.907 a Data set codes: RP, raw polar; CP, centroid translated polar; TA, tangent angle; EF, elliptic Fourier; EFO, EF with zeroth harmonic left out; and NEF, normalized elliptic Fourier. tive contribution of harmonics that do not explain much of the variance (e.g., the odd harmonics in Fig. 3). This increases the ef- fective dimensionality of the space so that the first three axes explain a smaller pro- portion of the now larger total variance. The relationships indicated by unstan- dardized data seem to correspond more to one's subjective visual impression of sim- ilarity. This can be seen by examining par- ticular pairs of species for which very dif- ferent relationships are indicated when the data are standardized. For example, based on unstandardized data from the RP ma- trix, species PS42 and AE81 are indicated to be very similar to each other but rather different from the other species. When the matrix was standardized, species PS42 is indicated to be very different from all oth- er species. If one superimposes the wing outlines for these two species one finds that there is a general overall agreement in the outlines. There are, however, some small differences (that would result in dif- ferent values for the higher order har- monics). Species AN3, AN11, and AE55 present another interesting example of the effect of standardization. Based on stan- dardized data AN 11 and AE55 are indicat- ed to be much more similar to each other than either is to AN3. Using unstandard- ized data one finds that AN3 and AE55 now are more similar than either is to AN11. These latter similarities agree closely with the visual impression one ob- tains by looking at the actual wing out- lines. The explanation for the unsatisfactory effects of standardization appear to be just that it gives equal weight to all harmonics and, thus, results in much greater empha- sis on differences among species in the higher order harmonics (which are much more sensitive to small irregularities in the outlines and to measurement errors). This is clearly seen in the results of principal components analyses where many more of the higher order harmonics had high loadings on the first few axes. Therefore, unstandardized data and variance-covari- ance matrices were used in all analyses re- ported below. The overall degree of similarity across the different data sets for the relationships among the 127 species is shown in Table 2. This table gives matrix correlations be- tween pairs of average taxonomic distance matrices from the various data sets. The matrices all are positively correlated, in- dicating some degree of general agree- ment. Data set TA has much lower corre- lations with the other data sets. Leaving out the zeroth elliptical harmonic appar- ently had little effect in this study since EF and EFO had a correlation of 1.000 (the two matrices are not, however, identical). The NEF distance matrix is highly corre- lated with those of the EF and EFO data sets. Data sets CP and EF indicate quite similar relationships among the 127 species since the correlation between their dis- tance matrices was 0.798. The similarity among these data sets is a measure of the consistency with which the original wings were manually aligned. Table 3 gives some comparisons among the results of UPGMA cluster analysis ap- plied to the various data sets. The TA data set resulted in the largest cophenetic cor- relation (due in part to the fact that one species was indicated to be isolated from the other species). Pairs of phenograms were compared by computing the CIC con- sensus index from a strict consensus tree (Rohlf, 1982). This very stringent criterion indicates that the phenograms are by no means identical for the different data sets. The largest value, 0.752, indicates, for ex- ample, that the EF and the EFO pheno- grams have 75.2% of their clusters repre- sented identically in both trees. The lowest values are for comparisons of data set TA 1984 FOURIER METHODS AND WING SHAPES 311 Table 3. Comparisons among UPGMA phenograms based on distance coefficients from unstandardized data. Coefficients are the CIC consensus coefficients based on strict consensus trees. Cophenetic correlation coefficients are also given for each data set. Larger values correspond to closer agreement. Data set _ Cophenetic Data set 1 RP 2 CP 3 TA 4 EF 5 EFO correlation 1 RP 0.636 2 CP 0.048 0.720 3 TA 0.088 0.008 0.786 4 EF 0.088 0.080 0.016 0.704 5 EFO 0.088 0.088 0.024 0.752 0.690 6 NEF 0.016 0.096 0.016 0.128 0.128 0.708 with the others. Figure 7 shows the phe- nograms from the CP and NEF data set (the others are not furnished in order to save space). While similar wings usually cluster together, the clusters so found do not cor- respond well to conventional taxonomic groupings. Table 4 summarizes the results of com- parisons of pairs of three-dimensional principal components solutions for the various data sets. The values given are Gower's (1971) M2 coefficient. Again, data sets EF and EFO agree closely and data set TA gives the least agreement with the oth- ers. The percentages of variance ex- plained, matrix correlations with the orig- inal distance matrices and between the principal component solutions are also furnished in Table 4. Figure 8 shows the results of the prin- cipal components analysis for the CP and the NEF data sets. In order to properly ap- preciate the relationship among the points, it is essential that the first and second axes are plotted to the same scale (this impor- tant point is often overlooked, even in re- cent papers). It is apparent from these fig- ures that, while there is little evidence of taxonomic structure in the data, similarly- shaped wings are plotted close together in the figures. There are also consistent changes in the shapes of the wings as one goes from left to right and from top to bottom across the figures. These trends can best be visualized by reconstructing wings corresponding to points exactly along the principal component axes. Figure 9 de- picts wings that would be plus or minus three standard deviations from the mean along axes I and II and zero along all other axes. This shows very clearly that axis I has narrower wings at the right and wider at the left. Axis II is harder to describe but it clearly contrasts wings that are wider near the base versus near the apex. This technique is potentially a very powerful tool since it allows one to visualize what an organism would look like if it were to be located in a particular region of a prin- cipal component (or even a canonical vari- ate) space, even if no observed points fell into that particular region. The ordination for the TA data set (not shown) is more similar to the results given by the other data sets than is suggested by the values given in Table 4. Most of the differences seemed to be due to the fact that species CXI 16 is indicated to be an extreme out- lier in the plot for the TA data set, but not for the other data sets. DISCUSSION Bookstein et al. (1982) pointed out a number of limitations on the biological in- terpretations one can give to the coeffi- cients of Fourier functions (but see also the rebuttal by Ehrlich et al., 1983). Book- stein et al. pointed out, for example, that the methods used for the RP and CP data sets are sensitive only to differences in shape not to differences in interpretation of homology between radii at different points along an outline (however, if one's goal is to measure shape per se, this can be considered an advantage). Another limitation is that a change in part of the shape (a 'Tocal change") may result in changes in the values of many of the coef- 1984 FOURIER METHODS AND WING SHAPES 313 Table 4. Comparisons of projections of the 127 species onto the first three principal component axes based upon unstandardized data. Gower's M2 coefficient is given for each pair of data sets compared. The percent- ages of variance explained and matrix correlations are also provided for each data set. Data set Data set 1 RP 2 CP 3 TA 4 EF 5 EFO explained correlation 1 RP 85.33 0.984 2 CP 0.692 88.25 0.990 3 TA 1.460 0.989 79.32 0.985 4 EF 0.706 0.469 0.968 92.04 0.998 5 EFO 0.706 0.469 0.968 0.000 92.15 0.996 6 NEF 0.878 0.466 0.975 0.288 0.288 88.60 0.993 ficients (not just in coefficients corre- sponding to the high frequency compo- nents as one might have hoped). This makes interpretation of the coefficients more difficult. However, these problems are not unique to Fourier methods. The coefficients of linear, quadratic, cubic, etc., polynomials fitted to some data need not be biologically interpretable if the under- lying biological process that generated the relationship under study corresponds to some other functional relationship (such as exponential or sinusoidal). This is an inherent problem in curve fitting when one fits functions that do not correspond to the actual underlying biological pro- cesses. Traditional morphometric ap- proaches have an analogous problem since the usual arbitrary suites of morphometric measurements need not correspond di- rectly to the underlying developmental, evolutionary factors either. Principal com- ponents and canonical variates analyses have the same type of problem in inter- pretation. The most important limitation of our approach is that we deal only with overall shape—not with changes in dis- tances between homologous points. Book- stein (1982) considered the latter to be of more fundamental importance in morpho- metries. Work in progress by one of us (FJR) uses information on the shape and position of each wing vein. The branching points of the veins and their points of in- tersection with the margin of the wing will provide homologous reference points. The fact that the individual Fourier coef- ficients will almost never be morphologi- cally meaningful might actually be con- sidered an advantage—since it may help an investigator to resist the temptation to create complex interpretations of the meaning of the individual coefficients (Rohlf and Ferson, 1983). We would not expect, for example, that a relatively high contribution of the second spherical har- monic in the average mosquito wing (see Fig. 3) would lead an investigator to con- clude that there was a growth field that had the form of a quatrefoil. The investi- gator is thereby forced to deal with overall shape differences by using suites of coef- ficients. After comparing several methods, the question as to which method is "best" nat- urally arises. All the methods are mathe- matically equivalent in that given enough harmonics they can encode the shape of the wing outlines exactly. For objects with more complex outlines there may not be a single "center" that enables the radius function to be single-valued. This elimi- nates the two polar methods from consid- eration as general methods (even though they may work well for many shapes such as those of insect wings). The method Fig. 7. UPGMA phenograms based on average taxonomic distance from unstandardized Fourier coeffi- cients. (A) CP data set. (B) NEF data set. 314 SYSTEMATIC ZOOLOGY VOL. 33 Fig. 8. Plots of the first two dimensions from principal components analyses of unstandardized Fourier coefficients. (A) CP data set. (B) NEF data set. 1984 FOURIER METHODS AND WING SHAPES 315 based on tangent angles has two serious problems. The first is that, as pointed out by Zahn and Roskies (1972), reconstruc- tions with relatively few harmonics do not usually result in figures with closed con- tours (see Fig. 4). This makes it more dif- ficult to interpret trends in reconstructed images (such as shown in Fig. 9 for the CP and NEF methods) since they may not look realistic. The second problem is that the estimation of the coefficients seems to be more sensitive to "noise" in the image than the other methods (as mentioned above species CXI 16 was considered to be an outlier only by this method). This leaves the methods based on elliptic Fourier coef- ficients as the most generally useful. The use of the normalized coefficients, NEF, is a convenience since it saves the investi- gator from having to worry about aligning the images in a standard fashion. The only disadvantage is that a small amount of ad- ditional computation is required (which may be a problem if the analyses are being performed on a small microcomputer without floating-point hardware). Our original (somewhat naive) hope was that an analysis of wing shape might yield classifications comparable to conventional taxonomic classification, but even in the absence of such congruence we feel that the present application is quite successful. The ordinations enabled us to summarize and display the diversity of shapes of wings found among the North American Culicidae—even though the wings are fairly similar, and some of the differences rather subtle and difficult to describe by conventional means. The differences in wing shape must have functional and ecological significance, but we did not find any obvious associations with published accounts indicating that a given species was a "strong" or "weak" flier, was found in exposed or sheltered habitats, etc. The detection of such asso- ciations may require additional quantita- tive information about the behavior and physiology of the various species of mos- quitoes. These results are in contrast to the study A Fig. 9. Reconstructed wing outlines for points along the principal component axes shown in Figure 8. (A) CP data set. (B) NEF data set. by Plowright and Stephen (1973) who found good agreement between their morphometric analysis of bee wings and the conventional classification. Their characters described features of the vena- tion rather than the outline shape of the wings as in the present study. Brown (1979) and Brown and Shipp (1977, 1978) also used wing vein lengths to examine phenetic relationships within two groups of Australian Diptera and found that the results varied in their consistency with re- sults from previous studies. The previous studies had used standard taxonomic char- acters including genitalic characters as well as geographic distribution for erecting species groups. Simon (1983) used a series of 49 characters involving wing vein length and was able to establish differ- ences among year classes (broods) and life cycle forms of periodical cicadas (Magici- cada) where no differences had previously been established. Whether differences in wing outline shape will provide useful taxonomic information when more di- verse wing shapes are analyzed remains to be seen. Minor wing shape differences may also be expected between closely related 316 SYSTEMATIC ZOOLOGY VOL. 33 species or populations if these differences are important for ecological / functional dif- ferences in flight characteristics. One of us (JWA) is examining both wing shape and vein length differences among popula- tions and species of Hawaiian picture- winged Drosophila where drastic wing shape differences between species and sexes are apparent. The ready availability of inexpensive coordinate digitizers attached to micro- computers (or to terminals) makes these methods practical for general use. The ac- tual process of digitizing outlines is tedi- ous, however. It takes several minutes to position an image, carefully trace its out- line, and then store the information in a computer. This process can be greatly speeded up for conveniently-sized high- contrast images (where one does not have to adjust a microscope and the background can clearly be distinguished from the ob- ject) by the use of a television camera at- tached to a computer-controlled digitizer. Such a setup in our laboratory (see Ferson et al., in prep.) enables large numbers of specimens to be processed since only a few seconds are required for each image. The principal problem is the fact that the pre- cision of the resulting measurements is somewhat limited. The scanner in our lab- oratory has a resolution of maximally 480 by 640 pixels compared to the 1,700 by 1,700 field on the digitizer used in the present study. Most video units presently available for microcomputers have much less resolution (undoubtedly only a tem- porary limitation). Charge-coupled device (CCD) cameras are also available that have a resolution of 1,700 by 2,200 (as much res- olution as a digitizer), but the processing of an image with 3.74 x 106 bytes requires much more computational effort. The de- termination of which device should be used will depend upon the complexity of the outline of the image and how many harmonics are needed to describe it satis- factorily. This will have to be determined in each study. ACKNOWLEDGMENTS This paper represents contribution number 489 from Graduate Studies in Ecology and Evolution, State University of New York at Stony Brook. This work was supported in part by a grant (BSR8202269) from the National Science Foundation. Some of the pre- liminary computations were performed while one of us (FJR) was on leave at the IBM Thomas J. Watson Research Center. Mr. Scott Ferson assisted this proj- ect in many ways. His help was greatly appreciated. Barbara Thomson assisted in the preparation of the final figures. REFERENCES Anstey, R. L., and D. A. Delmet. 1973. Fourier anal- yses of zooecial shapes in fossil tubular bryozoans. Geol. Soc. Am. Bull., 84:1753-1764. Bookstein, F. L. 1982. Foundations of morphomet- ries. Annu. Rev. Ecol. Syst., 13:451-470. Bookstein, F. L., R. E. Strauss, J. M. Humphries, B. Chernoff, R. L. Elder, and G. R. Smith. 1982. A comment on the uses of Fourier methods in sys- tematics. Syst. Zool., 31:85-92. Brown, K. R. 1979. Multivariate assessment of phe- netic relationships within the tribe Luciliini (Dip- tera: Calliphoridae). Aust. J. Zool., 27:465-477. Brown, K. R., and E. Shipp. 1977. Wing morpho- metries of Australian Luciliini (Diptera: Calliphor- idae). Aust. J. Zool., 25:765-777. Brown, K. R., and E. Shipp. 1978. Wing morpho- metry analysis of Australian Sarcophaginae (Dip- tera: Sarcophagidae). Syst. Entomol., 3:179-188. Carpenter, S. J., and W. J. LaCasse. 1955. Mosqui- toes of North America (North of Mexico). Univ. California Press, Los Angeles. Crovello, T. J. 1969. Numerical taxonomy: Its val- ue to mosquito taxonomy. Mosquito Syst., 1:63-67. Ehrlich, R., R. B. Pharr, Jr., and N. Healy-Williams. 1983. Comments on the validity of Fourier de- scriptors in systematics: A reply to Bookstein et al. Syst. Zool., 32:302-306. Gower, J. 1971. Statistical methods of comparing different multi-variate analyses of the same data. Pages 138-149 in Mathematics in the archaeologi- cal and historical sciences (F. R. Hodson, D. G. Kendall, and P. Tautu, eds.). Edinburgh Univ. Press, Edinburgh. Hendrickson, J. A., and R. R. Sokal. 1968. A nu- merical taxonomic study of the genus Psorophora (Diptera: Culicidae). Ann. Entomol. Soc. Am., 61: 385-392. Hu, M. K. 1962. Visual pattern recognition by mo- ment invariants. IRE Trans, on Information The- ory, 8:179-187. Kaesler, R. L., and J. A. Waters. 1972. Fourier anal- ysis of the ostracod margin. Bull. Geol. Soc. Am., 83:1169-1178. Kincaid, D. T., and R. B. Schneider. 1983. Quanti- fication of leaf shape with a microcomputer and Fourier transform. Can. J. Bot., 61:2333-2342. Kuhl, F. P., and C. R. Giardina. 1982. Elliptic Fou- rier features of a closed contour. Computer Graph- ics and Image Processing, 18:236-258. Moellering, H., and J. N. Rayner. 1982. The dual axis Fourier shape analysis of closed cartographic forms. Cartographic J., 19:53-59. Moellering, H., and J. N. Rayner. 1983. The har- 1984 FOURIER METHODS AND WING SHAPES 317 monic analysis of spatial shapes using dual axis Fourier shape analysis (DAFSA). Geograph. Anal- ysis, 13:64-77. Moss, W. W., W. A. Steffan, N. L. Evenhuis, and D. L. Manning. 1979. The genus Toxorhynchites (Diptera: Culicidae); Analysis of T. splendens and allies using techniques of numerical taxonomy. Mosquito Syst., 11:258-274. Nielsen, L. T. 1969. A critique on numerical tax- onomy. Mosquito Syst., 1:23-25. Plowright, R. C., and W. P. Stephen. 1973. A nu- merical taxonomic analysis of the evolutionary re- lationships of Bombus and Psithyrus (Apidae: Hy- menoptera). Can. Entomol., 105:733-743. Ralston, A. 1965. A first course in numerical anal- ysis. McGraw Hill, New York. Roberts, D. McL., A. Waren, and C. R. Curds. 1983. Morphometric analysis of outline shape applied to the peritrich genus Vorticella. Syst. Zool., 32:377- 388. Rohlf, F. J. 1963. Classification of Aedes by numer- ical taxonomic methods (Diptera: Culicidae). Ann. Entomol. Soc. Am., 56:798-804. Rohlf, F. J. 1967. Correlated characters in numeri- cal taxonomy. Syst. ZooL, 16:109-126. Rohlf, F. J. 1972. Empirical comparison of three ordination techniques in numerical taxonomy. Syst. Zool., 21:271-280. Rohlf, F. J. 1977. Classification of Aedes mosquitoes using statistical methods. Mosquito Syst., 9:372-388. Rohlf, F. J. 1982. Consensus indices for comparing classifications. Math. Biosci., 59:131-144. Rohlf, F. J., and S. Ferson. 1983. Image analysis. Pages 583-599 in Numerical taxonomy (J. Felsen- stein, ed.). NATO ASI Series G, Ecological Sciences No. 1. Springer-Verlag, New York. Rohlf, F. J., and R. R. Sokal. 1967. Taxonomic structure from randomly and systematically scanned biological images. Syst. ZooL, 16:246-260. Simon, C. M. 1983. Morphological differentiation in wing venation among broods of 13- and 17-year periodical cicadas. Evolution, 37:104-115. Simon, C. M., W. A. Steffan, W. W. Moss, and N. L. Evenhuis. 1982. The genus Toxorhynchites (Dip- tera: Culicidae); Numerical phylogenetic analysis of Tx. splendens and allies with phenetic compari- sons. Mosquito Syst., 14:221-261. Sneath, P. H. A., and R. R. Sokal. 1973. Numerical taxonomy. W. H. Freeman and Co., New York. Steward, C. C. 1968. Numerical classification of the Canadian species of the genus Aedes, Syst. Zool., 17:426-437. Zahn, C. T., and R. Z. Roskies. 1972. Fourier de- scriptors for plane closed curves. IEEE Trans, on Computers, 21:269-281. Received 7 February 1984; accepted 1 March 1984.