The Deformation Method for Grid Generation and Image Registration

Liao, Guojun Department of Mathematics University of Texas at Arlington

10. 2007 University of Houston

Outline of My Talk

1. The Deformation Method for Grid Generation and Numerical PDEs on Moving Grid 2. Image Registration as an Optimal Control Problem

Joint work with D. Anderson, P. Bochev, X. Cai, H. Chen, de la Pena, D. Fleitas, M. Gunzburger, B. Jiang, Z. Lei, F. Liu, J. Liu, S. Osher, T. Pan, B. Semper, J. Su, J. Xue

Grid Generation and Adaptation

To solve differential equations numerically by finite difference, finite element or finite volume methods, a good grid is needed for discretization of the physical domain.

Adaptation to physical variables;

Adaptation to geometric changes.

Approaches of Grid Generation

1. Cartesian/Body-Fitted 2. Structured/Unstructured/Multi-Blocks 3. Local Refinement/Moving Grid 4. Hybrid/Overlapping For the intended applications, we focus on a moving grid method, which redistributes a fixed number of nodes.

Invertible Transformations

• Existence and numerical construction of invertible transformations by differential equations is an interesting and challenging problem. • Harmonic maps approach works in two dimensions only (see [Rado 1926], [Melas 1993])

Harmonic maps are widely used to generate a body-fitted structured grid on a general domain.

In 2D, due to the Rado’s Theorem [Rado 1926], a harmonic map into a convex domain is injective if its restriction to the boundary is one-to-one and onto between the boundaries.

But the 3D extension of the Rado’s theorem is false [Melas 1993]. (The grid folding problem)

The Deformation Method (Version 3):

Construct a transformation T with prescribed

Jacobian determinant J(T) = f (T(x, t), t) > 0.

Use a div-curl-ode system to solve nonlinear problems:

Problem 1: Adaptive mesh generation with prescribed J. (Direct Problem)

Problem 2: Reconstruction of invertible transformation. (Inverse Problem)

Problem 3: Non-rigid image registration.

The Deformation Method (Version 3)

• Let O (t) be a family of moving domains. • Given a monitor (size) function f (x, y, z, t) > 0, we normalize it so that . (1/f) = |O (0)|. • Step 1: Determine a node velocity v satisfying div (v/f) = - (d/dt)(1/f). • Step 2: Find the new position T(x, y, z, t) of each node (x, y, z) at t = 0 from the ordinary differential equation: dT/dt = v (T(x, y, z, t), t) with T (x, y, z, 0) = (x, y, z).

Theorem: Steps 1 and 2 guarantee that J(T) = f (T(x, y, z, t), t)

Proof 1: Show directly that Steps 1 and 2 imply d/dt (J/f) = 0. Therefore J/f = constant = 1; and thus J = f.

Begin with the product rule (denoting h = 1/f) d/dt (J/f) = d/dt (J h) = h dJ/dt + J dh/dt.

By the Abel’s lemma (which says if dM/dt = AM, then d/dt det(M)= trace(A) det (M)), with M = grad h; A = grad v; and trace(A) = div v, we get dJ/dt = (div v) J.

We get

${\displaystyle {\begin{aligned}{\frac {d}{dt}}({\frac {J}{f}})&=h{\frac {dJ}{dt}}+J{\frac {dh}{dt}}\\&=hJ(\nabla \cdot {\textbf {v}})+J({\frac {\partial h}{\partial t}}+(\nabla h)\cdot {\frac {dT}{dt}})\\&=hJ(\nabla \cdot {\textbf {v}})+J({\frac {\partial h}{\partial t}}+(\nabla h)\cdot {\textbf {v}})\\\end{aligned}}}$

Proof 2 of J = f is based on the Transport Formula

The integrand is zero due to Step 1 for h = 1/f. Thus the first term of the following equations is zero. After performing change of variables, we get

Implementations of the Deformation Method

1. (For fixed Domains only) Solve a potential w from the Poisson equation . w = - ./.t (1/f), Neumann condition on

the boundary. Then set v/f = grad w. Check: Indeed we have div (v/f) = div grad w = .w = - ./.t (1/f).

For Domains with Fixed or Moving Boundary

2. Solve the div-curl system div (v/f) = - ./.t (1/f) curl (v/f) = g with the Neumann condition on fixed boundary and Dirichlet condition on moving boundary.

A Cube Deforms to a Cylinder

