A parallel 3D image segmentation method for Coronary CT Angiography

. Coronary artery (CA) disease is one of the major cardiovascular diseases that has been proved to be the leading cause of human death in world. In this paper, a 3D image segmentation algorithm based on Lattice Boltzmann (LB) is proposed for 3D CA segmentation. After investigating the behavior of Boundary Treatment schemes, a non-equilibrium extrapolation scheme is applied to keep the stability of the computation, and reduce the frequency of the re-initialization. The denoising and clipping method are also proposed for the segmentation refining. The segmentation result shows our LB model for 3D CAD segmentation is effective. More importantly, the LB model has natural parallelism. Our model can run on massively parallel architectures, ranging from inexpensive embedded FPGAs and DSPs up to GPUs.


Introduction
Coronary artery (CA) disease is one of the major cardiovascular diseases that has been proved to be the leading cause of human death in world.CA disease is caused by plaque buildup in the walls of the arteries that supply blood to the heart (called coronary arteries) and other parts of the body.Plaque buildup causes the inside of the arteries to narrow over time, which can partially or totally block the blood flow.Because of the fact that coronary artery disease often develops over decades, patients might not notice a problem until they have a significant blockage or a heart attack.
For those at risk of CA disease, Coronary CT angiography (CTA) is the major method to examine CA disease.CTA is convenient, noninvasive and effective to examine the stenosis or blockage of blood vessels.For this reason, CTA has become an important tool in screening, diagnosis and treatment of cardiovascular diseases.For many years, many researchers have carried out extensive research on quantitative description of cardiovascular disease and vascular dynamics based on CTA images.The key step of these studies is the two-dimensional (2D) or three-dimensional (3D) segmentation of CTA images.
However, there are challenges in CTA image segmentation, such as complex structure of coronary artery, variant vessel diameter, strong noise especially the motion artifacts caused by heart beating and low intensity contrast between the coronary artery and surrounding tissues.Therefore, CTA segmentation tends to segment CAD and the surrounding tissues as a whole which results in over segmentation.
In past few decades, the commonly used image segmentation methods are based on histogram threshold, active contour, region growing, artificial neural network, and etc.The methods on histogram threshold and region growing are prone to oversegment in the low intensity contrast region, while the artificial neural network needs huge workload on sample collection and labeling.Considering the image dimensions , compared to the traditional methods based on 2D image, 3D image segmentation is more accurate and reliable.This is because 3D image segmentation can directly use the depth information to separate the foreground and background, and can easily extract the target.3D image segmentation has always been a classic challenge in the field of image processing [1][2] [3].The curve evolution methods based on active contour [4][5] and level sets [6][7] model are widely concerned, for example the CV model [8] proposed by Chan and Vese based on Munford-Shah segmentation model [9].However, because of the huge amount of data in threedimensional images, the computation of solving CV model is very large, which consumes a lot of time.
In this paper, a 3D image segmentation algorithm based on Lattice Boltzmann (LB) model [10] is proposed for 3D CAD segmentation.A a novel numerical tool, LB model is simple and efficient.More importantly, the LB model has natural parallelism.Our model can run on massively parallel architectures, ranging from inexpensive embedded FPGAs and DSPs up to GPUs.

CV Model
A CV model can be expressed in the summation of three energy terms. Where In formula (1), Curve C is a curve that separates image domain  into two parts: Embedding curve C as the zero level set of distance function  in domain  , Energy E can be rewritten as Where H is Heaviside function: The gradient of the Heaviside function H is:  is so called Dirac function.Thus formula (1) is converted from line integral to surface integral in entire image domain.
Minimize formula (2) leads to the following partial differential equation.
The first term in the right side of equation ( 3) is a regular term that make equation ( 3) to be well-posed.The Dirac function is to limit the evolution in a narrow band nearby curve C. Practically, the regular term can be replaced by a simpler curvature term: where,  H ,   are the regularized approximations of the Heaviside function H and  respectively.

The basic lattice Boltzmann model
The basic idea of the LB method is the construction on simplified discrete dynamics for simulating the macroscopic model which are described by partial differential equations (PDE) using densities of particles moving on a regular lattice.The general Lattice Boltzmann modeling consists of two steps: a translation step in which particles move from node to node on a lattice and a collision step in which particles are redistributed at each node.The two steps can be represented by the discrete lattice Boltzmann equation (LBE).In this section, a D2Q9 model [11] is introduced, where "D2" stands for "2 dimensions", while "Q9" for "9 speeds".The discrete LBE is given by: where  is the relaxation factor, and where  is the velocity of the particles.
The right hand side of the equation ( 6) represents the collision term, and the constant quantity  is the single relaxation time, which controls the rate of approaching to equilibrium.The equilibrium distribution functions in direction i e  is only determined by the total density  which is defined in terms of the particle distribution functions: where eq i f and i f must be positive, and c is the regulator for the diffusion speed.Formula ( 6) -( 10) described a D2Q9 (2 dimensions 9 directions) LB model.Using Chapman-Enskog expansion and Taylor developing both sided of equation ( 6), we can obtain the macroscopic equation of the LB model in the form of PDE as follows (The derivation can be referred to [10]): 4

