A Survey on Non-rigid 3D Shape Analysis
Hamid Laga






















































A Survey on Non-rigid 3D Shape Analysis
Laga et al.: A Survey on Non-rigid 3D Shape Analysis






































Shape is an important physical property of natural and manmade 3D objects that characterizes their external appearances. Understanding differences between shapes and modeling the variability within and across  shape classes, hereinafter referred to as shape analysis,   are fundamental problems to many applications,  ranging from  computer vision and computer graphics to biology and medicine.  This chapter provides an overview of some of the recent techniques that studied the shape of 3D objects that undergo non-rigid deformations including bending and stretching. Recent surveys that covered some aspects such classification, retrieval, recognition, and rigid or nonrigid registration,  focused on methods that use shape descriptors.   Descriptors, however, provide abstract representations that do not enable the exploration of shape variability. In this chapter, we focus on recent techniques that  treated the shape of 3D objects as points in some high dimensional space where  paths  describe deformations.  Equipping the space with a suitable  metric enables the quantification of the range of deformations of a given shape, which in turn  enables (1) comparing and classifying 3D objects based on their shape, (2) computing smooth deformations, geodesics, between pairs of objects, and (3) modeling and exploring continuous shape variability in a collection of  3D models.  This article surveys and classifies recent developments in this field, outlines fundamental issues, discusses their potential applications in computer vision and graphics,  and highlights opportunities for future research. Our primary goal is to bridge the gap between various techniques that have been often independently proposed by different communities including mathematics and statistics, computer vision and graphics, and medical image analysis.








Shape statistics, Riemannian manifolds and metrics, elastic registration, Karcher mean, bending and stretching.




A Survey on Non-rigid 3D Shape Analysis
    Hamid Laga 
Hamid Laga  is with the School of Engineering and IT, Murdoch University, Australia and with the Phenomics and Bioinformatics Research Center, University of South Australia.


E-mail: H.Laga@murdoch.edu.au.
Manuscript received April 19, 2005; revised December 27, 2012.
    December 30, 2023
=============================================================================================================================================================================================================================================================================================





























§ INTRODUCTION



IEEEexample:BSTcontrol
Shape is an important physical property of natural and manmade objects that characterizes their external appearances. Understanding differences between shapes and modeling the variability within and across  shape classes  are fundamental problems and building blocks to many applications,  ranging from  computer vision and computer graphics to biology and medicine.  
In medicine, studying shapes of 3D anatomical structures and modeling their growth patterns are of  particular interest since many diseases can be linked to alterations of these shapes <cit.>.   In computer graphics, 3D  modeling  using low level primitives (polygons)  is a tedious and time-consuming process. With the growing availability of 3D models  via online repositories,  recent research efforts have shifted their focus towards  data-driven techniques for generating  3D shape variations from existing ones.  Data-driven techniques are also gaining momentum in the field of 3D reconstruction and modeling from images, range scans and noisy point clouds <cit.>, which are ill-posed problems. By leveraging shared information across multiple objects, data-driven techniques aim at discovering patterns, which can be  used  as  
  geometric, structural,  and semantic priors for driving the 3D reconstruction and modeling processes.  














When dealing with numbers, statistics is the science for studying  and modeling variabilities. Geometric shapes are fundamentally different. Research in this area, driven by  the computer vision, medical imaging, computer graphics and statistics communities, seeks to develop mathematical tools for  manipulating shapes in the same way we manipulate numbers.  This includes:

    	
  * Comparing the shapes of two objects,    saying whether two objects have similar shape or not, and more importantly, quantifying and localizing  the similarities and differences.  
	
  * Computing summary statistics and prototypes, such as  the mean shape, modes of variations and  higher order statistics, of a collection of 3D models. 
	
  * Mathematical modeling of shape variations,    using, for example, probability distributions as generative models. 
		These models form priors for random sampling and for  statistical inferences, and 
	
  * Exploring  the shape variations for synthesizing  valid   shapes and performing  interpolations,   extrapolations, statistical inferences, and regressions.

Implementing these ideas, for 3D models,  requires solving several challenges, each one has been the subject of important research and contributions. 
The first one is about the representation  of shape.  When selecting any representation, it is important to  ensure that two different shapes cannot have the same representation and that a given representation can always be associated to a shape.

Second, almost any shape analysis task requires some measure of dissimilarity, hereinafter referred to as a metric,  for comparing shapes. The challenge when designing a suitable metric comes from the fact that 3D objects have arbitrary pose, scale, and parameterization. They also bend, stretch, and change their topology and structure.  Some of these variations, translation,  pose, scale, and parameterization, are shape-preserving, while others  affect the shape. Thus, the metric  should  be invariant to transformations that are shape-preserving  (translation, scale, rotation and re-parameterization) and should quantify  those that change the shape (bending, stretching, and structural changes). 


The third difficulty in  the analysis of 3D shapes  is registration, matching  points across objects.  Any shape metric requires a registration component to help decide which point on one object is compared to which point on the other.   This  is often a difficult problem to solve, especially across objects with large  deformations and pose variability. Thus, it is not surprising that a large body of literature is specifically dedicated to this problem, see for example <cit.> for a detailed survey of this topic.




Finally, there are several computation and numerical challenges that one needs to address. For instance, several shape representations are non-linear. This results in shape spaces and metrics that are non-Euclidean. As a consequence, standard techniques, such as Principal Component Analysis (PCA), for statistical analysis cannot be directly applied. Moreover, analytical solutions, for correspondence and geodesics, often  do not exist and one has to rely on optimizing non-convex potential functions. This makes the problem even challenging and one often has to compromise between accuracy and computational efficiency.




The purpose of this chapter is to provide a comprehensive survey of the state-of-the-art techniques in elastic 3D shape analysis. By elastic we refer to 3D shapes that undergo  elastic deformations, those that bend and stretch. This already forms a large class of natural 3D objects  with a significant large body of literature and techniques dedicated to the topic.  We structure the survey into four main sections.  Section <ref> provides   a general formulation of the problem and some background materials on which the subsequent topics are built upon.  
Section <ref> focuses on the concepts of shape spaces and metrics and their relation to physical deformations.
Section <ref> treats the elastic correspondence / registration problem and its closely related problem of computing geodesics between pairs of shapes. Section <ref>  discusses how these building blocks are used for computing mean shapes, studying the modes of variations and characterizing shape collections with statistical models. Section <ref> discusses a few concrete examples where these concepts are used. Finally, we conclude in Section <ref> by summarizing the main lessons learned in this survey and discussing future directions for research.















§ GENERAL FORMULATION



We provide in this section a general formulation of the elastic shape analysis problem. We first overview the basic representations of 3D objects, which are essential for the subsequent sections of this chapter. Then, we formulate the elastic shape analysis problem and give the taxonomy that we will use for classifying state-of-the-art techniques.





 §.§ Representations



Commonly, shape is interpreted as the boundary of a 3D object. Mathematically, it can be defined as a function f: D →^d, which maps every   point of the domain D to  a point in ^d. Here d is equal to 1 or 3 depending on the choice of the representation. 

A common choice of the function f, and thus the domain D, is the signed distance function, f: D ⊂^3 →, whose values at a point s ∈ D is the signed distance of that point to the nearest point on the boundary of the object being represented. Positive (respectively negative) values of f correspond to regions in space that are outside (respectively inside) the object.  Subsequently, points s such that f(s) = 0 form the boundary of the object. This representation is commonly known as level-set representation since for any  α∈, the points s such that f(s) = α form a surface called level-set.   This  representation is implicit since one has to solve the equation f(s) = 0 in order to identify the points that belong to the object of interest.


In contrast to implicit representations, one can represent the boundaries of an object using an explicit function f: D →^3, which assigns to every point s ∈ D a three-dimensional point f(s) = (x(s), y(s), z(s)). The set of points f(s) then forms the 3D object, which is often referred to as an embedding of the canonical domain D. Although this seems to be the most natural representation, the choice of the domain D remains an open challenge since there is no single domain D that is suitable for every type of objects. For instance, open surfaces such as 3D human faces can be easily embedded on a planar domain D ⊂^2. D can be, for example,  a unit disk. Similarly, closed surfaces of genus-0 can be embedded onto a unit sphere . 



But how is it done in practice? A 3D model is, in general,  produced either by scanning a real 3D object, using the various types of 3D sensors that are currently available, or by using a 3D modeling software. In any case, a 3D model is represented as an 
unorganized point cloud or a uniformly-sampled set of points. To facilitate rendering and visualization processes, the point clouds are converted into polygonal meshes since many rendering algorithms as well as graphics processing units are optimized for processing  planar, particularly triangular, polygons. As such, the boundary f of a 3D object is represented as a piecewise linear surface  in the form of a  triangulated mesh M, a set of vertices P = {p_1, ⋯, p_n} and triangles 𝒯 = {T_1, ⋯, T_m} such that the triangles  intersect only at common vertices or edges <cit.>.  


Finally, the parameterization process takes a triangulated mesh M, finds a suitable continuous domain , referred to as the parameterization domain,  and then estimates the underlying continuous surface f: →. If M is a disk-like surface, then  can be chosen as a simply-connected region  of the plane, the unit disk.  When dealing with closed genus-0 surfaces then the unit sphere  is the natural parameterization domain. We refer the reader to the paper of Floater et al. <cit.> for a detailed survey of different parameterization techniques.

In what follows, we assume, unless explicitly specified,  that = and thus s = (u, v) is composed of the longitude and latitude coordinates.



 §.§ Invariance requirements



