JSS Journal of Statistical Software

February 2014, Volume 56, Issue 13. [Link]

Momocs: Outline Analysis Using R

Vincent Bonhomme Sandrine Picq

French Institute of Pondicherry UMR CBAE, Montpellier

Cédric Gaucherel Julien Claude

French Institute of Pondicherry University of Montpellier II

Abstract
We introduce here Momocs, a package intended to ease and popularize modern mor-
phometrics with R, and particularly outline analysis, which aims to extract quantitative
variables from shapes. It mostly hinges on the functions published in the book entitled
Modern Morphometrics Using R by Claude (2008). From outline extraction from raw
data to multivariate analysis, Momocs provides an integrated and convenient toolkit to
students and researchers who are, or may become, interested in describing the shape and
its variation. The methods implemented so far in Momocs are introduced through a
simplistic case study that aims to test if two sets of bottles have different shapes.

Keywords: modern morphometrics, comparison of shapes, Fourier transform, outline analysis,

thin plate splines, R.

1. Introduction

1.1. The aim and purpose of morphometrics

The link, if there were one, between the form and the function of objects, living or inert, has
been one of the most enduring questions in the realm of science. In many situations, analyzing
the shape variation among objects can bring deep insights into their functioning and to the
underlying mechanism leading to their variation in shape. For instance, an evolutionary biol-
ogist may be interested in testing the proximal link between the shape and the photosynthetic
capacities, hydraulic properties, or the way leaves develop from a bud. Similarly, distal causes
such as local adaptations in comparing different populations, or the diversification of foliar
organs along the evolutionary lineage of a species, can be investigated. In a very different
2 Momocs: Outline Analysis Using R

field, geographers would like to test if differences in the shape of cities can be attributed to
external factors such as highways crossing them or oceans bordering them. In addition, mor-
phometrics, the so-called quantitative study of form, and derived techniques such as shape
clustering and recognition, can concern a wide spectrum of disciplines.
Morphometrics aims at analyzing the variation and covariation of the size and shape of objects,
defining altogether their form. Shape and form might be confusing words, used as synonyms
in many languages. Hereafter, we will use the definition of shape proposed by Kendall (1989)
and Small (1996) that it is “the total of all information invariant under translations, rotations,
and isotropic rescaling”. What remains if we draw a heart, and then rotate the page, or change
its size uniformly, or move the page about on the desk, is the shape of the heart. Shape by
essence is better described in three dimensions, except for flat objects (see Kendall 1989,
for a review). In this paper, however, only two dimensional shape will be considered, e.g.,
three-dimensional objects will be viewed from one side and considered as represented by their
projections on a plane.
How to quantify the shape of the heart that was described above, and how to compare a
given number of hearts to show variations among themselves? This paper aims to provide
a didactic introduction to the shape concept, and particularly a branch of morphometrics
relying on outline analysis. Finally, a case study of outline analysis is presented; it uses
Momocs, the R package presented here, whose name stands for modern morphometrics, and
is intended to simplify such morphometrical analyses.

1.2. How to compare shapes?

The everyday approach to describing shapes is to use words, such as “round”, “narrow”, “heart-
shaped”, and “symmetric”. Such usage points to a very limited vocabulary access: shapes can
be too complex, and differences between them too subtle for words. More importantly, they
are ad hoc descriptors. In describing the shape of a different object, one will probably not
use the same vocabulary thus making comparison hard to appreciate.
A quantitative framework was introduced by traditional morphometrics which measures dis-
tances, areas, etc. and compares them in an uni- or multivariate framework (see Rohlf and
Marcus 1993). For instance, in describing and comparing several human faces, the lengths of
the ears and the noses, the interpupillary distances, and other lengths or ratios of lengths,
taken homogeneously between the individuals, are used to test differences between genders
or the covariances between parts of the body. In his seminal book On Growth and Form,
Thompson (1917) compared shapes in what could be called today a morphometrics approach,
and offered new ways of understanding their variations: i) some changes in the developmental
processes of living organisms, while minor, can lead to dramatic morphological changes, and
ii) physical constraints such as growing mechanics are of first importance in the final form of
organisms (Figure 1). Morphometrics contributed towards raising development from a simple
bridge between genes and organisms to a central catalyst of evolutionary change between
species. Many attempts have also been made outside biology, notably using scalar indices
to describe and compare shapes: compactness, elongation, fractal dimension, etc. are such
examples used to compare boundaries of political or physical geographical objects such as
cities, states, or watersheds, etc. (Moellering and Rayner 1981; Wentz 2000; Miller and Wentz
2003). In particular Momocs and the R routines, on which it hinges, have already been used
to study the influence of height and body mass on human bodies outlines (Courtiol, Ferdy,
 Journal of Statistical Software 3

•
•
•
•
•
•
•
•

Figure 1: These popular D’Arcy Thompson’s fishes illustrate his central thesis: “An organism
is so complex a thing, and growth so complex a phenomenon, that for growth to be so uniform
and constant in all the parts as to keep the whole shape unchanged would indeed be an unlikely
and an unusual circumstance. Rates vary, proportions change, and the whole configuration
alters accordingly” (cited from Thompson 1917). From a morphological point of view, these
two fishes seem very different, but a little change in the growth rate of the caudal (right) part,
here illustrated by the grid and four particular points taken on it, leads to dramatic change.
This example serves to demonstrate that, by providing a rigorous method of describing shape,
morphometrics can help a better understanding of how entities in nature function. Redrawn
after Thompson (1917).

Godelle, Raymond, and Claude 2010), to provide complementary information in the shape
description of watersheds where only scalar indices are classically used (Bonhomme, Frelat,
and Gaucherel 2013a) and the first evidence of intraspecific variability in the shape of pollen
grains of anemophilous species (Bonhomme, Prasad, and Gaucherel 2013b).
At the other end of the XXth century, computerized data acquisition and treatment arose
synchronously with an array of new methodological developments and these techniques of
modern morphometrics revolutionized the historical scope of morphometrics (Rohlf and Slice
1990; Rohlf and Marcus 1993; Adams, Rohlf, and Slice 2004).