A Cube Deforms to a Cylinder Intermediate Time

A Cube Deforms to a Cylinder: Final Time

Stress Concentration Problem for a Plate with a Small Hole

A plate occupying the domain [-10,10]x[-10,10] \ a circle of radius 1, centered at (0,0) is subjected to a uniform load in the x-direction on its edges x = -10 and x = 10.

Only a quarter of the domain needs to be considered due to the symmetry of the problem.

Unstructured Mesh with 695 nodes, 640 linear quadrilateral elements X Y 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10

Plane Stress Equations

The first three are the stress-strain equations; the last two are the equilibrium conditions.

u and v are displacements in x and y direction respectively;

sx, sy and txy are stress components;

fx and fy are the body force components,

The material properties are entered in terms of the Young modulus and the Poisson ratio.

Least squares finite element method for the stress analysis problem

• We are looking for that minimizes the functional

where A is a first order differential operator, which is defined by the differential equations.

A is a first order differential operator determined by

The Matrices Are

A theoretical result of stress in the y- direction for infinite plate at (0,y) is

Using 695 nodes, 640 linear quadrilateral elements and reduced integration with one Gaussian point, the algebraic system is under-determined since

l= Nelem ´ NG ´ Neq - (Nnode ´ m - Nbc) = = 640 ´ 1´ 5 - (695 ´ 5 - 108) =- 167

The system will be extremly over -determined if two Gaussian points are used, since

l= N elem ´ N G ´ N eq - (N node ´ m - N bc ) = = 640´ 4´ 5 - (695 ´ 5 - 108) = 9433

Results with Two Gaussian Points are Poor

Sx 0 1 2 3 LSFEM solution 1 2 3 4 5 6 7 8 910

Y

Adapted Grid for Two Gaussian Points X Y 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10

Results with Two Gaussian Points on Adapted Grid are still not Good

1 2 3 4 5 6 7 8 9 10 1 2 3 exact LSFEM Y

Creation of Overlapping Elements

After adding 39 elements, we get a slightly over-determined system on the overlapped grid. With one Gaussian point, we have X Y 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 679 1 5 (695 5 108) 28

Results with one Gaussian point on the overlapped grid are good, but the max is still well below 3 due to insufficient resolution near (0,1).

Sx 0 1 2 3 LSFEM solution 1 2 3 4 5 6 7 8 910

Y

Adapted Overlapping Grid X Y 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10

Results with one Gaussian point on the adapted grid (refined near (0,1)) with overlapping. Now max = 3.07, which is in excellent agreement with the theoretical result: max = 3.00 (for infinite plate).

1 2 3 4 5 6 7 8 9 10 1 2 3 exact LSFEM Y

Propagating Shocks in Compressible Flows (Cylindrical Implosion)

Propagating Shocks

Propagating Shocks

Propagating Shocks

Shock-interactions of Hypersonic Flows • Left: Adaptive Grid; Right: Mach Contours

Turek’s Version (2) Ref[9]

• The moving mesh problem: constructing atransformation u, from the computational space to the physical space. The moving mesh method employedbelongs to the velocity-type class, which is based on Liao’s work and Moser’s work. • This method has several advantages: Only linear Poisson equation on fixed meshes are to be solved; monitor functions can be obtained directly from error distributions or distance functions; mesh tangling can be prevented; and the data structure

from the starting mesh is preserved.

Deformed Mesh Refined Near a Circle [9]

Moving Mesh Refined near 120 Particles [9]

Image Registration

• Image registration is the process of establishing acommon geometric reference frame between twoimages which could be of the same object (with different sensor pose and illumination) or different objects, could be from the same or different imagingmodalities (CT, MRI, PET, f MRI, etc), and possibly taken at different times. • The process usually consists of two components: (1) Rigid (or affine) registration through global rotation, translation, projection, and scaling. (2) Non-rigid registration, which accounts for local, nonlinear deformation.

Medical Image Registration

Template image Reference Image

Challenging Tasks

• Image registration is a key task in image fusion, segmentation, change detection, super-resolution imaging, robotic vision, and in building image information systems, among many others. • Although a lot of work has been done, automatic registration of images with complex nonlinear and local distortions, multimodal registration, and registration of N-D images (where N > 2) are still very challenging at this moment.

Regularizing Methods

