An Improved SIMPLEC Scheme for Fluid Registration

. The image registration is always a strongly ill-posed problem, a stable numerical approach is then desired to better approximate the deformation vectors. This paper introduces an efficient numerical implementation of the Navier Stokes equation in the fluid image registration context. Although fluid registration approaches have succeeded in handling large image deformations, the numerical results are sometimes inconsistent and unexpected. This is related, in fact, to the used numerical scheme which does not take into consideration the different properties of the continuous operators. To take into account these properties, we use a robust numerical scheme based on finite volume with pressure correction. This scheme, which is called by the Semi-Implicit Method for Pressure-Linked Equation-Consistent (SIM-PLEC), is known for its stability and consistency in fluid dynamics context. The experimental results demonstrate that the proposed method is more efficient and stable, visually and quantitatively, compared to some classical registration methods


Introduction
Image registration is a widely used technique in image processing; its main principle is to compute geometrical correspondences between two or more given images [7]. It is considered as a useful tool in numerous applications: including astronomy, robotics and especially in bio-medical imaging [19,25].
Generally, image registration problems can be classified into two types: intensity-based methods and feature-based methods. The first type methods are based on the image's gray spaces [2,3,20], while the second ones exploit the difference features existing between the given images, such as regions and edges, to find the correspondences in the image's feature spaces [23]. In most cases, the image registration task is formulated as an optimization problem involving a distance measure to evaluate the similarity. However, due to different illuminations (grey-levels) of the image regions, the pairing cannot be achieved successfully, which makes the image registration an ill-posed problem. An additive regularization term is then added, motivated by the nature of the transformation. In fact, each regularizer produces a different registration model, and this choice is very vital for the solution's existence and uniqueness. One of the classical regularizers is the elastic regularization [4]. Even if this model has been addressed in several image registration problems [24], it is still limited, since it does not handle largely deformed data sets [25]. Another weakness of the linear elasticity regularization is that it does not take in consideration discontinuous displacement. A more convincing choice for preserving discontinuities of the displacement was the total variation-based regularization (TV), see [28]. The TV model regularizes the velocity/displacement/mapping, in contrast, it gives unsatisfactory registration results when the displacements are very smooth. A more robust regularization is proposed in [10], which takes into account the advantages of the stopped diffusion of the total variation registration models. Recently, an improved discontinuity-preserving model was proposed in [30] motivated by several regularization techniques used in vectorvalued image denoising, which gives promising results; but didn't tackle large deformations. In the context of non-parameteric registration, Christensen [6] developed one of the most effective non parametric registration called the fluid registration [8,9]. The purpose of the fluid image registration is to compute the deformation u(x, t) at time t for a given force field f , while the deformed image is considered in a viscous fluid. The deformations are then governed by the Navier Stokes equations for momentum conservation (where the pressure is neglected). This equation is given by (1.1) The component ∆v is called the viscous term, while the component term ∇(∇ · v(x, t)) allows for contraction and expansion of the fluid. The second part of (1.1) defining material derivative of the displacement u. The constants µ and ν are the Lamé coefficients. Although the fluid image registration appears easy to solve, in reality it is quite complicated. One of the key complications is the choice of demons force [15]. Moreover, fluid image registration is not based on a direct optimization approach. Also, since we deal with a set of nonlinear partial differential equations (PDE), great care must be taken in the choice of the discretization method. In addition, when displacement u is discontinuous, the fluid image registration gives poor results since it skips the discontinuity regions. Another well-known drawback of this method is the high computational cost. Indeed, once the force field f is fixed, the first part of Equation (1.1) is solved for v(x, t) using the successive over-relaxation (SOR) scheme [6], which computes an accurate fluid model at the expense of a large computational time. Then, an explicit Euler scheme is used to advance u in time t. The method requires at each iteration the computation of the Jacobian matrix of the displacement field, which is also computationally very expensive. This framework is thus time-consuming, which motivates the search of faster implementations. Recently, a new fluid approach was proposed using the Semi-Implicit Method for Pressure-Linked Equation (SIMPLE) scheme, where the pressure is well defined in the context of image registration [1]. However, if the registered images are not with the same modalities and in the presence of noise, this method becomes inconsistent with the apparition of some artefacts. Hence the introduction of another more accurate scheme.
The contributions in this paper are summarized below. Firstly, we define the pressure term, in the image registration context, as the effect of each region with respect to the nearest one, which has not been treated by Christensen in the previous model (1.1).
Secondly, we propose a stable, consistent and fast numerical implementation of the Navier Stokes equations for an incompressible fluid to resolve the fluid image registration problem. Which allows to introduce ideas from computational fluid dynamics in the image analysis context [29].
The use of the finite difference approximation in the image registration is very limited since it does not take into account the continuous and discontinuous operator properties. An alternative and robust discretization scheme is then needed. In this paper, we propose a finite volume scheme type to handle different properties of the fluid registration in the discretization process of the proposed Navier Stokes equation [13]. The use of this scheme can efficiently take into account the defined pressure term in the registered image. After the success of the Semi-Implicit Method for Pressure Linked Equations [26] algorithm in solving the Navier Stokes equations and many other problems in computational fluid dynamics (CFD) [14,18,26], the Semi-Implicit Method for Pressure Linked Equations-Consistent (SIMPLEC) [12] is also proposed with good consistency. However, no such attempt has been proposed to solve the fluid image registration problem. Then, in this work we introduce the SIM-PLEC algorithm as a discretization approach to solve the Navier-Stokes equation applied to image registration problems. We can see that the numerical results confirm that our proposed scheme is more efficient compared with other classical methods.
This paper is organized as follows: Section 2 explains the main concepts of the image registration and introduces the proposed Navier Stokes equation for solving the fluid registration problem. In Section 3, we propose the SIM-PLEC discretization scheme. Finally, numerical results and comparisons of our approach with classical registration methods are presented in Section 4.

