Shape representation as distribution of geometric features: from - - PowerPoint PPT Presentation

1/46

Workshop on Statistical modeling for shapes and imaging

Shape representation as distribution of geometric features: from currents to oriented varifolds

Irène Kaltenmark

Inria - Institut de Mathématiques de Bordeaux Université de Bordeaux

joint work with Nicolas Charon and Benjamin Charlier March 2019

2/46

Outline

Shape representation as distribution of geometric features Context of computational anatomy Distributions of points Current and varifold representations Functional shapes Conclusion Growth model for computational anatomy

3/46

Shape Spaces [Grenander, Miller, Trouvé, Younes, Beg, ..., Arguillère]

Define a metric on a space of shapes

d(S0, S1) . = inf{dG(Id, g) | g · S0 = S1} . Consider G . = {ψv

1 | v

∈ L2([0, 1], V )} where V is a space of vector fields on a fixed ambient space, e.g. Rn, and ψv

t = Id +

t

vs ◦ ψv

s ds .

G is equipped with a distance: dG(Id, ψ) . = inf{

1

|vt|2

V dt | ψv 1 = ψ} .

G acts on S(Rn) : ψ · S = ψ(S) and induces an infinitesimal action of V on S(Rn) : ξ(S) : V → TSS(Rn), v → v · S = v(S) . Strength: this approach handles a large type of shape data and combines definition of distances with shape registration

4/46

From Shape Spaces to inexact matching

How to model the shapes ?

Approximate

Inexact matching problem: E(v) = 1 2

1

|vt|2

V dt + λ

2 D(ψ(S0), S1)2 Data attachment term: currents and varifolds model shapes as distributions of their local tangent or normal vectors [Glaunès 05, Charon 13]. Reproducing kernels of Hilbert spaces lead to various types of data fidelity metrics between shapes that do not depend on the parametrization. Framework presented here: the choice of the metric reduces to one or two scalar functions. Related work: normal cycles [Roussillon 16], optimal transport as data term [Feydy 17], square-root velocity (SRV) representation [Srivastava 09, 12]

5/46

Back to simple functions

How to evaluate f ?

Properties of the test functions : ◮ Scale (e.g. for Gaussian distributions) ◮ Shape of the tail

6/46

Shapes as distributions

• X

ϕ(||x − x0||)dvolX(x) ≈ mass of X around x0 ,

• X

ϕ(||x − y||)dvolX(x) = 0 .

7/46

Shapes duality

Functions Shapes RKHS approach If ∈ D′ ϕ ∈ D f : R → R µX ∈ D′ ω ∈ D X (submanifold) µX ∈ H′ ϕX ∈ H X (submanifold) If (ϕ) =

R f ϕ

µX(ω) =

X ω(x)dvolX(x)

General kernel: X ← → ϕX =

• X

ϕx Degenerate kernel (Dirac) : X ← → 1 1X =

• X

δx

8/46

Shapes duality

RKHS approach If ∈ D′ ϕ ∈ D f : R → R µX ∈ D′ ω ∈ D X (submanifold) µX ∈ H′ ϕX ∈ H X (submanifold) If (ϕ) =

R f ϕ

µX(ω) =

X ω(x)dvolX(x)

For two shapes X, Y , we have µX, µY H′ = ϕX, ϕY H . For two points x0, y0 ∈ Rn δx0, δy0H′ = ϕ(||y0 − x0||Rn) .

Note that here ϕ : R → R (equivariance to rigid motion), e.g. ϕ(u) = exp(−u2/σ2).

9/46

Scalar product and norm

µX, µY H′ =

• Y
• X

ϕ(||y − x||)dvolX(x)dvolY (y) = µY (ϕX) . µY acts linearly on ϕx0 : it integrates the mass of Y in a neigh- borhood of x0 µY (ϕx0) =

• Y

ϕ(||y − x0||)dvolY (y) . X and Y are similar around x0 if µY (ϕx0) ≈ µX(ϕx0) . . Finally, d(X, Y )2 . = ||µX − µY ||2

H′ = ||µX||2 − 2µX, µY + ||µY ||2 .

10/46

Discretization

µXi , µYj =

• Yj
• Xi

ϕ(||y − x||)dvolX(x)dvolY (y) ≈ mX

i

• Yj

ϕ(||y − cX

i ||)dvolY (y)

≈ mX

i mY j ϕ(||cY j

− cX

i ||)

µX, µY ≈

• i,j

mX

i mY j ϕ(||cY j

− cX

i ||)

µX0, µX ≈ 1 × ϕ(0) µX0, µX′ ≈

1

2 + 1 2

• × ϕ(0)

