16.2.1 Fourier series

Several approaches have been suggested for use on objects that do not have landmarks. The shape of outlines in a plane can be considered using Fourier series representations (Zahn and Roskies 1972). A Fourier series decomposition of a polar representation or a tangent angle representation of the x–y coordinates around the outline leads to a new set of variables (the Fourier coefficients). In particular, for a 2D outline a polar representation of the radius r(θ) as a function of angle θ with the horizontal axis is given by:

where a_j, b_j are the Fourier coefficients. If the object has been translated and rotated into a standard reference frame, and normalized to unit size in some way, then the Fourier coefficients represent the shape of the object. The Fourier coefficients can be estimated for each outline using least squares techniques, and an important consideration is to choose p. Standard multivariate analysis on the Fourier coefficients provides a suitable way of analysing such data. There are numerous applications of this technique in the literature including in engineering and biology – for some examples see Zahn and Roskies (1972); Rohlf and Archie (1984); Johnson et al. (1985); Rohlf (1990); Stoyan and Stoyan (1994); Renaud (1995); and Lestrel (1997). This Fourier approach is most often used when there is some approximate underlying correspondence between parts of the outlines. The technique can be extended to three dimensions using spherical harmonics (Chung et al. 2007; Dette and Wiens 2009). If the outlines have discontinuities or sharp edges then a wavelet representation is more appropriate. Other techniques include that of Parui and Majumder (1983) who consider shape similarity measures for open curves. Note that open curves (with distinct start and end points) are a little simpler to describe than closed curves (where the start point is an unknown parameter that needs to be factored out).

16.2.2 Deformable template outlines

In modelling the outline of an object, Grenander and co-workers (Grenander 1994; Grenander and Miller 2007) specify a series of points around the outline connected by straight line segments. In their method, variability in the template is introduced by pre-multiplying the line segment by random scale-rotation matrices. Let us assume that the outline is described by the parameters θ^T = {θ^T₁, θ₂^T}, where θ₁ denotes a vector of similarity parameters (i.e. location, size and rotation) whereas θ₂ denotes say k landmarks of the outlines. For θ₁, we can construct a prior with centre of gravity parameters, a scale parameter and a rotation by an angle α. Conditional on θ₁, we can construct a prior distribution of θ₂ as follows. Let η_j and e_j be the template edge and object edge, respectively, in complex coordinates, j = 1, …, k, where the edges are formed from the vertices (landmark points) by taking differences e_j = v_j − v_{j − 1}, with v₀ = v_k and the edges satisfy the constraint ∑^k_{j = 1}e_j = 0. A simple model is the Gaussian Markov Random Field (MRF) of order 1

where {t_j} are random Gaussian rotations and

The process {θ₂} is an improper Gaussian MRF of second order. The prior on {θ₂} is invariant under location but not scale or rotation. Kent et al. (1996b) and Mardia (1996a) give full details and also provide a link with the snake (see Section 17.3.5).

Grenander and Manbeck (1993) use edge models in an application involving automatic testing of defective potatoes. Using real notation consider a template with edges g⁰₁, …, g_k⁰ with ∑^k_{j = 1}g_j⁰ = 0. Define the random edges g_j by

where r is global size, ψ is global rotation and ψ_j are individual rotations. We can take the distribution of r to be log normal, ψ to be uniform on (0, 2π), the joint distribution of ψ₁, …, ψ_k, to be multivariate von Mises such as (Mardia 1975):

where 0 < ψ_j < 2π, j = 1, …, k; ψ_{k + 1} = ψ₁, and β₁, β₂ > 0. Thus β₁ and β₂ are deformation parameters, β₁ allows similar neighbouring rotations whereas β₂ maintains fidelity to the template. Since ∑g_j = 0, there are constraints on {ψ_j}. Let us take r = 1 and ψ = 0. The constraints ∑g_j = 0 and ∑g⁰_j = 0 imply the following constraints on ψ₁, …, ψ_k:

where ∑g⁰_jx = ∑g⁰_jy = 0. These can be simplified further for small values of {ψ_j} leading to linear constraints on ψ_j. The field of view encloses a single position and therefore location parameters do not matter. Grenander and Manbeck (1993) also use information from a colour model in the defective potato testing application.

16.2.3 Star-shaped objects

Consider star-shaped objects in which can be represented by a radial function R(u), where R(u) > 0 and u ∈ S^{d − 1} are suitable spherical coordinates. An example of a star-shaped object suitable for such analysis is the leaf outline of Figure 16.1, which was modelled by Mardia and Qian (1995).