Mathematical model of fluid registration
The registration problem is usually based on a minimization problem between two images, template image and reference image [21]. Given two images T , R : Ω ⊂ R n → R, compactly supported on Ω (a bounded convex domain), where n denotes spatial dimension of the given images (in our case n = 2). The purpose of registration is to look for a transformation u : Ω ⊂ R n such that ideally T (u(x)) looks like R(x) as much as possible for all x ∈ Ω. In summary, the desired transformation u is achieved by minimizing a so-called distance measure D. Since this problem is ill-posed, an appropriate regularization S is used. The variational formulation of the image registration problem consists of finding a minimizer u of where A denotes the set of admissible transformations. The regularization term S is in general based on the gradient of u, noted ∇u. There is in fact many choices of the regularization term, the widely used are: diffusion registration [22] and elastic registration [21]. One of the more successful regularization choices was the fluid one [11], which is formulated such as an elastic potential of ∂ t u. To solve this problem, the Euler-Lagrange identity is used, which coincides with the resolution of the Navier-Lamé equation. There are also other regularizations derived from the diffusion one, such as the curvature term and the hyperelastic energy [5], which do not have good physical motivations.
Based on the proprieties of this fluid regularization, we propose a dynamic equation for the incompressible Newtonian fluids (2.2), which are governed by the Navier Stokes equations. This equation coupled the velocity vector field v(x, t) to a scalar pressure p(x, t) such as 3) The main analogy followed in this paper is the parallel between the incompressible Newtonian fluid and the image velocity of each pixel v(x, t) under the image registration concept. The introduced Navier Stokes equation (2.2) is well posed in the image registration task. In fact, the parameter ν is supposed to be the factor of the diffusion in the imaging problems. The pressure p is also modelled in the imaging task, which represents the effect of each region on the nearest one in the image during the registration process. In fact, since the pressure is given in gradient form, it may represent the contours. We suppose in the following that p is an external force. The term ∇ · v = 0 is well posed since the pixels are incondensable. Finally, f is supposed to be the external forces obtained by the gradient of the distance D between the two images; T and R. A classical choice for this distance is the sum square difference (SSD) measure defined as Hence, the external forces are computed such as: For this choice, we can effectively assure the existence of a solution to the problem (2.1) using techniques in [11]. There are many advantages in the use of the Navier Stokes equation. Firstly, we don't have any problem with the existence of a solution, since it is well-developed in the literature [29]. Secondly, there are many and stable numerical approaches that we can use to resolve this equation derived from a classical example of fluid dynamics. For a mathematical transparency of the Navier Stokes equation (2.2), a convenient way is to consider it nondimensional form [16]. This is obtained by introducing a reference length L * and a reference time T * , then we set where V * is the reference velocity defined as V * = L * T * . By a substitution into (2.2), we obtain for v ′ , P ′ , f ′ the following equation where R e is a nondimensional constant called the Reynolds number defined as 1 Re = ν L * V * , and to avoid the complexity of the notations, we keep the same notation in (2.2), i.e., we substitute v ′ , P ′ , f ′ by v, p, f . Then, we rewrite the Equation (2.5) as follows: In the following section, we discuss the proposed discretization scheme for the Navier Stokes equation (2.6), in the image registration context.

