- Email: [email protected]
Generalized regularization term for non-parametric multimodal image registration
ARTICLE IN PRESS
Signal Processing 87 (2007) 2837–2842 www.elsevier.com/locate/sigpro
Fast communication
Generalized regularization term for non-parametric multimodal image registration$ Jorge Larrey-Ruiz, Juan Morales-Sa´nchez, Rafael Verdu´-Monedero Departamento de las Tecnologı´as de la Informacio´n y las Comunicaciones, Universidad Polite´cnica de Cartagena, Spain Received 5 December 2006; received in revised form 22 May 2007; accepted 23 May 2007 Available online 8 June 2007
Abstract In the field of non-rigid medical image registration, many regularizers based on first- or second-order derivatives have been studied. In this paper, a new regularization term based on fractional order derivatives is proposed for the registration of multimodal (e.g., medical) images. It can be seen as a generalization of the diffusion and curvature smoothing terms, but with this approach it is possible to obtain better registration results in terms of both similarity of the images and smoothness of the transformation. This registration scheme is tested on two realistic medical imaging scenarios, comparing the obtained results with the optimally registered diffusion and curvature cases. r 2007 Elsevier B.V. All rights reserved. Keywords: Image registration; Variational methods; Biomedical imaging
1. Introduction Image registration is the process of finding an optimal geometric transformation that aligns corresponding points in two views of an object [1]. Particularly, in medical imaging there are several applications that require a non-rigid (i.e., nonlinear) registration to correct the local differences between the images [2]. Non-rigid image registration can be either parametric, where the transformation is expanded in terms of some parameters of basis functions, or non-parametric (commonly referred to $ This work is partially supported by the Spanish Ministerio de Ciencia y Tecnologı´ a, under grant TEC2006-13338/TCM, and by the Consejerı´ a de Educacio´n y Cultura de Murcia, under grant 03122/PI/05. Corresponding author. Tel.: +34 968 33 8861; fax: +34 968 32 5973. E-mail address: [email protected] (J. Larrey-Ruiz).
as elastic registration), which is based on the regularized minimization of a distance measure, until an equilibrium state (described by a set of partial differential equations) is achieved [3]. The non-rigid registration is an ill-posed problem, because it does not satisfy Hadamard’s criterion (i.e., uniqueness of the solution, if one exists, is not guaranteed). A common method to solve this problem is to add some prior knowledge of the displacement field (e.g., to assume that a physical plausible deformation for medical images should be regular and smooth). A regularization term is then used to preferentially obtain more likely solutions. For instance, so-called elastic registration techniques deal with a regularization term given by the energy of an elastic deformation [4]. Another popular way to add the constraints is to use a Tikhonov regularization term [5] (i.e., to minimize the energy of first-order derivatives of the solution),
0165-1684/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2007.05.028
ARTICLE IN PRESS 2838
J. Larrey-Ruiz et al. / Signal Processing 87 (2007) 2837–2842
which leads to a diffusion operator that makes the displacement field smooth, since it penalizes oscillating deformations. This diffusion term was one of the first regularizers proposed to solve the dense image matching problem and has been widely used in the literature [6,7]. One relatively recent regularizer is the curvature term, based on second-order derivatives [8,9], which can make the non-linear registration more robust against a suboptimal preregistration step, i.e., if the regularization term uses only second-order derivatives of the displacement field on the boundary of the data sets (as addressed in [10]), then affine linear deformations are corrected automatically by this smoother. During the past two years, a great interest in finding new regularization terms suitable for nonparametric medical image registration has arisen. Some recent examples are the minimal curvature [10], the membrane energy [11] or the consistent symmetric [12] approaches, all of them based on the diffusion and/or curvature smoothing terms. This paper contributes with the generalized formulation of a new regularization term based on fractionalorder derivatives, which is proposed for the registration of multimodal images. The diffusion and curvature smoothing terms can be seen as particular cases, but with this regularizer it is possible to obtain better registration results in terms of both similarity of the images and smoothness of the transformation. An efficient DCT-based implementation is also proposed for it. The paper is organized as follows. We start out with the mathematical setting of the general registration problem, the introduction of the new regularization term and the deduction of its ‘‘weak’’ formulation within a non-parametric variational framework. Next, the new registration scheme is tested on two realistic examples, comparing the obtained results with the optimally registered diffusion and curvature cases. Finally, the main issues presented throughout the paper are discussed. 2. The variational framework Given two data sets, a reference R RðxÞ and a template T TðxÞ, with x 2 X 0; 1½d , the aim of image registration is to find a global and/or local transformation from T onto R in such a way that the transformed template matches the reference. Then the purpose of the registration is to determine a displacement field u uðxÞ such that T u Tðx uðxÞÞ is similar to R in the geometrical sense. It turns
out that this problem may be formulated in terms of a variational approach [6,12]. To this end we introduce the joint energy functional to be minimized: J½u ¼ D½R; T; u þ a S½u,
(1)
where D (related to the external forces) represents a distance measure and S (internal constraints) is the regularization term, which determines the smoothness of the displacement field u, penalizing arbitrary transformations that may lead to cracks, foldings, or other unwanted deformations. The parameter a40 controls the strength of the smoothness of the displacement vectors versus the similarity of the images. When dealing with medical images with different intensity levels (i.e., in a multimodal registration scenario), mutual information (MI) is typically used to measure the distance between them. Assuming that the transformation class is sufficient to align the images, the MI is maximum if T u matches R (i.e., if the images are perfectly aligned) [13]. D is obtained by writing the MI in terms of the Kullback–Leibler divergence, as addressed in [14]. 2.1. Proposed regularization term At this point, the regularization term to be considered in this work is introduced: d Z 1X S¼ hDs=2 ul ; Ds=2 ul iRd dx 2 l¼1 X d Z 1X ¼ kDs=2 ul k2 dx, ð2Þ 2 l¼1 X with s 2 ½1; 2, where D r2 r> r is the ddimensional Laplace operator, and where h; iRd denotes the dot (or inner) product in Rd . Note that if s ¼ 1, then Ds=2 ¼ r, and a diffusion registration is being performed; if s ¼ 2, then Ds=2 ¼ D, and a curvature registration is in this case performed. A precise definition of the operator Ds=2 can be found e.g., in [15], and references therein. According to the calculus of variations, a displacement field u minimizing Eq. (1) is necessarily a solution of the Euler–Lagrange equation: fðx; uÞ þ aA½uðxÞ ¼ 0,
(3)
subject to appropriate boundary conditions [9]. The force field fðx; uÞ, related to the distance measure, drives the deformation; its expression can be found e.g., in [14]. On the other hand, A½u is a partial
ARTICLE IN PRESS J. Larrey-Ruiz et al. / Signal Processing 87 (2007) 2837–2842
differential operator which can be deduced from the directional (or Gaˆteaux) derivative of the smoother: d Z X dS½u; v ¼ hDs=2 ul ; Ds=2 vl iRd dx l¼1
Z ¼
X
X
hA½u; viRd dx.
ð4Þ
Using the fact that Ds=2 ul is a locally summable function with respect to the Lebesgue measure in R Rd , i.e., x kDs=2 ul kdxo þ 1 for every bounded domain x X, and imposing Neumann boundary conditions, i.e., Ds=2 ul jqX ¼ 0, where qX is the boundary of X, we can use the formula of partial integration (see e.g., [16]) to get the equation of the generalized derivative of order s=2, in the Sobolev sense, of Ds=2 ul : d Z X l¼1
X
hDs=2 ul ; Ds=2 vl iRd dx ¼ ð1Þs
d Z X l¼1
ðDs=2 Þ2 ul vl dx. X
(5) Finally, from (4) and (5) we can obtain the partial differential operator !s d X s A½u ¼ ðDÞ u ¼ qxk xk u, (6) k¼1
2839
where ~ denotes the rearrangement of a d-dimensional matrix into a single column vector, x is the current iteration and t is the time-step. The identity matrix I and the matrix representation of A½u, A, are two-dimensional matrices of size N N, where N ¼ n1 n2 nd is the cardinal of the data sets to be registered. Usually, ~ uð0Þ is initialized to zero. The l numerical solution of Eq. (7) allows for different implementations (e.g., through the discrete approximation of the spatial fractional derivatives). Particularly, an efficient, OðN log NÞ implementation based on the discrete cosine transform (DCT) can be easily deduced from [8], taking as a starting point the matrix representation of the partial differential operator of the diffusion registration, Adiff , which is defined recursively from the matrix stencils for the discrete Laplace operator. At this point, it should be noted that A ¼ ðAdiff Þs for every order s 2 ½1; 2 [3]. Next, the term ðI þ taAÞ1 has to be pre- and postmultiplied by the product of the unitary cosine matrices that diagonalize the structured matrix A and implement the DCT and its inverse. As a result, the iterative semi-implicit (7) can be rewritten, without rearranging the d-dimensional matrices, as ðxÞ ulðxþ1Þ ¼ IDCTfð1 þ ta lj1 ;...;jd Þ1 DCTfuðxÞ l tf l gg,
(8)
which includes the diffusion and curvature cases: if s ¼ 1, then A½u ¼ Du is the partial differential operator for diffusion registration; if s ¼ 2, then A½u ¼ D2 u is the partial differential operator for curvature registration. For the more general case of s 21; 2½, the fractional power of the Laplacian might seem meaningless in the spatial domain, but it can be seen as an operator that partially includes features of both diffusion and curvature registration schemes, and which can be easily implemented in another domain (e.g., the Laplace or Fourier domains). One common way to solve a non-linear partial differential equation like (3) is to introduce an artificial time t and to compute the steady-state solution of the time-dependent partial differential equation via a time-marching algorithm, letting the time derivative qt uðx; tÞ equal to the negative Gaˆteaux derivative of the joint energy functional [12]. The resulting expression can be discretized in time and space, and then solved by employing the following semi-implicit iterative scheme: ðxÞ
~ ~ ulðxþ1Þ ¼ ðI þ taAÞ1 ð~ uðxÞ l tf l Þ;
l ¼ 1; . . . ; d, (7)
where l ¼ 1; . . . ; d, j l ¼ 1; . . . ; nl , and !s d X ðj l 1Þp 2 cos . lj1 ;...;jd ¼ 2d nl l¼1
(9)
Finally, it should be noted that the matrix inversion and the product of matrices in (7) have become a pointwise inversion and a pointwise product in (8). 3. Results In this section, the proposed regularization term is tested on two relevant medical imaging applications: the atlas matching (or scene-to-model) scenario and the multimodal (or intra-subject) registration scenario. In the first case, the aim is to register an axial slice obtained from a computed tomography (CT, Fig. 1(b)) and an image corresponding to an atlas or model (Fig. 1(a)). The MI before the registration process is MI ¼ 1:14. In the second simulation, the objective is to register an axial slice obtained from a positron emission tomography (PET) scan (Fig. 2(b)) and an axial slice from a CT of a human thorax (Fig. 2(a)). The MI before the registration is in this case MI ¼ 0:98.
ARTICLE IN PRESS 2840
J. Larrey-Ruiz et al. / Signal Processing 87 (2007) 2837–2842
Fig. 1. Scene-to-model registration: (a) reference; (b) template; (c) T u and uðxÞ, s ¼ 1; (d) regular grid deformed, s ¼ 1; (e) T u and uðxÞ, s ¼ 2; (f) regular grid deformed, s ¼ 2; (g) T u and uðxÞ, s ¼ 1:85; (h) regular grid deformed, s ¼ 1:85
Fig. 2. Multimodal registration. (a) reference image; (b) template; (c) T u ðxÞ and uðxÞ, s ¼ 1; (d) regular grid deformed, s ¼ 1; (e) T u ðxÞ and uðxÞ, s ¼ 2; (f) regular grid deformed, s ¼ 2; (g) T u ðxÞ and uðxÞ, s ¼ 1:15; (h) regular grid deformed, s ¼ 1:15.
In both experiments, exists at least one value of s in the interval 1; 2½ that outperforms, in a lower number of iterations of the registration algorithm, the optimal results obtained with the diffusion and curvature schemes in terms of both similarity of the images and smoothness of the computed transformation. These values of s were obtained by sweeping the fractional order of derivation from 1 to 2, in steps of 0:05. The optimal registration procedure consisted of two sequential steps, as addressed in [17]: (1) Initial estimation of the simulation parameters: ^ the value of For a small number of iterations x, the regularization parameter a^ that minimizes the joint energy functional (1) is obtained. The
proportionality r between x^ and a^ is then ^ a. computed and fixed, r ¼ x=^ (2) Optimal parameters computation: Keeping constant the value of r, the aim is to find the optimal registration parameters, ao and xo , computed as ao ¼ k a^ and xo ¼ k x^ ¼ r ao , where the multiplying parameter k is obtained as the minimum value from which the slope of the joint energy functional is close to zero (i.e., the algorithm has converged). The computed parameters allow for the best trade-off between the energies of the joint functional (1). In particular, it should be noted that with less iterations the algorithm could not have converged yet, and with more iterations D could possibly be lower, but at the expense of a higher
ARTICLE IN PRESS J. Larrey-Ruiz et al. / Signal Processing 87 (2007) 2837–2842 Table 1 Summary of the parameters involved in the first simulation (Fig. 1) Registration scheme
ao
xo
MI
S
Diffusion ðs ¼ 1Þ Curvature ðs ¼ 2Þ Proposed ðs ¼ 1:85Þ
300 20 103 16 103
650 400 320
1:39 1:42 1:44
0:046 0:002 0:001
2841
curvature regularizers can accomplish the minimization of (1) in a more straightforward way (i.e., the convergence mode is faster). However, an analytical, detailed study of the convergence has not been tackled by the authors yet (actually, in the literature this study has not been rigourously carried out for a general non-parametric variational case). 4. Conclusion
Table 2 Summary of the parameters involved in the second simulation (Fig. 2) Registration scheme
ao
xo
MI
S
Diffusion ðs ¼ 1Þ Curvature ðs ¼ 2Þ Proposed ðs ¼ 1:15Þ
450 175 103 300
800 350 300
1:60 1:26 1:62
0:079 0:081 0:076
S, and therefore the registration would not be optimal from a variational point of view (i.e., J would not be minimum). Tables 1 and 2 summarize the simulation parameters ao (i.e., optimal regularization parameter) and xo (i.e., optimal number of iterations), as well as the computed measure of similarity (i.e., MI) and the regularization energy S, for each registration scheme. In these simulations, the highest value of the similarity measure and the lowest value of the regularization energy correspond to the novel approach. As can be appreciated in Fig. 2(h) and Fig. 2(d), an advantage of using s ¼ 1:15 versus s ¼ 1 is that in this way global rigid transformations are partially corrected (and thus the resulting grid is smoother), because the regularizer includes a slight component of curvature registration. In the same way, the use of s ¼ 1:85 instead of s ¼ 2 avoids local curvatures that could be unlikely for the tissues under consideration, because in this case the regularizer includes a small component of diffusion registration, see Fig. 1(h) versus Fig. 1(f). Finally, it should be noted that several experiments over different types of medical images show that the proposed hybrid approach can perform the optimal registration in the lowest number of iterations xo . This is also true for the experiments under consideration, as shown in the tables. This fact can be intuitively understood if one thinks that the partial inclusion of features of both diffusion and
In this paper, a regularization term based on fractional-order derivatives is proposed for the optimal registration of multimodal (e.g., medical) images. It can be seen as a generalization of the diffusion and curvature smoothing terms which includes features of both regularizers, thus providing better registration results in terms of both similarity of the images and smoothness of the transformation. In other words, the novel approach generalizes these PDE-based non-parametric registration methods to a continuum of fractional derivative-based PDE techniques, providing a gradual and better implementation of the combined diffusion–curvature approach than using P a Rcomposite regularizer defined as e.g., S ¼ 12 dl¼1 X gkrul k2 þ ð1 gÞðDul Þ2 dx, where g 2 ½0; 1; for an instance, if g ¼ 0:5 in this latter example, the result of the registration is quite close to the diffusion behavior, but the proposed regularizer with s ¼ 1:5 delivers a registration result that can actually be considered as in between both diffusion and curvature approaches. However, it should be noted that neither first-order (diffusion), second-order (curvature), nor the novel (fractionalorder) regularization terms are entirely physically appropriate. In fact, the situation is rather that no constraint is appropriate everywhere in the images, because these typically may represent tissues with very different properties, which may therefore need to be treated differently.
References [1] B. Zitova´, J. Flusser, Image registration methods: a survey, Image Vision Comput. 21 (2003) 997–1000. [2] J. Maintz, M. Viergever, A survey of medical image registration, Med. Image Anal. 2 (1) (1998) 1–36. [3] J. Modersitzki, Numerical Methods for Image Registration, Oxford Science, Oxford, 2003. [4] R. Bajcsy, S. Kovacic, Multiresolution elastic matching, Comput. Vision Graphics Image Process. 46 (1) (1989) 1–21. [5] A.N. Tikhonov, V.Y. Arsenin, Solutions of Ill-posed Problems, Winston & Sons, Washington, DC, 1977.
ARTICLE IN PRESS 2842
J. Larrey-Ruiz et al. / Signal Processing 87 (2007) 2837–2842
[6] Y. Amit, A nonlinear variational problem for image matching, SIAM J. Sci. Comput. 15 (1) (1994) 207–224. [7] B. Fischer, J. Modersitzki, Fast diffusion registration, In: Contemporary Mathematics, vol. 313, Inverse Problems, Image Analysis, and Medical Imaging, American Mathematical Society, Providence, RI, 2002, pp. 117–129. [8] B. Fischer, J. Modersitzki, A unified approach to fast image registration and a new curvature based registration technique, Linear Algebra Appl. 308 (2004) 107–124. [9] U.-D. Braumann, J.-P. Kuska, Influence of the boundary conditions on the results of non-linear image registration, IEEE Int. Conf. Image Process. I ( (2005) 1129–1132. [10] S. Henn, K. Witsch, Image registration based on multiscale energy information, Multiscale Modelling Simulation 4 (2) (2005) 584–609. [11] V. Noblet, et al., Retrospective evaluation of a topology preserving non-rigid registration method, Med. Image Anal. 10 (2006) 366–384.
[12] Z. Zhang, et al., Consistent multi-modal non-rigid registration based on a variational approach, Pattern Recognition Lett. 27 (2006) 715–725. [13] P. Viola, W.M. Wells, Alignment by maximization of mutual information, Int. J. Comput. Vision 24 (1997) 137–154. [14] G. Hermosillo, et al., A variational approach to multi-modal image matching, Techicnal Report 4117, INRIA, February 2001. [15] M.M. Meerschaert, et al., Fractional vector calculus for fractional advection–dispersion, Physica A 367 (2006) 181–190. [16] I.N. Bronstein, et al., Handbook of Mathematics, Springer, Berlin, 2004. [17] J. Larrey-Ruiz, J. Morales-Sa´nchez, Optimal parameters selection for non-parametric image registration methods, Lecture Notes in Computer Science, Springer, Berlin, Vol. 4179, 2006, pp. 564–575.
Signal Processing 87 (2007) 2837–2842 www.elsevier.com/locate/sigpro
Fast communication
Generalized regularization term for non-parametric multimodal image registration$ Jorge Larrey-Ruiz, Juan Morales-Sa´nchez, Rafael Verdu´-Monedero Departamento de las Tecnologı´as de la Informacio´n y las Comunicaciones, Universidad Polite´cnica de Cartagena, Spain Received 5 December 2006; received in revised form 22 May 2007; accepted 23 May 2007 Available online 8 June 2007
Abstract In the field of non-rigid medical image registration, many regularizers based on first- or second-order derivatives have been studied. In this paper, a new regularization term based on fractional order derivatives is proposed for the registration of multimodal (e.g., medical) images. It can be seen as a generalization of the diffusion and curvature smoothing terms, but with this approach it is possible to obtain better registration results in terms of both similarity of the images and smoothness of the transformation. This registration scheme is tested on two realistic medical imaging scenarios, comparing the obtained results with the optimally registered diffusion and curvature cases. r 2007 Elsevier B.V. All rights reserved. Keywords: Image registration; Variational methods; Biomedical imaging
1. Introduction Image registration is the process of finding an optimal geometric transformation that aligns corresponding points in two views of an object [1]. Particularly, in medical imaging there are several applications that require a non-rigid (i.e., nonlinear) registration to correct the local differences between the images [2]. Non-rigid image registration can be either parametric, where the transformation is expanded in terms of some parameters of basis functions, or non-parametric (commonly referred to $ This work is partially supported by the Spanish Ministerio de Ciencia y Tecnologı´ a, under grant TEC2006-13338/TCM, and by the Consejerı´ a de Educacio´n y Cultura de Murcia, under grant 03122/PI/05. Corresponding author. Tel.: +34 968 33 8861; fax: +34 968 32 5973. E-mail address: [email protected] (J. Larrey-Ruiz).
as elastic registration), which is based on the regularized minimization of a distance measure, until an equilibrium state (described by a set of partial differential equations) is achieved [3]. The non-rigid registration is an ill-posed problem, because it does not satisfy Hadamard’s criterion (i.e., uniqueness of the solution, if one exists, is not guaranteed). A common method to solve this problem is to add some prior knowledge of the displacement field (e.g., to assume that a physical plausible deformation for medical images should be regular and smooth). A regularization term is then used to preferentially obtain more likely solutions. For instance, so-called elastic registration techniques deal with a regularization term given by the energy of an elastic deformation [4]. Another popular way to add the constraints is to use a Tikhonov regularization term [5] (i.e., to minimize the energy of first-order derivatives of the solution),
0165-1684/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2007.05.028
ARTICLE IN PRESS 2838
J. Larrey-Ruiz et al. / Signal Processing 87 (2007) 2837–2842
which leads to a diffusion operator that makes the displacement field smooth, since it penalizes oscillating deformations. This diffusion term was one of the first regularizers proposed to solve the dense image matching problem and has been widely used in the literature [6,7]. One relatively recent regularizer is the curvature term, based on second-order derivatives [8,9], which can make the non-linear registration more robust against a suboptimal preregistration step, i.e., if the regularization term uses only second-order derivatives of the displacement field on the boundary of the data sets (as addressed in [10]), then affine linear deformations are corrected automatically by this smoother. During the past two years, a great interest in finding new regularization terms suitable for nonparametric medical image registration has arisen. Some recent examples are the minimal curvature [10], the membrane energy [11] or the consistent symmetric [12] approaches, all of them based on the diffusion and/or curvature smoothing terms. This paper contributes with the generalized formulation of a new regularization term based on fractionalorder derivatives, which is proposed for the registration of multimodal images. The diffusion and curvature smoothing terms can be seen as particular cases, but with this regularizer it is possible to obtain better registration results in terms of both similarity of the images and smoothness of the transformation. An efficient DCT-based implementation is also proposed for it. The paper is organized as follows. We start out with the mathematical setting of the general registration problem, the introduction of the new regularization term and the deduction of its ‘‘weak’’ formulation within a non-parametric variational framework. Next, the new registration scheme is tested on two realistic examples, comparing the obtained results with the optimally registered diffusion and curvature cases. Finally, the main issues presented throughout the paper are discussed. 2. The variational framework Given two data sets, a reference R RðxÞ and a template T TðxÞ, with x 2 X 0; 1½d , the aim of image registration is to find a global and/or local transformation from T onto R in such a way that the transformed template matches the reference. Then the purpose of the registration is to determine a displacement field u uðxÞ such that T u Tðx uðxÞÞ is similar to R in the geometrical sense. It turns
out that this problem may be formulated in terms of a variational approach [6,12]. To this end we introduce the joint energy functional to be minimized: J½u ¼ D½R; T; u þ a S½u,
(1)
where D (related to the external forces) represents a distance measure and S (internal constraints) is the regularization term, which determines the smoothness of the displacement field u, penalizing arbitrary transformations that may lead to cracks, foldings, or other unwanted deformations. The parameter a40 controls the strength of the smoothness of the displacement vectors versus the similarity of the images. When dealing with medical images with different intensity levels (i.e., in a multimodal registration scenario), mutual information (MI) is typically used to measure the distance between them. Assuming that the transformation class is sufficient to align the images, the MI is maximum if T u matches R (i.e., if the images are perfectly aligned) [13]. D is obtained by writing the MI in terms of the Kullback–Leibler divergence, as addressed in [14]. 2.1. Proposed regularization term At this point, the regularization term to be considered in this work is introduced: d Z 1X S¼ hDs=2 ul ; Ds=2 ul iRd dx 2 l¼1 X d Z 1X ¼ kDs=2 ul k2 dx, ð2Þ 2 l¼1 X with s 2 ½1; 2, where D r2 r> r is the ddimensional Laplace operator, and where h; iRd denotes the dot (or inner) product in Rd . Note that if s ¼ 1, then Ds=2 ¼ r, and a diffusion registration is being performed; if s ¼ 2, then Ds=2 ¼ D, and a curvature registration is in this case performed. A precise definition of the operator Ds=2 can be found e.g., in [15], and references therein. According to the calculus of variations, a displacement field u minimizing Eq. (1) is necessarily a solution of the Euler–Lagrange equation: fðx; uÞ þ aA½uðxÞ ¼ 0,
(3)
subject to appropriate boundary conditions [9]. The force field fðx; uÞ, related to the distance measure, drives the deformation; its expression can be found e.g., in [14]. On the other hand, A½u is a partial
ARTICLE IN PRESS J. Larrey-Ruiz et al. / Signal Processing 87 (2007) 2837–2842
differential operator which can be deduced from the directional (or Gaˆteaux) derivative of the smoother: d Z X dS½u; v ¼ hDs=2 ul ; Ds=2 vl iRd dx l¼1
Z ¼
X
X
hA½u; viRd dx.
ð4Þ
Using the fact that Ds=2 ul is a locally summable function with respect to the Lebesgue measure in R Rd , i.e., x kDs=2 ul kdxo þ 1 for every bounded domain x X, and imposing Neumann boundary conditions, i.e., Ds=2 ul jqX ¼ 0, where qX is the boundary of X, we can use the formula of partial integration (see e.g., [16]) to get the equation of the generalized derivative of order s=2, in the Sobolev sense, of Ds=2 ul : d Z X l¼1
X
hDs=2 ul ; Ds=2 vl iRd dx ¼ ð1Þs
d Z X l¼1
ðDs=2 Þ2 ul vl dx. X
(5) Finally, from (4) and (5) we can obtain the partial differential operator !s d X s A½u ¼ ðDÞ u ¼ qxk xk u, (6) k¼1
2839
where ~ denotes the rearrangement of a d-dimensional matrix into a single column vector, x is the current iteration and t is the time-step. The identity matrix I and the matrix representation of A½u, A, are two-dimensional matrices of size N N, where N ¼ n1 n2 nd is the cardinal of the data sets to be registered. Usually, ~ uð0Þ is initialized to zero. The l numerical solution of Eq. (7) allows for different implementations (e.g., through the discrete approximation of the spatial fractional derivatives). Particularly, an efficient, OðN log NÞ implementation based on the discrete cosine transform (DCT) can be easily deduced from [8], taking as a starting point the matrix representation of the partial differential operator of the diffusion registration, Adiff , which is defined recursively from the matrix stencils for the discrete Laplace operator. At this point, it should be noted that A ¼ ðAdiff Þs for every order s 2 ½1; 2 [3]. Next, the term ðI þ taAÞ1 has to be pre- and postmultiplied by the product of the unitary cosine matrices that diagonalize the structured matrix A and implement the DCT and its inverse. As a result, the iterative semi-implicit (7) can be rewritten, without rearranging the d-dimensional matrices, as ðxÞ ulðxþ1Þ ¼ IDCTfð1 þ ta lj1 ;...;jd Þ1 DCTfuðxÞ l tf l gg,
(8)
which includes the diffusion and curvature cases: if s ¼ 1, then A½u ¼ Du is the partial differential operator for diffusion registration; if s ¼ 2, then A½u ¼ D2 u is the partial differential operator for curvature registration. For the more general case of s 21; 2½, the fractional power of the Laplacian might seem meaningless in the spatial domain, but it can be seen as an operator that partially includes features of both diffusion and curvature registration schemes, and which can be easily implemented in another domain (e.g., the Laplace or Fourier domains). One common way to solve a non-linear partial differential equation like (3) is to introduce an artificial time t and to compute the steady-state solution of the time-dependent partial differential equation via a time-marching algorithm, letting the time derivative qt uðx; tÞ equal to the negative Gaˆteaux derivative of the joint energy functional [12]. The resulting expression can be discretized in time and space, and then solved by employing the following semi-implicit iterative scheme: ðxÞ
~ ~ ulðxþ1Þ ¼ ðI þ taAÞ1 ð~ uðxÞ l tf l Þ;
l ¼ 1; . . . ; d, (7)
where l ¼ 1; . . . ; d, j l ¼ 1; . . . ; nl , and !s d X ðj l 1Þp 2 cos . lj1 ;...;jd ¼ 2d nl l¼1
(9)
Finally, it should be noted that the matrix inversion and the product of matrices in (7) have become a pointwise inversion and a pointwise product in (8). 3. Results In this section, the proposed regularization term is tested on two relevant medical imaging applications: the atlas matching (or scene-to-model) scenario and the multimodal (or intra-subject) registration scenario. In the first case, the aim is to register an axial slice obtained from a computed tomography (CT, Fig. 1(b)) and an image corresponding to an atlas or model (Fig. 1(a)). The MI before the registration process is MI ¼ 1:14. In the second simulation, the objective is to register an axial slice obtained from a positron emission tomography (PET) scan (Fig. 2(b)) and an axial slice from a CT of a human thorax (Fig. 2(a)). The MI before the registration is in this case MI ¼ 0:98.
ARTICLE IN PRESS 2840
J. Larrey-Ruiz et al. / Signal Processing 87 (2007) 2837–2842
Fig. 1. Scene-to-model registration: (a) reference; (b) template; (c) T u and uðxÞ, s ¼ 1; (d) regular grid deformed, s ¼ 1; (e) T u and uðxÞ, s ¼ 2; (f) regular grid deformed, s ¼ 2; (g) T u and uðxÞ, s ¼ 1:85; (h) regular grid deformed, s ¼ 1:85
Fig. 2. Multimodal registration. (a) reference image; (b) template; (c) T u ðxÞ and uðxÞ, s ¼ 1; (d) regular grid deformed, s ¼ 1; (e) T u ðxÞ and uðxÞ, s ¼ 2; (f) regular grid deformed, s ¼ 2; (g) T u ðxÞ and uðxÞ, s ¼ 1:15; (h) regular grid deformed, s ¼ 1:15.
In both experiments, exists at least one value of s in the interval 1; 2½ that outperforms, in a lower number of iterations of the registration algorithm, the optimal results obtained with the diffusion and curvature schemes in terms of both similarity of the images and smoothness of the computed transformation. These values of s were obtained by sweeping the fractional order of derivation from 1 to 2, in steps of 0:05. The optimal registration procedure consisted of two sequential steps, as addressed in [17]: (1) Initial estimation of the simulation parameters: ^ the value of For a small number of iterations x, the regularization parameter a^ that minimizes the joint energy functional (1) is obtained. The
proportionality r between x^ and a^ is then ^ a. computed and fixed, r ¼ x=^ (2) Optimal parameters computation: Keeping constant the value of r, the aim is to find the optimal registration parameters, ao and xo , computed as ao ¼ k a^ and xo ¼ k x^ ¼ r ao , where the multiplying parameter k is obtained as the minimum value from which the slope of the joint energy functional is close to zero (i.e., the algorithm has converged). The computed parameters allow for the best trade-off between the energies of the joint functional (1). In particular, it should be noted that with less iterations the algorithm could not have converged yet, and with more iterations D could possibly be lower, but at the expense of a higher
ARTICLE IN PRESS J. Larrey-Ruiz et al. / Signal Processing 87 (2007) 2837–2842 Table 1 Summary of the parameters involved in the first simulation (Fig. 1) Registration scheme
ao
xo
MI
S
Diffusion ðs ¼ 1Þ Curvature ðs ¼ 2Þ Proposed ðs ¼ 1:85Þ
300 20 103 16 103
650 400 320
1:39 1:42 1:44
0:046 0:002 0:001
2841
curvature regularizers can accomplish the minimization of (1) in a more straightforward way (i.e., the convergence mode is faster). However, an analytical, detailed study of the convergence has not been tackled by the authors yet (actually, in the literature this study has not been rigourously carried out for a general non-parametric variational case). 4. Conclusion
Table 2 Summary of the parameters involved in the second simulation (Fig. 2) Registration scheme
ao
xo
MI
S
Diffusion ðs ¼ 1Þ Curvature ðs ¼ 2Þ Proposed ðs ¼ 1:15Þ
450 175 103 300
800 350 300
1:60 1:26 1:62
0:079 0:081 0:076
S, and therefore the registration would not be optimal from a variational point of view (i.e., J would not be minimum). Tables 1 and 2 summarize the simulation parameters ao (i.e., optimal regularization parameter) and xo (i.e., optimal number of iterations), as well as the computed measure of similarity (i.e., MI) and the regularization energy S, for each registration scheme. In these simulations, the highest value of the similarity measure and the lowest value of the regularization energy correspond to the novel approach. As can be appreciated in Fig. 2(h) and Fig. 2(d), an advantage of using s ¼ 1:15 versus s ¼ 1 is that in this way global rigid transformations are partially corrected (and thus the resulting grid is smoother), because the regularizer includes a slight component of curvature registration. In the same way, the use of s ¼ 1:85 instead of s ¼ 2 avoids local curvatures that could be unlikely for the tissues under consideration, because in this case the regularizer includes a small component of diffusion registration, see Fig. 1(h) versus Fig. 1(f). Finally, it should be noted that several experiments over different types of medical images show that the proposed hybrid approach can perform the optimal registration in the lowest number of iterations xo . This is also true for the experiments under consideration, as shown in the tables. This fact can be intuitively understood if one thinks that the partial inclusion of features of both diffusion and
In this paper, a regularization term based on fractional-order derivatives is proposed for the optimal registration of multimodal (e.g., medical) images. It can be seen as a generalization of the diffusion and curvature smoothing terms which includes features of both regularizers, thus providing better registration results in terms of both similarity of the images and smoothness of the transformation. In other words, the novel approach generalizes these PDE-based non-parametric registration methods to a continuum of fractional derivative-based PDE techniques, providing a gradual and better implementation of the combined diffusion–curvature approach than using P a Rcomposite regularizer defined as e.g., S ¼ 12 dl¼1 X gkrul k2 þ ð1 gÞðDul Þ2 dx, where g 2 ½0; 1; for an instance, if g ¼ 0:5 in this latter example, the result of the registration is quite close to the diffusion behavior, but the proposed regularizer with s ¼ 1:5 delivers a registration result that can actually be considered as in between both diffusion and curvature approaches. However, it should be noted that neither first-order (diffusion), second-order (curvature), nor the novel (fractionalorder) regularization terms are entirely physically appropriate. In fact, the situation is rather that no constraint is appropriate everywhere in the images, because these typically may represent tissues with very different properties, which may therefore need to be treated differently.
References [1] B. Zitova´, J. Flusser, Image registration methods: a survey, Image Vision Comput. 21 (2003) 997–1000. [2] J. Maintz, M. Viergever, A survey of medical image registration, Med. Image Anal. 2 (1) (1998) 1–36. [3] J. Modersitzki, Numerical Methods for Image Registration, Oxford Science, Oxford, 2003. [4] R. Bajcsy, S. Kovacic, Multiresolution elastic matching, Comput. Vision Graphics Image Process. 46 (1) (1989) 1–21. [5] A.N. Tikhonov, V.Y. Arsenin, Solutions of Ill-posed Problems, Winston & Sons, Washington, DC, 1977.
ARTICLE IN PRESS 2842
J. Larrey-Ruiz et al. / Signal Processing 87 (2007) 2837–2842
[6] Y. Amit, A nonlinear variational problem for image matching, SIAM J. Sci. Comput. 15 (1) (1994) 207–224. [7] B. Fischer, J. Modersitzki, Fast diffusion registration, In: Contemporary Mathematics, vol. 313, Inverse Problems, Image Analysis, and Medical Imaging, American Mathematical Society, Providence, RI, 2002, pp. 117–129. [8] B. Fischer, J. Modersitzki, A unified approach to fast image registration and a new curvature based registration technique, Linear Algebra Appl. 308 (2004) 107–124. [9] U.-D. Braumann, J.-P. Kuska, Influence of the boundary conditions on the results of non-linear image registration, IEEE Int. Conf. Image Process. I ( (2005) 1129–1132. [10] S. Henn, K. Witsch, Image registration based on multiscale energy information, Multiscale Modelling Simulation 4 (2) (2005) 584–609. [11] V. Noblet, et al., Retrospective evaluation of a topology preserving non-rigid registration method, Med. Image Anal. 10 (2006) 366–384.
[12] Z. Zhang, et al., Consistent multi-modal non-rigid registration based on a variational approach, Pattern Recognition Lett. 27 (2006) 715–725. [13] P. Viola, W.M. Wells, Alignment by maximization of mutual information, Int. J. Comput. Vision 24 (1997) 137–154. [14] G. Hermosillo, et al., A variational approach to multi-modal image matching, Techicnal Report 4117, INRIA, February 2001. [15] M.M. Meerschaert, et al., Fractional vector calculus for fractional advection–dispersion, Physica A 367 (2006) 181–190. [16] I.N. Bronstein, et al., Handbook of Mathematics, Springer, Berlin, 2004. [17] J. Larrey-Ruiz, J. Morales-Sa´nchez, Optimal parameters selection for non-parametric image registration methods, Lecture Notes in Computer Science, Springer, Berlin, Vol. 4179, 2006, pp. 564–575.