Shape Correspondence and Functional Maps SGP 2017 course Maks - - PowerPoint PPT Presentation

Shape Correspondence and Functional Maps

SGP 2017 course

Maks Ovsjanikov

Based on joint work with: E. Corman, A. Chambolle, M. Ben-Chen, O. Azencot, A. Butscher, R. Rustamov, J. Solomon, F. Chazal, L. Guibas, D. Nogneng, ….

Laboratoire d’Informatique de l’École polytechnique

General Overview

Overall Objective: Create tools for computing and analyzing mappings between geometric objects.

General Overview

Rather than comparing points on objects it is often easier to compare real-valued functions defined on them. Such maps can be represented as matrices. Overall Objective: Create tools for computing and analyzing mappings between geometric objects.

Course Overview

Course Notes: [Related] Course Website:

http://www.lix.polytechnique.fr/~maks/fmaps_course/

Linked from the website. Or use:

Sample Code:

See Sample Code link on the website.

• r http://bit.do/fmaps

http://bit.do/fmaps_notes

Course Overview

Motivation and Problem Taxonomy Rigid Matching: ICP Functional Map representation, properties Open problems, Q&A Basic pipeline for non-rigid matching Extensions, Improvements

What is a Shape?

Discrete: a graph embedded in 3D (triangle mesh). Continuous: a surface embedded in 3D.

• Connected.
• Manifold.
• Without Boundary.

Common assumptions:

What is a Shape?

5k – 200k triangles

Shapes from the FAUST, SCAPE, and TOSCA datasets

Discrete: a graph embedded in 3D (triangle mesh). Continuous: a surface embedded in 3D.

Overall Goal

Given two shapes, find correspondences between them.

Overall Goal

Given two shapes, find correspondences between them. Finding the best map between a pair of shapes.

Problem Taxonomy

Local vs. Global

refinement (e.g. ICP) | alignment (search) . .

Rigid vs. Deformable

rotation, translation | general deformation.

Semi vs. Fully Automatic

given landmarks | a priori model

Learning-Based vs. Direct

known examples | unseen data

Problem Taxonomy

Local vs. Global

refinement (e.g. ICP) | alignment (search) . .

Rigid vs. Deformable

rotation, translation | general deformation.

Semi vs. Fully Automatic

given landmarks | a priori model

Learning-Based vs. Direct

known examples | unseen data

Why Shape Matching?

Given a correspondence, we can transfer:

texture and parametrization segmentation and labels deformation

Other applications: shape interpolation, reconstruction ...

Sumner et al. ‘04.

12

Rigid Shape Matching

• The unknowns are the rotation/translation

parameters of the source onto the target shape.

• Given a pair of shapes, find the optimal Rigid

Alignment between them.

Iterative Closest Point (ICP)

• Classical approach: iterate between finding

correspondences and finding the transformation:

example in 2D

M

N

Given a pair of shapes, and , iterate:

• 1. For each

find nearest neighbor .

• 2. Find optimal transformation

minimizing:

arg min

R,t

X

i

kRxi + t yik2

2

M

N

xi ∈ M yi ∈ N

Iterative Closest Point

• Classical approach: iterate between finding

correspondences and finding the transformation:

M

N

Given a pair of shapes, and , iterate:

• 1. For each

find nearest neighbor .

• 2. Find optimal transformation

minimizing:

arg min

R,t

X

i

kRxi + t yik2

2

M

N

xi ∈ M yi ∈ N

Iterative Closest Point

• Classical approach: iterate between finding

correspondences and finding the transformation:

M

N

Given a pair of shapes, and , iterate:

• 1. For each

find nearest neighbor .

• 2. Find optimal transformation

minimizing:

arg min

R,t

X

i

kRxi + t yik2

2

M

N

xi ∈ M yi ∈ N

Iterative Closest Point

• Classical approach: iterate between finding

correspondences and finding the transformation:

M

N

Given a pair of shapes, and , iterate:

• 1. For each

find nearest neighbor .

• 2. Find optimal transformation

minimizing:

arg min

R,t

X

i

kRxi + t yik2

2

M

N

xi ∈ M yi ∈ N

