Physically-Based Modeling for Medical Image Registration Huai-Ping - - PowerPoint PPT Presentation

Physically-Based Modeling for Medical Image Registration

Huai-Ping Lee COMP 768 Presentation

What is Image Registration?

A moving image Im(x) and a fixed image If(x), can be 2D or 3D Rigid registration: find an affine transformation such that is closest to If(x)

Find a matrix to do the job, relatively easy

Non-rigid registration: find a deformation field to do the job

F F(Im) F

2

Medical Image Registration

Human body is always deforming

3 Fluid-flow-based deformable registration Rigid translation based on bony landmarks

http://titan.radonc.unc.edu/~foskey/prost_web/warp_movies/

Medical Image Registration

Note that we want to match not only the contour but also internal voxels

Correspondence of voxels between daily scanned images and the first image is important for calculating total radiation dose received by each part of the patient

a) b) c)

Optical Flow FEM Registration Displacement field

4

Outline

Classical Methods for Deformable Registration Viscous Fluid Model Liner Elasticity Model

• Finite Element Method
• Mesh Generation
• Boundary Condition (External Force)

5

Classical Methods

Image based methods

Key point matching + interpolation [Shen02] Optical flow [Thirion98] Energy minimization [Foskey05] [Lu04]

Model based methods

Mesh adaption + interpolation [Kaus07]

6

Viscous Fluid Model [Christensen 96]

The entire human body is modeled as compressible fluid, governed by PDE u: displacement field v: velocity field, given by Convective derivative: rate of change of material at point x, given by chain rule µ∇2v(x) + (µ + λ)∇(∇· v(x)) = f(x, u(x)) Diffusion Change of flux external force v = ∂u ∂t + (v· ∇) u du dt = ∂u ∂t + ∂u ∂x ∂x ∂t = ∂u ∂t +

d

• i=1

∂ui ∂xi vi

7

Viscous Fluid Model [Christensen 96]

External force is computed with a similarity measure By taking the derivative, they defined force at each dΩ Boundary condition is v = 0 The PDE is solved with finite difference method f(x, u(x, t)) = −α(IT (x − u(x, t)) − IS(x))∇T

• x−u(x,t)

8

C(IT (x), ISu(x, t)) = α 2

• Ω

|IT (x − u(x, t)) − IS(x)|2 dΩ

Viscous Fluid Model [Christensen 96]

Moving Image Target Image Deformed Image 9 Displacement-x Displacement-y

Viscous Fluid Model

Advantage: can deal with large deformation Disadvantage: does not consider different properties of different organs Human body is not totally fluid... Finite difference method does not fit non-regular domain very well

10

Linear Elasticity Model

Model organs as elastic objects Hooke’s Law says For multi-dimensional problems, it becomes

σ: stress, measure of force per unit area within a body ε: strain, deformation of materials caused by stress

D: elasticity matrix

σ = Dε F = −kx

11

Linear Elasticity

σ =         σx σy σz τxy τyz τxz         σ = Dε ε =         

∂u ∂x ∂u ∂y ∂u ∂z ∂u ∂x + ∂u ∂y ∂u ∂y + ∂u ∂z ∂u ∂x + ∂u ∂z

         =         

∂u ∂x ∂u ∂y ∂u ∂z ∂u ∂x ∂u ∂y ∂u ∂y ∂u ∂z ∂u ∂x ∂u ∂z

         = Lu

12

Strain and Stress

From Wikipedia

13

Linear Elasticity

Potential energy of elastic body Minimize E by driving its derivative to zero, and solve for u with finite element method

E = 1 2

• Ω

σT εdΩ +

• Ω

FudΩ E = 1 2

• Ω

uT LT DLudΩ +

• Ω

FudΩ

14

Finite Element Method

The domain is subdivided into finite number of elements that consist of nodes For each element, approximate function u with shape (basis) functions Ni and node values ui

u =

Nnodes

• i=1

Niui

15

FEM−Basis Function

One-dimensional, linear basis functions

From Wikipedia

FEM−Basis Function

Two-dimensional, linear basis functions

a) b)

17

Finite Element Method

Potential energy for a node in an element To minimize E, take

18

E

• uel

1 , . . . , uel Nnodes

• =

1 2

• Ω

Nnodes

• i=1

Nnodes

• j=1

uel