An important challenge, when designing a framework for shape analysis,  is to discard effects or variables  that do not affect the shape of an object and  to account only for those that affect shape.  
Kendall <cit.> described shape as a property of an object once its rotation, translation,  scale, are parameterization are removed from the representation.
In the field of shape analysis, the transformations associated with translation, rotation, and uniform scaling are termed shape-preserving. They are nuisance variables that should be discarded. Translation and scaling are probably the easiest ones to deal with.  For instance, let  be the space of all surfaces f: →^3.  One can discard translation by first translating f so that its center of mass will be located at the origin: 
f(s) → f(s) - ∫_a(s)f(s)/∫_a(s).
Here a(s) is the local surface area at s. It is defined as the norm of the normal vector to the surface at s. In the case of spherically parameterized surfaces, i.e. = then s=(u, v) and n(s)= ∂f/∂u×∂f/∂v. The scale component can be also discarded by normalizing the surface f in such a way that its overall surface area is one:
f(s) →f(s)/√(∫_D a(s)   ).
The space   of all origin-centered and unit-surface-area surfaces is called pre-shape space.

In addition to scale and translation, 3D objects undergo two other shape-preserving transformations, which are rotation and re-parameterization. Rotations, denoted by O, are elements of SO(3), which are  3× 3 matrices. They transform each point f(s) ∈^3 into Of(s). Here, we denote by  Of the rotated surface. 

Re-parameterization is a diffeomorphism γ: →, which transforms a surface f into f∘γ. Let Γ denote the space of all such diffeomorphisms. Diffeomorphisms are shape preserving transformations; that is, the surfaces  f and f∘γ have the same shape, albeit having different parameterizations. Note that re-parameterization is important in 3D shape analysis because it provides registration. Consider two surfaces f_1 and f_2 where f_1 is, for example, the surface of a cat  while f_2 is the surface of a horse, eventually in a different pose, see Figure <ref>. Let s  ∈ such that f_1(s) corresponds to the nose tip of f_1. If  f_1 and f_2 are arbitrarily parameterized,  which is often the case since they  have been parameterized independently from each other, f_2(s) may refer to any other location on f_2, the ear tip for example. Thus f_1 and f_2 are not in correct correspondence.  Putting f_1 and f_2 in correspondence is equivalent to reparameterizing f_2 with a diffeomorphism ∈ such that for every s ∈, f_1(s) and f_2((s)) = (f_2 ∘)(s) point to the same feature, nose tip to nose tip, etc..





The simplest way of removing the rotation variability is to make sure that all shapes being analyzed are  aligned, in a consistent fashion, to a common coordinate system. One way to do this is to compute the principal axis of the objects being analyzed, using Principal Component Analysis (PCA), and align these axis to the axis of the coordinate frame. This approach, however, is very sensitive to outliers and is sub-optimal especially in the presence of large elastic deformations. 


Alternatively, invariances to rotations and re-parameterizations can be dealt with  algebraically. The idea, in the case of rotations for example, is that a surface f and any other surface obtained by rotating f have the same shape. That is, ∀∈ SO(3), f and f are equivalent. Similarly, a surface f and any of its possible reparameterizations f ∘ are equivalent.   Thus the set [f] = { f ∘, ∈ SO(3), ∈} forms an equivalence class under the action of the rotation and reparameterization groups.  The set =  / (SO(3)×) of all equivalence classes is called shape space.









 §.§ Problem statement and taxonomy
 

We are given a collection  of origin-centered and scale-normalized 3D objects {f_i ∈, i=1, ⋯, n}. For the generality of formulation, we assume that each object f is represented as an embedding: f: →^3, which assigns to every point s ∈ a 3D point f(s) ∈^3.  When talking about shape analysis,  we are interested in (1)  quantifying similarities and differences in their shape, and (2) characterizing shape variability in the collection, which in turn can be done by computing summary statistics such as mean shape and modes of variation, and fitting probability distributions to the population in the same way as it is normally done with numbers. 


Intuitively, the mean shape  can be defined as the object whose shape is as close (similar) as possible to all the shapes of the other objects in the collection. This is known as the Karcher mean. Mathematically,  we seek to find a shape  and a set of optimal rotations ^*_i and diffeomorphisms ^*_i, i=1, ⋯, n that align every surface f_i onto the computed mean shape :

    (, ^*_i, ^*_i )  = _f, _i, _i∑_i=1^n d_(f, _i f_i ∘_i),
  
where d_(·, ·) is a certain measure of  dissimilarity (or distance) that quantifies shape differences. This equation can be simplified by writing it in terms of orbits or elements of the shape space  as follows:

    []= _[f] ∑_i=1^n d_([f],  [f_i]).
  
where

    d_([f_1],  [f_2]) = min_,  d_(f_1,  f_2 ∘).



Similarly, one can define  the variance as  the expected value of the squared distance from the mean <cit.>:

    σ^2 = 1/n∑_i=1^n d_([], [f_i])^2.

Implementing this simple general formulation requires solving the following challenging fundamental problems;

	
  * Choice of the metric d_.  Computing shape statistics depends on the choice of the metric, or measure dissimilarity d_. The challenge is in designing a metric that captures shape differences while being invariant to the set of transformations that preserve shape. 
	
  * The registration or correspondence problem.  Naturally, the solution to Eqn. (<ref>) depends on the way the input surfaces are parameterized, or registered, the process of finding the optimal rotations O_i and diffeomorphisms _i. 
	
  * Computational issues  that arise when using complex metrics that are physically motivated. One often deals with large collections of 3D models and thus the process of evaluating the metric, computing geodesics and shape statistics should be computationally efficient.

The following sections discuss these challenging problems and present different solutions that have been proposed in the literature.








§ SHAPE SPACES AND METRICS




In what follows, we represent the boundary of a 3D object as an embedding f: →^3. Here, we assume that the parameterization domain is the unit sphere.     It maps every point s = (u, v) ∈ to a point f(s) ∈^3.  For clarity purposes, we assume at this stage that the given surfaces have been normalized for translation and scale, thus they are elements of .



In shape analysis, a   dissimilarity measure quantifies  the amount of deformations, or energy, that one needs to apply to one shape in order to align it onto the other.  Consider the example of Fig. <ref> where f_1 is a straight cylinder and f_2 is a bended one. The deformation of f_1 onto f_2 results in a sequence of m intermediate  shapes F_i: →^3, i=1, …, m such that F_0 =  f_1,  F_m = f_1, and F_i = F_i-1 + V_i-1.
 Here, V_i-1: →^3 is a (deformation) vector field such that F_i(s) = F_i-1(s) + V_i-1(s).  The sequence F = {F_i, i=1,⋯, m} can be seen as a path or a curve in the preshape space . Its length is given by:

    L(F) = ∑_i=1^m√(( F_i-1,  F_i)),
 
where  is a certain measure of distance between the surfaces F_i-1 and F_i. When m is sufficiently large, F can be seen as a continuous deformation of f_1 onto f_2. It can be interpreted as a parameterized path F: [0, 1] →  such that F(0) = f_1, F(1) = f_2. Its length is given by:

    d_ (f_1, f_2) = L(F) = ∫_0^1  d F/dτ(τ)   d F/dτ(τ) ^1/2 dτ.

Here, d F/dτ(τ) is in fact an infinitisemal  vector field that deforms the surface F(τ). It can be also interpreted as a small perturbation of the surface.  The inner product, ··, also known as the metric, measures the strength  of this vector field.  Thus, integrating over τ provides the length of the path.

In general,  a surface f is treated as an element, or a point, in a pre-shape space . A vector field that deforms f is then a vector δ f that is tangent to  at f. Let T_f() denote the tangent space to  at f.  A metric is an inner product on T_f(). It takes two tangent vectors v_1 and v_2 ∈ T_f() and returns their inner product v_1v_2. The norm of a vector v ∈ T_f(), according to the metric ··, is given by vv^1/2. Since there are many paths that deform f_1 onto f_2, we are particularly interested in the shortest one, under the metric, which is called the geodesic in the pre-shape space :

    F^* = min_FL(F).
 
The length of the geodesic path is the dissimilarity, in the preshape space, between the two surfaces:

    d_ (f_1, f_2) = inf L(F) = L(F^*).

Since almost every 3D shape analysis task includes a step in which shapes are compared based on some measure of dissimilarity, the choice of the metric is very critical for the subsequent analysis tasks.  







 §.§ Kendall's shape space




The fundamentals  of  statistical analysis of shapes have been laid by Kendall <cit.>, as early as 1977, and advanced by many others <cit.>.  The basic idea is to   represent  the shape of a 3D object using  a finite set  p of  ordered points, called landmarks,  p = {p_i ∈ℝ^3, i=1, …}.  These points are sampled from the shape boundary and put in correspondence across shapes.  

In the case of our setup where a surface f is treated as a mapping from a parameterization domain  to ^3, Kendall's representation is equivalent to having a fixed finite sampling of the parameterization domain . That is D = {s_i, i=1, …, }. As such, the space  of all surfaces that can be represented with  landmarks is ^3.  The pre-shape space, , the space of surfaces that are normalized for translation and scale, is a hypersphere ^3-4, which is a subspace of ^3.  To make the representation invariant to rotations, Kendall  defines the group action SO(3) ×^3-4, (O, f) = Of whose orbits are  given by [f] = {(O, f) | O∈ SO(3)}. The  shape space  can then be viewed as  the quotient space ^3n-4 / SO(3).  The rotation-invariant geodesic distance between two surfaces f_1 and f_2 is given as the geodesic distance between their orbits on , :

    d_([f_1], [f_2]) = min_O ∈ SO(3) d_(Of_1, f_2) = min_O ∈ SO(3) d_(f_1, Of_2),

where d_ is the length of the arc on  that  connects f_1 to Of_2.  A geodesic path is then a great circle on    between  f_1 and O^*f_2 where
O^* = _O ∈ SO(3) d_(f_1, Of_2).
Kendall's shape space,  introduced in 1977, and its subsequent developments <cit.>, are considered as the foundation that lead to the major recent developments in statistical shape analysis.  








  §.§.§ Morphable models





In the case of planar objects, simply treating  Kendall's representation as a vector space, assuming that the preshape space  is Euclidean, and ·· is  an  metric,  results in the  2D  Active Shape Models (ASM)  of  Cootes et al. <cit.>.  This formulation has  been later used for the analysis of the shape of 3D objects that are in static poses, leading to what is now known as  morphable models. Morphable models have been  originally introduced for the analysis of the 3D shape of human faces <cit.>. They have been later used for the analysis of the 3D shape of  human bodies <cit.> as well as the shape of objects originating from various domains including  archaeology, astronomy, morphometrics,  and medical diagnosis <cit.>.

