Highly-Expressive Spaces of Well-Behaved Transformations: Keeping It Simple
Oren Freifeld
MIT CSAIL

Søren Hauberg
DTU Compute

Kayhan Batmanghelich
MIT CSAILJohn W. Fisher III
MIT CSAIL

freifeld@csail.mit.edu

sohau@dtu.dk

kayhan@csail.mit.edu

fisher@csail.mit.edu

Abstractv θ (x)nn

We propose novel finite-dimensional spaces of R → R
transformations, n ∈ {1, 2, 3}, derived from (continuouslydefined) parametric stationary velocity fields. Particularly, we obtain these transformations, which are diffeomorphisms, by fast and highly-accurate integration of continuous piecewise-affine velocity fields; we also provide an exact solution for n = 1. The simple-yet-highly-expressive
proposed representation handles optional constraints (e.g.,
volume preservation) easily and supports convenient modeling choices and rapid likelihood evaluations (facilitating
tractable inference over latent transformations). Its applications include, but are not limited to: unconstrained
optimization over monotonic functions; modeling cumulative distribution functions or histograms; time warping; image registration; landmark-based warping; real-time diffeomorphic image editing. Our code is available athttps://github.com/freifeld/cpabDiffeo

6
5
4
3
2
1
0
−1
−2
−3
−4
10

8

6

φθ (x; t)t=0
t = 0.5
t = 1.0
t = 1.5
t = 8.0

4

2

0
0

2

4

6

8

10

x

(a) Ω = [0, 10] ⊂ R(b) Ω ⊂ R2 (c) Ω ⊂ R2

Figure 1: An Ω → Ω map, φθ (·, t) : x 7→ φθ (x, t), derived
from a CPA velocity field, v θ : Ω → Rn . (a) A 1D example. Note v θ vanishes at fixed points of φθ (·, t). (b) Top:A 2D v θ shown in select locations. Middle: Rather thana Middlebury-style visualization, we use one emphasizing
the CPA property. The horizontal component, vhθ (left), andthe vertical component, vvθ (right), are shown as heat maps
whose colors range from blue = −λ, via green = 0, to
red = λ where λ = maxx∈Ω max(|vhθ (x)|, |vvθ (x)|). Bottom: Isrc ◦ φθ (·, 1) (where Isrc appears in Fig. 2). (c) a 2D
example with boundary and volume-preserving constraints.

1. Introduction
Spaces of well-behaved transformations, particularly diffeomorphisms, play a key role in computer vision. Unfortunately, current state-of-the-art representations of highlyexpressive diffeomorphisms are overly complicated. Thus,
despite their potential power and mathematical beauty, theirapplicability is somewhat limited, especially when large
datasets are involved or when computing resources and time
are limited. Moreover, owing to their complexity, using
powerful inference tools in such spaces remains challenging. Motivated by not only the practicalities of probabilistic
modeling and statistical inference but also a desire to makediffeomorphisms broadly accessible, we propose a representation that combines simplicity, expressiveness, and efficiency. Particularly, we propose novel spaces of transfor-

mations that are based on (fast, highly-accurate) integration of Continuous Piecewise-Affine (CPA) velocity fields.
Existing spaces offer only subsets of the benefits of the
proposed spaces: 1) high expressiveness; 2) finite dimensionality; 3) ease of implementation; 4) modest mathematical preliminaries; 5) convenient modeling choices; 6) easeof handling optional constraints; 7) fast and highly-accuratecomputations that lead to fast likelihood evaluations. Thesein turn facilitate the use of inference tools – typified by thecase of analysis-by-synthesis methods – that are usually too
expensive for rich transformation spaces. Possible applications are numerous, as we demonstrate with: image editing;
optimization over monotonic functions; modeling Cumulative Distribution Functions (CDFs) and histograms withorder-preserving geometry; time warps; landmark-based
image warping and animation. Another related potential application is image registration.

Acknowledgement: O.F. and J.W.F. are partially supported by
U.S. Office of Naval Research MURI program, award N000141110688,
and VITALITE, which receives support from U.S. Army Research Office
MURI, award W911NF-11-1-0391. S.H. was funded by the Danish Council for Independent Research, Natural Sciences.

1

Isrc
(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)
th

