A Parallel Image Registration Algorithm Based on a Lattice Boltzmann Model

Abstract

Image registration is a key pre-procedure for high level image processing. However, taking into consideration the complexity and accuracy of the algorithm, the image registration algorithm always has high time complexity. To speed up the registration algorithm, parallel computation is a relevant strategy. Parallelizing the algorithm by implementing Lattice Boltzmann method (LBM) seems a good candidate. In consequence, this paper proposes a novel parallel LBM based model (LB model) for image registration. The main idea of our method consists in simulating the convection diffusion equation through a LB model with an ad hoc collision term. By applying our method on computed tomography angiography images (CTA images), Magnet Resonance images (MR images), natural scene image and artificial images, our model proves to be faster than classical methods and achieves accurate registration. In the continuity of 2D image registration model, the LB model is extended to 3D volume registration providing excellent results in domain such as medical imaging. Our method can run on massively parallel architectures, ranging from embedded field programmable gate arrays (FPGAs) and digital signal processors (DSPs) up to graphics processing units (GPUs).

1. Introduction

Lattice Boltzmann (LB) method was first introduced as a powerful numerical tool for solving partial differential equations (PDE) in computational fluid dynamics (CFD) [1,2]. The biggest difference between the LB method and the traditional numerical methods is due to the fact that the LB method is based on the microscopic description of physical systems. By simulating the microscopic behavior of physical systems, the macroscopic equation of the LB method can match certain PDEs. The idea of the LB method is to transpose discrete dynamics for simulating the macroscopic model described by a PDE, using densities of particles colliding and streaming on a regular lattice [3]. The density of each lattice node follows the same simple evolution rules, in which only each node of the lattice and its neighbors are involved. Intrinsically adapted to parallel calculation, the LB method is a good candidate for massive parallel computation.

Some approaches can be used to design a parallel algorithm for image processing. One approach is to try to convert a sequential algorithm to a parallel algorithm. If a serial algorithm already exists for a certain image processing problem, then the inherent parallelism in that algorithm may be implemented in parallel. Inherent parallelism is parallelism that occurs naturally within an algorithm, not as a result of any special effort on the part of the algorithm or machine designer. For example, in [3], the three-dimensional data registration is improved by parallel programming for solving the simultaneous localization and mapping problem of the robot motion model. The main idea is to parallize the NNS (the nearest neighbor search) computation, which is the first step of three-dimensional data registration. In [4], the iterative digital image correlation is divided into three sequential stages. Each stage is implemented by parallel computing. It should be noted that exploiting inherent parallelism in a sequential algorithm might not always lead to an efficient parallel algorithm. Another approach is to design a totally new parallel algorithm that is more efficient than the existing one. For example, the operator splitting (OS) method [5,6] and the lattice Boltzmann method. The OS method decomposes high-dimensional PDEs into several low-dimensional ones. Therefore, the difficulty of solving PDEs is reduced. Both the OS and LB methods are widely used in parallel computation fields. However, compared with the OS method, the LB method also offers higher computational efficiency and stability besides the parallelism [7]. Because of these advantages, variant LB methods have been proposed as an alternative to the conventional numerical solvers of PDEs in CFD, e.g., diffusion equations, convection diffusion equations, reaction diffusion equations, etc. Inspired by the application of PDEs in image processing, LB models have been applied as a novel method to the field of PDE-based image processing [8] in the past decade, e.g., image segmentation [9], image denoising [10] and image enhancement [11] by elaborating specific LB models.

Image registration is a key procedure for many image processing tasks in which the information is extracted by mapping two individual images (the reference and template images) obtained at different times or by different devices. For example, as a critical medical image fusion preprocessing, multiple medical images should be aligned before they are fused together in order to gather different important information into single one image. The imaging process is always affected by many complex factors. This will cause non-rigid deformation between two images to be registered. Therefore, the registration method is required to have the ability to deal with non-rigid deformation, especially in the fields of medical imaging, remote sensing and so on, which require high registration accuracy. Among the large number of non-rigid image registration methods, the method based on PDEs has become one of the most popular registration methods because of its easy extensibility, intuitive construction and the ability to express continuous change. Over the last three decades, many PDE-based models of physics phenomena have been successfully applied in both digital image processing and computer vision domains. The main idea is to map the image processing problem to some physical phenomenon. For example, the optical flow model is used to recover image motion at each pixel from spatio-temporal image brightness variations. As a consequence, the image registration is achieved.

At present, various registration methods have found important applications in many areas such as target recognition and tracking [12,13], image mosaic [14], image-guided surgery, image-guided radiotherapy [15], and so on. Generally speaking, registration methods can be divided into two types: feature-based registration and intensity-based registration. Rather than the total information contained in the images, feature-based image registration focuses on image features such as the edges, corners [16], etc. Therefore, the algorithm can be significantly reduced from the computational complexity point of view. In recent decades, many feature-based registration algorithms have been proposed [17,18,19]. However, feature-based methods require sufficient features, which should be properly detected and matched, otherwise the registration accuracy will decrease. For the cases low on features, a good choice is to use intensity-based registration, which is formulated as a PDE. The intensity-based image registration method could be adopted to deal with nonlinear transformation; the full information of the images is used, and does not require feature detection and matching steps. However, the main problem is that the intensity-based method is time consuming.

