FLUID IMAGE REGISTRATION USING A FINITE VOLUME SCHEME OF THE INCOMPRESSIBLE NAVIER STOKES EQUATION

. This paper proposes a stable numerical implementation of the Navier-Stokes equations for ﬂuid image registration, based on a ﬁnite volume scheme. Although ﬂuid registration methods have succeeded in handling large deformations in various applications, they still suﬀer from perturbed solutions due to the choice of the numerical implementation. Thus, a robust numeri- cal scheme in the optimization step is required to enhance the quality of the registration. A key challenge is the use of a ﬁnite volume-based scheme, since we have to deal with a hyperbolic equation type. We propose the classical Patankar scheme based on pressure correction, which is called Semi-Implicit Method for Pressure-Linked Equation (SIMPLE). The performance of the proposed algorithm was tested on magnetic resonance images of the human brain and hands, and compared with the classical implementation of the ﬂuid image registration [13], in which the authors used a successive overrelaxation in the spatial domain with Euler integration in time to handle the nonlinear viscous. The obtained results demonstrate the eﬃciency of the proposed approach, vi- sually and quantitatively, using the diﬀerences between images criteria, PSNR and SSIM measures.


1.
Introduction. Image registration is one of today's most challenging problems in image processing, its main object being to find geometrical correspondences between two images or more, which we call the reference and the template [9,39,31,10,19,45,44]. Image registration is an important tool in various areas of applications such as astronomy, robotics and especially in bio-medical imaging [32,40,38,30,20,23,2]. These images contain a similar but not identical object. Indeed, scalar intensities are known while a transformation vector is to be calculated. To compute this vector, a natural way is to formulate the image registration 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. An additive regularization term is then needed, motivated by the nature of the transformation. Broit introduced the so-called elastic regularization [5], this term is based on linear elasticity model. Even if this model has been addressed in several image registration problems [24,11,29], it has some limitations such as it does not handle largely deformed data sets [32]. A more robust term is constructed by nonlinear elasticity which is addressed in many works, see [26,22,33,34]. The advantage of nonlinear models is that they preserve topology much better and also allow for intuitive deformations when the displacements are large. In contrast, the nonlinear models still limited in the treatment of irregular deformations. Recently, a non-local topology-preserving segmentation-guided registration method is performed to handle large deformations but it is specified to the smooth ones where the shapes to be matched are viewed as hyperelastic materials. To treat large transformation, Christensen [14] developed one of the most effective models used in the non parametric registration called the fluid registration scheme and Thirion [12,13], other derived fluid registration are also proposed in [4,8,27]. The purpose of the fluid image registration is to compute the displacement u for a given force field f , while the deforming image is considered to be embedded in a viscous fluid whose motion is governed by the Navier-Stokes equations for conservation of momentum (where the pressure is neglected). This equation is given by µ∆v(x, t) + (λ + µ)∇(∇ · v(x, t)) = f (x, t, u(x, t)), v(x, t) = ∂ t u(x, t) + v(x, t) · ∇u(x, t).
The component ∆v is called the viscous term because it constrains the velocity field spatially, while the component term ∇(∇ · v(x, t)) allows for contraction and expansion of the fluid. The second part of (1) defining material derivative of the displacement u, nonlinearity relates the velocity v and the displacement vector field. Constants µ and λ are the Lam coefficients. Although the fluid problem appears easy to solve, in reality it is quite complicated. One of the key complications is the choices of demons force [6,17,7]. Moreover, the fluid image registration is not based on an optimization approach, since we deal with a set of nonlinear partial differential equations (PDE), great care must be taken in the choice of the discretization method. Another well-known drawback of this method is the computational cost. Indeed, numerically, fixing the force field f , the first part of equation (1) is solved for v using the successive over-relaxation (SOR) scheme [14], 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. 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 for faster implementations.
The focus of this paper is on a stable and fast numerical implementation of the Navier-Stokes equations for an incompressible fluid to resolve the fluid image registration problem, introducing ideas from computational fluid dynamics into problems in image analysis [42,18,15].
The use of the finite difference approximation is very limited because it does not take into account the continuous operator properties. To overcome this problem, we choose an alternative and appropriate discretization scheme for the continuous operator that we will descritize. In this paper, we propose a finite volume type scheme to handle different properties of the fluid registration in the discretization process of the proposed Navier-Stokes equation [28], precisely, we define the pressure term in the registred image. After the success of the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) [35] algorithm (introduced by Brian Spalding and his student Suhas Patankar [37]) in solving the Navier-Stokes equations and many other problems in computational fluid dynamics (CFD) [16,35,41], no such attempt has been proposed to solve the fluid image registration problem. Because of this, in our work we use the SIMPLE algorithm as a discretization approach to solve the Navier-Stokes equation applied to image registration problem. The numerical results confirm that our proposed scheme is more efficient compared with the finite difference discretization used in previous literature [31].
This paper is organized as follows: Section 2 explains the main concepts of the image registration and the proposed Navier-Stokes equation for solving the fluid registration problem. In Section 3, we propose the discretization scheme. Finally, simulations and comparisons of our result with the finite difference method are presented in Section 4.

