A variational joint segmentation and registration framework for multimodal images

Abstract

Image segmentation and registration are closely related image processing techniques and often required as simultaneous tasks. In this work, we introduce an optimization-based approach to a joint registration and segmentation model for multimodal images deformation. The model combines an active contour variational term with mutual information (MI) smoothing fitting term and solves in this way the difficulties of simultaneously performed segmentation and registration models for multimodal images. This combination takes into account the image structure boundaries and the movement of the objects, leading in this way to a robust dynamic scheme that links the object boundaries information that changes over time. Comparison of our model with state of art shows that our method leads to more consistent registrations and accurate results.

Introduction

In the last decades, the development of imaging devices significantly increased image data collections and consecutively the need for human interpretation. To reduce labour hours of specialists new image processing techniques are developed. Two main image processing techniques are image registration and segmentation. Image registration provides an understanding of the data behavior in two or more different scenarios taken at different times. Meanwhile, image segmentation is the process of partitioning objects or features considering the characteristics of each pixel. In a wide range of problems, such as comparison of images taken at different time or modalities, image registration depends on image segmentation and vice-versa. Erdt et al. ¹ indicates that more than 20% of scientific research in medical imaging area requires a combined registration and segmentation scheme. The combined registration and segmentation models can be divided into two different categories: simultaneous registration and segmentation or joint segmentation and registration which lastly has been referred to as regmentation. In the last years, regmentation functional models have been introduced on monomodal imaging frameworks. These models are variational rigid registration based models,^2,3 variational nonrigid registration based models,^4,5 or atlas based models. ⁶ Recently, Ibrahim-Rada-Chen ⁷ present a linear curvature regmentation variational based model which gets profit on the global and local deformation and easily can cope with images which contain more than one object.

In difference with monomodal images, multimodal image processing requires spatial correspondence estimation or extract certain information from images with different modalities (or protocols). Such example can be from imaging using CT and MRI protocols, which brings different image perspective into the medical field of clinical and pre-clinical diagnostics. Even though assessing image similarity becomes challenging for multimodal images, multimodal images are desired as they bring different image perspectives and properties. In difference with monomodal images registration, where simple image similarity measures between image pairs can provide a good spatial transformation, dealing with multimodalities is substantially more difficult. Registering multimodal images, which are acquired through different mechanisms and have distinct modalities, involves the alignment of both images in terms of shapes and salient components while preserving the modality of one given image. For a comprehensive overview of registration techniques in a systematic manner, the work of Sotiras et al. ⁸ is referred. By adding the active contour without edges model into the curvature joint registration based on MI proposed by Modersitzki, ⁹ a modified joint image segmentation and registration delivers results which are shown in Figure 1. The same images are processed with the new proposed joint regmentation model and the results shown in the last section indicate that the combination of segmentation and registration tasks into a single framework is relevant as avoids the errors produced from simultaneous tasks. The existing registration models for mono modal-images can be classified as unsupervised and supervised methods. Among different registration similarity measure, the most popular one is mutual information (MI), introduced by Viola et al. ¹⁰ Recently, learning-based approaches measuring image similarity have been proposed as well. Techniques involving KL-divergence, ¹¹ max-margin structured learning, ¹² boosting, ¹³ or deep learning ¹⁴ have been investigated. However, their large training sets make them limited by the lack of generalization to previously unseen object classes and the registration performance depends on numerical optimization of the optimal registration parameters. There are a few unsupervised works to solve the joint segmentation and registration problem. The work proposed by Wang and Vemuri ¹⁵ introduces a registration driven by segmentation of a reference image without varying degrees of non-rigidity. The model applies cross cumulative residual entropy as a distance measure ¹⁵ and piecewise constant Chan-Vese model for segmentation. ¹⁶ Although, the algorithm can accommodate image pairs having very distinct intensity distributions but fails in deforming large objects. Droske et. al ¹⁷ presented variational method based on a local energy density which uses a generalized motion of image morphology for multimodal image registration. A similar idea was lastly proposed by Aganj and Fischl ¹⁸ where the segmentation score was applied to register two multimodal brain images. This model is not designed for general image deformation as the model applies rigid image transformation. Recently, a joint segmentation and registration model has been proposed by Debroux and Le Guyader, ¹⁹ which is formulated based on the combination of nonlocal total variation and nonlocal shape descriptors. The fitting term of this model uses a weighted sum of squared difference distance measure, similar to ⁷ and the model depends on the many parameters such as window size, patch size, number of neighbors pixels, etc. To overcome the difficulties mentioned above we propose a new variational model which combines linear curvature, known for its ability to generate smooth transformations,^9,20 and MI as a distance measure^10,21 for the spatial transform.

