Physically-Based Modeling for Medical Image Registration Huai-Ping Lee COMP 768 Presentation
What is Image Registration? A moving image I m ( x ) and a fixed image I f ( x ), can be 2D or 3D Rigid registration: find an affine transformation such that is F ( I m ) F closest to I f ( x ) Find a matrix to do the job, relatively easy Non-rigid registration: find a deformation field to do the job F 2
Medical Image Registration Human body is always deforming Rigid translation based on bony landmarks Fluid-flow-based deformable registration http://titan.radonc.unc.edu/~foskey/prost_web/warp_movies/ 3
Medical Image Registration a) b) c) Optical Flow FEM Registration Displacement field 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 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 µ ∇ 2 v ( x ) + ( µ + λ ) ∇ ( ∇· v ( x )) = f ( x , u ( x )) Diffusion Change of flux external force u : displacement field v = ∂ u ∂t + ( v · ∇ ) u v : velocity field, given by Convective derivative: rate of change of material at point x , given by chain rule d d u dt = ∂ u ∂t + ∂ u ∂ x ∂t = ∂ u ∂u i � v i ∂t + ∂ x ∂x i i =1 7
Viscous Fluid Model [Christensen 96] External force is computed with a similarity measure C ( I T ( x ) , I S u ( x , t )) = α � | I T ( x − u ( x , t )) − I S ( x ) | 2 d Ω 2 Ω By taking the derivative, they defined force at each d Ω � � f ( x , u ( x , t )) = − α ( I T ( x − u ( x , t )) − I S ( x )) ∇ T � � x − u ( x ,t ) Boundary condition is v = 0 The PDE is solved with finite difference method 8
Viscous Fluid Model [Christensen 96] Moving Image Target Image Deformed Image Displacement-x Displacement-y 9
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 F = − k x For multi-dimensional problems, it becomes σ = D ε σ : stress, measure of force per unit area within a body ε : strain, deformation of materials caused by stress D : elasticity matrix 11
Linear Elasticity ∂ u ∂ u 0 0 σ x ∂ x ∂ x ∂ u ∂ u 0 0 σ y ∂ y ∂ y ∂ u ∂ u 0 0 σ z ∂ z ∂ z = Lu ε = = σ = ∂ u ∂ x + ∂ u ∂ u ∂ u 0 τ xy ∂ y ∂ x ∂ y ∂ u ∂ y + ∂ u ∂ u ∂ u 0 τ yz ∂ z ∂ y ∂ z ∂ u ∂ x + ∂ u ∂ u ∂ u τ xz 0 ∂ z ∂ x ∂ z σ = D ε 12
Strain and Stress From Wikipedia 13
Linear Elasticity Potential energy of elastic body E = 1 � � σ T εd Ω + Fu d Ω 2 Ω Ω E = 1 � � u T L T DLu d Ω + Fu d Ω 2 Ω Ω Minimize E by driving its derivative to zero, and solve for u with finite element method 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 N i and node values u i N nodes � N i u i u = i =1 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 x B el i = L i N el i n (see Eq. 2.1) N nodes N nodes 1 � T B el T DB el � � u el 1 , . . . , u el u el j u el � � = j d Ω E N nodes i i 2 Ω i =1 j =1 N nodes � � F N el i u el + i d Ω Ω i =1 u el 1 , . . . , u el � � ∂ E To minimize E , take N nodes = 0 ; i = 1 , . . . , N nodes ∂ u el i N nodes � � T DB el � B el j u el F N el j d Ω = − ; i = 1 , . . . , N nodes d Ω i i Ω Ω j =1 18
FEM − Matrix Form K el u el = − F el T DB el K el Ω B el � i,j = j d Ω ( Stiffness Matrix ) i F el Ω F N el � i = d Ω i 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
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 1 0 0 0 ν ν (1 − ν ) (1 − ν ) 1 0 0 0 ν ν (1 − ν ) (1 − ν ) 1 0 0 0 ν ν E (1 − ν ) (1 − ν ) (1 − ν ) D = 1 − 2 ν 0 0 0 0 0 (1 + ν )(1 − 2 ν ) 2(1 − ν ) 1 − 2 ν 0 0 0 0 0 2(1 − ν ) 1 − 2 ν 0 0 0 0 0 2(1 − ν ) E : Young’s modulus ν : Poisson ratio Determined empirically! 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 Left Right Back tetrahedra, prisms may appear Subdivide them Case� 1 Case� 5 further Case� 2 Case� 6 Case� 3 Case� 7 Case� 4 Case� 8 25
Mesh Generation [Ferrant99] a) b) c) 26
Mesh Generation [Ferrant01] Subdivide tetrahedra when necessary Multi-resolution scheme 2� tets 1� tet� +� 1� pyramid 4� tets 4� tets 1� tet� +� 1� prism 2� prisms 2� tets� +� 1 prism 4� tets� +� 2� pyramids 1� tet� +� 2� pyramids 2� tets� +� 2� pyramids +� 1 pyramid 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 � � σ T εd Ω + ( I 1 ( x + u ( x )) − I 2 ( x )) 2 d Ω E = Ω Ω And use first order approximation I 1 ( x + u ( x )) ≈ I 1 ( x ) + ∇ I 1 ( x ) · u ( x ) 30
Volumetric Matching [Ferrant99] After taking the derivatives, we get four equations for each element: 4 � � T DB el � � � B el j + N el i ∇ I T 1 ∇ I 1 N el u el ∇ I 1 N el � � j d Ω = I 1 − I 2 d Ω i j i Ω Ω j =1 (5.6 ( K + G ) u = F a) b) c) Optical Flow FEM Registration Displacement field of FEM 31
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 v t − v t − 1 + Kv t = − F v t − 1 τ 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
Surface Matching [Ferrant00] Define a distance map of the edges in the target image: � � re D I ( x ) oint x . S Intuitively, external force should push the edge in source image toward the edge in target image � �� � F ( x ) = S min G exp ∇ I ( x ) D
Surface Matching [Ferrant00] S min : make F point towards a point that is closer to the edge � ���� � � � � � +1 if D I ( x ) x + ∇ I ( x ) > D I D S min = � ���� � � � � � − 1 if D I ( x ) x + ∇ I ( x ) < D I D G exp : avoid sticking on a wrong edge � +1 if k n ∇ I ( x ) > 0 G exp = − 1 if k n ∇ I ( x ) < 0 n : normal of surface k =1 if the region to be matched has higher intensity than surroundings; k= -1 otherwise
Recommend
More recommend