Discretization
To solve the Navier Stokes equation (2.6), we use a finite volume discretization based on a staggered control volume illustrated in the Figure 1. To calculate the variables v 1 (the horizontal component of velocity), v 2 (the vertical component of velocity) and p (Pressure), we use three staggered mesh [27]. The control volumes for v 1 and v 2 are displaced with respect to the control volume for the continuity equation. In the Figure 1, "P " denotes the node at which the partial differential equation is approximated, and "E", "W", "N", "S" are its neighbours. Cell faces "e", "w" for v 1 and "n", "s" for v 2 are in the midway between the nodes.
Firstly, we rewrite the momentum equation of (2.6) in a differential form for both velocity components To solve this system, we discretize each equation separately and we use the SIMPLE method. The discretization momentum equation for v 1 is derived by integrating the first equation of (3.1) over the control volume corresponding to v 1 , using the points (E, P, n, s) shown in Figure 1 and over the time interval; from t to t + ∆t. Thus, using the fact that the velocity v 1 and v 2 do not depend on the vertical and horizontal components, respectively, under each volume control face, we have then with vol = ∆x∆y and v 0 1 is the previous iteration of v 1 (in general it represents v 1 at t = 0). To simplify the study of this problem, we introduce some new entities defined as If we substitute this new entities in the Equation (3.2), we have then which it is equivalent to Using the Equations (3.3) and (3.4), we find On the other hand, we suppose that the flow is unidirectional (since the registration is concentrated in a particular direction), laminar (since the pixels are not mixed) with a constant pressure on each grid (the effect of the neighbors pixels is neglected in a fixed grid). Thus, the first equation of (3.1) becomes Also, the second equation representing the vertical direction becomes To compute the entities (J E , J W , J N and J S ), we use the solution of the two Equations (3.6) and (3.7). Firstly, we have to rewrite the Equation (3.6) under each grid as follow (3.8) Also, the expression of the Equation (3.7) in each grid is given by The solution of the Equation (3.8) is calculated using this expression In addition, the solution of (3.9) is given by We assume that P e = (R e v 1 ) E ∆x, where P e is called the Peclet number, representing the local report on the boundary portion of the volume control for inertial and viscous forces. The Equation (3.10) is now described as We can now calculate the expression of J E using the Equation (3.10) injecting the expression of J E in (3.5), we have By the same way, we can find the others expressions (J i −v 1P ·F i ), i ∈ (W, N, S) given as follow Therefore, (3.5) amounts to