Mathematically, for a given pair of vectors v_1 and v_2 ∈ T_f(),  the  metric is defined as:

    v_1v_2 = v_1v_2,

where ·· is the standard inner product. The geometric interpretation of this metric is as follows; A surface f_1 can be optimally deformed to another surface f_2 by simply adding to each point f_1(s) a displacement vector (s) such that
f_2(s) = f_1(s) + (s).
The optimal path F^* that connects f_1 to f_2, induced by the metric of Eqn. <ref>,  is simply a straight line. Subsequently, the difference, or dissimilarity,  between f_1 and f_2 is given by 

    d(f_1, f_2) = ∫_ds  =   ∫_f_1(s) - f_2(s) ds.
 
Assuming that d(·, ·) is Euclidean, and under a fixed parameterization, finding the optimal rotation that aligns one surface onto another is equivalent to finding the  rotation ∈ SO(3) that minimizes f_1 - O f_2.  
This   results in the Procrustean metric on the shape space , which can be efficiently solved using Singular Value Decomposition (SVD) <cit.>. Now,  for a fixed rotation and reparameterization, the surfaces are already registered,  both the mean shape of Eqn. (<ref>) as well as the covariance matrix   have analytical solutions:

    = 1/∑_i=1^ f_i,   and  = 1/-1∑_i=1^ (f - )×  (f - )^T.

Thus, efficient parameterization of the shape variability in a collection of 3D shapes can  be performed using Principal Component Analysis (PCA). Let λ_k, k=1, …, d, be the leading eigenvalues of  K,  and v_k their corresponding eigenvectors, called modes of variation. A new statistically feasible shape f  is given by

f =     + ∑_i=1^dα_k v_k,  where α_k ∈ are coefficients that control the contribution of each mode of variation.  Arbitrary new shapes can be generated by varying the coefficient vector  α = (α_1, α_2, …, α_d) whose  probability follows a Gaussian distribution  defined as
-logPr(α ) = 1/2∑_i=1^dα_i^2/λ_i + const.  Thus, if one wants to synthesize a plausibe shape, a shape that is similar, but not identical, to those in the collection, one only has to generate arbitrary αs and retain those that have high probability.

This approach is known as morphable models and is an  extension of the Active Shape Models (ASM) developed by    Cootes et al. <cit.> for the analysis of planar objects.
 They have been first introduced to 3D shape analysis by Blanz and Vetter in their seminal work on the synthesis of 3D faces <cit.>.  Morphable models have   been later generalized to characterize the space of the entire human body shapes <cit.>.  Both in <cit.> and <cit.>, the objects under study are assumed to be in neutral poses and their landmarks are in a one-to-one correspondence. 
















  §.§.§ The non-linear nature of Kendall's shape space








Although the  metric on Kendall's shape space has been extensively used in a wide range of applications, archaeology, astronomy, morphometrics, medical diagnosis <cit.>, it is only suitable for computing geodesics, and thus statistics, between shapes that undergo small deformations, anatomical shapes,  3D faces, or human shapes in a neutral pose.   This main limitation is illustrated in Figure <ref> using two synthetic cylinders.  The example shows optimal deformation paths between surfaces f_1 and f_2, where f_1 is a straight cylinder and f_2 is a bent cylinder. Figure <ref>-(a) shows the linear path between f_1 and f_2 in  (after full registration) obtained using the  metric by connecting  each pair of corresponding points with a straight line. The intermediate shapes along this path shrink unnaturally. Figure <ref>-(b) shows a natural  geodesic path computed with a metric that captures the shape deformations (here using SRNF inversion, which will be described later). 

 
Kilian et al. <cit.>  build on Kendall's representation a new metric that does not suffer from the shrinkage issue of the  metric. The idea is to define the distance between a surface f_1 and f_2 as the strength of a vector field u ∈ T_f_1() that penalizes elastic deformations, enforcing as rigid-as-possible and as-isometric-as-possible deformations along geodesic paths.   That is, given a shape f ∈ = ℝ^3n,  define  the inner product between two vectors u and v in T_f()  as:

    uv =  uv^R + uv^I + λ∑_s=1^nu(s)v(s)a(s),

where λ = 0.001 and a(s) is the local Voronoi area at the vertex s. The first term penalizes non-rigid deformations. To define it, consider first that the deformation vector u can be written as u= u' + u” where u' is only composed of rigid transformations and u” is only composed of non-rigid deformations. To have a metric that is invariant to rigid transformations while quantifying non-rigid ones then one can set uu^R  =  u”u”. As such,  uu^R = 0 if and only if u is a rigid motion. This term  ensures invariance to translation, scale and rotation without normalization in  a pre-processing step.    

The second term of Eqn. (<ref>) favours bending over stretching.  It is defined as the strength of the deformations that remain after discarding isometric motions, motions that preserve geodesic distances along the surface of f.  In other words, it can be seen as a measure of stretch in terms changes in the length of curves along the surface f. Since Kilian et al. <cit.>'s approach works on triangulated meshes then

    uu^I = ∑_(s, t)u(s) - u(t)  f(s) - f(t)  ^2,
 
where s and t are two adjacent vertices on the triangulated mesh that represents the surface f.   Thus, uu^I = 0 if and only if  u just bends the shape f. In other terms, when deforming the straight cylinder of Figure <ref>, the metric will favor bending over stretching and thus it will minimize the shrinkage that is observed when using the standard  metric. 


The last term is a regularizer that makes the metric Riemannian.  It also ensures that the geodesic  between two isometric shapes is unique.  The main limitation of this metric is the choice of the parameter λ in order to  ensures that the metric is Riemannian. Large values of λ will give high weight to the Euclidean term of Eqn. (<ref>). Thus,  the metric may suffer  from the shrinkage issue observed in morphable models (see Figure <ref>). Note that Kilian et al. empirically set this parameter to 0.001. 








 §.§ Metrics that capture physical deformations



Instead of using the  metric on , one would like to explicitly capture the type of deformations that one needs to apply to f_1 in order  to align it onto f_2. Such deformations can be of two types; bending and stretching.  This requires redefining Eqn. (<ref>) in terms of an energy function that penalizes these deformations:

    d_^2(f_1, f_2)   =   (f_1, f_2) 
       =    ∫_0^1  {α_s ( F(τ) )   +    α_b ( F(τ)  ) } d τ,

where _s and _b are respectively the stretching and bending energies. The parameters α_s and β_s ∈ control the contribution of each of the terms to the energy function. Below we discuss two important pieces of work that implement this model. 

The first one (section <ref>) is inspired by the elasticity theory in physics <cit.> where surfaces are treated as thin shells, a thin three-dimensional material of thickness δ. Stretching in this case is caused by in-layer (tangential) shear or compression while bending is caused by friction due to transversal shear <cit.>. The second one (Section <ref>) showed that with specific choices of the bending and stretching terms as well as their weights, the deformation model of Eqn. (<ref>) reduces to an  metric in the space of Square Root Normal Fields (SRNF), a special representation of surfaces <cit.>. 



  §.§.§ The shape space of thin shells



Based on a physical deformation paradigm, a surface f can be viewed as the middle layer of a thin viscous material of some finite (albeit small) thickness δ <cit.>.    Wirth et al. <cit.>, Heeren et al. <cit.>,  Berkels et al. <cit.>, and Zhang et al. <cit.>  used the deformation model of Eqn. (<ref>), with α_s = 1 and α_b proportional to the squared thickness δ^2. The stretch energy, referred to as the membrane term, is given as (f_1, f_2) = ∫_((s)) ds, where
	
    ((s) )    =   μ/2((s)) + λ/4((s)) - 
          (  μ/2  +  λ/4)  log  ((s)) - μ -λ/4.
	
The log term penalizes material compression. λ and μ are the Lamé constants of the tangential Newtonian dissipation measure.   The Cauchy-Green strain tensor,  =   _1^-1_2, where g_i is the first fundamental form or metric of the surface f_i,  accounts for changes in the metric  between the source and target surfaces.  In particular, the trace of the Cauchy-Green strain tensor, (G),  accounts for local changes in curve lengths while the determinant (G) accounts for local changes in surface area <cit.>.  Thus,  this can be seen as the generalization of the stretch metric of Kilian et al. <cit.>, which only considered changes in the length of curves. 

The bending energy should account for out-of-plane bending and changes in curvature. Since such information is encoded in the shape operator, then the most natural choice of measure for quantifying bending is variation in the shape operator <cit.>:

    ( f_1, f_2  ) =  ∫__2(s)) -  _1(s) _F^2 ds,

where ·_F denotes the Frobenius norm. The energy of Eqn. (<ref>) takes into account  full change in the second fundamental form.  Heeren et al. <cit.> and Windheuser et al. <cit.> simplified this metric by only considering changes in the mean curvature, which is the trace of the shape operator ;
	
    ( f_1, f_2  ) =  ∫_(_2(s)) -  (_1(s)) ^2 ds.

This energy is also known as the Willmore energy. In the discrete setup, Heeren et al. <cit.> approximated the mean curvature by the dihedral angle, the angle between the normals of two adjacent faces on the triangulated mesh.


Since this deformation model acts on the Cauchy-Green tensor  and mean curvature  then one can define the preshape space as the space of all pairs  (,   ).   An important observation that has been made in <cit.> is that the Hessian of the elastic deformation energy results in a proper Riemannian metric on the space of shells, modulo rigid body motions. That is, for an infinitesimal deformation  v ∈ T_(f), ()(v, v) = D^2 (v, v) is a proper metric that measures the strength of this deformation. (Here, D^2 (v, v) is the second derivative of  in the direction of v.)   Moreover,  Heeren et al. <cit.> showed that ()(v, v)  = 0 if and only if v induces an infinitesimal rigid motion. This last property is highly desirable in shape analysis where, as described in Section <ref>, one seeks a framework that is invariant to rigid motions of shapes. This properties guarantees this invariance without normalization for translation and rotation of the shapes being analyzed. Shapes, however, need to be normalized for scale. 



















  §.§.§ The shape space of square-root representations