Figure 1.

Simultaneous segmentation and linear curvature registration based on MI model noisy synthetic image and CT image. First column segmented template image, T & ${\varphi }_0(\mathbf{x})$ . Second column reference image R. Third column the deformation field. Forth column the transformed template image R & ${\varphi }_0(\mathbf{x}+u(\mathbf{x})$ . Last column the template image after transformation, $T(\mathbf{x}+u(\mathbf{u}))$ .

The outline of this paper is as follows. In the next section, we review the existing models for monomodal image regmentation and multimodal image registration. Then, we introduce our proposed new model for multimodal image regmentation. In the last section, we show some numerical tests including comparison.

Related work

Deformable multimodal registration attempts to find spatial correspondence between image pairs that are acquired with distinct scanning protocols. In the registration process, these correspondences are estimated by determining the spatial transformation which is applied to a template image to make similar to the reference image. Generally, image registration is driven by the chosen similarity measure and the chosen deformation model. For images taken from the same modality, the image similarity usually is considered to be high, thus simple and direct image similarity measures are used, for instance, the sum of squared intensity differences (SSD). We start this section with an overview of the work introduced by Ibrahim-Rada-Chen ⁷ for monomodal regmentation of two given images. The model jointly merges the segmentation and registration into a variational function represented in a level set. This method improves the Guyader-Vese model ²² by replacing the nonlinear elastic term with linear curvature model and adding a weighted Heaviside sum of SSD term.

Let R be a reference image and T template image in a given domain $\Omega \subset R^d$ , where d is the spatial dimensional of the given data. The registration problem is formulated as a transformation vector $\varphi (\mathit{x})$ such that $T(\varphi (\mathit{x}))=R(\mathit{x})$ . In the variational approach, image registration, the transformation is $\varphi (\mathit{x})=\mathit{x}+\mathit{u}(\mathit{x})$ , with $\mathit{u}(\mathit{x})$ the displacement vector field $\mathit{u}(\mathit{x})=(u_1(\mathit{x}),u_2(\mathit{x}))$ . The aim of the Ibrahim-Rada-Chen ⁷ model is to match the contour of the template image and at the same time segment the reference image to show the deformation of displacement field lead by the segmentation process. The segmentation of the template image T is represented by the zero level line ${\varphi }_x:\Omega \to ℝ$ which is the template contour Γ, denoted as in. ⁷ The joint functional proposed by Ibrahim-Rada-Chen ⁷ monomodal image regmentation is given as follows:

$$\begin{array}{l}\underset{c_1,c_2,\mathit{u}(\mathit{x})}{\min }\mathcal{J}=\mu {\mathcal{D}}^{\text{SSDH}}(T,R,{\varphi }_0(\mathit{x}),\mathit{u}(\mathit{x}))\\ +{\lambda }_1{\displaystyle {\int }_{\Omega }|}R(\mathit{x})-c_1{|}^2{H}_{\mathrm{\epsilon }}({\varphi }_0(\mathit{x}+\mathit{u}(\mathit{x})))\text{d}\mathit{x}\\ +{\lambda }_2{\displaystyle {\int }_{\Omega }|}R(\mathit{x})-c_2{|}^2(1-H_{\mathrm{\epsilon }}({\varphi }_0(\mathit{x}+\mathit{u}(\mathit{x}))))\text{d}\mathit{x}+\alpha {\mathcal{S}}^{\text{LC}}(\mathit{u})\end{array}$$

(1)

where λ₁, λ₂, μ, α are numerical constants. $\mathit{u}(\mathit{x})=(u_1(\mathit{x}),u_2(\mathit{x}))$ is the displacement vector field, ${\varphi }_0$ is a zero-level-set, ${\mathcal{D}}^{\text{SSDH}}$ is the sum of squared difference distance measure

$$\begin{array}{l}{\mathcal{D}}^{\text{SSDH}}={\displaystyle {\int }_{\Omega }(T(}\mathit{x}+\mathit{u}(\mathit{x}))-R(\mathit{x}){)}^2{H}_{\mathrm{ϵ}}({\varphi }_0(\mathit{x}+\mathit{u}(\mathit{x})))\text{}\text{d}\mathit{x}\hfill \end{array}$$

(2)

which is weighted by the regularised Heaviside function

$$H_{\mathrm{ϵ}}({\varphi }_0(\mathit{x}+\mathit{u}(\mathit{x})))=\frac{1}{2}(1+\frac{2}{\pi }\arctan \frac{{\varphi }_0(\mathit{x}+\mathit{u}(\mathit{x}))}{\mathrm{ϵ}}).$$

(3)

The term ${\mathcal{S}}^{\text{LC}}(\mathit{u})$ in (1) is the curvature regulariser term

$${\mathcal{S}}^{\text{LC}}(\mathit{u})={\displaystyle {\int }_{\Omega }{(\Delta u_1)}^2}+{(\Delta u_2)}^2\text{ d}\mathit{x}.$$

(4)

The values of c₁ and c₂ in equation (1) present the average intensity values inside and outside the boundary ${\varphi }_0(\mathit{x})$ in the reference image. By minimizing the functional (1) containing the SSD term and linear curvature term, we get:

$$\begin{array}{l}{c}_1=\frac{{\displaystyle {\int }_{\Omega }R}(\mathit{x})H_{\mathrm{\epsilon }}({\varphi }_0(\mathit{x}+\mathit{u}))\text{}\text{d}\mathit{x}}{{\displaystyle {\int }_{\Omega }{H}_{\mathrm{\epsilon }}}({\varphi }_0(\mathit{x}+\mathit{u}))\text{}\text{d}\mathit{x}},\\ c_2=\frac{{\displaystyle {\int }_{\Omega }R}(\mathit{x})(1-H_{\mathrm{\epsilon }}({\varphi }_0(\mathit{x}+\mathit{u})))\text{}\text{d}\mathit{x}}{{\displaystyle {\int }_{\Omega }1}-H_{\mathrm{\epsilon }}({\varphi }_0(\mathit{x}+\mathit{u}))\text{}\text{d}\mathit{x}}.\end{array}$$

(5)

In multilevel representation, the registration problem is solved using quasi-Newton method by updating the deformation field $\mathit{u}(\mathit{x})$ on each level. Even though the model shows good performance, it cannot cope with multimodal images.

Referring to registration for multimodality, different techniques have been introduced.^23–25 A significant work, proposed by Modersitzki, ⁹ based on mutual information and curvature regularizer, competes the other registration models. The joint functional for multimodal registration is given as follows:

$$J(\mathit{x}+\mathit{u}(\mathit{x}))=\mathcal{D}[T,R]+\alpha {\mathcal{S}}^{\text{LC}}$$

(6)

where ${\mathcal{S}}^{\text{LC}}$ is the curvature regularization term, α is a numerical constant and $\mathcal{D}[T,R]$ presents the distance measures considered as functionals in multimodal image pair T and R. Modersitzki ⁹ suggested to use mutual information (MI) (shortly reviewed below as CASE 1) or Normalized Gradient Field (NGF) (shortly reviewed below as CASE 2) to measure the similarity distance between two multimodal images. CASE 1: The joint functional based on MI, $\mathcal{D}[T,R]={\mathcal{D}}^{\text{MI}}$

$$J(\mathit{x}+\mathit{u}(\mathit{x}))={\mathcal{D}}^{MI}+\alpha {\mathcal{S}}^{\text{LC}}$$

(7)

where ${\mathcal{D}}^{MI}$ presents mutual information as a distance measure between the given images T and R and it is given with the formula:

$${\mathcal{D}}^{\text{MI}}={\displaystyle {\int }_{\Omega }{\rho }_{[T]}}\log {\rho }_{[T]}dt+{\displaystyle {\int }_{\Omega }{\rho }_{[R]}}\log {\rho }_{[R]}dr-{\displaystyle {\int }_{{\Omega }^2}{\rho }_{[T,R]}}\log {\rho }_{[T,R]}d(t,r)$$

(8)

where ${\rho }_{[T]}(t)$ and ${\rho }_{[R]}(r)$ are marginal densities which are expressed as:

$${\rho }_{[T]}(t)={\displaystyle {\int }_{\Omega }{\rho }_{[T,R]}}(t,r)dr,{\rho }_{[R]}(r)={\displaystyle {\int }_{\Omega }{\rho }_{[T,R]}}(t,r)dt$$

(9)

and ${\rho }_{[T,R]}(t,r)$ presents joint gray value distributions. Its discretized form is expressed as:

$${\rho }_{[T,R]}(t,r)=\frac{1}{m}{\displaystyle \sum _{j=1}^{m}K}(t-T(\mathit{x}+\mathit{u}(\mathit{x})),\sigma )K(r-R(\mathit{x}),\sigma )$$

(10)

with σ that stands for the width of Parzen density estimator and K for kernel and m presents the size of the sample image. The kernel positions in template and reference image are represented with t and r respectively in the gray value space for j = 1. CASE 2: The joint functional based on NGF, $\mathcal{D}[T,R]={\mathcal{D}}^{\text{NGF}}$

$$J(\mathit{x}+\mathit{u}(\mathit{x}))={\mathcal{D}}^{\mathit{NGF}}+\alpha {\mathcal{S}}^{\text{LC}},$$

(11)

${\mathcal{D}}^{\mathit{NGF}}$ presents Normalized Gradient Force distance measure for two given multimodal images, T and R which is defined as follows

$$\begin{array}{l}{\mathcal{D}}^{\text{NGF}}={\displaystyle {\int }_{\Omega }1}-{(n{[T(\mathit{x}+\mathit{u}(\mathit{x}))]}^{T}n[R(\mathit{x})])}^{2}d\mathit{x}\hfill \end{array}$$

(12)

where

$$\begin{array}{l}n{[T(\mathit{x}+\mathit{u}(\mathit{x}))]}^{\text{T}}n[R(\mathit{x})]\\ \approx {\left(\frac{\nabla T_i}{\sqrt{\Vert \nabla T_i{\Vert }^2+{\eta }^2}}\right)}^{\text{T}}\left(\frac{\nabla R_i}{\sqrt{\Vert \nabla R_i{\Vert }^2+{\eta }^2}}\right)\end{array}$$

(13)

where η is an important constant parameter because it determines the edge and what has to be within the specified noise level. Moreover, the terms $\nabla T_i$ and $\nabla R_i$ present gradient intensity of template image and reference image, respectively which are computed as

$$\nabla T_i={[{({\partial }_1^{h}T(\mathit{x}+\mathit{u}(\mathit{x})))}_i,{({\partial }_2^{h}T(\mathit{x}+\mathit{u}(\mathit{x})))}_i]}^T,$$

(14)

$$\nabla R_i=[{({\partial }_1^{h}R(\mathit{x}))}_i,{({\partial }_2^{h}R(\mathit{x}))}_i]$$

(15)

The NGF can be considered as the L2 norm of R, the residual of the alignment of the normalised gradients of two images at a pixel position x ,

$$r_{\mathit{x}}^h=1-{(n{[T(\mathit{x}+\mathit{u}(\mathit{x}))]}^{T}n[R(\mathit{x})])}^2$$

(16)

for discrete images T and R of size N × N using finite difference method. The images are discretised on a uniform mesh using vertex centred discretisation where ${\mathit{x}}_{i,j}$ denotes the pixel position or on a non-uniform mesh with finite difference method. The gradient is calculated using

$${\partial }_{x_1}{T}^h({\mathit{x}}_{i,j})=\frac{T^h({\mathit{x}}_{i+1,j})-T^h({\mathit{x}}_{i-1,j})}{2h},{\partial }_{x_2}{T}^h({\mathit{x}}_{i,j})=\frac{T^h({\mathit{x}}_{i,j+1})-T^h({\mathit{x}}_{i,j-1})}{2h},$$

where the first order central finite difference scheme is used to approximate the first order derivatives. The discretized form of NGF distance measure:

$${\mathcal{D}}^{\text{NGF}}=h^2{\displaystyle \sum _{i=1}^{N^2}{r}_i}$$

(17)

with the supposition that the spacial discretization is $h_1=h_2=h$ . In short words, NGF is based on the alignment of the edges in the reference and template images. The gradient is normalised with its magnitude.

The proposed new joint regmentation model for multimodality

Our joint functional model for multimodal image regmentation combines an active contour without edges, a curvature regulariser, and a mutual information distance measure. Our proposed regmentation functional is given with the formula:

$$\begin{array}{l}\underset{c_1,c_2,\mathit{u}(\mathit{x})}{\min }\mathcal{J}={\lambda }_1{\displaystyle {\int }_{\Omega }|}R(\mathit{x})-c_1{|}^2{H}_{\mathrm{ϵ}}({\varphi }_0(\mathit{x}+\mathit{u}(\mathit{x})))\text{}\text{d}\mathit{x}\hfill \\ +{\lambda }_2{\displaystyle {\int }_{\Omega }|}R(\mathit{x})-c_2{|}^2(1-H_{\mathrm{ϵ}}({\varphi }_0(\mathit{x}+\mathit{u}(\mathit{x}))))\text{}\text{d}\mathit{x}+\gamma {\mathcal{S}}^{\text{LC}}(\mathit{u})\hfill \\ +\mu {\mathcal{D}}^{\text{MIH}}(T,R,{\varphi }_0(\mathit{x}),\mathit{u}(\mathit{x}))\hfill \end{array}$$

(18)

where ${\mathcal{S}}^{\text{LC}}$ is the linear curvature, ${\mathcal{D}}^{\text{MIH}}$ represents the mutual information distance term which is weighted by the parameter μ. The ${\mathcal{D}}^{\text{MIH}}$ term is evaluated as the mutual information between foreground of the template and the reference, $T{H}_{\mathrm{ϵ}}$ and $R{H}_{\mathrm{ϵ}}$ , respectively:

$${\mathcal{D}}^{\text{MIH}}={\displaystyle {\int }_{\Omega }{\rho }_{[T{H}_{\mathrm{ϵ}}]}}\log {\rho }_{[T{H}_{\mathrm{ϵ}}]}dt+{\displaystyle {\int }_{\Omega }{\rho }_{[R{H}_{\mathrm{ϵ}}]}}\log {\rho }_{[R{H}_{\mathrm{ϵ}}]}dr-{\displaystyle {\int }_{{\Omega }^2}{\rho }_{[T{H}_{\mathrm{ϵ}},R{H}_{\mathrm{ϵ}}]}}\log {\rho }_{[T{H}_{\mathrm{ϵ}},R{H}_{\mathrm{ϵ}}]}d(t,r)$$

The joint functional is solved with respect to the displacement field using discretise then optimise approach based on the quasi-Newton method in multilevel framework for faster implementation. The grid points are located at the center of the cell ${\Omega }^h=\left\{{\mathit{x}}_{i,j}=(x_{1,i},x_{2,j})=((i-0.5)h,(j-0.5)h)\right\},1\le i,j\le N$ . We re-define the solution vector $\mathit{U}={\left[\begin{array}{c}{\mathit{u}}_1\\ {\mathit{u}}_2\end{array}\right]}_{2N^2\times 1},\mathit{x}={\left[\begin{array}{c}{\mathit{x}}_1\\ {\mathit{x}}_2\end{array}\right]}_{2N^2\times 1},$ where ${\mathit{u}}_1,{\mathit{u}}_2$ present displacement vectors and ${\mathit{x}}_1,{\mathit{x}}_2$ position vectors. Thus, the discretised form of the joint functional in (18), by a finite difference method is:

$$\begin{array}{l}\underset{c_1,c_2,\mathit{U}}{\min }\mathcal{J}={\lambda }_1{\displaystyle \sum _{i,j=1}^N|}R\left({\mathit{x}}_{i,j}\right)-c_1{|}^2{H}_{\mathrm{\epsilon }}\left({\varphi }_0\right({\mathit{x}}_{i,j}+\mathit{u}\left({\mathit{x}}_{i,j}\right)\left)\right)\\ +{\lambda }_2{\displaystyle \sum _{i,j=1}^N|}R\left({\mathit{x}}_{i,j}\right)-c_2{|}^2(1-H_{\mathrm{\epsilon }}({\varphi }_0({\mathit{x}}_{i,j}+\mathit{u}({\mathit{x}}_{i,j}\left)\right)\left)\right)\\ \\ +\mu {\eta }_t{\eta }_r{\displaystyle \sum _{i=1}^{n_t}{\displaystyle \sum _{j=1}^{n_r}{\rho }_{i,j}}}\log ({\rho }_{i,j}+\mathrm{\epsilon })\\ +\gamma {\displaystyle \sum _{l=1}^2{\displaystyle \sum _{i,j=1}^N(}}-4u_l\left({\mathit{x}}_{i,j}\right)+u_l\left({\mathit{x}}_{i+1,j}\right)\\ +u_l\left({\mathit{x}}_{i-1,j}\right)+u_l\left({\mathit{x}}_{i,j+1}\right)+u_l\left({\mathit{x}}_{i,j-1}\right){)}^2\end{array}$$

(19)

where ${\eta }_t=(t_n-t_0)/n_t$ and ${\eta }_r=(r_n-r_0)/n_r$ , n_t and n_r present the number of grid. The ϵ is just a small numerical constant to prevent extra considerations such as $″0log{0}^{″}$ . Furthermore, we are using homogeneous Neumann boundary conditions approximated by one side differences $u_l({\mathit{x}}_{i,1})=u_l({\mathit{x}}_{i,2}),u_l({\mathit{x}}_{1,j})=u_l({\mathit{x}}_{2,j}),u_l({\mathit{x}}_{i,N-1})=u_l({\mathit{x}}_{i,N}),u_l({\mathit{x}}_{N-1,j})=u_l({\mathit{x}}_{N,j}),l=1,2.$

Minimizing equation (19) brings a system of nonlinear equation with unknown U:

$$△\mathcal{J}=0$$

(20)

where

$$\begin{array}{l}\Delta \mathcal{J}={\lambda }_1{\displaystyle \sum _{i,j=1}^N|}R({\mathit{x}}_{i,j})-c_1{|}^2{H}_{\mathrm{ϵ}}({\varphi }_0({\mathit{x}}_{i,j}+\mathit{u}({\mathit{x}}_{i,j}))\hfill \\ \hspace{0.17em}+{\lambda }_2{\displaystyle \sum _{i,j=1}^N|}R({\mathit{x}}_{i,j})-c_2{|}^2(1-H_{\mathrm{ϵ}}({\varphi }_0({\mathit{x}}_{i,j}+\mathit{u}({\mathit{x}}_{i,j})))\hfill \\ \hspace{0.17em}+\mu {\mathcal{D}}^{\text{MIH}}+\gamma {\mathbf{u}}^T\mathbf{B}\mathbf{u}\hfill \end{array}$$

(21)

The gradient of the regularization term is computed as the multiplication of B matrix (constant matrix) of size $2N^2\times 2N^2$ that contains the coefficients of U. The matrix B is written as

$$\mathbf{B}=\left[\begin{array}{cc}{L}^{T}L& 0\\ 0& L^{T}L\end{array}\right],$$

(22)

where L is a block tridiagonal matrix from the regularization term. For finding the solution of our minimization problem, we start with zero initial guess, $\mathit{U}=0,$ solving

$$\mathit{H}\delta \mathit{U}=-\mathit{G}$$

(23)

where

$$\mathit{G}={\nabla }_u\mathcal{J}$$

(24)

for $\delta \mathit{U}$ and update $\mathit{U}\leftarrow \mathit{U}+\tau \delta \mathit{U}$ with τ as the Armijo line search parameter. ²⁶ H and G are the Hessian and gradient matrix for the functional $\mathcal{J}$ in equation (19) with respect to the displacement vector U . The algorithm for the proposed regmentation model based on MI is given in Algorithm 1 where the multilevel approach is adapted for fast computation of the model.

Algorithm 1 Joint Regmentation Method Algorithm for Multimodal Images: $({\mathcal{D}}^{\text{MIH}},{\mathcal{S}}^{\text{LC}},\varphi (\mathbf{x}),u(\mathbf{x}))\leftarrow \mathit{JRMI}(R,T,{\lambda }_1,{\lambda }_2,\mu ,\mathit{U},{\varphi }^0$ )

Step 1. Initialize the level set function $\varphi (\mathbf{x})$ to be a binary function as follows:

$${\varphi }^0=\varphi (\mathbf{x},t=0)=\{\begin{array}{l}-\rho \text{if }\mathbf{x}\in {\Omega }_0-\partial {\Omega }_0\hfill \\ 0\text{if }\mathbf{x}\in \partial {\Omega }_0\hfill \\ \rho \text{if }\mathbf{x}\in \Omega -{\Omega }_0\hfill \end{array}$$

where $\rho >0$ is a constant, ${\Omega }_0$ is a subset in the image domain Ω and $\partial {\Omega }_0$ is the boundary of ${\Omega }_0$ , and zero inital guess for the displacement field

$$\mathit{U}=0.$$

Step 2. Finding deformation field using multilevel strategy:

for $\mathit{level}=\mathit{coarseLevel}:\mathit{finestLevel}$ do

* Compute the deformation field on each level

  ${\mathit{U}}^{(\mathit{level})}\leftarrow \mathit{Registration}(R^{(\mathit{level})},T^{(\mathit{level})},{\mathit{U}}^{(\mathit{level})},{\varphi }^0$ ),

using Quasi-Newton method as follows:

$•$ • update c₁ and c₂ using equation (5)

$•$ • update ${\mathit{U}}^{(\mathit{level})}$ and solve equation (23)

$•$ • continue until convergence criterion is satsified

* ${\mathit{U}}^{(\mathit{level})}$ is interpolated to each next finer level after registration

end for

Numerical results and conclusions

In this section, we present several examples for synthetic as well as real data in comparison with the linear curvature model based on MI applied on a synthetic image and a set of real images. In addition, the results of our proposed joint regmentation are compared to image registration applying mutual information as a distance measure and image segmentation simultaneously. The image pairs used in all our experiments have significantly different intensity profiles, deformation and presence of noise. In each iteration we compute the Jacobian matrix of the transformation,

$$J=\left[\begin{array}{cc}1+\frac{\partial u_1}{\partial x_1}& \frac{\partial u_1}{\partial x_2}\\ \frac{\partial u_2}{\partial x_1}& 1+\frac{\partial u_2}{\partial x_2}\end{array}\right]$$

(25)

and the minimal value of Jacobian matrix $\mathcal{F}=\min (\det (J))$ is also calculated to make sure there is no folding or cracking in the deformation field if its value is greater than 0. Figures 2 and 3 present the results of the proposed model in comparison with Ibrahim et al. ⁷ regmentation model, and linear curvature image registration model based on MI and NGF proposed by Modersitzki et al.²⁸ The results show that the Ibrahim et al. ⁷ model, which uses the SSD similarity measure, is not capable to register and segment the given synthetic multimodal images whereas the proposed model can successfully cope with it. The fail of the method proposed by Ibrahim et. al ⁷ is expected as the model is not designed for multimodal images. On the other hand, the linear curvature registration model which applies MI, proposed by Modersitzki et al.,²⁸ provides larger deformation than the same model using NGF, referring to Figure 3. Even though, this model still has a poor registration in comparison with the proposed method, shown in Figure 2 last three columns.

Figure 2.

Comparison between monomodal image regmentation ⁷ and the proposed model. The monomodal joint regmentation fails to register with parameters ${\lambda }_1={\lambda }_2=\mu =1$ and α = 5 (in (c) and (d)). The proposed model successfully deforms the given template image with ${\lambda }_1={\lambda }_2=\gamma =0.0005$ and $\mu =2.0e-05$ (in (e), (f) and (g)). The MI distance difference of our proposed model before is $D^{\mathit{MIH}}=-0.22806$ which successfully minimizes the distance after to $D^{\mathit{MIH}}=-0.53424$ . (a) T and ${\varphi }_0(\mathbf{x})$ , (b) R, (c) SSD $\mathcal{F}=-0.0994$ , (d) SSD $T(\mathbf{x}+u(\mathbf{x}))$ , (e) MIH $\mathcal{F}=0.27148$ , (f) MIH R and ${\varphi }_{(\mathbf{x}+u(}\mathbf{x}))$ , (g) MIH $T(\mathbf{x}+u(\mathbf{x}))$ .

Figure 3.

Comparison between the linear curvature models for synthetic images using MI and NGF²⁸ (refer to (a),(b) and (c),(d) for MI and NGF results of derformation field and registered template image, respectively). The proposed model is compared to both of these models and shows that successfully deformes the given template image. In contrast, curvature model using MI gets stuck in local minima which is demonstrated by the distance measurements before registration $D^{MI}=-1.0303$ and after the registration $D^{MI}=-1.0426$ . Additionally, the curvature model using NGF does not show large deformation with regularizer parameter α = 2. MI $\mathcal{F}=-0.5835$ (a) MI $T(\mathbf{x}+u(\mathbf{x}))$ (b) NGF $\mathcal{F}=-0.7465$ (c) NGF $T(\mathbf{x}+u(\mathbf{x}))$ (d) MIH $\mathcal{F}=0.27148$ (e) MIH R and ${\varphi }_{(\mathbf{x}+u(}\mathbf{x}))$ (f) MIH $T(\mathbf{x}+u(\mathbf{x}))$ (g)

Figures 4 and 5 show the results on real multimodal data images. In the first row of both figures, we show comparison results between the proposed model and linear curvature model²⁸ of two chest images (T1 and T2 weighted images). The second and the third row have the thorax (PET and CT images²⁸) and brain images with high deformation and presence of noise of a T2 weighted CT image for the template and a MRI image for the reference. We clearly see that the proposed model delivers good results, referring to the fifth, sixth, and seventh column of Figure 4, whereas the linear curvature model²⁸ shown in the third and fourth column stuck in local minima due to highly non convexity of MI functional. Our multimodal joint regmentation provides accurate segmentation and registration in comparison also to linear curvature based on NGF model.²⁸ These results are shown in Figure 5. The gradient field distance measure involves second order image derivatives, as a consequence there can be a problem dealing with noisy images. Table 1 shows the parameters and distances measurements before and after the performance of linear curvature model²⁸ and our model. We note that the parameter α, λ₁, λ₂, μ, and γ are tuned according to the paired modalities. It can be noticed a similarity for those parameters in between the chest image case where T1 and T2 weighted CT images are involved and the image of the brain with T2 weighted CT image for the template and an MRI image for the reference, whereas those parameters change for the thorax image where PET image for the template is paired with a CT image for the reference. This drawback is expected similar to all the other inverse problems. In conclusion, the proposed model is suitable for jointly segmenting and registering images with intensity difference, severe deformation, and presence of noise. The model avoids in this way the need of pre-or post-processes such as pre-registration step which might be required for the segmentation task or vice-versa.

Figure 4.

Comparison between linear curvature model using MI²⁸ and the proposed model for real images data. T & ${\varphi }_0(\mathbf{x})$ (a) R (b) $\mathcal{F}=-0.3094$ (c) $T(\mathbf{x}+u(\mathbf{x}))$ (d) $\mathcal{F}=0.53174$ (e) R & ${\varphi }_0(\mathbf{x}+u(\mathbf{x}))$ (f) $T(\mathbf{x}+u(\mathbf{x}))$ (g) T & ${\varphi }_0(\mathbf{x})$ (h) R (i) $\mathcal{F}=-0.4770$ (j) $T(\mathbf{x}+u(\mathbf{x}))$ (k) $\mathcal{F}=0.56764$ (l) R & ${\varphi }_0(\mathbf{x}+u(\mathbf{x}))$ (m) $T(\mathbf{x}+u(\mathbf{x}))$ (n) T & ${\varphi }_0(\mathbf{x})$ (o) R (p) $\mathcal{F}=-0.3856$ (q) $T(\mathbf{x}+u(\mathbf{x}))$ (r) $\mathcal{F}=0.46249$ (s) R & ${\varphi }_0(\mathbf{x}+u(\mathbf{x}))$ (t) $T(\mathbf{x}+u(\mathbf{x}))$ (u)

Figure 5.

Comparison between linear curvature model using NGF²⁸ and the proposed model for real images data. T & 0 f (x) (a) R (b) F = –0.3094 (c) T(x + u(x)) (d) F = 0.56233 (e) R & 0 f (x+u(x)) (f) T(x + u(x)) (g) T & 0 f (x) (h)R (i) $\mathcal{F}=-0.3926$ (j) $T(\mathbf{x}+u(\mathbf{x}))$ (k) $\mathcal{F}=0.56764$ (l) R & ${\varphi }_0(\mathbf{x}+u(\mathbf{x}))$ (m) $T(\mathbf{x}+u(\mathbf{x}))$ (n) T & ${\varphi }_0(\mathbf{x})$ (o) R (p) $\mathcal{F}=-0.5284$ (q) $T(\mathbf{x}+u(\mathbf{x}))$ (r) $\mathcal{F}=0.46249$ (s) R & ${\varphi }_0(\mathbf{x}+u(\mathbf{x}))$ (t) $T(\mathbf{x}+u(\mathbf{x}))$ (u)

Table 1.

The first, second and third columns show the linear curvature model ¹⁷ parameters and distance measurements before and after, while the other columns are referring to our model, shown in Figure 4.

Figure 4	α	before ${\mathcal{D}}^{\text{MI}}$	after ${\mathcal{D}}^{\text{MI}}$	λ 1	λ 2	μ	γ	before ${\mathcal{D}}^{\text{MIH}}$	after ${\mathcal{D}}^{\text{MIH}}$
Chest	0.05	–1.4938	–1.3204	0.01	0.01	0.01	0.00045	–1.1677	–1.6247
Thorax	0.5	–0.2478	–0.2608	2	2	0.1	2	–0.40906	–0.88796
Brain	0.003	–1.1994	–1.2234	0.002	0.002	0.002	0.000035	–1.4477	–1.5911