Finally, the discretization equation can be written as
; ; P n =(R e v 2 ) n ∆y; A S = F S exp(P s ) exp(P s − 1) ; P s = (R e v 2 ) s ∆y, The momentum equations for the other direction is obtained with the same way, see [26] for more details. So the discretization for v 2 is given as where b 2 = volf 2 + vol ∆t v 0 2P . Since we have calculated the different momentum equations, we need to correct the pressure in each step to verify that ∇·v(x, t) = 0.

The pressure correction equation
The main idea of this step is to improve the guessed pressure p * such that the velocity field will progressively converge to the solution of the continuity equation For any guessed pressure p * , the velocities v * 1 and v * 2 , obtained by solving the momentum equations (3.11)-(3.12), satisfy (3.14) The velocities v 1 and v 2 are obtained through (3.11)-(3.12) using the correct pressure p, which satisfy the continuity condition. The correction of the guessed pressure by p ′ = p − p * is therefore necessary to correct the velocities v * 1 and v * 2 by (3.13). The relation between p ′ and v ′ is then obtained by subtracting the Equation (3.14) from (3.11)-(3.12) (3. 15) In the SIMPLEC algorithm, the terms nb A nb v ′ 1P and nb A nb v ′ 2P are subtracted from both sides of Equation (3.15). This yields To introduce a consistent approximation, the terms nb Using the fact that v 1P and v 2P satisfy the Equation (3.17), if we replace these expressions in Equation (3.4), we obtain the following correction pressure equation: where A E = ∆yR e d 1 , A W = ∆yR e d 1 , A N = ∆xR e d 2 , A S = ∆xR e d 2 , To solve these equations we use the algebrical equation form to obtain the fields verifying the conservation equations [26]. Finally, we summarize the proposed method in the Algorithm 1.

Algorithm 1
The proposed finite volume algorithm for image registration Input: v * 1 , v * 2 and the pressure field p * , the Reynolds number. Output: The velocities v 1 , v 2 and the pressure p. Solve the momentum equations (3.14) to obtain v * 1 and v * 2 .

4:
Calculate p by adding p ′ to p * .

5:
Calculate v 1 and v 2 from their previous values using the velocitycorrection formula (3.17).

6:
Treat the corrected pressure p as a new guessed pressure p * , return to step 2.

7:
Repeat the whole procedure until convergence is reached.

Approximation of displacement u
The next step is to compute the displacement u = (u 1 , u 2 ) from the associated velocities v 1 and v 2 , calculated through Algorithm 1, based on (2.3) and an Euler scheme [1]. For each grid point x i ∈ R 2 with a fixed index i = (i 1 , i 2 ) ∈ N 2 we have u k (x i ) ∈ R 2 , and we set where u k,1 i and u k,2 i are, respectively, the first and second component approximations of u. We also set which is the velocity approximation, where v k,1 i is the first component and v k,2 i is the second one.
As the first step, to approximate the ∇u term, a centred finite difference approximation is used to compute the Jacobian matrix J k i of u k i . Secondly, for the partial time derivative ∂ t u, we use a forward finite difference approximation. Therefore, the displacement and the velocity are connected through the following Euler scheme which is performed for all i The proposed algorithm to compute the transformations u is finally summarized in Algorithm 2. Compute of the Jacobian Matrix J k i .

4:
If the relative change of the distance measure is brought below a user supplied tolerance tol = 10 −5 , the iteration is stopped,

5:
If else, calculate the new v k+1 1 and v k+1 2 from the Algorithm 1, then return to step 2.

