Assimilation of image structures
Nicolas Papadakis1,2,
Vincent Chabot2, Alexandros Makris2, M¨
aelle Nodet2, Arthur Vidard2
1 CNRS, Institut de Math´
ematiques de Bordeaux
2 ´
Equipe Moise, Laboratoire Jean Kuntzmann/Inria Grenoble
GRETSI
September 5th 2013
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard)

Scientiﬁc context
Theme: Mathematical and numerical methods for modelling and
understanding the evolution of geophysic ﬂuids: ocean, atmosphere,
ice...
Applications: short and long term forecasting, risk analysis...
Today’s subject: Use of satellite images to monitor numerical models of
the ocean: data assimilation
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 1/1

Scientiﬁc context
Theme: Mathematical and numerical methods for modelling and
understanding the evolution of geophysic ﬂuids: ocean, atmosphere,
ice...
Applications: short and long term forecasting, risk analysis...
Today’s subject: Use of satellite images to monitor numerical models of
the ocean: data assimilation
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 1/1

Overview
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard)

Overview
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard)

Data assimilation in geophysics
Three main types of information
Numerical model: temporal evolution of physical quantities (velocity,
temperature, sea height, pressure...)
Observations: obtained from stations, ballons, driftbuoys (sparse/dense
in space/time)
Uncertaincies: conﬁdence in the background, the model and the
observations (modelling of mesurement errors)
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 2/1

Data assimilation in geophysics
Basic concepts of data assimilation
Estimate the physic state X(t) in a time interval T = [0; 1]
Use of a background prior X0 and observations Y(t) available in T
Minimization of a functional:
J(X) = ||X(0) − X0||2
B
+
T
||Y(t) − H(X(t), t)||2
R(t)
subject to an evolution law ∂tX = M(X), the model
B and R(t): covariance error matrices of background and observations
H and M can be non linear operators: non convex problem
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 3/1

Data assimilation in geophysics
Basic concepts of data assimilation
Estimate the physic state X(t) in a time interval T = [0; 1]
Use of a background prior X0 and observations Y(t) available in T
Minimization of a functional:
J(X) = ||X(0) − X0||2
B
+
T
||Y(t) − H(X(t), t)||2
R(t)
subject to an evolution law ∂tX = M(X), the model
B and R(t): covariance error matrices of background and observations
H and M can be non linear operators: non convex problem
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 3/1

Data assimilation in geophysics
Basic concepts of data assimilation
Estimate the physic state X(t) in a time interval T = [0; 1]
Use of a background prior X0 and observations Y(t) available in T
Minimization of a functional:
J(X) = ||X(0) − X0||2
B
+
T
||Y(t) − H(X(t), t)||2
R(t)
subject to an evolution law ∂tX = M(X), the model
B and R(t): covariance error matrices of background and observations
H and M can be non linear operators: non convex problem
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 3/1

Sequential data assimilation [Kalman, 1960]
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 4/1

Sequential data assimilation [Kalman, 1960]
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 4/1

Sequential data assimilation [Kalman, 1960]
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 4/1

Variational data assimilation [Le Dimet, 1982]
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 5/1

Variational data assimilation [Le Dimet, 1982]
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 5/1

Variational data assimilation [Le Dimet, 1982]
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 5/1

Overview
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard)

Image observation
Objectives
Assimilation of high-resolution satellite images of the ocean
Altimetry SST Chlorophyllis
Observation of sub-mesoscale physics through the evolution of structures
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 6/1

Image observation
Problems
Inhomogeneous spatio-temporal resolution, occlusions
SST image Temporal variation
Source: [B´
er´
eziat and Herlin, 2011]
Dimension of the observation space
Modeling of observation errors
Observation of quantities not deﬁned in the numerical models
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 7/1

Image observation
Problems
Inhomogeneous spatio-temporal resolution, occlusions
SST image Temporal variation
Source: [B´
er´
eziat and Herlin, 2011]
Dimension of the observation space
Modeling of observation errors
Observation of quantities not deﬁned in the numerical models
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 7/1

