Anthropomorphic image reconstruction via optimal control and hypoelliptic diffusion

22c721aa043f752b3b6e3299df04b306?s=47 GdR MOA 2015
December 03, 2015

Anthropomorphic image reconstruction via optimal control and hypoelliptic diffusion

by U. Boscain


GdR MOA 2015

December 03, 2015


  1. Anthropomorphic Image Reconstruction via Optimal Control and Hypoelliptic Diffusion Ugo

    Boscain CNRS, CMAP, Ecole Polytechnique, Paris, France and INRIA team GECO Jean-Paul Gauthier, LSIS, Toulon Dario Prandi, Paris Dophine Alexei Remizov, CMAP, Ecole Polytechnique Roman Chertovskih, Universidade do Porto work supported by the ERC grant GECOMETHODS December 3, 2015
  2. An algorithm of image reconstruction based on a model of

    the mammals visual cortex V1 65% of corruption
  3. 85% of corruption

  4. the Sub-Riemannian model of the visual cortex V1 Hubel and

    Wiesel (1959), Hoffman (’89), Petitot (’99), Citti and Sarti (2003), our research group. I will mainly present the final version of the model studied in [1] U. Boscain, J. Duplaix, J.P. Gauthier, F. Rossi, Anthropomorphic image reconstruction via hypoelliptic diffusion. SIAM J. CONTROL OPTIM.Vol. 50, No. 3, pp. 1309-1336, 2012. [2] U. Boscain, R. Chertovskih, J.P. Gauthier, A. Remizov. Hypoelliptic diffusion and human vision: a semi-discrete new twist. SIAM Journal on Imaging Sciences 2014, Vol. 7, No. 2, pp. 669-695 very interesting mathematically sub-Riemannian geometry (optimal control problems linear in the control and with quadratic cost) theory of hypoelliptic diffusion (in the sense of H¨ ormander)
  5. →Our original purpose was to write a model in a

    coherent mathematical language to make prediction that compare with psycological experiments could validate/invalidate the model →But the results are now at the state of the art (and beyond) in image processing
  6. Plan of the talk Plan of the talk why the

    visual cortex V1 has a structure of sub-Riemannian manifold reconstructing level sets via geodesics reconstruction of complex images: the hypoelliptic diffusion model new ideas to reconstruct deeply corrupted images (dynamic restoration)
  7. Two ideas coming from neurophyology of the visual cortex V1

    A. Hubel and Wiesel (1959) observed that there are (groups of) neurons sensitive to positions and directions. Hence the visual cortex lifts an image on the PT R2 = R2 × P1 (experimental fact) B. an image is reconstructed by minimizing the energy necessary to excite groups of neurons that are not excited by the image in PT R2 (postulate)
  8. A1. The lift in PTR2 Hubel and Wiesel (Nobel prize

    1981) observed that there are (groups of) neurons sensitive to directions.
  9. A2. The lift in PTR2 and sensible to "aligned" orientations

    curve visual cortex V1 activation activation orientation columns Plane of the image connections among orientation columns in the same hypercolumn (vertical) hypercolumns connections among orientation columns belonging to different hypercolumns (horizontal)
  10. A3. The lift in PTR2 the visual cortex stores an

    image as a set of points and tangent directions, i.e. it makes a lift to PT R2 = R2 × P1. The projective tangent bundle of R2 (or bundle of directions of the plane). π/2 x2 α ∈ [0, π]/ ∼ α x1 P1 −π/2 PT R2 can be seen as a fiber bundle whose base is R2 and whose fiber at the point (x1, x2) is the set of straight lines (i.e. directions without orientation) d(x1,x2) attached to (x1, x2).
  11. A2. The lift of a curve (x1(t), x2(t)) of class

    C1 curve in R2, such that ( ˙ x1(t), ˙ x2(t)) ̸= 0 ∀ t lift ↓ (x1(t), x2(t), α(t)), curve in R2 × P1 s.t. α(t) = arctan ˙ x2(t) ˙ x1(t) ∈ [−π/2, π/2]/ ∼ Example: (cos(t), sin(t)): -1 -0.5 0 0.5 1 x 0 0.2 0.4 0.6 0.8 1 y 0 1 2 3 Θ -1 -0.5 0 0.5 x 0 0.2 0.4 0.6 0.8 y →every C1 submanifold of R2 has a lift. →not all curves in R2 × P1 are lift of planar curves
  12. A3. Which curves are lift of planar curves? A curve

    in (x1(t), x2(t), α(t)) in R2 × P1 is the lift of a planar curve if α(t) = arctan ˙ x2(t) ˙ x1(t) (1) Here (x1(t), x2(t)) should be of class C1 with ( ˙ x1(t), ˙ x2(t)) ̸= 0 ∀ t and α(t) should be continuous It is better to enlarge a little the class of curves writing, instead than (1), the following condition ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∃ u(.) and v(.) ∈ L∞ ˙ x1(t) = u(t) cos(α(t)) ˙ x2(t) = u(t) sin(α(t)) ˙ α(t) =: v(t) α ∈ [−π/2, π/2]/ ∼ (2) Here (x1(t), x2(t), α(t)) should be of class Lipschitz.
  13. i.e. writing x = (x1, x2, α) if ˙ x

    = uX1 + vX2, X1 = ⎛ ⎝ cos(α) sin(α) 0 ⎞ ⎠ , X2 = ⎛ ⎝ 0 0 1 ⎞ ⎠ , in other words ˙ x ∈ (x) := Span{X1(x), X2(x)} even if in each point ˙ x belongs to a 2-D space, there are curves going everywhere in PT R2 since: [X1, X2](x) / ∈ (x) and dim(Span{X1, X2, [X1, X2](x)) = 3 completely non-integrable distribution ↔ H¨ ormander condition ⇓ (Chow theorem) for each pair of points there exists a trajectory joining them
  14. B1. How V1 reconstruct an interrupted curve 2 γ (

    ) γ ( ) b c 0 0 x x 1 Consider a smooth curve γ0 : [a, b] ∪ [c, d] → R2, interrupted in ]b, c[. We want to complete γ0 by a curve γ : [b, c] → R2 that is: γ(b) = γ0(b), γ(c) = γ0(c) ˙ γ(b) ∼ ˙ γ0(b) ̸= 0, ˙ γ(c) ∼ ˙ γ0(c) ̸= 0. we assume γ(b) ̸= γ(c), ˙ γ0(b) ̸= 0, ˙ γ0(c) ̸= 0
  15. B1. What to minimize? IDEA: Given an orientation column that

    is already active, it is easy to make activation of orientation columns that are: -) close to it, -) sensitive to a similar direction i.e. close in R2 × P1.
  16. The most natural cost for lift of planar curves on

    PTR2 Riemannian length: c b ˙ x2 1 + ˙ x2 2 + β2 ˙ α2 ds = c b u2 + β2v2 ds → min on all curves in PT R2 that are lift of planar curves (non-holonomic constraint). Then we get a problem of sub-Riemannian geometry (on PT R2): ˙ x = uX1 + vX2, c b u2 + β2v2ds → min, x = (x1, x2, α), X1 = ⎛ ⎝ cos(α) sin(α) 0 ⎞ ⎠ , X2 = ⎛ ⎝ 0 0 1 ⎞ ⎠ , initial and final positions are fixed in PT R2.
  17. digression: sub-Riemannian manifold Definition A sub-Riemannian manifold is a triple

    S = (M, , g), where (i) M is a connected smooth manifold of dimension n ≥ 3; (ii) is a smooth distribution of constant rank k < n, i.e. a smooth map that associates to q ∈ M a k-dimensional subspace q of Tq M satisfying the H¨ ormander condition. span{[X1, [. . . [Xj−1, Xj ]]](q) | Xi ∈ , j ∈ N} = Tq M, ∀ q ∈ M, (3) where denotes the set of horizontal smooth vector fields on M, i.e. = {X ∈ Vec(M) | X(q) ∈ q ∀ q ∈ M} . (iii) gq is a Riemannian metric on q , that is smooth as function of q.
  18. A Lipschitz continuous curve γ : [0, T ] →

    M is said to be horizontal if ˙ γ(t) ∈ γ(t) for a.e. t ∈ [0, T ]. Given an horizontal curve γ : [0, T ] → M, the length of γ is l(γ) = T 0 g(˙ γ, ˙ γ) dt. (4)
  19. horizontal curve (q) The distance induced by the sub-Riemannian structure

    on M is the function d(q0, q1) = inf{l(γ) | γ(0) = q0, γ(T ) = q1, γ horizontal}. (5) The hypothesis of connectedness of M and the H¨ ormander ⇒ →the finiteness and the continuity of d(·, ·) with respect to the topology of M (Chow’s Theorem) →The function d(·, ·) is called the Carnot-Caratheodory distance and gives to M the structure of metric space.
  20. Orthonormal frame Locally, the pair ( , g) can be

    given by assigning a set of k smooth vector fields (called a local orthonormal frame) spanning and that are orthonormal for g, i.e. q = span{X1(q), . . . , Xk (q)}, gq (Xi (q), Xj (q)) = δij . (6) The problem of finding the curve minimizing the length between two given points q0, q1, become the optimal control problem ˙ q = k i=1 ui Xi (q) T 0 k i=1 u2 i → min ↔ T 0 k i=1 u2 i → min ↔T → min, k i=1 u2 i ≤ 1 q(0) = q0, q(T ) = q1 you can think that a sub-Riemannian manifold is like a Riemannian manifold in which the inverse metric has some zero eigenvalues
  21. basic features in SRG t2 t Spheres are highly non-isotropic

    there are geodesics loosing optimality close to the starting point ⇒ spheres are never smooth even for small time (cut and conjugate loci are adiacent to the starting point) The Hausdorff dimension of (M, d) is always larger than its topological dimension (Mitchell theorem) geodesics are computed via Hamiltonian formalism (the Pontryagin Maximum Principle), but in general there are also “abnormal geodesics” (as it occurs in constrained optimization)
  22. Lot of open questions for instance regularity of minimizers: in

    (n − 1, n) they are smooth, but in larger co-dimension it is not known topology of small balls: it is not know if they are homeomorphic to euclidean balls see the book “A tour of sub-Riemannian geometry” by R. Montgomery the lecture notes “Introduction to Riemannian and sub-Riemannian geometry” Agrachev, Barilari, B.
  23. Remarks on this cost 1) The factor β can be

    eliminated with the transformation (x1, x2) → (βx1, βx2), i.e. by a “dilation of the initial conditions”. As a consequence, there is only one sub-Riemannian cost on PT R2 invariant by rototranslations of the plane, modulus dilations (observed by Agrachev) 2) c b u2 + β2v2ds → min ∼ c b u2 + β2v2 ↑ ↑ connec. among connnec. among hypercolumns orient.columns ds → min good model for the energy necessary to activate orientation columns which are not directly activated by the image
  24. 3) It is a compromise between length and curvature of

    the planar curve. Let γ = (x, y): c b ˙ x2 1 + ˙ x2 2 + β2 ˙ α2 ds = c b u2 + β2v2 ds = c b ∥˙ γ∥2 + β2∥˙ γ∥2K2 ds 4) there is existence of minimizers in the natural functional space D := {γ ∈ C2([b, c], R2) | ∥˙ γ(t)∥2 + β2∥˙ γ(t)∥2K2 γ (t) ∈ L1([b, c], R), γ(b) = x0, γ(c) = x1, ˙ γ(b) ≈ v0, ˙ γ(c) ≈ v1}. (7) →minimizers are analytic functions (there are no abnormal minimizers) 5) Since it is a sub-Riemannian cost, there is a natural hypoelliptic diffusion equation that can be used to reconstruct images
  25. Computation of optimal trajectories for curve-reconstruction step (a) compute candidate

    optimal trajectories with the Pontryagin Maximum Principle (they can be computed explicitly in terms of elliptic functions) step (b) evaluate their optimality (very difficult point) the local behavior of optimal trajectory is very complicated cut cut−conjugate conjugate →Yuri Sachkov and collaborators [2010-2013]
  26. This program has been essentially concluded by Y. Sachkov on

    SE(2) (double covering of PTR2) Cut Locus Part of the cut locus due to topological reasons seen as a full torus with no boundary Id R2 seen as an open disc S1 SE(2) ∼ R2 × S1
  27. Preliminary results of reconstruction of level sets by Yuri Sachkov

    and his group Original
  28. Preliminary results of reconstruction of level sets by Yuri Sachkov

  29. Preliminary results of reconstruction of level sets by Yuri Sachkov

  30. Complex images (not just a simple contour): ? ? ?

  31. An idea to reconstruct images all possible paths are activated

    as a Brownian motion For instance: dx = X1dW1 + X2dW2, → ∂t ψ(t, x) = (X1 2 + X2 2)ψ(t, x) X1 2 + X2 2 = (cos(α)∂x1 + sin(α)∂x2 )2 + β2∂2 α (sub-elliptic Heat equation, under H¨ ormander condition ⇒ solutions are smooth) The diffusion is highly non isotropic ∂t ψ(t, x) = (X1 2 + X2 2)ψ(t, x) −t log Pt (x, y) → − d(x, y)2 4t Varadhan-Leandre
  32. PLAN: Let I : R2 → [0, 1] be the

    image 0) smoothing the image with a Gaussian (it is made by the eyes) to get well defined level sets fc = I ∗ G(σx , σy ) ∈ L2(R2, R) ∩ C∞(R2, R), generically we get a Morse function. See: [1] Boscain, J. Duplaix, J.P. Gauthier, F. Rossi, SIAM SICON 2012 1) lifting the image to PT R2 and using it as an initial condition for the hypoelliptic heat eq. In the same paper we proved that the lift of a Morse function is a smooth surface in PT R2 2) computing the hypoelliptic diffusion ← 3) projecting down the image
  33. STEP 1) lifting an image to a distribution in PTR2

    Let us lift the Morse function fc in PT R2. This is made by associating with every point (x1, x2) of R2 the direction α ∈ R/ ∼ of the level set of fc at the point (x1, x2). L(fc ) = {(x1, x2, α) ∈ R2 c × P1 s.t. ∇fc (x1, x2) · (cos(α), sin(α)) = 0}, c level sets of f All directions are associated α c f
  34. Define the distribution on R × P1: ˆ fc (x1,

    x2, α) = fc (x1, x2)δ(d((x1, x2, α), L(f)) α 1 P smoothing Image lift of the image x x1 (distribution supported in L(f)) 2 corrupted part Ω
  35. STEP 2: hypoelliptic evolution There is no agreement on how

    to compute numerically the hypoelliptic diffusion. It is difficult due to two different scales finite difference schemes (Citti and Sarti group) numerical implementation of convolution kernels (Duits group) numerical integration of the Generalized Fourier Transform of the hypoelliptic diffusion equation on SE(2) that is a double covering of PT R2 The solution of the hypoelliptic diffusion equation can then be written as solutions of Mathieu type equations, ∂t Φ(t, θ) = (β2 d2 dθ2 + λ2 cos2(θ))Φ(t, θ) (8) [B, Duplaix, Gauthier, Rossi 2012]
  36. 3) Projecting down non isotropic diffusion 1 2 x x

    Image lift of the image α 1 P choice 1: fr (t, x1, x2) = P 1 ψ(x1, x2, α) dα choice 2: fr (t, x1, x2) = maxα∈P 1 {ψ(x1, x2, α)}←
  37. →This algorithm does not use the knowledge of where the

    image is corrupted (it is not an inpainting problem)
  38. Dynamic restoration Can we do more by using the information

    of where the image is corrupted to work on images as ? (85% of corrupted pixels)
  39. STEP0 we divide the pixels in “bad (corrupted)” and “good”

    (non corrupted) STEP1 we make a diffusion for δt using the previous method STEP2 each good point is “averaged” with the original point STEP3 bad points close to good points that got a sufficient mass are “upgraded” to good points. STEP4 repeat from STEP1
  40. None
  41. None
  42. None
  43. None
  44. Petitot observation: the visual cortex is probably doing “pure” hypoelliptic

    diffusion →V1 is able to reconstruct an image like: (pure hypoell. diffusion) →but not an image like (hyp. diff. + dyn. restor.) (pure hypoelliptic diffusion will not give good results on this immage)
  45. The END