Instead of using the Hessian of the deformation energy, Jermyn et al. <cit.> showed that by representing a surface f using the pair (,  ),  one can define a full elastic metric on the product space of metrics (first fundamental forms)  and unit normal fields   as follows:

    (δ, δ) (δ g, δ)  = ∫_ ds  √( |g| ){α(  (g^-1δ   g)^2   )  + 
    β(g^-1δ   g )^2      +   μδñ(s)δñ(s) }.

The first term of Eqn. (<ref>)   accounts for changes in the shape of a local patch that preserve patch area. The second term quantifies changes in the area of the patches. Thus, the first two terms capture full stretch. The last term measures changes in the direction of the surface normal vectors and thus captures bending. The parameters  α, β, and μ are positive weights that one can use to penalize the different  types of deformations, if one wants to favor stretch over bending then the third term should be given a higher weight.  As such, instead of working in the original space,  and , one can define a mapping  such that (f) = (, ) where (s) is the the first fundamental form of f at s and  (s) is the unit normal vector to the surface at s.  The  preshape space  is then the product space of metrics and unit normal vectors. The shape space is the  
quotient space =  /(SO(3)×), via the action of the rotation and re-parameterization groups.


Jermyn et al. <cit.> defined a new representation of surfaces called the square-root normal fields (SRNF), which are essentially   surface normals scaled by the square root of the local area: 

    Q(f ) (s) = q(s) = (s)/√( |  (s)  | )

where (s) = ∂ f/∂ u×∂ f/∂ v and s = (u, v) ∈. The space of SRNFs, hereinafter denoted by , has very nice properties that are relevant for shape analysis. In particular, Jermyn et al. <cit.> showed that the   metric on  is a special case of the full elastic metric given by  Eqn. (<ref>). It corresponds to the case where α = 0, β = 1/4, and μ = 1:

    δ qδ q =  ∫_ ds  √( |g| ){1/4(g^-1δ   g )^2      +    δñ(s)δñ(s)}.

Since √( |g| ) is just |n|, the Eqn. (<ref>) is equivalent to

    δ qδ q =    1/4∫_ ds {δ a(s) δ a(s)/ a(s)          +     δñ(s)δñ(s) a(s)   }.

Thus, the  metric in the space of SRNFs is equivalent to a weighted sum of surface area change and surface bending. This property of SRNFs makes them very promising as a representation of surfaces for elastic shape analysis.  If geodesics, mean shapes, PCA, etc. could be computed in  under the  metric, and then mapped back to  or  then there would be large gains in computational efficiency with respect to other metrics such as those defined on the space of thin-shells. 
 Unfortunately, there is no analytical expression for Q^-1 for arbitrary points in . Moreover, the injectivity and surjectivity of Q remain to be determined, meaning that for a given q∈, there may be no f∈  such that Q(f) = q, and if such an f does exist, it may not be unique.
 
If one cannot invert the representation, one can always pull the  metric back to  under Q and perform computations there, as in  <cit.>, but this is computationally expensive, and rather defeats the purpose of having an  metric in the first place.  Although an analytical solution to the inversion problem remains an open problem, Laga et al. <cit.> showed that one can always estimate, for a given q ∈ a surface f∈ whose SRNF representation is as close as possible q. 


SRNF representation is a generalization of the Square-Root Velocity Function (SRVF) <cit.> proposed for elastic analysis of planar curves and used in many applications including biology and botany <cit.>. It is also  a special case of the family of square-root maps (SRM). Another example of SRM is the Q-maps introduced  by Kurtek et al. <cit.>, which, similar to SRNF maps,  is used for the analysis of the shapes of parameterized surfaces using the  metric.  It is defined a mapping Q: →, where q = Q(f) such that  q(s) = √(|n(s)|)f(s).  




Properties. 
SRNF and Q-maps share very nice mathematical properties that make them suitable for shape analysis. In particular, both SRNFs and Q-maps are a subset of 𝕃^2(𝕊^2, ℝ^3). This means that one can use the standard  metric in the space of SRNFs or Q-maps for the analysis of shapes instead of the complex and expensive pullback elastic metrics. Furthermore, if a surface is re-parameterized according to f → f ∘γ, γ∈Γ,  its corresponding Q-map or SRNF is given by (q, γ) = (q∘γ)√(J_γ), where J_γ is the determinant of the Jacobian of γ.  


Another  main motivation  behind the use of  Q-maps and SRNFs in surface analysis is that the action of Γ  is by isometries under the  metric, q_1 - q_2 = ( q_1, γ) - ( q_2, γ),  ∀γ∈Γ.  This property is very important for joint rotation and re-parameterization-invariant comparison, and thus elastic registration, of surfaces as will be discussed in Section <ref>. As such, Q-maps and SRNF maps have been used for optimal re-parameterization, and thus registration,  of spherical and quadrilateral surfaces <cit.> that undergo isometric as well as elastic deformations,  and surfaces that contain missing parts. They have been also used for computing geodesic paths between  such surfaces <cit.>. The strength of these two approaches is that shape comparison, registration and geodesic computations  are performed under the same elastic Riemannian metric.  

The main difference between Q-maps and SRNFs is that (1) the former lacks a physical interpretation;  it does not have a clear relationship to an underlying elastic metric. The representation was solely devised for convenience of being able to compare the shapes of parameterized surfaces using the  metric, and (2) unlike Q-maps, SRNF representation is invariant to translation since it is computed using the surface normals.







 §.§ Transformation-based representations




The approaches described in Sections <ref>  and <ref> represent a shape as an embedding f: →^3 while the metric quantifies the deformations the shape undergoes. Another approach is to model transformations that 3D objects undergo rather than the transformed objects themselves. These representations derive from Grenander's pattern theory <cit.> in which a shape is not represented as such but as a deformation of  another (fixed) shape, T,  called template <cit.>. This paradigm has been introduced by Grenander <cit.>, developed in <cit.> and applied to various problems in computational anatomy and more recently to the analysis of human body shapes <cit.>. 























Let's assume that we are given a set of registered surfaces represented as triangulated meshes with the same graph topology as the template T.  We assume that the template has m triangular faces, T={T_k, k=1… m}. Each  triangular face can be represented with a matrix of its three vertices,  [v_0, v_1, v_2] ∈ℝ^3×3, or with its edge matrix e=[v_1 - v_0, v_2 - v_0] ∈ℝ^3×2.  Any surface   X_i = {X^i_k, k=1… m} of  m triangles  can be represented with a set of transformation matrices P_i = {P^i_k, k=1… m }  that deform the template T, X^i_k = P^i_k T_k where    P^i_k  is a 3×3  matrix that encodes various types of deformations that  T^i_k can undergo.  Hereinafter, we refer to  P as the deformation model.    

 
Existing methods  differ in the choice of (1) the  template T,   (2) the deformation model P, and (3) the metric used to measure distances in the space of deformations.   We will discuss these aspects in the following subsections.





  §.§.§ Choice of the template T 
 

When dealing with 3D models of the same class, human body shapes, one can define the template T as a typical 3D model in a neutral pose. This has been used in the SCAPE model <cit.> and in many other papers <cit.> for the analysis of human body shapes.  In general, the template T can be any surface, a sphere, of same topology and triangulation as the surfaces being analyzed.  In practice, however, domain specific templates are desirable to improve robustness and computation time. Alternatively,   one can define T as a single canonical triangle in the xy-plane <cit.>.  An exception to these  representations is the SCAPE model <cit.> where, in addition to  the template mesh T, an articulated skeleton is used to capture articulated motion.  These approaches operate on  meshes that are in full correspondence with exactly the  same connectivity. 







  §.§.§ Deformation models


The deformation model P encodes  different deformations that a triangular face X^i_k∈ X_i undergoes with respect   to its corresponding triangle T_k on T.  These  can be  rotations due to articulated motion (bending an arm of a human body), stretching due to shape differences between say two individuals in the same pose, or pose induced deformations (muscle stretch due to arm  bending).  


 SCAPE <cit.> used two deformation models, one acting on a template mesh and another on an articulated skeleton (which captures articulated motion). This  greatly simplifies the mathematical formulation;   A given mesh  X_i is represented with  three sets of transformation matrices: Q^i = {Q^i_k }, S^i = {S^i_k }, and R^i = {R^i_ l[k]}  for  pose induced deformations,  body shape deformations, and articulated motions, respectively.   
  Here,  R^i acts on an articulated skeleton used to capture articulated motion where R^i_ l[k] refers to the rotation of the l-th part of the skeleton to which the k-th face belongs to. 
  
Using  separable deformation models assumes that intrinsic  shape deformations and  pose variations are independent. This, however,   prevents one  from capturing phenomena where there is a strong correlation between body pose and local deformations (near the joints or at the muscles). For example, the  surface deformation generated by the motion performed by an athletic person exhibits different properties than the same motion carried out by a person with less pronounced  muscles <cit.>. Separable models cannot capture such  correlation between pose and body shape.

To overcome this limitation, Hasler et al. <cit.>  introduced a bi-linear model of pose and shape, without using a skeleton, where   each triangular face X^i_k is defined as X^i_k = R^i_k S^i_k T. The matrices R^i_k and S^i_k, which represent the pose and shape components,  respectively, are affine transforms applied to the template  T.  Note that, while  SCAPE  requires a set of 3D shapes of the same  subject  in different poses to learn the space of pose deformations, Hasler et al. <cit.>  learn   both spaces using  a random mix of 3D models from different subjects  in various poses.  Both  methods have shown detailed  representation of the human body shape. However, using strong priors on a human body model, manually defined skeleton and body parts, makes these methods  difficult to extend to other, especially free-form, object categories without redesigning the representation <cit.>.