Figure 2: Samples from the prior (§ 3). (a-h) The top 3 rows echo those in Fig. 1b and Fig. 1c. The 4 row shows adeformed grid overlaid on the image. The 5th row shows select trajectories. Note that the trajectories and transformations
are differentiable (hence continuous) but not piecewise affine. The tessellation in (e-h) is a refinement of the one in (a-d). (a)& (e): T θ ∈ M . (b) & (f): T θ ∈ M vp . (c) & (g): T θ ∈ M ∂ . (d) & (h): T θ ∈ M ∂,vp . See § 3 for details.
Henceforth, n ∈ {1, 2, 3} and Ω is either the whole of
Rn or an interval/rectangle/box (in which case, Ω ( Rn ).
CPA velocity fields. We base our method on spaces of
CPA velocity fields (Fig. 1). The term ‘piecewise’ is w.r.t.
P, a tessellation (§ 3) of Ω. Let VΩ,P be such a space.
While VΩ,P depends on Ω and P, we will usually notationally suppress these dependencies, and will just write V. One
appeal of these spaces is that they are finite-dimensional andlinear (but note that their elements, i.e. the velocity fields,
are usually nonlinear). Let d = dim(V). The spaces Rd
and V are identified with each other, where every θ ∈ Rd
is identified with exactly one element of V, denoted by v θ ,and vice versa. A finer P implies a higher d and richer velocity fields (Figs. 2 and 3). Appealingly, easily-imposed
optional constraints yield useful linear subspaces of V.
From CPA velocity fields to trajectories. Modulo a
detail (addressed in § 3.2) related to the case Ω ( Rn , any
continuous velocity field, whether Piecewise-Affine (PA) ornot, defines differentiable R → Ω trajectories. Particularly,
if x ∈ Ω then v θ ∈ V defines a trajectory, t 7→ φθ (x, t),such that φθ (x, 0) = x; see Fig. 2. This trajectory is thesolution to the integral equation
Z tθv θ (φθ (x, τ )) dτ , v θ ∈ V . (1)φ (x, t) = x +0Remark 1. Equation (1), whose solution, φθ (x, ·), appears
both inside and outside the integral, should not be confused with the piecewise-quadratic Ω → Rn map, y 7→Ryv θ (x) dx. The latter, a popular tool in computervision [40] and numerical analysis, is unrelated to our work.
0n×1

CPA-Based (CPAB) transformations. Modulo that detail, any continuous velocity field, whether PA or not, defines a transformation; i.e., a map whose input and output are viewed as points in Ω. Particularly, letting x varyand fixing t, x 7→ φθ (x, t) is an Ω → Ω transformation. Without loss of generality we set t = 1 and defineT θ (·) , φθ (·, 1). Since we integrate CPA velocity fields,
we coin our transformations CPA-Based (CPAB). We write
T θ = exp(v θ ) to indicate the relation between v θ and T θ ;the rational for this symbol is explained in our supplemental
material (henceforth referred to as Sup. Mat.). We let
M , exp(V) ,{exp(v θ ) : v θ ∈ V}(2)denote the space of CPAB transformations, again notationally suppressing the dependencies on Ω and P; i.e., formally
we should write MΩ,P . Importantly, M is nonlinear.
Remark 2. CPAB transformations are not CPA (except indegenerate cases): while T θ is continuous, it is not PA.
CPAB transformations are “nice”. E.g., they are diffeomorphisms, (T θ )−1 ∈ M , and the inversion is easy.
Remark 3. x 7→ φθ (x, t = 1) should not be confused with(the non-diffeomorphism, parametric optical-flow-like representation) x 7→ x + v θ (x), the latter being only a Taylor
approximation of the former. Thus, the way we utilize flexible parametrized vector fields is different from, e.g., [31].

(a) P1(b) P2(c) P3(d) P4(e) P5Figure 3: Tessellating a 2D region in several resolutions.
Integration of CPA velocity fields. Integral equations
usually lack analytic solutions. Since CPA velocity fields
are Lipschitz-continuous and almost-everywhere smooth,generic integration solvers are quite effective for them.
However, we can do even better. One of our contributions
is showing that integration of such fields is given in either
closed form (when n = 1) or almost closed form (n > 1).
Besides the obvious pluses (accuracy, computing time), thissolution makes it easier to interpret the resulting trajectories/transformation and is key to the theorems in § 3.
A specialized numerical solver. In practice, that solution has one shortcoming. It requires tedious bookkeeping(if n > 1) and invoking certain routines (easy if n = 1 but
cumbersome if n > 1). This is especially a hurdle in GPU
implementations. We thus propose a practical alternative, aspecialized solver for integrating CPA velocity fields, which
is faster and more accurate than non-specialized solvers.
Convenient modeling, tractable inference. Smoothness priors on M are easy to build and use. As is common with nonlinear spaces of nonlinear transformations, thenonlinearity of (θ, x) 7→ T θ (x) prohibits closed-form posterior or maximum likelihood calculations; however, since(θ, x) 7→ T θ (x) is evaluated fast, so is θ 7→ p(data|T θ ).
This facilitates the use of inference methods that rely on
multiple likelihood evaluations such as most Markov Chain
Monte Carlo (MCMC) methods. Lastly, CPAB transformations support coarse-to-fine analysis.

2. Related Work
CPA maps. In addition to fields such as numerical analysis and control theory, CPA maps are also used in computer vision [7, 40]. Our use of them differs from suchworks in that we integrate them as velocity fields to obtain transformations. The ODE equivalent to Eqn. (1) isdφθ (x,t)= v θ (φθ (x, t)). Unlike us, who solve it and exdt
ploit its link to transformations, Lin et al. [26] use it only asa regression model for fitting patterns in velocity data.
Pattern theory. Representing objects via transformations acting on them is a cornerstone in Grenander’s patterntheory. Our work is influenced by the impressive works
initiated by Grenander and continued in the geometryoriented subcommunities of computer vision and medical
imaging. Due to space limits, we mention only a few:
[15, 16, 8, 24, 42, 4, 1, 18, 17, 43, 30, 45, 9]. However, rather than focusing on complicated, possibly ∞dimensional spaces (whose both representations and asso-

P

Nc

Nv

D = 6Nc

d = 2Nv

P1
P2
P3
P4
P5

4
16
64
256
1024

5
13
41
145
1025

24
96
384
1536
6144

10
26
82
290
1090

′
Table 1: Values of Nc , Nv , D = dim(VΩ,P), and d =dim(VΩ,P ) for the P’s shown in Fig. 3. See § 3.

ciated integration are in practice often discretized and/or
approximated), or the geometry of transformation spaces,
here we take a different approach, emphasizing simplicity
and practicality, in order to provide the larger computervision community with a simple and powerful tool. Thus,
our text avoids the direct use of differential geometry, modulo an optional discussion in the Sup. Mat. The ability
to rapidly evaluate transformations accurately and to capture a wide range of deformations while keeping complexity low sets this work apart from other works on diffeomorphisms. Also, in applications involving N landmark pairs,
most methods can deal with only small values of N , as theircomplexity is O(N p ), p > 1 (e.g., some methods invert adense nN ×nN matrix). Thus, e.g., they cannot utilize tools
for dense-correspondence extraction (e.g., [25]), highlighting a disconnect with the larger computer-vision community. Ours is O(N ) and most computations are parallelized.
Finitely-many affine building blocks. The PA transformations from [7] should be confused with neither PA/CPA
velocity fields nor CPAB transformations; see Remark 2. In
fact, PA transformations are neither differentiable (or evencontinuous) nor invertible. Closer to ours is the work of Arsigny et al. [3] who, like us, use finitely-many affine building blocks to build flexible stationary velocity fields. Whilethey spatially average the blocks, we use a CPA constraint
on them. A field in our space is usually not in their space(so their integration is inapplicable to it) and vice versa, preventing a direct comparison of integrations. Unlike in our
integration method, they use an approximation throughout;i.e., they must approximate the entire trajectory. Their approximation differs from others in that it uses a weighted
sum of matrix exponentials. As the approximation holds
only near zero, they must divide the field by a large number (which is tied to #steps). To keep #steps small, theysmartly generalize the scaling-and-squaring method. Their
scheme is exact only for a single block; however, expressiveness requires more blocks, so the accuracy drops. It isthus unsurprising they focus on a small number of blocks,
and that to save time they use their method only in the laststage of inference. For similar dimensions of spaces of velocity fields, we can use larger (hence fewer) steps thanthem since most of our steps are exact. Unlike them, forn = dim(Ω) = 1 (regardless of the number of blocks)

1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
−0.2
−0.4−8 −6 −4 −2
0.8

0

2

4

6

8

0.6
0.4
0.2
0.0
−0.2
−0.4
−0.60.0

0.2

0.4

0.6

0.8

1.0

1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
−0.2
−0.4−8 −6 −4 −2 0 2
0.25
0.20
0.15
0.10
0.05
0.00
−0.05
−0.10
−0.150.0
0.2
0.4
0.6

4

6

8

1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
−0.2
−0.4−8 −6 −4 −2
1.2

0

2

4

6

8

1.0
0.8
0.6
0.4
0.2
0.0
0.8

1.0

−0.20.0

0.2

0.4

0.6

0.8

1.0

1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
−0.2
−0.4−8 −6 −4 −2 0 2
0.30
0.25
0.20
0.15
0.10
0.05
0.00
−0.05
−0.100.0
0.2
0.4
0.6

4

0.8

6

8

1.0

1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
−0.2
−0.4−8 −6 −4 −2 0 2
8
7
6
5
4
3
2
1
0
0.0
0.2
0.4
0.6

4

0.8

6

8

1.0

1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
−0.2
−0.4−8 −6 −4 −2 0 2
0.45
0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.05
0.000.0
0.2
0.4
0.6

velocitytransformed

1.0
0.8

0.00

0.6
−0.05
0.4
4

6

8

−0.10
0.0
0.15

0.2
0.2

0.4
0.6
0.8
velocity

1.0

0.10

1.0

−1
0
1
2
transformed−23transformed

0.8

0.6

0.6

0.4

0.4

0.2
−2

−1
0
1
derivative

2

3

0.00
1.0

0.2
5
10
15
transformed

20

0.00
1.0

0.8

0.8

0.00

0.6

0.8

0.6

0.6

−0.05

0.4

0.4

0.4

−0.200.0

0.6

histogram

1.0

0.8

1.0

5

10
15
histogram

20

10

20

0.4
0.2

−0.15

we solve exactly, while for n > 1 we solve exactly mostof the trajectory, and use a numerical solver in only partsof it. Beyond the comparison of integrations, which favors
ours, our representation has additional advantages such that
it is simpler, suits better for a GPU implementation, handles constraints and several modeling choices more easily,and scales better with an increase in the number of blocks.
Discrete representations and approximations. Allassonnire et al. [1] efficiently approximate diffeomorphisms
via non-differentiable maps. While we focus on a generalpurpose representation, Diffeomorphic demons [43] is a fastregistration tool based on discretely-defined fields, which as
its authors note, may be inconsistent with a diffeomorphic
framework and may not preserve orientation. They also
cannot impose volume preservation. Lastly, a computer representation of a sequence of discrete fields requires plenty ofmemory. Such issues can be obviated by adapting demons
to our compact and continuously-defined fields. More generally, approximations based on discretely-defined velocity
fields and/or discrete diffeomorphisms are widely used, e.g.,
in medical imaging [43]. Unlike these works, both our fields
and transformations are continuously-defined.
Statistics on manifolds and tangent spaces. Like others, e.g. [3, 42, 11, 29, 13, 20], we deal with the nonlinearityof the space of interest by utilizing the linear tangent spaceat the identity. Though not explored here, other types ofstatistics on manifolds [32, 12, 34, 38, 36, 21, 19, 41, 28,
44, 10, 23, 14] may also be applicable.
CDF/histogram modeling. Working with CDFs is preferred to working with densities since the latter might notexist and, moreover, the Lp (or a sphere-based) distance between two densities might be arbitrarily large even if theirprobability measures are essentially the same. One principled approach to representing CDFs uses p-Wassersteinspaces, usually p ∈ {1, 2}, the p = 1 case is related to EarthMover’s Distance [33]. A limitation of this approach is thatit requires bounded pth -moments and non-trivial computations. 1-Wasserstein methods also lack a full geometric

1.0

1.2

−0.10

Figure 4: Monotonic regression on synthetic data (top) by
inferring latent CPA fields (bottom), defined over 100 equallength cells. In the 4 leftmost columns we used a Gaussian
likelihood; in the 2 rightmost columns we used a robust(Geman-McClure) one. Columns 1,3,5 show ML solutions.
The rest show MAP solutions with a smoothness prior. Regardless of outliers/likelihood/prior, the function is increasing (this is obtained effortlessly from the representation).

derivative

0.8

0.05
0.8

0.0−3
1.0

0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.05
0.00−3

0.2

0.4

0.6

(a)

0.8

1.0

0.0−3

0.2

0.2
−2

−1

0

(b)

1

2

3

0.0−3

−2

−1

0

(c)

1

2

3

0.00

0.2
5

10

(d)

15

20

0.00

5

15

(e)

Figure 5: CDF/histogram representation. 1st row: Nc = 3.
2nd row: Nc = 100. F0 , not shown, is a CDF of a uniform
distribution on J = [−3, 3]. (a) v θ ∈ V ∂ ([0, 1], P). (b)dF . Note it is not piecewise constant.
F = T θ ◦ F0 . (c) dx
H0 , not shown, is a cumsum of a 20-bin uniform histogram.
(e) H = T θ ◦ H0 , a cumsum of a new histogram h (f) h.
space whose elements are always valid histograms/CDFs,
complicating synthesis of new points. Methods that map ahistogram/density to a sphere (e.g., [35]) have such geometry but suffer from two issues. First, they do not respect theordering of the bins or the real line. While in some applications of histograms bin ordering is immaterial, it is strange
to ignore the ordering of R. Second, large moves on thesphere lead to CDFs/histograms with negative values. The
problems above do not exist in our method.
Image warping. Also related to our work is the problem
of image warping [5, 2, 27, 6]. Unlike most methods, ours
is diffeomorphic and handles additional constraints easily.
Stationary Velocity Fields. Here we focus on Stationary Velocity Fields (SVF). While not all diffeomorphisms
can be defined via SVF, this popular approach (e.g., [3])to reduce computations and simplifying inference still captures a wide range of real-world deformations. Our CPA
fields reduce computations and simplify inference even further. While it is possible to use non-stationary CPA fields(i.e., θ = func(t)), we are unaware of a specialized integration for these, so one may have to resort to standard solvers.
Either way, there may still be a practical gain since the dimension of the inference problem will remain finite (i.e., dtimes the number of time steps) and our arguments related
to priors, constraints, and coarse-to-fine analysis still hold.

3. The Proposed Mathematical Representation
The lemmas/theorems below are proved in the Sup. Mat.

3.1. CPA Velocity Fields
A tessellation of Ω, denoted by P, is a set of Nc closed
Nc
subsets of Ω, {Uc }c=1
, also called cells, such that theirunion is Ω and the intersection of any two adjacent cells
is just their shared border. An Nc ′ -cell tessellation P ′ ,Nc ′ > Nc , is a refinement of P, a relation denoted P ≺ P ′ ,if each cell in P ′ is a subset of a cell in P. Unless statedotherwise, it is assumed that the cells are either intervals

cInput: P = {Cc }N
c=1 ; (A1,θ , . . . , ANc ,θ ); U = a sequenceof Npts points in Ω; t > 0; Nsteps , nsteps ∈ Z+
Output: φθ (U, t) , {φθ (x, t)}x∈U ⊂ Ω
1 ∆t ← t/Nsteps ; δt ← ∆t /nsteps
2 for c ∈ {1, .., Nc } do in parallelfi
TAc,θ ,∆t ← exp(∆t A
3c,θ )4 for x ∈ U do in parallel5x0 ← x6
for i ∈ {1, .., Nsteps } doxetemp = [ xTtemp 1 ]T ← TAγ(xi−1 ),θ ,∆t x
ei−1
7

8

9

10
11
12

if γ(xi−1 ) == γ(xtemp ) thenxi ← xtemp // analytic updateelse// inter-cell bdry crossedxi ← genericODEsolver( ini. con.=xi−1 , step
size=δt , nsteps , v(·) : x 7→ Aγ(x),θ xe)φθ (x, t) ← xi

Algorithm 1: Integrating v ∈ V. See Sup. Mat. fordetails on the generic solver and a discussion on the
hypothetical concern related to powers of TAc,θ ,∆t .
(if n = 1), triangles (n = 2), or tetrahedra (n = 3). Let
Pvert = {ξ : ξ is a vertex of some Uc ∈ P} and let Nv denote its cardinality; e.g., see Fig. 3 and Table 1. Rather thanregarding Pvert as a set, we regard it as an ordered Nv -tuple(the ordering is arbitrary, but assumed fixed).
Fix P and let x ∈ Ω. We define the membership function
γ : Ω → {1, . . . , Nc } , γ : x 7→ min {c : x ∈ Uc } ; (3)i.e., if x is not on an inter-cell border then γ(x) = c ⇐⇒x ∈ Uc , while a border point is assigned to the cell of thelowest index. A map, f : Ω → Rn , is called PA (w.r.t. P) ifNce whereare affine; i.e., f (x) = Aγ(x) x{f |Uc }c=1
e , [ x1 ] .
x

our proposed solver
a generic solver

5

25

45

65

85

105

0.06
0.33

0.08
0.43

0.08
0.41

0.09
0.43

0.11
0.47

0.12
0.48

Table 2: Timing comparison (1D). Results, in [sec], are averages of timings of integration of random CPA fields. The
tessellations, of [0, 1], consist of Nc intervals. Npts = 105 ,Nsteps = 100 and nsteps = 10 (see Algorithm 1).
n × (n + 1) matrices. If A ∈ Rn×(n+1) , then vec(A) ∈2
Rn +n is its row-by-row flattening to a column vector.
T
Likewise, vec(A) , [ (vec(A1 ))T ... (vec(ANc ))T ] ∈ RD .
′An inner product on V is defined by hvA1 , vA2 i =vec(A1 )T vec(A2 ). We now explain how to build an orthonormal basis for V, how its elements are parametrizedby θ, and how PA velocity fields are projected onto V.
Lemma 2. An element of V is any vA ∈ V ′ such thatvec(A) satisfies a linear system of constraints, denoted by
L vec(A) = 0; see Sup. Mat. for the construction of L.

θAc ∈ Rn×(n+1) ∀c ∈ {1, . . . , Nc }

Nc

(4)

A map, f : Ω → Rn , is called CPA if it is continuous and
PA. If f is CPA w.r.t. P it is CPA w.r.t. P ′ ≻ P. A vector
field v (on Ω) is an Ω → Rn map viewed as the mappingof points to vectors; i.e., if x ∈ Ω then we view v(x) as ann-dimensional “arrow”. The terms velocity field and vector
field will be used interchangeably. A vector-field space V
is a linear space whose elements are vector fields; dim(V)′may be finite or infinite. Let VΩ,Pand VΩ,P denote thespaces of PA and CPA velocity fields on Ω w.r.t. P. Note′
VΩ,P ⊂ VΩ,P. For short, we will also write V ′ and V, suppressing the dependency on Ω and P.
Lemma 1. V ′ and V are linear spaces, D , dim(V ′ ) =(n2 + n) × Nc , d , dim(V) = n × Nv , and d ≤ D withequality if and only if Nc = 1.
See, e.g., Table 1. A generic element of V ′ is denotedby vA where A , (A1 , . . . , ANc ) groups its Nc associatedLet B = [ B1 ... Bd ] ∈ RD×d denote a specific orthonormal basis of null(L): the one obtained via SVD of L. LetTθ = [ θ1 ... θd ] ∈ Rd . If vec(A) = Bθ then vA is CPA.
Particularly, vvec−1 (Bj ) is CPA for every j ∈ {1, . . . , d},dand {vvec−1 (Bj ) }j=1 is a basis for V. For a visualization ofthis basis, see Sup. Mat. Regardless of whether vA is CPA
or not, vvec−1 (BB T vec(A)) is CPA. When earlier we denoted a generic element of V by v θ we meant that θ standsdfor the coefficients w.r.t. the {vvec−1 (Bj ) }j=1 basis:e=v θ (x) = Aγ(x),θ x

Xd

j=1

θj vvec−1 (Bj ) (x)

(5)

ÄP
ä
d
where Aθ , (A1,θ , . . . , ANc ,θ ) , vec−1
j=1 θj Bj .
This basis is orthonormal w.r.t. the inner product h·, ·iB :V × V → R,(v θ1 , v θ2 ) 7→ θ1T θ2 .
Letθθvvert , {v (ξ)}ξ∈Pvert denote the ordered Nv -tuple of ndimensional values v θ takes at the Nv points in Pvert .
θLemma 3. Lθ7→vvert, the map sending θ to vvert, is a linearθ
bijection and its associated d × d matrix is given in closed
form. We denote its inverse by Lvvertθ7→θ .

If Ω ( Rn , let ∂Ω denote the boundary of Ω, and let
n : ∂Ω → S n−1 denote the unit normal to ∂Ω.
Lemma 4. Optional constraints such as v θ (x), n(x) =Nc
0 on ∂Ω and/or {tr(Ac ) = 0}c=1, are linear and can thus
extend L to have more rows. The null space of the extended
L is a linear subspace (whose dimension is denoted by somed′ < d) of the null space of the original L. Reusing the′symbol B = [ B1 ... Bd′ ] ∈ RD×d to denote the orthonormal basis of this new null space obtained via SVD of theextended L, we get a d′ -dimensional basis of CPA velocityd′fields, {vvec−1 (Bj ) }j=1 , that satisfy the new constraint(s).

P

P1

P2

P3

P4

P5

Npts

Nsteps

P1

P2

P3

P4

P5

our proposed solver
a generic solver

0.06
0.22

0.09
0.23

0.15
0.25

0.24
0.26

0.27
0.28

64 × 64
64 × 64
256 × 256
256 × 256512 × 512
512 × 512
1024 × 1024
1024 × 1024

10
100
10
100
10
100
10
100

.001
.002
.004
.023
.011
.079
.031
.255

.001
.003
.006
.028
.016
.129
.051
.414

.002
.003
.009
.046
.025
.185
.083
.519

.002
.005
.011
.064
.036
.287
.146
1.005

.006
.014
.014
.083
.041
.300
.143
1.070

Table 3: Timing comparison (2D). Results, in [sec], are averages of timings of integration of random CPA fields. The
tessellations are as in Fig. 3, the image size is 512 × 512,
Nsteps = 100 and nsteps = 10 (see Algorithm 1).
The corresponding subspaces are denoted V ∂ , V tr andV= V ∂ ∩ V tr . We still denote an element of any of′these by v θ , with the understanding that then θ ∈ Rd . We
also define, similarly to Eqn. (2), subsets of M : M ∂ ,exp(V ∂ ), M vp , exp(V tr ) and M ∂,vp , exp(V ∂,tr ).

Table 4: Timings (in [sec]) for varying Npts and Nsteps .
The setup is as the one described in in Table 3.

∂,tr

3.2. CPA-Based Transformations
Via Eqn. (1), v θ ∈ V implies T θ ∈ M , a relation denoted T θ = exp(v θ ). The detail alluded to in § 1 is that for
Eqn. (1) to be well defined, φθ (x, τ ) must always be in Ω,
the domain of v θ . This holds if Ω = Rn . It also holds ifΩ ( Rn and v θ ∈ V ∂ . If Ω ( Rn and v θ ∈ V \ V ∂ , it
might not. Also, some of the results below require both thecontinuity and PA properties of v θ . We need a lemma.
Lemma 5. If Ω ( Rn , we can extend P and constrain thevelocity fields such that we obtain a linear subspace of V
whose elements extend to CPA fields on Rn . In which case,
we also redefine Ω to be equal to (the whole of) Rn .
Henceforth, if Ω ( Rn and we work with V (or V tr ),as opposed to V ∂ (or V ∂,tr ), we will assume the procedure from Lemma 5 has been applied. Thus, T θ is always an Ω → Ω map, Eqn. (1) is well defined and v θ isindeed CPA. Note that ∀x ∈ Ω, φθ (x, 0) = x and that
{x : v θ (x) = 0} are fixed points of T θ (see Fig. 1a). Let
n×(n+1)
us define, let
î a usefuló building block. If A ∈ R
A(n+1)×(n+1)
eA = 01×(n+1) ∈ R. Let expm denote thee is [ 01×n 1 ].
matrix exponential. The last row of expm(A)ee isMoreover, det expm(A) > 0. Particularly, expm(A)einvertible. If n = 1 or A has a special structure, expm(A)has a closed form. Otherwise, since n is small and due toe a generic expm routine approximates it
the structure of A,twell. Let t ∈ R and define ψθ,c: Ω → Rn viaóî tψθ,c (x)fie , TAc,θ ,t , expm(tA, TAc,θ ,t x
(6)
c,θ )
1

(a solution to an ODE with an affine velocity field, ξ 7→
e We note that the solution to Eqn. (1) is the compoAc,θ ξ).
sition of a finite number, denoted by m, of such solutions:tmt2
t1
φθ (x, t) = (ψθ,c◦ . . . ◦ ψθ,c◦ ψθ,c)(x) .
m
2
1(7)

As mentioned earlier, T θ is defined via T θ (x) = φθ (x, 1).
The compact form of Eqn. (7) hides an ugly truth: the numm
ber of the trajectory segments, m, their durations, {ti }i=1 ,mand the indices of the cells involved, {ci }i=1 (where a cell
may appear more than once), all depend on x. Except thefirst cell, c1 = γ(x), they also depend on θ and t. Thus, amore precise and cumbersome notation would be:tm(ψθ,cx,θ,tm(x,θ,t)x,θ,tt2 (x,θ,t)t1 (x,θ,t)(x,θ,t) ◦ . . . ◦ ψθ,c2 (x,θ,t) ◦ ψθ,γ(x) )(x) .

Remark 4. While invertible affine matrices form a group,φθ (·, t) : Ω → Ω is not affine, exactly because of the xdependencies mentioned above. In fact, since the above
quantities vary with x even within a cell, the continuous
x 7→ φθ (x, t) is not PA, hence not CPA.
A map, T : Ω → Ω, is a called a diffeomorphism (onΩ) if T −1 exists and both T and T −1 are differentiable. Let
G be the space of (orientation-preserving) diffeomorphismson Ω. Let H ⊂ G be its restriction to diffeomorphismsthat can be obtained from differentiable stationary velocity
fields via integration; i.e., T (·) = φ(·, 1) where φ(x, t) =Rt
x+ 0 v(φ(x, τ )) dτ with a differentiable v : Ω → Ω. Both
G and H are ∞−dimensional nonlinear spaces.
Theorem 1. (i) If α ∈ R then v θ = v αθ /α and φθ (·, t) =φαθ (·, t/α). (ii) (T θ )−1 exists, is in M , and equals toT −θ . (iii) M ⊂ G. (iv) If T θ ∈ M vp then T θ is volumepreserving (hence the vp). (v) M is a d-dimensional nonlinear space; similar statements hold for M ∂ , M vp andM ∂,vp , with their respective values of d′ . (vi) If P1 ≺ P2 ≺. . . is a tessellation sequence such that eventually all them→∞cells become arbitrarily small, then MΩ,Pm −−−−→ H.
Corollary 1. If n = 1, T θ is increasing. Let J be either Ror an interval. If F0 : J → Ω is either non-decreasing, increasing, right-continuous, continuous, differentiable, a diffeomorphism, or a step function (7 non-mutually-exclusive
cases), then so is F θ , T θ ◦ F0 : J → Ω. If, in addition,Ω = [0, 1], v θ ∈ V ∂ and F0 is a CDF, then F θ is a CDF.
This gives an unconstrained parametrization of large
classes of increasing functions and CDFs/histograms via afinite-dimensional linear space (Figs. 4 and 5). In the CDF
case, F0 need have neither a density nor finite moments.
Equation (7) is very useful. It justifies and leads to a simple proof of Theorem 1, provides an insight into the struc-

0

0

src

100

0

dst

100

0

src

100

dst

200

200

200

200

300

300

300

300

400

400

400

400

500

500

500

600
00

100

200

300

400

500

600

700

src

100

600
00800

100

300

400

500

600

700

dst

100

200

200

600
00800

500

100

200

300

400

500

600

700

src

100

200

600
00800

dst

300

300

300

400

400

500

200

300

400

(a) Isrc

500

600

700

600
0800

500

100

200

300

400

(b) Idst

500

600

700

600
0800

300

400

500

600

700

800

400

500

600

700

800

200

400

100

200

T θ (src)
dst

200

300

500

100

100

400

600
0

T θ (src)

100

dst

500

100

200

300

400

500

600

700

600
0800

(c) v θ

100

200

300

(d) Isrc ◦ T θ(e) (Isrc ◦ φθ )(·; t) and Idst

Figure 6: Landmark-based warping. One example in each row. (a) Isrc (and src landmarks). (b) Isrc (and src/dst landmarks).
(c) The inferred v θ . (d) A new image is synthesized by warping Isrc using T θ (and dst & T θ (src) landmarks). (d) Animation.
ture of t 7→ φθ (x, t), and, when n = 1, leads to a closedform solution for x 7→ T θ (x). Lastly, when n ∈ {2, 3}, itsuggests a practical, essentially-exact solution.
mx,θ,tTheorem 2. If n = 1 then mx,θ,t , {ti (x; θ, t)}i=1, andmx,θ,t{ci (x; θ, t)}i=2 have closed forms. Thus, so does T θ (x).

If n > 1 then T θ (x) is given essentially in closed fromin the sense that Eqn. (7) still holds, but a generic expm
is needed (a mild issue) and, more importantly, to find thequantities mentioned in Theorem 2 one needs to solve, forevery x, a number (depending on x) of sequential problems
e 6= γ(ξ) where thee ξ)of the form arg mint>0 γ(expm(tA)
pair (ξ, A) ∈ Ω × Rn×(n+1) changes in each of these problems (while the set of all A’s is shared across different x’s,
we cannot know a-priori which A will appear in each problem). Doing this for n = 1 is easy, fast, and accurate. But
for n ≥ 2, partly since φθ (x, ·) may reenter a cell, this requires tedious bookkeeping, multiple expm calls, and multiple invocations of a numerical solver. This leads to slow
approximated solutions that are hard to implement in GPU.
Inference over latent transformations often requires evaluating T θ for multiple values of θ (e.g., 10,000), and ahigh Npts , where Npts is number of points to be transformed (e.g., the number of pixels in an image). We thus
prefer a more practical alternative and propose a specialized solver for integrating CPA fields. This solver, summarized in Algorithm 1, alternates between the analytic solution (with neither bookkeeping nor need to solve explicitly
mx,θ
mx,θ
for mx,θ,1 , {ci (x; θ, 1)}i=2
and {ti (x; θ, 1)}i=1) and ageneric solver. Thus, our solver is both faster and more accurate than a generic solver since most of the trajectory issolved exactly, while only in small portions of the trajectory does our solver resort to the generic one. Importantly,
the number of expm evaluations required is constant, Nc ,and grows with neither Npts nor t. In fact, since the algorithm is highly accurate, fast, and simple, and since usually
Nc ≪ Npts , we prefer to use it even when n = 1.

Construction of Gaussian smoothness priors. It iseasy to build a Gaussian smoothness prior on V ′ (e.g., we
use a D × D covariance whose correlations decay withinter-cell distances), denoted N (0D×1 , ΣPA ), and then the
V ′ → V map, A 7→ θ = B T vec(A), induces a prior on V:p(θ) = N (0d×1 , ΣCPA ) where ΣCPA = B T ΣPA B. See
Figs. 2 as well as Sup. Mat. for samples from p(θ).

4. Applications and Results
We focus here on 1D/2D; see Sup. Mat. for 3D.
Integration Accuracy. By construction, Algorithm 1 ismore accurate than any generic solver used as its subroutine. We also verified empirically, by comparing the integration error (that we can compute this error at all, is by thevirtue of our closed-form solution) averaged over random
CPA velocity fields. See Sup. Mat. for more details.
Timings. We timed our implementation on an Nvidia
GTX 780 card. Tables 2 and 3 compare our solver and ageneric one. Our implementation can be further optimized;
e.g., we wrapped CUDA calls via python and not, say, C++.
Table 4 shows how Nsteps and Npts affect the timings. Empirically, Nsteps = 10 typically sufficed for good results.
Real-Time diffeomorphic image editing. Our code includes a GUI where a user chooses velocities at some/all
vertices of P. Using Lemma 3, the Gaussian priors on v θ ,and the fast integration, we compute, in real time, the conditional warp, exp(E(v θ |user’s choice)). A video with ademo is included in our Sup. Mat.
Expressiveness. While simple, CPAB transformations
capture a wide range of diffeomorphisms (obviously, likeany other work on diffeomorphisms, we make no claims
about capturing non-diffeomorphisms such as occlusions).
One way to see this is by inspecting samples from p(θ).
M also contains what is sometimes called large deformations (an overloaded term, used in some texts for elementsof G outside its identity component), e.g., a rectangle-tohorseshoe morphing; see Sup. Mat. Expressiveness is alsoSignals

1.5

Unwarping the signals

1.5

s0 1.00.5−6σ1 1.0
0.5

s

u

−4σ

s
s
s
s

0.0
1.5
1.0
1 0.5
0.0
1.5
1.0
2 0.5
0.0
1.5
1.0
3 0.5
0.0
1.5
1.0
4 0.5
0.0
1.5
1.0
5 0.5
0.0
1.5
1.0
0.5
0.0

u
u
u
u

s̄

0.0
1.5
1.0
1 0.5
0.0
1.5
1.0
2 0.5
0.0
1.5
1.0
3 0.5
0.0
1.5
1.0
4 0.5
0.0
1.5
1.0
5 0.5
0.0
1.5
1.0
0.5
0.0

−2σ

100

200

300

400

500

600

700

0.0
1.5
1.0
1 0.5
0.0
1.5
1.0
1 0.5
0.0
1.5
1.0
0.5
0.0
1.5
1.0
1 0.5
0.0
1.5
1.0
1 0.5
0.0
1.5
1.0
1 0.5
0.0

0
+2σ
+4σ
+6σ

ū
0

Eigen warp applied to ū

1.5

s0 1.0
0.5

0

100

200

300

Time

400

500

600

700

Time

0

100

200

300

400

500

600

700

Time

5

Figure 7: Nonlinear time warping in real MoCap Data. Left: a reference signal, s0 , 5 other signals, {si }i=1 , and their mean,P55s̄ = 15 5i=1 si . The mean is smeared since {si }i=1 are misaligned. Center: having inferred the warps, {T θi }i=1 , we unwarp
5
the signals, setting ui = si ◦ T −θi . Note {ui }i=1 are better aligned, as is evident by the details preserved in their mean,P
1 P5
θ̄+jσ1 ξ1ū = 5 i=1 ui . Right: warping ū by T
where θ̄ = 51 5i=1 θi , j ∈ {−6, −4, −2, 0, 2, 4, 6}, ξ1 is the first principal
5component of {θi − θ̄}i=1 , and σ1 is the corresponding standard deviation.
exemplified in the experiments below (Figs. 6 and 7).
Inferring latent transformations. Reasoning over realworld signals requires three oft-confounded stages: a representation; a model (e.g., a likelihood); inference. Whilethese are related it is crucial to distinguish between them;e.g., Horn and Schunk’s representation+model of opticalflow had been long regarded highly inaccurate, till modern
inference rendered it almost state-of-the-art [39]. Focusing
here on representations, we avoid advocating a particularlikelihood or inference type, and instead use, in a coarseto-fine manner, simple choices such as Gaussian likelihoodand either conjugate-gradient optimization or the Metropolis algorithm (where, in the proposal, we use either global
or local moves, utilizing Lemma 3). Without claiming thesechoices are optimal, these turned out to be effective.
Latent transformations arise in two typical settings: withand without known correspondences. In the first, also called(exact/inexact) diffeomorphic point matching [24, 18], one
Npts
has landmark pairs {(xi , yi )}i=1
and seeks T : Ω → Ω
Npts
such that {T (xi ) ≈ yi }i=1 . We note that for n = 1, by
Corollary 1, our ability to easily solve this problem using
CPAB transformations means we can solve a monotonic
regression problem, as we demonstrate in Fig. 4. Particularly, by adding boundary constraints and simulated annealing, we can well fit a CPAB transformation to any reasonable CDF; see also our Sup. Mat. More generally,
but still for n = 1, our method can be used to solve inference/optimization problems over monotonic functions,
CDFs and (cumulative sums of) histograms/filters (while respecting bin ordering). The case where n ∈ {2, 3} appears,
e.g., in landmark-based warping. As shown in Fig. 6, using real images and landmarks from a hands dataset ([37],www2.imm.dtu.dk/pubdb/p.php?403), we infer atransformation between different hands in different poses.
We then apply it to the whole source image, creating newhands, whose appearance belongs to the source hands, withgeometry akin to the destination hands. We also animate thesource image by evaluating φθ for different t’s (including
outside of [0, 1]); e.g., Fig. 6e shows animated frames forthe example in the second row of Fig. 6a-6d. Our Sup. Mat.
contains more examples, including ones that show we canhandle a very large Npts . In the no-correspondences case,there are two subsets of Ω, denoted Usrc and Udst , and two
feature maps, Isrc : Usrc → F and Idst : Udst → F, where
F is a feature space (e.g., color), and one seeks T : Ω → Ωsuch that Idst ◦ T −1 ≈ Isrc (or Isrc ◦ T ≈ Idst ); e.g., thisis the case for time warps and landmark-free image registration. In Fig. 7 we show inferred time warps for MoCap
data [22]. A similar process can be applied to the images inFig. 6, though one may prefer information-theoretic measures to pixelwise likelihood. Lastly, we show an example
for statistics on the inferred transformations. Particularly,
Fig. 7 (right) shows warping ū (defined in the figure’s caption) along the 1st eigen warp (computed using PCA on theinferred velocity fields); i.e., we use Euclidean statistics inthe linear V and map the results to the nonlinear M .

5. Conclusion
We proposed new well-behaved transformations that are
simple, expressive, efficient, and have many applications.
We showed the transformations are computed fast and accurately and that they are highly expressive. Like the spacein [3], our spaces lack closure under composition; but sinceT θ2 ◦ T θ1 ∈ G, our representation can still be used withinframeworks that apply consecutive transformations [43].
While this work focused on representation it was primarily motivated by modeling and inference: e.g., the ability
to quickly run MCMC inference over diffeomorphisms is ahallmark of the method, showing that simplicity pays off.

References[1] S. Allassonnière, A. Trouvé, and L. Younes. Geodesic shooting and diffeomorphic matching via textured meshes. In
EMMCVPR. Springer, 2005. 3, 4
[2] N. Arad, N. Dyn, D. Reisfeld, and Y. Yeshurun. Image warping by radial basis functions: Application to facial expressions. CVGIP, 1994. 4[3] V. Arsigny, O. Commowick, X. Pennec, and N. Ayache. A
log-euclidean polyaffine framework for locally rigid or affine
registration. In BIR. Springer, 2006. 3, 4, 8
[4] M. F. Beg, M. I. Miller, A. Trouvé, and L. Younes. Computing large deformation metric mappings via geodesic flows ofdiffeomorphisms. IJCV, 2005. 3[5] F. L. Bookstein. Principal warps: Thin-plate splines and thedecomposition of deformations. TPAMI, 1989. 4[6] T.-L. Chen and S. Geman. Image warping using radial basis
functions. JAS, 2014. 4[7] T. F. Cootes, C. J. Twining, V. Petrovic, R. Schestowitz, andC. J. Taylor. Groupwise construction of appearance models
using piecewise affine deformations. In BMVC, 2005. 3[8] P. Dupuis, U. Grenander, and M. I. Miller. Variational
problems on flows of diffeomorphisms for image matching.
QAM, 1998. 3[9] S. Durrleman, S. Allassonnière, and S. Joshi. Sparse adaptive
parameterization of variability in image ensembles. IJCV,
2013. 3[10] P. Fletcher. Geodesic regression and the theory of leastsquares on Riemannian manifolds. IJCV, 2012. 4[11] P. Fletcher, C. Lu, and S. Joshi. Statistics of shape via principal geodesic analysis on Lie groups. In CVPR, 2003. 4[12] P. Fletcher, C. Lu, S. Pizer, and S. Joshi. Principal geodesic
analysis for the study of nonlinear statistics of shape. IEEE
TMI, 2004. 4[13] O. Freifeld and M. Black. Lie bodies: A manifold representation of 3D human shape. In ECCV. Springer, 2012. 4[14] O. Freifeld, S. Hauberg, and M. Black. Model transport:
Towards scalable transfer learning on manifolds. In CVPR,
2014. 4[15] U. Grenander. General pattern theory: A mathematical studyof regular structures. Clarendon Press, 1993. 3[16] U. Grenander and M. Miller. Representations of knowledge
in complex systems. JRSS, Series B, 1994. 3[17] U. Grenander and M. Miller. Pattern theory: from representation to inference. Oxford University Press, 2007. 3[18] H. Guo, A. Rangarajan, and S. Joshi. Diffeomorphic point
matching. In Handbook of Mathematical Models in Computer Vision. Springer, 2006. 3, 8
[19] M. T. Harandi, M. Salzmann, and R. Hartley. From manifold
to manifold: geometry-aware dimensionality reduction for
SPD matrices. In ECCV. Springer, 2014. 4[20] S. Hauberg, O. Freifeld, and M. Black. A geometric take on
metric learning. In NIPS, 2012. 4[21] J. Hinkle, P. T. Fletcher, and S. Joshi. Intrinsic polynomials
for regression on riemannian manifolds. JMIV, 2014. 4[22] C. Ionescu, D. Papava, V. Olaru, and C. Sminchisescu. Human3.6m: Large scale datasets and predictive methods for
3d human sensing in natural environments. TPAMI, 2014. 8[23] S. Jayasumana, R. Hartley, M. Salzmann, H. Li, and M. Harandi. Kernel methods on the riemannian manifold of symmetric positive definite matrices. In CVPR, 2013. 4[24] S. C. Joshi and M. I. Miller. Landmark matching via large
deformation diffeomorphisms. IEEE TIP, 2000. 3, 8
[25] J. Kim, C. Liu, F. Sha, and K. Grauman. Deformable spatial
pyramid matching for fast dense correspondences. In CVPR,
2013. 3[26] D. Lin, E. Grimson, and J. Fisher. Modeling and estimating
persistent motion with geometric flows. In CVPR, 2010. 3[27] C. Liu, H.-Y. Shum, and C.-S. Zhang. A two-step approachto hallucinating faces: global parametric model and local
nonparametric model. In CVPR, 2001. 4[28] M. Lorenzi and X. Pennec. Geodesics, parallel transport &one-parameter subgroups for diffeomorphic image registration. IJCV, 2012. 4[29] E. G. Miller and C. Chefd’hotel. Practical non-parametricdensity estimation on a transformation group for vision. In
CVPR, 2003. 4[30] D. Mumford and A. Desolneux. Pattern theory: the stochastic analysis of real-world signals. AK Peters, 2010. 3[31] T. Nir, A. M. Bruckstein, and R. Kimmel.
Overparameterized variational optical flow. IJCV, 2008. 2[32] X. Pennec. Intrinsic statistics on riemannian manifolds: Basic tools for geometric measurements. JMIV, 2006. 4[33] Y. Rubner, C. Tomasi, and L. Guibas. Metric for distributions
with applications to image databases. In CVPR, 1998. 4[34] S. Sommer, F. Lauze, S. Hauberg, and M. Nielsen. Manifold
valued statistics, exact principal geodesic analysis and theeffect of linear approximation. ECCV, 2010. 4[35] A. Srivastava, I. Jermyn, and S. Joshi. Riemannian analysis
of probability density functions with applications in vision.
In CVPR, 2007. 4[36] A. Srivastava, S. H. Joshi, W. Mio, and X. Liu. Statisticalshape analysis: Clustering, learning, and testing. TPAMI,
2005. 4[37] M. B. Stegmann and D. D. Gomez. A brief introduction tostatistical shape analysis, 2002. 8[38] R. Subbarao and P. Meer. Nonlinear mean shift over riemannian manifolds. IJCV, 2009. 4[39] D. Sun, S. Roth, and M. J. Black. Secrets of optical flow
estimation and their principles. In CVPR, 2010. 8[40] R. Szeliski. Computer vision: algorithms and applications.
Springer Science & Business Media, 2010. 2, 3[41] S. Taheri, P. Turaga, and R. Chellappa. Towards viewinvariant expression analysis using analytic shape manifolds.
In AFGR, 2011. 4[42] M. Vaillant, M. I. Miller, L. Younes, and A. Trouvé. Statistics on diffeomorphisms via tangent space representations.
NeuroImage, 2004. 3, 4[43] T. Vercauteren, X. Pennec, A. Perchant, and N. Ayache. Diffeomorphic demons: Efficient non-parametric image registration. NeuroImage, 2009. 3, 4, 8
[44] Q. Xie, S. Kurtek, H. Le, and A. Srivastava. Parallel transportof deformations in shape space of elastic surfaces. In ICCV,2013. 4[45] L. Younes. Spaces and manifolds of shapes in computer vision: An overview. IVC, 2012. 3