To override the disadvantages in both the feature and intensity-based methods, an LB model is proposed for intensity-based image registration. The main idea of our method consists in simulating the convection diffusion equation by establishing an LB model. Based on the efficiency of the LB model in dealing with the PDE, our model outperforms the classical nonlinear transformation methods, such as Wangs demons [20], optical flow model [21] and B-spline model [22], in terms of MSE (mean square error), RE (relevant error) and NMI (normalized mutual information). Because of the natural parallelism of the LB model [23,24], our model is efficient and parallel. It has the potential to be efficiently executed on massively parallel architectures on embedded field programmable gate arrays (FPGAs) and digital signal processors (DSPs) as well as graphics processing units (GPUs).

To explain our theoretical approach, this paper is structured as follows. Section 2 introduces our method including the basic theory of the LB model and how the two-dimensional (2D) and three-dimensional (3D) LB models for image registration are constructed step by step. Then, the experiments and discussion are presented in Section 3 and Section 4, and finally we reach the conclusion in Section 5.

2. Method

2.1. Basic Theory of the Lattice Boltzmann Model

The basic idea of the LB model is the construction of simplified discrete dynamics for simulating the macroscopic model, which is described by partial differential equations (PDEs) using densities of particles moving on a regular lattice. On the lattice, each node follows the same parallel evolution rules. The general lattice Boltzmann modeling consists in two steps: a translation step, during which particles move from node to node on a lattice, and a collision step, during which particles are redistributed at each node. The two steps can be represented by the discrete lattice Boltzmann equation (LBE). In this paper, a D2Q9 mode is used, where “D2” stands for “2 dimensions”, while “Q9” stands for “9 speeds”. The discrete LBE is given by [1]:

$$f_i\left(\stackrel{\rightharpoonup }{r}+{\sigma }_h{\stackrel{\rightharpoonup }{e}}_i,t+{\sigma }_t\right)-f_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)=-\frac{1}{\tau }\left(f_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)-f_i^{(0)}\left({\stackrel{\rightharpoonup }{r}}_i,t\right)\right),\stackrel{\rightharpoonup }{r}\in R^2$$

(1)

where $i=0,1,2,\dots ,8$ and $f_i\left(\stackrel{\rightharpoonup }{r}+{\sigma }_h{\stackrel{\rightharpoonup }{e}}_i,t+{\sigma }_t\right)$ is the particle density (also named density distribution function) on node $\left(\stackrel{\rightharpoonup }{r}+{\sigma }_h{\stackrel{\rightharpoonup }{e}}_i\right)$, at time $t+{\sigma }_t$ with time step ${\sigma }_t$ and designating lattice spacing ${\sigma }_h$ in direction ${\stackrel{\rightharpoonup }{e}}_i$:

$${\stackrel{\rightharpoonup }{e}}_i=\{\begin{array}{l}(0,0),i=0\hfill \\ (1,0),(0,1),(-1,0),(0,-1),i=1,3,5,7\hfill \\ (1,1),(-1,1),(-1,-1),(1,-1),i=2,4,6,8\hfill \end{array}$$

(2)

$${\sigma }_h=v{\sigma }_t$$

(3)

where $v$ is the velocity of the particles.

The right-hand side of Equation (1) represents the collision term, and the constant quantity $\tau$ is the single relaxation time, which controls the rate of approach to equilibrium. The equilibrium distribution functions $f_i^{(0)}\left(\stackrel{\rightharpoonup }{r},t\right)$ in direction ${\stackrel{\rightharpoonup }{e}}_i$ is only determined by the total density $\rho$ which is defined in terms of particle distribution functions:

$$f_i^{(0)}={\beta }_i\rho =\{\begin{array}{l}\rho /\left(9c\right),i=1,2,3,\dots ,8\hfill \\ \rho -8\rho /\left(9c\right),i=0\hfill \end{array}$$

(4)

where

$${\displaystyle \sum _{i=0}^8{f}_i=}{\displaystyle \sum _{i=0}^8{f}_i^{(0)}=}\rho$$

(5)

$f_i^{(0)}$ and $f_i$ must be positive, and $c$ is the regulator for the diffusion speed. Similar to [25], the macroscopic equation can be derived as follows.

$$\frac{\partial \rho }{\partial t}=\frac{{\sigma }_h{}^2}{18{\sigma }_t}div\left[\nabla \left(\left(\tau -0.5\right)\rho \right)\right]$$

(6)

Equation (6) is an anisotropic diffusion equation, also called Perona–Malik diffusion equation, which has an important application to the image denoising while preserving the image edges. In Section 2.2, a more general LB model will be developed in detail.

2.2. LB-Based 2D Image Registration Model

Suppose a reference image $R$ and a template image $T$ are given, the aim of image registration is to search a locale or global nonlinear optimal transformation $F$ which satisfies $F\left(T\right)=R$. This means the transformed template $T$ should match the reference $R$. In other words, we should find the displacement vector $\left(X,Y\right)$ which must satisfy $T\left(x-X,y-Y\right)=R\left(x,y\right)$ or $T\left(x-X,y-Y\right)\approx R\left(x,y\right)$. This problem can be expressed by the so-called sum of squared differences (SSD):

$$E=\frac{1}{2}{\displaystyle \underset{\mathsf{\Omega }}{\iint }{\left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]}^{2}dxdy}$$

(7)

where $\mathsf{\Omega }$ is the image spatial domain. Thus, the problem changes to find the optimal $\left(X,Y\right)$ to minimize energy $E$:

$$\underset{X,Y}{{\displaystyle MIN}}\left\{{\displaystyle \underset{\mathsf{\Omega }}{\iint }{\left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]}^{2}dxdy}\right\}$$