How to compare images and models?
Indirect observation or Pseudo-observation
Extraction of information from images (contours, motion) that can be
compared to the variables of the model [P. and M´
emin, 2007; Michel, 2011...]
⇒ Complex error modelling
Direct observation
Deﬁne a suitable norm for comparing structures of model X and image
data Y:
J(X) = ||X(0) − X0||2
B
+
T
Y(t) − H(X(t), t)2
R(t)
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 8/1

How to compare images and models?
Indirect observation or Pseudo-observation
Extraction of information from images (contours, motion) that can be
compared to the variables of the model [P. and M´
emin, 2007; Michel, 2011...]
⇒ Complex error modelling
Direct observation
Deﬁne a suitable norm for comparing structures of model X and image
data Y:
J(X) = ||X(0) − X0||2
B
+
T
Y(t) − H(X(t), t)2
R(t)
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 8/1

How to compare images and models?
Indirect observation or Pseudo-observation
Extraction of information from images (contours, motion) that can be
compared to the variables of the model [P. and M´
emin, 2007; Michel, 2011...]
⇒ Complex error modelling
Direct observation
Deﬁne a suitable norm for comparing structures of model X and image
data Y:
J(X) = ||X(0) − X0||2
B
+
T
||Y(t) − H(X(t), t)||2
R(t)
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 8/1

How to compare images and models?
Indirect observation or Pseudo-observation
Extraction of information from images (contours, motion) that can be
compared to the variables of the model [P. and M´
emin, 2007; Michel, 2011...]
⇒ Complex error modelling
Direct observation
Deﬁne a suitable norm for comparing structures of model X and image
data Y:
J(X) = ||X(0) − X0||2
B
+
T
||Y(t) − H(X(t), t)||2
R(t)
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 8/1

Overview
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard)

Direct observation: Past works
Deﬁnition of a norm between image and variable of the model
Link of the temporal variation of images Y(t) with surface velocity w:
T
||∂tY + ∇Y · w||2
[P. and M´
emin, 2008; B´
er´
eziat and Herlin, 2011; Beyou, Cuzol, Gorthi and M´
emin, 2013]
Increase of the state space with a passive tracer Z representing the
visualized phenomena Y:
T
||H(Y) − H(Z)||2
H is the operator mapping images in a suitable space
Couple the tracer dynamics with the model velocity:
∂tZ = M(Z, w)
[Titaud, Vidard, Souopgui and Le Dimet, 2009]
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 9/1

Direct observation: Past works
Deﬁnition of a norm between image and variable of the model
Link of the temporal variation of images Y(t) with surface velocity w:
T
||∂tY + ∇Y · w||2
[P. and M´
emin, 2008; B´
er´
eziat and Herlin, 2011; Beyou, Cuzol, Gorthi and M´
emin, 2013]
Increase of the state space with a passive tracer Z representing the
visualized phenomena Y:
T
||H(Y) − H(Z)||2
H is the operator mapping images in a suitable space
Couple the tracer dynamics with the model velocity:
∂tZ = M(Z, w)
[Titaud, Vidard, Souopgui and Le Dimet, 2009]
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 9/1

Direct observation: Past works
Compare its amplitude to the image intensity (pixels, curvelets)
[Titaud, Vidard, Souopgui and Le Dimet, 2009]
Sequence of ﬂuorescein images
CORIOLIS experimental turntable (Grenoble, France)
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 10/1

Direct observation: Current works
Compare its position to the image level lines (gradients, normals)
[Ba and Fablet, 2010; Chabot, Nodet, P. and Vidard, 2013]
Sequence of ﬂuorescein images
CORIOLIS experimental turntable (Grenoble, France)
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 10/1

Experimental framework
State variable: X = {2D velocity w = (u, v), sea height h, tracer Z}
2D Shallow-water model:
∂tu − (f + ζ)v + ∂xB = −r∗u + κ∆u
∂tv + (f + ζ)u + ∂yB = −r∗v + κ∆v
∂th + ∂x(hu) + ∂y(hv) = 0.
Coupling the tracer evolution with the model velocities:
∂tZ + u∂xZ + v∂yZ = ν∆Z
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 11/1