If the diameters of the faces are small wrt to the scale of ϕ (how flat is ϕ around 0), the associ- ated metric is robust to parametrization.

11/46

Discretization

◮ This metric induces registration algorithms that do not require pointwise correspondences between shapes. For example, one can register a curve with 100 points on a curve with 150 points. ◮ Shapes union: if X = X1 ∪ X2, then µX = µX1 + µX2 − µX1∩X2 . ◮ This metric is quite robust to discontinuities. Conversely, normal cycles allow to discriminate topological properties [Roussillon 16].

12/46

Mass sensitivity

X Y

µY ≈ 2µX so that d(X, Y )2 = ||µX − µY ||2 ≈ ||µX||2

Application to brain registration

Inclusion detection ? ”d(X, Y )2 ≈ ||µY2||2” X and Y have the same length but ”||µX|| > ||µY ||”

This is not without consequences

13/46

Data compression

[Gori 17] Cross section of 3 fiber bundles Consider a fiber bundle B = n

i Xi such that for

any i, Xi is close to an average fiber ˜ X, i.e. µXi ≈ µ˜

X ,

then µB =

n

• i

µXi ≈ nµ˜

X .

14/46

Outline

Shape representation as distribution of geometric features Context of computational anatomy Distributions of points Current and varifold representations Functional shapes Conclusion Growth model for computational anatomy

15/46

Limitations

The tangent space distributions of shapes highlight natural correspondences... as well as the shapes’ orientation:

16/46

Tangent bundle

For curves and surfaces in R2

• r R3,
• riented tangent spaces can be represented

by a unique unit vector (tangent or normal). µX is now a distribution on a space of test functions ω : Rn × Sn−1 → R , (n = 2, 3) : µX(ω) =

• X

ω(x, − → tx )dvolX(x) δ

− → tx x , δ − → ty y H′ = ϕ(||y −x||Rn)γ(−

→ tx , − → ty Rn) .

Note that ϕ : R → R, γ : R → R (equivariance to rigid motion).

17/46

Representation choices

X, Y two shapes, µX, µY the respective distributions of their local tangent or normal vectors µX, µY H′ =

• X×Y

ϕ(||y − x||Rn)γ(− → tx , − → ty Rn)dvolX(x)dvolY (y) ≈

• i,j

mX

i mY j ϕ(||cY j

− cX

i ||)γ(

− → tX

i ,

− → tY

j )

Data type γ(− → tx , − → ty ) γ(u) Currents [Glaunès 05] − → tx , − → ty u (unoriented) Varifolds [Charon 13] − → tx , − → ty 2 u2 Oriented Varifolds* exp 2− → tx , − → ty /σ2

T

• exp

2u/σ2

T

• ∝ exp

−|− → tx − − → ty |2/σ2

T

• Remark: σT does not depend on the scale of the shapes.

*[Kaltenmark,Charlier,Charon 17]

18/46

Representation choices

X, Y two shapes, µX, µY the respective distributions of their local tangent or normal vectors µX, µY H′ =

• X×Y

ϕ(||y − x||Rn)γ(− → tx , − → ty Rn)dvolX(x)dvolY (y) ≈

• i,j

mX

i mY j ϕ(||cY j

− cX

i ||)γ(

− → tX

i ,

− → tY

j )

Data type γ(− → tx , −− → tx ) Currents [Glaunès 05] −γ(− → tx , − → tx ) = −1 (unoriented) Varifolds [Charon 13] γ(− → tx , − → tx ) = 1 Oriented Varifolds* γ(− → tx , −− → tx ) = 0 *[Kaltenmark,Charlier,Charon 17]

19/46

Current regularization

δ

− → tx x , µY H′ = 2

• i=1

mY

i ϕ(||yi − x||)−

→ tx , − → tY

i

≈ ϕ(0)− → tx , mY

1

− → tY

1 + mY 2

− → tY

2

≈ ˜ mϕ(||˜ y − x||)− → tx , − → t˜

y

≈ δ

− → tx x , ˜

mδ

− → t˜

y

˜ y H′

Mass cancellation effect: if mY

1

− → tY

1 + mY 2

− → tY

2 ≈ 0, then ||µY ||H′ ≈ 0 .

Counterpart: thin and sharp structures are not seen by the metric or at best their mass is underestimated, e.g. tails, tip of horns, leaf stem, etc.

20/46

Mass conservation with varifolds

Theorem (Charon 13)

If ϕ and γ are two continuous positive functions, ϕ(0) > 0, and γ(1) > 0, then there exists c > 0, such that for any k-dimensional compact submanifold X cvolk(X) ≤ ||µX||H′ .

21/46

Applications in the diffeomorphometric framework

