INTERPOLATION BETWEEN IMAGES BY CONSTRAINED
OPTIMAL TRANSPORT
Said Kerrache and Yasushi Nakauchi
Graduate School of Systems and Information Engineering, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Japan
Keywords:
Constrained optimal transport, Monge-Kantorovich, Interpolation between images.
Abstract:
In this paper, the recently proposed technique of constrained optimal transport is used to interpolate between
images under specified constraints. The intensity values in both images are considered as mass distributions,
and a flow of minimum kinetic energy is computed to transport the initial distribution to the final one, while
satisfying specified constraints on the intermediate mass as well as the the velocity or the momentum field.
As an application, the proposed method is used for interpolating between images with a common unchanged
part, as well as under constraint on the volume expansion and contraction. The latter is achieved by imposing
bounds on the divergence of the velocity field of the flow. This constraint is discretized then integrated into the
problem Lagrangian using the augmented Lagrangian method. A variation of the solution is also presented,
where the constraint is decoupled into two constraints coordinated by an additional Lagrange multiplier. This
allows a considerable speedup, though numerical robustness decreases in certain cases. Constrained optimal
transport can potentially be used in image registration. In particular, the proposed method for controlling the
volume change has potential application in registration of images under volume change constraints as it is the
case for medical images depicting muscle movements or those with contrast enhancing structures.
1 INTRODUCTION
Interpolation between images is a basic task used in
many areas of computer graphics and computer vision
(Beier and Neely, 1992; Schaefer et al., 2006; Stich
et al., 2008; Manning and Dyer, 1999; Seitz and Dyer,
1996) . Given two images, the goal of image interpo-
lation is to find a continuous sequence of images that
transforms one image to the other. There is a large
number of interpolation techniques, the most simple
of which is linear interpolation. In practice, however,
more sophisticated techniques are used, notably the
thin-plate spline interpolation method (see (Wolberg,
1998) for a survey). Image interpolation is tightly
linked to image registration. Indeed, many techniques
for image registration achieve the desired transforma-
tion by gradually deforming one image to the other
(Zitova, 2003). Like image registration, the interpo-
lation problem is ill-posed, and consequently, there
may be a large (even infinite theoretically) number of
interpolations that can match two images. However,
it is important to consider only those interpolations
that cause a reasonable deformation to the image. For
instance, transformations that cause extreme distor-
tions of the image, folding for instance, are excluded.
This can be achieved using a number of techniques,
of which the method of optimal transport is of partic-
ular interest to this work . In this method, the images
are considered as mass distributions, and a map that
minimizes the transport cost between the two images
is sought (Haker et al., 2004).
Optimal transport is a field that is attracting an
increasing number of researchers and is increasingly
applied in various fields in mathematics, science and
technology (Villani, 2009). Within the field of image
registration, optimal transport belongs to the fluid reg-
istration category (Bro-Nielsen and Gramkow, 1996),
and it is in many respects related to the optical flow
technique (Beauchemin and Barron, 1995). It has
been applied successfully for registering various types
of images (Haker et al., 2004; Museyko et al., 2009).
Despite this, optimal transport suffers from a number
of shortcomings. First, it offers no control on the in-
termediate images, which may sometimes lead to un-
desired interpolations . Indeed, depending on the field
of application, often there may be some physical (or
medical) considerations that limit the type of applica-
ble transformations (Haber and Modersitzki, 2004b;
Rohlfing et al., 2003). For instance, interpolations
that lead to objects colliding or occupying the same
space may be rejected (see Figure 5 and Section 4 for
an example). This is an example where intermediate
75
Kerrache S. and Nakauchi Y..
INTERPOLATION BETWEEN IMAGES BY CONSTRAINED OPTIMAL TRANSPORT.
DOI: 10.5220/0003375200750084
In Proceedings of the International Conference on Computer Vision Theory and Applications (VISAPP-2011), pages 75-84
ISBN: 978-989-8425-47-8
Copyright
c
 2011 SCITEPRESS (Science and Technology Publications, Lda.)

