1.4.1 Biology: Mouse vertebrae

In an experiment to assess the effects of selection for body weight on the shape of mouse vertebrae, three groups of mice were obtained: Large, Small and Control. The Large group contains mice selected at each generation according to large body weight, the Small group was selected for small body weight and the Control group contains unselected mice. The bones form part of a much larger study and these bones are from replicate E of the study (Falconer 1973; Truslove 1976; Johnson et al. 1985, 1988; Mardia and Dryden 1989b).

We consider the second thoracic vertebra T2. There are 30 Control, 23 Large and 23 Small bones. The aims are to assess whether there is a difference in size and shape between the three groups and to provide descriptions of any differences. Each vertebra was placed under a microscope and digitized using a video camera to give a grey level image, see Figure 1.4. The outline of the bone is then extracted using standard image processing techniques (for further details see Johnson et al. 1985) to give a stream of about 300 coordinates around the outline. Six landmarks were taken from the outline using a semi-automatic procedure described by Mardia (1989a) and Dryden (1989, Chapter 5), where an approximate curvature function of the smoothed outline is derived and the mathematical landmarks are placed at points of extreme curvature as measured by this function. In Figure 1.6 we see the six landmarks and also in between each pair of landmarks, nine equally spaced pseudo-landmarks are placed.

The dataset is available in the R package shapes and the three groups can be accessed by typing:

library(shapes)
data(mice)

The dataset is stored as a list with three components: mice$x is an array of the coordinates in two dimensions of the six landmarks for each bone; mice$group is a vector of group labels; and mice$outlines is an array of 60 points on each outline in two dimensions, containing the landmarks and pseudo-landmarks. To print the k × m × n array of landmarks in R, type mice$x and to print the group labels type mice$group (‘c’for Control, ‘l’ for Large and ‘s’ for Small). In order to plot the landmark data we can use:

par(mfrow=c(1,3))
joins<-c(1,6,2:5,1)
plotshapes(mice$x[,,mice$group=="c"],joinline=joins)
title("Control")
plotshapes(mice$x[,,mice$group=="l"],joinline=joins)
title("Large")
plotshapes(mice$x[,,mice$group=="s"],joinline=joins)
title("Small")

Here the plotshapes function plots 2D (x, y) coordinates of each object, and lines are drawn between the landmarks given in the joinline option. Here lines are drawn from landmark 1 to 6 to 2 to 3 to 4 to 5 and finally back to 1 on each object. The plot is given in Figure 1.7.

In order to plot the outline data we can use:

par(mfrow=c(1,3))
joins<-c(1:60,1)
plotshapes(mice$outlines[,,mice$group=="c"],joinline=joins,col=2)
title("Control")
plotshapes(mice$outlines[,,mice$group=="l"],joinline=joins,col=2)
title("Large")
plotshapes(mice$outlines[,,mice$group=="s"],joinline=joins,col=2)
title("Small")

and the result is given in Figure 1.8. Here the points are drawn in red (col=2) and the lines are drawn to connect points 1 through 60 and back to 1.

Note that the coordinates in mice$x are also available in the shapes package individually by group: qcet2.dat, qlet2.dat, and qset2.dat, which can be useful for short-cuts in coding.

It is of interest to examine size and shape differences in the three groups, and how shape is related to size.

1.4.2 Image analysis: Postcode recognition

A random sample of handwritten British postcodes was collected and digitized (Anderson 1997), and an example digit ‘3’ is shown in Figure 1.9. It is of interest to classify each of the handwritten characters so that mail can be automatically sorted. The problem is a classic one in image analysis and many methods have been suggested, with varying degrees of success (e.g. see Hull 1990). The location and size of the characters are not so important for recognition but orientation information may be crucial, for example an ‘M’ must not be confused with a ‘W’. Some successful attempts at reading handwritten numbers include Simard et al. (1993); Hastie and Tibshirani (1994); and Hastie and Simard (1998). A survey of relevant work is given by Plamondon and Srihari (2000), and a related topic is hand-drawn gesture recognition (e.g. see Mardia et al. 1993).

Anderson (1997) obtained mathematical landmarks and pseudo-landmarks on the digital images by hand, and in particular for the digit 3 there were 13 landmarks, as shown in Figure 1.9. It is of interest to examine the average shape and variability in shape in the data, which can then be used as a prior model for digit recognition from images of handwritten postcodes.

The landmark data are given the dataset digit3.dat in the shapes package. The data are displayed in Figure 1.10 using the R command:

plotshapes(digit3.dat,joinline=1:13)

1.4.3 Biology: Macaque skulls