Inexact matching problem: deform a shape S0 towards a target shape S1 E(v) = 1 2

1

|vt|2

V dt + λ

2 |µψ1(S0) − µS1|2

H′

This energy is minimized by a gradient descent. Each step consists in: ◮ Computing the gradient of the data attachment term ◮ Pulling it backward over time ◮ Updating t → vt ◮ Updating ψ1(S0) (end point of the geodesic generated by v)

22/46

Denoising with currents

Source: blue. Target: red.

23/46

Denoising with currents

Source: green. Target: red.

24/46

Comparisons

Source: blue. Target: red. Currents Unoriented varifolds Oriented varifolds

25/46

Unoriented varifolds and multi-attachment terms

Source: blue. Target: red.

26/46

Restoration with oriented varifolds

Source: blue. Target: red. Initial shape is an ellipsoid

27/46

Functional shapes

Consider two shapes X, Y ⊂ Rn , respectively equipped with a signal f X : X → Rd, f Y : Y → Rd, then one can define µ(X,f X ),µ(Y ,f Y )H′ =

• X×Y

ϕ(||y − x||Rn)γ(− → tx , − → ty Sn−1)ϕf (||f Y (y) − f X(x)||Rd )dvolX(x)dvolY (y) [Charon 14] Clustering fibers under the constraint: take into account the ending points of the fibers [Gori 17] Examples of signal: thickness/pressure measurements, brain activity map, etc.

28/46

Functional shapes

◮ Computing the previous metric for functional shape is straighforward ◮ Plug your new favorite metric in a deformation framework:

• Do the deformations act on the feature ? (yes: dΦ(x) · −

→ tx , no: labels from segmentation)

• Do you need an independant action on the feature/signal ? (metamorphosis).

For registration allowing signal deformation, see Fshapes Toolkit [Charon, Charlier 15,18] Matching of 2 membranes in the retina with thickness measurements.

29/46

Mass sensitivity

Inclusion detection ? ”d(X, Y )2 ≈ ||µY2||2” X and Y have the same length but ”||µX|| > ||µY ||”

This is not without consequences

29/46

Mass sensitivity

Inclusion detection ? ”d(X, Y )2 ≈ ||µY2||2” X and Y have the same length but ”||µX|| > ||µY ||”

This is not without consequences

Interlude

Jean Feydy

CMLA, ENS de Cachan

30/46

Conclusion - Why two metrics ?

On Growth and Form, D’Arcy Thompson, 1917.

31/46

Conclusion - What global model for computational anatomy ?

Assumption: the geometry gives enough prior to retrieve the biological homologies between shapes Kernel metrics+LDDMM OT+LDDMM [Feydy et al.] Low variability Local changes Medium variability High variability Complex phenomena Priors on deformations [Gris] Elastic model for thickness estimation [Younes] Growth model with priors on the biological process [Kaltenmark] Etc. False True

32/46

Conclusion

Currents and varifolds metrics can be used in many applications outside the diffeomorphometry framework ! For example, to replace the Dice coefficient c = 2|X ∩ Y | |X| + |Y | . Top: red and blue have same area. Bottom: blue is twice larger than red. Code: Fshapes Toolkit (Matlab) [Charlier, Charon]. You can customize very simply the functions ϕ and γ ! The next talk will introduce the KeOps library [Charlier, Feydy, Glaunès]!

33/46

Expedition in a non Diffeomorphic World

Shape representation as distribution of geometric features Context of computational anatomy Distributions of points Current and varifold representations Functional shapes Conclusion Growth model for computational anatomy

34/46

Growth and homology issues

Longitudinal data

1 2 3 4 2 3 4 5 2 3 4 5 1 6 1 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

Type I Type I & II

(a) Plant growth (b) Bone growth On the macroscopic scale, one can identify two main types of growth process: ◮ Type I: a growth homogeneously distributed. ◮ Type II: a growth process that involves new material on specific areas (e.g. plant growth, crystal growth or mineralized tissues as bone, horn, mollusc shells, tendon, cartilage, tooth enamel).

35/46

Individual evolution model

Creation processes call for inner partial mappings. Def: a growth scenario is a curve of shapes, equipped with a flow of diffeomorphisms that ◮ define partial relations of homology between the ages of the shapes, ◮ allow the creation of new coordinates.

36/46

Population model

t

Subject 2 Template Biological Coordinate System Subject 1

Shape spaces can be generalized for these growth scenarios. The metric results from the action of a group of spatio-temporal deformations. New structuring hypothesis: a shape space of growth scenarios models a population

• f shape curves that share a common creation process.

This creation process is encoded by a biological coordinate system.

37/46

Time-dependent infinitesimal action