2.
Mathematical model of fluid registration. The registration problem is usually based on a minimization problem between two images, template image and reference image [31,25]. Given are two images T , R : Ω ⊂ R d → R, compactly supported on Ω, where d = 2. The goal is to find a transformation x+u(x) : Ω ⊂ R d or a deformation such that ideally T (x + u(x)) ≈ R(x) for all x ∈ Ω. This goal is achieved by minimizing a so-called distance measure D. As this problem is ill-posed, an appropriate regularization S is inevitable.
A variational formulation of the image registration problem is to find a minimizer u of where A denotes the set of admissible transformations and α is a regularization parameter.
The regularization term S is in general based on the gradient of u, denoted by ∇u. Important choices for this regularizer include the diffusion registration [31] (3) and the elastic registration [31,25] (4) where d is the dimension of the image, λ and µ are the Lam coefficients. One of the most successful regularization choices was the fluid one [31], given 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. Based on the properties of this regularization, we propose a dynamic equation for the incompressible Newtonian fluids (5), which are governed by the Navier-Stokes equations, this equation coupling/pairing the velocity vector field v to a scalar pressure p such that ∂v ∂t + v · ∇v = −∇p + ν∆v + f ( . , u( . )), where the velocity v is calculated according to the displacement u such that The main analogy followed in this paper is the parallel between the incompressible Newtonian fluid and the image velocity of each pixel v(x) under the image registration concept. To do so, a Navier-Stokes equation (5) is introduced, where ν is supposed to be the factor of the diffusion. Also, the pressure p is modelled as the effect of each region on the nearest one in the image during the registration process. We suppose in the following that p is an external force. On the other hand, ∇ · v = 0 is well posed since the pixels are incondensable. The analogy is summarized in Table 1. Finally, f is supposed to be the external force obtained by  Table 1. The analogy between the incompressible Newtonian fluid and the image registration.
the gradient of the distance D between the two images; T and R. A typical choice for this distance is the sum square difference (SSD) measure defined as (7) D SSD (T , R; u) = T (· + u(·)) − R(·) 2 L 2 (Ω) . Hence the external force f is computed by We can effectively assure the existence of a solution to problem (5) using techniques in [42], since we deal with a two dimensional problem. 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 [42] while the uniqueness is assured for only the 2D case. Secondly, there are many stable numerical approaches that we can use to resolve this equation derived from a classical example of fluid dynamics. Finally, we are able to implement this method efficiently and we have a sufficient theoretical framework in which we can rely on to understand the obtained results.
For a mathematical transparency of the Navier-Stokes equation (5), a convenient way is to consider its dimensionless form [18], which is obtained by introducing a reference length L * and a reference velocity V * , then we set where T * is the reference time defined as T * = L * V * . We note that concerning the choice of these parameters, we initialise the reference velocity by V * = 1 and the reference length L * by the length of the image domain. By a substitution into (5), we obtain for v ′ , P ′ , f ′ the following equation where R e is a dimensionless constant called the Reynolds number defined as In practice, the Reynolds number is always taken in the interval [2,100]. There is in fact a relation between the choice of the Reynold number and the deformations type. The Reynold term becomes higher when the deformations between template and reference images are larger. To avoid the complexity of the notations, we keep the same notation in (5), i.e. we substitute v ′ , P ′ , f ′ by v, p, f . Then, we rewrite the equation (9) as follows In the following section, we discuss the proposed discretization scheme for the Navier-Stokes equation (11), appropriate to the image registration problem.