Experimental framework
State variable: X = {2D velocity w = (u, v), sea height h, tracer Z}
2D Shallow-water model:
∂tu − (f + ζ)v + ∂xB = −r∗u + κ∆u
∂tv + (f + ζ)u + ∂yB = −r∗v + κ∆v
∂th + ∂x(hu) + ∂y(hv) = 0.
Coupling the tracer evolution with the model velocities:
∂tZ + u∂xZ + v∂yZ = ν∆Z
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 11/1

Experimental framework
State variable: X = {2D velocity w = (u, v), sea height h, tracer Z}
2D Shallow-water model:
∂tu − (f + ζ)v + ∂xB = −r∗u + κ∆u
∂tv + (f + ζ)u + ∂yB = −r∗v + κ∆v
∂th + ∂x(hu) + ∂y(hv) = 0.
Coupling the tracer evolution with the model velocities:
∂tZ + u∂xZ + v∂yZ = ν∆Z
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 11/1

Observation space and error modelling
Study of different observation spaces for comparing images:
Data compression: wavelet transform T (Z) and coefﬁcient thresholding
Structure position: gradients ∇Z, normals ∇Z/||∇Z||
Twin experiments to study the robustness to data noise
Modeling of additive or multiplicative noise into the observation
covariance matrices R
Groundtruth Additive noise Multiplicative noise
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 12/1

Observation space and error modelling
Study of different observation spaces for comparing images:
Data compression: wavelet transform T (Z) and coefﬁcient thresholding
Structure position: gradients ∇Z, normals ∇Z/||∇Z||
Twin experiments to study the robustness to data noise
Modeling of additive or multiplicative noise into the observation
covariance matrices R
Groundtruth Additive noise Multiplicative noise
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 12/1

Results
Observation space vs amplitude of data noise
Wavelets Pixels Gradients
14.8 dB 60.1% 60.8% 34.0%
20.8 dB 28.5% 26.2% 17.8%
26.8 dB 17.1% 15.6% 12.4%
Table : Decrease of the velocity RMS error w.r.t the background [Chabot, Nodet,
P. and Vidard, 2013]
More details to improve the results of the above table:
Please go to the wonderful poster of Vincent Chabot tomorrow
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 13/1

Results
Observation space vs amplitude of data noise
Wavelets Pixels Gradients
14.8 dB 60.1% 60.8% 34.0%
20.8 dB 28.5% 26.2% 17.8%
26.8 dB 17.1% 15.6% 12.4%
Table : Decrease of the velocity RMS error w.r.t the background [Chabot, Nodet,
P. and Vidard, 2013]
More details to improve the results of the above table:
Please go to the wonderful poster of Vincent Chabot tomorrow
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 13/1

Overview
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard)

Optimal transport
Let ρ0 and ρ1 be two positive densities
M is the the set of transport maps pushing forward ρ0 to ρ1:
ρ0(x) = ρ1(M(x))| det(∂M(x))|
Deﬁne a cost C(M) associated to each transport map M ∈ M:
C(M) = ||x − M(x)||2
Wasserstein distance W2(ρ0, ρ1): minimal cost C(M∗)
Optimal transport map M∗: unique
Associated geodesic path: trajectories in straight lines
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 14/1

Optimal transport
Let ρ0 and ρ1 be two positive densities
M is the the set of transport maps pushing forward ρ0 to ρ1:
ρ0(x) = ρ1(M(x))| det(∂M(x))|
Deﬁne a cost C(M) associated to each transport map M ∈ M:
C(M) = ||x − M(x)||2
Wasserstein distance W2(ρ0, ρ1): minimal cost C(M∗)
Optimal transport map M∗: unique
Associated geodesic path: trajectories in straight lines
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 14/1

Optimal transport
Let ρ0 and ρ1 be two positive densities
M is the the set of transport maps pushing forward ρ0 to ρ1:
ρ0(x) = ρ1(M(x))| det(∂M(x))|
Deﬁne a cost C(M) associated to each transport map M ∈ M:
C(M) = ||x − M(x)||2
Wasserstein distance W2(ρ0, ρ1): minimal cost C(M∗)
Optimal transport map M∗: unique
Associated geodesic path: trajectories in straight lines
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 14/1