The biological coordinate system is the set of all the coordinates required to parametrize an individual growth scenario. These coordinates have a birth tag that allows to anticipate their emergence in the ambient space and to explicit the gradual embedding of the ages of the shape through the flow. The action of a vector field is filtered by the birth tag ∀t ∈ [0, 1], ∀x ∈ X, ˙ qt(x) = 1 1τ(x)≤t vt(qt(x)) . filter standard infinitesimal action

38/46

Reconstitution of a scenario

We consider a population of scenarios, all aligned on the time interval [0, 1]. This population admits a common biological coordinate system (X, τ). ◮ The final shape Star of an individual’s scenario is given. ◮ The initial position q0 : X → Rd is known. Aim: We want to retrieve the complete scenario of this individual.

?

39/46

LDDMM: Gradient descent on the vector field

Inexact matching problem

Our aim is to find the simplest vector field vopt ∈ L2

V such that q1(X) approximates

the target shape Star as well as possible. General minimization problem: E : B × L2

V → R

E(q0, v) = 1 2

1

|vt|2

V dt + A(q1)

where ˙ qt = ξ(qt,t)(vt) . Any minimizer must satisfy vopt

t

= K V J (qt, pt, t) where J : B × B∗ × [0, 1] − → V ∗ (q, p, t) − → ξ∗

(q,t) · p

and p1 = −dA(q1) ∈ B∗ , ˙ pt = −∂qξ(qt,t)(vt)∗ · pt.

40/46

Construction of a diffeomorphism with a RKHS

Example of deformation when the space of vector fields V is a Reproducing Kernel Hilbert Space with a Gaussian kernel. The vector field below is generated by 3 control points (δpi

xi )i=1,2,3, pairs of positions and directions.

⇒ Large Deformation Diffeomorphic Metric Mapping (LDDMM) methods

41/46

Application

Superimposed trajectories of the solution and the target: The continuous trajectory is completely determined by the initial position q0 and the initial momenta p0 (right fig- ure).

42/46

The end

Merci pour votre attention !

A penny for your thoughts: What generic term for currents, varifolds, oriented varifolds, and normal cycles ?

43/46

Momentum and Momentum Map with the growth dynamic

With the growth dynamic, the momentum map builds the sum of the control points δp(x)

q(x) indexed by the active coordinates. The time variable t specifies the actual

support Xt ⊂ X. J (qt, pt, t) =

• x∈Xt

δpt(x)

qt(x) ,

vopt

t

= K V J (qt, pt, t) =

• x∈Xt

kV (qt(x), · )pt(x) . Examples: vopt

t

=

• x∈Xt

e− |qt (x)− · |2

2σ2

pt(x)

• r

vopt

t

=

• x∈Xt

pt(x) .

44/46

Variants of the minimization problem

For a large class of solutions (q, p) ∈ C([0, 1], B × B∗), t → |vopt

t

|V is strictly

• increasing. However, when the deformation are generated by global rotations and

translations, we would like this norm to be rather constant. Flop.png ◮ Adapted norm: E(q0, v) = 1 2

1

αt |vt|2

V dt + A(q1) .

The minimizers of this energy satisfy vopt

t

= 1 αt K V J (qt, pt, t). ◮ Constrained norm: consider the classic energy E(q0, v) = 1 2

1

|vt|2

V dt + A(q1) ,

and minimize it under the constraint that for any t ∈ [0, 1], |vt|2

V = ct

where c : [0, 1] → R+ is known. (augmented Lagrangian method)

45/46

Non constant growth

Comparison of the norms of the vector fields between the target (red) and the solution (blue):

Norm of At*(blue) 0.5 1 0.5 1 1.5 2 2.5 3 3.5 Norm of Nt*(blue) 0.5 1 0.5 1 1.5 2 2.5 3 3.5 Norm of At*(blue) 0.5 1 0.5 1 1.5 2 2.5 3 3.5 Norm of Nt*(blue) 0.5 1 0.5 1 1.5 2 2.5 3 3.5

Adapted norm Constrained norm

46/46

Hamiltonian framework

For any minimizer vopt ∈ L2

V , the trajectory (q, p) ∈ C([0, 1], B × B∗) is solution of

the system

  

˙ qt = ∂Hr

∂p (qt, pt, t)

˙ pt = − ∂Hr

∂q (qt, pt, t)

where Hr(qt, pt, t) = 1 2 |K V J (qt, pt, t)|2

V .

The novelty here is that the Hamiltonian function depends on the time. Yet, it can be controlled by the initial condition (q0, p0). It leads to a new minimization problem: ˆ E(q0, p0) = E(q0, v) where vt = K V J (qt, pt, t) , and to a new algorithm (so-called shooting method on the initial momentum).