Sage Journals: Discover world-class research

Abstract

Variational methods for image registration basically involve a regularizer to ensure that the resulting well-posed problem admits a solution. Different choices of regularizers lead to different deformations. On one hand, the conventional regularizers, such as the elastic, diffusion and curvature regularizers, are able to generate globally smooth deformations and generally useful for many applications. On the other hand, these regularizers become poor in some applications where discontinuities or steep gradients in the deformations are required. As is well-known, the total (TV) variation regularizer is more appropriate to preserve discontinuities of the deformations. However, it is difficult in developing an efficient numerical method to ensure that numerical solutions satisfy this requirement because of the non-differentiability and non-linearity of the TV regularizer. In this work we focus on computational challenges arising in approximately solving TV-based image registration model. Motivated by many efficient numerical algorithms in image restoration, we propose to use augmented Lagrangian method (ALM). At each iteration, the computation of our ALM requires to solve two subproblems. On one hand for the first subproblem, it is impossible to obtain exact solution. On the other hand for the second subproblem, it has a closed-form solution. To this end, we propose an efficient nonlinear multigrid (NMG) method to obtain an approximate solution to the first subproblem. Numerical results on real medical images not only confirm that our proposed ALM is more computationally efficient than some existing methods, but also that the proposed ALM delivers the accurate registration results with the desired property of the constructed deformations in a reasonable number of iterations.

Keywords

Image registration augmented Lagrangian method nonlinear multigrid method total variation regularization

Introduction

One of important tasks in image analysis is image registration. Generally speaking, it is the process of matching two or more images of the same or similar object obtained from different times, perspectives, and/or imaging sensors. This process can be done mathematically by computing a spatial geometric transformation of deformations (also known as displacement fields) that maps each point in one image onto a corresponding point in the other image with an optimal or meaningful manner. Image registration has played an important role in several areas of applications. Especially in medical applications, for example, it has been used routinely in medical diagnosis, treatment guidance and monitoring for providing complementary information. A good survey of the medical applications can be seen in literature^1–4 and the references therein.

Variational methods for image registration have been actively and extensively studied and applied in the field of image analysis. The main idea is to find spatial correspondences between two given images, a so-called reference image $R : Ω \subset ℝ^{2} \to V \subset ℝ$ and a so-called template image $T : Ω \subset ℝ^{2} \to V \subset ℝ$ , where the image domain Ω is assumed to be a rectangle. This is usually done by finding an optimal deformation $u : ℝ^{2} \to ℝ^{2}, u : x \mapsto u (x) = {(u_{1} (x), u_{2} (x))}^{⊤},$ such that the transformed version of the template image $T_{u} (x) = T (x + u (x))$ and the reference image R are spatially matched, according to an image-to-image dissimilarity measure, $D (u)$ . Without loss of generality, this work assumes that $Ω = {[0, 1]}^{2} \subset ℝ^{2}$ and $V = [0, 1]$ for 2 D gray intensity images. Note that the choice of Ω is somewhat arbitrary. The unit square $Ω = {[0, 1]}^{2}$ is adopted and used in all numerical sections throughout this work in order to make the performance comparison of the proposed numerical method with that of literature.⁵

If the image intensities of the given images R and T are comparable, one may have various choices for $D$ . Probably the most popular choice for the dissimilarity measure is provided by the so-called sum of squared differences (SSD)

D (u) = \frac{1}{2} \int_{Ω} {(T_{u} (x) - R (x))}^{2} d x

The registration task is then to solve the minimization problem

\min_{u} {D (u)}

(1)

where u is searched over a set of admissible functions

U

minimizing

D

. As is well-known, the minimization of

D

does not have a unique minimizer and it becomes necessary to impose a constraint on the solution u via a regularizer

R

. By modifying (1), the desired deformation u is a minimizer of the variational problem

\min_{u} {J_{α} (u) = D (u) + α R (u)}

(2)

where

α > 0

is a positive constant that compromises the quality of the similarity between

T_{u}

and R, measured by

D

, and the level of penalty for unwanted deformations, measured by

R

In this work, we focus on the regularizer of the form

R (u) = \sum_{l = 1}^{2} R^{TV} (u_{l})

(3)

where

R^{TV} (u_{l}) = \int_{Ω} | \nabla u_{l} | d x = \int_{Ω} \sqrt{u_{l, x}^{2} + u_{l, y}^{2}} d x

with

u_{l, x} = \frac{\partial u_{l}}{\partial x}

and

u_{l, y} = \frac{\partial u_{l}}{\partial y}

to ensure that the constructed deformation u is unique and preserves discontinuities of the deformations. This term is known as the total variation (TV) regularizer.^5–7 In contrast to other traditional regularizers generating globally smooth deformations, e.g. elastic⁸ diffusion⁹ and curvature regularizers^10,11 the TV regularizer is able to produces locally non-smooth deformations, which are required in matching several moved objects or partially occluded objects in medical applications, particularly at organ boundaries during the breathing induced organ motion.

In order to get a numerical solution of the variational problem (2), the standard gradient descent method can be applied. We first embed the associated Euler-Lagrange (EL) equation of (2) into a dynamic equation and drive it to a steady state. This yields the explicit scheme as given by

{\begin{matrix} \frac{u_{1}^{[k + 1]} - u_{1}^{[k]}}{τ} = α \nabla \cdot (\frac{\nabla u_{1}^{[k]}}{\sqrt{{(u_{1, x}^{[k]})}^{2} + {(u_{1, y}^{[k]})}^{2} + β}}) – f_{1} (u^{[k]}), \\ \frac{u_{2}^{[k + 1]} - u_{2}^{[k]}}{τ} = α \nabla \cdot (\frac{\nabla u_{2}^{[k]}}{\sqrt{{(u_{2, x}^{[k]})}^{2} + {(u_{2, y}^{[k]})}^{2} + β}}) – f_{2} (u^{[k]}), \end{matrix}

(4)

where

f_{l} (u^{[k]}) = (T_{u^{[k]}} - R) \partial_{u_{l}} T_{u^{[k]}}, k = 0, 1, 2, \dots

is the time step,

τ > 0

is the step size, and

β > 0

is a small real parameter to avoid zero division; see more details in literature.^12–14 This numerical scheme is easy to implement, but very slow to converge due to the constraint of stability conditions in the step size τ. In order to speed up the convergence of (4), we may linearize nonlinear ‘coefficients’ in the associated system and define the iteration step as follows:

{\begin{matrix} \frac{u_{1}^{[k + 1]} - u_{1}^{[k]}}{τ} - α \nabla \cdot (\frac{\nabla u_{1}^{[k + 1]}}{\sqrt{{(u_{1, x}^{[k]})}^{2} + {(u_{1, y}^{[k]})}^{2} + β}}) \\ = - f_{1} (u^{[k]}), \\ \frac{u_{2}^{[k + 1]} - u_{2}^{[k]}}{τ} - α \nabla \cdot (\frac{\nabla u_{2}^{[k + 1]}}{\sqrt{{(u_{2, x}^{[k]})}^{2} + {(u_{2, y}^{[k]})}^{2} + β}}) \\ = - f_{2} (u^{[k]}) \end{matrix}

(5)

The idea of the linearized gradient descent (LGD) method in (5) is to linearize and solve the EL equation via a fixed-point (FP) iteration in a similar way to the so-called Lagged-diffusivity method¹⁵ or Quasi-Newton scheme.^16,17 For each iteration, a linear system needs to be solved. As can be seen, both gradient descent methods in (4) and (5) share two drawbacks. At first, these gradient descent methods provide only the the approximate solutions of the original problem (2), since the TV regularizer $R (u)$ in (3) is replaced with

R_{β} (u) = \sum_{l = 1}^{2} R_{β}^{TV} (u_{l})

(6)

where

R_{β}^{TV} (u_{l}) = \int_{Ω} {| \nabla u_{l} |}_{β} d x = \int_{Ω} \sqrt{u_{l, x}^{2} + u_{l, y}^{2} + β} d x

is used to avoid non-differentiability of

R^{TV} (u_{l})

and thus approximated to get (4) and (5). On the second, the choice of β will effect on the computational efficiency of the numerical methods and the smoothness of the constructed deformations. Larger the β, more efficient the methods are, whereas more smooth the constructed deformations will be. Therefore the registration performance in smoothing deformations by these gradient descent methods is very sensitive to the parameter β. Indeed, there are a few different methods existing in the literature for solving the associated EL equation of (2). For example, the use of the nonlinear multigrid (NMG) methods can be found in literature.^5–7,18 However, these existing NMG methods also suffer from difficulties related to the non-differentiability and non-linearity of

R^{TV} (u_{l})

. To develop efficient numerical methods for TV-based image registration is still a challenging task and has been an active research area so far.

Recently, the variable-splitting methods are the well-known techniques in the field of image restoration for solving variational models, which require the minimization of nonlinear and non-differentiable functionals. To the best of the author’s knowledge, the variable-splitting methods for TV-based image registration model in (2) are still missing in the literature. In this paper, the proposed method is based on the so-called augmented Lagrangian method (ALM). We compare the performance of the proposed ALM with the LGD method in (5) and the NMG method developed by literature.⁵ Numerical results show that our proposed ALM is more computationally efficient than these existing methods

The rest of this paper is organized as follows. In the second section, we derive our ALM algorithm to solve the variational problem (2), following its numerical implementation in the third section. In the fourth section, numerical comparisons are carried out to confirm the effectiveness of the proposed ALM. Finally, some concluding remarks are made in the last section.

The proposed ALM

Before we derive the proposed ALM, we first rewrite (2) as

\min_{u} {J (u) = D (u) + α \sum_{l = 1}^{2} \int_{Ω} | \nabla u_{l} | d x}

(7)

Next, we introduce an auxiliary variable $w = (w_{1}, w_{2})$ such that $w_{l} = \nabla u_{l}$ and reformulate the minimization problem (7) to the the following constrained minimization problem

\min_{u, w} J (u, w)

(8)

J (u, w) = D (u) + α \sum_{l = 1}^{2} \int_{Ω} | w_{l} | d x

s . t . w_{l} = \nabla u_{l} for l=1, 2

For the minimization of (8), this work proposes to use the ALM and rewrite the constrained minimization problem (8) into an unconstrained minimization problem. We define the augmented Lagrangian functional for the above constrained minimization as follows

L (u, w; λ_{1}, λ_{2}) = J (u, w) + \sum_{l = 1}^{2} \frac{θ_{l}}{2} \int_{Ω} {| w_{l} - \nabla u_{l} |}^{2} d x + \sum_{l = 1}^{2} \int_{Ω} λ_{l} \cdot (w_{l} - \nabla u_{l}) d x

(9)

where

θ_{1}, θ_{2}

are the positive penalty parameters and

λ_{1}, λ_{2}

are the Lagrange multipliers.

Algorithm 1. ALM for TV-based image registration

1) Initialization: set m = 0, choose $θ_{1}, θ_{2} > 0$ and $λ_{1}^{(0)}, λ_{2}^{(0)}$ .

2) Repeat