(8)

Equation (8) can be resolved by the variational method [26]. We can obtain the Euler–Lagrange equation of the variational problem (8) as follows [27]:

$$\{\begin{array}{l}\frac{\partial X}{\partial t}=\left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]\frac{\partial T}{\partial x}\\ \frac{\partial Y}{\partial t}=\left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]\frac{\partial T}{\partial y}\end{array}$$

(9)

However, Equation (9) is ill-posed because of the none-unique transformation from $T$ onto $R$. To make it well-posed, the regular term is always added to Equation (8), for example, a diffusion term:

$$\{\begin{array}{l}\frac{\partial X}{\partial t}=div\left(\frac{\nabla X}{\left|\nabla X\right|}\right)+\left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]\frac{\partial T}{\partial x}\\ \frac{\partial Y}{\partial t}=div\left(\frac{\nabla Y}{\left|\nabla Y\right|}\right)+\left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]\frac{\partial T}{\partial y}\end{array}$$

(10)

The first term of the right-hand side of the two equations in (10) is the regular term, while the second term is the SSD term, called image force term. Equation (10) has been proposed by Fisher et al. [27] for curvature-based image registration. In view of the PDE, Equation (10) comes to be a kind of convection diffusion equation which can be parallelly simulated by the LB method [28], but we need to modify the discrete LBE defined in Section 2.1 as follows:

$$\{\begin{array}{l}f{x}_i\left(\stackrel{\rightharpoonup }{r}+{\sigma }_h{\stackrel{\rightharpoonup }{e}}_i,t+{\sigma }_t\right)-f{x}_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)=-\frac{1}{\tau }\left(f{x}_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)-f{x}_i^{(0)}\left({\stackrel{\rightharpoonup }{r}}_i,t\right)\right)+\epsilon F_x\\ F_x=\frac{\alpha }{9}\left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]\frac{\partial T}{\partial x}\end{array}$$

(11)

and

$$\{\begin{array}{l}f{y}_i\left(\stackrel{\rightharpoonup }{r}+{\sigma }_h{\stackrel{\rightharpoonup }{e}}_i,t+{\sigma }_t\right)-f{y}_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)=-\frac{1}{\tau }\left(f{y}_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)-f{y}_i^{(0)}\left({\stackrel{\rightharpoonup }{r}}_i,t\right)\right)+\epsilon F_y\\ F_y=\frac{\alpha }{9}\left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]\frac{\partial T}{\partial y}\end{array}$$

(12)

where

$$\begin{array}{l}f{x}_i^{(0)}={\beta }_{i}X=\{\begin{array}{l}X/\left(9c\right),i=1,2,3,\dots ,8\hfill \\ X-8X/\left(9c\right),i=0,\hfill \end{array}\\ {\displaystyle \sum _{i=0}^{8}f{x}_i=}{\displaystyle \sum _{i=0}^{8}f{x}_i^{(0)}=}X\end{array}$$

(13)

and

$$\begin{array}{l}f{y}_i^{(0)}={\beta }_{i}Y=\{\begin{array}{l}Y/\left(9c\right),i=1,2,3,\dots ,8\hfill \\ Y-8Y/\left(9c\right),i=0,\hfill \end{array}\\ {\displaystyle \sum _{i=0}^{8}f{y}_i=}{\displaystyle \sum _{i=0}^{8}f{y}_i^{(0)}=}Y\end{array}$$

(14)

Without loss of generality, Equations (11) and (12) can be simplified as:

$$f_i\left(\stackrel{\rightharpoonup }{r}+v{\sigma }_t{\stackrel{\rightharpoonup }{e}}_i,t+{\sigma }_t\right)-f_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)=-\frac{1}{\tau }\left(f_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)-f_i^{(0)}\left({\stackrel{\rightharpoonup }{r}}_i,t\right)\right)+\epsilon F$$

(15)

Note that $\epsilon$ is a small value. $\alpha$ is a positive one, and $f_i^{(0)}$ is defined by Equations (4) and (5). Expanding the left side of Equation (15) to the Taylor series at $\left(\stackrel{\rightharpoonup }{r},t\right)$, we have:

$${\displaystyle \sum _{n=1}^{+\infty }\frac{{\sigma }_t{}^n}{n!}{\left(v{\stackrel{\rightharpoonup }{e}}_i\cdot \nabla +{\partial }_t\right)}^n{f}_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)}=-\frac{1}{\tau }\left(f_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)-f_i^{(0)}\left({\stackrel{\rightharpoonup }{r}}_i,t\right)\right)+\epsilon F$$

(16)

Introducing the Chapman–Enskog expansion [25] into the LB model, we have:

$$\{\begin{array}{l}{f}_i=f_i{}^{\left(0\right)}+\epsilon f_i{}^{\left(1\right)}+{\epsilon }^2{f}_i{}^{\left(2\right)}+\dots +{\epsilon }^n{f}_i{}^{\left(n\right)}\\ {\partial }_t=\epsilon {\partial }_{t1}+{\epsilon }^2{\partial }_{t2}\\ {\partial }_x=\epsilon {\partial }_{\alpha }\\ {\partial }_y=\epsilon {\partial }_{\beta }\end{array}$$

(17)

According to Equation (17), we obtain

$${\displaystyle \sum _i{f}_i}={\displaystyle \sum _i\left\{f_i{}^{\left(0\right)}+{\displaystyle \sum _{n=1}^{+\infty }{\epsilon }^n{f}_i{}^{\left(n\right)}}\right\}}=\rho$$

(18)

Considering formula (5), we conclude that

$${\displaystyle \sum _i{f}_i{}^{\left(n\right)}}=0,n\ne 0$$

(19)

Substituting (17) into (16) and comparing the coefficients of $\epsilon$ and ${\epsilon }^2$ on both sides, the following equations are obtained:

$${\epsilon }^1:{\sigma }_{t}v{\stackrel{\rightharpoonup }{e}}_i\cdot {\nabla }_{\alpha \beta }{f}_i^{(0)}+{\sigma }_t{\partial }_{t1}{f}_i^{(0)}=-\frac{1}{\tau }{f}_i^{(1)}+F$$

(20)

$${\epsilon }^2:{\sigma }_{t}v{\stackrel{\rightharpoonup }{e}}_i\cdot {\nabla }_{\alpha \beta }{f}_i^{(1)}+{\sigma }_t{\partial }_{t2}{f}_i^{(0)}+\frac{1}{2}{\sigma }_t{}^2{\left(v{\stackrel{\rightharpoonup }{e}}_i\cdot {\nabla }_{\alpha \beta }+{\partial }_{t1}\right)}^2{f}_i^{(0)}=-\frac{1}{\tau }{f}_i^{(2)}.$$

(21)

Equation (20) can be rewritten as

$$f_i^{(1)}=\tau F-\tau {\sigma }_{t}v{\stackrel{\rightharpoonup }{e}}_i\cdot {\nabla }_{\alpha \beta }{f}_i^{(0)}-\tau {\sigma }_t{\partial }_{t1}{f}_i^{(0)})$$

(22)

Summing up both sides of Equation (22) from $i$ = 0 to 8 leads to

$${\partial }_{t1}\rho =\frac{9F}{{\sigma }_t}$$

(23)

Summing up both sides of Equation (21) from $i$ = 0 to 8, we have

$${\sigma }_{t}v{\displaystyle \sum _i{\stackrel{\rightharpoonup }{e}}_i\cdot {\nabla }_{\alpha \beta }{f}_i^{(1)}}+{\sigma }_t{\displaystyle \sum _i{\partial }_{t2}{f}_i^{(0)}}+\frac{1}{2}{\sigma }_t{}^2{\displaystyle \sum _i{\left(v{\stackrel{\rightharpoonup }{e}}_i\cdot {\nabla }_{\alpha \beta }+{\partial }_{t1}\right)}^2{f}_i^{(0)}}=0$$

(24)

Substituting (22) into Equation (24) leads to

$${\partial }_{t2}\rho =\frac{2}{3c}{v}^2{\sigma }_t{\nabla }_{\alpha \beta }\left[\left(\tau -0.5\right){\nabla }_{\alpha \beta }\rho \right]$$

(25)

Thus, according to Equations (23), (25) and (17), we can obtain

$$\frac{\partial \rho }{\partial t}=\frac{2{\sigma }_h{}^2}{3c{\sigma }_t}div\left[\left(\tau -0.5\right)\nabla \rho \right]+\frac{9\epsilon }{{\sigma }_t}F$$

(26)

Let

$$\{\begin{array}{l}\rho =X\hfill \\ \tau =0.5+\frac{1}{\left|\nabla X\right|}\hfill \\ F=\left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]\frac{\partial T}{\partial x}\hfill \end{array}$$

(27)

We have

$$\frac{\partial X}{\partial t}=\frac{2{\sigma }_h{}^2}{3c{\sigma }_t}Div\left(\frac{\nabla X}{\left|\nabla X\right|}\right)+\frac{\alpha \epsilon }{{\sigma }_t}\left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]\frac{\partial T}{\partial x}$$

(28)

Similarly, we can obtain