Discretization.
To solve the Navier-Stokes equation (11), we use a finite volume discretization based on a staggered control volume illustrated in Fig. 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 meshes [21,36]. The use of a staggered mesh has a number of advantages for solving incompressible flow problems. In particular, this type of grid overcomes numerical problems associated with pressurevelocity coupling, which is the case when a collocated grid is used (the spurious pressure oscillations). Indeed, the representation of curved boundaries on a staggered Cartesian grid avoids some complexities which are inevitable when a non-staggered one is used. The control volumes for v 1 and v 2 are displaced with respect to the control volume for the continuity equation. In Fig. 1, 'P' denotes the node at which the partial differential equation is approximated, and 'E', 'W', 'N', 'S' are its neighbours. Cell faces 'e' and 'w' for v 1 and 'n' and 's' for v 2 are in the midway between the nodes.
Firstly, we rewrite the momentum equation of (11) in a differential form for both velocity components To solve this system of equations, 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 (12) over the control volume corresponding to v 1 , using the different point (E, P, n, s) shown in Fig. 1 and over the time interval from t to t + ∆t. Thus, using the fact that the velocities v 1 and v 2 are not depending on the vertical and horizontal component respectively under each volume control face, we have 1)∆t, x)), and vol = ∆x∆y. To simplify the study of this problem, we introduce some new entities defined as If we substitute these new entities in equation (14), we then have Also, the integration of the third continuity equation in (12) over the corresponding control volume and time gives which is equivalent to where (15) and (16), we find On the other hand, to calculate the entities (J E , J W , J N , J S ) in each grid, we need a approximation method. For this reason, we suppose that the flow is unidirectional, laminar with a constant pressure and without the external force (f = 0) on each grid [35]. Thus, the first equation of (12) in the ox direction becomes and the second one representing the equation in the vertical direction oy In order to respect the conservation of the velocities v 1 and v 2 , the quantities must be continuous across the common boundary of two control volumes. Therefore, we add to the equations (18) and (19) the boundary conditions which depend on the control volume corresponding to velocity components v 1 and v 2 , given by the following systems.
The first one (20) using the grid {P, E, s, n}, and the second one (21) using the grid {w, e, P, N }, see Fig. 1.
To calculate the entities (J E , J W , J N and J S ), we use the solution of the equations (20) and (21). The solution of equation (20) is given by the following expression In addition, the solution of (21) is given by Here v 1e and v 1n are only notations to make the difference between the two solutions of equations 20 and 21. To simplify the notations, we assume that where P e is called Peclet number, representing the local report on the boundary portion of the volume control for inertial forces and viscous forces. Equation (22) is now described as We can now calculate the expression of J E using the equation (25) injecting the expression of J E in (17), we have By the same way, we can find the other expressions ( given as follows

Finally, the discretization equation can be written as
where ; P e = (R e v 1 ) E ∆x, ; P n = (R e v 2 ) n ∆y, ; P s = (R e v 2 ) s ∆y, . The momentum equations for the other direction is obtained by the same way [35]. So the discretization for v 2 is given as where Since we have calculated the different momentum equations, we need to correct the pressure in each step to verify that ∇ · u = 0.
3.1. 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, last equation of (12).
Generally, the correct pressure p is obtained using where p ′ is the pressure correction. Then, we have to compute the new velocity components corresponding to the new corrected pressure p ′ . The two corresponding velocity corrections, denoted by v ′ 1 and v ′ 2 respectively, are obtained as Using the fact that v * 1 satisfies equation (27), if we replace the expression of v 1 in equation (27), we obtain an algebrical problem that is easy to solve following the same steps as in [35].
Finally, the correction pressure equation is given by where , To solve these equations, we use the algebraic equations to obtain the fields verifying the conservation equations [35]. We summarize the proposed method in the algorithm below. Once we have calculated v 1 and v 2 , the computation of Data: v * 1 , v * 2 and the pressure field p * , the Reynolds number 1 ≤ R e ≤ 1000. Result: The velocities v 1 , v 2 and the pressure p repeat the whole procedure until convergence is reached.
1. Solve the momentum equations (27) and (28)  Algorithm 1: The proposed finite volume algorithm for image registration the deformations u is given through a finite difference scheme. In the following subsection, we detail the computation of the deformations u.

Approximation of displacement u.
To compute the displacement u = (u 1 , u 2 ) from the associated velocities v 1 and v 2 , calculated through algorithm 1, we use a Euler scheme to resolve equation (6). 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 for the first step, to approximate the ∇u term, a centred finite difference approximation is used. For a fixed iteration k and for i = (i 1 , i 2 ), we have which represents the Jacobian matrix with 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 Our algorithm to compute the transformations u is finally summarized in Algorithm 2. 3. If the relative change of the distance measure is brought below a user supplied tolerance tol = 10 5 , the iteration is stopped, and v k+1 2 from the algorithm 1, return to step 1.

Algorithm 2:
The proposed algorithm to compute the displacement u 4. Results and discussion. In this section, we evaluate the performance of the proposed registration model. We have tested our algorithm on a benchmark of approximately one hundred MRI images. We present only five of them chosen with different deformations type and grey-level histogram. To measure the robustness of the proposed algorithm, we compare it with the classical implementation of the fluid registration model [13], for which the authors used a finite difference scheme. We approved that neither human participants nor animals are involved in this research.
In the first experiment, we choose firstly two images of human hands used for the first time in [1] and downloaded from the website 1 . In Fig. 2, we present the reference and template image.
Our aim is to find an optimal geometric transformation between corresponding images, and finally finding the reference image by registering the template one. In Fig. 3, we represent the registered image using the proposed approach compared with the obtained one by the fluid image registration proposed in [13]. In Fig. 4, we present the grid of deformations that are applied to the template image and the Jacobian determinants of deformation field. To show the evolution of the deformations with respect to the iterations, we present in Fig. 5 the application of the deformation to a rectangular grid with a rotation of 90 • degrees. We note that in some cases we apply a zoom to a region of the grid to better see the deformations. Finally, in the Fig. 6, we compare our algorithm with [13] using the error between the two images (reference and template). In the second example we take two MRI images of human brain downloaded from the web site 2 . We use the same process followed in the first example to present the next ones and keeping the same order of the figures. Fig. 7 shows the reference and template images of brain, Fig. 8 illustrates the obtained results, Fig. 9 shows deformation field for all the grid and the Jacobian determinants, while Fig. 10 shows the evolution of the deformations in a rectangular grid and the error between the deformed template and the reference images using the fluid image registration [13] and the proposed approach are presented in Fig. 11. In the third experiment, we apply the proposed approach to another slice of human brain. Again, we compare our model to the fluid image registration model [13] in Fig. 13 and, as in the previous tests, the deformation field and the Jacobian determinants in Fig. 14, while the evolution of the deformation in a rectangle grid and the error between template and reference image are presented, respectively, in Fig. 15 and Fig. 16. We follow the same steps as in previous tests for the human head image Fig. 17 and we present the results obtained in Fig. 18,  Fig. 19, Fig. 20 and Fig. 21. Finally, we show the obtained results in Fig. 23,  Fig. 24, Fig. 25 and Fig. 26 associated to the EPI slice image of human brain (Fig. 22) downloaded from the web site 3 . Note that the chosen parameters in the implementation of the proposed approach are presented in Table 2. For the Jacobian determinant, it is well known from calculus that the determinant of the Jacobian matrix det(J = I + ∇u) (here, I = diag(1, ..., 1) ∈ R d×d ) can be used to assess the invertibility of the deformation mapping x → x + u(x) as well as local volume change. Indeed, the Jacobian determinant value changes in each iteration and for all the tests. In general, the Jacobian determinant values are in the interval [0.005, 0.5]. For example, in the first test (human hands) min(det(J)) = 0.0105 and max(det(J)) = 0.4962, while in the fourth test (human brain) min(det(J)) = 0.0051 and max(det(J)) = 0.4988.
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 classical implementation of the fluid image registration [13]. However, we remark that the grey-levels of the deformed template are different from the original template ones in all tests approximately. There is in fact many reasons for the grey-levels loss. The widespread reason is due to the interpolation process, but in our problem, the main reason comes from the use of the Navier-Stokes equation (velocity-pressure (v, p) version). In fact, the NavierStokes equations can be simplified by introducing the stream function ψ satisfying ∇ ⊥ ψ = v and vorticity ω which satisfies ω = ∇ · v as dependent variables. Based on the work proposed by Bertalmio et al. [3], the image intensity I is related to the vorticity ω by a Poisson's equation ∆I = ω and according to this equation, I is more regular than ω (the smoothing effect of the Laplace operator). Indeed, the image intensity I defines the velocity field by v = ∇ ⊥ I. Due to this high regularity of the image I, intensities are perturbed. Then, a grey-level measure is needed to evaluate the ability of the proposed approach in avoiding this drawback. We then 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 [43]. In Table 3, the SSIM and PSNR values are calculated for all the used MRI images to confirm the robustness of the proposed approach in preserving the grey-level values compared to the fluid registration scheme and Thirion [13]. The best results are represented by a bold number. Always, the proposed method outperforms the others in terms of both PSNR and SSIM. On the other hand, we can also see some topological changes, especially in Fig. 13 where the two central disconnected components are merged, this is in fact related to the diffusion effect caused by the term 1 Re ∆v. We can also see the diffusion process in the figures that come after. To see the effect of the diffusion process, we have added two tests where we show the evolution of the registered image with respect to the Reynold's coefficient R e using images with same and different modalities (see Figs. 27 and 28). We precise that in all the above experiments, we have chosen the coefficient R e manually with respect to the low error between the deformed template and the reference image.       Figure 11. Difference error between template and reference images of human brain 1 using fluid registration (on the left) and the proposed approach (on the right).  Figure 16. Difference error between template and reference images of human brain 2 using fluid registration (on the left) and the proposed approach (on the right).          The full codes for the proposed model are implemented in MATLAB 2013. Typically, the execution of the main implemented programme requires an average of 2 ∼ 4 minutes on a 2.4 GHz Pentium Quad core computer for 256 × 256 grey-scale images; for the color and large-size images the computation time becomes more considerable. While the execution of the fluid image registration takes about 3 ∼ 7 minutes in the same conditions.

Conclusion.
A new numerical approach for the fluid image registration based on Navier-Stokes equation was introduced to developing the quality of the registration. A finite volume-based scheme was used to handle the discretization part errors. The performance of this new approach has been assured by the efficiency of the results, both visually and using the error criteria.
List of abbreviation. SIMPLE: Semi-Implicit Method for Pressure-Linked Equation.
Compliance with ethical standards.
• Funding: This research was entirely funded by institution of the authors.
• Conflict of interest: The authors declare that they have no conflict of interest.
• Neither human participants nor animals are involved in this research.