1.3. Modern morphometrics: Landmark configuration and outline analysis

Modern morphometrics considers shape as a whole, taking into account all the geometrical
relationships of the input data. The two main approaches in use are: the study of landmark
configurations, and outline and surface analyses. Both of them preserve the geometrical
information, i.e., relative positions between all points are kept (Moellering and Rayner 1981;
Kuhl and Giardina 1982). This allows shape reconstruction from their numerical signature,
a fact which is of great interest since we can then define the most frequent shapes, the rare
or the impossible ones, infer intermediate shapes, etc. in other words, these approaches can
reveal some functional links between the shape and its variation with the processes leading
to it.
Configuration of landmarks can be summarized as follows: the relative positions of a set of
points, called landmarks, are considered globally e.g., by using a matrix of their pairwise
4 Momocs: Outline Analysis Using R

• •
•
2

•
•
3 •
• •
1
3

• •
1

•
2
• • •

Figure 2: Two biological shapes that can be explored with modern morphometrics. The mouse
jaw on the left, redrawn after Claude (2008), exhibits many potential landmarks. Some of
them are (1) frontiers between bones, (2) sharp angles on a single bone, and (3) the region
with high curvature. The same approach cannot be easily applied with the Ginko biloba leaf
on the right since less landmarks can be identified or are not consistently present on Ginko
leaves, for instance the separation between the two lobes. Outline analysis can become a
solution in cases such as this.

euclidean distances (Richtsmeier, Cheverud, and Lele 1992; Richtsmeier, Burke Deleon, and
Lele 2002). The landmarks need to be structurally similar (ideally homologous) between
individuals. For instance, if the variation of the vertebrate skull is considered, the bones
that share a common origin are defined as homologous and can be part of a configuration
of landmarks along with any other geometrical feature that can be unambiguously identified
(Figure 2, see also Macleod 1999, for a general discussion on homology).
Other objects, for instance a leaf, do not display any or display too few landmarks. Methods
for analyzing configurations of landmarks cannot be easily applied here, and outline analysis
will be used instead. This second approach considers outline as a whole. Outline is defined
here as the closed polygon formed by the (x; y) coordinates of pixels defining it. One popular
outline analysis approach, on which Momocs is focused so far, uses Fourier series to describe
shapes, and is detailed below. Finally, three dimensional surfaces or outlines are an other
area of morphometrics, but will not be considered here.

1.4. Momocs: Analysis of outline variation using R

