1: Start with an initial guess |
2: Set |
3: for |
4: if |
5: |
6: |
7: end if |
8: Compute |
9: Compute |
10: Compute |
11: Solve |
|
for |
12: |
13: end for |
This paper proposes a numerical method for solving a non-rigid image registration model based on optimal mass transport. The main contribution of this paper is to address two issues. One is that we impose a proper periodic boundary condition, such that when the reference and template images are related by translation, or a combination of translation and non-rigid deformation, the numerical solution gives the underlying transformation. The other is that we design a numerical scheme that converges to the optimal transformation between the two images. As an additional benefit, our approach can decompose the transformation into translation and non-rigid deformation. Our numerical results show that the numerical solutions yield good-quality transformations for non-rigid image registration problems.
Citation: |
Figure 2. Standard 7-point stencil and wide stencil discretizations. (a) The stencil points of the discretization (15). (b) The stencil points of the discretization (17). (c) Wide stencil discretization: We apply a local coordinate rotation at the grid point $\mathbf{x}_{i, j}$ by the angle $\theta_{i, j}$. The grey dashed lines are the orthogonal axes $\{(\mathbf{e}_z)_{i, j}, (\mathbf{e}_w)_{i, j}\}$. The stencil length is $\sqrt{h}$ ($\sqrt{h}>h$). The grey stars are the stencil points $\mathbf{x}_{i, j}\pm\sqrt{h}(\mathbf{e}_z)_{i, j}$ and $\mathbf{x}_{i, j}\pm\sqrt{h}(\mathbf{e}_w)_{i, j}$. The unknowns at these stencil points are approximated by the bilinear interpolation from the neighboring points (black dots)
Figure 3. Optimal versus non-optimal transformations. (a) Constant images $R$ and $T$, where $\rho^T = \rho^R = 1$. (b) The optimal transformation $\phi^*$ obtained by our monotone numerical scheme. It is an identity mapping. The figure shows the deformed image of a square grid under $\phi^*$, which is the square grid itself. (c) The non-optimal transformation $\phi$ obtained by a non-monotone finite difference scheme in [3]. The figure shows that a square grid is severely deformed under $\phi$
Figure 4. Image registration using the periodic and Neumann boundary conditions, where $T$ and $R$ are related by translation. (a) Template image $T$. (b) Reference image $R$. (c) Transformed image $T_{\phi^*}$ under the periodic boundary condition. (d) Displacement of pixels from $T$ to $T_{\phi^*}$ under the periodic boundary condition, which is a pure translation. (e) A deformed grid obtained by applying the transformation $(\phi^*)^{-1}$ on a square grid. $(\phi^*)^{-1}$ is computed under the periodic boundary condition. The thick black lines show where the boundary of $\Omega = [0, 1]\times [0, 1]$ is moved to under $(\phi^*)^{-1}$. The color bar is the morphing magnitude $\mu$. The intensity of the color shows the degree of morphing effect under $(\phi^*)^{-1}$. (f) Displacement of pixels from $T$ to $T_{\phi^*}$ under the Neumann boundary condition. (g) A deformed grid obtained by applying the transformation $(\phi^*)^{-1}$ on a square grid. $(\phi^*)^{-1}$ is computed under the Neumann boundary condition
Figure 5. Image registration using the periodic and Neumann boundary conditions, where $T$ and $R$ are related by a combination of translation and dilation. (a) Template image $T$. (b) Reference image $R$. (c) Transformed image $T_{\phi^*}$ under the periodic boundary condition. (d1) Displacement of pixels from $T$ to $T_{\phi^*}$ under the periodic boundary condition. (d2) Decomposition of the displacement into a combination of translation component (green) and dilation component (red). (e) A deformed grid obtained by applying the transformation $(\phi^*)^{-1}$ on a square grid. $(\phi^*)^{-1}$ is computed under the periodic boundary condition. The thick black lines show where the boundary of $\Omega = [0, 1]\times [0, 1]$ is moved to under $(\phi^*)^{-1}$. The color bar is the morphing magnitude $\mu$. The intensity of the color shows the degree of morphing effect under $(\phi^*)^{-1}$. (f) Displacement of pixels from $T$ to $T_{\phi^*}$ under the Neumann boundary condition. (g) A deformed grid obtained by applying the transformation $(\phi^*)^{-1}$ on a square grid. $(\phi^*)^{-1}$ is computed under the Neumann boundary condition
Figure 6. Image registration using the periodic and Neumann boundary conditions, where $T$ and $R$ are related by a combination of translation and rotation. (a) Template image $T$. (b) Reference image $R$. (c) Transformed image $T_{\phi^*}$ under the periodic boundary condition. (d1) Displacement of pixels from $T$ to $T_{\phi^*}$ under the periodic boundary condition. (d2) Decomposition of the displacement into a combination of translation component (green) and local rotation component (red). (e) A deformed grid obtained by applying the transformation $(\phi^*)^{-1}$ on a square grid. $(\phi^*)^{-1}$ is computed under the periodic boundary condition. The thick black lines show where the boundary of $\Omega = [0, 1]\times [0, 1]$ is moved to under $(\phi^*)^{-1}$. The color bar is the morphing magnitude $\mu$. The intensity of the color shows the degree of morphing effect under $(\phi^*)^{-1}$. (f) Displacement of pixels from $T$ to $T_{\phi^*}$ under the Neumann boundary condition. (g) A deformed grid obtained by applying the transformation $(\phi^*)^{-1}$ on a square grid. $(\phi^*)^{-1}$ is computed under the Neumann boundary condition
Figure 7. Mass transport registration under periodic boundary condition. (a) Template image $T$. (b) Reference image $R$. (c) Transformed image $T_{\phi^*}$. (d) Difference between transformed image $T_{\phi^*}$ and $R$. (e) Pre-specified underlying transformation between $T$ and $R$. (f) Transformation given by the numerical scheme, which is a good approximation to the pre-specified underlying transformation in (e). (g) A deformed grid obtained by applying the transformation $(\phi^*)^{-1}$ on a square grid
Figure 8. Empirical two-step registration, where $T$ and $R$ are the same as Figure 7(a)-(b). The registration is implemented by FAIR package [38]. (a) Transformed image $T_\phi$. (b) Difference between transformed image $T_\phi$ and $R$. (c) Transformation given by the empirical approach, consisting of a rigid pre-registration (green arrows) and a non-rigid elastic deformation (red arrows). (d) A deformed grid obtained by applying the transformation $\phi^{-1}$ on a square grid
Figure 9. Medical image registration using the periodic boundary condition. (a) Template image $T$. (b) Reference image $R$. (c) Transformed image $T_{\phi^*}$. (d) Decomposition of the displacement into a combination of translation component (green) and non-rigid deformation component (red). (e) A deformed grid obtained by applying the transformation $(\phi^*)^{-1}$ on a square grid. $(\phi^*)^{-1}$ is computed under periodic boundary condition
Figure 10. Medical image registration using the periodic boundary condition. (a) Template image $T$. (b) Reference image $R$. (c) Transformed image $T_{\phi^*}$. (d) Decomposition of the displacement into a combination of translation component (green) and non-rigid deformation component (red). (e) A deformed grid obtained by applying the transformation $(\phi^*)^{-1}$ on a square grid. $(\phi^*)^{-1}$ is computed under periodic boundary condition
Figure 11. Image registration from Lena to Tiffany using the periodic boundary condition. (a) Template image $T$. (b) Reference image $R$. (c) Transformed image $T_{\phi^*}$. (d) Decomposition of the displacement into a combination of translation component (green) and non-rigid deformation component (red). (e) A deformed grid obtained by applying the transformation $(\phi^*)^{-1}$ on a square grid. $(\phi^*)^{-1}$ is computed under periodic boundary condition
Figure 12. Image registration using the periodic and Neumann boundary conditions, where $T$ and $R$ are given in Figure 5(a) and Figure 5(b). (a) Difference between $T$ and $R$. (b) Difference between $T_{\phi}$ and $R$, where $T_{\phi}$ is a transformed image under rigid registration only (or, $T_{\phi}$ is a translation of $T$ to align with $R$). (c) Difference between unmorphed transformed image $T_{\phi^*}^{unmorph}$ and $R$ under the periodic boundary condition. (d) Difference between unmorphed transformed image $T_{\phi^*}^{unmorph}$ and $R$ under the Neumann boundary condition
Table 1. Modified Levenberg-Marquardt algorithm
1: Start with an initial guess |
2: Set |
3: for |
4: if |
5: |
6: |
7: end if |
8: Compute |
9: Compute |
10: Compute |
11: Solve |
|
for |
12: |
13: end for |
Table 2.
Sum of the squared differences before registration
Examples | Example 2 | Example 3 | Example 4 | Example 5 | Example 6 | Example 7 |
| 0.47 | 0.52 | 0.61 | 0.64 | 0.69 | 0.59 |
| 0 | | | | | |
Table 3. A summary of morphing effect at a point (or an infinitesimal element)
net flow of mass/pixels | area change of a square element | change of mass/pixels intensity | morphing magnitude | color of a square element |
zero | invariance | invariance | | white |
inflow | compressed | increase | | red |
outflow | expanded | decrease | | blue |
Table 4. The errors of the motion fields (26). The values are computed for Examples 2-5 in Sections 6.3 and 6.4. In each example, the mass transport registration under the periodic boundary condition is compared against either Neumann boundary condition or elastic registration
Examples | Example 2 | Example 3 | Example 4 | Example 5 |
| Periodic: | Periodic: | Periodic: | Mass transport, periodic: |
Neumann: | Neumann: | Neumann: | Two-step empirical: |
Table 5.
Number of steps for convergence (residual tolerance
Example | Example 3 | |||
Image size | 100x100 | 200x200 | 400x400 | 800x800 |
Number of steps for convergence | 5 | 3 | 3 | 3 |
CPU time for corrections of translation kernels (sec) | 1.0 | 4.6 | 30 | 259 |
CPU time for the primary nonlinear solver (sec) | 3.1 | 7.3 | 58 | 1083 |
Total CPU time (sec) | 4.1 | 11.9 | 88 | 1342 |
Table 6.
Number of steps for convergence (residual tolerance
Examples | Example 3 | Example 4 | Example 5 | Example 6 | Example 7 |
Image size | 600x600 | ||||
Number of steps for convergence | 3 | 3 | 10 | 10 | 19 |
CPU time for the primary nonlinear solver (sec) | 147 | 152 | 668 | 627 | 1613 |
Table 7.
Residual versus the number of iteration
The number of iteration | 1 | 2 | 3 | 4 | 5 |
Residual | 1382 | 195 | 2.32 | 1.71 | 0.131 |
The number of iteration | 6 | 7 | 8 | 9 | 10 |
Residual | 0.0236 | 0.00492 | 9.50x10 | 3.42x10 | 9.41x10 |
G. Barles
and P. E. Souganidis
, Convergence of approximation schemes for fully nonlinear second order equations, Asymptotic Anal., 4 (1991)
, 271-283.
![]() ![]() |
|
J. -D. Benamou, Y. Brenier and K. Guittet, The Monge-Kantorovitch mass transfer and its
computational fluid mechanics formulation, Internat. J. Numer. Methods Fluids, 40 (2002),
21-30, ICFD Conference on Numerical Methods for Fluid Dynamics (Oxford, 2001).
doi: 10.1002/fld.264.![]() ![]() ![]() |
|
J.-D. Benamou
, B. D. Froese
and A. M. Oberman
, Two numerical methods for the elliptic Monge-Ampère equation, M2AN Math. Model. Numer. Anal., 44 (2010)
, 737-758.
doi: 10.1051/m2an/2010017.![]() ![]() ![]() |
|
K. Böhmer
, On finite element methods for fully nonlinear elliptic equations of second order, SIAM J. Numer. Anal., 46 (2008)
, 1212-1249.
doi: 10.1137/040621740.![]() ![]() ![]() |
|
S. C. Brenner
, T. Gudi
, M. Neilan
and L.-y. Sung
, $C^0$ penalty methods for the fully nonlinear Monge-Ampère equation, Math. Comp., 80 (2011)
, 1979-1995.
doi: 10.1090/S0025-5718-2011-02487-7.![]() ![]() ![]() |
|
C. Broit, Optimal registration of deformed images.
![]() |
|
L. G. Brown
, A survey of image registration techniques, ACM Computing Surveys (CSUR), 24 (1992)
, 325-376.
doi: 10.1145/146370.146374.![]() ![]() |
|
P. A. Browne
, C. J. Budd
, C. Piccolo
and M. Cullen
, Fast three dimensional r-adaptive mesh redistribution, J. Comput. Phys., 275 (2014)
, 174-196.
doi: 10.1016/j.jcp.2014.06.009.![]() ![]() ![]() |
|
C. J. Budd
and J. F. Williams
, Moving mesh generation using the parabolic Monge-Ampère equation, SIAM J. Sci. Comput., 31 (2009)
, 3438-3465.
doi: 10.1137/080716773.![]() ![]() ![]() |
|
L. A. Caffarelli
, Boundary regularity of maps with convex potentials. II, Ann. of Math., 144 (1996)
, 453-496.
doi: 10.2307/2118564.![]() ![]() ![]() |
|
K. Y. Chan and J. W. Wan, Reconstruction of missing cells by a killing energy minimizing
nonrigid image registration, in Engineering in Medicine and Biology Society (EMBC), 2013
35th Annual International Conference of the IEEE, IEEE, 2013,3000-3003.
doi: 10.1109/EMBC.2013.6610171.![]() ![]() |
|
R. Chartrand
, B. Wohlberg
, K. R. Vixie
and E. M. Bollt
, A gradient descent solution to the Monge-{K}antorovich problem, Appl. Math. Sci. (Ruse), 3 (2009)
, 1071-1080.
![]() ![]() |
|
Y. Chen and J. W. Wan, Monotone mixed narrow/wide stencil finite difference scheme for Monge-Ampère equation, arXiv preprint, arXiv: 1608.00644.
![]() |
|
G. E. Christensen,
Deformable Shape Models for Anatomy, PhD thesis, Washington University Saint Louis, Mississippi, 1994.
![]() ![]() |
|
M. G. Crandall
, H. Ishii
and P.-L. Lions
, User's guide to viscosity solutions of second order partial differential equations, Bull. Amer. Math. Soc. (N.S.), 27 (1992)
, 1-67.
doi: 10.1090/S0273-0979-1992-00266-5.![]() ![]() ![]() |
|
M. G. Crandall
and P.-L. Lions
, Viscosity solutions of Hamilton-Jacobi equations, Trans. Amer. Math. Soc., 277 (1983)
, 1-42.
doi: 10.1090/S0002-9947-1983-0690039-8.![]() ![]() ![]() |
|
E. J. Dean
and R. Glowinski
, Numerical methods for fully nonlinear elliptic equations of the Monge-Ampère type, Comput. Methods Appl. Mech. Engrg., 195 (2006)
, 1344-1386.
doi: 10.1016/j.cma.2005.05.023.![]() ![]() ![]() |
|
P. Dupuis
, U. Grenander
and M. I. Miller
, Variational problems on flows of diffeomorphisms for image matching, Quarterly of Applied Mathematics, 56 (1998)
, 587-600.
doi: 10.1090/qam/1632326.![]() ![]() ![]() |
|
X. Feng
, R. Glowinski
and M. Neilan
, Recent developments in numerical methods for fully nonlinear second order partial differential equations, SIAM Rev., 55 (2013)
, 205-267.
doi: 10.1137/110825960.![]() ![]() ![]() |
|
X. Feng
and M. Neilan
, Vanishing moment method and moment solutions for fully nonlinear second order partial differential equations, J. Sci. Comput., 38 (2009)
, 74-98.
doi: 10.1007/s10915-008-9221-9.![]() ![]() ![]() |
|
B. Fischer
and J. Modersitzki
, Fast inversion of matrices arising in image processing, Numerical Algorithms, 22 (1999)
, 1-11.
doi: 10.1023/A:1019194421221.![]() ![]() ![]() |
|
P. A. Forsyth and G. Labahn, Numerical methods for controlled Hamilton-Jacobi-Bellman PDEs in finance,
Journal of Computational Finance, 11 (2007), 1pp.
doi: 10.21314/JCF.2007.163.![]() ![]() |
|
B. D. Froese
, A numerical method for the elliptic Monge-Ampère equation with transport boundary conditions, SIAM J. Sci. Comput., 34 (2012)
, A1432-A1459.
doi: 10.1137/110822372.![]() ![]() ![]() |
|
B. D. Froese
and A. M. Oberman
, Convergent finite difference solvers for viscosity solutions of the elliptic Monge-Ampère equation in dimensions two and higher, SIAM J. Numer. Anal., 49 (2011)
, 1692-1714.
doi: 10.1137/100803092.![]() ![]() ![]() |
|
S. K. Godunov
, A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics, Mat. Sb. (N.S.), 47 (1959)
, 271-306.
![]() ![]() |
|
A. A. Goshtasby,
2-D and 3-D Image Registration: For Medical, Remote Sensing, and Industrial Applications, John Wiley & Sons, 2005.
![]() |
|
S. Haker and A. Tannenbaum, Optimal mass transport and image registration, in Variational
and Level Set Methods in Computer Vision, 2001. Proceedings. IEEE Workshop on, IEEE,
2001, 29-36.
doi: 10.1109/VLSM.2001.938878.![]() ![]() |
|
S. Haker
, L. Zhu
, A. Tannenbaum
and S. Angenent
, Optimal mass transport for registration and warping, International Journal of Computer Vision, 60 (2004)
, 225-240.
doi: 10.1023/B:VISI.0000036836.66311.97.![]() ![]() |
|
D. L. Hill, P. G. Batchelor, M. Holden and D. J. Hawkes, Medical image registration,
Physics in Medicine and Biology, 46 (2001), R1.
doi: 10.1088/0031-9155/46/3/201.![]() ![]() |
|
M. Irani
and S. Peleg
, Improving resolution by image registration, CVGIP: Graphical Models and Image Processing, 53 (1991)
, 231-239.
doi: 10.1016/1049-9652(91)90045-L.![]() ![]() |
|
M. Knott
and C. S. Smith
, On the optimal mapping of distributions, Journal of Optimization Theory and Applications, 43 (1984)
, 39-49.
doi: 10.1007/BF00934745.![]() ![]() ![]() |
|
N. V. Krylov
, The control of the solution of a stochastic integral equation, Teor. Verojatnost. i Primenen., 17 (1972)
, 111-128.
![]() ![]() |
|
K. Levenberg
, A method for the solution of certain non-linear problems in least squares, Quart. Appl. Math., 2 (1944)
, 164-168.
doi: 10.1090/qam/10666.![]() ![]() ![]() |
|
P. -L. Lions, Hamilton-Jacobi-Bellman equations and the optimal control of stochastic systems,
in Proceedings of the International Congress of Mathematicians, Vol. 1, 2 (Warsaw, 1983),
PWN, Warsaw, 1984,1403-1417.
![]() ![]() |
|
J. A. Maintz
and M. A. Viergever
, A survey of medical image registration, Medical Image Analysis, 2 (1998)
, 1-36.
doi: 10.1016/S1361-8415(01)80026-8.![]() ![]() |
|
D. W. Marquardt
, An algorithm for least-squares estimation of nonlinear parameters, J. Soc. Indust. Appl. Math., 11 (1963)
, 431-441.
doi: 10.1137/0111030.![]() ![]() ![]() |
|
J. Modersitzki,
Numerical Methods for Image Registration, Numerical Mathematics and Scientific Computation, Oxford University Press, New York, 2004, Oxford Science Publications.
![]() ![]() |
|
J. Modersitzki,
FAIR: Flexible Algorithms for Image Registration, vol. 6 of Fundamentals of Algorithms, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2009.
doi: 10.1137/1.9780898718843.![]() ![]() ![]() |
|
O. Museyko
, M. Stiglmayr
, K. Klamroth
and G. Leugering
, On the application of the Monge-Kantorovich problem to image registration, SIAM J. Imaging Sci., 2 (2009)
, 1068-1097.
doi: 10.1137/080721522.![]() ![]() ![]() |
|
A. M. Oberman
, Wide stencil finite difference schemes for the elliptic Monge-Ampère equation and functions of the eigenvalues of the Hessian, Discrete Contin. Dyn. Syst. Ser. B, 10 (2008)
, 221-238.
doi: 10.3934/dcdsb.2008.10.221.![]() ![]() ![]() |
|
V. I. Oliker
and L. D. Prussner
, On the numerical solution of the equation $(\partial^ 2z/\partial x^ 2)(\partial^ 2z/\partial y^ 2)-((\partial^ 2z/\partial x\partial y))^ 2 = f$ and its discretizations. Ⅰ, Numer. Math., 54 (1988)
, 271-293.
doi: 10.1007/BF01396762.![]() ![]() ![]() |
|
S. Osher
and C.-W. Shu
, High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations, SIAM J. Numer. Anal., 28 (1991)
, 907-922.
doi: 10.1137/0728049.![]() ![]() ![]() |
|
K. Rohr,
Landmark-based Image Analysis: Using Geometric and Intensity Models, vol. 21, Springer Science & Business Media, 2001.
doi: 10.1007/978-94-015-9787-6.![]() ![]() |
|
L.-P. Saumier
, M. Agueh
and B. Khouider
, An efficient numerical algorithm for the $L^ 2$
optimal transport problem with periodic densities, IMA J. Appl. Math., 80 (2015)
, 135-157.
doi: 10.1093/imamat/hxt032.![]() ![]() ![]() |
|
A. Sotiras
, C. Davatzikos
and N. Paragios
, Deformable medical image registration: A survey, IEEE Transactions on Medical Imaging, 32 (2013)
, 1153-1190.
doi: 10.1109/TMI.2013.2265603.![]() ![]() |
|
P. Thevenaz
, U. E. Ruttimann
and M. Unser
, A pyramid approach to subpixel registration based on intensity, IEEE Transactions on Image Processing, 7 (1998)
, 27-41.
doi: 10.1109/83.650848.![]() ![]() |
|
J.-P. Thirion
, Image matching as a diffusion process: An analogy with maxwell's demons, Medical Image Analysis, 2 (1998)
, 243-260.
doi: 10.1016/S1361-8415(98)80022-4.![]() ![]() |
|
E. F. Toro,
Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, Third edition. Springer-Verlag, Berlin, 2009.
![]() ![]() |
|
U. Trottenberg, C. W. Oosterlee and A. Schüller,
Multigrid, Academic Press, Inc., San Diego, CA, 2001, With contributions by A. Brandt, P. Oswald and K. Stüben.
![]() ![]() |
|
A. Trouvé
, Diffeomorphisms groups and pattern matching in image analysis, International Journal of Computer Vision, 28 (1998)
, 213-221.
![]() |
|
P. Viola
and W. M. Wells Ⅲ
, Alignment by maximization of mutual information, International Journal of Computer Vision, 24 (1997)
, 137-154.
![]() |
An example of image registration using the Neumann boundary condition.
(a) Template image
Standard 7-point stencil and wide stencil discretizations.
(a) The stencil points of the discretization (15).
(b) The stencil points of the discretization (17).
(c) Wide stencil discretization: We apply a local coordinate rotation at the grid point
Optimal versus non-optimal transformations.
(a) Constant images
Image registration using the periodic and Neumann boundary conditions, where
Image registration using the periodic and Neumann boundary conditions, where
Image registration using the periodic and Neumann boundary conditions, where
Mass transport registration under periodic boundary condition.
(a) Template image
Empirical two-step registration, where
Medical image registration using the periodic boundary condition.
(a) Template image
Medical image registration using the periodic boundary condition.
(a) Template image
Image registration from Lena to Tiffany using the periodic boundary condition.
(a) Template image
Image registration using the periodic and Neumann boundary conditions, where