INFLUENCE OF THE BOUNDARY CONDITIONS ON THE ... - CiteSeerX
INFLUENCE OF THE BOUNDARY CONDITIONS ON THE RESULT OF NON-LINEAR IMAGE REGISTRATION Ulf-Dietrich Braumann and Jens-Peer Kuska University of Leipzig Interdisciplinary Center for Bioinformatics Härtelstraße 16-18, D-04107 Leipzig ABSTRACT The focus of many non-parametric image registration algorithms lies on the solution of non-linear partial differential equations. We offer a simple solution procedure therefor based on discrete Fourier transform. Boundary conditions can strongly influence the result of the registration. The issue is investigated on the example of non-linear curvature-based registration. 1. INTRODUCTION Registration is a fundamental task in image processing and is used to match two or more pictures [1]. It has various applications in particular in medical image processing. The non-linear registration of images is widespread in various variants [2] for twodimensional (and in analogy for three-dimensional) image data. The goal of such non-linear procedures is to find a transformation u(x) : W ® W with x = (x1 , x2 )¨ Î W Ì R2 so that the reference R(x) and the transformed image T Ix - u(x)M have a minimum absolute difference, and the transformation u(x) fulfills certain smoothness constraints. The smoothness of the transformation is obtained by some physically motivated differential operators. Two of them are the elastic smoothness constraint [3, 4] and the fluid constraint [5, 6]. In the present paper we restrict our interest to the curvature-based image registration [7, 8, 9, 2, 10]. In contrast to the solution technique in [10] we offer a simple multiresolution algorithm based on discrete Fourier transform [11, 12]. The simple formulation allows the inclusion of different boundary conditions by choosing appropriate variants of discrete Fourier transforms. The importance of the boundary conditions depends on the class of images that has to be registered. For many cases the objects that should be aligned can be surrounded by a sufficient number of background pixels, as e. g. in CT or MR. For large images [13, 14, 15] without a sharp object boundary like histological sections as obtained from transmitted light microscopy this is not efficient. 2. METHOD 2.1. Curvature-based registration The curvature-based image registration bases on the variational problem [10] min ID(u) + ΑScurv (u)M (1)
defining a joint registration criterion consisting of the sum of squared differences D(u) =
2 1 JT Ix - u(x)M - R(x)N âx à 2 W
,
(2)
and a smoothing term Scurv (u) =
1 (Du)¨ × (Du) âx 2 àW
(3)
which was introduced in [7]1 and recently studied in [9], wherein D denotes the two-dimensional Laplace operator. According to the calculus of variations, a function u(x) minimizing equ. (1) necessarily should be a solution for the Euler-Lagrange equation -ΑD2 u(x) + f Ix, u(x)M = 0
(4)
f (x, u(x)) = JR(x) - T Ix - u(x)MN × ÑT Ix - u(x)M
(5)
with
and
¶4 u ¶4 u ¶4 u (x) + 2 2 2 (x) + 4 (x). (6) ¶x14 ¶x1 ¶x2 ¶x2 This equation (4) is highly non-linear and typically cannot be solved directly. One common way to get a smooth convergence for the solution is to extend u(x) towards a time dependent function u(x, t) with limt®¥ u(x, t) = u(x), and consequently to solve D2 u(x) =
¶u (x, t) = -ΑD2 u(x, t) + f Ix, u(x, t)M. ¶t
(7)
However, the transformation found by equ. (7) for u is strongly influenced by its boundary conditions. 2.2. Boundary conditions Due to the fourth order for the two spatial dimensions one has to choose a total of four boundary conditions for u(x Î G, t), ¶n u(x Î G, t), ¶2n u(x Î G, t) and ¶3n u(x Î G, t) where G is the boundary of the image and n is the normal direction on the boundary. For the image registration the spatial coordinate x is restricted to the rectangular area x Î [0, N] ´ [0, M] with N and M denoting the image extent in x1 and x2 direction, respectively, and one needs the boundary conditions on the four edges of this area. We consider three types of boundary conditions for the fourth order equation (4) in the space variable x:
u
German Research Foundation (DFG), grant BIZ – 6 1/2
1 The equation in [7] uses a regularized form of the Laplacian (D + Ε) instead of D to avoid the division by zero in the inverse operator.
(a) periodic boundary conditions (b) vanishing first and third derivative in normal direction to the boundary (c) no displacement on the boundary and vanishing second derivative on the boundary. The cases (a) and (b) have the obvious advantage, that a rigid translation of the T , i. e. u(x) = t with a constant vector t, does not violate the boundary conditions and fulfills the homogenous equation (4). None of the boundary conditions is consistent with a rigid transformation that includes a rotation around a point c within the image, i. e. u(x) = R(x-c)+t, with the rotation matrix R. The only reason for this choice of boundary conditions is that a fast solution method exists for them.
Laplace operator on every level is approximated by the equation 25.3.33 from [16]. For the computation of V (l) on every level l one has to solve a system of linear equations for the components of V (l) . The best way to solve the system of equations for periodic boundary conditions is to use the discrete Fourier coefficients [12] of M (l) -1 N (l) -1
V (l) m,n =
For vanishing first and third derivative on the boundary the discrete cos-transform (type-II DCT) with M (l) -1 N (l) -1
V (l) m,n = 2.3. Solution procedure Since we operate on digital images we have to determine the discrete solution U in a regular mesh with U (l) (t) = u(xi, j , t) for i = 1, . . . N (l) , j = 1, . . . , M (l) and l = 0, . . . , L. The superscript l denotes a certain discretization level, where l = 0 is the coarsest grid and l = L the finest grid. We will use N (l+1) = 2N (l) and M (l+1) = 2M (l) and introduce the operator R AU (l+1) (t)E = U (l) (t) for the restriction of solution U to the coarse grid and the operator P AU (l) (t)E = U (l+1) (t) for the prolongation to the finer level. It should be noticed, that P and R are not only doing a resampling on the finer and coarser grids, respectively. Both operators include a scaling operation because the arguments of the operators describe displacements on the mesh, e. g. in our case P has to scale the arguments by a factor of two, and R has to divide the the arguments by two. The time dependent equation (7) on the level l = 0 is discretized with the time step h by an implicit mid-point rule for the linear operator D2 and by an explicit Euler rule for the non-linear term f : Α h 2 (0) V (0) (t + h) - U (0) (t) = ID V (t + h) + D2U (0) (t)M 2 +h F (0) (U (0) (t)),
(8)
Α h 2 (l) ID V (t + h) + D2U (l) (t)M 2 h F (l) + IP[V (l-1) ](t + 1) + U (l) (t)M (9) 2
This equation uses the implicit mid-point rule also in the nonlinear term of F, but instead of an iteration of V (l) the prolongation of V (l-1) from the finer level is used as approximation to V (l) , of course the equation (9) may be iterated, but this is not done here. When V (L) (t + h) on the finest level is computed, the new U (l) (t + h) are determined by U (L) = V (L) U (l) = R AU (l+1) E with l = L - 1, . . . , 0.
4 1 Πm (l) â âVˆ cos C (l) KΜ + OG 2 N (l) M (l) Μ=0 Ν=0 Μ,Ν M ´ cos C
1 Πn KΝ + OG 2 N (l)
(12)
and for zero-displacement at the boundary and vanishing third derivative the discrete sin-transform (type-II DST) with (l)
V (l) m,n =
(l)
M -1 N -1 Π (m + 1) IΜ + 12 M 4 ˆ (l) sin C V G â â Μ,Ν N (l) M (l) Μ=0 Ν=0 M (l)
´ sin C
Π (n + 1) IΝ + 12 M G (13) N (l)
is used. (l) We only list the iteration equation for Vˆ Μ,Ν with periodic boundary conditions, so from equation (9) one gets (l) (l) Vˆ Μ,Ν (t + h) = -Uˆ Μ,Ν (t) +
(l) (l) 6 J2 Uˆ Μ,Ν (t) + h Fˆ Μ,Ν N (l) q IΜ, Ν; Ω(l) Μ , ΩΝ M
(14)
q(Μ, Ν; ΩΜ , ΩΝ ) = 6 - 4 Β :23 + 10 cos(Μ ΩΜ ) cos(Ν ΩΝ )
where F (0) is the discrete approximation of f on the mesh points of level l = 0, and V (l) is a temporary approximation for U (l) . Than, the finer levels l > 0 are computed by the equation V (l) (t + h) - U (l) (t) = -
mΜ nΝ 1 (l) Vˆ exp C2 Π ä K (l) + (l) OG . (11) (l) (l)â â Μ,Ν N M Μ=0 Ν=0 M N
77 A cos(Μ ΩΜ ) + cos(Ν ΩΝ )E 4 7 + A cos(2 Μ ΩΜ ) + cos(2 Ν ΩΝ )E 2 1 - A cos(3 Μ ΩΜ ) + cos(3 Ν ΩΝ ) 4 + cos(Μ ΩΜ - 2 Ν ΩΝ ) + cos(2 Μ ΩΜ - Ν ΩΝ ) -
+ cos(2 Μ ΩΜ + Ν ΩΝ ) + cos(Μ ΩΜ + 2 Ν ΩΝ )E> (15) 4
with Β = Α h/ ∆(l) and Ω(l) = 2 Π/ N(l) and Ω(l) = 2 Π/ M(l) , Ν Μ where ∆(l) is the space step used for the discretization in x- and y-direction. With respect to the discrete cos- and sin-transforms (as taken for zero-derivative and zero-displacement boundary conditions, resp.) one obtains similar equations. Only the dependence on q(Μ, Ν; ΩΜ , ΩΝ ) changes. For the discrete cos- and sin-transform one has to replace q(Μ, Ν; ΩΜ , ΩΝ ) with q(Μ+1/ 2, Ν+1/ 2; ΩΜ / 2, ΩΝ / 2) in equation (14).
(10) 3. RESULTS AND DISCUSSION
When the cycle U(t) S® V(t + h) S® U(t + h) is complete, the convergence of U(t + h) is checked and the time t is incremented Å Å by h until ÅÅÅU (L) (t + h) - U (L) (t)ÅÅÅ < Ε with Ε 1. The square of the
We illustrate the registration performance based on a synthetic image which was non-linearly distorted as shown in fig. 1. We
Fig. 1. The reference image R (left) and the deformed image T (right)
apply a simple procedure for solving equ. (7) based on discrete Fourier transform. This easily allows tests with different boundary conditions for u(x), namely periodic boundary conditions (cf. fig. 2), zero-derivative boundary conditions (cf. fig. 3), and zero-displacement boundary conditions (cf. fig. 4). The visualization of the results in figs. 2–4 show the registered image on the left, and on the right a regular mesh under the transformation x - u(x) printed over a line integral convolution [17] of the vector field u(x) as background. All registration iterations use Α = 64 and h = 64 and a total of 512 iterations over four resolution levels of the 256 ´ 256 sample images. Since the time integration scheme is unconditional stable, for all time steps h larger time steps and a fewer number of iterations may be used. In many cases, the periodic boundary conditions give the best results for the registration process. For special cases the zerodisplacement boundary conditions may be used. However, in this particular case with image contents close to the margin, the zeroderivative boundary condition gives the best result. A more realistic example is shown in fig. 5 using x-ray radiographs taken from [8]. Here the parameters Α = 5 and h = 64 with 1024 iterations are used. This case shows that the appropriate choice of boundary conditions is essential for the quality of the result. Inappropriate boundary conditions, like the zero-displacement condition in this example can hinder the convergence of the solution. 4. CONCLUSIONS
Fig. 2. Registration result (left) and displacement field (right) for periodic boundary conditions
We have presented a simple non-linear curvature-based image registration algorithm that allows various boundary conditions. The proper choice of the boundary conditions depends on the kind of images that should be registered. In our experiments the periodic boundary conditions turns out to be the most flexible one. The real Fourier transforms as for zero-displacement and zero-derivative boundary conditions are nearly twice as fast as the variant for periodic boundary conditions. It should be noted that we always have found a single transformation without the need of restarting the iteration. If required by the images, other variants of boundary conditions in vertical and horizontal direction may be used by the combination of different forms of discrete Fourier transforms in the space directions. 5. ACKNOWLEDGMENTS
Fig. 3. Registration result (left) and displacement field (right) for zero-derivative boundary conditions
The authors would like to thank the referees for their comments, especially for bringing reference [7] to our attention. 6. REFERENCES [1] Lisa Gottesfeld Brown, “A survey of image registration techniques,” ACM Computing Surveys, vol. 24, no. 4, pp. 325– 376, 1992. [2] Jan Modersitzki, Numerical Methods for Image Registration, Oxford University Press, New York, NY, 2004. [3] Gary E. Christensen, Deformable Shape Models for Anatomy, Ph.D. thesis, Electrical Engineering D. Sc. Dissertation, Washington University, St. Louis, Missouri, August 1994.
Fig. 4. Registration result (left) and displacement field (right) for zero-displacement boundary conditions
[4] Gary E. Christensen, Michael I. Miller, Ulf Grenander, and Michael W. Vannier, “Individualizing neuroanatomical at-
Fig. 5. Registration of two x-ray radiographs of different human hands, the first column shows the ‘deformed’ image T (top) and the reference image R (bottom). The three remaining columns illustrate the respective registration results (bottom) and their corresponding displacement fields (top) for the following cases: periodic boundary conditions (second column), zero derivative boundary conditions (third column), and zero displacement boundary conditions (last column).
[5]
[6]
[7]
[8]
[9]
[10]
[11]
lases using a massively parallel computer,” IEEE Computer, pp. 32–38, January 1996. Gary E. Christensen an Richard D. Rabbit and Michael I. Miller, “Deformable templates using large deformation kinematics,” IEEE Transactions on Image Processing, vol. 10, no. 5, pp. 1435–1447, 1996. Morten Bro-Nielsen and Claus Gramkow, “Fast fluid registration of medical images,” in VBC’96: Proceedings of the 4th International Conference on Visualization in Biomedical Computing. 1996, pp. 267–276, Springer-Verlag. Yali Amit, “A nonlinear variational problem for image matching,” SIAM Journal on Scientific Computing, vol. 15, no. 1, pp. 207–224, Jan. 1994. Bernd Fischer and Jan Modersitzki, “Curvature based registration with applications to MR-mammography,” in Computational Science - ICCS 2002, P.M.A. Sloot et al., Ed. 2002, pp. 203–206, Springer LNCS 2331. Bernd Fischer and Jan Modersitzki, “Curvature based image registration,” Journal of Mathematical Imaging and Vision, vol. 18, pp. 81–85, Jan. 2003. Bernd Fischer and Jan Modersitzki, “A unified approach to fast image registration and a new curvature based registration technique,” Linear Algebra and its Applications, vol. 380, pp. 107–124, 2004. Matteo Frigo and Steven G. Johnson, “FFTW: An adaptive software architecture for the FFT,” in Proc. 1998 IEEE Intl. Conf. Acoustics Speech and Signal Processing. 1998, vol. 3, pp. 1381–1384, IEEE.
[12] Matteo Frigo and Steven G. Johnson, “The design and implementation of FFTW3,” Proceedings of the IEEE, vol. 93, no. 2, 2005, special issue on "Program Generation, Optimization, and Adaptation". [13] Ulf-Dietrich Braumann, Jens-Peer Kuska, Jens Einenkel, Lars-Christian Horn, and Michael Höckel, “Quantification of tumour invasion fronts using 3D reconstructed histological serial sections,” in Bildverarbeitung für die Medizin 2004 – Algorithmen, Systeme, Anwendungen. Mar. 2004, Informatik aktuell, pp. 70–74, Springer-Verlag. [14] Ulf-Dietrich Braumann, Jens-Peer Kuska, Jens Einenkel, Lars-Christian Horn, and Michael Höckel, “How to quantify cervical carcinoma invasion fronts in 3D?,” in International Journal of Gynecological Cancer, Sept. 2004, vol. 14, Suppl. 1, p. 12. [15] Ulf-Dietrich Braumann, Jens-Peer Kuska, Jens Einenkel, Lars-Christian Horn, and Michael Höckel, “3D-reconstruction and quantification of cervical carcinoma invasion fronts from histological serial sections,” IEEE Transactions on Medical Imaging, accepted for publication. [16] Milton Abramowitz and Irene A. Stegun, Pocketbook of Mathematical Functions, Verlag Harri Deutsch, Frankfurt/Main, 1984. [17] Brian Cabral and Leith C. Leedom, “Imaging vector fields using line integral convolution,” Computer Graphics (SIGGRAPH ’93 Proceedings), vol. 27, pp. 263–272, 1993.
defining a joint registration criterion consisting of the sum of squared differences D(u) =
2 1 JT Ix - u(x)M - R(x)N âx à 2 W
,
(2)
and a smoothing term Scurv (u) =
1 (Du)¨ × (Du) âx 2 àW
(3)
which was introduced in [7]1 and recently studied in [9], wherein D denotes the two-dimensional Laplace operator. According to the calculus of variations, a function u(x) minimizing equ. (1) necessarily should be a solution for the Euler-Lagrange equation -ΑD2 u(x) + f Ix, u(x)M = 0
(4)
f (x, u(x)) = JR(x) - T Ix - u(x)MN × ÑT Ix - u(x)M
(5)
with
and
¶4 u ¶4 u ¶4 u (x) + 2 2 2 (x) + 4 (x). (6) ¶x14 ¶x1 ¶x2 ¶x2 This equation (4) is highly non-linear and typically cannot be solved directly. One common way to get a smooth convergence for the solution is to extend u(x) towards a time dependent function u(x, t) with limt®¥ u(x, t) = u(x), and consequently to solve D2 u(x) =
¶u (x, t) = -ΑD2 u(x, t) + f Ix, u(x, t)M. ¶t
(7)
However, the transformation found by equ. (7) for u is strongly influenced by its boundary conditions. 2.2. Boundary conditions Due to the fourth order for the two spatial dimensions one has to choose a total of four boundary conditions for u(x Î G, t), ¶n u(x Î G, t), ¶2n u(x Î G, t) and ¶3n u(x Î G, t) where G is the boundary of the image and n is the normal direction on the boundary. For the image registration the spatial coordinate x is restricted to the rectangular area x Î [0, N] ´ [0, M] with N and M denoting the image extent in x1 and x2 direction, respectively, and one needs the boundary conditions on the four edges of this area. We consider three types of boundary conditions for the fourth order equation (4) in the space variable x:
u
German Research Foundation (DFG), grant BIZ – 6 1/2
1 The equation in [7] uses a regularized form of the Laplacian (D + Ε) instead of D to avoid the division by zero in the inverse operator.
(a) periodic boundary conditions (b) vanishing first and third derivative in normal direction to the boundary (c) no displacement on the boundary and vanishing second derivative on the boundary. The cases (a) and (b) have the obvious advantage, that a rigid translation of the T , i. e. u(x) = t with a constant vector t, does not violate the boundary conditions and fulfills the homogenous equation (4). None of the boundary conditions is consistent with a rigid transformation that includes a rotation around a point c within the image, i. e. u(x) = R(x-c)+t, with the rotation matrix R. The only reason for this choice of boundary conditions is that a fast solution method exists for them.
Laplace operator on every level is approximated by the equation 25.3.33 from [16]. For the computation of V (l) on every level l one has to solve a system of linear equations for the components of V (l) . The best way to solve the system of equations for periodic boundary conditions is to use the discrete Fourier coefficients [12] of M (l) -1 N (l) -1
V (l) m,n =
For vanishing first and third derivative on the boundary the discrete cos-transform (type-II DCT) with M (l) -1 N (l) -1
V (l) m,n = 2.3. Solution procedure Since we operate on digital images we have to determine the discrete solution U in a regular mesh with U (l) (t) = u(xi, j , t) for i = 1, . . . N (l) , j = 1, . . . , M (l) and l = 0, . . . , L. The superscript l denotes a certain discretization level, where l = 0 is the coarsest grid and l = L the finest grid. We will use N (l+1) = 2N (l) and M (l+1) = 2M (l) and introduce the operator R AU (l+1) (t)E = U (l) (t) for the restriction of solution U to the coarse grid and the operator P AU (l) (t)E = U (l+1) (t) for the prolongation to the finer level. It should be noticed, that P and R are not only doing a resampling on the finer and coarser grids, respectively. Both operators include a scaling operation because the arguments of the operators describe displacements on the mesh, e. g. in our case P has to scale the arguments by a factor of two, and R has to divide the the arguments by two. The time dependent equation (7) on the level l = 0 is discretized with the time step h by an implicit mid-point rule for the linear operator D2 and by an explicit Euler rule for the non-linear term f : Α h 2 (0) V (0) (t + h) - U (0) (t) = ID V (t + h) + D2U (0) (t)M 2 +h F (0) (U (0) (t)),
(8)
Α h 2 (l) ID V (t + h) + D2U (l) (t)M 2 h F (l) + IP[V (l-1) ](t + 1) + U (l) (t)M (9) 2
This equation uses the implicit mid-point rule also in the nonlinear term of F, but instead of an iteration of V (l) the prolongation of V (l-1) from the finer level is used as approximation to V (l) , of course the equation (9) may be iterated, but this is not done here. When V (L) (t + h) on the finest level is computed, the new U (l) (t + h) are determined by U (L) = V (L) U (l) = R AU (l+1) E with l = L - 1, . . . , 0.
4 1 Πm (l) â âVˆ cos C (l) KΜ + OG 2 N (l) M (l) Μ=0 Ν=0 Μ,Ν M ´ cos C
1 Πn KΝ + OG 2 N (l)
(12)
and for zero-displacement at the boundary and vanishing third derivative the discrete sin-transform (type-II DST) with (l)
V (l) m,n =
(l)
M -1 N -1 Π (m + 1) IΜ + 12 M 4 ˆ (l) sin C V G â â Μ,Ν N (l) M (l) Μ=0 Ν=0 M (l)
´ sin C
Π (n + 1) IΝ + 12 M G (13) N (l)
is used. (l) We only list the iteration equation for Vˆ Μ,Ν with periodic boundary conditions, so from equation (9) one gets (l) (l) Vˆ Μ,Ν (t + h) = -Uˆ Μ,Ν (t) +
(l) (l) 6 J2 Uˆ Μ,Ν (t) + h Fˆ Μ,Ν N (l) q IΜ, Ν; Ω(l) Μ , ΩΝ M
(14)
q(Μ, Ν; ΩΜ , ΩΝ ) = 6 - 4 Β :23 + 10 cos(Μ ΩΜ ) cos(Ν ΩΝ )
where F (0) is the discrete approximation of f on the mesh points of level l = 0, and V (l) is a temporary approximation for U (l) . Than, the finer levels l > 0 are computed by the equation V (l) (t + h) - U (l) (t) = -
mΜ nΝ 1 (l) Vˆ exp C2 Π ä K (l) + (l) OG . (11) (l) (l)â â Μ,Ν N M Μ=0 Ν=0 M N
77 A cos(Μ ΩΜ ) + cos(Ν ΩΝ )E 4 7 + A cos(2 Μ ΩΜ ) + cos(2 Ν ΩΝ )E 2 1 - A cos(3 Μ ΩΜ ) + cos(3 Ν ΩΝ ) 4 + cos(Μ ΩΜ - 2 Ν ΩΝ ) + cos(2 Μ ΩΜ - Ν ΩΝ ) -
+ cos(2 Μ ΩΜ + Ν ΩΝ ) + cos(Μ ΩΜ + 2 Ν ΩΝ )E> (15) 4
with Β = Α h/ ∆(l) and Ω(l) = 2 Π/ N(l) and Ω(l) = 2 Π/ M(l) , Ν Μ where ∆(l) is the space step used for the discretization in x- and y-direction. With respect to the discrete cos- and sin-transforms (as taken for zero-derivative and zero-displacement boundary conditions, resp.) one obtains similar equations. Only the dependence on q(Μ, Ν; ΩΜ , ΩΝ ) changes. For the discrete cos- and sin-transform one has to replace q(Μ, Ν; ΩΜ , ΩΝ ) with q(Μ+1/ 2, Ν+1/ 2; ΩΜ / 2, ΩΝ / 2) in equation (14).
(10) 3. RESULTS AND DISCUSSION
When the cycle U(t) S® V(t + h) S® U(t + h) is complete, the convergence of U(t + h) is checked and the time t is incremented Å Å by h until ÅÅÅU (L) (t + h) - U (L) (t)ÅÅÅ < Ε with Ε 1. The square of the
We illustrate the registration performance based on a synthetic image which was non-linearly distorted as shown in fig. 1. We
Fig. 1. The reference image R (left) and the deformed image T (right)
apply a simple procedure for solving equ. (7) based on discrete Fourier transform. This easily allows tests with different boundary conditions for u(x), namely periodic boundary conditions (cf. fig. 2), zero-derivative boundary conditions (cf. fig. 3), and zero-displacement boundary conditions (cf. fig. 4). The visualization of the results in figs. 2–4 show the registered image on the left, and on the right a regular mesh under the transformation x - u(x) printed over a line integral convolution [17] of the vector field u(x) as background. All registration iterations use Α = 64 and h = 64 and a total of 512 iterations over four resolution levels of the 256 ´ 256 sample images. Since the time integration scheme is unconditional stable, for all time steps h larger time steps and a fewer number of iterations may be used. In many cases, the periodic boundary conditions give the best results for the registration process. For special cases the zerodisplacement boundary conditions may be used. However, in this particular case with image contents close to the margin, the zeroderivative boundary condition gives the best result. A more realistic example is shown in fig. 5 using x-ray radiographs taken from [8]. Here the parameters Α = 5 and h = 64 with 1024 iterations are used. This case shows that the appropriate choice of boundary conditions is essential for the quality of the result. Inappropriate boundary conditions, like the zero-displacement condition in this example can hinder the convergence of the solution. 4. CONCLUSIONS
Fig. 2. Registration result (left) and displacement field (right) for periodic boundary conditions
We have presented a simple non-linear curvature-based image registration algorithm that allows various boundary conditions. The proper choice of the boundary conditions depends on the kind of images that should be registered. In our experiments the periodic boundary conditions turns out to be the most flexible one. The real Fourier transforms as for zero-displacement and zero-derivative boundary conditions are nearly twice as fast as the variant for periodic boundary conditions. It should be noted that we always have found a single transformation without the need of restarting the iteration. If required by the images, other variants of boundary conditions in vertical and horizontal direction may be used by the combination of different forms of discrete Fourier transforms in the space directions. 5. ACKNOWLEDGMENTS
Fig. 3. Registration result (left) and displacement field (right) for zero-derivative boundary conditions
The authors would like to thank the referees for their comments, especially for bringing reference [7] to our attention. 6. REFERENCES [1] Lisa Gottesfeld Brown, “A survey of image registration techniques,” ACM Computing Surveys, vol. 24, no. 4, pp. 325– 376, 1992. [2] Jan Modersitzki, Numerical Methods for Image Registration, Oxford University Press, New York, NY, 2004. [3] Gary E. Christensen, Deformable Shape Models for Anatomy, Ph.D. thesis, Electrical Engineering D. Sc. Dissertation, Washington University, St. Louis, Missouri, August 1994.
Fig. 4. Registration result (left) and displacement field (right) for zero-displacement boundary conditions
[4] Gary E. Christensen, Michael I. Miller, Ulf Grenander, and Michael W. Vannier, “Individualizing neuroanatomical at-
Fig. 5. Registration of two x-ray radiographs of different human hands, the first column shows the ‘deformed’ image T (top) and the reference image R (bottom). The three remaining columns illustrate the respective registration results (bottom) and their corresponding displacement fields (top) for the following cases: periodic boundary conditions (second column), zero derivative boundary conditions (third column), and zero displacement boundary conditions (last column).
[5]
[6]
[7]
[8]
[9]
[10]
[11]
lases using a massively parallel computer,” IEEE Computer, pp. 32–38, January 1996. Gary E. Christensen an Richard D. Rabbit and Michael I. Miller, “Deformable templates using large deformation kinematics,” IEEE Transactions on Image Processing, vol. 10, no. 5, pp. 1435–1447, 1996. Morten Bro-Nielsen and Claus Gramkow, “Fast fluid registration of medical images,” in VBC’96: Proceedings of the 4th International Conference on Visualization in Biomedical Computing. 1996, pp. 267–276, Springer-Verlag. Yali Amit, “A nonlinear variational problem for image matching,” SIAM Journal on Scientific Computing, vol. 15, no. 1, pp. 207–224, Jan. 1994. Bernd Fischer and Jan Modersitzki, “Curvature based registration with applications to MR-mammography,” in Computational Science - ICCS 2002, P.M.A. Sloot et al., Ed. 2002, pp. 203–206, Springer LNCS 2331. Bernd Fischer and Jan Modersitzki, “Curvature based image registration,” Journal of Mathematical Imaging and Vision, vol. 18, pp. 81–85, Jan. 2003. Bernd Fischer and Jan Modersitzki, “A unified approach to fast image registration and a new curvature based registration technique,” Linear Algebra and its Applications, vol. 380, pp. 107–124, 2004. Matteo Frigo and Steven G. Johnson, “FFTW: An adaptive software architecture for the FFT,” in Proc. 1998 IEEE Intl. Conf. Acoustics Speech and Signal Processing. 1998, vol. 3, pp. 1381–1384, IEEE.
[12] Matteo Frigo and Steven G. Johnson, “The design and implementation of FFTW3,” Proceedings of the IEEE, vol. 93, no. 2, 2005, special issue on "Program Generation, Optimization, and Adaptation". [13] Ulf-Dietrich Braumann, Jens-Peer Kuska, Jens Einenkel, Lars-Christian Horn, and Michael Höckel, “Quantification of tumour invasion fronts using 3D reconstructed histological serial sections,” in Bildverarbeitung für die Medizin 2004 – Algorithmen, Systeme, Anwendungen. Mar. 2004, Informatik aktuell, pp. 70–74, Springer-Verlag. [14] Ulf-Dietrich Braumann, Jens-Peer Kuska, Jens Einenkel, Lars-Christian Horn, and Michael Höckel, “How to quantify cervical carcinoma invasion fronts in 3D?,” in International Journal of Gynecological Cancer, Sept. 2004, vol. 14, Suppl. 1, p. 12. [15] Ulf-Dietrich Braumann, Jens-Peer Kuska, Jens Einenkel, Lars-Christian Horn, and Michael Höckel, “3D-reconstruction and quantification of cervical carcinoma invasion fronts from histological serial sections,” IEEE Transactions on Medical Imaging, accepted for publication. [16] Milton Abramowitz and Irene A. Stegun, Pocketbook of Mathematical Functions, Verlag Harri Deutsch, Frankfurt/Main, 1984. [17] Brian Cabral and Leith C. Leedom, “Imaging vector fields using line integral convolution,” Computer Graphics (SIGGRAPH ’93 Proceedings), vol. 27, pp. 263–272, 1993.