• Optimize a similarity measure C-sim between the template (moving) and the reference (fixed) images. i.e., C-sim = the sum of squared differences (SSD). • Since the optimization problem for nonrigid registration is “ill posed”, a regularization term C-reg (penalty on the gradient or the Laplacian of displacements, or integrals based on continuum mechanics, etc) is added to C-sim. • The problem becomes: Optimize functional C=C-sim + k C-reg, where k is a parameter that represents the tradeoff between the similarity and regularity.

Issues with Regularity Terms

• If the weighing parameter k is too small, the algorithm will be unstable. If k is too large, the regularity will be too strong and the resulting transformations will not accurately optimize the similarity measure. • A major shortcoming of this approach is that the resulting transformation with any significant k distorts the similarity measure.

Optimal Control Approach

• Optimizing any chosen similarity measure subject to the constraint of a div-curl-ode system withoutregularizing term. • The method is an integration of two recent mathematical developments: 3.The deformation method for construction of differentiable and invertible transformations for mesh adaptation. The method is based on solving adiv-curl-ode system. 2. The theory and computational methods for optimalcontrol problems with partial differential equationsas constraints [13]. The right hand sides of the divcurl system are control variables.

Reconstruction of Transformation from a Div-Curl-Ode System

• The Direct Problem (Deformation Method Version 1): Given f and g, determine T(x, 1) such that J(T(x,1)) = f(x).

• Div u(x) = f(x) –1 • Curl u(x) = g(x) with u·n = 0, n = the outward normal on .D, • v (x, t) = u (x)/(t+(1– t)f(x)) for t in [0,1] • .T(x, t)/.t = v (T(x, t), t) with T(x, 0) = x.

The Inverse Problem

• Let S(x) be a given invertible transformation on D. Let f = (an approximation of) J(S). • Construct T(x, 1) by (an extension of [13]) Min . |T(x, 1) – S(x)|^2 dx over all g (with div g = 0) subject to • Div u(x)= f(x) –1 • Curl u(x) = g(x) with u ·n = 0, n = the outward normal on .D, • v (x, t) = u (x)/(t+(1– t)f(x)) for t in [0,1] • .T(x, t)/.t = v (T(x, t), t) with T(x, 0) = x.

Registration as an Optimal Control Problem

Similarity measure SSD = . (I1(x, y, z) – I2(f(x, y, z)))2 dD

minimize SSD = . (I1(x) – I2(T(x,1)))2 dD over all possible

controls f and g with f(x) > 0, . f dD = |D|, and div g = 0, subject to the constraints:

• Div u(x) = f(x) – 1, Curl u(x) = g(x) with u·n = 0, n = the outward normal on .D, • v(x, t) = u(x)/(t + (1– t) f(x)) for t in [0, 1] and x in D, • .T(x, t)/.t = v(T(x, t), t) for t in [0,1] and each fixed x in D with T(x, 0) = x.

Analysis of the Div-Curl-Ode Method

1. Existence, uniqueness, regularity, and injectivity of the solution to the direct problem; 2. How about the inverse problem? 3. The admissible space of the optimal control problem consists of all invertible, differentiable transformations; 4. Existence of the solution to the optimal control problem with div-curl-ode for image registration; 5. Co-state method and direct method.

Mesh Adaptation: 1/f is Proportional to the Gradient of a Mona Lisa image

Effect of Curl

• Left: Curl = 0 Right: Curl . 0

Local Stretching and Rotation x y 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Reconstruction by Div-Curl x y 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

3D Example 0 0.2 0.4 0.6 0.8 1 V3 0 0.2 0.4 0.6 0.8 1 V1 0 0.2 0.4 0.6 0.8 1 V2 X Y Z

Reconstruction (Max Error < 0.25 h, h =1/40) 0 0.2 0.4 0.6 0.8 1 Z 0 0.2 0.4 0.6 0.8 1 X 0 0.2 0.4 0.6 0.8 1 Y X Y Z

Constrained Optimization by the Lagrange Multiplier Method

Minimize an integral F(u) = . H(u(x)) dx subject to G(u, f, g) = 0, where u is in a function space V1, G is a mapping from V1 . V2, V2 is another function space.

Example: G(u, f, g) = (h, z) is defined by

div u – f = h and curl u – g = z. Reference:

1. Nonlinear Analysis and its Applications, vol. III, by E. Zeidler, Springer, 1988. 2. [13] by Gunzburger and Lee, 1999.

Existence of Lagrange Multipliers