Iterative Closest Point

• Classical approach: iterate between finding

correspondences and finding the transformation:

M

N

Given a pair of shapes, and , iterate:

• 1. For each

find nearest neighbor .

• 2. Find optimal transformation

minimizing:

arg min

R,t

X

i

kRxi + t yik2

2

M

N

xi ∈ M yi ∈ N

Iterative Closest Point

• Classical approach: iterate between finding

correspondences and finding the transformation:

M

N

Given a pair of shapes, and , iterate:

• 1. For each

find nearest neighbor .

• 2. Find optimal transformation

minimizing:

arg min

R,t

X

i

kRxi + t yik2

2

M

N

xi ∈ M yi ∈ N

Iterative Closest Point

• Classical approach: iterate between finding

correspondences and finding the transformation:

M

N

Given a pair of shapes, and , iterate:

• 1. For each

find nearest neighbor .

• 2. Find optimal transformation

minimizing:

arg min

R,t

X

i

kRxi + t yik2

2

M

N

xi ∈ M yi ∈ N

• 1. Finding nearest neighbors: can be done with space-

partitioning data structures (e.g., KD-tree).

• 2. Finding the optimal transformation

minimizing:

Iterative Closest Point

• Classical approach: iterate between finding

correspondences and finding the transformation:

Can be done efficiently via SVD decomposition.

arg min

R,t

X

i

kRxi + t yik2

2

M

N

Arun et al., Least- Squares Fitting of Two 3-D Point Sets

ICP: Optimal Transformation

Problem Formulation: 1. Given two sets points: in . Find the rigid transform: that minimizes: 2. Closed form solution: 1. Construct: , where 2. Compute the SVD of C: 1. If 2. Else 3. Set Note that C is a 3x3 matrix. SVD is very fast.

{xi}, {yi}, i = 1..n

R, t

N

X

i=1

kRxi + t yik2

2

Arun et al., Least-Squares Fitting of Two 3-D Point Sets

C = UΣV T

det(UV T ) = 1, Ropt = UV T

Ropt = U ˜ ΣV T , ˜ Σ = diag(1, 1, . . . , −1)

C = PN

i=1(yi − µY )(xi − µX)T

µX = 1

N

P

i xi, µ

topt = µY − RoptµX

, µY = 1

N

P

i yi,

Non-Rigid Shape Matching

Unlike rigid matching with rotation/translation, there is no compact representation to optimize for in non-rigid matching.

Non-Rigid Shape Matching

• What does it mean for a correspondence to be “good”?
• How to compute it efficiently in practice?

Main Problems:

Isometric Shape Matching

Good maps must preserve geodesic distances. Possible Model:

Geodesic: length of shortest path lying entirely on the surface.

dM(x, y)

dN (T(x), T(y))

M

N

Isometric Shape Matching

Approach:

Find the point mapping by minimizing the distance distortion: The unknowns are point correspondences.

Topt = arg min

T

X

x,y

kdM(x, y) dN (T(x), T(y))k

dM(x, y)

dN (T(x), T(y))

M

N

Isometric Shape Matching

Approach:

The space of possible solutions is highly non-linear, non-convex.

Problem:

Find the point mapping by minimizing the distance distortion:

Topt = arg min

T

X

x,y

kdM(x, y) dN (T(x), T(y))k

dM(x, y)

dN (T(x), T(y))

M

N

Functional Map Representation

We would like to define a representation of shape maps that is more amenable to direct optimization.

1. A compact representation for “natural” maps. 2. Inherently global and multi-scale. 3. Handles uncertainty and ambiguity gracefully. 4. Allows efficient manipulations (averaging, composition). 5. Leads to simple (linear) optimization problems.

Background: Laplace-Beltrami Operator

Given a compact Riemannian manifold without boundary, the Laplace-Beltrami operator :

∆ : C∞(M) ! C∞(M), ∆f = div rf

∆

M

Heat Equation on a Surface

Given a compact surface without boundary the evolution of heat is given by: ∂f

∂t = ∆f = div rf.

Laplace-Beltrami Operator

Given a compact surface without boundary, the Laplace-Beltrami operator :