2.1) Compute $(u^{(m)}, w^{(m)})$ as an (approximate) minimizer of the augmented Lagrangian functional with the Lagrange multipliers $λ_{1}^{(m - 1)}, λ_{2}^{(m - 1)}$ , i.e.,

(u^{(m)}, w^{(m)}) \approx \underset{u, w}{\arg \min} L (u, w; λ_{1}^{(m - 1)}, λ_{2}^{(m - 1)})

(10)

2.2) Update Lagrange multipliers

λ_{1}^{(m)} \leftarrow λ_{1}^{(m - 1)} + θ_{1} (w_{1}^{(m)} - \nabla u_{1}^{(m)})

(11)

λ_{2}^{(m)} \leftarrow λ_{2}^{(m - 1)} + θ_{2} (w_{2}^{(m)} - \nabla u_{2}^{(m)})

(12)

until a stopping rule for ALM method is satisfied.

We propose an iterative algorithm to solve the minimization of $L (u, w; λ_{1}, λ_{2})$ ; see Algorithm 1. Since u and w are coupled together in the minimization problem (9), it is very difficult to solve all variables simultaneously. We separate the minimization problem into two sub-problems and develop an alternating minimization procedure to approximate the solution. This process is repeated until one of the following stopping rules is satisfied:

{RelSSD}^{(m)} = \frac{D (u^{(m)})}{D (u^{(0)})} < ϵ_{1},

(13)

\max {\frac{∥ u_{1}^{(m)} - u_{1}^{(m - 1)} ∥}{∥ u_{1}^{(m)} ∥}, \frac{∥ u_{2}^{(m)} - u_{2}^{(m - 1)} ∥}{∥ u_{2}^{(m)} ∥}} < ϵ_{2},

(14)

\max {\frac{∥ w_{1}^{(m)} - \nabla u_{1}^{(m)} ∥}{| Ω |}, \frac{∥ w_{2}^{(m)} - \nabla u_{2}^{(m)} ∥}{| Ω |}} < ϵ_{3},

(15)

\max {\frac{∥ λ_{1}^{(m)} - λ_{1}^{(m - 1)} ∥}{∥ λ_{1}^{(m)} ∥}, \frac{∥ λ_{2}^{(m)} - λ_{2}^{(m - 1)} ∥}{∥ λ_{2}^{(m)} ∥}} < ϵ_{4},

(16)

m \geq ϵ_{5},

(17)

where ϵ₁, ϵ₂, ϵ₃ and ϵ₄ denote the predefined small positive numbers, and ϵ₅ is the maximum iteration. Here m denotes the index of the current iteration.

The two-subproblems are as follows: u-subproblem. Given $w = (w_{1}, w_{2}), λ_{1}$ and $λ_{2}$ , we search a minimizer $u = {(u_{1}, u_{2})}^{⊤}$ of the minimization problem:

\min_{u} J (u)

(18)

J (u) = D (u) + \sum_{l = 1}^{2} \frac{θ_{l}}{2} \int_{Ω} {| w_{l} - \nabla u_{l} |}^{2} d x - \sum_{l = 1}^{2} \int_{Ω} \nabla u_{l} \cdot λ_{l} d x

(19)

According to the calculus of variation, the solution of u -subproblem is determined by the associated EL equation