Another property of the representations above is that they have redundant degrees of freedom (DOF). For instance, Hasler et al. model deformations with 15 DoF and a non-linear encoding of triangle deformations that captures dependencies between pose and
shape.  Freifeld and  Black <cit.> showed that triangle deformations lie in a 6 dimensional nonlinear manifold.   More precisely, X_k^i = R^i_k O^-1_k A^i_k S^i_k  O_k T_k, where O_k is a rotation that maps the template triangle T_k to a canonical planar triangle. The scaling  S^i_k and transformation A^i_k are in-plane transformations. Finally,  R^i_k and  O^-1_k are rotations that map the planar triangle to the 3D space. Note that O_k is fixed since it depends on the template.  Freifeld and  Black also showed that R^i_k O^-1_k A^i_k S^i_k  O_k is invertible and has 6 DoF, which eliminates deformations that do not have physical meaning. A triangular mesh of m faces is  then represented with m-tuples (R, A, S), where R∈ SO(3), A is a 2×2 matrix, and S>0. 

















Finally, medial representations (M-reps), which are an extension of curve skeletons <cit.> and medial surfaces <cit.>,   have been  introduced by Blum <cit.> as a tool  for representing the geometry and structure of anatomical shapes <cit.>.   Unlike the deformation models described above, which act on the shape's surface, a medial surface is formed by the centers of all spheres that are interior to the object and tangent to the object boundary in two or more points.
Each medial point, called medial atom,   is a 4-tuple m = {x, r, n_0, n_1}, where x∈ℝ^3 is the location of the atom, r ∈ℝ^+ is the radius of the  maximum sphere that is tangent to the surface in two or more points, n_0 and n_1 ∈𝕊^2 are two unit spoke directions.  A medial atom m  is an element of  the manifold ℱ(1) = ℝ^3 ×ℝ^+ ×𝕊^2 ×𝕊^2.  Fletcher et al. <cit.> showed that by representing the atoms with the deformation, translation, scaling, and rotation, of a canonical atom, (0,0,1), they become elements of 𝒮(1) = ℝ^3 ×ℝ^+ × SO(3)/SO(2)  × SO(3)/SO(2), which is a Lie group. A grid of  medial atoms representing a 3D object is   an element of the Lie group 𝒮 = 𝒮(), the product of  copies of 𝒮(1). 


Since their introduction, M-reps have been  used   in various medical image processing tasks, including segmentation <cit.>, registration, shape discrimination <cit.>, and  brain-structure analysis <cit.>.  








  §.§.§ Metrics on the space of deformations 
  

Representing shapes with the deformations of a template leads to shape spaces that are non Euclidean  (rotations, for example, are elements of SO(3)).  Some methods   typically use a Euclidean representation of deformations and measure distances  using the  metric,  ignoring the geometry of the space of deformations. The SCAPE model <cit.>, for instance, learned a shape model using the Euclidean distance and PCA on the  deformation matrices S using only  the subjects that are in the same neutral pose. Similarly, Hasler et al. <cit.> learned, separately,  
two low-dimensional models of  shape and pose, and then combined them at a later stage. 

Unlike these works, both M-reps <cit.> and Lie bodies <cit.> exploit the manifold structure of the space of deformations where geodesics and statistics are computed using proper Riemannian metrics. An important property that makes these methods of practical interest is that, despite the non-linearity of the metric as well as the manifold,  geodesic distances and geodesic paths between a pair of shapes have a closed-form  formula, unlike other methods which require expensive optimization schemes, see <cit.> and <cit.> for the mathematical details.   This significantly simplifies the process of comparing shapes and computing geodesics and  summary statistics.

Note that these representations require   that all the shapes being analyzed to have the same connectivity and to be in one-to-one correspondence with each other and with the template.  Such correspondences need to be computed using a different approach, and more importantly, using an optimality criteria that  is often different from the one used for computing geodesics and statistics. 













§ REGISTRATION AND GEODESICS



Now that we have introduced the concepts of metrics, pre-shape spaces, and shape spaces, we turn our attention to the problem of computing geodesics, smooth optimal deformations, F^*,  between a source surface f_1 and a target surface f_2 ∈:

    F^* = min_(O, ) ∈ SO(3)×{min_F:[0, 1] →,
    F(0) = f_1, F(1) = O(f_2 ∘) L(F)  },

where L(F) is given by Eqn. (<ref>) and depends on the choice of representation and metric. This double optimization problem finds the shortest path, and thus the geodesic, in the shape space between [f_1] and [f_2]. The inner minimization problem assumes that the surfaces are in correspondence, fixes the rotation O and reparameterization γ, and finds the shortest path in  that connects f_1 to O(f_2 ∘).  The outer optimization is a registration step, which finds the optimal rotation  O and reparameterization  that bring  f_2   as close as possible to f_1 without changing its shape. 
Clearly, these two problems, registration and geodesics, are interrelated and should be solved jointly using the same metric (or optimality criteria). 
















 §.§ Registration