• 1. Is invariant under isometric deformations.
• 2. Has a countable eigendecomposition:

that forms an orthonormal basis for .

• 3. Characterizes the geodesic distances fully.

∆

∆φi = λiφi

L2(M) M

The Laplace-Beltrami operator has an eigendecomposition:

Laplace-Beltrami Eigenfunctions

∆

∆φi = λiφi

Ordered from low frequency (smoothest) to higher frequency (oscillating).

λ0 = 0 λ1 = 2.6 λ2 = 3.4 λ3 = 5.1

φ0

φ1

φ2

φ3

. . .

λi = Z

M

krφi(x)k2dµ(x)

Any (square-integrable) can be represented as a linear combination of the LB eigenfunctions.

Laplace-Beltrami Eigenfunctions

= a0 a1 + . . . +

f =

∞

X

i=0

aiφi

f : M → R

ai = Z

M

f(x)φi(x)dµ(x)

φ0 φ1

…. that forms an orthonormal basis for :

L2(M)

In the Discrete World

• Functions are defined at vertices of the mesh.
• Integration is defined with respect to a discrete

volume measure:

• diagonal matrix of area weights.
• Laplacian is discretized as a matrix

kfk2

2 = f T Af

A

L = A−1W

i j

αij βij t1 t2

Lij = 1 2A(j) (cot(αij) + cot(βij))

Can be derived from 1st order FEM.

In the Discrete World

• Computing the eigenfunctions of the Laplacian reduces

to solving the generalized eigenvalue problem:

• eigs function in Matlab
• Both A and W are sparse positive semidefinite.

Lφ = λφ ⇔ Wφ = λAφ

Number of triangles Computation time (in s) 5000 0.65 25000 2.32 50868 3.6 105032 10

Time to compute 100 basis functions.

Functional Approach to Mappings

The map induces a functional correspondence: TF (f) = g, where g = f T M

N

T T

36

Functional maps: a flexible representation of maps between shapes, O., Ben-Chen, Solomon, Butscher, Guibas, SIGGRAPH 2012

Given two shapes and a pointwise map T : N → M

Functional Approach to Mappings

f : M → R

TF

TF (f) = g : N → R

The map induces a functional correspondence: T

37

TF (f) = g, where g = f T

Functional maps: a flexible representation of maps between shapes, O., Ben-Chen, Solomon, Butscher, Guibas, SIGGRAPH 2012

Given two shapes and a pointwise map T : N → M

Functional Approach to Mappings

f : M → R

TF

TF (f) = g : N → R

The map induces a functional correspondence: T

38

TF (f) = g, where g = f T

Functional maps: a flexible representation of maps between shapes, O., Ben-Chen, Solomon, Butscher, Guibas, SIGGRAPH 2012

Given two shapes and a pointwise map T : N → M

Functional Approach to Mappings

The induced functional correspondence is linear:

f : M → R

TF

TF (f) = g : N → R

TF (α1f1 + α2f2) = α1TF (f1) + α2TF (f2)

39

Given two shapes and a pointwise map T : N → M

Functional Map Representation

The induced functional correspondence is complete.

f : M → R

TF

TF (f) = g : N → R

40

Given two shapes and a pointwise map T : N → M

Observation

Express both and in terms of basis functions:

f TF (f)

Since is linear, there is a linear transformation from to .

TF

{ai}

{bj}

M

f : M → R

g : N → R

TF

N

f = X

i

aiφM

i

g = TF (f) = X

i

biφN

i

Functional Map Definition

Functional map: matrix C that translates coefficients from to .

ΦM

ΦN

Functional Maps

Definition:

Cij : coefficient of

in the basis of .

TF (φM

j )

φN

i

Cij = Z

N

TF (φM

j )φN i dµ

In an orthonormal basis: For a fixed choice of basis functions , and a linear transformation between functions, a functional map is a matrix C, s.t. for any if , then:

{φM}, {φN }

TF

f = P

i aiφM i

T(f) = P

i biφN i

b = Ca

Example 1

Given two shapes with points and a map: If functions are represented as vectors (in the hat basis), the functional map is given by matrix-vector product:

matrix encoding the map T,