images are constrained. As another example, muscles
such as the heart can shrink and expand, but within a
certain limit (Haber et al., 2010). In this case, it is the
relative change in volume that must be constrained.
Furthermore, and as it is the case with many fluid
registration techniques, optimal transport can result in
extreme deformations leading to topological changes
or unreasonable expansions or contractions of the dif-
ferent image regions. Figure 1 shows the interpolation
between two synthetic images that differ by a simple
translation using optimal transport. Figure 2, on the
other hand, shows the interpolation between the same
two images using a divergence free flow. The two re-
sults seem only slightly different. However, a closer
look at the divergence of the optimal transport flow
in Figure 3 shows that the image is not translated as
in the case of the divergence-free flow, but instead,
the region of high intensity is expanded, whereas the
region of low intensity is shrunk. This is the result
of minimizing the kinetic energy: only the difference
in mass is transported, not the whole mass in the re-
gion. Here again, limiting the relative change in vol-
ume would avoid such unreasonable transformations.
In this paper, we use the method of constrained
optimal transport (Kerrache and Nakauchi, 2011) to
interpolate between images. Constraints on the in-
terpolation: intermediate images as well as relative
change in volume can be imposed. Mathematically,
the relative change in volume is simply the Jacobian
of the transformation. Imposing a bound on the Jaco-
bian can be achieved by bounding the divergence of
the velocity field as it shown in Section 3. Aside from
the interpolation itself, an additional motivation for
this work is constrained image registration. Indeed,
once the interpolation is computed it becomes possi-
ble to compute the map between the two images by
integrating the velocity field of the flow. This point
will be the subject of a further study. The remain-
der of this paper is organized as follows. Section 2
presents the notion of constrained optimal transport.
Section 3 presents the proposed method. Section 4
contains the experimental results of the proposed al-
gorithms. Finally, Section 5 concludes the paper and
gives some future research directions.
2 CONSTRAINED OPTIMAL
TRANSPORT
Optimal transport admits two distinct formulations.
In the time-independent formulation, the problem
data consist in two spaces, each supporting a proba-
bility measure. The goal is to find a mass- preserving
map between the two spaces that minimizes a certain
transport cost. The cost can be the distance traveled
by each particle of mass, a function of the distance
or otherwise. In the time-dependent formulation, the
supporting space is the same, and the task consists
in transporting an initial mass distribution to a final
configuration. The cost in the time-dependent case is
also related to the distance between source and desti-
nation, but it has a continuous description. Each curve
joining any two points is given a cost, and therefore,
the transport cost is dependent on the path followed
and not only the start and ending points. Under ap-
propriate assumptions (Villani, 2009), the transport
between two mass distributions can be interpreted as
traversing a path in the set of probability measures
over the space under consideration. Finding an op-
timal transport consists then in finding a minimizing
curve in the space of probabilities. This formulation is
used by (Kerrache and Nakauchi, 2011) to introduce
the concept of constrained optimal transport, where
certain curves in the space of probability measures are
declared infeasible and can not be used as transport
paths. This allows to formulate a number of interest-
ing problems. For instance, it is possible to eliminate
certain mass distributions, or eliminate transport plans
that cause certain changes to the initial density. This
is the case in this paper, where transport plans that
cause undesired intermediate images or deformations
are eliminated.
In (Benamou and Brenier, 2000), the problem of
optimal mass transport in a closed convex subset D
of R
d
with the squared Euclidean distance as cost
is transformed to an optimal control problem of a
potential flow (Cohen and Kundu, 2004). The ap-
proach consists in computing a flow, in the sens of
fluid dynamics, having the minimum kinetic energy
that moves the initial density to the final one. More
precisely, the problem is formulated as
inf
ρ,m
Z
1
0
Z
D
|
m(t,x)
|
2
2ρ(t,x)
dxdt, (1)
s.t. ∂
t
ρ + ∇ · m = 0, (2)
ρ(0,·) = ρ
0
, ρ(1,·) = ρ
1
, (3)
where ρ (t,x) is the density, m (t,x) is the momentum
of the flow, ρ
0
(x) and ρ
1
(x) are two bounded positive
density functions defined on D, such as
Z
D
ρ
0
(x)dx =
Z
D
ρ
1
(x)dx = 1
This problem is then transformed to the following
saddle point problem
inf
φ,q
sup
µ
L (φ, q, µ) = F (q) + G (φ) +
h
µ,∇φ − q
i
, (4)
where µ = (ρ,m), G(φ) =
R
D
φ(0,x)ρ
0
(x) −
φ(1,x)ρ
1
(x)dx, F is defined by
VISAPP 2011 - International Conference on Computer Vision Theory and Applications
76

