Statistical Geometry Processing Winter Semester 2011/2012 - - PowerPoint PPT Presentation

(Deformable) Shape Matching

Statistical Geometry Processing

Winter Semester 2011/2012

Rigid Shape Matching

3

Iterated Closest Points (ICP)

The main idea:

• Pairwise matching technique
• Registers two scans
• Multi-part matching is a different story (more on this later)
• We want to minimize the distance between the two parts
• We set up a variational problem
• Minimize distance “energy” by rigid motion of one part

Part A

(stays fixed)

Part B

(moves, rotation & translation)

4

Iterated Closest Points (ICP)

Problem:

• How to compute the distance
• This is simple if we know the corresponding points.
• Of course, we have in general no idea of what corresponds...
• ICP-idea: set closest point as corresponding point
• Full algorithm:
• Compute closest point points
• Minimize distance to these closest points by a rigid motion
• Recompute new closest points and iterate

5

Closest Points

Distances: Closest points distances:

Part A

(stays fixed)

Part B

(moves, rotation & translation)

Part A

(stays fixed)

Part B

(moves, rotation & translation)

6

Iteration

Part A Part B Part A Part B Part A Part B final result

7

Variational Formulation

Variational Formulation:

 

 

    

   

n i B i nearest A i SO B SO

d A dist