• ne 1 per column with zeros everywhere else.

g = TT f C = TT

M

N

T :

nM, nN

T : nN × nM

T : N → M

Example 2

If functions are represented in the reduced basis:

matrix of the first eigenfunctions of as columns. matrix of the first eigenfunctions of as columns.

The functional map matrix:

if if

ΦM : nM × kM

ΦN : nN × kN

kM

kN

∆N

ΦT

N ΦN = Id

C = Φ+

N TT ΦM

C = ΦT

N AN TT ΦM

ΦT

N AN ΦN = Id

∆M

C = ΦT

N TT ΦM

Given two shapes with points and a map:

matrix encoding the map T,

• ne 1 per column with zeros everywhere else.

nM, nN

T : nN × nM

T : N → M

: left pseudo-inverse.

+

Example Maps in a Reduced Basis

Triangle meshes with pre-computed ponitwise maps “Good” maps are close to being diagonal

Reconstructing from LB basis

Map reconstruction error using a fixed size matrix.

47

0.5 1 1.5 2 2.5 3 3.5 4 4.5

reconstruction error

Number of basis (eigen)-functions 27.9k vertices

Functional Map algebra

• 1. Map composition becomes matrix multiplication.
• 2. Map inversion is matrix inversion (in fact, transpose).
• 3. Algebraic operations on functional maps are possible.

E.g. interpolating between two maps with

C = αC1 +(1−α)C2.

Course Overview

Motivation and Problem Taxonomy Rigid Matching: ICP Functional Map Representation, properties Open Problems, Q&A Basic pipeline for non-rigid matching Extensions, Improvements

In practice we do not know C. Given two objects our goal is to find the correspondence. How can the functional representation help to compute the map in practice?

Shape Matching

?

Matching via Function Preservation

where Given enough pairs, we can recover C through a linear least squares system.

f = P

i aiφM i ,

g = P

i biφN i .

{a, b}

Suppose we don’t know C. However, we expect a pair of functions and to correspond. Then, C must be s.t.

Ca ≈ b

f : M → R

g : N → R

Function preservation constraint is general and includes:

• Attribute (e.g., color) preservation.
• Descriptor preservation (e.g. Gauss curvature).
• Landmark correspondences (e.g. distance to the point).
• Part correspondences (e.g. indicator function).

Map Constraints

Suppose we don’t know C. However, we expect a pair of functions and to correspond. Then, C must be s.t.

Ca ≈ b

f : M → R

g : N → R

Commutativity Constraints

Regularizations: Commutativity with other operators:

C

Note that the energy: is quadratic in C.

SM SN

CSM = SN C

kCSM SN Ck2

F

Regularization

Lemma 1:

The mapping is isometric, if and only if the functional map matrix commutes with the Laplacian:

Implies that exact isometries result in diagonal functional maps.

C∆M = ∆N C

Functional maps: a flexible representation of maps between shapes, O., Ben-Chen, Solomon, Butscher, Guibas, SIGGRAPH 2012

Lemma 2:

The mapping is locally volume preserving, if and only if the functional map matrix is orthonormal:

Map-Based Exploration of Intrinsic Shape Differences and Variability, Rustamov et al., SIGGRAPH 2013

CT C = Id

Regularization

Lemma 3:

If the mapping is conformal if and only if:

CT ∆1C = ∆2

Regularization

Map-Based Exploration of Intrinsic Shape Differences and Variability, Rustamov et al., SIGGRAPH 2013

Basic Pipeline

Given a pair of shapes :

• 1. Compute the first k (~80-100) eigenfunctions of the Laplace-

Beltrami operator. Store them in matrices:

• 2. Compute descriptor functions (e.g., Wave Kernel Signature)
• n . Express them in , as columns of :
• 3. Solve
• 4. Convert the functional map

to a point to point map T.

Copt =

diagonal matrices of eigenvalues

• f LB operator

M, N

ΦM, ΦN

ΦM, ΦN

∆M, ∆N :

Copt = arg min

C

kCA Bk2 + kC∆M ∆N Ck2

A, B

M, N

Conversion to point-to-point