This paper introduces Momocs ([Link] a package
to analyze outlines of shapes, using R (R Core Team 2013). The package originated from some
core functions published by one of us (Claude 2008) and reviewed by Bowman (2009). These
functions were turned into an integrated framework and a standalone R package. Below are
provided step-by-step guidelines for performing modern morphometrics, from which Momocs
derives its name, using R. The package’s vignette A Graphical Introduction to Momocs and
Outline Analysis Using R (Bonhomme 2012) also provides an extensive description of the
functions of the package.
While both outline analysis and R (to a much greater extent) have been used in increasing
measure, so far no dedicated tool has been aviable at CRAN. Momocs aims to fill this gap.
Other tools exist, but they focus on configuration of landmarks: shapes by Dryden (2012),
MorphoJ by Klingenberg (2011), and recently geomorph by Adams and Otarola-Castillo
 Journal of Statistical Software 5

(2012). Outline analyses can be performed by some standalone programs but often not under
an open license and only on certain operating systems. The Stony Brook University webpage
lists most of them ([Link] the SHAPE suite by Iwata, Ne-
sumi, Ninomiya, Takano, and Ukai (2002) being broadly used. Momocs is placed under the
GPL license and through this license, it may become a collaborative hub for other researchers
to explore new methods or implement existing approaches.

2. Mathematical background of Fourier-based outline analyses

Approaches for analyzing outlines rather than landmark configurations estimate parameters
of functions rather than the relative positions of landmarks after superimposition. One of the
strategies is to adjust Fourier series to some shape descriptors, and this is the main approach
followed, so far, by Momocs.

2.1. Fourier transformations and closed outlines

Outline analysis does necessarily require that the outline on which points are sampled to be
structurally defined: what is extracted is the geometrical information contained in the outline
itself and taken as a whole. However, when comparing shapes, outlines should correspond
to structurally similar features. In particular, Fourier-based approaches are powerful enough
to extract this geometric information. They have their basis in the idea of Fourier series:
to decompose a periodic function into a sum of more simple trigonometric functions such as
sine and cosine. These simple functions have frequencies that are integer multiples, i.e., are
harmonics, of one another. The lower harmonics provide approximations for the coarse-scale
trends in the original periodic function, while the high-frequency harmonics fit its fine-scale
variations.
Fourier series can be used in morphometrics, amongst many other derived applications since
closed outlines can be considered as periodic functions. If we start somewhere on the outline
and follow it, we will pass again and again by the same starting point and thus periodic
functions can describe this outline. Such functions are: the distance of any point on the
outline to the centroid of the shape, the variation of the tangent angle for any point, or the
(x; y) coordinates on the plane. A, or several, periodic functions are then obtained and can
be decomposed (and thus described) by Fourier series. These three different methods are
available in Momocs, and hereafter called “radius variation”, “tangent angle” and “elliptical
analysis” (Figure 3), and their comparison has been discussed extensively by Rohlf and Archie
(1984).
The principle of Fourier series described above applies to continuous functions. Since a shape
is based on a finite number of discrete points, typically coordinates on a plane (or a space), a
discrete equivalent of Fourier series is used in morphometrics. A given number of points called
pseudo-landmarks, have to be sampled along the outline before performing analysis of outline
variation. All Fourier decompositions result in an harmonic sum of trigonometric functions
weighted with harmonic coefficients. They are (usually) normalized to remove homothetic,
translational or rotational differences between shapes. Two or four coefficients, depending on
the approach used, are obtained for each harmonic calculated and can then be considered as
quantitative variables. Nyquist frequency precludes more harmonics than half the number
of points fitted, which is thus their upper limit. The geometrical information contained in
6 Momocs: Outline Analysis Using R

Tangent angle

Variation of the tangent angle (in radians)

2π
400 4 10 15 19
5 3
6
200 2 3π/2 3
7 1
8
10 9 20
0 11 π
9
12 14
13 19
−200 6 7 8
π/2 2 13 14
15 12
16 18 11 18
5 16 20
−400 4
17 0 1 17

−1000 −500 0 500 1000 0 5 10 15 20

Points along the outline

Radius variation Elliptical analysis

Deviation from the first point (in pixels)

500 250 4
1 3 dx
5 dy
Radius length (in pixels)

6
400 5 16 17 50 2
4 18 19 20 1 7
2 3 15 8
10 19 20
300 6 10 11 −150 2 9 20
12 13 11 18
12
200 7 −350 3 13 14 19
14 16 17
4 8 9 14 15
15 16 18
100 −550 6 7
8 9 13 17
5 10
0 −750 11 12

0 5 10 15 20 0 5 10 15 20
Points along the outline Points along the outline

Figure 3: Twenty equidistant points have been sampled, starting from the beak, counter-
clockwise and along the curvilinear abscissa of the dove shape (top left) outline, inspired by
Picasso’s drawing. This outline can be described using Fourier-based methods. Tangent angle
(top right) illustrates the variation of the tangent angle along the outline. Radius variation
(bottom left) illustrates the length of the radius, here considered as the distance between
the center of the shape (the cross within the dove outline) and the points along the outline.
Elliptical analysis (bottom right) shows the two curves corresponding to xn − x1 (in blue) and
yn − y1 (in red).

the outlines are thus quantified and can be analyzed with classical multivariate tools. The
following sections detail the core functions of Momocs, and an extensive description can be
found in Rohlf and Archie (1984) and Claude (2008).

2.2. Fourier radius variation

Zahn and Roskies (1972) stated that, given a closed outline, the radius r, taken as the distance
from the outline centroid and a given point of the outline, can be expressed as a periodic
function of the angle θ. Harmonics from 0 to k approximate the function:

k
1 X
r(θ) = a0 + an cos(wn θ) + bn sin(wn θ)
2
n=1

with:

p p
2X 2X
an = ri cos(nθi ), bn = ri sin(nθi ),
p p
i=1 i=1
 Journal of Statistical Software 7

and:
r p
2X
a0 = ri .
p
i=1

w refers to the pulse and p is the number of sampled points along the outline (equivalent here
to the number of sampled radii in this case). The an and bn harmonic coefficients, extracted
for every individual shape, can be used for multivariate analyses to compare a set of outlines.

2.3. Fourier tangent angle

Radius variation may fail to fit some complex outlines, in particular when a given radius
intercepts the outline twice, a situation that can arise when the outline presents convexities
and concavities. Zahn and Roskies (1972) proposed also another approach. The Fourier
tangent angle fits the cumulative change in the angle of a tangent vector (φ(t)), as a function
of the cumulative curvilinear distance t along the outline.
Given a closed outline, previously scaled to 2π, φ(t) can be expressed as:
φ(t) = θ(t) − θ(0) − t,
where t is the distance along the outline, θ(t) the angle of the tangent vector at t and θ(0)
the angle of the tangent vector taken for the first point. It can be removed for normalizing
the coefficients obtained. Two coefficients per harmonic can be estimated as follows:
p p
2X 2X
an = φ(t) cos(nθi ), bn = φ(t) sin(nθi ),
p p
i=1 i=1

and:
r p
2X
a0 = φ(t)
p
i=1

2.4. Elliptic Fourier analysis

The last approach presented here is due to Giardina and Kuhl (1977) and Kuhl and Giardina
(1982) who developed a method for fitting separately the x and y coordinates of an outline
projected on a plane. This method has become very popular since it has great advantages
over the other Fourier-based approaches: equally spaced points are not required, virtually
any outline can be fitted (see Rohlf and Archie 1984; Crampton 1995; Renaud and Michaux
2003) and the coefficients can be made independent of outline position and normalized for
size.
Let T be the perimeter of a given closed outline, here considered as the period of the signal.
One sets ω = 2π/T to be the pulse. Then, the curvilinear abscissa t varies from 0 to T . One
can express x(t) and y(t) as follows:
+∞
a0 X
x(t) = + an cos(nωt) + bn (sin nωt)
2
n=1
8 Momocs: Outline Analysis Using R

with:
ZT ZT
2 2
an = + x(t) cos(nωt)dt, bn = + x(t) sin(nωt)dt.
T T
0 0

Similarly,
+∞
c0 X
y(t) = + cn (cos nωt) + dn (sin nωt)
2
n=1
with:
ZT ZT
2 2
cn = + y(t) cos(nωt)dt, dn = + y(t) sin(nωt)dt.
T T
0 0

Since the outline contains a finite number of points given by k, one can calculate discrete
estimators for every harmonic coefficient of the nth rank:
k  
T X ∆xp 2πntp 2πntp−1
an = 2 2 cos − cos ,
2π n ∆tp T T
p=1
k  
T X ∆xp 2πntp 2πntp−1
bn = 2 2 sin − sin .
2π n ∆tp T T
p=1

∆x1 = x1 − xk and cn and dn are calculated similarly. a0 and c0 correspond to the estimate
of the coordinates of the centroid of the original outline and are estimated by:
p p
2X 2X
a0 = xi , c0 = yi .
T T
i=1 i=1

Intuitively, for all positive integers n, the sum of a cosine curve and a sine curve represent
the nth harmonic content of the x and y projections of the k-edged polygon, and for any
n, these two curves define an ellipse in the plane. Ferson, Rohlf, and Koehn (1985) noticed
that in the “time”, say one period, it takes the nth harmonic to traverse its ellipse n times,
the (n + 1)th harmonic has traversed its own ellipse n + 1 times (Figure 4). Ferson et al.
(1985) noticed that the reconstruction of the original polygon is done by vector-adding these
ellipses for all harmonics, which echoes (the ancient astronomer) Ptolemy’s epicycles, and the
reconstruction obtained from N harmonics is the best possible fit in a least-squares sense.
In elliptic Fourier analysis, four coefficients per harmonic are obtained, two for x and two for
y. We can use the first harmonic, the one that defines the best-fitting ellipse, to normalize the
harmonic coefficients and make them invariant to size and rotation. The harmonic coefficients
can also be normalized for the location of the first outline coordinate. When this is done, the
shapes are individually aligned according to their first fitted ellipse. If at least one homologous
point can be defined, one can rather use them to align the outlines. Normalized elliptic Fourier
coefficients, further symbolized by An , Bn , Cn and Dn , are obtained:
     
An B n 1 cos ψ sin ψ an bn cos nθ − sin θ
=
Cn D n λ − sin ψ cos ψ cn dn sin nθ cos nθ
 Journal of Statistical Software 9

Figure 4: For all positive integers, the sum of a cosine curve and a sine curve defines an
ellipse in the plane. Elliptic Fourier analysis is based on an harmonic sum of such ellipses
as in Ptolemy’s astronomical system with higher harmonic order ellipses “rolling” within all
lower order ellipses. Three harmonics are here shown at four locations on the original outline.

The scale λ is estimated as the magnitude of the semi-major axis of the ellipse as defined by
the first harmonic. The second right term corresponds to the orientation of the first ellipse,
with ψ being the rotation angle, the third to the original harmonic coefficient, and the last to
the rotation of the starting point to the end of the ellipse, with a rotation angle of θ. Ferson
et al. (1985) also supplied the following formulas with which to calculate these parameters:
2 × (a1 b1 + c1 d1 )
ψ = 0.5 arctan
a21 + c21 − b21 − d21
with first p
λ = a∗2 + c∗2
and
θ = arctan(c∗ /a∗ )
with
a∗2 = a1 cos ψ + b1 sin ψ
c∗2 = c1 cos ψ + d1 sin ψ

3. Momocs: Outline analysis using R

Preliminaries
Momocs is S4-oriented (Chambers 1998) which has many advantages in terms of usage and
10 Momocs: Outline Analysis Using R

whisky
beer

Figure 5: The bottles dataset included in Momocs and analyzed here.

programming: it prevents typing errors, provides validity checking, allows inheritance and
encapsulation, etc. (see Genolini 2008). In practical terms, lists of coordinates and matrices
of harmonic coefficients, are handled through ‘Coo’ and ‘Coe’ class objects respectively, to
which methods can be applied. For those not familiar with S4 objects, data stored in objects
can be retrieved and used as classical S3 objects in R: matrices, factors, etc. The Momocs
documentation provides an extensive description of these classes and the methods that can
be applied to them. The case study presented below will focus on the basic (and probably
typical) use of the package.

All the following examples are based on the bottles dataset from the package (see ?bot and
Figure 5). We want to test if whisky and beer bottles have different shapes. How to calibrate
outline analysis parameters and then obtain a matrix of normalized harmonic coefficients will
also be discussed. While elliptic Fourier analysis will be presented because it is one of the most
popular outline analysis tools today, the methodology employed will be equally valid for other
approaches as well. On the extracted harmonic coefficients, some multivariate analyses will
be presented: principal component analysis and morphological space, to illustrate the global
bottle diversity, and multivariate ANOVA, to test for shape difference between the two sets
of bottles. Then, linear discriminant analysis and hierarchical clustering will be introduced as
perspectives for Momocs and because they are common and helpful statistical tools for those
interested in multiple comparison. Finally, thin plate splines analysis will be introduced: this
is not only a tribute to D’Arcy Thompson’s work but it may also bring great insights into the
developmental differences underlying differences in the shapes compared.
 Journal of Statistical Software 11

3.1. Outline extraction

The input data for the morphometric analysis of outlines is, for two-dimension shapes, a set
of (x; y) pixel coordinates sampled on each outline. Momocs provides facilities to convert a
set of .jpg black-and-white images, e.g., as classically obtained with image analysis stations
when working on biological objects (see ?[Link]). Other image formats will be included
in subsequent releases. The algorithm presented in Rohlf (1990), and implemented in R by
Claude (2008) is used to extract the list of coordinates from a black mask that corresponds
to the shape on a white background. Lists of coordinates, stored in separate .txt files, will
be a convenient alternative when data are derived from GIS (see ?[Link]).
Outlines are finally included in a ‘Coo’ class object. Outlines can be visualized in a one page
graph (Figure 5). They can be centered, aligned, scaled and homologous landmarks can be
defined to perform a Procrustes alignment (see Friess and Baylac 2003) before an elliptical
Fourier analysis. When the outlines become rough due to artifacts during the digitization
process (for instance when automatic outlining produces noise around the outline), outlines
can be smoothed either when outlines are extracted from images, or before the calculation
of harmonic coefficients (see ?[Link] and ?eFourier for instance). In order to specify
explanatory variables going along with the coordinate or coefficient set, grouping factors or
covariates can be specified through a [Link], and then used to create subsets (see ?Coo).
Finally, we emphasize that Momocs supports import/export from/to other formats used in
morphometrics such as .chc, .nef, etc.

R> library("Momocs")
R> data("bot", package = "Momocs")
R> panel(bot, cols = rep([Link](2), each = 20))
R> bot

A Coo object (see ?Coo)

******************************

General
--------------------
- 40 outlines
- 1864 +/- 371 coordinates per outline
- No landmark defined
- All outlines are closed
- 1 grouping factor(s) defined

Coordinates: @coo
--------------------
[,1] [,2]
[1,] 37 561
[2,] 37 560
[3,] 37 559
[4,] 37 558

[...]
12 Momocs: Outline Analysis Using R

Factors: @fac
--------------------
type : beer whisky

3.2. Calibration of outline analysis

Fourier-based approaches can fit any outline provided that the number of harmonics is large
and the outline smooth enough between sampled points, while the signal/noise ratio can be
very low for high order harmonics. The latters describe details that may be due to many
things, i.e., digitalization artifacts or user bias, but not to real differences between shapes.
On the other hand morphometrics is also used when differences between shapes are subtle.
Conflicting situations such as those suggested by this familiar saying are actually a recurrent
issue in morphometrics: what is the right number of harmonics? Unfortunately, no objec-
tive criterion exists so far, and the criterion used usually depends on the scope of the study.
This might not be fully satisfactory to morphometrics newcomers but some approaches are
presented below that can help to choose the most appropriate number of harmonics. Further-
more, a recent approach by Claude (2013) allows to study the measurement-error depending
on the harmonic rank.

Through shape reconstruction

First, a ‘Coo’-object can be passed to [Link]() to observe the reconstructed shape for a
range of harmonics (Figure 6).

R> hqual(bot, method = "eFourier", id = 16, [Link] = 1:49,

+ palette = [Link], [Link] = "panel")

Through deviations
A qualitative approach would be of limited value, quite unlike the method we have presented
to quantify deviations. The idea is to define, for a given number of sampled points, the best
possible fit (i.e., obtained with half this number of points), and to then compare the euclidean
distances obtained with a lower number of harmonics for every point of this outline and the
best possible outline with these sampled points (Figure 7). One can for instance choose the
minimal number of harmonics that leads to an average deviation of 1 pixel.

R> hquant(bot, [Link] = c(12, 16, 20, 24, 32))

Through harmonic power

Finally, we can also estimate the number of harmonics after examining the spectrum of
harmonic Fourier power. The power is proportional to the harmonic amplitude and can
be considered as a measure of shape information. As the rank of a harmonic increases, the
power decreases and adds less and less information. We can evaluate the number of harmonics
 Journal of Statistical Software 13

1 2 3 4 5 6 7

8 9 10 11 12 13 14

15 16 17 18 19 20 21

22 23 24 25 26 27 28

29 30 31 32 33 34 35

36 37 38 39 40 41 42

43 44 45 46 47 48 49

Figure 6: The “pecheresse” beer bottle in the bot dataset reconstructed from different numbers
of harmonics. Twelve harmonics give a very satisfactory reconstruction and for 20 the result
is almost perfect.

Deviations along the outline

4
Number of harmonics
12
16
20
24
32

3
Deviation in pixels

0
10 20 30 40 50 60
Points sampled along the outline

Figure 7: The bottles dataset included in Momocs and analyzed here. The y-axis represents
the deviation in pixels between the best possible fit with a given number of harmonics (here
32) for every sampled points along the outline (on the x-axis). Standard deviations for the
40 shapes are displayed.
14 Momocs: Outline Analysis Using R

Fourier power spectrum

1.0
Quantiles
0
0.5
1

0.8
Cumulative harmonic power
0.6
0.4
0.2

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
Number of harmonics included
40 outlines

Figure 8: Cumulated harmonic Fourier power calculated from the bottles dataset. The 12
first harmonics gather nearly 100% of the harmonic power. Maxima and medians are also
plotted.

that we must select, so their cumulative power gathers 99% of the total cumulative power
(Crampton 1995, Figure 8). The power of a given harmonic is calculated as:

A2n + Bn2 + Cn2 + Dn2

P owern =
2

Figure 8 is obtained by typing:

R> hpow(bot)

Computing elliptic Fourier analysis

Once the right number of harmonics has been determined, elliptic Fourier analysis is performed
on the ‘Coo’-object using the eFourier method and a

[Number of Outlines] × [4 coefficients × Number of Harmonics]

matrix is obtained along with grouping factors, individual names, etc. and returned as a ‘Coe’
class object. By default, the obtained coefficients are normalized so that the first fitting ellipse
is re-aligned along the x-axis. Other options can be considered, for instance, one can also
normalize by performing a Procrustes alignment when landmarks can be defined (see Friess
and Baylac 2003, ?ProcGPAlign and the package’s vignette which includes an illustration of
this approach).

R> botF <- eFourier(bot, nb.h = 20)

R> botF
 Journal of Statistical Software 15

A matrix of harmonic coefficients obtained with

elliptical Fourier analysis (see ?Coe)
******************************

General
--------------------
- 40 outlines
- 20 harmonics
- 1 grouping factor defined

[...]

Before multivariate analysis can be performed, one may be interested in having a global view
of the elliptic Fourier analysis: which coefficients vary and what is the geometrical variation
they depict (Figure 9). For instance, bn and cn represent the asymmetry of the shapes that
can vary from coarse to fine scales, i.e., from lower to higher rank harmonics (Iwata, Niikura,
Matsuura, Takano, and Ukai 1998; Yoshioka, Iwata, Ohsawa, and Ninomiya 2004). Figure 9
is obtained by typing:

R> hcontrib(botF, [Link] = 1:8)

R> boxplot(botF)

3.3. Analyzing Fourier coefficients

Principal component analysis (PCA) and other multivariate approaches can be directly per-
formed on this ‘Coe’ class object (or directly on the matrix stored in the @coe slot) since all of
the harmonic coefficients can be considered as quantitative variables. Both pros and cons of
multivariate techniques apply to an analysis of the harmonic coefficients, just as they apply to
other types of data. The purpose here is not describe them extensively but rather to present
what is currently implemented in Momocs.

Principal component analysis

Momocs takes profit of the ade4 package by Dray and Dufour (2007). The pca method
can be used on a ‘Coe’ class object and performs a PCA with centering but no rescaling
by default. In other words, the small-amplitude coefficients will contribute less than the first
coefficients. It returns a ‘dudi’ object to which all suitable ade4 functions can also be applied.
An almost exhaustive wrapper that gathers graphical functions from ade4 such as the display
of eigenvalues, confidence ellipses and “stars”, individual labeling, neighboring graphs, etc. and
that also adds dedicated features such as the display of the morphological space is provided
by [Link]. See ?[Link] for an exhaustive description of this highly tunable function.
Plotting the PCA (Figure 10) is straightforward: below, we first compute elliptical Fourier
analysis with 20 harmonics, get out a ‘dudi’ object, and finally plot it.

R> botF <- eFourier(bot, nb.h = 20)

R> botD <- pca(botF)
16 Momocs: Outline Analysis Using R

Harmonic contribution

0
0.5
Amplification factor
1
2
5
10

1 2 3 4 5 6 7 8
Harmonic number

Variation of harmonic coefficients

0.10

Harmonic coefficients
A
B
C
D
0.05
Coefficient value

0.00
−0.05
−0.10
−0.15

2 3 4 5 6 7 8
Harmonic number

Figure 9: Analysis of harmonic coefficients after an elliptic Fourier computation on the bot
dataset. The top figure shows the effect of every harmonic on the shape reconstruction. Every
harmonic is represented on the x-axis; then the corresponding coefficients are multiplied by
values that illustrate their removal (0), the normal shapes they lead to (1) or their exaggerated
effect in shape reconstruction (2 and above). When the first harmonic is removed, on which all
other ellipses “roll on”, the reconstructed shape is obviously very bad. Exaggerated coefficients
help to understand their contribution to the final shape, e.g., the second harmonic multiplied
5 or 10 times indicates that it contributes to describe the bottleneck while the third and
fourth harmonic contribute to the constriction on the middle of the bottle. The bottom figure
illustrates the variation of every coefficient along the whole dataset. Here, both the Bn and
Cn coefficients are very small depicting the (bilateral) symmetry of the studied shapes.
 Journal of Statistical Software 17

Principal Component 2

Principal Component 1
Principal Component 2

beer whisky

Eigenvalues

Principal Component 1

Figure 10: Some examples of the factorial maps depicting morphological variation that can
be obtained with Momocs. The first two principal component axes are shown, (PC1 and PC2
are the x- and y-axis, respectively). The top figure displays the factorial map with no display
of the classes, a neighboring graph, shapes reconstructed from the factorial map using the
first two PC axes, and “rugs” along the axes; the bottom figure shows eigenvalues, confidence
ellipses for the two groups and another option for displaying the morphological space, as a
grid of bottles.

R> [Link](botD, [Link] = "circle", neighbors = TRUE)

R> [Link](botD, 1, pch = 20, eigen = TRUE,
+ [Link] = 0.5, [Link] = 1.5, rug = FALSE)
18 Momocs: Outline Analysis Using R

Multivariate analysis of variance

We can test for a difference between subsets of shapes using multivariate analysis of variance
(MANOVA), with every harmonic coefficient being considered as an homologous quantitative
variable measured for every shape of the dataset. This can be achieved with:

R> [Link](botF, "type")

The number or retained harmonics was not specified.

Analysis done with 9 harmonics
Df Hotelling-Lawley approx F num Df den Df Pr(>F)
fac 1 3826.1 318.84 36 3 0.000247 ***
Residuals 38
---
Signif. codes: 0 "***" 0.001 "**" 0.01 "*" 0.05 "." 0.1 " " 1

The outcome of this analysis shows that shapes of the whisky and beer bottles significantly
differ.

Hierarchical clustering
Momocs also includes a method to perform hierarchical clustering that hinges on dist and
hclust for calculation, and [Link] from the ape package for graphical output (see Par-
adis, Claude, and Strimmer 2004 and Paradis 2012). This can be achieved with the code
below (Figure 11):

R> clust(botF)

Call:
hclust(d = dist(botE@coe, method = method))

Cluster method : complete

Distance : euclidean
Number of objects: 40

Other methods, such as linear discriminant analysis, etc. exemplified in Claude (2013) will be
integrated in further versions of Momocs.

Thin plate splines

Deformation grids as those that contributed to popularizing D’Arcy Thompson’s ideas can
be obtained using thin plate splines mathematical formalization. The notion of thin plate
splines has been borrowed from mechanics and involves the bending of a thin sheet of metal
(see Bookstein 1991). The deformations required to pass from the mean shape to the extreme
points of the morphological space can be calculated and displayed on the PCA. One can also
perform thin plate splines analysis based on the harmonic coefficients and the reconstructed
shapes (Figure 12).
 Journal of Statistical Software 19

gie
wh

whisky

oran
isk

whis
yg

llan
oban

lenm
len

tamdhu

whisky jb
ky a

ga

ffe
dro
be

k
ky g

ou
a
b

rie
whisky
ee

ym
er

mru

ch
na

sk
f

rc
ra

an
whis

ce
ch

isk
be

nz
y

m
ke

i
er

ro

rn
is

de
ur

wh
k

ka

na
ing t er

lin
ild
alk

ne

be
fis
yw

er
ew

r
be he k i

be
er r is n
hn
bee can wh ky jo e
r ho ey is rous
ega
rdee w h
a m ousg
n ky f
beer we whis
stmalle hma
beer bra
beer pecheresse whisky bushmills
iler beer sie
beer jup rraneva
da
bee
ro r gri
r tau wh m b
bee s isk erge
n
tine be y ja
allan el er ckd
a
yb uv gu nie

w
s k r d ine ls

hi
i
e

wh e ss

be

sk
pp

be
be

y
er
tra

whis

yo
whisky
tr u

er
la

tan
ero

ic
re
er

en

hi
whisky chivas

chy

gle
mark
be

ky h
otr
almo
v
us

foo
ma
ldp
de

redbrea

ighla

t
makers
ky d

y
yo
er
be

isk

ndp
whis
wh

st

ark
whisky

Figure 11: An example of the hierarchical clustering and the graphical output that can be
obtained with Momocs on the matrix of coefficients calculated on the bot dataset.

R> botFg <- meanShapes(botF)

R> [Link](botFg$beer, botFg$whisky)
R> [Link](botFg$beer, botFg$whisky, amp = 2, [Link] = 500,
+ palette = [Link])
R> [Link](botFg$beer, botFg$whisky, [Link] = 2000, amp = 2)

4. Perspectives for Momocs

The R package Momocs is generic enough to become a gateway in analysis of outline variation
outside the developmental biology and evolutionary biology for which it has been developed.
Complementary techniques can easily be included, by the package’s developers upon request
or by third parties who have developed new mathematical approaches. Momocs stands for
modern morphometrics but, so far, only deals with outline analysis. For instance, we plan
alongside package updates, a better digitalization step including outline acquisition with the
help of Bezier curves, complementary approaches of outline sampling such as local over-
sampling and additional smoothing algorithms, the integration of 3D algorithms and more
gateways between Momocs and others morphometric programs.
Moreover, morphometric data, extracted as lists of coordinates are scarcely available on the
web, but may be very useful for meta-analysis, development, and as support for course ma-
terial. Momocs contains such datasets kindly provided by its users. We believe in such open
20 Momocs: Outline Analysis Using R

Figure 12: Some examples of the visualization of the deformation grids obtained using thin
plate splines in Momocs. From top to bottom: a simple deformation grid; isodeformation
lines; and a vector field that all depict the bindings required to pass from the average shapes
for beer and whisky.

data philosophy, as encouraged by the R Core Team and the many package developers, and
hope that Momocs will become a hub for such activity.

Acknowledgments
We thank the reviewers for all of their work and valuable comments made during the review
process. We warmly thank Allan Bailur for his invaluable help with the presentation in
English, along with Anna Seastrand on earlier versions of this manuscript. ANR-BDIV-07-008
supported VB and CG for this study; ANR-11-CEPL-002 and ANR-11-BCV-008 supported
JC for this study.
 Journal of Statistical Software 21

References

Adams D, Otarola-Castillo E (2012). geomorph: Geometric Morphometric Analysis of 2d/3d

Landmark Data. R package version 1.0, URL [Link]
geomorph.

Adams DC, Rohlf FJ, Slice DE (2004). “Geometric Morphometrics: Ten Years of Progress
Following the ‘Revolution’.” Italian Journal of Zoology, 71(1), 5–16.

Bonhomme V (2012). “A Graphical Introduction to Momocs and Outline Analysis Using R.”
Version 2012-10-17, URL [Link]

Bonhomme V, Frelat R, Gaucherel C (2013a). “Application of Elliptical Fourier Analysis to

Watershed Boundaries: A Case Study in Haiti.” Geomorphologie – Quantitative Hydro-
Geomorphology, 1, 17–26.

Bonhomme V, Prasad S, Gaucherel C (2013b). “Intraspecific Variability of Pollen Morphology

as Revealed by Elliptic Fourier Analysis.” Plant Systematics and Evolution, 299(5), 811–
816.

Bookstein FL (1991). Morphometric Tools for Landmark Data: Geometry and Biology. Cam-
bridge University Press.

Bowman A (2009). “Review of ‘Morphometrics with R’.” Journal of Statistical Software, Book
Reviews, 31(1), 1–2. Book author: Julien Claude, URL [Link]
b01/.

Chambers JM (1998). Programming with Data: A Guide to the S Language. Springer-Verlag.

Claude J (2008). Morphometrics with R. Springer-Verlag, New York.

Claude J (2013). “Log-Shape Ratios, Procrustes Superimposition, Elliptic Fourier Analysis:

Three Worked Examples in R.” Hystrix, the Italian Journal of Mammalogy, 24(1), 94–102.

Courtiol A, Ferdy JB, Godelle B, Raymond M, Claude J (2010). “Height and Body Mass
Influence on Human Body Outlines: a Quantitative Approach Using an Elliptic Fourier
Analysis.” American Journal of Physical Anthropology, 142(1), 22–29.

Crampton JS (1995). “Elliptical Fourier Shape Analysis of Fossil Bivalves: Some Practical
Considerations.” Lethaia: An International Journal of Palaeontology and Stratigraphy,
28(2), 179–186.

Dray S, Dufour AB (2007). “The ade4 Package: Implementing the Duality Diagram for
Ecologists.” Journal of Statistical Software, 22(4), 1–20. URL [Link]
org/v22/i04/.

Dryden I (2012). shapes: Statistical Shape Analysis. R package version 1.1-6, URL http:
//[Link]/package=shapes.

Ferson S, Rohlf FJ, Koehn RK (1985). “Measuring Shape Variation of Two-Dimensional

Outlines.” Systematic Biology, 34(1), 59–68.
22 Momocs: Outline Analysis Using R

Friess M, Baylac M (2003). “Exploring Artificial Cranial Deformation Using Elliptic Fourier
Analysis of Procrustes Aligned Outlines.” American Journal of Physical Anthropology,
122(1), 11–22.

Genolini C (2008). A (Not So) Short Introduction to S4. URL [Link]

org/doc/contrib/[Link].

Giardina CR, Kuhl FP (1977). “Accuracy of Curve Approximation by Harmonically Related

Vectors with Elliptical Loci.” Computer Graphics and Image Processing, 6(3), 277–285.

Iwata H, Nesumi H, Ninomiya S, Takano Y, Ukai Y (2002). “Diallel Analysis of Leaf Shape
Variations of Citrus Varieties Based on Elliptic Fourier Descriptors.” Breeding Science,
52(2), 89–94.

Iwata H, Niikura S, Matsuura S, Takano Y, Ukai Y (1998). “Evaluation of Variation of Root

Shape of Japanese Radish (Raphanus Sativus L.) Based on Image Analysis Using Elliptic
Fourier Descriptors.” Euphytica, 102(2), 143–149.

Kendall D (1989). “A Survey of the Statistical Theory of Shape.” Statistical Science, 4(2),
81–120.

Klingenberg CP (2011). “MorphoJ: An Integrated Software Package for Geometric Morpho-

metrics.” Molecular Ecology Resources, 11(2), 353–357.

Kuhl FP, Giardina CR (1982). “Elliptic Fourier Features of a Closed Contour.” Computer
Graphics and Image Processing, 18(3), 236–258.

Macleod N (1999). “Generalizing and Extending the Eigenshape Method of Shape Space
Visualization and Analysis.” Paleobiology, 25(1), 107–138.

Miller HJ, Wentz EA (2003). “Representation and Spatial Analysis in Geographic Information
Systems.” The Annals of the Association of American Geographers, 93(3), 574–594.

Moellering H, Rayner JN (1981). “The Harmonic Analysis of Spatial Shapes Using Dual Axis
Fourier Shape Analysis (DAFSA).” Geographical Analysis, 13(1), 64–77.

Paradis E (2012). Analysis of Phylogenetics and Evolution with R. 2nd edition. Springer-
Verlag, New York.

Paradis E, Claude J, Strimmer K (2004). “ape: Analyses of Phylogenetics and Evolution in

R Language.” Bioinformatics, 20(2), 289–290.

R Core Team (2013). R: A Language and Environment for Statistical Computing. R Founda-
tion for Statistical Computing, Vienna, Austria. URL [Link]

Renaud S, Michaux JR (2003). “Adaptive Latitudinal Trends in the Mandible Shape of

Apodemus Wood Mice.” Journal of Biogeography, 30(10), 1617–1628.

Richtsmeier JT, Burke Deleon V, Lele SR (2002). “The Promise of Geometric Morphometrics.”
American Journal of Physical Anthropology, 119(S35), 63–91.

Richtsmeier JT, Cheverud JM, Lele S (1992). “Advances in Anthropological Morphometrics.”

Annual Review of Anthropology, 21, 283–305.
 Journal of Statistical Software 23

Rohlf F, Archie J (1984). “A Comparison of Fourier Methods for the Description of Wing
Shape in Mosquitoes (Diptera: Culicidae).” Systematic Biology, 33(3), 302–317.

Rohlf FJ (1990). “An Overview of Image Processing and Analysis Techniques for Morpho-
metrics.” In Proceedings of the Michigan Morphometrics Workshop, Special Publication
Number 2, pp. 47–60. University of Michigan Museum of Zoology, Ann Arbor.

Rohlf FJ, Marcus LF (1993). “A Revolution in Morphometrics.” Trends in Ecology and

Evolution, 8(4), 129–32.

Rohlf FJ, Slice D (1990). “Extensions of the Procrustes Method for the Optimal Superimpo-
sition of Landmarks.” Systematic Zoology, 39(1), 40–59.

Small CG (1996). The Statistical Theory of Shape. Springer-Verlag.

Thompson DW (1917). On Growth and Form. Dover.

Wentz EA (2000). “A Shape Definition for Geographic Applications Based on Edge, Elonga-
tion, and Perforation.” Geographical Analysis, 32(2), 95–112.

Yoshioka Y, Iwata H, Ohsawa R, Ninomiya S (2004). “Analysis of Petal Shape Variation of

Primula sieboldii by Elliptic Fourier Descriptors and Principal Component Analysis.” The
Annals of Botany, 94(5), 657–664.

Zahn CT, Roskies RZ (1972). “Fourier Descriptors for Plane Closed Curves.” IEEE Transac-
tions on Computers, 21(3), 269–281.

Affiliation:
Vincent Bonhomme
Department of Ecology
French Institute of Pondicherry
600 001 Pondicherry, India
E-mail: [Link]@[Link]
Current affiliation:
Centre de Bio-Archéologie et d’Ecologie
UMR 5059 CBAE
34 000 Montpellier, France

Sandrine Picq
Centre de Bio-Archéologie et d’Ecologie
UMR 5059 CBAE
34 000 Montpellier, France
Current affiliation:
Université du Québec en Abitibi-Témiscamingue
Rouyn-Noranda J9X 5E4, Québec, Canada
24 Momocs: Outline Analysis Using R

Julien Claude
Institut des Sciences de l’Évolution
Université de Montpellier II
34 095 Montpellier, France

Cédric Gaucherel
Department of Ecology
French Institute of Pondicherry
600 001 Pondicherry, India

Journal of Statistical Software [Link]

published by the American Statistical Association [Link]
Volume 56, Issue 13 Submitted: 2012-12-24
February 2014 Accepted: 2013-09-07