A Model Problem by Div-Curl

I1= R(x, y) = 10 + 30 y2 I2= T(x, y) = 10 + 30y.

They agree on the top and bottom boundaries. But T is linear in y, while R is quadratic in y.

Exact solution: F(x, y) = (x, y2). The transformation pushes each horizontal mesh-line at y downward to y2.

Figure Deformation after one step Note: The new method produces an almost exact result in seconds. SSD is reduced from 29.5 to 5.2 in one step (SSD is reduced to 0.35 in two steps)

Two Ellipses by Div-Curl

Template image T Reference Image R

Deforming Ellipses

After 1 step After 10 steps

The Registration Error

Error after the first step Error after ten steps

The Aperture Problem

Template Image Reference Image Ground Truth

9 9 9 9 9 9 9 88 8 8 8 8 8 7 7 7 7 7 7 7 6 6 6 6 6 6 6 5 5 5 5 5 5 5

9 9 9 9 9 9 9 88 8 8 8 8 8 7 7 7 7 7 7 7 6 6 6 6 6 6 6 5 5 5 5 5 5 5

9 9 9 9 9 9 9 8 8 8 8 8 8 8 7 7 7 7 6 6 6 6 6 6 6 5 5 5 5 5 5 5 7 7 7 6

5

4

3

2

1

1 2 3 4 5 6 7 812345678 6

5

4

The Aperture Problem

3

2

1

1234 567 8

Fast Parametric Method The div-curl-ode Method

Template Image

Reference Image

Deformed (from Template) Image

SSD / Iteration

Summary

• The new approach is based on “solid” mathematical foundation. It accounts for local volume changes through the divergence of the transformation; and accounts for local rotation through the curl vector. • It is based on a linear partial differential system, and thus its numerical implementation is well understood. • The method may be used in other optimization problems that involve motion or displacement estimation.

References

1. J. Moser, Volume elements of a Riemannian manifold, Trans. AMS, 120, 1965 2. A. D. Melas, An Example of a Harmonic Map Between Euclidean Balls, Proceedings of the American Mathematical Society, Vol. 117, No. 3. 1993 3. T-W. Pan, J. Su and G. Liao, A numerical grid generator based on Moser's deformation method, Numerical Methods for Partial Differential Equations, Vol. 10, p. 20-31, 1994 4. F. Liu, S. Ji, and G. Liao, An adaptive grid method and its application to steady Euler flow calculations, SIAM J. Sci. Comput. 20, 811-825, 1998 5. G. Liao, F. Liu, C. de la Pena, D. Pang, and S. Osher, Level set based deformation methods for adaptive grids, Journal of Computational Physics, vol. 159, 2000

References (Continues)

6. Z. Lei and G. Liao, Adaptive grids for resolution enhancement, Shock Waves, vol. 12, p.153-156, 2002 7. X. Han, C. Xu, and J. L. Prince, A 2D Moving Grid Geometric Deformable Model, IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2003 8. D. Fleitas, X. Cai, G. Liao, B. Jiang, The Least-Squares Finite Element Method on Overlapping Elements, Journal of Computational Information Systems 1:2, 2005 12. S. Turek and D. Wan, Fictitious Boundary and Moving Mesh Methods for the Numerical Simulation of Particulate Flow, Journal of Computational Physics 222 (2007) 28–56 10. G. Liao et. al. Optimal Control Approach to Data Set Alignment, to appear in Applied Math. Letters

References (Continues)

11. Analysis and computation of adaptive moving grids by deformation, P. Bochev, G. de la Pena, and G. Liao, Numerical Methods for Partial Differential Equations, Vol. 12, pp. 489-506, 1996

12. A note on harmonic maps, H. Liu and G. Liao, Appl. Math. Lett., Vol. 9, No. 4,1996 13. Gunzburger and Lee, Analysis and approximation of optimal control problems for first-order elliptic systems in three dimensions, Appl. Math. Comp., vol. 100, 1999

Related Books

1. Numerical Grid Generation, by J. Thompson et. al. 1985 2. The Fundamentals of Grid Generation, by Knupp, P. and S. Steinberg, 1993 3. Computational Grid Generation, Adaptation, and Solution Strategies, by G. Carey, 1997 4. Least Squares Finite Element Method, by B. Jiang, Springer 1998 5. Grid Generation Methods, by V. A, Leseikin, Springer 1999