Optimal transport
Wasserstein distance for comparing images
L2 interpolation W2 interpolation
+ A pertinent distance for comparing image structures
- Computational cost for signal of dimension > 1
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 15/1

Optimal transport
Wasserstein distance for comparing images
L2 interpolation W2 interpolation
+ A pertinent distance for comparing image structures
- Computational cost for signal of dimension > 1
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 15/1

Wassersein distance for image assimilation
A symmetric alternative to previous optical ﬂow-based approaches
Images seen as densities, but can be generalized [Maitre and Lombardi, 2013]
Example: Estimation of velocity (Non realistic application of the
Kalman Filter...)
Groundtruth scenario: a 1D density is translated with a constant speed wc
Model state: a density Z(x, t) and its velocity w(x, t)
Dynamics: Z(x + w(x, t), t + 1) = Z(x, t)
Observations: snapshots Y(x, t) of the groundtruth density
Kalman ﬁltering of (Z, w) using Y to recover the velocity value wc:
L2 distance ||Y − Z||2
OT distance W(Y, Z)2
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 16/1

Wassersein distance for image assimilation
A symmetric alternative to previous optical ﬂow-based approaches
Images seen as densities, but can be generalized [Maitre and Lombardi, 2013]
Example: Estimation of velocity (Non realistic application of the
Kalman Filter...)
Groundtruth scenario: a 1D density is translated with a constant speed wc
Model state: a density Z(x, t) and its velocity w(x, t)
Dynamics: Z(x + w(x, t), t + 1) = Z(x, t)
Observations: snapshots Y(x, t) of the groundtruth density
Kalman ﬁltering of (Z, w) using Y to recover the velocity value wc:
L2 distance ||Y − Z||2
OT distance W(Y, Z)2
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 16/1

Wassersein distance for image assimilation
A symmetric alternative to previous optical ﬂow-based approaches
Images seen as densities, but can be generalized [Maitre and Lombardi, 2013]
Example: Estimation of velocity (Non realistic application of the
Kalman Filter...)
Groundtruth scenario: a 1D density is translated with a constant speed wc
Model state: a density Z(x, t) and its velocity w(x, t)
Dynamics: Z(x + w(x, t), t + 1) = Z(x, t)
Observations: snapshots Y(x, t) of the groundtruth density
Kalman ﬁltering of (Z, w) using Y to recover the velocity value wc:
L2 distance ||Y − Z||2
OT distance W(Y, Z)2
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 16/1

Wassersein distance for image assimilation
Sequential assimilation
Perfect data
Observations Result L2 Result OT
Velocity error ||w(., t) − wc)|| along time
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 17/1

Wassersein distance for image assimilation
Sequential assimilation
Noisy data
Observations Result L2 Result OT
Velocity error ||w(., t) − wc)|| along time
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 18/1

Wassersein distance for image assimilation
Sequential assimilation
Very noisy data
Observations Result L2 Result OT
Velocity error ||w(., t) − wc)|| along time
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 19/1

A last advantage of Optimal Transport
2D OT
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 20/1

A last advantage of Optimal Transport
2D OT in a complex domain [P., Peyr´
e and Oudet, 2013]
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 21/1

Conclusion and Perspectives
Preliminary works on image distances
Taking into account the structures contained in the images:
⇒ New kind of observation for the forecasting of geophysic ﬂuids
Speed up of OT algorithms
Test on non-experimental data
Operational models: sQG, NEMOVAR, ROMS
Compatibility of image data with classical obervations
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 22/1

Conclusion and Perspectives
Preliminary works on image distances
Taking into account the structures contained in the images:
⇒ New kind of observation for the forecasting of geophysic ﬂuids
Speed up of OT algorithms
Test on non-experimental data
Operational models: sQG, NEMOVAR, ROMS
Compatibility of image data with classical obervations
Assimilation of image structures (N. Papadakis, V. Chabot, A. Makris, M. Nodet, A. Vidard) 22/1