The function quantifies components of shape variation for a set of specimens as described by their patterns of symmetry
 and asymmetry. Here, shape variation is decomposed into variation among individuals, variation among sides (directional 
 asymmetry), and variation due to an individual x side interaction (fluctuating symmetry). These components are then 
 statistically evaluated using Procrustes ANOVA. Statistical assessment of model effects for shape variation is accomplished using permutation procedures. 
 Methods for both matching symmetry and object symmetry can be implemented. Matching symmetry is when each object contains mirrored 
 pairs of structures (e.g., right and left hands) while object symmetry is when a single object is symmetric 
 about a midline (e.g., right and left sides of human faces). Details on general approaches for the study of symmetry in geometric 
 morphometrics may be found in: Mardia et al. 2000; Klingenberg et al. 2002.
As input, the function receives either an array (p x k x n) containing raw landmarks (requiring 
 GPA to be performed) or a gpagen object (if GPA has been previously performed). If one wishes 
 to incorporate semilandmarks, GPA should be performed first using gpagen. Otherwise, 
 bilat.symmetry can perform the initial GPA, assuming all landmarks are fixed. For 
 "object.sym = FALSE, landmarks should be of dimension (p x k x 2n), as each specimen is 
 represented by both left and right configurations.
Analyses of symmetry for matched pairs of objects is implemented when object.sym=FALSE. Here, a 3D array [p x k x 2n] 
 contains the landmark coordinates for all pairs of structures (2 structures for each of n specimens). Because the two sets of 
 structures are on opposite sides, they represent mirror images, and one set must be reflected prior to the analysis to 
 allow landmark correspondence. IT IS ASSUMED THAT THE USER HAS DONE THIS PRIOR TO PERFORMING THE SYMMETRY ANALYSIS. 
 Reflecting a set of specimens may be accomplished by multiplying one coordinate dimension 
 by '-1' for these structures (either the x-, the y-, or the z-dimension). A vector containing information on individuals 
 and sides must also be supplied. Replicates of each specimen may also be included in the dataset, and when specified will be 
 used as measurement error (see Klingenberg and McIntyre 1998).
Analyses of object symmetry is implemented when object.sym=TRUE. Here, a 3D array [p x k x n] contains the landmark 
 coordinates for all n specimens. To obtain information about asymmetry, the function generates a second set of objects 
 by reflecting them about one of their coordinate axes. The landmarks across the line of symmetry are then relabeled to obtain
 landmark correspondence. The user must supply a list of landmark pairs. A vector containing information on individuals 
 must also be supplied. Replicates of each specimen may also be included in the dataset, and when specified will be 
 used as measurement error.
Notes for geomorph 3.0
Compared to older versions of geomorph, some results can be expected to be slightly different.  Starting with geomorph 3.0,
 results use only type I sums of squares (SS) with either full randomization of raw shape values or RRPP (preferred with nested terms)
 for analysis of variance (ANOVA).  Older versions used a combination of parametric and non-parametric results, as well as a combination
 of type I and type III SS.  While analytical conclusions should be consistent (i.e., "significance" of effects is the same), these
 updates maintain consistency in analytical philosophy.  This change will require longer computation time for large datasets, but the trade-off
 allows users to have more flexibility and eliminates combining disparate analytical philosophies.
Note also that significance of terms in the 
 model are found by comparing F-values for each term to those obtained via permutation.  F-ratios and df are not strictly necessary (a ratio of SS would suffice), 
 but they are reported as is standard for anova tables. Additionally, users will notice that the df reported are based on the number of observations rather than 
 a combination of objects * coordinates * dimensions, as is sometimes found in morphometric studies of symmetry. However, this change has no effect 
 on hypothesis testing, as only SS vary among permutations (df, coordinates, and dimensions are constants). 
 
The generic functions, print, summary, and plot all work with bilat.symmetry.