Optimal Transport
in Imaging Sciences
Gabriel Peyré
www.numerical-tours.com
Slide 2
Slide 2 text
Other applications:
Texture synthesis
Image segmentation
Statistical Image Models
Input image Modified image
Source image (X)
Style image (Y)
Source image after color transfer
J. Rabin Wasserstein Regularization
Optimal transport framework Sliced Wasserstein projection Applications
Application to Color Transfer
Source image (X)
Sliced Wasserstein projection of X to style
image color statistics Y
Optimal transport framework Sliced Wasserstein projection Applications
Application to Color Transfer
Source image (X)
Style image (Y)
Sliced Wasserstein projection of X to style
image color statistics Y
Source image after color transfer
J. Rabin Wasserstein Regularization
Colors distribution: each pixel point in R3
Xi
Quotient space: A ne space:
Discrete measure:
Constant weights: pi
= 1/N. Fixed positions Xi
(e.g. grid)
Histogram
Point cloud
{(pi
)
i
\
i
pi
= 1}
RN d/ N
Discrete Distributions
Xi
Rd
i
pi
= 1
µ =
N 1
i=0
pi Xi
Slide 5
Slide 5 text
⇥ ⇥ argmin
N i
||Xi Y (i)
||p
Wasserstein distance:
Projection on statistical constraints:
Proj
C
(f) = Y
Wp
(µX, µY
)p =
i
||Xi Y (i)
||p
Metric on the space of distributions.
C = {f \ µf
= µY
}
Discrete distributions:
Optimal Assignments
Optimal assignment:
Xi
Y (i)
µX
µY
µX
=
1
N
N
i=1
Xi
Slide 6
Slide 6 text
Higher dimensions: combinatorial optimization methods
Hungarian algorithm, auctions algorithm, etc.
µ =
i
pi Xi
⇥ =
i
qi Yi
Wp
(µ, )p solution of a linear program.
Arbitrary distributions:
Computing Transport Distances
Xi
Yi
Explicit solution for 1D distribution (e.g. grayscale images):
sorting the values, O(N log(N)) operations.
O(N5/2 log(N)) operations.
intractable for imaging problems.
Slide 7
Slide 7 text
Probabilistic coupling:
Extension to:
= N2
i=1
qi Yi
weighted distributions
arbitrary number of points
µ,
= P RN1 N2 \ P 0, P1 = p, P 1 = q
If N1
= N2
, permutation matrix:
i
||Xi Y (i)
||p =
i,j
Pi,j
||Xi Yj
||p
Defined even if N1
= N2
.
Takes into account weights.
Linear objective.
Coupling Matrices
µ = N1
i=1
pi Xi
q
p
P = P = (
i (j)
)
i,j
Slide 8
Slide 8 text
Kantorovitch Formulation
Linear programming (Kantorovitch): W(µ, )p = P , C
P argmin
P µ,
P, C =
i,j
Ci,jPi,j
Optimal coupling P :
p q
q
p
P
P
P, C = cst
extremal points:
N , P = P
Theorem: If pi
= qi
= 1/N,
Slide 9
Slide 9 text
Example of 2-D Optimal Assignments
µ
Slide 10
Slide 10 text
Transport:
Equalization:
Optimal assignement:
Input color images: fi
RN 3.
min
N
||f0 f1
⇥ ||
T : f0
(x) R3 f1
( (i)) R3
˜
f0
= T(f0
) ˜
f0
= f1
Color Histogram Equalization
⇥i
=
1
N
x
fi(x)
T
Optimal transport framework Sliced Wasserstein projection Applications
Application to Color Transfer
Source image (X)
Sliced Wasserstein projection
image color statistics Y
Optimal transport framework Sliced Wasserstein projection Applications
Application to Color Transfer
Source image (X)
Style image (Y)
Sliced Wasserstein projection of X to style
image color statistics Y
Source image after color transfer
J. Rabin Wasserstein Regularization
Application to Color Transfer
Source image (X)
Style image (Y)
Sliced Wasserstein projection of X to style
image color statistics Y
Source image after color transfer
J. Rabin Wasserstein Regularization
f0
T(f0
)
µ1
f1
µ0 µ1
f0
T
(p = 2)
Theorem: E(X) = SW(µX, µY
)2 is of class C1 and
SW(µX, µY
)2 =
|| ||=1
W(µX , µY
)2d
E(X) = Xi Y⇥ (i)
, d .
Sliced-Wasserstein Distance
Key idea: replace transport in Rd by series of 1D transport.
Xi
Xi,
Sliced Wasserstein distance:
Projected point cloud: X = { Xi, ⇥}i
.
[Rabin, Peyr´
e, Delon & Bernot 2010]
where
N
are 1-D optimal assignents of X and Y .
Possible to use SW in variational imaging problems.
Fast numerical scheme : use a few random .
Slide 13
Slide 13 text
Stochastic gradient descent of E(X):
E (X) = W(X , Y )2
Step 1: choose at random.
Step 2:
X( ) converges to C = {X \ µX
= µY
}.
X( +1) = X( ) E (X( ))
Sliced Wasserstein Transportation
N , X = Y
Theorem:
X( ) Proj
C
(X(0))
Numerical observation:
Final assignment
X
is a local minimum of
E
(
X
) =
SW
(
µX , µY )
2
Slide 14
Slide 14 text
Solving is computationally untractable.
and is close to f0
f(0) = f0
(Stochastic) gradient descent:
At convergence: µ ˜
f0
= µf1
f( +1) = f( ) E(f( ))
˜
f0
solves min
f
E(f) = SW(µf , µf1
)
f( ) ˜
f0
Sliced Wasserstein Transfer
Approximate Wasserstein projection:
min
N
||f0 f1
⇥ ||
Slide 15
Slide 15 text
Input image f0
Target image f1
Transfered image ˜
f0
Color Transfer
Application to Color Transfer
Source image (X)
Style image (Y)
Sliced Wasserstein projection of X to style
image color statistics Y
Source image after color transfer
Application to Color Transfer
Source image (X)
Style image (Y)
Sliced Wasserstein projection of X to style
image color statistics Y
Source image after color transfer
Source image (X)
Style image (Y)
Sliced Wasserstein projection of X to style
image color statistics Y
Source image after color transfer
Slide 16
Slide 16 text
Transfered image ˜
f1
Color Exchange
Source image (X)
X ⇥ Y
Source image (X)
Style image (Y)
X ⇥ Y
Source image (X)
Style image (Y)
X ⇥ Y
Source image (X)
Style image (Y)
X ⇥ Y
Y ⇥ X
Input image f0
Target image f1
Transfered image ˜
f0
Slide 17
Slide 17 text
Transfer ˜
f0
= T(f0
).
T : R3 R3 not regular.
Color Transfer Artifacts
Application to Color Transfer
Source image (X)
Sliced Wasserstein projection of X to style
image color statistics Y
Source image after color transfer
Source image (X)
Style image (Y)
Sliced Wasserstein projection of X to style
image color statistics Y
Source image after color transfer
ansport framework Sliced Wasserstein projection Applications
plication to Color Transfer
Source image (X)
Sliced Wasserstein projection of X to style
image color statistics Y
Source image after color transfer
Input image f0
Target image f Transfered image ˜
f0
T amplifies noise.
µ exists and is unique.
µ
µ1
µ3
Theorem:
W2
(µ1
, µ )
W
2 (µ
2 ,µ )
W2
(µ3
,µ
)
If µi
=
Xi
, then µ =
X
X =
i
iXi
Generalizes Euclidean barycenter.
i
i
= 1
Barycenter of {(µi, i
)}L
i=1
:
µ argmin
µ
L
i=1
iW2
(µi, µ)2
if
µ1 does not vanish on small sets,
Wasserstein Barycenter
[Agueh, Carlier, 2010]
µ2
Slide 20
Slide 20 text
Barycenter:
| (k, ) \ Pk,⇥
⇥= 0| N0
+ N1
1.
0 1 t
µ0 µ1
µt
P
µ0
= N1
k=1
p0
(k)
Xk
P argmin
P µ0,µ1 k,
Pk,
||Xk Y ||2
µt =
X
k,`
P?
k,` (1 t)Xk+tY`
Special Case: 2 Distributions
Theorem: [Folklore]
µ0
µ1
µ1
= N2
=1
p1
( )
Y
Gradient descent:
Disadvantage:
Non-convex problem local minima.
Advantages:
µX
that solves
min
X
i
i SW(µX, µXi
)2
Ei
(X) = SW(µX, µXi
)2
X( +1) = X( ) ⇥
L
i=1
i Ei
(X( ))
Sliced Wasserstein Barycenter
µX
is a sum of N Diracs.
Sliced-barycenter:
Smooth optimization problem.
Raw image
sequence
Compute
Wasserstein
barycenter
Project on
the barycenter
Application: Color Harmonization
. Step 1: compute Sliced-Wasserstein Barycenter of color statistics;
. Step 2: compute Sliced-Wasserstein projection of each image onto the
Barycenter;
Raw image sequence
J. Rabin – GREYC, University of Caen Approximate Wasserstein Metric for Texture Synthesis and Mixing
. Step 1: compute Sliced-Wasserstein Barycenter of color statistics;
. Step 2: compute Sliced-Wasserstein projection of each image onto the
Barycenter;
Raw image sequence
J. Rabin – GREYC, University of Caen Approximate Wasserstein Metric for Texture Synthesis and Mixing
. Step 1: compute Sliced-Wasserstein Barycenter of color statistics;
. Step 2: compute Sliced-Wasserstein projection of each image onto the
Barycenter;
Raw image sequence
J. Rabin – GREYC, University of Caen Approximate Wasserstein Metric for Texture Synthesis and Mixing
Wasserstein Barycenter Sliced Wasserstein Barycenter Experimental results Applications Conclusion
Color transfer
Color harmonization of several images
. Step 1: compute Sliced-Wasserstein Barycenter of color statistics;
. Step 2: compute Sliced-Wasserstein projection of each image onto the
Barycenter;
asserstein Barycenter Sliced Wasserstein Barycenter Experimental results Applications Conclusion
Color transfer
Color harmonization of several images
. Step 1: compute Sliced-Wasserstein Barycenter of color statistics;
. Step 2: compute Sliced-Wasserstein projection of each image onto the
Barycenter;
serstein Barycenter Sliced Wasserstein Barycenter Experimental results Applications Conclusion
olor transfer
Color harmonization of several images
. Step 1: compute Sliced-Wasserstein Barycenter of color statistics;
. Step 2: compute Sliced-Wasserstein projection of each image onto the
Barycenter;
analysis synthesis
Probability
distribution
µ = µ(p)
Exemplar f0
Outputs f µ(p)
Gaussian models: µ = N(m, ), parameters p = (m, ).
Texture Synthesis
Problem: given f0
, generate f
“random”
perceptually “similar”
Slide 27
Slide 27 text
Input exemplar:
d = 1 (grayscale), d = 3 (color)
N1
N2
N3
Images Videos
Gaussian model:
m RN d, RNd Nd
X µ = N(m, )
Texture synthesis:
given (m, ), draw a realization f = X( ).
highly under-determined problem.
Factorize = AA (e.g. Cholesky).
Compute f = m + Aw where w drawn from N(0, Id).
Texture analysis: from f0
RN d, learn (m, ).
f0
RN d
Gaussian Texture Model
N1
N2
Slide 28
Slide 28 text
Stationarity hypothesis: X(· + ) X
(periodic BC)
Block-diagonal Fourier covariance:
i,j
=
1
N
x
f0
(i + x) f0
(j + x) Rd d
ˆ
y( ) = ˆ( ) ˆ
f( )
y = f computed as
where ˆ
f( ) =
x
f(x)e
2ix1⇥1
N1
+ 2ix2⇥2
N2
= 0, ˆ( ) = ˆ
f0
( ) ˆ
f0
( ) Cd d
is a spot noise = 0, ˆ( ) is rank-1.
MLE of :
Maximum likelihood estimate (MLE) of m from f0
:
i, mi
=
1
N
x
f0
(x) Rd
Spot Noise Model [Galerne et al.]
Slide 29
Slide 29 text
Cd C
Input f0
RN 3 Realizations f
Example of Synthesis
Synthesizing f = X( ), X N(m, ):
= 0, ˆ
f( ) = ˆ
f0
( ) ˆ
w( )
Convolve each channel with the same white noise.
w N(N 1, N 1/2Id
N
)
Slide 30
Slide 30 text
Input distributions (µ0, µ1
) with µi
= N(mi, i
).
E0 E1
W2
(µ0, µ1
)2 = tr (
0
+
1
2
0,1
) + ||m0 m1
||2,
T
0,1
= ( 1/2
1 0
1/2
1
)1/2
S = 1/2
1
+
0,1
1/2
1
T(x) = Sx + m1 m0
where
Ellipses: Ei
= x Rd \ (mi x) 1
i
(mi x) c
Theorem:
If
⇢
ker(⌃0) \ ker(⌃1)? = {0},
ker(⌃1) \ ker(⌃0)? = {0},
Gaussian Optimal Transport
Slide 31
Slide 31 text
The set of Gaussians is geodesically convex:
µt
= ((1 t)Id + tT) µ0
= N(mt, t
)
Gaussian Geodesics
µ1
mt
= (1 t)m0
+ tm1
t
= [(1 t)Id + tT]
0
[(1 t)Id + tT]
µ0
0,1
= ( 1/2
1 0
1/2
1
)1/2
T(x) = Sx + m1 m0
S = 1/2
1
+
0,1
1/2
1
Input distributions (µ0, µ1
) with µi
= N(mi, i
).
Slide 32
Slide 32 text
t
f[0]
f[1]
Geodesic of Spot Noises
0 1
Theorem:
i.e. ˆ
i
( ) = ˆ
f[i]( ) ˆ
f[i]( ) .
f[t] = (1 t)f[0] + tg[1]
ˆ
g[1]( ) = ˆ
f[1]( )
ˆ
f[1]( ) ˆ
f[0]( )
| ˆ
f[1]( ) ˆ
f[0]( )|
Then t [0, 1], µt
= µ(f[t])
Let for i = 0, 1, µi
= µ(f[i]) be spot noises,
Slide 33
Slide 33 text
OT Barycenters
µ1 µ3
µ2
Slide 34
Slide 34 text
Optimal transport
Euclidean Rao
2-D Gaussian Barycenters
Slide 35
Slide 35 text
Input
Spot Noise Barycenters
Slide 36
Slide 36 text
Dynamic Textures Mixing
Slide 37
Slide 37 text
Fast sliced approximation.
Wasserstein distance approach
Extension to a wide range
of imaging problems.
in out
Color transfert, segmentation, . . .
Conclusion
14 Anonymous
P(⇥⇥
) P(⇥ ) P(⇥ c
) P(⇥⇥
) P(⇥ ) P(⇥ c
)
Image modeling with statistical constraints
Colorization, synthesis, mixing, . . .
Source image (X)
Style image (Y)
Source image a
J. Rabin Wasserstein Regularization