Given a functional map C, we would like to convert to to a point-to-point map. Option 1: declare Problems: high computational complexity , low accuracy.

f : M → R

g : N → R

TF

T(x) = arg max

y

ΦN Cδx O(nMnN )

Conversion to point-to-point

Given a functional map C, we would like to convert to to a point-to-point map. Option 2: declare Advantages: computational complexity , higher accuracy (e.g., works with the identity map).

T(x) = arg min

y

kδy Cδxk2

f : M → R

g : N → R

TF

O(nM log nN )

Incorporating Orthonormality

In many practical situations we would expect a volume- preserving map, which implies:

CT C = Id

Option: use post-processing to enforce this constraint. Iterate: 1. Compute the point-to-point map T. 2. Solve for the functional map: Exactly the same objective as ICP, but in higher dimension. Can use the same method!

arg min

C, s.t. CT C=Id

X

x∈M

kCδx δT (x)k2

2

Results

A very simple method that puts together many constraints and uses 100 basis functions gives reasonable results:

61 Functional maps: a flexible representation of maps between shapes, O., Ben-Chen, Solomon, Butscher, Guibas, SIGGRAPH 2012

Results

radius 0.025 radius 0.05

Functional maps: a flexible representation of maps between shapes, O., Ben-Chen, Solomon, Butscher, Guibas, SIGGRAPH 2012

A very simple method that puts together many constraints and uses 100 basis functions gives reasonable results:

62

Segmentation Transfer without P2P

To transfer functions we do not need to convert functional to pointwise maps. E.g. we can also transfer segmentations: for each segment, transfer its indicator function, and for each point pick the segment that gave the highest value.

Course Overview

Motivation and Problem Taxonomy Rigid Matching: ICP Functional Map Representation, properties Open Problems, Q&A Basic pipeline for non-rigid matching Extensions, Improvements

Basic Pipeline Extensions

Maps in Collections Manifold Optimization Promoting Pointwise Maps

Informative Descriptor Preservation via Commutativity for Shape Matching, Nogneng, O., Eurographics 2017 Functional map networks for analyzing and exploring large shape collections Huang, Wang, Guibas, SIGGRAPH 2014 Image Co-Segmentation via Consistent Functional Maps Huang, Wang, Guibas, CVPR 2013 MADMM: A generic algorithm for non-smooth optimization on manifolds. Kovnatsky, Glashoff, M. Bronstein, ECCV, 2016.

Basic Pipeline

Key Optimization Step: Can solve in (at least) two ways:

Copt = arg min

C

kCA Bk2

F + kC∆M ∆N Ck2 F

• 1. Rewrite in vector form:

Use standard LS solver. e.g. in Matlab.

copt = arg minc kAc bk2

F

c = A\b

• 2. Use a generic non-linear solver. For the gradient use:

∂tr(CA) ∂C

= AT ,

∂tr(CT A) ∂C

= A

advantages: less memory, often faster, more flexible.

Basic Pipeline

For area-preserving maps:

Copt = arg min kCA Bk2

F + kC∆M ∆N Ck2 F

• No-longer a convex problem. However, can use

manifold optimization techniques.

CT C = Id

Basic Pipeline

For area-preserving maps:

Copt = arg min kCA Bk2

F + kC∆M ∆N Ck2 F

• Manifold of orthonormal matrices:

Stiefel manifold.

• Can use manopt package. Only

needs extrinsic gradient.

• Can be extended to handle non-smooth
• bjective functions.

CT C = Id

CT C = Id

MADMM: A generic algorithm for non-smooth optimization on manifolds. Kovnatsky, Glashoff, M. Bronstein, ECCV, 2016.

Basic Pipeline

Using a robust norm:

MADMM: A generic algorithm for non-smooth optimization on manifolds. Kovnatsky, Glashoff, M. Bronstein, ECCV, 2016.

Copt = arg min

CT C=Id

kCA Bk2,1 + kC∆M ∆N Ck2

F

Correspondence computed with data containing 10% outliers

Least Squares Robust MADMM

Consistent Maps in Shape Collections

Given a collection of shapes, we expect the maps to satisfy the loop closure property.