Figure 1: Interpolation using optimal transport.
Figure 2: Interpolation using a divergence-free flow.
Figure 3: The divergence of the optimal transport flow.
INTERPOLATION BETWEEN IMAGES BY CONSTRAINED OPTIMAL TRANSPORT
77

F (q) =

0 if q ∈ K,
+∞ otherwise,
with
K =
n
(a,b) : R × R
d
→ R × R
d
,
a +
|
b
|
2
2
≤ 0 pointwise
)
,
and
h
·,·
i
is the inner product defined by
h
u,v
i
=
Z
1
0
Z
D
u · v. (5)
This formulation allowed the authors to devise an
algorithm to compute the optimal flow by an Uzawa
type algorithm, the description of which can be
found in (Fortin and Glowinski, 1983) and (Glowin-
ski and Tallec, 1989). The maximization over µ,
which is a free variable in Problem (4), implies that
|
∇φ − q
|
must be zero at the saddle point. This al-
lowed the authors to augment Lagrangian by the term
r/2
|
∇φ − q
|
2
, where r is a positive constant. Problem
(4) is thus transformed to
inf
φ,q
sup
µ
L (φ, q, µ) = F (q)+ G (φ) +
h
µ,∇φ − q
i
+
r
2
|
∇φ − q
|
2
. (6)
The authors then present a numerical scheme for solv-
ing Problem (6), where φ, q and µ are updated iter-
atively. The computation of φ consists in solving a
Poisson equation, q is obtained by solving a pointwise
optimization problem, whereas µ is updated using a
gradient-type rule.
In (Kerrache and Nakauchi, 2011), Problem (1)
is solved under the additional constraint that µ =
(ρ,m) ∈ U, where U is a closed convex set. With
this additional constraint,
|
∇φ − q
|
is not guaranteed
to reach 0. In fact,
|
∇φ − q
|
not reaching zero indi-
cates that the constraint is active. Because of this, it is
not possible to use the same approach as previously,
and the problem is rewritten by introducing new de-
coupling variables and Lagrange multipliers. Further-
more, the Lagrangian is augmented with additional
quadratic terms that improve the quality of the numer-
ical algorithms. The new form is (see (Kerrache and
Nakauchi, 2011)):
inf
φ,q,p,b
sup
µ∈U,ν,η
L
r,s
(φ,q, p,b,µ,ν,η) = F (q) + G (φ)
+
h
µ, p − b
i
+
h
ν,∇φ − p
i
+
h
η,b − q
i
+
r
2
|
∇φ − p
|
2
+
r
2
|
b − q
|
2
−
s
2
|
µ − ν
|
2
−
s
2
|
µ − η
|
2
. (7)
Here r and s are positive numbers, p and b are decou-
pling variables and ν, η are the new Lagrange multi-
pliers. The authors then present a number of Uzawa
type algorithms to solve Problem (7).
3 THE PROPOSED METHOD
The interpolation between two images is achieved by
computing a flow of minimum kinetic energy that
transports the intensity values of the first image to
match those in the second one. This of course as-
sumes that both images have the same total amount
of intensity values, which can be achieved by nor-
malizing both images (Haker et al., 2004). Using
constrained optimal transport, it is possible to impose
constraints on the intermediate interpolation. For in-
stance, it is possible to fix a part of the image that is
common to both the initial and the final images during
the interpolation process. An example exploiting this
idea is presented in the section of experiment results.
Similarly, it is possible to isolate different parts of the
image and interpolate them independently by forcing
the normal components of the momentum to be null
at the boundaries of each portion.
A more complex type of constraints can be im-
posed on the divergence of the velocity field, which
translates into a constraint on the rate of expansion
and contraction of the image content. Indeed, the
transformation created by the fluid flow can be com-
puted by integrating in time the path that each particle
follows. A particle initially at (x
0
,y
0
) follows a path
that verifies the following system of ordinary differ-
ential equations:

dx/dt = v
x
(x (t) , y (t),t)
dy/dt = v
y
(x (t) , y (t),t)
with initial conditions (x
0
,y
0
). Here v
x
and v
y
are
the components of the velocity field of the flow. The
map that associates to each point (x, y) the point
(x (t) , y (t)) solution to the above system of equations
is called the flow map and denoted by ϕ(x,t) or ϕ
t
(x).
The transformation caused by the flow is therefore ϕ
1
,
that is the map that associates to each point its final
destination. A classical result in fluid dynamics is that
the Jacobian J (x,t) of ϕ
t
evolves under the relation
(see (Chorin and Marsden, 1993)):
∂
∂t
J (x,t) = ∇ · v (ϕ(x,t),t)J (x,t),
where v is the velocity field of the flow. If we impose
the condition
α ≤ ∇ · v (x,t) ≤ β, ∀x ∈ D, (8)
VISAPP 2011 - International Conference on Computer Vision Theory and Applications
78

and assume that J (x,t) is strictly positive, then it fol-
lows that
α ≤
Z
1
0
∂
∂t
J (x,t)
J (x,t)
dt ≤ β.
Since at t = 0, ϕ is the identity map, we have that
exp(α) ≤ J (x,1) ≤ exp(β). Therefore, controlling the
divergence of the velocity field allows the control of
the Jacobian of the transformation. The advantage of
this approach is that instead of a strong non-linear
constraint on the determinant of the flow map, we
have a much simpler one, which is linear in the veloc-
ity field. Notice that the constraint is enforced point-
wise, which is much stronger than a constraint on the
average.
An important observation is that the mass preser-
vation condition and the constraint on the volume ex-
pansion and contraction can be contradictory. Indeed,
preserving and matching the mass, which in the case
of images represents the intensity, may require a min-
imum change in the volume. Therefore, the bounds
on the Jacobian must be reasonably chosen. Further-
more, the constraint (8) is non-convex as a function
of µ, which implies that local optimization algorithms
as those presented in this paper may find a local mini-
mum only, depending on the initial values of the vari-
ables.
3.1 Optimization
There are two approaches that can be used to solve a
constrained transport problem. The first is to project
the variable µ onto the feasible set at each iteration.
This method is used in Algorithm 1 presented here-
after, which is proposed in (Kerrache and Nakauchi,
2011). In this algorithm and the following ones, V ,
H, H
0
, H
00
are Hilbert spaces of functions and vector
fields over the computational domain equipped with
inner products analogous to (5) and the associated
norms. For fixed µ, p and b, the two Lagrangians L
p,µ
and L
b,µ
appearing in the algorithms are defined by
L
p,µ
(φ,ν) = G (φ) +
h
ν,∇φ − p
i
+
r
2
|
∇φ − p
|
2
−
s
2
|
ν − µ
|
2
,
L
b,µ
(q,η) = F (q) +
h
η,b − q
i
+
r
2
|
q − b
|
2
−
s
2
|
η − µ
|
2
.
These Lagrangians are introduced merely for simpli-
fying notation.
Algorithm 1 (Kerrache and Nakauchi, 2011)
• Initialization:

p
0
,b
0
,µ
0
,ν
0
,η
0

∈ H × H ×U ×
H × H given arbitrarily.
• Repeat the following steps until convergence.
1. Compute φ
n
such that: L
p
n−1
,µ
n

φ
n
,ν
n−1

≤
L
p
n−1
,µ
n

φ,ν
n−1

, ∀φ ∈ V ;
2. p
n
= p
n−1
− ρ
r

µ
n
− ν
n−1
+ r

p
n−1
− ∇φ
n

;
3. ν
n
= ν
n−1
+ ρ
s

∇φ
n
− p
n−1
− s

ν
n−1
− µ
n

;
4. Compute q
n
such that: L
b
n−1
,µ
n

q
n
,η
n−1

≤
L
b
n−1
,µ
n

q,η
n−1

, ∀(q,η) ∈ H × H;
5. b
n
= b
n−1
− ρ
r

η
n−1
− µ
n
+ r

b
n−1
− q
n

;
6. η
n
= η
n−1
+ ρ
s

b
n−1
− q
n
− s

η
n−1
− µ
n

;
7. µ
n+1
= arginf
µ∈U
|
µ−
1
2

ν
n
+ η
n
+
1
s
(p
n
− b
n
)



2
.
As shown in (Kerrache and Nakauchi, 2011), comput-
ing φ
n
is equivalent to solving the Poisson equation
−r∆φ = ∇ ·

ν
n−1
− rp
n−1