{\begin{array}{l} - θ_{1} △ u_{1} + f_{1} (u) = g_{1}, \\ - θ_{2} △ u_{2} + f_{2} (u) = g_{2}, \end{array}

(20)

subject to the homogeneous Neumann boundary conditions

\frac{\partial u_{l}}{\partial n} = 0

, where n denotes the unit outward normal vector on the image boundary

\partial Ω, f_{l} (u) = (T_{u} - R) \partial_{u_{l}} T_{u}, g_{l} = - θ_{1} \nabla \cdot w_{l} - \nabla \cdot λ_{l}

, and l = 1, 2.w-subproblem. Given u ,

λ_{1}

and

λ_{2}

, the following minimization problem

\min_{w} J (w)

(21)

J (w) = α \sum_{l = 1}^{2} \int_{Ω} | w_{l} | d x + \sum_{l = 1}^{2} \frac{θ_{l}}{2} \int_{Ω} {| w_{l} - \nabla u_{l} |}^{2} d x + \sum_{l = 1}^{2} \int_{Ω} w_{l} \cdot λ_{l} d x,

(22)

is solved for

w = (w_{1}, w_{2})

. Despite the fact that the variables

w_{1}

and

w_{2}

do not decouple we can still explicitly solve this minimization problem for w using the generalized shrinkage formula as used in:

w_{l} = \max (| S_{l} | - \frac{α}{θ_{l}}, 0) \frac{S_{l}}{| S_{l} |}

(23)

with the convection

(\frac{0}{| 0 |} = 0)

, where

S_{l} = \nabla u_{l} - \frac{λ_{l}}{θ_{l}}

and l = 1, 2. Note that several numerical methods, such as iterative reweighted least squares and interior-point method, have been applied to the minimization of subproblems as represented in (21), but these methods cannot usually converge to satisfactory solutions, especially for large-scale problems. To guarantee the convergence of minimization of subproblem (21), as well as to reduce the costs of numerical computation, generalized shrinkage/thresholding functions as given in (23) are employed as the solvers; see literature.^19–21

Numerical implementation

In this section, we present the details of how to solve the equation (20) and update the variables $w_{1}, w_{2}, λ_{1}$ , and $λ_{2}$ for each iteration.

Finite difference discretization

In order to discretize the EL equation (20), let

{(u_{l}^{h})}_{i, j} = u_{l}^{h} (x_{1_{i}}, x_{2_{j}}), {(u^{h})}_{i, j} = {({(u_{1}^{h})}_{i, j}, {(u_{2}^{h})}_{i, j})}^{⊤}

{(w_{l, 1}^{h})}_{i, j} = w_{l, 1}^{h} (x_{1_{i}}, x_{2_{j}}), {(w_{l, 2}^{h})}_{i, j} = w_{l, 2}^{h} (x_{1_{i}}, x_{2_{j}}),

{(w_{l}^{h})}_{i, j} = ({(w_{l, 1}^{h})}_{i, j}, {(w_{l, 2}^{h})}_{i, j}), {(g_{l}^{h})}_{i, j} = g_{l}^{h} (x_{1_{i}}, x_{2_{j}}),

denote the grid functions for l = 1, 2 with the grid spacing

h = (h_{1}, h_{2}) = (1 / n_{1}, 1 / n_{2})

, where the integers

n_{1} = 1 / h_{1}

and

n_{2} = 1 / h_{2}

are the number of uniform intervals in the x₁ and x₂ coordinate directions. Each grid point x in the discretized domain

Ω_{h} \subset Ω

is given by

x = {(x_{1_{i}}, x_{2_{j}})}^{⊤} = {(\frac{(2 i - 1) h_{1}}{2}, \frac{(2 j - 1) h_{2}}{2})}^{⊤}

for

1 \leq i \leq n_{1}

and

1 \leq j \leq n_{2}

Applying the finite difference approximations with (20), the discrete EL equation at a grid point $(x_{1_{i}}, x_{2_{j}})$ over the discrete domain Ω _h is given by

{\begin{matrix} ︸ \underset{N_{1}^{h} {(u^{h})}_{i, j}}{- θ_{1} △^{h} {(u_{1}^{h})}_{i, j} + f_{1}^{h} {(u^{h})}_{i, j}} = {(g_{1}^{h})}_{i, j}, \\ ︸ \underset{N_{2}^{h} {(u^{h})}_{i, j}}{- θ_{2} △^{h} {(u_{2}^{h})}_{i, j} + f_{2}^{h} {(u^{h})}_{i, j}} = {(g_{2}^{h})}_{i, j}, \end{matrix}

(24)

with the following notation

δ_{x_{1}}^{-} {(u_{l}^{h})}_{i, j} = ({(u_{l}^{h})}_{i, j} - {(u_{l}^{h})}_{i - 1, j}) / h_{1},

δ_{x_{1}}^{+} {(u_{l}^{h})}_{i, j} = ({(u_{l}^{h})}_{i + 1, j} - {(u_{l}^{h})}_{i, j}) / h_{1},

δ_{x_{2}}^{-} {(u_{l}^{h})}_{i, j} = ({(u_{l}^{h})}_{i, j} - {(u_{l}^{h})}_{i, j - 1}) / h_{2},

δ_{x_{2}}^{+} {(u_{l}^{h})}_{i, j} = ({(u_{l}^{h})}_{i, j + 1} - {(u_{l}^{h})}_{i, j}) / h_{2},

\nabla^{h} {(u_{l}^{h})}_{i, j} = {(δ_{x_{1}}^{+} {(u_{l}^{h})}_{i, j}, δ_{x_{2}}^{+} {(u_{l}^{h})}_{i, j})}^{⊤},

\nabla^{h} \cdot {(w_{l}^{h})}_{i, j} = δ_{x_{1}}^{-} {(w_{l, 1}^{h})}_{i, j} + δ_{x_{2}}^{-} {(w_{l, 2}^{h})}_{i, j},

- △^{h} {(u_{l}^{h})}_{i, j} = - \nabla^{h} \cdot (\nabla^{h} {(u_{l}^{h})}_{i, j}),

= \frac{1}{h_{1}^{2}} ({(Σ)}_{i, j}^{h} {(u_{l}^{h})}_{i, j} - {({\bar{Σ}}^{h})}_{i, j} {(u_{l}^{h})}_{i, j}),

{(Σ)}_{i, j}^{h} = 2 (1 + γ^{2}), γ = h_{1} / h_{2},

{(\bar{Σ})}_{i, j}^{h} {(u_{l}^{h})}_{i, j} = {(u_{l}^{h})}_{i + 1, j} + {(u_{l}^{h})}_{i - 1, j} + γ^{2} {(u_{l}^{h})}_{i, j + 1} + γ^{2} {(u_{l}^{h})}_{i, j - 1},

f_{1}^{h} {(u^{h})}_{i, j} = f_{1}^{h} {(u_{1}^{h}, u_{2}^{h})}_{i, j}

= (T_{i, j}^{h^{*}} - R_{i, j}^{h}) \frac{T_{i + 1, j}^{h^{*}} - T_{i - 1, j}^{h^{*}}}{2 h_{1}},

f_{2}^{h} {(u^{h})}_{i, j} = f_{2}^{h} {(u_{1}^{h}, u_{2}^{h})}_{i, j}

= (T_{i, j}^{h^{*}} - R_{i, j}^{h}) \frac{(T_{i, j + 1}^{h^{*}} - T_{i, j - 1}^{h^{*}})}{2 h_{2}},

T_{i, j}^{h^{*}} = T^{h} (i + {(u_{1}^{h})}_{i, j}, j + {(u_{2}^{h})}_{i, j}),

{(g_{l}^{h})}_{i, j} = - δ_{x_{1}}^{-} (θ_{l} {(w_{l, 1}^{h})}_{i, j} + {(λ_{l, 1}^{h})}_{i, j})

- δ_{x_{2}}^{-} (θ_{l} {(w_{l, 2}^{h})}_{i, j} + {(λ_{l, 2}^{h})}_{i, j}) .

We note that all finite difference approximations need to be adjusted at the image boundary $\partial Ω_{h}$ using the approximations of the boundary conditions

{(u_{l}^{h})}_{i, 1} = {(u_{l}^{h})}_{i, 2}, {(u_{l}^{h})}_{i, n_{2}} = {(u_{l}^{h})}_{i, n_{2} - 1},

{(u_{l}^{h})}_{1, j} = {(u_{l}^{h})}_{2, j}, {(u_{l}^{h})}_{n_{1}, j} = {(u_{l}^{h})}_{n_{1} - 1, j} .

In the following subsections the symbols ‘h’ and ‘ ${(\cdot)}_{i, j}$ ’ will sometimes drop for simplicity.

Nonlinear multigrid method for solving u-subproblem

The most difficult part of the proposed ALM is the solution of the EL equation for u -subproblem in (20). This subproblem is nonlinear and impossible to obtain a closed-form solution. It therefore requires one to approximately solve the corresponding nonlinear discrete system via some iterative method (e.g, a FP method and a Newton-type method). As is well-known, the computational costs of iterative methods are too expensive and we are motivated to develop an efficient numerical method for reducing the computational work.

In this section, we propose a nonlinear multigrid (NMG) method due to literature.²² The basic idea is to accurate the convergence of some basic iterative method on the finest grid by relying on the complementary interplay of smoothing and coarse-grid correction principles. For a more comprehensive treatment of MG methods in the area of image registration, we refer literature^{5,11,18,23–28} and references therein.

Let us denote the nonlinear discrete system in (24) using the following notation:

{\begin{matrix} N_{1}^{h} (u^{h}) = g_{1}^{h}, \\ N_{2}^{h} (u^{h}) = g_{2}^{h} \end{matrix}

(25)

where

u^{h}

denotes the exact solution. Let

{\bar{u}}^{h} = {({\bar{u}}_{1}^{h}, {\bar{u}}_{2}^{h})}^{⊤}

be an approximate solution of

u^{h}

resulting from applying a few iterations of the smoother on a fine-grid Ω _h (pre-smoothing step). Define the algebraic error as

e^{h} = u^{h} - {\bar{u}}^{h}

. Therefore, the residual equation is given by

N_{l}^{h} ({\bar{u}}^{h} + e^{h}) - N_{l}^{h} ({\bar{u}}^{h}) = g_{l}^{h} - N_{l}^{h} ({\bar{u}}^{h}) = r_{l}^{h} (l = 1, 2) .

To correct ${\bar{u}}^{h}$ numerically on Ω _h , $e^{h}$ is required to compute, but this computation is considerably and prohibitively expensive on Ω _h . Since the smoother in the pre-smoothing step removes the high frequency component of the error, one can represent the nonlinear system to the coarse grid Ω _H as given by

︸ \underset{N_{l}^{h} (u^{h})}{N_{l}^{h} ({\bar{u}}^{h} + e^{h})} = ︸ \underset{g_{l}^{h}}{r_{l}^{h} + N_{l}^{h} ({\bar{u}}^{h})}

↓

(26)

︸ \underset{N_{l}^{H} (u^{H})}{N_{l}^{H} ({\bar{u}}^{H} + e^{H})} = ︸ \underset{g_{l}^{H}}{r_{l}^{H} + N_{l}^{H} ({\bar{u}}^{H})},

where H is the index for the new cell size

H_{1} \times H_{2}

with

H_{1} \geq h_{1}

and

H_{2} \geq h_{2}

. Next, the nonlinear residual equation (27) on Ω _H has to be solved by an efficient method (the coarsest grid solver). The coarse-grid correction

e^{H} = u^{H} - {\bar{u}}^{H}

on Ω _H is then interpolated back to that on the fine grid

e^{h}

, which can now be used to update

{\bar{u}}^{h}

on Ω _h by using

{\bar{u}}_{new}^{h} =

{\bar{u}}^{h} + e^{h}

(coarse-grid correction step). Finally, the smoother is again performed with a few iterations on Ω _h to remove high frequency parts of the interpolated error (post-smoothing step). This procedure is known as a two-grid cycle, and it can be extended to a MG method with recursive application.

Below, we provide the technical details for our MG components.

The standard coarsening method is used in Ω _H by doubling the grid spacing in each space direction — i.e. $h_{l} \to 2 h_{l} = H_{l}$ for l = 1, 2.

The restriction operator $I_{h}^{H}$ and interpolation operators $I_{H}^{h}$ are used as the intergrid transfer operators and determined by the averaging and bilinear interpolation techniques — cf. literature^29–33 for more details.

The discretization coarse grid approximation method is used to compute the coarse-grid operators.

The coarsest grid solver is the LGD method.

The MG cycle is $V (ν_{1}; ν_{2})$ -cycle, where ν₁ and ν₂ denote respectively the numbers of pre- and post-smoothing steps.

The MG smoother to be discussed shortly in the next section is obtained from a coupled outer-inner iteration method using a relaxation parameter ω.

The implementation of the proposed NMG method can be summarized as follows: ${\bar{u}}^{h} \leftarrow NMGCYC ({\bar{u}}^{h}, g_{1}^{h}, g_{2}^{h}, R^{h}, T^{h}, ν_{1}, ν_{2}, θ_{1}, θ_{2}, ω, Siter)$

$1)$ If $Ω_{h} = coarsest$ grid $(| Ω_{h} | = 4 \times 4)$ , solve (25) using the LGD method and then stop.

Else continue with following step.

$2)$ Pre-smoothing:

For k = 1 to ν₁,

$[{\bar{u}}^{h}] \leftarrow Smoother ({\bar{u}}^{h}, g_{1}^{h}, g_{2}^{h}, R^{h}, T^{h}, θ_{1}, θ_{2}, ω, Siter)$

(Siter represents the maximum number of the inner iterations)

$3)$ Restriction to the coarse grid:

${\bar{u}}_{1}^{H} \leftarrow I_{h}^{H} {\bar{u}}_{1}^{h},$ ${\bar{u}}_{2}^{H} \leftarrow I_{h}^{H} {\bar{u}}_{2}^{h}, R^{H} \leftarrow I_{h}^{H} R^{h}, T^{H} \leftarrow I_{h}^{H} T^{h}$

$4)$ Set the initial solution for the coarse-grid problem: $[{\tilde{u}}_{1}^{H}, {\tilde{u}}_{2}^{H}] \leftarrow [{\bar{u}}_{1}^{H}, {\bar{u}}_{2}^{H}]$

$5)$ Compute the new right-hand side for the coarse-grid problem:

$g_{1}^{H} \leftarrow I_{h}^{H} (g_{1}^{h} - N_{1}^{h} ({\bar{u}}^{h})) + N_{1}^{H} ({\bar{u}}^{H})$ ,

$g_{2}^{H} \leftarrow I_{h}^{H} (g_{2}^{h} - N_{2}^{h} ({\bar{u}}^{h})) + N_{2}^{H} ({\bar{u}}^{H})$

$6)$ Implement the NMG method on the coarse-grid problem:

${\bar{u}}^{H} \leftarrow NMGCYC ({\bar{u}}^{H}, g_{1}^{H}, g_{2}^{H}, R^{H}, T^{H}, ν_{1}, ν_{2}, θ_{1}, θ_{2},$

$ω, Siter)$

$7)$ Add the coarse-grid corrections:

${\bar{u}}_{1}^{h} \leftarrow {\bar{u}}_{1}^{h} + I_{H}^{h} ({\bar{u}}_{1}^{H} - {\tilde{u}}_{1}^{H}),$ ${\bar{u}}_{2}^{h} \leftarrow {\bar{u}}_{2}^{h} + I_{H}^{h} ({\bar{u}}_{2}^{H} - {\tilde{u}}_{2}^{H})$

$8)$ Post-smoothing:

For k = 1 to ν₂,

$[{\bar{u}}^{h}] \leftarrow Smoother ({\bar{u}}^{h}, g_{1}^{h}, g_{2}^{h}, R^{h}, T^{h}, θ_{1}, θ_{2}, ω, Siter)$

To solve u -subproblem using (25) numerically, the proposed NMG method is stopped if the maximum number of the MG cycles $ε_{1}$ is reached or the mean of the relative residuals obtained from the discrete EL equation (25) is smaller than a small prescribed number $ε_{2} > 0$ .

Finally, the pseudo-code implementation of the proposed NMG method can be given in Algorithm 2.

Algorithm 2 The NMG algorithm for u -subproblem

\bar{u} h \leftarrow NMG ({\bar{u}}^{h}, g_{1}^{h}, g_{2}^{h}, R^{h}, T^{h}, ν_{1}, ν_{2}, θ_{1}, θ_{2}, ω, Siter, \vec{ε})

$1)$ Initialization: select $\vec{ε} = {(ε_{1}, ε_{2})}^{⊤}$ , set K = 0,

${[{\bar{u}}^{h}]}^{[K]} = {\bar{u}}^{h}$ and $RelRes = ε_{2} + 1$ .

$2)$ Repeat

$2.1)$ ${[{\bar{u}}^{h}]}^{[K + 1]} \leftarrow NMGCYC ({[{\bar{u}}^{h}]}^{[K]}, g_{1}^{h}, g_{2}^{h}, R^{h}, T^{h},$

$ν_{1}, ν_{2}, θ_{1}, θ_{2}, ω, Siter)$ .

$2.2)$ Compute

$RelRes = mean {\frac{| | g_{l}^{h} - N_{l}^{h} ({[{\bar{u}}^{h}]}^{[K + 1]}) | |}{| | g_{l}^{h} - N_{l}^{h} ({[{\bar{u}}^{h}]}^{[0]}) | |} | l = 1, \dots, 2}$ .

$2.3)$ Set $K = K + 1$ .

until a stopping rule for NMG method is met.

The MG smoother

Following literature^{5,11,18,23–27} the local Fourier analysis can be used to guarantee that there exists an efficient point-wise smoother within a MG method for solving the discrete nonlinear system (25). To obtain a high-potential point-wise smoother, this work proposes a coupled outer-inner iteration method in a FP framework.

Let ν denote the index for the outer step. We start the outer iteration by introducing the iterative scheme as follows