In an investigation into sex differences in the crania of a species of macaque Macaca fascicularis (a type of monkey), random samples of 9 male and 9 female skulls were obtained by Paul O’Higgins (Hull-York Medical School) (Dryden and Mardia 1993). A subset of seven anatomical landmarks was located on each cranium and the three-dimensional (3D) coordinates of each point were recorded.

It is of interest to assess whether there are any size and shape differences between the sexes. If there are any differences, then a description of the differences is required. An artist’s impression of the 3D skull with the anatomical landmarks is given in Figure 1.11.

The data are obtained by typing:

library(shapes)
data(macaques)

The 3D landmarks are available in the 7 × 3 × 18 dimensional array macaques$x, and the genders are in macaques$group.

We plot the data in Figure 1.12 using the command shapes3d as follows:

joins<-c(1,2,5,2,3,4,1,6,5,3,7,6,4,7)
colpts<-rep(1:7,times=18)
shapes3d(macaques$x,col=colpts,joinline=joins)

The command shapes3d uses the rgl library in R, which in turn uses OpenGL graphics. The 3D plots can be easily rotated and moved in the graphics window by clicking and moving the mouse, thus giving a good idea of the 3D geometry of the configuration.

Note that the coordinates in macaques$x are also available in the shapes package individually by group: macf.dat, and macm.dat.

1.4.4 Chemistry: Steroid molecules

Dryden et al. (2007) and Czogiel et al. (2011) analyse a dataset of steroids, which are small molecules with a wide variety of uses. The dataset consists of between 42 and 61 atoms for each of 31 steroid molecules. The three-dimensional coordinates, atom type, van der Waals radius and partial charges of each atom are given. The collection of steroids has been considered by a number of authors, including Wagener et al. (1995). This particular version of the data was constructed by Jonathan Hirst and James Melville (School of Chemistry, University of Nottingham). The steroids have different binding affinities to the corticosteroid binding globulin (CBG) receptor, and so each molecule has an activity class of either ‘1’ high, ‘2’ intermediate or ‘3’ low binding affinity. It is of interest to examine how the shape (‘steric’) properties of the molecules are related to activity class. This dataset is quite challenging in that the molecules have different numbers of atoms, and the correspondence between atoms (labelling) is not given. However, the 17 carbon atoms in the three cyclohexane rings and one cyclopentane ring are common to all the steroids, and these do correspond in a sensible way. The carbon rings are plotted in Figure 1.13 using the following commands.

data(steroids)
joins<-c(1:6,1,6,5,4,7:10,5,4,7,11:14,8,14:17,13)
shapes3d(steroids$x[,,],col=rep(1:17,times=31),joinline=joins)

1.4.5 Medicine: Schizophrenia magnetic resonance images

Bookstein (1996b) considers 13 landmarks taken on near midsagittal 2D slices from magnetic resonance (MR) brain scans of 14 control volunteers and 14 schizophrenia patients. It is of interest to study any shape differences in the brain between the two groups, either in average shape or in shape variability. If shape differences between the two groups can be established, then this should enable researchers to gain an increased understanding about the condition. In Figure 1.14 we see the 13 landmarks on a 2D slice from a scan of a schizophrenia patient. The landmarks are: 1, splenium, posteriormost point on corpus callosum; 2, genu, anteriormost point on corpus callosum; 3, top of corpus callosum, uppermost point on arch of callosum (all three landmarks registered to the diameter of the callosum); 4, top of head, a point relaxed from a standard landmark along the apparent margin of the dura; 5, tentorium of cerebellum at dura; 6, top of cerebellum; 7, tip of fourth ventricle; 8, bottom of cerebellum; 9, top of pons, anterior margin; 10, bottom of pons, anterior margin; 11, optic chiasm; 12, frontal pole, extension of a line from 1 through 2 until it intersects the dura; and 13, superior colliculus.

The data are plotted in Figure 1.15 using the following commands.

data(schizophrenia)
plotshapes(schizophrenia$x,symbol=as.integer(schizophrenia$group))

1.4.6 Medicine and law: Fetal alcohol spectrum disorder

Another important application of shape analysis is in the assessment of fetal alcohol spectrum disorders (FASDs). An MR image from the corpus callosum of a prisoner with a landmark (Rostrum) and 39 semi-landmarks is displayed in Figure 1.16. The shape of the corpus callosum has been used in court cases in expert witness testimony to help assess whether or not a defendant had been affected by FASD. Statistical shape analysis has been used successfully to help waive the death penalty for many defendants. For further details see Mardia et al. (2013a) and Section 7.10.

1.4.7 Pharmacy: DNA molecules