LB based two-dimensional CV model
Equation ( 11) is an anisotropic diffusion model while the relaxation factor  is variable value.In our model,  is defined as: . Therefore, the global description of the LB model introduced in section 3 is written as: Compared with equation ( 4), a new LB model is modified based on formula (6): where F is diffusion source which is defined as: Similar to section 2, we can obtain the macroscopic equation:

LB method based three-dimensional CV model
The two-Dimensional CV model can be directly modified to Three-Dimensional.
The equilibrium distribution function is as follows: Direction i e  (as shown in Fig. 1) is defined as: ( The macroscopic equation is shown as following: For simplicity, h  and t  are set to 1: 5 c do not changes.

Parameter settings
It should be noted that the value of the distance function  must be greater than 0 to ensure the stability of the LB model.That is because the LB will get unstable if the particle density is negative.For simplicity, we initialize the distance function as following formula: where to equation (11), relaxation factor  should be larger than 0.5, otherwise equation (11) will get ill-posed.
Diffusion coefficient c 3 2 not only determines the contribution of the regular term to the curve evolution,but also affects the behavior of the segmentation.Obviously, A smaller c means increased rigidness of the curve, but more stability of the model, and vice versa.In Fig. 2(a)-(c), the CAD segmentation with c=50, 5, 1 is presented.With the decreasing of c, some small tissues are missed in the segmentation, and the segmentation is of lower accuracy.This is mainly because the curve are getting more rigid, in other words, the curvature of the curve is decreased.However, a larger c will reduce the contribution of the regular term, which will cause the instability of the model.In the experiments, we suggest set c to 10.

Boundary treatment
Periodic scheme and bounce-back scheme [12] are two common methods to deal with simple boundaries in lattice Boltzmann model.However, both the schemes have poor universality, for example, they can not deal with the boundary containing gradient information.In the experiments, these two boundary treatment schemes may degrade the distance function quickly, as a consequence, re-initialization of the distance function should be implemented frequently.
In our research, we use an non-equilibrium extrapolation scheme [13] which is more stable and easy to implement.As shown in Fig. 3 where To investigate the performance of the non-equilibrium extrapolation scheme, two indexes are used, such as the mean value and variance of gradient   denoted by MVG and VG respectively.It is worth noting that, we set F to 0 so that the distance function  is not effected by the diffusion source which can change the distance function dramatically in the image segmentation.The indexes MVG and VG over iterations based on nonequilibrium extrapolation and periodic scheme are plotted in Fig. 4 and Fig. 5 respectively.Fig. 4 illustrate that, compared to the non-equilibrium extrapolation scheme, the MVG based on periodic scheme decreased much rapidly, and almost equals to 0 after about 2700 iterations.This means the distance function  has degraded to a constant function.In Fig. 5, the VG based on periodic scheme increase to the maximum after 700 iterations, and then, decreased to 0 after 2700 iterations.Actually this is due to the fact that, in the first few iterations, the distance function degrades to 0 only nearby the image boundary.And then, because of the diffusion procedure, the degradation expands quickly to the whole image domain until the distance function equals to constant 0. This indicates that the periodic scheme is more unstable which need to re-initialize the distance function more frequently than non-equilibrium extrapolation scheme.

Conclusion
In this paper, we proposed a 3D CA segmentation in coronary CT Angiography.The algorithm is simple and efficient.More importantly, the LB model has natural parallelism.We investigate the behavior of the affect of Boundary Treatment to the distance function.
The value and variance of   over iterations show different boundary treatments may cause different effects on the degrade speed of the distance function.In our research, the non-equilibrium extrapolation scheme is applied to keep the stability of the computation, and reduce the frequency of the re-initialization.The segmentation result shows our LB model for 3D CAD segmentation is effective.The algorithm is simple and efficient.More importantly, the LB model has natural parallelism.Especially, our method can run on massively parallel architectures, ranging from inexpensive embedded FPGAs and DSPs up to GPUs.
This work was supported by Jiangsu industry university research cooperation project (BY2020650).

1 (inside C) and 2 
(outside C).I is the image to be segmented.Variable Fig. 1.Direction i e  .

. 6 ITM
the contour as a circle centered at point The distance function is positive outside the contour and negative inside.Parameter B should satisfy radius B  , so that distance function radius   .According Web of Conferences 45, 01036 (2022) CSCNS2021 https://doi.org/10.1051/itmconf/20224501036

,
, assume points E, O and A locate on the boundary of the image; B,C,D locate inside the image while F,G,H outside the image.Take point O for instance, the particle density   t O f i , is divided into two parts: and non-equilibrium part

) 7 ITMFig. 3 .
Fig. 3. Points E, O and A locate on the boundary of the image; B,C and D locate inside the image while F, G and H outside the image.

10 ITMFig. 9 .
Fig. 9.The first refining.a) The plane used to clip the segmentation; b) the processing result after clipping and denoising.

Table 1 .
Parameter setting in the refinings of our experiment.