fast fluid registration with dirichlet boundary conditions - RIT - People
FAST FLUID REGISTRATION WITH DIRICHLET BOUNDARY CONDITIONS: A TRANSFORM-BASED APPROACH Nathan D. Cahill1,2, J. Alison Noble2, David J. Hawkes3, and Lawrence A. Ray1 1
Research & Innovation, Carestream Health, Inc., Rochester, NY 14652, USA Wolfson Medical Vision Laboratory, University of Oxford, Oxford OX1 3PJ, UK 3 Centre for Medical Image Computing, University College London, London WC1E 6BT, UK 2
ABSTRACT Fluid registration is an example of a nonrigid image registration algorithm that uses a deformation field to define the transformation between two images. The velocity of the deformation field is governed by the Navier-Lamé equations, which can be discretized and solved numerically via fixed-point iteration. The fixed-point iteration generates a succession of linear PDE systems, which can be solved quickly via discrete Fourier transform (DFT) techniques, as shown in the prior art. The major drawback of this approach is that it is only applicable when the boundary conditions of the velocity field are assumed to be periodic. This paper shows that by considering the adjoint of the Navier-Lamé operator, the succession of linear PDE systems can be solved quickly via discrete sine transform (DST) techniques, generating velocity fields that satisfy Dirichlet boundary conditions (where the velocities are zero on the boundaries). 1. INTRODUCTION Any image registration algorithm relies on defining three components: a similarity measure indicating the similarity between two images, a class of geometric transformations relating the two images, and an optimization technique providing the mechanism for optimizing the similarity measure over the space of valid transformations. Nonrigid registration algorithms define transformations as being parametric (e.g., local rigid, local affine, thin-plate splines, B-splines, radial basis functions, and finite element method models) or nonparametric (diffeomorphic deformation fields). Modersitzki [1] provides an overview of nonparametric image registration methods and illustrates their flexibility. Nonparametric image registration techniques typically define the similarity measure as a combination of two terms: a term that indicates the actual image similarity, and a term that regularizes the deformation field in order to ensure smoothness and prevent folds or tears. One of the earliest nonparametric image registration techniques, elastic
1424406722/07/$20.00 ©2007 IEEE
712
registration [2], defines a regularizing term based on the linearized elastic potential of the deformation field. Fluid registration, originally presented in [3, 4], uses a similar regularizing term as elastic registration, but introduces a time-dependence, operating on the velocity of the deformation field instead of the displacement. Fluid registration has gained more traction in recent research, as it has shown [1] to be more flexible in dealing with complicated nonrigid deformations. Computational solutions to the fluid registration problem involve solving linear systems of partial differential equations (PDEs). Early techniques [3, 4] used iterative linear systems solutions that can be used with a variety of boundary conditions; however, such techniques can be slow to converge. More recent techniques [1, 5, 6] solve the PDE systems directly, by posing them as deconvolution problems that are solved via Fourier-based methods or inverse filtering. Although these techniques are much faster, the first two are only applicable when the boundary conditions are periodic, and the third only provides an approximate solution. In this paper, we show that the fluid registration problem can be posed in such a way that a fast transform-based approach can be used to find the deformation field subject to Dirichlet boundary conditions; i.e., the velocities of the resulting deformation field are zero on the boundaries. The remainder of this paper is organized as follows: Section 2 provides the mathematical background needed to understand the fluid registration problem and current solution schemes. Section 3 describes the alternative formulation of the fluid registration problem along with the proposed fast transform-based approach for its solution. Section 4 illustrates the application of this approach to the deformable registration of simulated mammograms. Finally, Section 5 presents some conclusions and suggestions for future work. 2. BACKGROUND This section presents a brief introduction to the fluid
ISBI 2007
Authorized licensed use limited to: Rochester Institute of Technology. Downloaded on August 12,2010 at 20:19:00 UTC from IEEE Xplore. Restrictions apply.
registration problem; next, it follows with a summary of current solution techniques.
Using the calculus of variations [7], nonparametric image registration techniques are transformed from functional optimization problems into systems of partial differential equations (PDE) known as Euler-Lagrange equations. This allows us to present the fluid registration technique in the following way. Given a reference image I ref , a floating image I float , both defined on the closure of the open set : , and given a deformation ĭ defined by: ĭ x : x u x , (1) where u x is the displacement of the deformation field, we define I ufloat to be the deformed image I float ĭ . Under this notation, the fluid registration problem involves solving the following system of partial differential equations, known as the Navier-Lamé equations: P'v O P v b x, u x , (2) subject to suitable boundary conditions (typically Dirichlet, Neumann, or periodic). The vector field v is the velocity of the deformation, related to u by:
v x, t :
d dt
u x, t w t u x, t >u x, t @v x, t .
(3)
The constants P and O are known as the Lamé constants, and are typically set to P 1 and O 0 . The force vector b relates to the Gâteaux derivative of the sum of squared differences (SSD) between the reference image and the deformed floating image:
b x, u x
I
ref
x I ufloat x I ufloat x .
(4)
The Navier-Lamé equations are semi-linear, so a convenient way to solve (2) is by fixed-point iteration. Given an initial estimate u 0 of the displacement, the
corresponding initial velocity v 0 is computed from (3). The velocity field v k 1 is implicitly defined by:
(5) P'v k 1 O P v k 1 b x, u k x , where u k is computed from Euler integration of (3). Equation (5) is a linear PDE system that can be solved for successive values of k until the iteration converges on a fixed velocity field v f . An alternative way to express (5) is with the NavierLamé operator L :
Lv k 1
b x, u k x ,
d1 , d 2 T
(6)
w , w T wx1 wx2
.
(8)
2.2. Current Solution Techniques A variety of techniques for numerically approximating the solution of the linear PDE system (6) have arisen in the literature. All techniques have involved discretizing (7) by making finite difference approximations to the gradient, divergence, and Laplacian. Early techniques [3, 4] used iterative linear system solvers, such as successive overrelaxation, to solve the discrete version of (6). A number of different authors [1, 5, 6] have recognized that (6) can be solved much more quickly than iterative linear system solvers would allow. Modersitzki [1] and BroNielsen and Gramkow [5] derive the eigenfunctions of the Navier-Lamé operator, and Modersitzki [1] shows how, under periodic boundary conditions, the knowledge of these eigenfunctions leads to a Fourier transform-based approach to the solution of (6). Xu and Dony [6] construct an approximate FIR inverse to the FIR filter that implements the discrete version of (6). The approximate inverse filter is convolved with the right-hand side of (6), yielding an approximate solution. Although this approach could theoretically be used with any boundary conditions, it would still yield only an approximate solution, rather than an exact solution (to within numerical precision).
3. PROPOSED SOLUTION TECHNIQUE In order to describe the technique for quickly solving the fluid registration problem subject to Dirichlet boundary conditions, we first provide an alternate formulation of (6). Next, we identify eigenfunctions of the linear operator in the alternative formulation and show how the knowledge of these eigenfunctions leads to a solution based on the discrete sine transform.
3.1. Alternative Formulation of Navier-Lamé Equations Given the Navier-Lamé operator (7), it is easy to compute its adjoint, L† :
L† P d T d I O P dd T . †
(9)
Multiplying (6) on both sides by the adjoint operator yields: (10) L†Lv k 1 det L v k 1 L†b x, u k x .
The determinant of L manipulation:
where
(7)
I is the 2 u 2 identity matrix, and d is given by: d
2.1. The Fluid Registration Problem
v
L P d T d I O P dd T ,
can be found by algebraic
713 Authorized licensed use limited to: Rochester Institute of Technology. Downloaded on August 12,2010 at 20:19:00 UTC from IEEE Xplore. Restrictions apply.
§ O 2 P d12 Pd 22 det L det¨¨ © O P d1d 2
defined on an M u N grid of pixels, then the discrete adjoint Navier-Lamé equations are given by: S 0,0 S 0,0 v k 1 ,1 S 0,1 b k ,1 S 0, 2 b k , 2 , (14)
O P d1d 2 ·¸ Pd12 O 2 P d 22 ¸¹
P O 2P d Td . 2
S
(11)
Hence, (10) can be written as:
2
(12)
or, alternatively, with Laplacian notation, as: P O 2P '2 v k 1 L†b x, u k x .
S
0, 0
v
k 1 , 2
S 1,1 b k ,1 S 1, 2 b k , 2 ,
where
L†b x, u k x ,
P O 2P d Td v k 1
0, 0
§0 ¨0 ©
(13)
This system of equations, which we will refer to as the adjoint Navier-Lamé equations, when considered in conjunction with Dirichlet boundary conditions, has eigenfunctions that can be exploited to form a fast discrete sine transform-based solver.
§ 0 ¨ O 2 P ¨ 0 ©
S
1, 2 T
S 0,1
1
0·
1
0¹
P O 2 P ¨ 1 4 1 ¸ ,
S 0, 0
P
S 1,1
0
·
2 O 3 P O 2 P ¸ , ¸ P 0 ¹
§ 1 0 ¨0 0 4 ¨ ©1 0 and is the convolution operator. S 0, 2
(15)
¸
O P
1
·
0¸ ¸ 1 ¹
(16)
,
(17)
3.2. Eigenfunctions In a manner similar to that of Modersitzki [1], we show that products of sine functions are eigenfunctions of the operator (11) with respect to Dirichlet boundary conditions.
Theorem 2
1:
0,1 2 .
:
Let
a, b N 0
For
3.4. Numerical Solution Scheme If we use the eigenfunctions as building blocks to construct a velocity field, we find that the velocity field can be represented by an expansion in sine waves:
v mk,n1
with
2
a b z 0 , let Z a ,b x1 , x2 : sin aSx1 sin bSx2 .
Then Z a,b is an eigenfunction of the operator det L (cf., equation (11)) with corresponding eigenvalue
2
with respect to the Dirichlet boundary conditions.
Proof: For fixed a and b , we define Z 'Z
2
2
2
4 MN
Z a,b . Hence
where
S a b Z ,
and therefore,
E i, j
Hence,
2
det L Z P O 2P S 4 a 2 b 2 Z NZ . Also, it is clear that Z x1 ,0 = Z 0, x2 = Z x1 ,1 = Z 1, x2 2
= 0 for all x1 , x2 >0,1@ . Ŷ
3.3. Finite Difference Approximation Using the knowledge of the eigenfunctions of (11), we can derive a numerical solution scheme for the adjoint Navier-Lamé equations. However, we first must derive the discrete form of (13) using a finite difference approximation. If we consider v k 1 and
b x, u k x
b k ,1 , b k , 2
T
v
k 1 ,1
, v k 1 , 2
Snj
(18)
i 0 j 0
¦ ¦ >vˆ i ,kj1 ,t E i, j @sin MSmi1 sin N 1 ,
M 1 N 1
Snj
i 0 j 0
4 P O 2 P cos MSi 1 cos
2 , Sj N 1
2
(19)
(20)
and where t 1,2 . Because the discrete sine transform is its own inverse [8], we can now present a fast algorithm for solving the linear PDE system (5) with Dirichlet boundary conditions:
' 'Z S 4 a 2 b 2 Z .
'2Z
¦ ¦ vˆ i ,kj1 sin MSmi1 sin N 1 .
M 1 N 1
This is, in fact, the discrete sine transform of the velocity field. With some algebraic manipulation (which is not included here because of space considerations), it can be shown that the left-hand sides of (14) are given by: S 0,0 S 0,0 v mk,n1 ,t
S 4 a 2 b 2 P O 2P
N a ,b
4 MN
Algorithm: 1. Compute the right-hand side of (14) by direct convolution. 2. Compute the sine transform of the result of step 1. 3. Divide each component of the result of step 2 by the corresponding E i, j of (20), yielding vˆ k 1 .
4. Compute the discrete sine transform of vˆ k 1 , yielding the velocity field v k 1 .
T
The only caveat to this algorithm is that the term E 0,0 is
to be vector images
identically zero, making step 3 undefined for this
714 Authorized licensed use limited to: Rochester Institute of Technology. Downloaded on August 12,2010 at 20:19:00 UTC from IEEE Xplore. Restrictions apply.
component. As the corresponding (DC) component of v k 1 is in the null space of det L anyway, vˆ 0k,01 ,t can
simply be assigned a value of zero. 4. BREAST IMAGE REGISTRATION EXAMPLE
In order to illustrate the results of fluid registration subject to Dirichlet boundary conditions, we consider a pair of simulated mammograms and contrast the resulting deformation field with one that was subject to periodic boundary conditions. The mammograms are simulated using the method of Hipwell et al. [9], and are provided by the authors of that article. They correspond to two different rotations and compressions of a 3-D MR breast image with fatty tissue suppressed. The top left image of Figure (1) is a mammogram simulated at 30% compression and 10° rotation, and the top right image is a mammogram simulated at 40% compression and 0° rotation of the 3-D MR breast image. The bottom images illustrate the deformation fields resulting from fluid registration of the left mammogram to the right mammogram, subject to Dirichlet boundary conditions (left) and periodic boundary conditions (right).
Figure (1): (top left) first mammogram, (top right) second mammogram, deformation field with Dirichlet (bottom left) and periodic (bottom right) boundary conditions.
Experiments verify that the proposed fast algorithm yields deformation fields identical to those found using successive overrelaxation [4], up to the precision of the convergence tolerance. 5. CONCLUSION
By analyzing an alternate form of the Navier-Lamé equations, we have derived a fast algorithm for fluid registration subject to Dirichlet boundary conditions. Breast images were used to illustrate the deformation field resulting from this technique, in contrast with a deformation field that was subject to periodic boundary conditions. One potential direction for future work is to extend this analysis to derive a fast algorithm for fluid registration
subject to Neumann boundary conditions. This should be possible within the same framework of using the adjoint Navier-Lamé equations, but with an expansion of the velocity field in cosine waves. In addition, the extension of both fast algorithms to 3-D fluid registration should be investigated, and compared to other fast 3-D fluid registration algorithms, such as the anisotropic multigrid method described by Crum [10]. 10. ACKNOWLEDGMENTS
N. D. Cahill thanks Bill Crum at the Centre for Medical Image Computing, University College London, for fruitful discussions on fluid registration, and John Hipwell and Christine Tanner, also at the Centre for Medical Image Computing, for providing the simulated mammograms. 11. REFERENCES [1] J. Modersitzki, Numerical Methods for Image Registration, Oxford University Press, December 2003. [2] C. Broit, “Optimal Registration of Deformed Images,” Ph.D. Thesis, Computer and Information Sciences, University of Pennsylvania, 1981. [3] G. Christensen, R. Rabbit, and M. Miller, “3D Brain Mapping Using a Deformable Neuroanatomy,” Physics in Medicine and Biology, Vol. 39, pp. 609-618, 1994. [4] G. Christensen, R. Rabbit, and M. Miller, “Deformable Templates Using Large Deformation Kinematics,” IEEE Transactions on Image Processing, Vol. 5, No. 10, pp. 1435-1447, October 1996. [5] M. Bro-Nielsen and C. Gramkow, “Fast Fluid Registration of Medical Images,” Proceedings of the 4th International Conference on Visualization in Biomedical Computing, LNCS 1131, pp. 267276, 1996. [6] X. Xu and R. Dony, “Fast Fluid Registration Using Inverse Filtering for Non-rigid Image Registration,” Proc. 3rd International Symposium on Biomedical Imaging: Macro to Nano, pp. 470-473, April 2006. [7] S. Fomin and I. Gelfand, Calculus of Variations, Dover Publications, 2000. [8] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes in C++: The Art of Scientific Computing, 2nd Ed., Cambridge University Press, 2002. [9] J. Hipwell, C. Tanner, W. Crum, and D. Hawkes, “X-ray Mammogram Registration: A Novel Validation Method,” Proc. 8th International Workshop on Digital Mammography, LNCS 4046, pp. 197-204, 2006. [10] W. Crum, C. Tanner, and D. Hawkes, “Anisotropic Multiscale Fluid Registration: Evaluation in Magnetic Resonance Breast Imaging,” Physics in Medicine and Biology, Vol. 50, pp. 51535174, 2005.
715 Authorized licensed use limited to: Rochester Institute of Technology. Downloaded on August 12,2010 at 20:19:00 UTC from IEEE Xplore. Restrictions apply.
Research & Innovation, Carestream Health, Inc., Rochester, NY 14652, USA Wolfson Medical Vision Laboratory, University of Oxford, Oxford OX1 3PJ, UK 3 Centre for Medical Image Computing, University College London, London WC1E 6BT, UK 2
ABSTRACT Fluid registration is an example of a nonrigid image registration algorithm that uses a deformation field to define the transformation between two images. The velocity of the deformation field is governed by the Navier-Lamé equations, which can be discretized and solved numerically via fixed-point iteration. The fixed-point iteration generates a succession of linear PDE systems, which can be solved quickly via discrete Fourier transform (DFT) techniques, as shown in the prior art. The major drawback of this approach is that it is only applicable when the boundary conditions of the velocity field are assumed to be periodic. This paper shows that by considering the adjoint of the Navier-Lamé operator, the succession of linear PDE systems can be solved quickly via discrete sine transform (DST) techniques, generating velocity fields that satisfy Dirichlet boundary conditions (where the velocities are zero on the boundaries). 1. INTRODUCTION Any image registration algorithm relies on defining three components: a similarity measure indicating the similarity between two images, a class of geometric transformations relating the two images, and an optimization technique providing the mechanism for optimizing the similarity measure over the space of valid transformations. Nonrigid registration algorithms define transformations as being parametric (e.g., local rigid, local affine, thin-plate splines, B-splines, radial basis functions, and finite element method models) or nonparametric (diffeomorphic deformation fields). Modersitzki [1] provides an overview of nonparametric image registration methods and illustrates their flexibility. Nonparametric image registration techniques typically define the similarity measure as a combination of two terms: a term that indicates the actual image similarity, and a term that regularizes the deformation field in order to ensure smoothness and prevent folds or tears. One of the earliest nonparametric image registration techniques, elastic
1424406722/07/$20.00 ©2007 IEEE
712
registration [2], defines a regularizing term based on the linearized elastic potential of the deformation field. Fluid registration, originally presented in [3, 4], uses a similar regularizing term as elastic registration, but introduces a time-dependence, operating on the velocity of the deformation field instead of the displacement. Fluid registration has gained more traction in recent research, as it has shown [1] to be more flexible in dealing with complicated nonrigid deformations. Computational solutions to the fluid registration problem involve solving linear systems of partial differential equations (PDEs). Early techniques [3, 4] used iterative linear systems solutions that can be used with a variety of boundary conditions; however, such techniques can be slow to converge. More recent techniques [1, 5, 6] solve the PDE systems directly, by posing them as deconvolution problems that are solved via Fourier-based methods or inverse filtering. Although these techniques are much faster, the first two are only applicable when the boundary conditions are periodic, and the third only provides an approximate solution. In this paper, we show that the fluid registration problem can be posed in such a way that a fast transform-based approach can be used to find the deformation field subject to Dirichlet boundary conditions; i.e., the velocities of the resulting deformation field are zero on the boundaries. The remainder of this paper is organized as follows: Section 2 provides the mathematical background needed to understand the fluid registration problem and current solution schemes. Section 3 describes the alternative formulation of the fluid registration problem along with the proposed fast transform-based approach for its solution. Section 4 illustrates the application of this approach to the deformable registration of simulated mammograms. Finally, Section 5 presents some conclusions and suggestions for future work. 2. BACKGROUND This section presents a brief introduction to the fluid
ISBI 2007
Authorized licensed use limited to: Rochester Institute of Technology. Downloaded on August 12,2010 at 20:19:00 UTC from IEEE Xplore. Restrictions apply.
registration problem; next, it follows with a summary of current solution techniques.
Using the calculus of variations [7], nonparametric image registration techniques are transformed from functional optimization problems into systems of partial differential equations (PDE) known as Euler-Lagrange equations. This allows us to present the fluid registration technique in the following way. Given a reference image I ref , a floating image I float , both defined on the closure of the open set : , and given a deformation ĭ defined by: ĭ x : x u x , (1) where u x is the displacement of the deformation field, we define I ufloat to be the deformed image I float ĭ . Under this notation, the fluid registration problem involves solving the following system of partial differential equations, known as the Navier-Lamé equations: P'v O P v b x, u x , (2) subject to suitable boundary conditions (typically Dirichlet, Neumann, or periodic). The vector field v is the velocity of the deformation, related to u by:
v x, t :
d dt
u x, t w t u x, t >u x, t @v x, t .
(3)
The constants P and O are known as the Lamé constants, and are typically set to P 1 and O 0 . The force vector b relates to the Gâteaux derivative of the sum of squared differences (SSD) between the reference image and the deformed floating image:
b x, u x
I
ref
x I ufloat x I ufloat x .
(4)
The Navier-Lamé equations are semi-linear, so a convenient way to solve (2) is by fixed-point iteration. Given an initial estimate u 0 of the displacement, the
corresponding initial velocity v 0 is computed from (3). The velocity field v k 1 is implicitly defined by:
(5) P'v k 1 O P v k 1 b x, u k x , where u k is computed from Euler integration of (3). Equation (5) is a linear PDE system that can be solved for successive values of k until the iteration converges on a fixed velocity field v f . An alternative way to express (5) is with the NavierLamé operator L :
Lv k 1
b x, u k x ,
d1 , d 2 T
(6)
w , w T wx1 wx2
.
(8)
2.2. Current Solution Techniques A variety of techniques for numerically approximating the solution of the linear PDE system (6) have arisen in the literature. All techniques have involved discretizing (7) by making finite difference approximations to the gradient, divergence, and Laplacian. Early techniques [3, 4] used iterative linear system solvers, such as successive overrelaxation, to solve the discrete version of (6). A number of different authors [1, 5, 6] have recognized that (6) can be solved much more quickly than iterative linear system solvers would allow. Modersitzki [1] and BroNielsen and Gramkow [5] derive the eigenfunctions of the Navier-Lamé operator, and Modersitzki [1] shows how, under periodic boundary conditions, the knowledge of these eigenfunctions leads to a Fourier transform-based approach to the solution of (6). Xu and Dony [6] construct an approximate FIR inverse to the FIR filter that implements the discrete version of (6). The approximate inverse filter is convolved with the right-hand side of (6), yielding an approximate solution. Although this approach could theoretically be used with any boundary conditions, it would still yield only an approximate solution, rather than an exact solution (to within numerical precision).
3. PROPOSED SOLUTION TECHNIQUE In order to describe the technique for quickly solving the fluid registration problem subject to Dirichlet boundary conditions, we first provide an alternate formulation of (6). Next, we identify eigenfunctions of the linear operator in the alternative formulation and show how the knowledge of these eigenfunctions leads to a solution based on the discrete sine transform.
3.1. Alternative Formulation of Navier-Lamé Equations Given the Navier-Lamé operator (7), it is easy to compute its adjoint, L† :
L† P d T d I O P dd T . †
(9)
Multiplying (6) on both sides by the adjoint operator yields: (10) L†Lv k 1 det L v k 1 L†b x, u k x .
The determinant of L manipulation:
where
(7)
I is the 2 u 2 identity matrix, and d is given by: d
2.1. The Fluid Registration Problem
v
L P d T d I O P dd T ,
can be found by algebraic
713 Authorized licensed use limited to: Rochester Institute of Technology. Downloaded on August 12,2010 at 20:19:00 UTC from IEEE Xplore. Restrictions apply.
§ O 2 P d12 Pd 22 det L det¨¨ © O P d1d 2
defined on an M u N grid of pixels, then the discrete adjoint Navier-Lamé equations are given by: S 0,0 S 0,0 v k 1 ,1 S 0,1 b k ,1 S 0, 2 b k , 2 , (14)
O P d1d 2 ·¸ Pd12 O 2 P d 22 ¸¹
P O 2P d Td . 2
S
(11)
Hence, (10) can be written as:
2
(12)
or, alternatively, with Laplacian notation, as: P O 2P '2 v k 1 L†b x, u k x .
S
0, 0
v
k 1 , 2
S 1,1 b k ,1 S 1, 2 b k , 2 ,
where
L†b x, u k x ,
P O 2P d Td v k 1
0, 0
§0 ¨0 ©
(13)
This system of equations, which we will refer to as the adjoint Navier-Lamé equations, when considered in conjunction with Dirichlet boundary conditions, has eigenfunctions that can be exploited to form a fast discrete sine transform-based solver.
§ 0 ¨ O 2 P ¨ 0 ©
S
1, 2 T
S 0,1
1
0·
1
0¹
P O 2 P ¨ 1 4 1 ¸ ,
S 0, 0
P
S 1,1
0
·
2 O 3 P O 2 P ¸ , ¸ P 0 ¹
§ 1 0 ¨0 0 4 ¨ ©1 0 and is the convolution operator. S 0, 2
(15)
¸
O P
1
·
0¸ ¸ 1 ¹
(16)
,
(17)
3.2. Eigenfunctions In a manner similar to that of Modersitzki [1], we show that products of sine functions are eigenfunctions of the operator (11) with respect to Dirichlet boundary conditions.
Theorem 2
1:
0,1 2 .
:
Let
a, b N 0
For
3.4. Numerical Solution Scheme If we use the eigenfunctions as building blocks to construct a velocity field, we find that the velocity field can be represented by an expansion in sine waves:
v mk,n1
with
2
a b z 0 , let Z a ,b x1 , x2 : sin aSx1 sin bSx2 .
Then Z a,b is an eigenfunction of the operator det L (cf., equation (11)) with corresponding eigenvalue
2
with respect to the Dirichlet boundary conditions.
Proof: For fixed a and b , we define Z 'Z
2
2
2
4 MN
Z a,b . Hence
where
S a b Z ,
and therefore,
E i, j
Hence,
2
det L Z P O 2P S 4 a 2 b 2 Z NZ . Also, it is clear that Z x1 ,0 = Z 0, x2 = Z x1 ,1 = Z 1, x2 2
= 0 for all x1 , x2 >0,1@ . Ŷ
3.3. Finite Difference Approximation Using the knowledge of the eigenfunctions of (11), we can derive a numerical solution scheme for the adjoint Navier-Lamé equations. However, we first must derive the discrete form of (13) using a finite difference approximation. If we consider v k 1 and
b x, u k x
b k ,1 , b k , 2
T
v
k 1 ,1
, v k 1 , 2
Snj
(18)
i 0 j 0
¦ ¦ >vˆ i ,kj1 ,t E i, j @sin MSmi1 sin N 1 ,
M 1 N 1
Snj
i 0 j 0
4 P O 2 P cos MSi 1 cos
2 , Sj N 1
2
(19)
(20)
and where t 1,2 . Because the discrete sine transform is its own inverse [8], we can now present a fast algorithm for solving the linear PDE system (5) with Dirichlet boundary conditions:
' 'Z S 4 a 2 b 2 Z .
'2Z
¦ ¦ vˆ i ,kj1 sin MSmi1 sin N 1 .
M 1 N 1
This is, in fact, the discrete sine transform of the velocity field. With some algebraic manipulation (which is not included here because of space considerations), it can be shown that the left-hand sides of (14) are given by: S 0,0 S 0,0 v mk,n1 ,t
S 4 a 2 b 2 P O 2P
N a ,b
4 MN
Algorithm: 1. Compute the right-hand side of (14) by direct convolution. 2. Compute the sine transform of the result of step 1. 3. Divide each component of the result of step 2 by the corresponding E i, j of (20), yielding vˆ k 1 .
4. Compute the discrete sine transform of vˆ k 1 , yielding the velocity field v k 1 .
T
The only caveat to this algorithm is that the term E 0,0 is
to be vector images
identically zero, making step 3 undefined for this
714 Authorized licensed use limited to: Rochester Institute of Technology. Downloaded on August 12,2010 at 20:19:00 UTC from IEEE Xplore. Restrictions apply.
component. As the corresponding (DC) component of v k 1 is in the null space of det L anyway, vˆ 0k,01 ,t can
simply be assigned a value of zero. 4. BREAST IMAGE REGISTRATION EXAMPLE
In order to illustrate the results of fluid registration subject to Dirichlet boundary conditions, we consider a pair of simulated mammograms and contrast the resulting deformation field with one that was subject to periodic boundary conditions. The mammograms are simulated using the method of Hipwell et al. [9], and are provided by the authors of that article. They correspond to two different rotations and compressions of a 3-D MR breast image with fatty tissue suppressed. The top left image of Figure (1) is a mammogram simulated at 30% compression and 10° rotation, and the top right image is a mammogram simulated at 40% compression and 0° rotation of the 3-D MR breast image. The bottom images illustrate the deformation fields resulting from fluid registration of the left mammogram to the right mammogram, subject to Dirichlet boundary conditions (left) and periodic boundary conditions (right).
Figure (1): (top left) first mammogram, (top right) second mammogram, deformation field with Dirichlet (bottom left) and periodic (bottom right) boundary conditions.
Experiments verify that the proposed fast algorithm yields deformation fields identical to those found using successive overrelaxation [4], up to the precision of the convergence tolerance. 5. CONCLUSION
By analyzing an alternate form of the Navier-Lamé equations, we have derived a fast algorithm for fluid registration subject to Dirichlet boundary conditions. Breast images were used to illustrate the deformation field resulting from this technique, in contrast with a deformation field that was subject to periodic boundary conditions. One potential direction for future work is to extend this analysis to derive a fast algorithm for fluid registration
subject to Neumann boundary conditions. This should be possible within the same framework of using the adjoint Navier-Lamé equations, but with an expansion of the velocity field in cosine waves. In addition, the extension of both fast algorithms to 3-D fluid registration should be investigated, and compared to other fast 3-D fluid registration algorithms, such as the anisotropic multigrid method described by Crum [10]. 10. ACKNOWLEDGMENTS
N. D. Cahill thanks Bill Crum at the Centre for Medical Image Computing, University College London, for fruitful discussions on fluid registration, and John Hipwell and Christine Tanner, also at the Centre for Medical Image Computing, for providing the simulated mammograms. 11. REFERENCES [1] J. Modersitzki, Numerical Methods for Image Registration, Oxford University Press, December 2003. [2] C. Broit, “Optimal Registration of Deformed Images,” Ph.D. Thesis, Computer and Information Sciences, University of Pennsylvania, 1981. [3] G. Christensen, R. Rabbit, and M. Miller, “3D Brain Mapping Using a Deformable Neuroanatomy,” Physics in Medicine and Biology, Vol. 39, pp. 609-618, 1994. [4] G. Christensen, R. Rabbit, and M. Miller, “Deformable Templates Using Large Deformation Kinematics,” IEEE Transactions on Image Processing, Vol. 5, No. 10, pp. 1435-1447, October 1996. [5] M. Bro-Nielsen and C. Gramkow, “Fast Fluid Registration of Medical Images,” Proceedings of the 4th International Conference on Visualization in Biomedical Computing, LNCS 1131, pp. 267276, 1996. [6] X. Xu and R. Dony, “Fast Fluid Registration Using Inverse Filtering for Non-rigid Image Registration,” Proc. 3rd International Symposium on Biomedical Imaging: Macro to Nano, pp. 470-473, April 2006. [7] S. Fomin and I. Gelfand, Calculus of Variations, Dover Publications, 2000. [8] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes in C++: The Art of Scientific Computing, 2nd Ed., Cambridge University Press, 2002. [9] J. Hipwell, C. Tanner, W. Crum, and D. Hawkes, “X-ray Mammogram Registration: A Novel Validation Method,” Proc. 8th International Workshop on Digital Mammography, LNCS 4046, pp. 197-204, 2006. [10] W. Crum, C. Tanner, and D. Hawkes, “Anisotropic Multiscale Fluid Registration: Evaluation in Magnetic Resonance Breast Imaging,” Physics in Medicine and Biology, Vol. 50, pp. 51535174, 2005.
715 Authorized licensed use limited to: Rochester Institute of Technology. Downloaded on August 12,2010 at 20:19:00 UTC from IEEE Xplore. Restrictions apply.