Molecular dynamics simulations are a widely used and powerful method of gaining an understanding of the properties of molecules, particularly biological molecules such as DNA. The simulations are undertaken with a computer package, for example AMBER (Salomon-Ferrer et al. 2013), and involve a deterministic model being specified for the molecule. The model consists of point masses (atoms) connected by springs (bonds) moving in an environment of water molecules, also treated as point masses and springs. At each time step the equations of motion are solved to provide the next position of the configuration in space. The simulations are very time-consuming to run – for example several weeks of computer time may be needed to generate a few nanoseconds of data.

We consider the statistical modelling of a specific DNA molecule configuration in water. In particular, we concentrate on the simple case of 22 phosphorous atoms, where the k = 22 atom locations are recorded in angstroms (in m = 3 dimensions) and are observed over a period of time. The temporal data are highly correlated. There are many questions of interest, for example, can we describe the main features of geometric variability, can we estimate the full configuration space of the molecule, and can we simulate the molecule over much longer time scales using fast statistical techniques. In the shape package there is a very small set of DNA data with n = 30 observations, which can be obtained using:

data(dna.dat)

A plot of the data is given in Figure 1.17.

1.4.8 Biology: Great ape skulls

In an investigation to assess the cranial differences between the sexes of ape, 29 male and 30 female adult gorillas (Gorilla gorilla), 28 male and 26 female adult chimpanzees (Pan), and 30 male and 24 female adult orangutans (Pongo) were studied. The data are described in detail by O’Higgins (1989) and O’Higgins and Dryden (1993). Eight landmarks are chosen in the midline plane of each skull as shown in Figure 1.18. The landmarks are anatomical landmarks and are located by an expert biologist.

The dataset is available in the R package shapes and can be accessed by typing:

data(apes)

The dataset is stored as a list with two components: apes$x is an array of coordinates in eight landmarks in two dimensions for each skull; and apes$group is a vector of group labels. To print the k × m × n array of landmarks in R, type apes$x and to see the group labels type apes$group (‘gorf’for female gorillas, ‘gorm’ for male gorillas, ‘panf’ for female chimpanzees, ‘panm’ for male chimpanzees, ‘pongof’ for female orangutans, and ‘pongom’ for male orangutans). In order to plot the landmark data we can use:

par(mfcol=c(2,3))
plotshapes(apes$x[,,apes$group=="gorf"],col=1)
title("Female Gorillas")
plotshapes(apes$x[,,apes$group=="gorm"],col=2)
title("Male Gorillas")
plotshapes(apes$x[,,apes$group=="panf"],col=3)
title("Female Chimpanzees")
plotshapes(apes$x[,,apes$group=="panm"],col=4)
title("Male Chimpanzees")
plotshapes(apes$x[,,apes$group=="pongof"],col=5)
title("Female Orang Utans")
plotshapes(apes$x[,,apes$group=="pongom"],col=6)
title("Male Orang Utans")

Again the plotshapes function plots 2D (x, y) coordinates of each object, and here we have drawn each group using a different colour. The plot is given in Figure 1.19.

It is of interest to assess whether there is a size difference between the sexes and whether there are any shape differences between the sexes in the face and braincase regions. A biologist would also be interested in geometrical descriptions of the shape difference, and how shape relates to size and other covariates.

Note that the coordinates in apes$x also available in the shapes package individually by group: gorf.dat, gorm.dat, panf.dat, panm.dat, pongof.dat and pongom.dat.

1.4.9 Bioinformatics: Protein matching

A protein is a sequence of amino acids, of which there are 20 types, and each amino acid has a one-letter code: A, C, D, E, F, G, H, I, K, L, M, N, P, Q, R, S, T, V, W, Y. A protein structure is the 3D configuration of atoms determined by the sequence and the 3D size and shape of the protein structure determines its function. The structure is summarized using key atoms, such as α-carbon atoms in each amino acid, and other secondary structure features such as α-helices and β-sheets (e.g. Mardia 2013b). Predicting the 3D structure of the protein from the amino acid sequence (protein folding) is one of the grand challenges in science.

Green and Mardia (2004, 2006) considered matching two proteins from the PDB data bank http://www.rcsb.org/pdb (Berman et al. 2000). The particular pair of proteins that were compared had PDB codes 1a27 (Human Type I 17Beta-Hydroxysteroid Dehydrogenase) (Mazza 1997) with 63 points; and 1cyd (Mouse lung carbonyl reductase) (Tanaka et al. 1996) with 40 points. The points are the active sites of the proteins and consist of the coordinates of the centres of gravity of the amino acids that make up the nicotinamide adenine dinucleotide phosphate (NADP) binding sites of two proteins. The 20 amino acid types of each active site are labelled in one of four groups as follows:

• Group 1: Hydrophobic. A, F, I, L, M, P, V.
• Group 2: Charged. D, E, K, R.
• Group 3: Polar. C, H, N, Q, S, T, W, Y.
• Group 4: Glycine. G.

Representations of the full proteins are given in Figure 1.20.

It is of interest to align the two molecules in order to find common similar geometrical parts. The size and shape of the proteins are of interest, and scale invariance is not required here. In this application the landmarks are unlabelled, and estimation of the correspondence between subsets of the active sites of each protein is required.

1.4.10 Particle science: Sand grains

A dataset of the curved outlines of sand grain profiles is available. There are 24 sea sand and 25 river sand grain profiles in two dimensions. The original data were provided by Dietrich Stoyan (Stoyan and Stoyan 1994; Stoyan 1997). On each particle outline there are 50 points, which were were extracted at approximately equal arc lengths by the method described in Kent et al. (2000, section 8.1). The sea particles are from the Baltic Sea and the river particles from the Caucasian River Selenchuk. It is of interest to describe the differences in shape variability of the sand grains between the two groups. Here we are interested in the shapes of the continuous closed curves. The points on each outline are unlabelled – there is no natural correspondence between points. The data are displayed in Figure 1.21.

data(sand)
plotshapes(sand$x[,,sand$group=="sea"],
   sand$x[,,sand$group=="river"],joinline=c(1:50,1))

1.4.11 Biology: Rat skull growth

We consider a well-known dataset of landmarks located on X-rays of rat skulls as they grow. The data are described in Bookstein (1991) and studied by several other authors including Goodall and Lange (1989); Monteiro (1999); Le and Kume (2000a); Kent et al. (2001); and Kenobi et al. (2010). The rats were carefully X-rayed at ages 7, 14, 21, 30, 40, 60, 90 and 150 days, and there are 18 rats with complete sets of eight landmarks in two dimensions at each age. The X-rays were taken of the skulls of the same rats recorded throughout their lifetimes, and it is of interest to describe the size and shape changes as the rats undergo growth.

The data for the 18 rats with complete data are displayed in Figure 1.22 and are available in the shapes package using:

data(rats)

with the landmark coordinates in rats$x, identification number in rats$no and the time in rats$time.

1.4.12 Biology: Sooty mangabeys

Twelve landmarks are taken from the midline of the skulls of a type of monkey, sooty mangabey (Cercocebus atys), in a study described by O’Higgins and Dryden (1992), see Figure 1.23. The specimens ranged from young juveniles to an adult female and an adult male. The objective is to describe the size and shape differences in the individuals in the series from the young juveniles to the older juveniles, and then to the adults. A further problem is to examine whether the individuals can be modelled by a regression line in shape space.

The data are available in the shapes package using:

data(sooty)

The data is an array of size 10 × 2 × 7, with the first two configurations being centred, rotated coordinates for the smallest juvenile and the male adult, and then the last five observations contain the original data for the three juveniles, then the female adult and the male adult.

1.4.13 Physiotherapy: Human movement data

Kume et al. (2007) consider an application in the study of human movement data, consisting of k = 4 landmarks (lower back, shoulder, wrist and index finger) moving in time. The landmark coordinates are obtained by recording the 3D locations of small reflective markers using a system of seven video cameras. Each individual in the study sat at a table and was asked to move his or her index finger towards a target point which was positioned straight, or to the right or or to the left at different angles. The data were collected by Dr James Richardson, Université Paris Sud, France. We are interested in the shapes of the landmark configurations, and we concentrate on the shapes of the configurations in the plane of the table and this subset was considered by Kume et al. (2007). We shall concentrate on a subset of the data where five curves are available, which are labelled a, b, c, d, e. The dataset consists of the projected view of the movements in the plane of the table at which the subject is sitting. For each of the individuals we have 10 equally spaced time points (after initially linearly transforming the times to [0, 1] for each curve). The data are available as a four-dimensional array in the shapes package using:

data(humanmove)

The array is of dimension 4 × 2 × 10 × 5 and represents the coordinates of the k = 4 landmarks in m = 2 dimensions with 10 time points for each of the 5 movements.

1.4.14 Genetics: Electrophoretic gels

A technique for the identification of proteins involves the comparison of electrophoretic gel images (Horgan et al. 1992). Two examples of such images are gel A and gel B shown in Figure 1.24. The images were obtained from particular strains of parasites which carry malaria. The objective is to use the gel image to identify the strain of parasite.