In general the radial function R(u) contains size, shape, location and rotation information of the surface. However, location and rotation information is not present if the object can be regarded as being in a fixed location and orientation, for example from preliminary registration. In some applications, it is of interest to restrict u to part of the sphere . The size of object the could be taken as the geometric mean

and, assuming the registration is fixed, the shape of surface is taken as X(u) = R(u)/S. A possible model for log R(u) is a Gaussian process with mean log μ(u) and covariance function .

An alternative procedure is to take the size as with shape X_A(u) = R(u)/A. A suitable model here could be a Gaussian process for R(u) (Mancham and Molchanov 1996). The two procedures will be similar for small variations about the mean μ(u).

In practice a discrete set of points will be available, and we assume that the points are identified on the surface for k values of (e.g. on regularly spaced rays from the origin). Assuming an initial registration has been determined, the collection of radii R = (R₁, …, R_k)^T measures the size and shape of the surface. Let X* = (log R₁, …, log R_k)^T be the vector of log-radii. The size of the surface is measured by:

and the shape of the surface is measured by the vector

(16.3)

A possible model for log-shape is a (singular) multivariate normal distribution, and PCA can be used to reduce the dimensionality for inference.

Example 16.1 Dryden (2005b) and Brignell et al. (2010) consider the brain surface data of Section 1.4.15. Estimates of the mean surface shape were obtained using (16.3), and shape variability was investigated using PCA. In addition surface symmetry measures were calculated, and the schizophrenia patients were found to be more symmetrical than the controls – an effect known as brain torque.   □

16.2.4 Featureless outlines

In some situations there is no sensible correspondence between outlines, for example in the sand particles of Section 1.4.10. Miller et al. (1994) develop an approach where the complex edge transformations {t_i} are assumed to have a circulant Toeplitz covariance matrix in constructing a complex multivariate normal model. The model has been applied successfully to images of mitochondria (Grenander and Miller 1994). This model with circulant symmetry has also been studied by Kent et al. (2000), who consider additional constraints to make {θ₂} invariant under the similarity transformations. Kent et al. (2000) also explore the eigenstructure of the circulant covariance matrix. Kent et al. (2000) and Hobolth et al. (2002) apply their circulant models to the sand particles of Section 1.4.10, where the mean is a circle and the variability about the mean is of particular interest. In this application the shape variability of the river particles is significantly higher than the sea particles. The circulant models build on those developed by Grenander and Miller (1994), and Hobolth et al. (2002) make an explicit connection between edge-based and vertex-based models, with a focus on Markov random field models around the outline. Studies where the objects under study are featureless (i.e. no obvious correspondence in parts of the outline) are common, for example Grenander and Manbeck (1993); Grenander and Miller (1994); Mardia et al. (1996f); Stoyan (1997); Rue and Syversveen (1998); Hobolth and Vedel Jensen (2000); Kent et al. (2000); Hobolth et al. (2003) and Gardner et al. (2005) study potatoes, mitochondria, sand grains and ceramic particles, mushrooms, cells, cell nuclei, and tumour classification. Several of the models have continuous representations using stochastic processes, and these models are not dependent on the particular discretization of the curves. Parameter estimation is often carried out using likelihood-based methods (e.g., Kent et al. 2000; Hurn et al. 2001; Hobolth et al. 2002).

16.4.1 SRVF and quotient space for size-and-shape

Srivastava et al. (2011a) have introduced a very promising method for the size-and-shape analysis of elastic curves, building on earlier work that started with Klassen et al. (2003). The method involves the analysis of curve outlines or functions and is invariant under re-parameterization of the curve.

Let f be a real valued differentiable curve function in the original space, . From Srivastava et al. (2011a) the square root velocity function (SRVF) or normalized tangent vector of f is defined as , where

and ||f(t)|| denotes the standard Euclidean norm. After taking the derivative, the q function is now invariant under translation of the original function. In the 1D functional case the domain t ∈ [0, 1] often represents ‘time’ rescaled to unit length, whereas in two and higher dimensional cases t represents the proportion of arc length along the curve.

Let f be warped by a re-parameterization γ ∈ G, that is f○γ, where γ ∈ G : [0, 1] → [0, 1] is a strictly increasing differentiable warping function. Note that γ is a diffeomorphism. The SRVF of f○γ is then given as:

using the chain rule. We would like to consider a metric that is invariant under re-parameterization transformation G. If we define the group G to be domain re-parameterization and we consider an equivalence class for q functions under G, which is denoted as [q], then we have the equivalence class [q] ∈ Q, where Q is a quotient space after removing arbitrary domain warping. Inference can be carried out directly in the quotient space Q.

An elastic distance (Srivastava et al. 2011a) defined in Q is given as:

where ||q||₂ = {∫¹₀q(t)²dt}^1/2 denotes the L² norm of q.

There are several reasons for using the q representation instead of directly working with the original curve function f. One of the key reasons is that we would like to consider a metric that is invariant under re-parameterization transformation G. The elastic metric of Srivastava et al. (2011a) satisfies this desired property,

although it is quite complicated to work with directly on the functions f₁ and f₂. However, the use of the SRVF representation simplifies the calculation of the elastic metric to an easy-to-use L² metric between the SRVFs, which is attractive both theoretically and computationally. In practice q is discretized using p points.

If q₁ can be expressed as some warped version of q₂, that is they are in the same equivalence class, then d([q₁], [q₂]) = 0 in quotient space. This elastic distance is a proper distance satisfying symmetry, non-negativity and the triangle inequality. Note that we sometimes wish to remove scale from the function or curve, and hence we can standardize so that

(16.4)

and then the SRVF measures shape difference.

In the m ≥ 2 dimensional case it is common to also require invariance under rotation of the original curve. Hence we may also wish to consider an elastic distance (Joshi et al. 2007; Srivastava et al. 2011a) defined in Q given as:

The m = 2 dimensional elastic metric for curves was first given by Younes (1998).

16.4.2 Quotient space inference

Inference can be carried out directly in the quotient space Q, and in this case the population mean is most naturally the Fréchet/Karcher mean μ_Q. Given a random sample [q₁], …, [q_n] we obtain the sample Fréchet mean by optimizing over the warps for the 1D function case (Srivastava et al. 2011b):

In addition for the m ≥ 2 dimensional case (Srivastava et al. 2011a) we also need to optimize over the rotation matrices Γ_i where

This approach can be carried out using dynamic programming for pairwise matching, then ordinary Procrustes matching for the rotation, and the sample mean is given by:

Each of the parameters is then updated in an iterative algorithm until convergence.

16.4.3 Ambient space inference

Cheng et al. (2016) introduce a Bayesian approach for analysing SRVFs, where the prior re-parameterization for γ is a Dirichlet distribution. This prior distribution is uniform when the parameters in the Dirichlet are all a = 1, and a larger value of a leads to a transformation more concentrated on (i.e. translations). In the m ≥ 2 dimensional case, we consider a Gaussian process for the difference of two vectorized q functions in a relative orientation Γ, that is {vec(q₁ − q*₂)|γ, Γ} ∼ GP, where . The matrix Γ ∈ SO(m) is a rotation matrix with parameter vector θ. If we assign a prior for rotation parameters (Eulerian angles) θ corresponding to rotation matrix Γ, then the joint posterior distribution of (γ([t]), θ), given (q₁([t]), q₂([t])) is:

where γ, θ, κ are independent a priori. Cheng et al. (2016) carried out inference by simulating from the posterior using MCMC methods. For large samples estimation is consistent, which can be shown using the Bernstein–von Mises theorem (van der Vaart 1998, p. 141). Cheng et al. (2016) provide a joint model for the SVRF and the warping function, and hence it is in the ambient space.

Example 16.2 Consider a pairwise comparison of two mouse vertebrae from the dataset in Section 1.4.1, where there are 60 landmarks/pseudo-landmarks on each outline. We use Cheng et al. (2016)’s MCMC algorithm for pairwise matching with 50 000 iterations. The q-functions are obtained by initial smoothing, and then normalized so that ||q||₂ = 1. The registration is carried out using rotation through an angle θ about the origin, and a warping function γ. The original and registered pair (using a posterior mean) are shown in Figure 16.2 and the point-wise correspondence between the curves and a point-wise 95% credibility interval for γ(t) are shown in Figure 16.3. The start point of the curve is fixed and is given by the left-most point on the curve in Figure 16.3 that has lines connecting the two bones. The narrower regions in the credibility interval correspond well with high curvature regions in the shapes. We also applied the multiple curve registration, as shown in Figure 16.4.   □

There are many applications where the SRVF methods can be used, for example Srivastava et al. (2012) develop techniques for 2D and 3D shape analyses in gait-based biometrics, action recognition, and video summarization and indexing. The shape of curves that are invariant to re-parameterization has been the subject of deep mathematical work on different types of quotient spaces, with connections to elastic metrics (Michor and Mumford 2006, 2007; Younes et al. 2008a).

The choice as to whether to use quotient space inference or ambient space inference was discussed by Cheng et al. (2016), and in practice the two methods are often similar due to a Laplace approximation to the posterior density. This situation is similar to that of Chapter 14 in the analysis of unlabelled shape, where the choice between optimization (quotient space model) or marginalization (ambient space model) needs to be considered (see Section 14.2.2).