1 2 ) ( ) ( ) ( ), 3 ( 2 ), 3 (

3 3

min arg ) , ( min arg p t Rp x t Rx

t R t R  

Variables: Orthogonal matrix R, translation vector t

8

Numerical Solution

Question: How to minimize this energy?

• The energy is quadratic
• There is only one problem...
• Constraint optimization
• We have to use an
• rthogonal matrix...
• This problem can (still) be solved exactly.

9

Solution

First step: computing the translation

• Easy to see: average translation is optimal

(c.f. total least squares)

• This is independent of the rotation

Second step: compute the rotation

• (2a) Compute optimal linear map
• (2b) Orthogonalize





 

n i B i nearest A i

n

1 ) ( ) ( ) (

1 p p t

10

Optimal Linear Map

First:

• Subtract translation from points pi

(A) = pi (A) – t

• Then: Solve an unconstrained least-squares problem
• Finally: compute the orthogonal matrix R that is

closest to M.

) ( ) ( ) (

~ : .. 1

B i nearest A i

n i p p M   

~

) ( ) ( ) ( 3 , 3 3 , 2 3 , 1 2 , 3 2 , 2 2 , 1 1 , 3 1 , 2 1 , 1

~ : .. 1

B i nearest A i

m m m m m m m m m n i p p             

unknowns (9 variables)

11

Least-Squares Optimal Rotation

How to compute a least-squares (Frobenius norm)

• rthogonal matrix that fits a general matrix:
• Compute the SVD: M = UDVT
• The least-squares orthogonal fit is: R = UVT

(just set all singular values to one)

• We can compute this in one step:
• Solve the least-squares matrix fitting problem using SVD
• Omit the diagonal matrix straight ahead

12

Generalizations

Convergence speed:

• Convergence of basic “point-to-point” ICP is not so great
• Typically: 20-50 iterations for simple examples
• Problem: Zero-th order method

(flip point correspondences in each step)

• Improvement: “point-to-plane” ICP
• First order approximation
• Match points to tangential planes rather than points
• Converges much faster (3-5 iterations for similar examples)

13

Implementation

Part A Part B

 

    

    

n i B i nearest B i nearest A i SO B SO

d A nearest A nearest

1 2 ) ( ) ( ) ( ) ( ) ( ), 3 ( 2 ), 3 (

, min arg )) ( ( ), ( min arg

3 3

n p t Rp x n t Rx

t R t R  

14

Implementation

Implementation:

• We need normals for each point (unoriented)  kNN+PCA
• Compute closest point, project distance vector to its

normal

• Minimize the sum of all such distances:

Part A Part B



  

 

n i B i nearest B i nearest A i SO 1 2 ) ( ) ( ) ( ) ( ) ( ), 3 (

, min arg

3

n p t Rp

t R 

15

Comparison

Point-to-point: 19 iterations Point-to-plane: 3 iterations

(accuracy problems) (much more accurate result)

16

Implementation

Problem:

• No closed form solution for the optimal rotation with

point-to-plane correspondences

Solution:

• Numerical solution
• Setup non-linear optimization problem (rotation,

translation = 6 parameters)

• Use non-linear optimization technique
• Remaining problem: Parametrization of the rotations
• Trouble with singularities (spherical topology)

17

Local Linearization

Standard technique: local linearization

• Transformation: T(x) = Rx + t
• Linearize rotations:

x I x x                                                                               1 ) ( ) cos( ) sin( ) sin( ) cos( 1 ) cos( ) sin( 1 ) sin( ) cos( 1 ) cos( ) sin( ) sin( ) cos(

, , , , , , , ,

y T x T y T T                      

           

18

Local Linearization

Standard technique: local linearization

• Numerical solution: iterative solver
• We have a current rotation R(i – 1) from the last iteration:
• Taylor expension at R(i – 1):
• Solve for t, , ,  (linear expressison  quadratic opt.)

x R I x

) 1 ( ) ( , ,

1 ) (



                        

i i

y T     

  



  

 

n j B j nearest B j nearest A j i 1 2 ) ( ) ( ) ( ) ( ) ( ) ( , ,

, min arg

3

n p t p R

t     

19

Local Linearization

Then:

• Project R(i) back on the manifold of orthogonal matrices.

(for example using the SVD-based algorithm discussed before)

• Then iterate, until convergence.

Why does this work?

• The parametrization is non-degenerate
• For large , , , the norm of the matrix increases arbitrarily

(i.e.: the object size increases, away from the data)

• Therefore, the least-squares optimization will perform a number
• f small steps rather than collapse.

20

More Tricks & Tweaks

ICP Problems:

• Partial matching might lead to distortions / bias
• Remove outliers (M-estimator, delete “far away points”, e.g.

20% percentile in point-to-point distance)

• Remove normal outliers

(if connection direction deviates from normal direction)

• Sampling problems
• Problem: for example flat surface with engraved letters
• No convergence in that case
• Improvement: Sample correspondence points with distribution

to cover unit sphere of normal directions as uniformly as possible

Deformable Shape Matching

22

Example

?

What are the Correspondences?

23

What are we looking for?

Problem Statement: Given:

• Two surfaces S1, S2  ℝ3

We are looking for:

• A reasonable deformation function f: S1  ℝ3

that brings S1 close to S2

?

S1 S2

f

24

Example

?

Correspondences? no shape match too much deformation

• ptimum

25

This is a Trade-Off

Deformable Shape Matching is a Trade-Off:

• We can match any two shapes

using a weird deformation field

• We need to trade-off:
• Shape matching (close to data)
• Regularity of the deformation field (reasonable match)

26

Components: Matching Distance: Deformation / rigidity:

Variational Model

27

Variational Model

Variational Problem:

• Formulate as an energy minimization problem:

) ( ) ( ) (

) ( ) (

f E f E f E

r regularize match

 

28

Part 1: Shape Matching

Assume:

• Objective Function:
• Example: least squares distance
• Other distance measures:

Hausdorf distance, Lp-distances, etc.

• L2 measure is frequently used (models Gaussian noise)

S2 f(S1)

 

2 1 2 , 1 ) (

), ( ) ( S S f dist f E match 







1 1

1 2 2 1 ) (

) , ( ) (

S x match

d S dist f E x x

29

Point Cloud Matching

Implementation example: Scan matching

• Given: S1, S2 as point clouds
• S1 = {s1

(1), …, sn (1)}

• S2 = {s1

(2), …, sm (2)}

• Energy function:
• How to measure ?
• Estimate distance to a point sampled surface

si

(2)

fi(S1)  







m i i match

S dist m S f E

1 2 ) 2 ( 1 1 ) (

, | | ) ( s

 

x ,

1

S dist

30

Surface approximation

Solution #1: Closest point matching

• “Point-to-point” energy

 







m i i S in i match

s NN s dist m S f E

1 2 ) 2 ( ) 2 ( 1 ) (

) ( , | | ) (

1

si

(2)

f(S1)

31

Surface approximation

Solution #2: Linear approximation

• “Point-to-plane” energy
• Fit plane to k-nearest neighbors
• k proportional to noise level, typically k  6…20

si

(2)

f(S1)

32

Variational Model

Variational Problem:

• Formulate as an energy minimization problem:

) ( ) ( ) (

) ( ) (

f E f E f E

r regularize match

 

33

What is a “nice” deformation field?

• Isometric “elastic” energies
• Extrinsic (“volumetric deformation”)
• Intrinsic (“as-isometric-as

possible embedding”)

• Thin shell model
• Preserves shape (metric plus curvature)
• Thin-plate splines
• Allowing strong deformations, but keep shape

Part II: Deformation Model

) (

) (

f E

r regularize

34

Elastic Volume Model

Extrinsic Volumetric “As-Rigid-As Possible”

• Embed source surface S1 in volume
• f should preserve 3  3 metric tensor (least squares)

 



   

1

2 T ) (

) (

V r regularize

dx f f f E I

first fundamental form I (ℝ3×3)

S1 V1 f S2 f f (V1)

ambient space

35

Volume Model

Variant: Thin-Plate-Splines

• Use regularizer that penalizes curved deformation

second derivative (ℝ3×3)

S1 V1 f S2 Hf =(f ) f (V1)

ambient space





1

2 ) (

) ( ) (

V f r regularize

dx x H f E

36

How does the deformation look like?

• riginal

as-rigid-as possible volume thin plate splines

37

Intrinsic Matching (2-Manifold)

• Target shape is given and complete
• Isometric embedding

 



   

1

2 T ) (

) (

S r regularize

dx f f f E I

first fund. form (S1, intrinsic)

Isometric Regularizer

S1 f S2 f

tangent space

38

“Thin Shell” Energy

• Differential geometry point of view
• Preserve first fundamental form I
• Preserve second fundamental form II
• …in a least least-squares sense
• Complicated to implement
• Usually approximated

Elastic “Thin Shell” Regularizer

S1 S2 f

I II I II

39

Deformable ICP

How to build a deformable ICP algorithm

• Pick a surface distance measure
• Pick an deformation model / regularizer

) ( ) ( ) (

) ( ) (

f E f E f E

r regularize match

 

40

Deformable ICP

How to build a deformable ICP algorithm

• Pick a surface distance measure
• Pick an deformation model / regularizer
• Initialize f (S1) with S1 (i.e., f = id)
• Pick a non-linear optimization algorithm
• Gradient decent (easy, but bad performance)
• Preconditioned conjugate gradients (better)
• Newton, Gauss Newton (even better, but more work)
• LBGFS (quick & effective)
• Always use analytical derivatives!
• Run optimization

Animation Reconstruction

42

Real-time Scanners

space-time stereo courtesy of James Davis, UC Santa Cruz color-coded structured light courtesy of Phil Fong, Stanford University motion compensated structured light courtesy of Sören König, TU Dresden

43

Animation Reconstruction

Problems

• Noisy data
• Incomplete data (acquisition holes)
• No correspondences

noise holes missing correspondences

44

Animation Reconstruction

Remove noise, outliers

Fill-in holes (from all frames) Dense correspondences

Urshape Factorization Approach

46

“Factorization”

t = 0 t = 1 t = 2

data urshape

S f f f

deformation

47

Components

Variational Model

• Given an initial estimate,

improve urshape and deformation

Numerical Discretization

• Shape
• Deformation

Domain Assembly

• Getting an initial estimate
• Urshape assembly

48

Components

Variational Model

• Given an initial estimate,

improve urshape and deformation

Numerical Discretization

• Shape
• Deformation

Domain Assembly

• Getting an initial estimate
• Urshape assembly

49

Energy Minimization

Energy Function

E(f, S) = Edata + Edeform + Esmooth

Components

• Edata(f, S)

– data fitting

• Edeform(f) – elastic deformation, smooth trajectory
• Esmooth(S)

– smooth surface

Optimize S, f alternatingly

50

Data Fitting

Data fitting

• Necessary: fi(S) ≈ Di
• Truncated squared distance

function (point-to-plane)

S Di fi Di fi(S)

Edata(f, S)

51

Edeform(f)

Elastic Deformation Energy

S Di f

Regularization

• Elastic energy
• Smooth trajectories

52

Surface Reconstruction

Data fitting

• Smooth surface
• Fitting to noisy data

S Di

Esmooth(S)

fi

• 1(Di)

S fi

• 1(Di)

S

f

53

Factorization

t = 0 t = 1 t = 2

data urshape

S f f f

deformation

54

Summary: Variational Model

) ( ) , )( ( ) , , ( ) , , (

urshape n deformatio data

S E S E E E E d S E d S E

smooth velocity accel volume rigid match

                               f f f



   

) ( 2

) , ( ) , ( ) ( ) , (

S V F T rigid rigid

dx t t x S E I x f x f f

x x



 



  

) ( 2

1 ) , ( ) ( ) , (

S V vol volume

dx t x S E x f f

x





          

S acc accel

dx t t x S E

2 2 2

) , ( ) ( ) , ( x f f 



        

S velocity velocity

dx t t x S E

2

) , ( ) ( ) , ( x f f 

 



 

S uv smooth smooth

dx x s x S E

2 2

) ( ) ( ) ( 



 



T t n i i match

t

S f d dist trunc d f S E

1 1 2)

)) ( , ( ( ) , , (

55

Components

Variational Model

• Given an initial estimate,

improve urshape and deformation

Numerical Discretization

• Shape
• Deformation

Domain Assembly

• Getting an initial estimate
• Urshape assembly

56

geometry

Discretization

Sampling:

• Full resolution geometry
• Subsample deformation

deformation

57

geometry

Discretization

Sampling:

• Full resolution geometry
• High frequency
• Subsample deformation
• Low frequency

deformation

58

geometry

Discretization

Sampling:

• Full resolution geometry
• High frequency, stored once
• Subsample deformation
• Low frequency, all frames ⇒ more costly

deformation

59

Shape Representation

Shape Representation:

• Graph of surfels (point + normal + local connectivity)
• Esmooth – neighboring planes should be similar
• Same as the bunny exercise...

S

60

...but how about Neighborhoods?

Topology estimation

• Domain of S, base shape (topology)
• Here, we assume this is easy to get
• In the following
• k-nearest neighborhood graph
• Typically: k = 6..20

Limitations

• This requires dense enough sampling
• Does not work for undersampled data

61

Volumetric Deformation Model

• Surfaces embedded in “stiff” volumes
• Easier to handle than “thin-shell models”
• General – works for non-manifold data

Deformation

geometry “thick shell”

f S V

62

Deformation

Deformation Energy

• Keep deformation gradients ∇f as-rigid-as-possible
• This means: ∇fT∇f = I
• Minimize: Edeform = ∫T ∫V||∇f(x,t)T∇f(x,t) – I||2 dx dt

geometry

f S

“thick shell”

∇f

V

63

Additional Terms

More Regularization

• Volume preservation: Evol = ∫T ∫V|

|det(∇f) – 1| |2

• Stability
• Acceleration:

Eacc = ∫T ∫V| |∂t

2 f|

|2

• Smooth trajectories
• Velocity (weak):

Evel = ∫T ∫V| |∂t f| |2

• Damping

64

Discretization

How to represent the deformation?

• Goal: efficiency
• Finite basis:

As few basis functions as possible

geometry deformation

65

Discretization

Meshless finite elements

• Partition of unity, smoothness
• Linear precision
• Adaptive sampling is easy

geometry deformation

66

Topology:

• Separate deformation

nodes for disconnected pieces

• Need to ensure
• Consistency
• Continuity
• Euclidean / intrinsic

distance-based coupling rule

• See references for details

Meshless Finite Elements

67

Adaptive Sampling

Adaptive Sampling

• Bending areas
• Decrease rigidity
• Decrease thickness
• Increase sampling density
• Detecting bending areas:

residuals over many frames

68

Components

Variational Model

• Given an initial estimate,

improve urshape and deformation

Numerical Discretization

• Deformation
• Shape

Domain Assembly

• Getting an initial estimate
• Urshape assembly

69

Urshape Assembly

Adjacent frames are similar

• Solve for frame pairs first
• Assemble urshape step-by-step

frame 11 frame 12 frame 13 frame 14 frame 15 frame 16

[data set courtesy of C. Theobald, MPC-VCC]

70

Hierarchical Merging

S f(S) data f

71

Hierarchical Merging

S f(S) data f

72

Initial Urshapes

S f(S) data f

73

Initial Urshapes

S f(S) data f

74

Alignment

S f(S) data f

75

Align & Optimize

S f(S) data f

76

Hierarchical Alignment

S f(S) data f

77

Hierarchical Alignment

S f(S) data f

Results

79

80

81

82