i T Bel i T DBel j uel j dΩ

+

• Ω

Nnodes

• i=1

FN el

i uel i dΩ

x Bel

i = LiN el i

n (see Eq. 2.1)

• Ω

Nnodes

• j=1

Bel

i T DBel j uel j dΩ = −

• Ω

FN el

i

dΩ ; i = 1, . . . , Nnodes

∂E

• uel

1 , . . . , uel Nnodes

• ∂uel

i

= 0 ; i = 1, . . . , Nnodes

FEM−Matrix Form

Build the linear system for each element Use numerical integration methods such as Gaussian Quadratures Combine all local systems into one global system Solve the system with iterative solvers such as conjugate gradient

19

   Kel

i,j =

• Ω Bel

i T DBel j dΩ

Fel

i =

• Ω FN el

i

dΩ

Keluel = −Fel

(Stiffness Matrix)

FEM−Remaining Issues

How to decide the material properties?

The matrix D

How to create nodes and elements?

Tetrahedralization of the image (organs)

How to decide the right-hand-side?

External force term F

20

Material Properties

E: Young’s modulus ν: Poisson ratio Determined empirically!

D = E(1 − ν) (1 + ν)(1 − 2ν)          1

ν (1−ν) ν (1−ν) ν (1−ν)

1

ν (1−ν) ν (1−ν) ν (1−ν)

1

1−2ν 2(1−ν) 1−2ν 2(1−ν) 1−2ν 2(1−ν)

        

21

Mesh Generation

Delaunay Tetrahedralization

Needs to generate polygonal mesh in advance (e.g., with marching cubes [Brock05]) Provides well-shaped elements (near optimal aspect ratio)

Directly subdivide image domain into tetrahedra

Each voxel is labeled to be inside an organ or not [Ferrant99], [Ferrant01]

22

Mesh Generation [Brock05]

Run marching cubes algorithm on down- sampled labeled image Run Laplacian smoothing several times

23

Mesh Generation [Ferrant99]

Divide the domain into grids, and subdivide each grid into six tetrahedra Each tetrahedron is subdivided according to the image labels

24

Mesh Generation [Ferrant99]

When subdividing tetrahedra, prisms may appear Subdivide them further

Left Right Back Case 1 Case 4 Case 2 Case 3 Case 5 Case 8 Case 6 Case 7

25

Mesh Generation [Ferrant99]

a) b) c)

26

Mesh Generation [Ferrant01]

Subdivide tetrahedra when necessary Multi-resolution scheme

1 tet + 2 pyramids 2 tets + 2 pyramids 1 tet + 1 prism 2 prisms 1 tet + 1 pyramid 2 tets + 1 + 1 prism pyramid 4 tets + 2 pyramids 2 tets 4 tets 4 tets

27

Mesh Generation [Ferrant99, 01]

Advantage: simple Disadvantage: quality of tetrahedra is not as good as quality of those generated by Delaunay tetrahedralization

a) b) c) d) e) f)

28

Boundary Condition

Ku = F, where F is external force on each node Volumetric matching [Ferrant99] Surface matching [Ferrant00] (So what happened after 2000? nothing special...)

29

Volumetric Matching [Ferrant99]

Instead of computing F separately, use a local image similarity constraint in the energy function And use first order approximation

30

E =

• Ω

σT εdΩ +

• Ω

(I1(x + u(x)) − I2(x))2dΩ I1(x + u(x)) ≈ I1(x) + ∇I1(x) · u(x)

Volumetric Matching [Ferrant99]

After taking the derivatives, we get four equations for each element:

(K + G)u = F

a) b) c)

Optical Flow FEM Registration Displacement field of FEM 31

• Ω

4

• j=1
• Bel

i T DBel j + N el i ∇IT 1 ∇I1N el j

• uel

j dΩ =

• Ω
• I1 − I2
• ∇I1N el

i

dΩ (5.6

Results [Ferrant99]

a) b) c) a) b) c)

32

Surface Matching [Ferrant00]

Deform the boundary surface of an object in moving image toward the boundary in fixed image For boundary elements, change the linear system to F is a decreasing function of the gradient, minimized at the edges of the fixed image Solve for locations of boundary nodes Same method can be used for predefined landmarks, but F has to be defined differently

33

vt − vt−1 τ + Kvt = −Fvt−1