, with Neumann boundary
conditions in the time dimension given by r∂
t
φ(0,·) =
ρ
0
− ν
n
0
(0,·) + rp
n−1
0
(0,·) and r∂
t
φ(1,·) = ρ
1
−
ν
n
0
(1,·) + r p
n−1
0
(1,·), whereν
0
and p
0
are the first el-
ements of the vectors ν and p respectively. Similarly,
finding q
n
amounts to solving
inf
q∈K




b
n−1
+
η
n−1
r
− q




2
.
The projection method is efficient for problems
with simple constraints. However, for more complex
constraints that can be written as a set of a reason-
able number of inequalities, a faster method is to use
the augmented Lagrangian technique to enforce the
constraints gradually. The cost of each iteration is
significantly less than that in the projection method,
since only a non-constrained, or bound-constrained,
quadratic minimization is needed instead of a con-
strained one. The use of augmented Lagrangian has
also the advantage of not requiring a feasible start-
ing point, as it is the case for interior point methods.
This is an important advantage, since finding a feasi-
ble point can be very difficult in some problems, for
instance the present image interpolation task. This
method is presented in Algorithm 2. Here, the La-
grangian is augmented with
−
w
2
|
C (µ)
|
2
,
where w is a positive number. For the problem of con-
strained divergence, C (µ) is taken to be ∇ · v − ω for
some specified constant value ω. Notice that the algo-
rithm solves the problem with an equality constraint
C (µ) = 0. Inequality constraints can be handled eas-
ily by adding slack variables, in which case, Step 8
corresponds to a bound-constraint minimization.
INTERPOLATION BETWEEN IMAGES BY CONSTRAINED OPTIMAL TRANSPORT
79

Algorithm 2
• Initialization:

p
0
,b
0
,µ
0
,ν
0
,η
0
,λ
0

∈ H × H ×
H × H ×H × H
0
given arbitrarily.
• Repeat the following steps until convergence.
1. Compute φ
n
such that L
p
n−1
,µ
n

φ
n
,ν
n−1

≤
L
p
n−1
,µ
n

φ,ν
n−1

, ∀φ ∈ V ;
2. p
n
= p
n−1
− ρ
r

µ
n
− ν
n−1
+ r

p
n−1
− ∇φ
n

;
3. ν
n
= ν
n−1
+ ρ
s

∇φ
n
− p
n−1
− s

ν
n−1
− µ
n

;
4. Compute q
n
such that: L
b
n−1
,µ
n

q
n
,η
n−1

≤
L
b
n−1
,µ
n

q,η
n−1

, ∀(q,η) ∈ H × H;
5. b
n
= b
n−1
− ρ
r

η
n−1
− µ
n
+ r

b
n−1
− q
n

;
6. η
n
= η
n−1
+ ρ
s

b
n−1
− q
n
− s

η
n−1
− µ
n

;
7. µ
n+1
= argsup
µ∈H
h
h
µ, p
n
− b
n
i
−
s
2
|
µ − ν
n
|
2
−
s
2
|
µ − η
n
|
2
+
h
λ
n
,C (µ)
i
−
w
2
|
C (µ)
|
2
i
;
8. λ
n+1
= λ
n
− ρ
w
C

µ
n+1

.
In the algorithm, λ denotes the Lagrange multiplier
for the new constraint C (µ) = 0. The numbers ρ
r
,
ρ
s
and ρ
w
are stepping parameters that are chosen by
the user. The choice of these parameters is important.
If they are chosen too large, the algorithm diverges,
whereas if chosen to small, the algorithm takes un-
necessarily long time to converge.
To apply this algorithm, the divergence operator is
first discretized (in the current implementation, using
finite difference), then the minimization is performed
on the resulting discretization. The constraint on the
divergence is applied to the velocity field, which is
related to µ by a nonlinear transformation:
v =
m
ρ
.
As a result, the problem of minimizing µ in Step 8 of
Algorithm 2 can become time consuming when solv-
ing large problems. Indeed, the Hessian is not con-
stant and the number of its non-zero elements grows
quadratically with the size of discretization of the di-
vergence operator. A solution to this problem is to
use the Lagrange multiplier technique to decouple the
relation linking the velocity and µ from the bound
constraints. Accordingly, a new variable γ is added,
which represents the velocity field, and a new con-
straint is added to the problem: m = ργ (remember
that µ = (ρ,m)). This constraint is enforced using
a Lagrange multiplier ζ, and the Lagrangian is aug-
mented with the the quadratic term
−
u
2
|
m − ργ
|
2
,
where u is a positive constant. The new problem can
be solved using the following algorithm.
Algorithm 3
• Initialization:

p
0
,b
0
,µ
0
,ν
0
,η
0
,λ
0
,ζ
0
,γ
0

∈ H ×
H ×H × H ×H ×H
0
× H
00
× H
00
given arbitrarily.
• Repeat the following steps until convergence.
1. Compute φ
n
such that L
p
n−1
,µ
n

φ
n
,ν
n−1

≤
L
p
n−1
,µ
n

φ,ν
n−1

, ∀φ ∈ V ;
2. p
n
= p
n−1
− ρ
r

µ
n
− ν
n−1
+ r

p
n−1
− ∇φ
n

;
3. ν
n
= ν
n−1
+ ρ
s

∇φ
n
− p
n−1
− s

ν
n−1
− µ
n

;
4. Compute q
n
such that: L
b
n−1
,µ
n

q
n
,η
n−1

≤
L
b
n−1
,µ
n

q,η
n−1

, ∀(q,η) ∈ H × H;
5. b
n
= b
n−1
− ρ
r

η
n−1
− µ
n
+ r

b
n−1
− q
n

;
6. η
n
= η
n−1
+ ρ
s

b
n−1
− q
n
− s

η
n−1
− µ
n

;
7. µ
n+1
= argsup
µ∈H
h
h
µ, p
n
− b
n
i
−
s
2
|
µ − ν
n
|
2
−
s
2
|
µ − η
n
|
2
+
h
ζ
n
,m − ργ
n
i
−
u
2
|
m − ργ
n
|
2
i
;
8. γ
n+1
= argsup
γ∈H
0

ζ
n
,m
n+1
− ρ
n+1
γ

−
u
2


m
n+1
− ρ
n+1
γ


2
+
h
λ
n
,C (γ)
i
−
w
2
|
C (γ)
|
2
i
9. λ
n+1
= λ
n
− ρ
w
C

γ
n+1

.
10. ζ
n+1
= ζ
n
− ρ
u

m
n+1
− ρ
n+1
γ
n+1

.
In Algorithm 3, ρ
u
denotes the stepping parameter for
the new Lagrange multiplier, which once again must
be chosen properly to strike a balance between the
need for convergence and convergence speed. With
this new form, the optimization problem for µ is de-
composed into two optimization problems. However,
both are much simpler than the original one. Indeed,
both problems are quadratic, and in the case of µ, the
optimization can be done pointwise. This has the ef-
fect of improving the execution time.
4 EXPERIMENTAL RESULTS
A first test of interpolation by constrained optimal
transport is preformed using the images of Figure 4.
Figure 4: Test images: a building with moving clouds in the
background.
VISAPP 2011 - International Conference on Computer Vision Theory and Applications
80

Image at t=0
Image at t=0.14 Image at t=0.28
Image at t=0.42
Image at t=0.57
Image at t=0.71 Image at t=0.85
Image at t=1
Image at t=0
Image at t=0.14 Image at t=0.28
Image at t=0.42
Image at t=0.57
Image at t=0.71 Image at t=0.85
Image at t=1
Figure 5: Interpolation using free optimal transport.
Image at t=0
Image at t=0.14
Image at t=0.28 Image at t=0.42
Image at t=0.57 Image at t=0.71 Image at t=0.85 Image at t=1
Image at t=0
Image at t=0.14
Image at t=0.28 Image at t=0.42
Image at t=0.57 Image at t=0.71 Image at t=0.85 Image at t=1
Figure 6: Interpolation using constrained optimal transport.
Figure 7: Test images: MRIs of a female breast, left: during the wash-in phase. right: during the wash-out phase (Haber and
Modersitzki, 2004a).
In these two images, the building is a common el-
ement, whereas the background changes. Figure 5
shows the corresponding interpolation using free op-
timal transport. Notice that the building is distorted
by the motion of the clouds. In Figure 6, a more rea-
sonable interpolation is presented, where the build-
ing is constrained to remain constant. The solution is
computed using Algorithm 1, with periodic boundary
conditions.
A test of the constrained divergence interpolation
method is conducted on the test images shown in Fig-
ure 7
1
. These are magnetic resonance images of a
female breast. This kind of images results from in-
jecting a marker into the body that causes an addi-
tional contrast to appear in the final image. Because
1
Figure reproduced with permission, courtesy of IOP
Publishing Ltd and the authors. Original paper available
at http://dx.doi.org/10.1088/0266-5611/20/5/018
INTERPOLATION BETWEEN IMAGES BY CONSTRAINED OPTIMAL TRANSPORT
81