In shape analysis, registration quality is a major factor in determining the accuracy of the subsequent analysis steps such as computing optimal deformations (geodesics) or  summary statistics (mean shapes and modes of variations).  It has been extensively studied in the past two decades, see <cit.> and <cit.> for comprehensive surveys and different taxonomies of the state-of-the-art. In general, the registration problem can be mathematically stated as the problem of finding the optimal  transformation, T, and a re-parameterization γ: → such that T (f_2 ∘γ) is as close as possible to f_1, ∫ d^2(f_1(s),   T f_2(γ(s) ) ds is minimized. Here  d(·, ·) is a certain measure of closeness in the space of surfaces.  We have seen in Section <ref> how to remove the variabilities due to translation and scale. In this section, we will focus on the remaining variables, which are rotations and reparameterizations. Thus, we are given two surfaces f_1 and f_2 ∈ and we seek to solve the following optimization problem:

    (^*, ^*) = _(, ) ∈×  d_(f_1,    (f_2 ∘) ),


The optimization problem of Eqn. (<ref>) is a complex problem that can be solved by iterating between the optimization over rotations O∈ SO(3), finding the optimal rotation, given known correspondences (fixed ), and then finding the optimal correspondences or reparameterization  given fixed rotation O. 

If correspondences are given,  finding the best   that aligns f_2 onto f_1  is straightforward. Under the  metric,  the best rotation  can be analytically found using Singular Value Decomposition (SVD) leading to the Procrustes analysis of shapes <cit.>, which forms the basis of  2D Active Shape Models (ASM) <cit.> and  3D morphable models <cit.>. Similarly, if one is using the partial elastic metric of Jermyn et al. <cit.> (Section <ref>), then one can find the optimal rotation O in the space of SRNFs and apply it to the original surface since the SRNF map of Of is exactly Oq. 


In practice, correspondences are   unknown and better results are obtained if  one  solves simultaneously for the best correspondences and for the optimal rigid  alignment. 












  §.§.§ Landmark-based elastic registration


In this class of methods, 3D shapes are treated as sets of landmarks sampled from their surfaces. The problem of finding correspondences between two surfaces f_1 and f_2 is then reduced to the problem of matching landmarks across the surfaces using some optimality criteria.
The sampling of landmarks can be either uniform or adaptive to cope with the varying surface complexity <cit.>. For instance, in the case of triangulated meshes, each vertex can be treated as a landmark. In the case of parameterized surfaces, f:→^3, then one can sample points s_i from the domain   and use f(s_i) = f_i as landmarks. 

A popular approach for putting in correspondence the landmarks is the Iterative Closest Point (ICP) algorithm, proposed in <cit.>,  which iterates between two steps: (1) a matching step where correspondences are estimated either by nearest-neighbor search  using  point-to-point <cit.> or point-to-plane <cit.> distances, or using some local descriptors, and (2) an optimal alignment step using the estimated correspondences. Since  its introduction,  many  variants of the ICP algorithm have been proposed. They aimed at improving various aspects, such as  speed and quality of  convergence,   of the original algorithm, by, for example, combining multiple distance measures  with some local descriptors <cit.>. 
Nevertheless, ICP algorithm and its variants provide good results when the poses of the shapes being aligned are close to each other. 
When the difference in  pose  is large,  the search space for  correspondences is very large. In theory, however,  three pairs of corresponding points are sufficient to define a rigid transformation. This fact has been  used in <cit.> to derive RANSAC-based solutions. Their complexity is, however, high, of order O(^3), where  is the number of points being aligned. Aiger et al. <cit.> introduced the four point congruent sets (4PCS) algorithm, which reduced the complexity of RANSAC algorithms to O(^2). It has been later optimized to achieve a complexity of order O() <cit.>.


Another approach for finding correspondences is to characterize the shape's geometry around the landmarks using some descriptors <cit.> and then match each landmark to its closest one in terms of a distance measure in the descriptor space <cit.>. The main difficulty, however, is that  in the vast majority of applications, one is required to put in correspondence  two or more 3D objects that undergo non-rigid deformations that affect the geometry of the local regions. When restricting the analysis to surfaces that only undergo isometric deformations, bending,  which is a subclass of non-rigid deformations, one  can design methods or descriptors that only capture the intrinsic properties of shapes, the properties that remain unchanged  after changes due to bending are factored out. Examples of such approaches include Generalized Multidimensional Scaling (GMDS) <cit.>, spectral analysis based on Laplace-Beltrami eigenfunctions <cit.>, and Heat Kernel Signatures (HKS) <cit.>.  Several  papers have attempted to relax the isometry constraint so that 3D objects can be put in correspondence even when they undergo stretching. Zhang et al. <cit.>, for example, put two objects  in correspondence by finding the optimal bending and stretching that is needed to align one shape onto another. In this work, bending is measured in terms of changes in the direction of the normal vectors while stretching is measured in terms of changes in edge lengths. Thus, the approach uses a metric that is similar to those presented in Section <ref>.






  §.§.§ Elastic registration as a reparameterization problem



Instead of using landmarks, some previous works tried to find the optimal diffeomorphism ^* that puts f_2 in full correspondence with f_1. Since the search space for optimal re-parameterizations is large,
some  of the previous works  such as  SPHARM <cit.> and  SPHARM-PDM <cit.>, used a fixed parameterization that is analogous to the arc-length parameterization on curves.  Others  restricted the type of re-parameterizations to a specific subset of the space of diffeomorphisms.  Lipman et al. <cit.>, for example,  showed that  the  optimal correspondence between nearly-isometric genus-0 surfaces, conformally embedded onto a sphere,  can be found by searching for the best Möbius transform that aligns them. The approach is, however, restricted to nearly isometric surfaces. In addition, it only finds a sparse set of correspondences since complex surfaces cannot be conformally embedded onto a sphere with minimal distortion. Zeng et al. <cit.> extended Lipman et al.'s  approach and finds a dense matching in the presence of large elastic deformations.  Kim et al. <cit.>  showed that by blending multiple maps, each one is a Möbius transform, one can match surfaces under relatively large elastic deformations. Nevertheless, the correspondences with these methods are not one-to-one. The inter-surface  maps they produce are not guaranteed to be a diffeomorphism. Thus they are not suitable for statistical shape analysis tasks.  In addition, these methods are limited to surfaces that can be conformally parameterized onto a sphere. Although, in theory, any closed genus-0 surface can be conformally embedded into a sphere, in practice, however, the distortion that results from conformal embedding of complex surfaces that have elongated extruding parts, human or animal body shapes, is often very large. As a consequence, several important shape parts are being under-sampled  and thus cannot be correctly matched across surfaces.




Using a fixed set of landmarks  or a fixed parameterization of the surfaces being analyzed imposes constraints on the correspondences that can be found and often may lead to undesirable results especially in the presence of large elastic deformations. 
Fig. <ref>-(a) shows an example where  fixed parameterization fails. Better results can be  achieved by   treating the boundaries of objects as continuous surfaces, rather than discretizing them into point sets at the outset, and solving directly Eqn. (<ref>) for finding the optimal  sampling (reparameterization ) that better matches features across shapes as proposed in <cit.>, see Fig. <ref>-(b) for a 2D example. Continuous  representations make it possible to elastically model and optimize correspondences.  The process of matching two 3D shapes  is, in fact, reduced to the problem of optimally re-parameterizing one shape so that it is as close as possible to the other.  Thus, correspondence becomes a re-parameterization problem, which can be mathematically formulated using diffeomorphisms. 

In order to do that, however, the action of  on the space of surfaces should be by isometries, under the metric d_(·, ·)  
    ∀∈, d_(f_1∘γ, f_2 ∘γ) = d_(f_1, f_2).

Kurtek et al. <cit.> showed that when using the  metric in the space of surfaces , the condition of Eqn. (<ref>) is only valid when  is an area-preserving diffeomorphism. This is, however, too restrictive especially when dealing with elastic registration because elastic deformations correspond to changes in surface area.


Kurtek et al. <cit.> and Jermyn et al. <cit.> showed that in the space of Square Root Maps (SRM) and in the space of Square-Root Normal Fields (SRNF), respectively, the action of Γ is by isometries under the  metric. Given that the  metric in the space of SRNFs is equivalent to the partial elastic metric in the space of , one can 

	
  * Option 1: perform the analysis in  using the pullback of the partial elastic metric since the pullback metric inherits all the properties of the metric, or  
	
  * Option 2: map the surfaces being analyzed to the space of SRNFs, perform the analysis there, registration and geodesic computation, and finally map, if needed, the results back to the space of surfaces .


In practice, option 2 is computationally very efficient; We are given two surfaces f_1 and f_2, and their SRNF representations q_1= Q(f_1) and q_2=Q(f_2). We  perform  elastic registration by solving for the optimal rotation O and optimal re-parameterization γ such that  the  distance between q_1 and O(q_2, γ) is minimized.  
If  O^* and γ^* is the solution to this optimization problem, then  f̅_2= O^*(f_2∘γ^*) and its corresponding SRNF is given by q̅_2 = O^*(q_2, γ^*). 
Since this process uses the  metric, geodesics correspond to straight lines in the space of SRNFs and thus one can use all the standard computational tools available when dealing with Euclidean spaces.  

















 §.§ Geodesics




For clarity purposes, let's assume for now that we are given two surfaces f_1 and f_2 ∈ with a fixed parameterization. In other terms, a one-to-one correspondence between the surfaces is given and we are assuming  that this correspondence is correct.    The goal is to solve the inner optimization problem of Eqn. (<ref>), find a path ^*: [0, 1] → that minimizes Eqn. (<ref>). Solutions to this optimization problem depend on the choice of the metric and  the representation (or nature of the pre-shape space).  When using the  metric on , the solution is straightforward and in fact it can be expressed analytically as F^*(τ) = f_1 + τ (f_2 - f_1), with τ∈ [0, 1].
As shown in Figure <ref>-(a), using the  metric in  leads to undesirable distortions.  In this section, we will present the  main  strategies for computing geodesic paths when using elastic metrics that are physically motivated. We will also discuss details that depend on the choice of the metric and the representation.






  §.§.§ Geodesics using pullback metrics
  



The first class of approaches find geodesics by solving the optimization problem of Eqn. (<ref>) using one of the elastic metrics defined in Section <ref>. In the simplest case, we are given two surfaces f_1 and f_2 and we seek to find another surface f^* on the geodesic path F and which is exactly at the mid-point between f_1 anf f_2, f^* = F(0.5). Mathematically, f^* is the minimizer of the discrete path energy:

    f^* = min_f L(F) = (f_1, f) + (f, f_2).

which corresponds to solving for f the Euler-Lagrange equation

    ∇(f_1, f, f_2) = ∂(f_1, f)  /∂ f  + ∂(f, f_2)  /∂ f  = 0.

Heeren et al. <cit.> solve this non-linear optimization problem using the   Quasi-Newton method. That is, one can start with an initial estimate f^0  of f, f^0 = 0.5 (f^1 + f^2), and then iteratively update the estimate via a Newton iteration with step size control:
f^n+1 = f^n + ϵΔ f 
 where

    Δ f = -[  (f_1, f^n, f_1)   ]^-1∇.

The optimal step size ϵ can be determined using the David-Fletcher-Powell (DFP) formula <cit.>. Using this formulation, one can  solve for the entire time-discrete geodesic path recursively; once the mid-point f between f_1 and f_2 is computed, one can apply the same procedure to solve for the mid-point between f_1 and f, and between f and f_2, etc. <cit.>. Alternatively, one can directly find the entire discrete path F = { f_1= f^1, f^2, f^3, ⋯, f^K= f_2}, by optimizing

    (f^i, i= 1, ⋯, K) = ∑_i=1^K-1(f^i, f^i+1),

over  f^i, i= 2, ⋯, K-1. This is called a path straightening procedure. The optimization can start with an initial path, which is then iteratively straightened, updated so that its length with respect to the elastic metric is minimized.    In the space of shells <cit.>  this is done by   solving the set of Euler-Lagrange equations:

    ∂/∂ f^i = 0, i=2, ⋯, K-1.

Instead of solving a non-linear optimization problem, Kurtek et al. <cit.>  define the energy of a continuous path F as:

    E_path(F) = ∫_0^1  d F/dτ(τ)   d F/dτ(τ)     dτ.
 
The shortest path is then found by minimizing the path energy using a gradient descent approach where the current estimate of the geodesic path is updated according to F = F - ϵ∇ E_path. Let 𝒜_0 be the set of all smooth paths that start at f_1 and end at f_2 and let T_F(𝒜_0) be the tangent space to   𝒜_0 at a given point F.  T_F(𝒜_0) corresponds to perturbations of the path F. 
Kurtek at al. <cit.> describes the space of path perturbations using orthonormal basis ℬ = {α_1, ⋯, α_B } so that the gradient of the path energy can be written as

    ∇ E_path = ∑_i ∇ E_path (δα_i) δα_i.
 
We refer the reader to <cit.> for the full derivation of the directional derivatives ∇ E_path (δα_i) as well as the procedure for constructing the orthonormal basis ℬ.



 































  §.§.§ Geodesics in the space of SRNFs


 
The SRNF representation, compared to others,  is very promising as a representation of surfaces for elastic shape analysis. If geodesics, mean shapes, PCA, etc. could be computed in  and  under the  metric, and then mapped back to ,  there would be large gains in computational efficiency with respect to the optimization-based methods described in the previous sections and those introduced in  <cit.>.  Unfortunately,  unlike the curve case, there is no analytical expression for Q^-1, the inverse of SRNF map Q, for arbitrary points in . Moreover, the injectivity and surjectivity of Q remain to be determined, meaning that for a given q∈, there may be no f∈ such that Q(f) = q, and if such an f does exist, it may not be unique. 
If one cannot invert the representation, one can transfer geodesics and statistical analyses conducted in  back to  by pulling the  metric back to   under and perform computations there, as in <cit.>, but this is computationally expensive, and rather defeats the purpose of having an  metric in the first place.



To solve this inversion problem, Xie et al. <cit.> showed that in the case of star-shaped genus surfaces, an analytical solution to the inversion problem exists. Laga et al. <cit.>, on the other hand,    developed a method that, given q∈, finds an f∈ such that Q(f) = q, if one exists, or an f whose image Q(f) is the closest, in terms of the elastic metric, to q if it does not. 
This is achieved by formulating SRNF inversion as an optimization problem: find an element f ∈ whose image Q(f) is as close as possible to the given q∈ under the  norm.  This optimization problem is efficiently solved using a gradient descent approach <cit.>. To avoid undesirable solutions, since gradient descent procedures converge to local minima if not initialized appropriately, Laga et al. <cit.> 
 carefully engineered an orthonormal basis of , and used a spherical-wavelet based multiresolution and multiscale representation of the elements of . The basis can be generic, spherical harmonics, or domain-specific such as PCA basis in the presence of training datasets.  With this approach,  SRNF maps can be inverted with high accuracy and efficiency, with robustness to local minima, and with low computational complexity.  

Finally, the shapes on the geodesic between two star-shaped surfaces are not necessarily star-shaped, and thus,  the analytic form may not be the correct inversion. However, one can use the analytical solution   as an initial guess for the inverse, thereby better initializing the reconstruction-by-optimization problem.

By mapping back to  geodesics computed in , which are straight lines since they are equipped with the  metric, the  computational cost is reduced by  an order of magnitude compared to methods that use path straightening under the pullback metric.  
See Section <ref>  for a few examples of geodesics between pairs of complex genus-0 surfaces.


































  §.§.§ Comparison and discussion

As the reader can see, many shape spaces and metrics have been defined in the literature. To help the reader navigate in this rich literature, we summarize in this section their properties and discuss their similarities and differences.

The first  aspect to consider when choosing a metric and a shape space is  the invariance properties of the metric, which defines the type of pre-processing requirements. For  instance, as summarized in Table <ref>,  the metrics of Kilian et al. <cit.> as well as those defined in the space of thin-shells <cit.>  are invariant to translation, scale, and rotation and thus one does not need to normalize the input 3D models in a pre-processing step. On the other hand, diffeomorphisms do not act by isometries in these shape spaces under the defined metrics. As a consequence, these metrics are not invariant to re-parameterization and thus are not suitable for joint registration and geodesic computation.  On the other hand, the  metric on the space of square-root maps such as the SRNF <cit.> and Q-maps <cit.> are invariant to re-parameterization. Thus the correspondence and registration problem can be solved jointly, using the same metric, with the geodesic computation. The SRNF requires normalization for scale while the Q-maps require normalization for translation and scale. Finally, methods based on M-reps <cit.> and Lie bodies <cit.> require normalization for all rigid transformations as well as the pre-computation of correspondences using a different approach.  



The other aspect to consider is the types of deformations that are quantified, and subsequently,  the complexity of the metric and the representation, which affect the processing time. Ideally, the process of evaluating the metric should be fast so that geometry processing and shape analysis tasks can be performed in a computationally reasonable time. 
From this aspect, Euclidean shape spaces and Euclidean metrics are probably the most convenient to use, since evaluating the  metric, computing short paths, and performing statistics under the metric are straightforward. Moreover, when working in Euclidean spaces one can  use all the computational tools developed for vector spaces, interpolations and extrapolations, Principal Component Analysis (PCA), and probabilistic models. As such, the SRNF representation <cit.> is very convenenint since all these tasks can be performed in the SRNF space using the  metric, which is equivalent to the partial elastic metric in , and map the results back to  for visualization <cit.>. 







§ STATISTICAL ANALYSIS UNDER ELASTIC METRICS



At this stage we have all the necessary computational ingredients, which are a shape  space, a metric,  and a mechanism for computing  correspondences and geodesics. In this section, we will show how to use these ingredients for performing full statistical  analysis of the shape of 3D objects, given an input collection of 3D models, we want to

	
  * compute the mean shape and the modes of variations. This is often referred to as the the shape atlas,  
	
  * characterize a population of 3D shapes with a probability distribution, and finally 
	
  * synthesize arbitrary shapes by sampling from the fitted probability distribution. 

We have already seen in Section <ref> that when the metric is Euclidean, in the case of morphable models for example, computing such quantities is straightforward using standard tools in vector spaces. In this section, we focus on the case where the shape space and the metric are non-linear. 






 §.§ Statistical analysis using non-Euclidean metrics



When using non-Euclidean metrics,  Eqn. (<ref>), which is equivalent to Eqn. (<ref>), cannot be solved analytically. Several previous works have proposed 
gradient-based approaches  for finding the Karcher mean of Eqn. (<ref>), see for example <cit.>. For convenience, the general approach is  repeated here as  Algorithm <ref>.




In addition to computing the mean shape of a collection of 3D models, we are usually interested in understanding (and modelling) how the shape varies within the collection.  Let {f_1, ⋯, f_n} be a collection of shapes and  be their mean. For simplicity of presentation, we assume that all the surfaces have been registered to . Let α_i be the geodesic path, with respect to the chosen metric, between  and f_i. Fletcher et al. <cit.> defined the following concepts that are essential for representing variability of data that lie on general manifolds:

	
  * Variance. Flechet <cit.> defined variance as  the expected value of the squared Riemannian distance from the mean.
	
    σ^2 = 1/n∑_i=1^n d_ ([], [f_i])^2.
	
	
  * Geodesic subspaces. The lower-dimensional subspaces in PCA are linear subspaces. For general manifolds, it is natural to use geodesic curves as one-dimensional subspace, which is equivalent to principal directions in PCA.  Flechet <cit.> showed that  submanifold geodesics at the mean  are the generalization of the linear subspaces of PCA.	
	
  * Projection. The projection of a point f ∈ onto a geodesic submanifold H is defined as the point on H that is the nearest to f in terms of  the Riemannian distance. Thus, the projection
operator π_H(f): → H  can be defined as:
	
    π_H(f) = _[y] ∈ H d_([f], [y] )^2.


Given the non-linearity of the shape space, the modes of variation can be efficiently computed by first projecting the data to the tangent space to  at the Karcher mean using the exp map, performing standard PCA in that linear space, and then mapping the results back to  using the Log map. Algorithm <ref> summarizes this procedure. 
 

Note that these two algorithms for computing the Karcher mean and the modes of variation are generic and can be used with any representation of surfaces and metrics. They both involve co-registration, geodesics computation, and the computation of the exponential and log maps.  The former can be omitted when the surfaces are pre-registered using a different procedure and different optimality criteria. This is however suboptimal.  The computation of geodesics, exponential maps and log maps depend on the metric and the shape space.  When using M-reps or Lie bodies of Section <ref>,  geodesics and geodesic distances have analytical solutions.  In other cases, when using the metrics in the space of shells or when using the square-root representations, one has to use the optimization procedures described in Section <ref>.










 §.§ Statistical analysis by SRNF inversion




Computing shape summaries on non-linear shape spaces using complex elastic metrics, or deformation energies, requires computing geodesics many times. For instance, at each iteration of  Algorithm <ref>, n geodesics should be computed where  n is the number of shapes in the collection to average. This is computationally very expensive as discussed in <cit.>. SRNF representations offer an elegant and a computationally efficient alternative solution; 
Instead of using the pullback metric, one can simplify the iterative procedure given in Algorithm <ref> by first mapping all surfaces to the SRNF space, , resulting in {_1,…,_n}. Then, the mean shape in SRNF space, denoted by , is computed by iterating between co-registering all _is to , and subsequently updating . This is computationally very efficient since  is an Euclidean space. Finally, at the last step, the mean shape ∈[] is computed by SRNF inversion, finding a shape  such that  (f̅) =. The entire procedure is summarized in Algorithm <ref> and detailed in <cit.>.




Next, let  denote the SRNF of the final average shape. With a slight abuse of notation, let _i, i=1,…,n denote the SRNFs of the surfaces in the sample optimally registered to this average. Then,  since the SRNF space is Euclidean, the modes of variation can be computed in a standard way, using singular value decomposition (SVD) of the covariance matrix. Let u_q^k denote the k-th principal direction. Then, one can explore the variability in this direction around the mean using + λ u_q^k, where λ∈. To visualize the principal directions of variation in , we only need to compute, by SRNF inversion, f^k_λ such that (f^k_λ) = + λ u_q^k.  


Finally, given a sample of observed surfaces f_1,…,f_n ∈, one wants to fit a probability model to the data and sample from it new 3D shape instances. Estimating probability models on nonlinear spaces like  is difficult: classical statistical approaches, developed for vector spaces, do not apply directly. By first computing the SRNF representations of the surfaces, one can use all the usual statistical tools in the vector space , and then map the results back to  using the  SRNF inversion algorithm.  Let G(q) be the model probability distribution fitted to {q_1,…,q_n}.  A random sample q_s can be generated from G. We then find f_s such that Q(f_s) = q_s.  Note that G can be an arbitrary distribution, parametric or non-parametric. Here, we use a wrapped truncated normal distribution, which can be learned using Principal Component Analysis on .   We caution the reader that the distribution is defined and analyzed in SRNF space; a distribution on  induces a distribution on , but it  has not been derived explicitly. 


In terms of computational complexity, computing the mean shape only requires the inversion of one  SRNF, while   techniques based on the pullback metric such as  <cit.> are iterative, and require the computation of n geodesics per iteration, with n being the number of shapes to average. This is computationally very expensive especially since each geodesic is computed using the expensive pullback metric.






§ EXAMPLES AND APPLICATIONS




We consider in this section a few examples of elastic 3D shape analysis tasks using the tools discussed in this chapter. Note that this section is not intended to provide a thorough performance evaluation and comparison of the methods described in this chapter but instead it aims to show their potential and the range of shape analysis problems that can be efficiently solved with these methods. 

	
	


 §.§ Registration and geodesic deformations

The first example we consider is the joint  registration and geodesic deformation of 3D objects that undergo complex elastic deformations. Here, we would like to emphasize for the reader that elastic (or non-rigid) deformation is composed of bending and stretching. Works based on (generalized) multi-dimensional scaling <cit.> and (heat-)diffusion geometry <cit.> only consider bending, which preserves intrinsic properties of surfaces. Here, we consider full elastic deformations. 

In Figure <ref>, we consider  two human body shapes of the same subject but in significantly different poses and compute a joint registration and geodesic deformation. We first compare the quality of the registration when computed with a re-parameterization invariant elastic metric (Figure <ref>-(b) and Figure <ref>-(c)) and when computed with other methods, functional maps of <cit.>. As one can see in Figure <ref>-(a), the later method does not produce correct one-to-one correspondences compared to the methods that formulate registration as a re-parameterization problem. Next, we compare geodesics computed using the  metric in  vs. geodesics computed using elastic metrics by SRNF inversion. This is shown in Figures <ref>-(b), (c) and Figure <ref>-(a), (b), where shapes undergo complex bending and stretching.   
 








 §.§ Elastic co-registration of 3D shapes


Standard registration methods are often pairwise, they put in correspondence a pair of 3D objects. Here, in statistical shape analysis, we are dealing with multiple objects that should be registered to each other in a globally optimal way. Moreover, as discussed in Sections <ref> and <ref>, using a different criteria for registration, geodesics, and mean computation leads to suboptimal solutions.  An important by-product of the Karcher mean computation using a re-parameterization invariant metric is the simultaneous co-registration of multiple objects.  This is illustrated with the examples of Figures <ref> and <ref>. In the first example, we take multiple 3D human body models in different poses and belonging to different subjects, put them in correspondence while simultaneously computing their Karcher mean as described in Section <ref>.  We perform the same experiment on the example of Figure <ref> this time taking 3D models of the hamate bone in the hand wrest, compute the mean shape while simultaneously registering all the 3D models in the corpus.
















 §.§ Classification

The length of a geodesic, or deformation energy, between a pair of surfaces in the shape space is a measure of
dissimilarity that can be used for unsupervised as well as supervised classification of 3D objects. To demonstrate this, we use the SHREC07 watertight 3D model benchmark <cit.>, which is composed of 400 watertight 3D models evenly divided into 20 shape classes. We only consider the 13 classes that are composed of genus-0, triangulated meshes. First, we compute their spherical parameterizations <cit.>, and then normalize them  for translation and scale. Thus, these surfaces become elements of the preshape space . Next, we map the surfaces to the space of SRNFs,  and compute their pairwise distance matrix. That is, for a given pair of surfaces f_1 and f_2 and their respective SRNF maps q_1 and q_2, their dissimilarity is given as 

    d((f_1, f_2) = d_([f_1], [f_2]) =  min_O,  q_1 - O(q_2 ∘) ^2.

For comparison, we also compute the pairwise distances between the shapes using (1) the  distance in , i.e. without registration, and (2) the  metric in  after elastic registration using the SRNF framework. Figure <ref> shows the multidimensional scaling (MDS) plot of these three pairwise distance matrices. As one would expect,  when using the  metric in , the 3D models are all spread across the space (Figure <ref>-(a)). Clusters start to emerge when using the  metric in  with elastic registration performed using the SRNF framework (Figure <ref>-(b)). When using the  metric in , which is equivalent to partial elastic metric (Figure <ref>-(c)), we can clearly see that the 3D models are clustered by shape classes. This suggests that the SRNF framework is suitable for 3D shape classification and clustering. We refer the reader to <cit.> for a  detailed  analysis of the performance of  SRNFs in  classification and clustering. 









 §.§ Random 3D model synthesis


Finally, we consider the problem of synthesizing arbitrary 3D models that are similar to, but note exactly the same as, a given collection of 3D models.  We  consider a collection of 398 human shapes <cit.> composed of multiple subjects in different poses. The first subject has 35 different poses including a neutral one. All other subjects have a neutral pose and a few other poses.  We first compute the SRNF representations of these surfaces, perform statistical analysis in the SRNF shape space , using standard linear statistics such  as Principal Component Analysis, and finally map the results to the surface space, F, using the proposed SRNF inversion algorithm. Formally, if q̅ is the mean and V_i, i=1, ⋯, k are the leading eigenvectors in the space of SRNFs, then any arbitrary SRNF map q can be written as:

    q = q̅ + ∑_i=1^k λ_i V_i,

where λ_i are real coefficients. An arbitrary shape can then be synthesized by finding a surface f such that its SRNF is q. This can be efficiently computed using the SRNF inversion procedure described in Section <ref> and detailed in <cit.>.


Figure <ref>-(a) shows the mean shape and the first five modes of variation. Figure <ref>, on the other hand, shows 10 randomly syntehsized human body shapes. Observe that the statistical model can generate shapes of arbitrary subjects in arbitrary poses.











 §.§ Other applications

Elastic registration,  geodesics and statistical summaries are building blocs to many applications in computer graphics, computer vision, and medical image analysis. The methods presented in this chapter are able to perform most of  these tasks.   As a consequence, the framework has been used in various applications.  In medical applications,  Kurtek et al. <cit.> used the Euclidean distance in the space of Q-maps for shape-based analysis of 3D subcortical structures in the brain in order to detect attention-deficit hyperactivity disorder (ADHD).  Joshi et al. <cit.>, on the other hand,  used the SRNF representation for shape-based hippocampal modeling in Alzheimer's disease.  
In particular, Joshi et al. showed that the elastic metric preserves important anatomical features when matching hippocampal structures and thus is potentially useful for accurately capturing the biological homology. The approach naturally provides a statistical framework for computing means and covariances in the shape space, which can be used for modeling the biological variability of subcortical structures. They observed a slight overlap of the statistical group difference results, obtained using SRNF, when comparing with the SPHARM approach <cit.>. Joshi et al.'s results, however, reveal additional effects due to disease or progression in different subregions of the hippocampal anatomy. One possible explanation is that the deformations generated using re-parameterization-invariant elastic metrics also account for the variability due to parameterizations, and thus may contain a richer description about the shape variability of the anatomical structures.  Samir et al. <cit.> used the SRNF framework for elastic shape analysis of cylindrical surfaces and applied it to  3D/2D registration in endometrial tissue characterization and also to simulate the deformation of such tissues <cit.>. 

Elastic shape analysis, in particular SRM-based elastic shape analysis <cit.>  has been later extended to include landmark constraints in the analysis process <cit.>. The landmarks can be used either to guide  registration in the presence of large and complex elastic deformations <cit.>, or as anatomical features that one needs to preserve when  generating patient-specific atlas of anatomical organs such as carpal bones <cit.>.  








§ SUMMARY AND PERSPECTIVES



We have reviewed  in this chapter the fundamentals as well as the state-of-the-art methods for non-rigid 3D shape analysis. We have seen that there are three  fundamentally inter-related problems, which are  correspondence (registration), geodesic deformation, and statistical summarization of shapes that undergo elastic deformations. Thus, these problems are better solved jointly using the same metric that captures (and quantifies) physical deformations, which are bending and stretching.  We have particularly seen   three types of metrics, mainly (1) the Euclidean distance or  metric in the space of shapes , (2) the pullback of elastic metrics, including the elastic metric in the space of shells, and (3) the  metric in the space of Square-Root Normal Fields (SRNF). While metrics in the first  class  are easy to use and have lead to some of the major developments in statistical shape analysis including the active shape models  as well as the morphable models, they are limited to 3D objects that do not undergo large bending and stretching. Moreover, they are not reparameterization invariant and thus the correspondence problem should be solved in a pre-processing step using a different approach and a different optimality criteria.  


The second class of methods capture the true deformations of shapes and thus can generate natural geodesics and shape summaries. They are, however, computationally very expensive; computing one geodesic may take hours in some cases since they require solving a complex optimization problem. In fact, only a few methods (M-reps and Lie bodies) have closed-form formulas for computing geodesics and geodesic distances. Other methods, compute them by solving expensive optimization problems. 


Square Root representations, on the other hand, are particularly promising since they map 3D shapes onto an Euclidean space  where the  metric  is equivalent to the partial elastic metric in the preshape and shape spaces  and . Thus geodesics   become straight lines.  One can then perform all the computations in the space of SRNFs using standard vector calculus and only map the results back to the space of surfaces for visualization purposes.  Moreover, diffeomorphisms act by isometries on the elements of the SRNF space under the  metric. Thus, one can solve jointly the elastic registration problem under the same metric. The main limitations of these methods is that the  metric on the space of SRNFs is only equivalent to the partial elastic metric with fixed weights for the bending and stretching terms. The existence of such nice representation for arbitrary weights is in fact unknown and can be a promising direction for future research. 






















Finally, despite the rich literature in this field, there are several fundamental problems and challenges that require further research. Below we summarize some of the open challenges that we believe require further investigation. 


Topological and structural variabilities.   
While a significant progress has been made in studying geometric variability (driven by the computer vision and medical imaging communities),    analyzing structural variability of man-made 3D shapes, which is receiving a growing interest by the graphics community,  is still at its infant stage.  Structural variability  is highly related to the semantics and functionality of  shapes <cit.>. Methods that studied it  (driven by the computer graphics community and  motivated by the need to easily  create 3D models and enrich collections),  do not come with a concept of structure space and a metric for performing true statistics on structure. On the other hand, discrete and  continuous representations (driven by the computer vision, medical and biological imaging communities) lead to well defined shape spaces and metrics.  These methods, however,  perform poorly on 3D shapes that are composed of multiply interacting components, subcortical structures in the brain.  




Multiply-interacting shapes.  Objects are not isolated in nature. Their shape is often influenced by their surroundings and by objects with which they interact. Examples include   anatomical surfaces, subcortical structures in human brains,  or man-made objects whose shape is influenced by their surroundings and by the actions to which they are subject to. Graphics papers that explored structural variability in man-made 3D shapes capture some of these interactions. Extending these techniques  to natural objects (subcortical structures) is an interesting direction of research to explore. We believe that developing  a  mathematical representation of the space of structures  and capturing and modeling the way shapes of multiple objects co-vary  are  open challenges that  require  further investigation. Solving them  will open perspectives for new applications in medical image analysis, biology and 3D modeling in computer graphics.









Correspondence.   Despite the large amount of research on this particular topic, correspondence between 3D shapes still remains an open challenging problem, particularly in the presence of large elastic, topological and structural variations. Methods where correspondence and analysis are solved jointly under the same metric provide promising results, but still   fail when analyzing 3D shapes that are semantically or functionally similar but significantly differ in their geometry or structure.  This problem will be  simplified if semantics and functionality of shapes and their parts can be  incorporated into the analysis. Recent works that relate structure to functionality <cit.> provide a step forward in this direction. However, the main challenge when dealing with man-made 3D objects comes from the fact that correspondences are not one-to-one and thus cannot be represented with diffeomorphisms or other nice mappings.  Recent works on graph and tree-graph statistics <cit.> may provide part of the solution, although having a joint model that captures both geometric and structural variabilities remains an open problem.










































§ ACKNOWLEDGMENTS



§ ACKNOWLEDGMENT

The 3D human shape models used in this chapter have been kindly provided by Nils Hasler.  The other models are from the TOSCA <cit.> and SHREC07 datasets <cit.>. 






















IEEEtran