Functional map networks for analyzing and exploring large shape collections Huang, Wang, Guibas, SIGGRAPH 2014 Image Co-Segmentation via Consistent Functional Maps Huang, Wang, Guibas, CVPR 2013

…

Consistent Maps in Shape Collections

Given a collection of shapes, we expect the maps to satisfy the loop closure property. In functional language we expect: difficult (non convex) constraint. Y

ij∈cycle

Cij = Id

Consistency via Latent Space Optimization

Find the optimal latent space by solving: If this reduces to an eigenvalue problem. Given ‘s, solve for to enforce consistency. Restart.

Cij = YjY −1

i

Given a collection of maps, enforce loop closure by creating a “latent” shape:

min

Y kCijYi Yjk2 F

YT Y = Id

Yi

Cij

Consistency via Latent Space Optimization

Attention: In step 1, can help to make Given a collection of maps , and the latent functions

Cij

E(C, Y) = kCijYi Yjk2

F Find the optimal C, Y by alternating:

• 1. minYT Y=Id E(C, Y)
• 2. minC E(C, Y)

eigenvalue problem least squares Yi ∈ Rk×k2, k2 < k Attention: Can add other terms in step 2, e.g., kCij Cold

ij k2 F

Consistency via Latent Space Optimization

Application to Co-segmentation:

Functional map networks for analyzing and exploring large shape collections, Huang, Wang, Guibas, SIGGRAPH 2014

Consistency via Latent Space Optimization

Application to Co-segmentation:

Image Co-Segmentation via Consistent Functional Maps Wang, Huang, Guibas, CVPR 2013

Basic Pipeline Extensions

Maps in Collections Manifold Optimization Promoting Pointwise Maps

Informative Descriptor Preservation via Commutativity for Shape Matching, Nogneng, O., Eurographics 2017 Functional map networks for analyzing and exploring large shape collections Huang, Wang, Guibas, SIGGRAPH 2014 Image Co-Segmentation via Consistent Functional Maps Huang, Wang, Guibas, CVPR 2013 MADMM: A generic algorithm for non-smooth optimization on manifolds. Kovnatsky, Glashoff, M. Bronstein, ECCV, 2016.

Making Functional Maps Point-to-Point

(Known) Theoretical result: A functional map is point-to-point iff it preserves pointwise products of functions:

C(fh) = C(f)C(h) ∀ f, h (fh)(x) = f(x)h(x)

Question: When does a linear functional mapping correspond to a pull-back by a point-to-point map?

• J. von Neumann, Zur operatoren methode in der klassichen Mechanik, Ann. of Math.(2) 33 (1932)

Making Functional Maps Point-to-Point

(Known) Theoretical result: A functional map is point-to-point iff it preserves pointwise products of functions. Approach Represent descriptor functions via their action on functions through multiplication. C

Informative Descriptor Preservation via Commutativity for Shape Matching, Nogneng, O., Eurographics 2017

Sfk Sgk

CSfk = SgkC

Extended Basic Pipeline

Given a pair of shapes :

• 1. Compute the multi-scale bases for functions on the two
• shapes. Store them in matrices:
• 2. Compute descriptor functions (e.g., Gauss curvature) on

. Express them in , as columns of :

• 3. Solve
• 4. Convert the functional map

to a point to point map T.

Copt =

M, N

ΦM, ΦN

ΦM, ΦN Copt = arg min

C

kCA Bk2 + kC∆M ∆N Ck2

A, B

M, N

Informative Descriptor Preservation via Commutativity for Shape Matching, Nogneng, O., Eurographics 2017

+ X

k

kCSfk SgkCk2

Results with extended basic pipeline

Incorporating multiplicative operators improves results significantly.

Informative Descriptor Preservation via Commutativity for Shape Matching, Nogneng, O., Eurographics 2017

before after before after

Results with extended basic pipeline

The problem is well-constrained even for larger number of basis functions.

Informative Descriptor Preservation via Commutativity for Shape Matching, Nogneng, O., Eurographics 2017

Some Open Problems

What is the optimal choice of basis? How to guarantee a continuous pointwise map? What are better deformation models? How much can we do without converting to p2p?

Thank you! Questions?