Image at t=0 Image at t=0.14
Image at t=0.28
Image at t=0.42
Image at t=0.57 Image at t=0.71 Image at t=0.85 Image at t=1
Image at t=0 Image at t=0.14
Image at t=0.28
Image at t=0.42
Image at t=0.57 Image at t=0.71 Image at t=0.85 Image at t=1
Figure 8: Interpolation using constrained optimal transport.
Figure 9: Divergence of the velocity field in the free case.
of this, traditional intensity based methods can not
register properly these images. They cause the vol-
ume where the contrast appears to shrink considerably
(see (Haber and Modersitzki, 2004a) for more details
about these images and this phenomenon). The inter-
polation is done on these images with resolution 64 ×
64. For the sake of this example, the absolute value
of the divergence of the velocity field is bounded by
1.1, which limits the Jacobian, and hence the relative
change of volume, between approximately 3 and 1/3
pointwise. Figure 8 shows the resulting interpolation
computed by Algorithm 2. The figure shows that it
is possible to achieve low volume change rates and
mass preservation, while at the same time keeping the
interpolation smooth. Figure 9 and 10 compare the
divergence of the velocity field for the free and con-
strained interpolation respectively. The decrease in
the divergence is clearly visible. Figure 11 and 12
show the evolution of the convergence criteria for the
two algorithms. It can be seen that indeed, the en-
ergy is minimized, the constraints are being satisfied
gradually and that they are active, since
|
∇φ − q
|
does
not approach zero. To compare the two algorithms,
a test is done using 32 × 32 versions of the above
images. The execution time is 152 seconds for Algo-
rithm 3 versus 257 seconds for Algorithm 3, resulting
from a parallel execution on a four-core Xeon pro-
cessor with enough memory. The tolerance is taken
to be 5e-5. Results show that Algorithm 3 is almost
twice as fast as Algorithm 2, however, it is less ro-
bust. The additional decoupling creates an oscillation
phenomenon, which may cause the algorithm to miss
VISAPP 2011 - International Conference on Computer Vision Theory and Applications
82

Figure 10: Divergence of the velocity field in the constrained case.
0 200 400 600 800 1000
10
−6
10
−4
10
−2
10
0
Iteration
Energy
| ∇ ⋅ µ|
|∇φ−q|
|C|
Figure 11: Trace of convergence criteria for Algorithm 2.
0 50 100 150
10
−6
10
−5
10
−4
10
−3
10
−2
10
−1
10
0
10
1
Iteration
Energy
|∇ ⋅ µ|
|∇φ−q|
|C|
|m − ρ v |
Figure 12: Trace of convergence criteria for Algorithm 3.
the saddle point and diverge when the constraint are
too tight, that is small values of the bound on the di-
vergence. Thus for problems with tight constraint, it
is better to use Algorithm 2, whereas for other prob-
lems, using Algorithm 3 results in a considerable gain
in time.
5 CONCLUSIONS
This paper presented a new framework for con-
strained image interpolation using optimal transport.
Algorithms were introduced to solve the problem and
numerical experiments were conducted that demon-
strate their working. The proposed approach of inter-
polation by constrained optimal transport can be gen-
eralized to other types of constraints, as long as they
can be handled by the augmented Lagrangian tech-
nique or the projection method. However, for convex
constraints, such as the constraint on the divergence
of the velocity field presented in this paper, the so-
lution obtained depends heavily on the starting point.
The algorithms may get stuck in a local minimum or
fail to find a feasible point. This practically limits
the minimum value of the absolute value of the diver-
gence that can be achieved. As a future work, we plan
on investigating methods to obtain good initial points
from which to start the optimization, and techniques
to escape local minima and locally infeasible regions.
Also as future direction of research, we intend to
use this method for image registration, especially in
medical applications. As previously explained, this
requires the computation of the flow map. Although
this can theoretically be easily done by integrating the
velocity field, it is a challenging problem in practice,
due to the compressible nature of the flows involved.
INTERPOLATION BETWEEN IMAGES BY CONSTRAINED OPTIMAL TRANSPORT
83