Results and discussion
In this section, we test the performance of the proposed image registration model. In fact, we have tested our algorithm on a large image registration benchmark. We present only four tests chosen with different size and deformation. To measure the robustness of the proposed algorithm, we compare it with some competitive approaches, including the fluid registration model R F luid [11], the modified total variation (TV) regularization R M T V [17], the SIMPLE model [1] and also the improved discontinuity-preserving model R IDP proposed in [30]. In the Figure 2, we present the reference and template image for the four used tests.
Our aim is to find the deformation between the template and the reference images, and recovering the reference image by registering the template one. In Figures 3-6, we represent the registered template image using the proposed approach compared with the other methods using the difference between the reference and template images and also the deformed grid for each test.
If we focus on what happens in the error between the template image and the reference one, in all examples, we can see that the proposed method is more efficient than the other image registration methods. To evaluate the performance of the proposed approach with respect to noise reduction, we use two measures such as peak signal to noise ratio (PSNR) and the structural similarity (SSIM). The PSNR is a popular metric used to measure the quality of the estimated image, while the SSIM is a complementary measure, which gives an indication of image quality based on known characteristics of the human visual system [31]. The PSNR measures signal strength relative to noise in the image and is defined by where the M SE is the mean squared error defined by The proposed approach The SMPLE model [1] R F luid [11] R IDP [30] R M T V [17] Figure 3. In each row we present the obtained result by each method for Example 1 ; left: the registered template image; the middle: difference between reference and obtained image and in the right column we plot the deformation field.

The proposed approach
The SIMPLE model [1] R F luid [11] R IDP [30] R M T V [17] Figure 4. In each row we present the obtained result by each method for Example 2 ; left: the registered template image; the middle: difference between reference and obtained image and in the right column we plot the deformation field.

The proposed approach
The SIMPLE model [1] R F luid [11] R IDP [30] R M T V [17] Figure 5. In each row we present the obtained result by each method for Example 3 ; left: the registered template image; the middle: difference between reference and obtained image and in the right column we plot the deformation field.

The proposed approach
The SIMPLE model [1] R F luid [11] R IDP [30] R M T V [17] Figure 6. In each row we present the obtained result by each method for Example 4 ; left: the registered template image; the middle: difference between reference and obtained image and in the right column we plot the deformation field.
The SSIM is calculated on multiple windows of given image, i.e. the measurement between two windows x and y of size N × N is defined by where the variables, respectively, defined for x and y as follows: µ x and µ y , mean; σ 2 x and σ 2 y , variance; cov xy , covariance; c 1 = (k 1 L) 2 , c 2 = (k 2 L) 2 are two stabilizing constants; and L the dynamics of the pixel values, 255 for 8-bit encoded image. This metric gives an indication on the quality of the image based on the known characteristics of human visual system. Moreover, to measure the quality of the registered images, the relative reduction of the dissimilarity rel · SSD is used where u is the current optimal value and D stop is the value of D(u) at u = 0. In Table 1, the SSIM, PSNR and rel · SSD values are calculated for all the used images. The best results are represented by a bold number. Always, the proposed method outperforms the other in terms of both PSNR and SSIM. The full code for the proposed model is implemented in the MATLAB 2013. Typically, the execution of the main implemented programme requires an average of 1 ∼ 5 minutes on a 3 GHz Pentium Quad core computer for 256 × 256 grey-scale images; the runtime increases for larger image sizes. While the execution of the fluid image registration takes about 2 ∼ 7 minutes in the same conditions. However, for the two other regularizations : R IDP and R M T V , the execution time is sometimes less compared to that recorded by our method.
Finally, we want to check the speed of the convergence during time iterations of the proposed approach compared to the SIMPLE model [1]. To do that, we run the two codes for 500 time iterations for the Example 3 registration process and we compute the relative error in each iteration. The obtained results are presented in Figure 7 which shows that the approach based on SIMPLEC algorithm converge faster than the one based on SIMPLE [1].

Conclusions
A consistent numerical scheme for the fluid image registration based on Navier Stokes equation was introduced in the image registration context. A finite volume-based scheme using the SIMPLEC algorithm was performed to avoid the errors arising from the discretization part. The performance of this approach has been tested using different examples. The proposed approach outperforms some competitive ones both visually and using different criteria.