$$\frac{\partial Y}{\partial t}=\frac{2{\sigma }_h{}^2}{3c{\sigma }_t}Div\left(\frac{\nabla Y}{\left|\nabla Y\right|}\right)+\frac{\alpha \epsilon }{{\sigma }_t}\left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]\frac{\partial T}{\partial y}$$

(29)

where ε and ${\sigma }_t$ are small values and might as well assume ${\sigma }_t=\epsilon$.

Therefore, Equations (28) and (29) can be rewritten as

$$\{\begin{array}{l}\frac{\partial X}{\partial t}=\frac{2{\sigma }_h{}^2}{3c{\sigma }_t}Div\left(\frac{\nabla X}{\left|\nabla X\right|}\right)+\alpha \left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]\frac{\partial T}{\partial x}\\ \frac{\partial Y}{\partial t}=\frac{2{\sigma }_h{}^2}{3c{\sigma }_t}Div\left(\frac{\nabla Y}{\left|\nabla Y\right|}\right)+\alpha \left[T\left(x-X,y-Y\right)-R\left(x,y\right)\right]\frac{\partial T}{\partial y}\end{array}.$$

(30)

Overall, Equations (15), (4) and (5) describe a simple LB model for computing each component of the displacement. In the 2D LB-based registration model, we need two such simple LB models to solve Equations (28) and (29), respectively.

2.3. LB-Based 3D Volume Registration

Our original model can easily be easily extended to a three-dimensional model for performing volume registration. A D3Q15 model is chosen and used in our study.

In the three-dimensional model, the displacement is a vector with three components $\left(X,Y,Z\right)$. Similar to the problem of 2D registration, three simple LB models are used to compute the displacement. According to the simplified LBE defined by Equation (15), the 3D LB model can be written as follows.

$$\{\begin{array}{c}f{x}_i\left(\stackrel{\rightharpoonup }{r}+{\sigma }_h{\stackrel{\rightharpoonup }{e}}_i,t+{\sigma }_t\right)-f{x}_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)=-\frac{1}{\tau }\left(f{x}_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)-f{x}_i^{(0)}\left({\stackrel{\rightharpoonup }{r}}_i,t\right)\right)+\epsilon Fx\\ f{y}_i\left(\stackrel{\rightharpoonup }{r}+{\sigma }_h{\stackrel{\rightharpoonup }{e}}_i,t+{\sigma }_t\right)-f{y}_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)=-\frac{1}{\tau }\left(f{y}_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)-f{y}_i^{(0)}\left({\stackrel{\rightharpoonup }{r}}_i,t\right)\right)+\epsilon Fy\\ f{z}_i\left(\stackrel{\rightharpoonup }{r}+{\sigma }_h{\stackrel{\rightharpoonup }{e}}_i,t+{\sigma }_t\right)-f{z}_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)=-\frac{1}{\tau }\left(f{z}_i\left({\stackrel{\rightharpoonup }{r}}_i,t\right)-f{z}_i^{(0)}\left({\stackrel{\rightharpoonup }{r}}_i,t\right)\right)+\epsilon Fz\end{array}$$

(31)

where