In each gel there are a number of black spots, where each spot can be one of two types – invariant or variant. The invariant spots are present for all parasites and the arrangement of variant spots enables identification of the parasite. A problem with the technique when used in the field is that the gels are prone to deformations and so the gel images first need to be ‘registered’ (transformed) so that direct comparisons can be made.

In this application, 10 invariant spots have been picked out by an expert, as shown in Figure 1.25. The spots are available in the shapes package using:

data(gels)

The invariant spots are used to match gel A to gel B, either by a similarity transformation or by a more complicated transformation. A question of interest is: can a matching procedure be made resistant to some outlier points, for example mislabelled points?

1.4.15 Medicine: Cortical surface shape

In a study of the shape of cortical surfaces of schizophrenia patients and normal volunteers, some structural MR images were taken of the brain (Brignell et al. 2010). The dataset consists of n = 68 3D MR images of the brain from 29 male healthy controls, 25 male schizophrenia patients, 9 female healthy controls, and 5 female schizophrenia patients. The MR images are proton density weighted images and were collected by Sean Flynn at the University of British Columbia, Canada. Each volunteer’s image consists of 256 × 256 × 256 voxels (3D pixels of size 1 mm³). The cortical surface was extracted from each dataset using image processing techniques [specifically the brain extraction tool of Smith (2002)]. Finally 62 501 points were located on the surface of the cortex above a transverse plane passing through the anterior and posterior commisures, with the brain midline plane being located in a sagittal plane. We see an example of the cortical surface points of one of the subjects in Figure 1.26. It is of interest to describe the cortical surface shape, and whether there are any differences in shape between the schizophrenia group and the control group.

A subset of the data are available in the shapes package:

data(cortical)

which contains the 250 points on the outline intersecting the axial slice containing the commisures (the bottom-most axial cross-section in Figure 1.26).

1.4.16 Geology: Microfossils

The microfossil Globorotalia truncatulinoides is a microscopic planktonic found in the ooze on the ocean bed. Lohmann (1983) published 21 mean outlines of the microfossil which were based on random samples of organisms taken at different latitudes in the South Indian Ocean. Figure 1.27 shows the three mathematical landmarks selected on each of the outlines, where one landmark has been placed at the origin. The coordinates of the landmarks are extracted from Figure 7 of Bookstein (1986). It is of interest to examine whether the size of the organisms is related to the shape, and whether size or shape are related to the covariate of latitude. A more basic problem would be to obtain an estimate of the average shape of the fossils and to describe the structure of the shape variability.

The data are available in the shapes package using:

data(shells)

1.4.17 Geography: Central Place Theory

Central Place Theory was postulated by Christaller (1933) and is the situation where towns are distributed on a regular hexagonal lattice over a homogeneous area (with towns at hexagon centres, see Figure 1.28). Mardia et al. (1977) consider this hypothesis for a map of 44 places in 6 counties in Iowa, namely Union, Ringgold, Clarke, Decatur, Lucas and Wayne.

In order to examine whether Central Place Theory holds, one could examine the shapes of the triangles formed by a town and its neighbours to see if they are more equilateral than expected under a randomness hypothesis. A convenient triangulation of the towns is a Delaunay triangulation (Mardia et al. 1977; Green and Sibson 1978; Okabe et al. 2000). In Figure 1.28 we see that Voronoi polygons for ideal central places would be hexagons and Delaunay triangles would be equilateral triangles. In Figure 1.29 we see a Delaunay triangulation and the Voronoi polygons for the Iowa data. An important question to ask is: are the Delaunay triangles more equilateral than expected by chance?

The points here are unlabelled (there is no correspondence in the vertices of the triangles). Also the triangles are correlated due to neighbouring triangles sharing points. Kendall (1983, 1989) also studied shape in Delaunay triangulations, in order to investigate the Central Place Theory hypothesis, following Mardia et al. (1977).

1.4.18 Archaeology: Alignments of standing stones

Consider a map of the 52 megalithic sites that form the ‘Old Stones of Land’s End’ in Cornwall, UK, given in Figure 1.30. It was proposed by Alfred Watkins, in the early 1920s, that these and other megalithic sites were placed in deliberate straight lines, called ley lines. One approach is to consider the shapes of all possible triangles and to see if there are more ‘flat’ triangles (triangles with the largest angle close to 180°) than expected under a randomness hypothesis. The points are unlabelled, and in this dataset there are triangles in two dimensions.

This dataset is particularly important in the history of shape analysis because it motivated D.G. Kendall’s pioneering work. Analysis of these data is considered by Broadbent (1980), Kendall and Kendall (1980), Small (1988) and Stoyan et al. (1995) among others.