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
Wiesel (1959), Hoﬀman (’89), Petitot (’99), Citti and Sarti (2003), our research group. I will mainly present the ﬁnal version of the model studied in [1] U. Boscain, J. Duplaix, J.P. Gauthier, F. Rossi, Anthropomorphic image reconstruction via hypoelliptic diﬀusion. SIAM J. CONTROL OPTIM.Vol. 50, No. 3, pp. 1309-1336, 2012. [2] U. Boscain, R. Chertovskih, J.P. Gauthier, A. Remizov. Hypoelliptic diﬀusion 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 diﬀusion (in the sense of H¨ ormander)
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
visual cortex V1 has a structure of sub-Riemannian manifold reconstructing level sets via geodesics reconstruction of complex images: the hypoelliptic diﬀusion model new ideas to reconstruct deeply corrupted images (dynamic restoration)
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)
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)
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 ﬁber bundle whose base is R2 and whose ﬁber at the point (x1, x2) is the set of straight lines (i.e. directions without orientation) d(x1,x2) attached to (x1, x2).
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.
= 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
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.
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 ﬁnal positions are ﬁxed in PT R2.
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 ﬁelds 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.
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)
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 ﬁniteness 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.
given by assigning a set of k smooth vector ﬁelds (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 ﬁnding 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
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 Hausdorﬀ 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)
(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.
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
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 diﬀusion equation that can be used to reconstruct images
optimal trajectories with the Pontryagin Maximum Principle (they can be computed explicitly in terms of elliptic functions) step (b) evaluate their optimality (very diﬃcult point) the local behavior of optimal trajectory is very complicated cut cut−conjugate conjugate →Yuri Sachkov and collaborators [2010-2013]
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
image 0) smoothing the image with a Gaussian (it is made by the eyes) to get well deﬁned 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 diﬀusion ← 3) projecting down the image
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
to compute numerically the hypoelliptic diﬀusion. It is diﬃcult due to two diﬀerent scales ﬁnite diﬀerence schemes (Citti and Sarti group) numerical implementation of convolution kernels (Duits group) numerical integration of the Generalized Fourier Transform of the hypoelliptic diﬀusion equation on SE(2) that is a double covering of PT R2 The solution of the hypoelliptic diﬀusion 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]
(non corrupted) STEP1 we make a diﬀusion 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 suﬃcient mass are “upgraded” to good points. STEP4 repeat from STEP1
diﬀusion →V1 is able to reconstruct an image like: (pure hypoell. diﬀusion) →but not an image like (hyp. diﬀ. + dyn. restor.) (pure hypoelliptic diﬀusion will not give good results on this immage)