Surface Matching [Ferrant00]

Define a distance map of the edges in the target image: Intuitively, external force should push the edge in source image toward the edge in target image

re D

• I(x)
• int x. S

F(x) = SminGexp∇

• D
• I(x)

Surface Matching [Ferrant00]

Smin: make F point towards a point that is closer to the edge Gexp: avoid sticking on a wrong edge

Smin =    +1 if D

• I(x)
• > D
• I
• x + ∇
• D
• I(x)
• −1 if D
• I(x)
• < D
• I
• x + ∇
• D
• I(x)
• Gexp =

+1 if kn∇I(x) > 0 −1 if kn∇I(x) < 0

n: normal of surface k=1 if the region to be matched has higher intensity than surroundings; k= -1 otherwise

Surface Matching [Ferrant00]

a) b) c) d) e) f) a) b) f)

External Force [Ferrant00]

Compute external force from moved surface Let be the vector representing the displacement on the boundary nodes F can be computed with Solve Ku = -F for u at all nodes Interpolate u for each voxel

˜ u F = K˜ u

37

Results [Ferrant00]

a) b) c) Surface of moving image

• verlaid on fixed image

Surface of moving image Deformed surface overlaid

• n fixed image

a) b) c) c)

Moving Image Fixed Image Deformed Image Difference of b & c 38

Results [Ferrant00]

a) b) c)

Deformed image Original Moving image Fixed image 39

Linear Elasticity Model

Advantage:Elastic model is more close to human body we can assign different material properties for different organs Finite element method models non-regular domains much better than finite difference method Disadvantage: linear model only works for small deformation! Need nonlinear model for large deformation And what about interaction between organs? Body fluid and gas between organs?

40

Reference-Physically Based Methods

[Christensen 96] Gary Christensen, Richard D. Rabbitt, and Michael I. Miller. Deformable templates using large deformation kinematics. IEEE Trans. Image Process., 5(10), October 1996. [Ferrant99] Matthieu Ferrant, Simon K. Warfield, Charles R. G. Guttmann, Robert

• V. Mulkern, Ferenc A. Jolesz, and Ron Kikinis. 3d image matching using a finite

element based elastic deformation model. In MICCAI, pages 202–209, 1999. [Ferrant00] Matthieu Ferrant, Simon K. Warfield, Arya Nabavi, Ferenc A. Jolesz, and Ron Kikinis. Registration of 3d intraoperative MR images of the brain using a finite element biomechanical model. In MICCAI, pages 19–28, 2000. [Ferrant01] Matthieu Ferrant. Physics-based deformable modeling of volumes and surfaces for medical image registration, segmentation and visualization. PhD Thesis, 2001. [Brock05] K. K. Brock, M. B. Sharpe, L. A. Dawson, S. M. Kim, and D. A.

• Jaffray. Accuracy of finite element model-based multi-organ deformable image
• registration. Med. Phys., 32(6):1647–1659, June 2005.

41

Reference-Classical Methods

[Shen02] Dinnan Shen and Christos Davatzikos. HAMMER: hierarchical attribute matching mechanism for elastic registration. IEEE Trans. Medical Imaging, 21 (11), November 2002. [Thirion98] J.-P. Thirion. Image matching as a diffusion process: an analogy with Maxwell’s demons. Medical Image Analysis, 2(3):243-260, 1998. [Foskey05] Mark Foskey, Brad Davis, Lav Goyal, Sha Chang, Ed Chaney, Nathalie Strehl, Sandrine Tomei, Julian Rosenman and Sarang Joshi. Large deformation three-dimensional image registration in image-guided radiation

• therapy. Phys. Med. Biol., 50: 5869-5892, December 2005.

[Lu04] Weiguo Lu, Ming-Li Chen, Gustavo H. Olivera, Kenneth J. Ruchala, and Thomas R. Mackie. Fast free-form deformable registration via calculus of

• variations. Phys. Med. Biol., 49: 3067-3087, June 2004.

[Kaus07] Michael R. Kaus, Kristy K. Brock, Vladimir Pekar, Laura A. Dawson, Alan M. Nichol, and David A. Jaffray. Assessment of a model-based deformable image registration approach for radiation therapy planning. Int. J. Radiation Oncology Biol. Phys., 68(2): 572-580, 2007. 42