• PDF
• Cite
• Share
Article Contents  Article Contents

# Numerical method for image registration model based on optimal mass transport

The second author is supported by Natural Sciences and Engineering Research Council of Canada (NSERC)

• 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.

Mathematics Subject Classification: Primary: 65N06, 65N22; Secondary: 35J96.

 Citation: • • Figure 1.  An example of image registration using the Neumann boundary condition. (a) Template image $T$. (b) Reference image $R$. (c) Underlying transformation between $T$ and $R$, which is a pure translation. (d) Transformation given by the Neumann boundary condition

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 . 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 . (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 $U^{(0)}=\frac{1}{2}(x^2+y^2)$. 2: Set $(c_1^{(-1)}, c_2^{(-1)})=(\infty, \infty)$, $\Pi \equiv [-\frac{1}{2}, \frac{1}{2}]\times[-\frac{1}{2}, \frac{1}{2}]$. 3: for $k = 0, 1, ...$ until convergence do 4:   if $(c_1^{(k-1)}, c_2^{(k-1)})\neq (0, 0)$ then 5:     $(c_1^{(k)}, c_2^{(k)}) = \underset{(c_1, c_2)\in \Pi} {\textrm{arg min}} \| R(U^{(k)} + c_1V_1 + c_2V_2) \|$. 6:     $U^{(k+\frac{1}{2})} = U^{(k)} + c_1^{(k)}V_1 + c_2^{(k)}V_2$. 7:   end if 8:   Compute $(a^{(k+\frac{1}{2})}, \theta^{(k+\frac{1}{2})})$ by (20). 9:   Compute $R^{(k+\frac{1}{2})}\equiv R(U^{(k+\frac{1}{2})})$ by (21). 10:   Compute $\mathbf{J}^{(k+\frac{1}{2})}\equiv\mathbf{J}[U^{(k+\frac{1}{2})}]$ by (22). 11:   Solve $\begin{array}{l} [\lambda I + (\mathbf{J}^{(k+\frac{1}{2})})^T \mathbf{J}^{(k+\frac{1}{2})}] E^{(k+\frac{1}{2})} \\ = -(\mathbf{J}^{(k+\frac{1}{2})})^T R^{(k+\frac{1}{2})}. \end{array}$ for $E^{(k+\frac{1}{2})}$. 12:   $U^{(k+1)} = U^{(k+\frac{1}{2})} + E^{(k+\frac{1}{2})}$. 13: end for

Table 2.  Sum of the squared differences before registration $\|\rho^T-\rho^R\|_{L_2(\Omega)}$, and after registration $\|\rho^{T_{\phi^*}}-\rho^R\|_{L_2(\Omega)}$. The values are computed for Examples 2-7 in Sections 6.3, 6.4 and 6.6. For each example, $\|\rho^R\|$ has been normalized to 1. The approach is the mass transport registration under the periodic boundary condition

 Examples Example 2 Example 3 Example 4 Example 5 Example 6 Example 7 $\|\rho^T-\rho^R\|_{L_2(\Omega)}$ 0.47 0.52 0.61 0.64 0.69 0.59 $\|\rho^{T_{\phi^*}}-\rho^R\|_{L_2(\Omega)}$ 0 $2\times 10^{-5}$ $1\times 10^{-3}$ $7\times 10^{-4}$ $9\times 10^{-4}$ $8\times 10^{-3}$

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 $\mu$ color of a square element zero invariance invariance $\mu=0$ white inflow compressed increase $\mu>0$ red outflow expanded decrease $\mu < 0$ 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 $\| \phi^*(\mathbf{x}) - \phi^*_{true}(\mathbf{x}) \|_{L_2(\Omega)}$ Periodic: $0$ Periodic: $0.0053$ Periodic: $0.066$ Mass transport, periodic: $0.0055$ Neumann: $0.056$ Neumann: $0.056$ Neumann: $0.088$ Two-step empirical: $0.011$

Table 5.  Number of steps for convergence (residual tolerance $10^{-4}$), and CPU time for Example 3 with different image sizes. Here (ⅰ) corrections of translation kernels refer to Line 4-7 of Table 1, and (ⅱ) the primary nonlinear solver refers to Line 8-12 of Table 1. The experiments are run in MATLAB

 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 $10^{-4}$), and CPU time for nonlinear solver only for Examples 3-7 with the same image sizes. The experiments are run in MATLAB

 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 $k$ for Example 5

 The number of iteration $k$ 1 2 3 4 5 Residual 1382 195 2.32 1.71 0.131 The number of iteration $k$ 6 7 8 9 10 Residual 0.0236 0.00492 9.50x10$^{-4}$ 3.42x10$^{-4}$ 9.41x10$^{-5}$
• Figures(12)

Tables(7)

## Article Metrics  DownLoad:  Full-Size Img  PowerPoint