REFERENCES
Beauchemin, S. S. and Barron, J. L. (1995). The computa-
tion of optical flow. ACM Comput. Surv., 27:433–466.
Beier, T. and Neely, S. (1992). Feature-based image meta-
morphosis. SIGGRAPH Comput. Graph., 26:35–42.
Benamou, J.-D. and Brenier, Y. (2000). A computational
fluid mechanics solution to the Monge-Kantorovich
mass transfer problem. Numer. Math., 84(3):375–393.
Bro-Nielsen, M. and Gramkow, C. (1996). Fast fluid
registration of medical images. In Proceedings of
the 4th International Conference on Visualization in
Biomedical Computing, pages 267–276, London, UK.
Springer-Verlag.
Chorin, A. J. and Marsden, J. E. (1993). A Mathematical In-
troduction to Fluid Mechanics (Texts in Applied Math-
ematics) (v. 4). Springer.
Cohen, I. M. and Kundu, P. K. (2004). Fluid Mechanics,
Third Edition. Academic Press.
Fortin, M. e. and Glowinski, R. e. (1983). Augmented La-
grangian methods: Applications to the numerical so-
lution of boundary-value problems. Studies in Math-
ematics and its Applications, 15. Amsterdam-New
York-Oxford: North-Holland. XIX, 340 p.
Glowinski, R. and Tallec, P. L. (1989). Augmented La-
grangian and operator-splitting methods in nonlinear
mechanics. SIAM.
Haber, E., Horesh, R., and Modersitzki, J. (2010). Numer-
ical optimization for constrained image registration.
Numerical Linear Algebra with Applications, 17:343–
359.
Haber, E. and Modersitzki, J. (2004a). Numerical meth-
ods for volume preserving image registration. Inverse
Problems, 20(5):1621.
Haber, E. and Modersitzki, J. (2004b). Volume preserv-
ing image registration. In Barillot, C., Haynor, D. R.,
and Hellier, P., editors, Medical Image Computing
and Computer-Assisted Intervention MICCAI 2004,
volume 3216 of Lecture Notes in Computer Science,
pages 591–598. Springer Berlin / Heidelberg.
Haker, S., Zhu, L., Tannenbaum, A., and Angenent, S.
(2004). Optimal mass transport for registration and
warping. Int. J. Comput. Vision, 60(3):225–240.
Kerrache, S. and Nakauchi, Y. (2011). Computing con-
strained energy-minimizing flows. In 3rd Interna-
tional Conference on Computer Research and Devel-
opment. Accepted.
Manning, R. A. and Dyer, C. R. (1999). Interpolating view
and scene motion by dynamic view morphing. In
Proc. Computer Vision and Pattern Recognition Conf.,
volume 1, pages 388–394.
Museyko, O., Stiglmayr, M., Klamroth, K., and Leuger-
ing, G. (2009). On the application of the monge-
kantorovich problem to image registration. SIAM J.
Img. Sci., 2(4):1068–1097.
Rohlfing, T., Maurer, C. R., Bluemke, D. A., and Jacobs,
M. A. (2003). Volume-preserving nonrigid registra-
tion of MR breast images using free-form deformation
with an incompressibility constraint. Medical Imag-
ing, IEEE Transactions on, 22(6):730–741.
Schaefer, S., McPhail, T., and Warren, J. (2006). Image
deformation using moving least squares. ACM Trans.
Graph., 25:533–540.
Seitz, S. M. and Dyer, C. R. (1996). View morphing. In
SIGGRAPH96, pages 21–30.
Stich, T., Linz, C., Wallraven, C., Cunningham, D., and
Magnor, M. (2008). Perception-motivated interpola-
tion of image sequences. In Proceedings of the 5th
symposium on Applied perception in graphics and vi-
sualization, APGV ’08, pages 97–106, New York, NY,
USA. ACM.
Villani, C. (2009). Optimal transport. Old and new.
Grundlehren der Mathematischen Wissenschaften
338. Berlin: Springer. .
Wolberg, G. (1998). Image morphing: A survey. The Visual
Computer, 14:360–372.
Zitova, B. (2003). Image registration methods: a survey.
Image and Vision Computing, 21(11):977–1000.
VISAPP 2011 - International Conference on Computer Vision Theory and Applications
84