{\begin{array}{l} - θ_{1} △ u_{1}^{[ν + 1]} + f_{1} (u^{[ν + 1]}) = g_{1}, \\ - θ_{2} △ u_{2}^{[ν + 1]} + f_{2} (u^{[ν + 1]}) = g_{2} . \end{array}

Since this iterative scheme is fully implicit, a linearization procedure of the nonlinear term $f_{l} (u^{[ν + 1]})$ is required. As pointed out in literature^{5,11,23,26–28} it is appropriate to apply a global linearization technique with the nonlinear terms $f_{l} (u^{[ν + 1]})$ . Thus, the resulting approximation is given by

f_{l} (u^{[ν + 1]}) = f_{l} (u_{1}^{[ν + 1]}, u_{2}^{[ν + 1]})

\approx f_{l} (u_{1}^{[ν]}, u_{2}^{[ν]}) + σ_{l 1}^{[ν]} δ u_{1}^{[ν]} + σ_{l 2}^{[ν]} δ u_{2}^{[ν]}

(27)

where

σ_{l 1} (u^{[ν]}) = \partial_{u_{1}} f_{l} (u_{1}^{[ν]}, u_{2}^{[ν]})

= (\partial_{u_{l}} T (u^{[ν]})) (\partial_{u_{1}} T (u^{[ν]}))

+ (T (u^{[ν]}) - R) (\partial_{u_{1} u_{l}} T (u^{[ν]})),

σ_{l 2} (u^{[ν]}) = \partial_{u_{2}} f_{l} (u_{1}^{[ν]}, u_{2}^{[ν]})

= (\partial_{u_{l}} T (u^{[ν]})) (\partial_{u_{2}} T (u^{[ν]}))

+ (T (u^{[ν]}) - R) (\partial_{u_{2} u_{l}} T (u^{[ν]}))

To obtain a simple and fast iterative scheme, we use the approximations for $σ_{l 1}$ and $σ_{l 2}$ as introduced by literature.^11,34–36 Therefore, we have

σ_{11} (u^{[ν]}) = (\partial_{u_{1}} T_{u^{[ν]}}) (\partial_{u_{1}} T_{u^{[ν]}})

+ (T_{u^{[ν]}} - R) (\partial_{u_{1} u_{1}} T_{u^{[ν]}})

\approx (\partial_{u_{1}} T_{u^{[ν]}}) (\partial_{u_{1}} T_{u^{[ν]}}),

σ_{22} (u^{[ν]}) = (\partial_{u_{2}} T_{u^{[ν]}}) (\partial_{u_{2}} T_{u^{[ν]}})

+ (T_{u^{[ν]}} - R) (\partial_{u_{2} u_{l}} T_{u^{[ν]}})

\approx (\partial_{u_{2}} T_{u^{[ν]}}) (\partial_{u_{2}} T_{u^{[ν]}}),

σ_{12} (u^{[ν]}) = σ_{21} (u^{[ν]}) = 0.

These approximations leads to the following linearized system

{\begin{array}{l} - θ_{1} △ u_{1}^{[ν + 1]} + σ_{11} (u^{[ν]}) u_{1}^{[ν + 1]} = \\ g_{1} - f_{1} (u^{[ν]}) + σ_{11} (u^{[ν]}) u_{1}^{[ν]}, \\ - θ_{2} △ u_{2}^{[ν + 1]} + σ_{22} (u^{[ν]}) u_{2}^{[ν + 1]} = \\ g_{2} - f_{2} (u^{[ν]}) + σ_{22} (u^{[ν]}) u_{2}^{[ν]} \end{array}

(28)

Next, we apply the finite difference discretization as discussed in Section 3.1 with (28) and solve the associated linear system for each outer step ν by the Gauss-Seidel (GS) method as the inner iteration. The kth step of the GS method defined at a grid point $(x_{1_{i}}, x_{2_{j}})$ is then given by

{(u^{[ν + 1]})}_{i, j}^{[k + 1]} = {(v)}_{i, j}^{[k + 1]}

(29)

where

\begin{array}{l} (v)_{i, j}^{[k + 1]} = (\begin{array}{l} F_{1}^{[ν]} \\ F_{2}^{[ν]} \end{array}), \\ F_{1}^{[ν]} = (\frac{1}{(θ_{1} / h_{1}^{2}) {(Σ)}_{i, j} + {(σ_{11} (u^{[ν]}))}_{i, j}}) \\ \times [{(g_{1})}_{i, j} - f_{1} {(u^{[ν]})}_{i, j} + {(σ_{11} (u^{[ν]}))}_{i, j} {(u_{1}^{[ν]})}_{i, j} \\ + (θ_{1} / h_{1}^{2}) {(\bar{Σ})}_{i, j} {(u_{1}^{[ν + 1]})}_{i, j}^{[k + 1 / 2]}] \end{array}

F_{2}^{[ν]} = (\frac{1}{(θ_{2} / h_{1}^{2}) {(Σ)}_{i, j} + {(σ_{22} (u^{[ν]}))}_{i, j}}) \times

= [{(g_{2})}_{i, j} - f_{2} {(u^{[ν]})}_{i, j} + {(σ_{22} (u^{[ν]}))}_{i, j} {(u_{2}^{[ν]})}_{i, j}

+ (θ_{2} / h_{1}^{2}) {(\bar{Σ})}_{i, j} {(u_{2}^{[ν + 1]})}_{i, j}^{[k + 1 / 2]}]

and

{(\bar{Σ})}_{i, j} {(u_{l}^{[ν + 1]})}_{i, j}^{[k + 1 / 2]} = {(u_{l}^{[ν + 1]})}_{i + 1, j}^{[k]} + {(u_{l}^{[ν + 1]})}_{i - 1, j}^{[k + 1]}

+ γ^{2} {(u_{l}^{[ν + 1]})}_{i, j + 1}^{[k]}

+ γ^{2} {(u_{l}^{[ν + 1]})}_{i, j - 1}^{[k + 1]}

Here the superscripts k, $k + 1 / 2$ , and k + 1 denote the current, intermediate and new approximations computed by the GS method.

In order to obtain more efficiency, one introduces a relaxation parameter $ω \in (0, 2)$ and iterates the successive over relaxation (SOR) steps by

{(u^{[ν + 1]})}_{i, j}^{[k + 1]} = (1 - ω) {(u^{[ν + 1]})}_{i, j}^{[k]} + ω {(v)}_{i, j}^{[k + 1]}

(30)

Finally, the implementation of the proposed MG smoother (30) on Ω _h is summarized in Algorithm 3.

Algorithm 3 The MG smoother for the proposed NMG algorithm (Algorithm 2)

{\bar{u}}^{h} \leftarrow Smoother ({\bar{u}}^{h}, g_{1}^{h}, g_{2}^{h}, R^{h}, T^{h}, θ_{1}, θ_{2}, ω, Siter)

$1)$ ] Use input parameters to compute ${(v)}_{i, j}$ for all $1 \leq i \leq n_{1}$ and $1 \leq j \leq n_{2}$ .

$2)$ ] Perform SOR steps for $k = 0, 1, 2, \dots, Siter$

$2.1)$ ] Update ${({\bar{u}}^{h})}_{i, j}^{[k + 1]}$ by (30) for all $1 \leq i \leq n_{1}$ and $1 \leq j \leq n_{2}$ .

