Optimal Transport in Imaging Sciences

Talks a MODE Conference, Rennes

Transcript

1. Optimal Transport in Imaging Sciences Gabriel Peyré www.numerical-tours.com
2. 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
3. Overview • Wasserstein Distance • Sliced Wasserstein Distance • Wasserstein Barycenter • Gaussian Texture Models
4. 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
5. ⇥ ⇥ 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
6. 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.
7. 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
8. 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,
9. Example of 2-D Optimal Assignments µ
10. 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
11. Overview • Wasserstein Distance • Sliced Wasserstein Distance • Wasserstein Barycenter • Gaussian Texture Models
12. (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 .
13. 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
14. 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 ⇥ ||
15. 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
16. 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
17. 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.
18. Overview • Wasserstein Distance • Sliced Wasserstein Distance • Wasserstein Barycenter • Gaussian Texture Models
19. µ 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
20. 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
21. t = 0 t = 1 µ0 µ1 Example of Displacement Interpolation µ0 t=0 t=0.2 t=0.4 t=0.6 t=0.8 t=1 t=0 t=0.2 t=0.4 t=0.6 t=0.8 t=1 t=0 t=0.2 t=0.4 t=0.6 t=0.8 t=1 t=0 t=0.2 t=0.4 t=0.6 t=0.8 t=1 t=0 t=0.2 t=0.4 t=0.6 t=0.8 t=1 t=0 t=0.2 t=0.4 t=0.6 t=0.8 t=1 µ1 t=0.6 t=0.8 t=1 t=0 t=0.2 t=0.6 t=0.8
22. 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.
23. Examples: Gaussian Mixtures µ1 µ2 µ3 µ2 µ1 µ3 Wassertstein Sliced-Wassertstein
24. 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;
25. Overview • Wasserstein Distance • Sliced Wasserstein Distance • Wasserstein Barycenter • Gaussian Texture Models
26. 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”
27. 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
28. 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.]
29. 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 )
30. 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
31. 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 ).
32. 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,
33. OT Barycenters µ1 µ3 µ2
34. Optimal transport Euclidean Rao 2-D Gaussian Barycenters
35. Input Spot Noise Barycenters
36. Dynamic Textures Mixing
37. 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