$$\{\begin{array}{c}Fx=\frac{\alpha }{15}\left[T\left(x-X,y-Y,z-Z\right)-R\left(x,y,z\right)\right]\frac{\partial T}{\partial x}\\ Fy=\frac{\alpha }{15}\left[T\left(x-X,y-Y,z-Z\right)-R\left(x,y,z\right)\right]\frac{\partial T}{\partial y}\\ Fz=\frac{\alpha }{15}\left[T\left(x-X,y-Y,z-Z\right)-R\left(x,y,z\right)\right]\frac{\partial T}{\partial z}\end{array}$$

(32)

Considering the particles have 15 streaming directions in the three-dimensional model (as shown in Figure 1), we define the equilibrium distribution functions as:

$$f{x}_i^{(0)}=\{\begin{array}{l}X/\left(15c\right),i=1,2,3,\dots ,14\hfill \\ X-14X/\left(15c\right),i=0\hfill \end{array}$$

(33a)

$$f{y}_i^{(0)}=\{\begin{array}{l}Y/\left(15c\right),i=1,2,3,\dots ,14\hfill \\ Y-14Y/\left(15c\right),i=0\hfill \end{array}$$

(33b)

$$f{z}_i^{(0)}=\{\begin{array}{l}Z/\left(15c\right),i=1,2,3,\dots ,14\hfill \\ Z-14Z/\left(15c\right),i=0\hfill \end{array}$$

(33c)

and

$$\{\begin{array}{c}{\displaystyle \sum _{i=0}^{14}f{x}_i=}{\displaystyle \sum _{i=0}^{14}f{x}_i^{(0)}=}X\\ {\displaystyle \sum _{i=0}^{14}f{y}_i=}{\displaystyle \sum _{i=0}^{14}f{y}_i^{(0)}=}Y\\ {\displaystyle \sum _{i=0}^{14}f{z}_i=}{\displaystyle \sum _{i=0}^{14}f{z}_i^{(0)}=}Z\end{array}$$

(34)

Directions ${\stackrel{\rightharpoonup }{e}}_i$ (as shown in Figure 1) are defined as:

$$\stackrel{\rightharpoonup }{e}=\left[{\stackrel{\rightharpoonup }{e}}_0,{\stackrel{\rightharpoonup }{e}}_1,{\stackrel{\rightharpoonup }{e}}_2,\dots ,{\stackrel{\rightharpoonup }{e}}_{14}\right]=\left[\begin{array}{lllllllllllllll}0\hfill & 1\hfill & -1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & -1\hfill & 1\hfill & -1\hfill & -1\hfill & 1\hfill & 1\hfill & -1\hfill \\ 0\hfill & 0\hfill & 0\hfill & 1\hfill & -1\hfill & 0\hfill & 0\hfill & 1\hfill & -1\hfill & -1\hfill & 1\hfill & -1\hfill & 1\hfill & -1\hfill & 1\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & -1\hfill & 1\hfill & -1\hfill & 1\hfill & -1\hfill & 1\hfill & -1\hfill & 1\hfill & 1\hfill & -1\hfill \end{array}\right].$$

(35)

As in Section 2.2, we obtain the macroscopic equations of the 3D LB model:

$$\{\begin{array}{c}\frac{\partial X}{\partial t}=\frac{2{\sigma }_h^2}{3c{\sigma }_t}Div\left(\frac{\nabla X}{\left|\nabla X\right|}\right)+\alpha \left[T\left(x-X,y-Y,z-Z\right)-R\left(x,y,z\right)\right]\frac{\partial T}{\partial x}\\ \frac{\partial Y}{\partial t}=\frac{2{\sigma }_h^2}{3c{\sigma }_t}Div\left(\frac{\nabla Y}{\left|\nabla Y\right|}\right)+\alpha \left[T\left(x-X,y-Y,z-Z\right)-R\left(x,y,z\right)\right]\frac{\partial T}{\partial y}\\ \frac{\partial Z}{\partial t}=\frac{2{\sigma }_h^2}{3c{\sigma }_t}Div\left(\frac{\nabla Z}{\left|\nabla Z\right|}\right)+\alpha \left[T\left(x-X,y-Y,z-Z\right)-R\left(x,y,z\right)\right]\frac{\partial T}{\partial z}\end{array}$$

(36)

3. Experiments and Discussion

3.1. Algorithm

According to Equations (4) and (10) to (14), the LB model can be applied to the image registration operation and achieves the same results as the curvature-based registration method proposed by Fisher et al. Take 3D registration for instance; the algorithm consists of the iterative following steps:

(1)

Initialize the values of $X$, $Y$ and $Z$ as zero;

(2)

Particles collision and streaming according to Equations (31) and (32);

(3)

Update $X$, $Y$ and $Z$ according to Equation (34);

(4)

Update equilibrium distribution functions according to Equation (33);

(5)

Go to step (2).

3.2. Experimental Results

To investigate the performance of our model, the registrations of 2D images and 3D volume data are shown in Figure 2. The 2D image registration testing was carried out on four types of test images: the registration of artificial images, computed tomography angiography (CTA) images, magnetic resonance (MR) images and video frames. The test images are given in Figure 2. The left column is the reference images and the right column, the template images. To measure the performance of the LB registration model, our method is compared with Wangs demons, B-spline registration model and optical flow registration model. In 3D volume registration, the MR volume data are random distorted as template data. In the experiments, the image/volume intensities are normalized between 0 and 1. All the experiments are implemented on a PC with Intel(R) Core(TM) i7-4702MQ CPU, 8G RAM, with MATLAB 2016.

3.2.1. 2D Image Registration

Parameter settings are presented in

Table 1. Parameters setting of the LB model.

$\mathit{\tau }$	${\mathit{\sigma }}_{\mathit{t}}$	${\mathit{\sigma }}_{\mathit{h}}$	$\mathit{\alpha }$	$\mathit{c}$
1	1	1	2.5	1.5

. Concerning the stability of the LB model, the relaxation factor $\tau$ should be larger than 0.5 because the diffusion coefficient should be positive, otherwise the diffusion function will be ill-posed. In the LB model, the lattice nodes correspond to the image pixels. Concerning the image coordinates system in which the pixel indices are integer values ranging from 1 to the length of the row or column, the grid spacing of the lattice should be 1. This means ${\sigma }_h=1$. Parameter $v$ and ${\sigma }_t$ should satisfy ${\sigma }_h=v{\sigma }_t$, otherwise the particles will stream into the place where the coordinate is not the integer in the next step. Therefore, we need to use interpolation technology to increase grid density. On one hand, the increased grid density can improve the computation accuracy, but on the other hand, the computational complexity is also increased. According to Equations (29) and (30), both parameter $c$ and ${\sigma }_t$ contribute to the diffusion coefficients. However, the influence of parameter ${\sigma }_t$ on the diffusion coefficients can be replaced by adjusting parameter $c$ which is used as a regulator of the diffusion coefficient. For simplicity, parameters $v$, ${\sigma }_t$ and ${\sigma }_h$ are set to be 1, while $\tau$, $c$ and $\alpha$ are chosen empirically. In the experiments, the values of $\tau$, ${\sigma }_t$, ${\sigma }_h$, $\alpha$ and $c$ are, respectively, 1, 1, 1, 2.5, and 1.5.

First, we subjectively compared the LB model with three popular methods as mentioned in Section 3.2 by exhibiting the errors $D=\left|T\left(x-X,y-Y\right)-R\right|$, and the displacement grid.

In experiment 1, a pair of synthetic checkboard images is chosen as test images (Figure 2a). The template image (right column of Figure 2a) is a distorted sine wave.

In experiment 2, the test images (Figure 2b) are extracted from a computed tomography angiography (CTA) image sequence. Note that the region of high gray-level displays high-intensity major vessels, which are a result of a weakened blood vessel wall.

In experiment 3, we register two MR images (Figure 2c). The template image (right column of Figure 2c) is generated by using a random distortion [29] on the reference image (left column of Figure 2c) acquired from the BrainWeb MRI Simulated Normal Brain Database [30].

In experiment 4, we register two video frames for tracking the moving objects (dragon dancer) recorded in our Chinese spring festival celebration.

The registration results of the four experiments are presented in Figure 3, in which the left column is the registration results, the middle column is the errors $D$, and the right column is the displacement grid.

Considering the difference between the registration results and reference images, the LB model produces smaller errors $D$. Notice that the evaluations in Figure 3 are mainly based on the subjective experience of the observer. For objective evaluation, four quantitative metrics can be used to evaluate the performance of the proposed model, namely, mean square error (MSE), relative error (RE), normalized mutual information (NMI) and execution time per iteration (ETPI). MSE is a popular metric to measure the absolute differences between two images. As a measure of precision, RE is the ratio of the errors $\left|T\left(x-X,y-Y\right)-R\right|$ before and after image registration. NMI evaluates the correlation between $T\left(x-X,y-Y\right)$ and $R$ in information theory. Higher similarity between two images implies smaller MSE and RE, and greater NMI. The formulas are as follows:

$$MSE=\frac{{\displaystyle \sum _y{\displaystyle \sum _x{\left|R-\widehat{T}\right|}^2}}}{N}$$

(37)

$$RE=\frac{{\displaystyle \sum _x{\displaystyle \sum _y\left|R-\widehat{T}\right|}}}{{\displaystyle \sum _x{\displaystyle \sum _y\left|R-T\right|}}}$$

(38)

$$NMI=\frac{H\left(R\right)+H\left(\widehat{T}\right)}{H\left(R,\widehat{T}\right)}$$

(39)

where

$$\{\begin{array}{l}H\left(R\right)=-{\displaystyle \sum _{x\in R}^{}p\left(x\right)\log p\left(x\right)}\\ H\left(R,\widehat{T}\right)=-{\displaystyle \sum _{x\in R}^{}{\displaystyle \sum _{y\in \widehat{T}}^{}p\left(x,y\right)\log p\left(x,y\right)}}\end{array}$$

(40)

In Equation (37), $\widehat{T}$ is the deformed version of template $T$. In Equation (39), $H\left(R\right)$ is the information entropy of image $R$, while $H\left(R,\widehat{T}\right)$ is the joint information entropy.

The experimental results are shown in

Table 2. Experimental results of experiments 1 to 4. Four quantitative metrics are used to evaluate the performance of the proposed model, namely, mean square error (MSE), relative error (RE), normalized mutual information (NMI) and execution time per iteration (ETPI).

Experiment	Method	MSE	RE	NMI	ETPI (s)
1	LB model	7.8517	0.0295	1.6129	0.0160
	Wangs demons	136.639	0.0889	1.4093	0.0237
	B-Spline	1035.5	0.2173	1.3303	73.8824
	Optical flow	278.0928	0.1348	1.3984	0.0253
2	LB model	71.1387	0.3127	1.1951	0.0067
	Wangs demons	73.8568	0.3286	1.1919	0.0157
	B-Spline	147.9227	0.5063	1.1495	6.1764
	Optical flow	69.3849	0.3051	1.2135	0.0186
3	LB model	38.6518	0.4376	1.2507	0.0091
	Wangs demons	115.7656	0.6347	1.2149	0.0193
	B-Spline	250.1226	1.0404	1.1789	13.6078
	Optical flow	43.3073	0.4484	1.2489	0.0158
4	LB model	318.5778	0.2767	1.2466	0.0213
	Wangs demons	1479.3864	0.4590	1.1907	0.0297
	B-Spline	922.9010	0.5030	1.1774	49.4118
	Optical flow	968.6714	0.4240	1.2263	0.0298

.

In Table 2, from the point view of registration accuracy, the LB model gains the smallest MSE and RE and the greatest NMI in experiments 1, 3 and 4. In experiment 2, there is only a tiny gap between the LB model and the optical flow registration model, which has the best registration accuracy. In all the four experiments, the LB model obtains the smallest index values of time cost per iteration. This implies that the LB model is relatively faster than the other three methods.

Figure 4 amplifies the displacement $\left(X,Y\right)$ given by the LB model in experiment 2.

Note that the region of high gray-level displays high-intensity major vessels, which are a result of a weakened blood vessel wall. The displacements around the vessel tends to shrink so that the template image T can be aligned with the reference. As shown in Figure 5a–c, the convergence of the LB model is illustrated by indices MSE, RE and NMI, which approach constants with the iterations.

3.2.2. 3D Volume Registration

To test the performance of our LB model for 3D volume registration, we used the BrainWeb MRI Simulated Normal Brain Database [30]. A random distortion was added to the reference data ($181\times 217\times 181$) with a maximum of 30 voxels of deformation, as shown in Figure 6a. Figure 6c shows the template data after registration, and Figure 6b, the reference data. Figure 7 gives the slices of original MR volume data, the distorted MR volume data and the registered results. The experiment takes 5132 seconds for 500 iterations. We can notice that there are many 3D LB models besidesD3Q15, e.g., D3Q7, D3Q19 and D3Q27. An LB model with more velocities can achieve higher calculation accuracy but is more time consuming. In the experiment, we used D3Q15 for considering the compromise between the calculation accuracy and time consumption.

4. Discussion

Parameter $\tau$ should be larger than 0.5 so that the diffusion coefficient is constrained to being positive. Otherwise, the LB model will be unstable, because a diffusion equation with a negative diffusion coefficient is ill-posed. Even when the diffusion coefficient $\tau -\frac{1}{2}$ is positive but very close to 0 (or $\tau$ very close to 0.5), the LB model is also unstable. As shown in Figure 8, the registration comes to be unstable when $\tau$ approaches 0.503 from 0.505.

In this study,$\tau =\frac{1}{2}+\frac{1}{\left|\nabla X\right|}$, this means that the gradients of $X$ and $Y$ should not be very large. However, the regular term $Div\left(\frac{\nabla X}{\left|\nabla X\right|}\right)$ and $Div\left(\frac{\nabla Y}{\left|\nabla Y\right|}\right)$ determine the smoothness of $X$ and $Y$ respectively, while the weight $\frac{2{\sigma }_h}{3c{\sigma }_t}$ controls the contribution of the regular term in Equations (25) and (26). In other words, a higher value of $\frac{2{\sigma }_h}{3c{\sigma }_t}$ leads to smaller gradients of $X$ and $Y$.

Parameter $\alpha$ and $c$ substantially affect the performance of our model. Parameter $\alpha$ and $\frac{1}{c}$ denote the SSD and the regular terms weight, respectively. For higher registration accuracy, we hope for large $\alpha$ and small $\frac{1}{c}$ values. However, this will result in a decrease in the regularization effect. With regard to the stability, we need relatively smaller $\alpha$ and larger $\frac{1}{c}$. Considering Equations (4) and (5), parameter $c\in \left[1,+\infty \right)$, the maximum of $\frac{1}{c}$ is 1. Parameter $\alpha \in \left(0,+\infty \right)$ because $\alpha$ is positive. In Figure 9, we give the relative errors with different $\alpha$ and $c$. Large values of $\alpha$ and $c$ can lead to the instability of the LB model, and a sharp increase in the RE value. Smaller values of $\alpha$ and $c$ can ensure the stability, but the accuracy is not the best.

In Figure 10, we present the relative errors $RE$ in domain $\left\{RE\left(\alpha ,c\right)|1\le c\le 6and0\le \alpha \le 5\right\}$. We find that the larger parameter $c$ and $\alpha$, the smaller the relative errors, but with a poorer stability. In domain $\left\{c,\alpha |c>6and\alpha >5\right\}$, the LB model is even unstable. Figure 11 illustrates an unstable example with c = 6.5 and α = 5.5. On the basis of the test, we suggest choosing $\left(c,\alpha \right)$ in domain $A=\left\{c,\alpha |1\le c\le 5.0and0\le \alpha \le 9.0\right\}$. Generally speaking, in domain $A$, if $\left(c,\alpha \right)$ is far away from (5.0, 9.0), the model is more stable, but less accurate, and vice versa.

5. Conclusions

Image registration is a key pre-procedure for high-level image processing. In this paper, an original parallel LB-based model for image registration is constructed. The main idea of the proposed method consists in simulating the convection diffusion equation by establishing an LB model. Experiments on CTA images and artificial images indicate that our model is much faster than classical methods and achieves accurate registration.

The stability is analyzed in the discussion, mainly focusing on the parameter $c$ and $\alpha$. High and $\alpha$ values can help to improve the registration accuracy, but may lead to instability. In this paper, we suggest $c\in \left[1,5.0\right]$ and $\alpha =\left[0,9.0\right]$. We also extend the LB model to 3D volume registration. Because of the high parallelism and high efficiency for programming and computation, the LB method is faster than the classical methods on a personal computer. Especially, our method can run on massively parallel architectures on large-scale systems such as embedded FPGAs, DSPs or GUPs.