Note that we applied the so-called local Fourier analysis (LFA) to test the smoothing performance of the proposed MG smoother with different values of ω (the relaxation parameter) and Siter (the maximum number of inner iterations) for the four representative registration problems given by Problems 1–4 (as shown respectively in Figure 1. Our results shows that the proposed FP smoother with $ω = 1.85$ and Siter = 2 is recommended; see its performance in leading the proposed ALM method to fast convergence from Tables 3 and 4. The theoretical details for applying LFA to analyze the proposed smoother are straightforward and similar to our previous work in (11)

Figure 1.

Four registration problems of the real medical images. (left column) Reference images; (right column) Template images; (top row) Problem 1; (2nd row) Problem 2; (3rd row) Problem 3; (bottom row) Problem 4.

The closed form expressions for $w_{1}, w_{2}, λ_{1}$ , and $λ_{2}$

Based on the formulation (23), we can get the closed form expressions for $w_{1}$ and $w_{2}$ as given by

{(w_{l})}_{i, j} = \max (| {(S_{l})}_{i, j} | - \frac{α}{θ_{l}}, 0) \frac{{(S_{l})}_{i, j}}{| {(S_{l})}_{i, j} |},

(31)

where

S_{l} = {(δ_{x_{1}}^{+} {(u_{l})}_{i, j} - \frac{{(λ_{l, 1})}_{i, j}}{θ_{l}}, δ_{x_{2}}^{+} {(u_{l})}_{i, j} - \frac{{(λ_{l, 2})}_{i, j}}{θ_{l}})}^{⊤}

for l = 1, 2.

Similarly, based on the formulations in (11) and (12), we may update all the Lagrangian multipliers by

{(λ_{1}^{(k + 1)})}_{i, j} = {(λ_{1}^{(k)})}_{i, j} + θ_{1} ({(w_{1})}_{i, j}^{(k + 1)} - \nabla^{+} {(u_{1})}_{i, j}^{(k + 1)})

(32)

{(λ_{2}^{(k + 1)})}_{i, j} = {(λ_{2}^{(k)})}_{i, j} + θ_{2} ({(w_{2})}_{i, j}^{(k + 1)} - \nabla^{+} {(u_{2})}_{i, j}^{(k + 1)})

(33)

Numerical experiments

In this section, we present a number of numerical experiments to

compare the overall performance of the proposed ALM with two related numerical solutions, which are the LGD method (5) and the NMG method by Chumchob (5)

assess the accuracy and efficiency of the proposed ALM with regard to parameter changes.

We note first that four registration problems consisting of four real medical images to be denoted as Problems 1–4 were selected for the experiments, as shown respectively in Figure 1. Second, we used $n_{1} = n_{2} = n$ (i.e. we assume the grid spacing $h_{1} = h_{2} = h = 1 / n$ ) in all numerical tests. Third, the bilinear interpolation was employed to compute the transformed template image $T_{u}$ once the deformations are found in all cases. We also note that all numerical algorithms for the proposed ALM, LGD and NMG methods were started with $u^{[0]} = 0$ and implemented under MATLAB R2018a and run on a machine configured with Intel(R) Core(TM) i7 Quad-core 4.2 GHz and 32GB of RAM.

Performance comparison with the other two methods

In this test, the performance of three different methods for TV-based image registration model is compared. We apply the relative SSD (a qualitative measure in the accuracy), time per iteration (in seconds), total CPU time (in seconds), and total iterations to perform our evaluation on the four registration problems shown in Figure 1.

For the LGD method, we chose $τ = 10^{- 5}$ and $β = 10^{- 6}$ and implemented a direct solver using MATLAB’s $\$ -operator (also known as mldivide operator) for solving each linear system of equations arising from applying the FP method with the nonlinear discrete EL system.

For the NMG method, we implemented with the FAS-NMG method of Chumchob (5) The details of our numerical implementation can be summarized as follows. We solved the nonlinear discrete EL system with $β = 10^{- 6}$ in a V-cycle framework, where the coarsest grid is 4 × 4. The LGD method was used to solve the coarsest problem with $τ = 10^{- 3}$ . We used the standard coarsening method and standard interpolation technique for transferring data between grids. We applied 15 pre- and 15 post-smoothing steps with a so-called FP-SOR smoother.

For the proposed ALM, we chose the penalty parameters $θ_{1} = θ_{2} = 1 / 10$ . We used the predefined numbers $ϵ_{1} = 10^{- 2}, ϵ_{2} = 10^{- 4}, ϵ_{3} = ϵ_{4} = 10^{- 8}$ and $ϵ_{5} = 200$ with Algorithm 1. The predefined numbers $ε_{1} = 20$ and $ε_{2} = 10^{- 6}$ were applied by the proposed NMG method in Algorithm 2 to solve u -subproblem. In all registration problems, the proposed NMG method with the MG smoother in Algorithm 3 was performed with the MG parameters, $ν_{1} = ν_{2} = 5$ (5 pre- and 5 post-smoothing steps), $ω = 1.85$ (the relaxation parameter), and $Siter = 2$ (2 iterations of SOR steps).

All numerical methods in this test used the same regularization parameter $α = 1 / 50$ with the grid spacing h = 1/256. Each method was stopped when the relative solution, $\max {\frac{∥ u_{1}^{(m)} - u_{1}^{(m - 1)} ∥}{∥ u_{1}^{(m)} ∥}, \frac{∥ u_{2}^{(m)} - u_{2}^{(m - 1)} ∥}{∥ u_{2}^{(m)} ∥}}$ , was brought below $10^{- 4}$ or the number of iterations increases beyond a predefined number. Here we used the predefined numbers $m_{max}^{LGD} = 200$ for the LGD method, $m_{max}^{NMG} = 50$ for the NMG method, and $m_{max}^{ALM} = 200$ for the proposed ALM.

From Table 1, one can see that only the proposed ALM converges, whereas the LGD and NMG methods are unable to converge in a reasonable number of iterations. As expected, we have found in this case that the small value of β has a significant effect on the accuracy of the registered images and the convergence of the LGD and NMG method. Moreover, Figure 2 shows that for each registration problem the proposed ALM yields the best value of the relative SSD. It is important to note that the smaller the value of the relative SSD is, the accurate registered image is achieved. This evidence ensures that the registration results by the proposed ALM are more reliable than the LGD and NMG methods. Particularly, Figure 3 shows that the proposed ALM is computationally efficient than the other two methods in delivering the high quality of the registered images. The proposed ALM takes only a few iterations (say 5 iterations) to drop the relative SSD below 0.20, which means that the dissimilarities between the reference and registered images have been reduced more than 80% within the first 5 iterations. This is a remarkable result to conclude that the computational performance of the proposed ALM in solving TV-based image registration is much more efficient than those of the other two methods.

Table 1.

Comparison of the relative SSD, time per iteration, total CPU time (s), and total iterations by three different numerical methods on the four registration problems in Figure 1 with $α = \frac{1}{50}$ and $h = \frac{1}{256}$ . * indicates that maximum number of iterations reached without convergence.

	Registration problem
Method	1	2	3	4
LGD
Relative SSD	0.8027	0.7853	0.7841	0.8250
Time per iteration (s)	0.2467	0.2458	0.2457	0.2463
Total CPU time (s)	49.3423	49.1558	49.1301	49.2524
Total iterations	$200^{*}$	$200^{*}$	$200^{*}$	$200^{*}$
NMG
Relative SSD	0.2465	0.1670	0.0710	0.0852
Time per iteration (s)	0.4941	0.4930	0.4814	0.4914
Total CPU time (s)	24.7045	24.6477	24.0701	24.5724
Total iterations	$50^{*}$	$50^{*}$	$50^{*}$	$50^{*}$
Proposed ALM
Relative SSD	$0.1836$	$0.1476$	$0.0664$	$0.0681$
Time per iteration (s)	0.5970	0.6893	0.7105	0.7138
Total CPU time (s)	28.4776	40.6674	37.6567	38.5466
Total iterations	41	59	53	54

Figure 2.

Qualitative comparison of registered images by three different methods on the four registration problems in Figure 1. (left) Registered images by LGD; (middle) registered images by NMG; (right) registered images by the proposed ALM. ‘RelSSD’ means the relative SSD defined in (13) to represent the dissimilarities between R and $T_{u}$ . Note that the small the value of RelSSD is, the accurate registered image is obtained.

Figure 3.

History of the relative SSD by three different methods in solving the four registration problems in Figure 1.

In terms of computation times per iteration, we can see that the proposed ALM requires to solve two sub-problems for u and w in each step. Thus, in accordance with our measurements, we expect that more computation time per iteration is required for the proposed ALM than those of the other two methods. As can be seen from Table 1, our numerical experiments reflect this expectation.

In terms of total CPU times, Table 1 illustrates the LGD method is the slowest method, while the NMG method is slightly faster than the proposed ALM. However, both LGD and NMG method are unable to converge. Thus we can conclude that the LGD and NMG methods are less computationally efficient than the proposed ALM.

In Figure 4 we present the constructed deformations by the proposed ALM for all registration problems. In Figure 5, we show the constructed surfaces of both components of the deformations in Figure 4. In Figure 6, we present the middle slices of the constructed surfaces of the components of the deformations in Figure 5. As expected from the use of TV regularizer, the visual inspection in Figure 4 shows that the proposed ALM delivers the visually appealing results in preserving the discontinuities of the constructed deformation as shown in the corresponding close-up regions. We can also see in Figure 4 that the constructed deformation yields multiple motions in the image to be registered and the motion discontinuities can be observed at the boundaries of the local regions. We see further that the corresponding deformation in each registration problem is neither smooth nor reflecting a homogeneous motion. For each local region, the motion tends to change smoothly, and the gradients of the constructed deformation in this area are small. At the boundaries of the local regions, the gradients of the constructed deformation are large, and the discontinuities of the constructed deformation can be observed. Moreover, we can see from Figures 5 and 6 that the constructed surfaces of both components of the deformations are non-smooth and their middle slices are almost piecewise constant for all registration problems.

Figure 4.

Constructed deformations by the proposed ALM for the four registration problems in Figure 1. Note that multiple motions in the corresponding image and motion discontinuities at the boundaries of the local regions determined by the constructed deformation can be observed in the corresponding close-up regions.

Figure 5.

Constructed surfaces of both components of the constructed deformations by the proposed ALM for the four registration problems in Figure 1. (left) the first component u₁; (right) the second component u₂; (top row) Problem 1; (2nd row) Problem 2; (3rd row) Problem 3; (bottom row) Problem 4.

Figure 6.

The middle slices of the corresponding surfaces of both components of the constructed deformations shown in Figure 5 by the proposed ALM.

To summarize we have successfully developed the efficient and effective numerical method for TV-based image registration model. Our registration results on the real medical applications shown in Figure 1 demonstrate that the proposed ALM is more computationally efficient and effective than the other two methods. The most attractive features of our proposed method is that it is able to produce the high quality of the registered images in a reasonable number of iterations while it satisfies the requirement in constructing the deformations by the TV regularizer.

Performance tests with the regard to parameter changes

We now present numerical results from several test cases, to assess the accuracy and efficiency of our proposed numerical techniques with the regard to parameter changes.

h-Independence test

One of the key properties of MG techniques is that their convergence does not depend on the number of grid points. Thus, in this test we designed our numerical experiments to investigate this property with the proposed ALM in Algorithm 1, and to back up our proposed NMG method in Algorithms 2.

We implemented the proposed ALM with the same regularization parameter $α = \frac{1}{250}$ and the penalty parameters θ₁ and θ₂ as shown in Table 2. In order to stop our proposed ALM, we used the same predefined numbers used in Section 4.1. We also used the same MG parameters to solve u -subproblem as given in the previous section.

Table 2.

Registration results for Problems 1–4 in Figure 1 by the proposed ALM in Algorithm 1 with $α = \frac{1}{250}$ and different values of grid spacing h. Note that ‘Ratio’ means the ratio of the total CPU times in increasing both image dimensions by a factor of 2.

Registration	Grid spacing (h)
problem	$\frac{1}{256}$	$\frac{1}{512}$	$\frac{1}{1024}$	$\frac{1}{2048}$
1 ( $θ_{1} = θ_{2} = 1$ )
Relative SSD	0.0968	0.0958	0.0951	0.0951
Total iterations	22	22	21	21
Total CPU time (s)	22.4489	74.1501	308.4895	$1, 308.3476$
Ratio	–	3.3031	4.1603	4.2411
Average MG cycles	6	6	6	6
2 ( $θ_{1} = θ_{2} = 1$ )
Relative SSD	0.1028	0.1027	0.1025	0.1002
Total iterations	16	16	15	15
Total CPU time (s)	16.1900	58.8500	257.0531	$1, 117.0938$
Ratio	–	3.6350	4.3679	4.4294
Average MG cycles	7	7	7	7
3 ( $θ_{1} = θ_{2} = \frac{1}{4}$ )
Relative SSD	0.0359	0.0359	0.0358	0.0357
Total iterations	36	34	33	33
Total CPU time (s)	38.8623	137.9971	577.3338	$2, 578.2640$
Ratio	–	3.5509	4.1837	4.4658
Average MG cycles	7	7	7	7
4 ( $θ_{1} = θ_{2} = \frac{1}{4}$ )
Relative SSD	0.0377	0.0375	0.0374	0.0373
Total iterations	31	31	30	30
Total CPU time (s)	30.9928	116.8559	500.9036	$2, 238.7849$
Ratio	–	3.7704	4.2865	4.4695
Average MG cycles	6	6	6	6

Table 3.

Registration results for Problems 1–4 in Figure 1 by the proposed ALM in Algorithm 1 with $h = \frac{1}{256}$ and different values of the regularization parameter α.

Registration	Regularization parameter (α)
problem	$\frac{1}{250}$	$\frac{1}{200}$	$\frac{1}{150}$	$\frac{1}{100}$	$\frac{1}{50}$	$\frac{1}{25}$
1 ( $θ_{1} = θ_{2} = 1$ )
Relative SSD	0.0968	0.1015	0.1080	0.1313	0.2099	0.3466
Total iterations	20	21	25	30	33	49
Average MG cycles	6	6	6	6	6	6
2 ( $θ_{1} = θ_{2} = 1$ )
Relative SSD	0.1028	0.1028	0.1069	0.1174	0.1675	0.2104
Total iterations	16	18	20	27	37	40
Average MG cycles	7	6	6	6	6	6
3 ( $θ_{1} = θ_{2} = \frac{1}{4}$ )
Relative SSD	0.0359	0.0385	0.0436	0.0536	0.0958	0.2251
Total iterations	36	44	49	65	90	125
Average MG cycles	7	7	6	6	5	5
4 ( $θ_{1} = θ_{2} = \frac{1}{4}$ )
Relative SSD	0.0377	0.0406	0.0449	0.0548	0.0863	0.1909
Total iterations	31	34	49	57	102	139
Average MG cycles	6	6	5	5	5	5

Table 4.

Registration results for Problems 1–4 in Figure 1 by the proposed ALM in Algorithm 1 with $h = \frac{1}{256}, α = \frac{1}{250}$ and different values of the penalty parameter θ. Recall that * indicates that maximum number of MG cycles reached without convergence.

Penalty	Registration problem
parameter	1	2	3	4
θ = 10
Relative SSD	0.8128	0.8128	0.4923	0.5056
Total iterations	3	3	3	3
Average MG cycles	$20^{*}$	15	18	18
θ = 8
Relative SSD	0.7721	0.7774	0.4271	0.4401
Total iterations	3	3	3	3
Average MG cycles	$20^{*}$	19	15	15
θ = 6
Relative SSD	0.7098	0.7138	0.3476	0.3595
Total iterations	3	3	3	3
Average MG cycles	$20^{*}$	$20^{*}$	14	14
θ = 4
Relative SSD	0.6033	0.6029	0.2615	0.2710
Total iterations	3	3	3	3
Average MG cycles	$20^{*}$	$20^{*}$	19	19
θ = 2
Relative SSD	0.1935	0.2048	0.2187	0.2467
Total iterations	8	7	10	9
Average MG cycles	12	14	$20^{*}$	$20^{*}$
θ = 1
Relative SSD	0.0968	0.1028	0.2057	0.2443
Total iterations	20	16	17	19
Average MG cycles	6	7	$20^{*}$	$20^{*}$
$θ = \frac{1}{2}$
Relative SSD	0.0693	0.0696	0.0670	0.0441
Total iterations	39	31	25	19
Average MG cycles	5	5	$20^{*}$	$20^{*}$
$θ = \frac{1}{4}$
Relative SSD	0.0617	0.0586	0.0359	0.0377
Total iterations	60	44	36	31
Average MG cycles	4	5	7	6
$θ = \frac{1}{6}$
Relative SSD	0.0600	0.0549	0.0338	0.0348
Total iterations	75	60	46	52
Average MG cycles	4	4	5	5
$θ = \frac{1}{8}$
Relative SSD	0.0584	0.0536	0.0328	0.0340
Total iterations	89	100	63	59
Average MG cycles	4	4	5	4
$θ = \frac{1}{10}$
Relative SSD	0.0579	0.0535	0.0326	0.0336
Total iterations	108	191	75	62
Average MG cycles	4	4	4	4

In the numerical results shown in Table 2, one can see five quantities: the relative SSD, the total iterations, the total CPU times (in seconds), the ratio of the total CPU times in increasing both image dimensions by a factor of 2, and the number of average MG steps used to solve u -subproblem by Algorithm 2 with different values of grid spacing h.

As expected from a numerical technique using a MG framework, Table 2 shows that Algorithm 1 not only converges within a few iterations, but it also provides the accurate registration results because the dissimilarities between the reference and registered images given have been reduced more than 90% for Problem 1, 89% for Problem 2, and 96% for Problems 3 and 4. Moreover, the NMG method in Algorithm 2 can reduce the mean of the relative residuals to $ε_{2} = 10^{- 6}$ within a few MG steps for different values of h. In addition, we see an increase in the total CPU time by a factor of approximately 4.5 when both image dimensions is increased by a factor of 2. This shows that the proposed ALM has an optimal efficiency resulting from the NMG method proposed in Algorithm $2$ , which should be of order $O (N \log (N))$ where $N = n^{2}$ . Then we can conclude that the proposed ALM has no effect on the convergence for different values of h.

α-Dependence test

Next we evaluate to show how our proposed ALM in Algorithm 1 is affected with varying the regularization parameter α.

To this end, we performed the proposed ALM on the four registration problems in Figure 1 with the same grid spacing $h = \frac{1}{256}$ and the penalty parameters θ₁ and θ₂ as given in Table 2. Our proposed ALM in Algorithm 1 was stopped with the same predefined numbers used in Section 4.1. We also applied the same predefined number and MG parameters as shown in Section 4.1 with Algorithms 2 and 3 to solve u -subproblem.

Table 3 presents the numerical results by the proposed ALM in Algorithm 1 with different values of the regularization parameter α. It contains three quantities: the relative SSD, the total iterations, and the number of average MG steps used to solve u -subproblem by Algorithm 2.

As presented in Table 3, decreasing the values of α has significant effects on the accuracy of the registered images and the convergence of the proposed ALM in Algorithm 1, whereas the value of α has only a small effect on the convergence of the proposed NMG method in Algorithm 2. We can see that large α is not needed as small ones give better registration results, typically $α = \frac{1}{250}, \frac{1}{200}, \frac{1}{150}$ , and $\frac{1}{100}$ . It is important to note that the process to select the optimal value of α is a separate but important issue because it is in general unknown a priori and it significantly affects on the accuracy of registered images and the ALM performance.

In order to find a suitable α automatically, the ‘cooling’ (‘continuation’) process suggested in literature^{11,35,37–39} is recommended for real applications. The basic idea is to start with a high initial value of α and then slowly reduce α such that the obtained solution can be used to be an excellent starting point for the next in order to decrease $J$ in (8).

$(θ_{1}, θ_{2})$ -dependence test

We now present numerical results from several test cases to evaluate the registration performance of our proposed ALM with different values of θ₁ and θ₂.

We performed the proposed ALM in Algorithm 1 on the four registration problems shown in Figure 1 using $h = \frac{1}{256}$ and $α = \frac{1}{250}$ . The predefined numbers and the MG parameters for Algorithms 1–3 are the same values used in Section 4.1.

We note that we take $θ_{1} = θ_{2} = θ$ in this numerical test to evaluate the performance of the proposed ALM. Table 4 presents the registration results for the four registration problems by the proposed ALM in Algorithm 1 with different values of θ. Three qualities for the performance tests are the relative SSD, the total iterations, and the number of average MG steps used to solve u -subproblem by Algorithm 2.

One sees from Table 4 that as θ is decreased from 10 to $\frac{1}{10}$ the dissimilarities between the given images decrease (the accuracy of the registered image increase), whereas the cost of the proposed ALM increases. We observe that the smaller θ is the better the proposed NMG method performs and an increase in the number of smoothing steps is required to achieve fast MG convergence for large θ. Next, we see that decreasing the value of θ leads to the best registration result at $θ = \frac{1}{10}$ . We see further that large values of θ are not required as small ones deliver better registration performance. Typically, θ should be between $\frac{1}{10}$ and $\frac{1}{2}$ to give the accurate registration results in a reasonable number of iterations.

Concluding remarks

In this paper, we first explained how standard methods solve the TV-based image registration model. Next, we discussed how we are needed to develop a new method. In order to efficiently solve the model, we therefore proposed to use augmented Lagrangian method. We separated the associated minimization problem into two subproblems. As a result, the first subproblem is a nonlinear problem and impossible to obtain exact solution, whereas the other one has a closed-form solution. Next, we developed an efficient NMG method to solve the associated discrete nonlinear system. In order to assess the efficiency and effectiveness of our new method, we tested, using four registration problems of real medical images, how our new method performs. We found in our first numerical test that the registered images by different methods are not identical. Therefore, it can be concluded that different numerical algorithms for TV-based image registration model have a significant effect on the accuracy of registered images. We also found that our new method outperforms the existing methods. Moreover we found by the performance comparison in this test that our new method is able to deliver the accurate registration results with the desired properties of the constructed deformations in a reasonable number of iterations. Next, we found in the second test that there are no effect on the convergence of the new method for different numbers of grid points. Moreover, we observed from the third and last tests that the choice of the regularization and penalty parameters is important for the quality of the registered images and the computational performance of the new method. The outlines in selecting these two parameters to obtain the accurate registration results were discussed. Future work will extend the proposed method to high-order variational models for image registration.

Footnotes

Acknowledgements

The authors would like to thank the anonymous reviewers for a number of useful suggestions to improve the quality of the paper

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research,authorship,and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research,authorship,and/or publication of this article: The first author work was partially supported by the Thailand Research Fund Grant ♯ MRG6080169.

ORCID iD

Noppadol Chumchob

References

Maintz

Viergever

A survey of medical image registration. Med Image Anal 1998; 2: 1–36.

Hajnal

Hill

Hawkes

Medical image registration. The Biomedical Engineering Series. Boca Raton, FL: CRC Press, 2001.

Mani

Arivazhagan

Survey of medical image registration. J Biomed Eng Technol 2013; 1: 8–25.

Sotiras

Davatzikos

Paragios

Deformable medical image registration: a survey. IEEE Trans Med Imaging 2013; 32: 1153–1190.

Chumchob

Vectorial total variation-based regularization for variational image registration. IEEE Trans Image Process 2013; 22: 4551–4559.

F-Schauf

Henn

Hömke

, et al. 2006. Total variation based image registration. In: Tai X-C, Lie K-A, Chan TF, et al. (eds) Proceedings of the international conference on PDE-Based image processing and related inverse problems series: Mathematics and visualization. Berlin: Springer Verlag, pp.305–323.

F-Schauf

Henn

Witsch

Multigrid based total variation image registration. Comput Visual Sci 2008; 11: 101–113.

Broit

Optimal registration of deformed images . PhD thesis, University of Pennsylvania, USA, 1981.

Fischer

Modersitzki

Fast diffusion registration. Contemp Math 2002; 313: 117–129.

10.

Fischer

Modersitzki

Curvature based image registration. J Math Imag Vision 2003; 18: 81–85.

11.

Chumchob

Chen

A robust multigrid approach for variational image registration models. J Comput Appl Math 2011; 236: 653–674.

12.

Acar

Vogel

Analysis of bounded variation penalty method. Inverse Probl 1994; 10: 1217–1229.

13.

Chan

Chen

On a nonlinear multigrid algorithm with primal relaxation for the image total variation minimization. Numer Algor 2006; 41: 387–411.

14.

Rudin

Osher

Fatemi

Nonlinear total variation based noise removal algorithms. Physica D 1992; 60: 259–268.

15.

Chan

Mulet

On the convergence of the lagged diffusivity fixed point method in total variation image restoration. SIAM J Numer Anal 1999; 36: 354–367.

16.

Vogel

Computational methods for inverse problems. Philadelphia: SIAM Publications, 2002.

17.

Vogel

Oman

Iterative methods for total variation denoising. SIAM J Sci Comput 1996; 17: 227–238.

18.

Chumchob

Chen

A variational approach for discontinuity-preserving image registration. East-West J Math Special 2010; 2010: 266–282.

19.

Goldstein

Osher

The split bregman method for l1-regularized problems. SIAM J Imag Sci 2009; 2: 323–343.

20.

Wang

Yang

Yin

, et al. A new alternating minimization algorithm for total variation image reconstruction. SIAM J Imaging Sci 2008; 1: 248–272.

21.

Wenqi

Jinming

Zhaowen

, et al. Implementation of high-order variational models made easy for image processing. Math Methods Appl Sci 2016; 39: 4208–4233.

22.

Brandt

Multi-level adaptive solutions to BVPs.

Math Comp 1977; 31: 333–390.

23.

Chumchob

Chen

Brito

A fourth order variational image registration model and its fast multigrid algorithm. Multiscale Model Simul 2011; 9: 89–128.

24.

Chumchob

Chen

An improved variational image registration model and a fast algorithm for its numerical approximation. Numer Meth Partial Differential Eq 2012; 28: 1966–1995.

25.

Chumchob

Chen

Brito

A new variational model for removal of combined additive and multiplicative noise and a fast algorithm for its numerical approximation. Int J Comput Math 2013; 90: 140–161.

26.

Jewprasert

Chumchob

Chantrapornchai

A fourth-order compact finite difference scheme for higher-order

PDE

-based image registration. East Asian J Appl Math 2016; 5: 361–386.

27.

Chumchob

Jewprasert

Chantrapornchai

Multigrid solution of the nonlinear

PDE

s arising in elastic image registration with application to a group of monomodal images. ScienceAsia 2016; 42: 415–422.

28.

Thompson

Chen

An effective diffeomorphic model and its fast multigrid algorithm for registration of lung

CT

images. Comput Meth Appl Math 2020; 20: 141–168.

29.

Briggs

Henson

McCormick

A multigrid tutorial. 2nd ed. Philadelphia, USA: SIAM Publications, 2000.

30.

Hackbusch

Multi-grid methods and applications. Berlin, Heidelberg, New York Springer-Verlag, 1985.

31.

31. Trottenberg

Oosterlee

Schuller

Multigrid. Cambridge: Academic Press, 2001.

32.

Wesseling

Multigrid methods. Philadelphia, PA

Edwards

, 2004

33.

Wienands

Joppich

Practical Fourier analysis for multigrid method. USA: Chapman & Hall/CRC, 2005.

34.

Modersitzki

Numerical methods for image registration. Oxford, 2004.

35.

Haber

Modersitzki

A multilevel method for image registration. SIAM J Sci Comput 2006; 27: 1594–1607.

36.

Köstler

Ruhnau

Wienands

Multigrid solution of the optical flow system using a combined diffusion- and curvature-based regularizer. Numer Linear Algebra Appl 2008; 15: 201–218.

37.

Henn

Witsch

Image registration based on multiscale energy information. Multiscale Model Simul 2005; 4: 584–609.

38.

Lelièvre

Oldenburg

Magnetic forward modelling and inversion for high susceptibility. Geophys J Inter 2006; 166: 76–90.

39.

Newman

Hoversten

GM.

Solution strategies for two- and three-dimensional electromagnetic inverse problems. Inverse Problems 2000; 16: 1357–1375.

An augmented Lagrangian method for solving total variation (TV)-based image registration model

Abstract

Keywords

Introduction

The proposed ALM

Numerical implementation

Finite difference discretization

Nonlinear multigrid method for solving u-subproblem

The MG smoother

The closed form expressions for w 1 , w 2 , λ 1 , and λ 2

Numerical experiments

Performance comparison with the other two methods

Performance tests with the regard to parameter changes

h-Independence test

α-Dependence test

( θ 1 , θ 2 ) -dependence test

Concluding remarks

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

ORCID iD

References

The closed form expressions for $w_{1}, w_{2}, λ_{1}$ , and $λ_{2}$

$(θ_{1}, θ_{2})$ -dependence test