S Balla-Arabé
email: balla_arabe_souleymane@ieee.org

Architecture-Driven Level Set Optimization: From Clustering to Sub-pixel Image Segmentation < 1

Keywords: Level set method, image segmentation, graphics processing units, hybrid CPU-GPU architecture

Thanks to their effectiveness, Active contour models (ACMs) are of great interest for computer vision scientists. The level set methods (LSMs) refer to the class of geometric active contours. Comparing with the other ACMs, in addition to sub-pixel accuracy, it has the intrinsic ability to automatically handle topological changes. Nevertheless, the LSMs are computationally expensive. A solution for their time consumption problem can be hardware acceleration using some massively parallel devices such as graphics processing units (GPUs). But the question is: which accuracy can we reach while still maintaining an adequate algorithm to massively parallel architecture? In this work, we attempt to push back the compromise between, speed and accuracy, efficiency and effectiveness, to a higher level, comparing with state-of-the-art methods. To this end, we designed a novel architecture-aware hybrid CPU-GPU LSM for image segmentation. The initialization step, using the well-known k-means algorithm, is fast although executed on a CPU, while the evolution equation of the active contour is inherently local and therefore suitable for GPU-based acceleration. The incorporation of local statistics in the level set evolution allowed our model to detect new boundaries which are not extracted by the used clustering algorithm. Comparing with some cutting-edge LSMs, the introduced model is faster, more accurate, less subject to giving local minima, and therefore suitable for automatic systems. Furthermore, it allows two-phase clustering algorithms to benefit from the numerous LSM advantages such as the ability to achieve robust and sub-pixel accurate segmentation results with smooth and closed contours. Intensive experiments demonstrate, objectively and subjectively, the good performance of the introduced framework both in terms of speed and accuracy.

I. INTRODUCTION

N computer vision and pattern recognition, image segmentation [START_REF] Chen | Multiple Kernel Fuzzy C-means based Image Segmentation[END_REF][START_REF] Fu | Symmetry Constraint for Foreground Extraction[END_REF][START_REF] Balla-Arabé | A fast and robust level set[END_REF][4][5][6][START_REF] Balla-Arabé | Multi-Kernel Implicit Curve Evolution for Selected Texture Regions Segmentation in VHR Satellite Images[END_REF] is a major process by which a given image is partitioned into a number of meaningful and homogeneous regions, such that the union of any two neighborhood regions yields a heterogeneous segment. Due to the fact that there is no a general framework which is effective for all kinds of images, the task is non-trivial and more challenging in presence of noise. One should choose the proper method according to the characteristics of the image in hand.

In recent years, optimization methods have attracted much attention as powerful and natural image segmentation tools. Their basic principle is to achieve segmentation by minimizing a given energy function designed from the image information. Optimization methods can be roughly classified into two important classes: spatially discrete and spatially continuous representations. In spatially discrete approaches, the image pixels are usually considered as the nodes of a graph, and the aim of segmentation is to find cuts of this graph which have a minimal cost [START_REF] Boycov | Fast Approximate Energy Minimization Via Graph Cuts[END_REF][START_REF] Malik | Contour and Texture Analysis for Image Segmentation[END_REF].

Active contour models belong to the spatially continuous approaches: the segmentation of the image plane is considered as a problem of infinite-dimensional optimization. The main idea is to evolve a given curve in the direction of negative energy gradient by means of an appropriate partial differential equation. The level set method (LSM) [START_REF] Papandreou | Multigrid Geometric Active Contour Models[END_REF][START_REF] Wang | A Nonlinear Adaptive Level Set for Image Segmentation[END_REF][START_REF] Li | Implicit Active Contours Driven by Local Binary Fitting Energy[END_REF][START_REF] Li | A Level Set Method for Image Segmentation in the Presence of Intensity Inhomogeneities With Application to MRI[END_REF] designates the class of ACMs which uses the Eulerian framework, i.e., the geometric representation of the evolving curve, instead of the parametric one, i.e., the Lagrangian framework [START_REF] Zhu | Region Competition: Unifying Snakes, Region Growing, and Bayes/MDL for Multiband Image Segmentation[END_REF][START_REF] Brun | A Nonconservative Lagrangian Framework for Statistical Fluid Registration-SAFIRA[END_REF][START_REF] Nakhmani | Self-crossing Detection and Location for Parametric active contours[END_REF][START_REF] Kass | Snakes: Active Contour Models[END_REF]. Fig. 1 displays the organization of the ACMs.

Comparing with parametric active contours, the LSM presents more advantages, such as the ability to easily handle complex shapes, and topological changes. Furthermore, it allows a straightforward passage from two-dimensional (2D) to 3D space. The original idea of the LSM stems from the Hamilton Jacobi approach, i.e., a time-dependent equation for a moving surface [START_REF] Osher | Fronts Propagating and Curvature Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulation[END_REF][START_REF] Malladi | A Topology Independent Shape Modeling Scheme[END_REF][START_REF] Caselles | On Geodesic Active Contours[END_REF]. In space, the LSM aims to evolve a given curve toward its interior or exterior normal until defining the boundary of the object of interest. The curve evolution is driven by the level set equation (LSE) which is a nonlinear partial differential equation formulated as follow

( ),      Vk t     ( 1
)
where  is the level set function (LSF), V the velocity field which guides the active curve toward the needed boundaries,

 and  are user-controlling parameters and k is a nonlinear curvature term expressed as

.      k   (2)
Two main approaches are usually used to stop the evolving curve on the contours of the desired object. The first uses an edge indicator depending on the gradient of the image, as in classical snake and active contour models [START_REF] Gao | A Relay Level Set Method for Automatic Image Segmentation[END_REF][START_REF] Kichenassamy | Gradient Flows and Geometric Active Contours Models[END_REF][START_REF] Cremers | A Review of Statistical Approaches to Level Set Segmentation: Integrating Color, Texture, Motion and Shape[END_REF][START_REF] Delgado-Gonzalo | CYB-E-2015-04-0460: Architecture-Driven Level Set Optimization: From Clustering to Sub-pixel Image Segmentation < 10 "Snakes with an Ellipse-Reproducing Property[END_REF][START_REF] Li | Distance Regularized Level Set Evolution and its Application to Image Segmentation[END_REF]. These models have a lower computational complexity and are suitable to parallel architecture, but they are sensitive to noise and initialization.

Architecture-Driven Level Set Optimization: From Clustering to Sub-pixel Image Segmentation The second uses some regional attributes to stop the evolving curve on the real boundary. These models are more robust against noise and can detect objects with weak edges. One of the most known region-based level set algorithms was proposed in [START_REF] Chan | Active Contours Without Edges[END_REF] where Chan and Vese introduced a level set formulation to minimize the Mumford and Shah functional [START_REF] Mumford | Optimal Approximations by Piecewise Smooth Functions and Associated Variational Problems[END_REF]. They converted the problem into a mean curvature flow problem just like the active contours. In [START_REF] Karantzalos | A region-based level set segmentation for automatic detection of man-made objects from aerial and satellite images[END_REF], the authors used this approach in the field of remote sensing imagery for automatic detection of man-made objects (roads, buildings, etc.) from aerial and satellite images. In [START_REF] Bazi | Unsupervised Change Detection in Multispectral Remotely Sensed Imagery with Level set Methods[END_REF], Bazi et al. presented an unsupervised change detection method for multispectral images by designing a region-based energy functional using the difference image. The level set framework was then used to minimize the designed functional.

The region-based level set algorithms give better results than those of the classical active contours because the stopping term did not depend on the gradient of the image, reducing the dependence on strong edges and improving the robustness against noise. However, most of the times they cannot deal with intensity inhomogeneities, and are definitely not suitable for parallel programming, because they are not local. Furthermore, in order to solve the LSE, most of these methods suggest the use of some computationally expensive finite difference, finite element or finite volume approximations and an explicit computation of the curvature [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF][START_REF] Vese | A Multiphase Level Set Framework for Image Segmentation Using the Mumford and Shah Model[END_REF].

In our previous work [START_REF] Balla-Arabé | GPU Accelerated Edge-Region Based Level Set Evolution Constrained By 2D Gray-scale Histogram[END_REF] and [START_REF] Balla-Arabé | Level Set Region Based Image Segmentation Using Lattice Boltzmann Method[END_REF], in order to handle the problem of computational expense, we used the lattice Boltzmann method (LBM) as an alternative approach for solving the LSE. Comparing with [START_REF] Chan | Active Contours Without Edges[END_REF] for example, these approaches give better results in terms of efficiency and accuracy. The problem of time consumption is better handled because the curvature is implicitly computed and the algorithm is simple. But these introduced methods have some limits. In [START_REF] Balla-Arabé | GPU Accelerated Edge-Region Based Level Set Evolution Constrained By 2D Gray-scale Histogram[END_REF], we used regional attributes of a given pixel to design an Unsigned Pressure Force (UPF) which acted as a stop function for the active curve. This method is greatly dependent to the initialization which limits its use in automatic systems. Furthermore, the UPF is not local and therefore not suitable to parallel architecture-based acceleration. In [START_REF] Balla-Arabé | Level Set Region Based Image Segmentation Using Lattice Boltzmann Method[END_REF], we used both local and regional statistics to design an energy function which had to be minimized using the level set framework, the method was more effective and robust against initialization than the one introduced in [START_REF] Balla-Arabé | GPU Accelerated Edge-Region Based Level Set Evolution Constrained By 2D Gray-scale Histogram[END_REF]. But as in [START_REF] Balla-Arabé | Level Set Region Based Image Segmentation Using Lattice Boltzmann Method[END_REF], we still need to compute the mean values of pixels intensity inside and outside the active curve at each iteration, which increases the communication between the processors of any massively parallel device like GPU, and therefore slowdown the method.

In this work, we introduce a novel architecture-driven framework which allows obtaining a highly effective and fully local level set equation. Therefore, the proposed algorithm is adequate to parallel hardware-based acceleration, while allowing the achievement of promising segmentation results. In addition, the proposed method is a general framework which allows to considerably improving the results of most two-phase clustering methods such as k-means, the fuzzy c-means and their derivatives. By using the two intensity means output of the clustering algorithms and some local statistics of the image, we design an effective and inherently local energy functional. A fully local LSE is then obtained using the gradient descent methods. This strategy effectively allows combining, in a straightforward manner, the advantages of clustering technics and level set methods. On the one hand, the speed and the low complexity of some clustering algorithms, such as k-means, make the methods faster in comparison to some state of the art LSMs; on the other hand, the method has the advantages of the LSMs such as the ability to 1-achieve sub-pixel accuracy [START_REF] Chan | Active Contours Without Edges[END_REF]; 2-allow the incorporation of various prior knowledge, for example, shape and intensity distribution in order to achieve more robust segmentation results [START_REF] Chen | Using prior shapes in geometric active contours in a variational framework[END_REF]; 3-provide smooth and closed contours which are ineluctable for further applications such as shape analysis and recognition [START_REF] Li | Minimization of Region-Scalable Fitting Energy for Image Segmentation[END_REF].

Although, most of two-phase clustering methods can be used, in this paper we use the output of the k-means algorithm, which has been initialized using the interior and the exterior intensity means of the initial active curve. The k-means algorithm, which is intrinsically fast and non-local, is executed on the CPU, while the level set evolution is executed on an NVIDIA GPU. The major contributions of the present work can be summarized as follows: Intensive experiments, using the Berkeley segmentation dataset BSDS300 [START_REF] Martin | A Database of Human segmented Natural Images and its Application to Evaluating Segmentation Algorithms and Measuring Ecological Statistics[END_REF], demonstrate, subjectively and objectively, the effectiveness and efficiency of the proposed segmentation framework.

The remainder of the paper is organized as follows. Section 2 introduces the LBM. The formulation of the proposed method is given in Section 3. Section 4 presents the experimental results. The final section is the conclusion.

II. BACKGROUND

This section gives a general idea of the lattice Boltzmann model. The LBM was first designed to simulate Navier-Stokes equations for an incompressible fluid [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond Numerical Mathematics and Scientific Computation[END_REF][START_REF] Zhao | Lattice Boltzmann Based PDE Solver on the GPU[END_REF][START_REF] He | Lattice Boltzmann Model for Incompressible Navier-Stokes Equation[END_REF]. Its evolution equation is

( , 1) ( , ) ,             i i i coll f f r e t f r t t (3)
which can be decomposed into two steps

Collision : ( , ) ( , ) ,        coll i coll f f r t f r t t (4) Streaming : ( , 1) ( , ),       coll ii f r e t f r t (5)
where i f is the particle distribution function and  r a spatial variable. In this paper,    coll ft is the Bhatnager-Gross-Krook (BGK) collision model [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-component Systems[END_REF][START_REF] Buick | Gravity in a lattice Boltzmann model[END_REF][42][43], with a body force 

F , expressed as 2 1 .[ ( , ) ( , )] . ,              eq i i i coll fD f r t f r t F e t bc  (6)
where D is the grid dimension, b the link number at each grid point, c the length of each link which is set to 1 in this paper and  represents the relaxation time. eq i f is the local Maxwell-Boltzmann equilibrium particle distribution function expressed in its continuous form as

2 32 () (2 ) exp{ }, 2      eq u f RT RT   (7)
where   is the particle velocity and  u the macroscopic velocity. The equilibrium distribution can be expressed in discrete form as follows when modeling typical diffusion phenomenon, ( ) ,

  eq i i i i f A with f    (8) 2 2
, f e 00 , fe where  is the macroscopic fluid density. By performing the Chapman-Enskog expansion [START_REF] Chapmann | The Mathematical Theory of Non-Uniform Gases: An account of the kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases[END_REF], the following diffusion equation can be recovered from the lattice Boltzmann evolution equation [START_REF] Zhao | Lattice Boltzmann Based PDE Solver on the GPU[END_REF],

( ) .      div F t   (9)
Substituting  by the signed distance function  in Eq. ( 9), the LSE can be recovered. In our model we use the D2Q5

( 2  D , 5  b
) LBM lattice structure. Fig. 2 displays a typical D2Q5 model where each link has its velocity vector ( , )  i e r t . The body force F acts as the link with image data for the LBM solver.

Reference [START_REF] Chen | Cellular Automata Based Level Set Method for Image Segmentation[END_REF] used another approach to perform the level set image segmentation. Eq. ( 3) is the general evolution equation of the LBM; however in level-set-based image segmentation a stop function ()  gr is necessary to stop the evolving curve or surface at the boundaries of the object. In order to introduce the stop function into the LBM, the authors considered a medium between the nodes of the lattice. The particles can pass through the medium with a possibility of ()  i gr, and will be pushed back to where they were with a possibility of 1

( )   i gr . The LBM evolution equation is modified as 1 , 1 )[ , [ ( , , ] ] (1 ) ( ( ) ( ) ( )
( ) ) ( ) , ,                  eq i i i i i i i i i f r e t g r f r t f r t f r t g r f r e t   ( 10
)
where  is the convection coefficient. The macroscopic fluid density  is set as a signed distance function. More details about this approach can be found in [START_REF] Chen | Cellular Automata Based Level Set Method for Image Segmentation[END_REF].

III. THE PROPOSED HYBRID LEVEL SET FRAMEWORK FOR TWO-PHASE IMAGE SEGMENTATION

This section is dedicated to the conception and analysis of the proposed algorithm. We first design a novel energy functional, which will be minimized using the level set framework. Let 𝜙: Ω → ℝ be an LSF defined on a domain  . In all our work  is a signed distance function positive inside the zero level contour  and positive outside. The introduced energy functional is defined as follow > CYB-E-2015-04-0460: Architecture-Driven Level Set Optimization: From Clustering to Sub-pixel Image Segmentation < 4

( ) ( ) ( ),  main Regularization E E E    (11)
where ()

main E
 is the main energy term which contains the regional information of the image to be segmented, ()

Regularization E
 is the constraint on the smoothness and the length of the evolving contour.

By using the gradient descent method, the LSE can be recovered from the above defined energy function ,



 

E t   (12)
where  E  is the Gâteaux derivative [START_REF] Aubert | Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations[END_REF] of  . According to Eq. ( 11), Eq. ( 12) is equivalent to the following evolution equation

() .          Regularization main E E t    (13)
A. Design of ()

main E 
The main energy term ()

main E
 has been designed so that:

1-It is inherently local in order to suit to massively parallel architecture, 2-Its minimization yields to an effective and accurate image segmentation even in the presence of outliers.

()

main E  is therefore defined as 2 11 2 2 2 ( ) (1 exp{ ( ( ) ) }) ( ) (1 exp{ ( ( ) ) })(1 ( )) ( ), main Outliers E I x P H dx I x P H dx E                       (14)
Where Outliers E is the energy to be minimized in order to make the method robust against outliers, H is the Heaviside function,  with respect to  can therefore be formulated as

x
2 1 2 2 ( )[ (1 exp{ ( ( ) ) }) (1 exp{ ( ( ) ) })] . main in Outliers out E I x P E I x P                     (15)
In [START_REF] Balla-Arabé | GPU Accelerated Edge-Region Based Level Set Evolution Constrained By 2D Gray-scale Histogram[END_REF], we had defined an efficient energy term based on a 2D gray-scale histogram which can make the LSM more robust by avoiding the evolving curve stopping at the pixels which are more likely to be outliers. This energy term was defined as follows

.

        m I x y H dxdy H dxdy   (18)
In this work, inspired by the energy term defined in Eq. ( 16), we define

Outliers E as follows 12 2 1 2 1 1 2 2 2 ( ) ( , ) exp{ ( )} ( ) ( , ) (1 exp{ ( ( ) ) })
(1 exp{ ( ( ) ) }),

Outliers mean

E A P P I I H dx

with A P P I x P

I x P                    (19)
A straightforward analysis of

Outliers E
shows that the effect of it minimization is to make the present method robust against outliers by avoiding the evolving curve stopping at pixels for which the intensity value is far different from their local intensity average, that is, high gradient isolate pixels. The derivative of

Outliers

E

with respect to  can therefore be written as follows 12 [ ( , )exp{ (

)}] ( ),

Outliers mean

E A P P I I           ( 20
)
By adding this term in the main energy, we finally get the following equation

12 ( ) ( , )[1 exp{ ( )}].       main mean E A P P I I      (21) B. Design of () Regularization E 
In traditional LSM, the regulation term used as a constraint on the area and the length of the active contour [START_REF] Chan | Active Contours Without Edges[END_REF] is expressed as

( ) ( ) ( ) ( ) ( ) , Area Length E E E H dx H dx               (22)
where  and  are user controlling parameters. Its derivative with respect to  is expressed as

( ) ( ) ( ). Length Area E E E div                      (23)
In this paper, we introduced a local statistics-based constraint in the classical regularization term. Therefore, the modified regularization term is expressed as follows

( ) ( ) ( ) ( ) ( ) ( ) ( exp{ }) (1 ( )) ,              Regularization E H dx H dx C with C sign I H I H dx           (24)
where  and  are positive parameters. The term () C  is designed in order to attract the evolving curve towards the sectors of the image domain where the gradient is medium or high. Since the LSF  is defined so that it is positive inside the zero level contour and positive outside. A straightforward analysis of () C  leads to the following conclusion,

1) 0 ( ) 0, exp{ } ( ) ( ) 2) 0 exp{ } ( ) 0.                   C if I C sign I if I C      
We can therefore notice that when the curve is evolving, () C  is different from zero and has the sign of I when the curve is on a pixel for which exp{ }  I  . This makes the curve expands or contracts in order to evolve toward higher gradient pixels.

But since ()

main E
 is designed so that the evolving curve cannot stop on high gradient pixels, in order to make the proposed model robust against outliers, the designed regularization term will increase the ability to detect objects delineated by medium edges and weak edges depending on the threshold parameter  . The derivative of the proposed regularization term with respect to  can be written as follows

( )[ ( ) ( ) ( exp{ })]. Regularization E div sign I H I                   (25)
Finally, we obtain the following LSE

12 2 1 2 1 1 2 2 2 ( )( ( , )[1 exp{ ( )}] ( ) ( ) ( exp{ })) ( , ) (1 exp{ ( ( ) ) })
(1 exp{ ( ( ) ) }).

I x P                                     (26)
In order to extend the evolution to all the level set of  , the gradient projection method [START_REF] Rosen | The Gradient Projection Method for Nonlinear Programming, II, Non-linear Constraints[END_REF] allows us to replace ()  by   . To maintain the suitability to parallel programming of our model, we use the local LBM to solve the obtained LSE. By setting  as a signed distance function, that is, 1  , and constraining it to stay like that, Eq. ( 26) can be written as follow

                     (27)
which is similar to Eq. ( 9) with the body force expressed as

The proposed LSE can therefore be solved using the following lattice Boltzmann evolution equation

2 12 ( , 1) ( , ) (1 )[ ( , ) ( , )] ( ( , )[1 exp{ ( )}] ( ) ( exp{ })),                     i i i i i mean f r e t f r t f r t f r t D bc A P P II sign I H I       (29)
without the necessity to explicitly calculate the computational expansive curvature term since it is implicitly handled by LBM. Fig. 3 illustrates the flowchart of the proposed level set based algorithm, where we can clearly distinguish which part is executed on the CPU and which one on the GPU.

The principal implementation steps of the proposed method are as follows:

Input:

Initial zero level contour  (signed distance function),  ,  ,  ,  ,  ,  and  .

Output: The final zero LSF contours  .

Algorithm

IV. EXPERIMENTAL RESULTS AND ANALYSIS

This section demonstrates the performance of the proposed architecture-driven level set image segmentation method (A-DLSO) both in terms of efficiency and effectiveness.

A. Experimental setup

The experiments environment is the parallel computing toolbox of Matlab R2013a installed on a laptop Aspire V3-571G with an Intel Core (TM) i3-2348M (2.3GHz, 3MB L3 cache) processor and possessing a GPU NVIDIA GeForce 710M with 2GB dedicated VRAM. The executive time in all the experiments invoking the graphics card includes the transfer time of data from the CPU to the graphics card and vice versa. In this work, we have empirically set the parameters. In all the experiments, the chosen parameters are those who give the highest F-measure average. For our model, we fixed The optimized Matlab function arrayfun is used to execute the A-DLSO, the CV and the HZ code on the GPU. For example, the body force of the A-DLSO is computed using the following instructions

1-Id = gpuArray(I); 2-Imeand = gpuArray(Imean); 3-Fd = arrayfun(@Body_force,Id,Imeand).

The first and the second instruction transfer respectively 𝐼 and from the CPU to the GPU, while the third instruction computes the body force on the GPU using the kernel function Body_force.m programmed according to Eq. [START_REF] Balla-Arabé | Level Set Region Based Image Segmentation Using Lattice Boltzmann Method[END_REF]. Since all the output arguments are stored in the GPU memory, the function gather is used to transfer them back from GPU to CPU.

To best use the GPU capability, we restructure our code in order to maximally reduce the number of loops by using code vectorization technic. A good explanation of this method can be found in Matlab online support.

In this work, The GPU task and thread assignments are automatically handled by Matlab. Nevertheless, we should notice that the introduced model could be considerably faster by using the CUDA framework with our own task and thread management strategies. Furthermore, since the curve evolution is fully local, the use of CUDA kernels can provide a big speed-up, for example the collision step of the LBM can be a good candidate.

The supervised objective evaluation is undertaken using the F-measure based on precision and recall. It measures the similarity between two images. The higher it is, the better the segmentation result. It is formulated as follows [START_REF] Wang | Relationships Between Diversity of Classification Ensembles and Single-class Performance Measures[END_REF][START_REF] He | Learning from Imbalanced Data[END_REF] In order to better allow the subjective evaluation, we presented the results of all the methods that we used in two forms, by binary images and by only contours. To obtain the binary representation, the interior of the obtained contour is represented by white pixels and the exterior by black pixels. The blue triangle is the initial contour. The dimensions used are 450 948 

, with the exception of Fig. 4 where the dimensions of the image are of 948 450  . The proposed method and all the methods used for comparison were run using the intensity information I of color images, which was obtained by performing a weighted sum of the R, G and B components as recommended by the 601 resolution of the International Commission on Illumination (CIE), 0.2989 0.5870 0.1140 .

      I R G B (31)

B. Experimental results

In this section we compare the introduced framework with the following methods 1-the recent parametric kernel graph cuts (PKGC) based image segmentation method introduced by Ben Salah and Mitiche in [START_REF] Ben Salah | Multiregion Image Segmentation by Parametric Kernel Graph Cuts[END_REF]; 2-the region-based ACM (CV) introduced by Chan and Vese in [START_REF] Chan | Active Contours Without Edges[END_REF]; 3-the two-phase fuzzy c-means (FCM) clustering method; 4-the two-phase k-means clustering algorithm; 5-the local ACM (HZ) introduced by Hagan and Zhao in [START_REF] Hagan | Parallel 3D Image Segmentation of Large Data Set on a GPU Cluster[END_REF]; 6-the ACM for image segmentation in the presence of intensity inhomogeneities introduced by Li et al. in [START_REF] Li | A Level Set Method for Image Segmentation in the Presence of Intensity Inhomogeneities With Application to MRI[END_REF]; 7-the fast ACM (GF) described by Gibou and Fedkiw in [START_REF] Gibou | A Fast Hybrid k-means Level Set Algorithm for Segmentation[END_REF].

Fig. 4 gives evidence of the way the proposed method improves clustering algorithms by allowing them to achieve sub-pixel accurate segmentation with smooth and closed contours. In this work, we used the mean values of the two final classes obtained by the k-means algorithm. Nevertheless, another clustering algorithm, such as the FCM or support vector machine (SVM), can be used. b) and (c), more useful contours are detected by our model. In Fig. 6 (d), we superimposed the k-means segmentation result with the final level set contour of the proposed algorithm represented in red. We can see that the contour is closed, smoother, and the F-measure, comparing the result with human segmentation, gives a score of 89.95 % for our introduced framework and only 77.09 % for the k-means. Thus, the result obtained by our model is far better in term of accuracy. This is an ineluctable proof on how it can greatly improve the quality performance of two-phase clustering algorithms.

Fig. 5 demonstrates the robustness of the proposed method against outliers. In this experiment, we simulated outliers using salt and pepper noise with different density values. The comparison with the k-means algorithm clearly shows the superiority of the proposed framework. The last row shows that even if the k-means fails, the A-DLSO is able to achieve a pretty good segmentation thanks to the use of local statistics, and the constraints on the area and the length introduced during the curve evolution.

From Fig. 6 to Fig. 8, we challenged the proposed algorithm using the PKGC, the CV, the FCM, the k-means, the HZ, the Li and GF methods. The results of the supervised objective evaluation are displayed in Table I and the 3D histogram of Fig. 9. The executive times of the serial implementation on the CPU of all the algorithms used in this section are presented in Table II. Finally, the executive times of the GPU-based implementation of the A-DLSO, CV and HZ are displayed in Table III.

According to Table I, it can be seen that in almost all the experiments the proposed algorithm has the highest average F-measure, i.e., the produced result is the closest to human segmentation used as ground-truth. The subjective evaluation also confirms the superiority of the our model since it extracted more useful, thin and non-discontinuous contours. While the CV method got trapped in local minima, like in Fig. 6 (b), Fig. 6 (c) and Fig. 7 (b), the proposed method gives very good results which confirms its robustness against initialization, and therefore its suitability to automatic systems. Comparing with the A-DLSO, the HZ method presents lower performance, for example in Fig. 6 (a) and Fig. 8 (b), where it gives over-segmentation results. In many cases, the Li method also gives over segmented results, as demonstrated by Fig. 6 (a), Fig. 6 (c) and Fig. 8 (c). The GF method is less effective than the A-DLSO, particularly when it comes to detecting steep corners, or in terms of global segmentation like in Fig. 6 (b), Fig. 6 (c) and Fig. 7 (b). In comparison with the ground-truth, the A-DLSO gives better results than the parametric kernel based method, PKGC. As with the k-means clustering algorithm, the FCM and its derivatives can also be improved by our framework.

When comparing the executive times displayed by Tables II and III, we can see that even when serially implemented on the CPU, the proposed LSM is faster than almost all the methods used in these experiments. Aside from the k-means algorithm, only the HZ method, which uses a very simple speed function, is faster than the introduced method, but the quality of its results is far lower. The A-DLSO is far faster than the CV and the Li based LSM. In almost all cases, it is more than one hundred times faster than the CV method, and more than eighty times faster than the Li method. A straightforward analysis of the executive times of the GPU-based implementation evinces that the curve evolution part of the introduced method has been accelerated more than ten times compared to its serial version on CPU. This is a strong proof of its suitability to massively parallel architecture. The CV method which is non-local, because of the computation of the interior and exterior means of the active curve at each looping, has been accelerated only roughly two fold. Overall, we can conclude that the introduced A-DLSO presents very attractive performance in terms of speed comparing with level set based models.

Fig. 10 proves the robustness of the A-DLSO to initialization. It can be seen that from different initial contours, the proposed methods quietly produces the same finally good result, while the CV method is almost always trapped into a local minimum.

V. CONCLUSIONS AND PERSPECTIVES

This work presents an effective and efficient level set based image segmentation framework suitable for heterogeneous CPU-GPU architecture. The method which is also fast when serially implemented on the CPU, allows considerably improving the quality of the segmentation results obtained using most two-phase clustering algorithms such as SVM, k-means, FCM and their derivatives. Thus, by combining the advantage of LSMs and clustering methods, the A-DLSO is fast, gives smooth contours with sub-pixel accuracy, and can easily handle complex shapes. Its robustness to initialization makes it a good candidate for automatic systems. Intensive experiments demonstrate the good performance of the proposed method. Future works will cover an extension of the A-DLSO to multiphase image segmentation, while conserving its intrinsic locality. The memory complexity aspect will also be considered, since it is one of the limits of the LBM used to solve the proposed LSE. Furthermore, some research on A-DLSO automatic parameterization will be carried out using some machine learning technics. The first row presents the original images with the initial contour, and the corresponding human segmentation as ground truth. The second row presents the results of the proposed algorithm. The third row presents the results of the k-means clustering method. The fourth row presents the results of the FCM clustering method. The fifth row presents the results of the PKGC method. The sixth row presents the results of the CV method. The seventh row presents the results of the HZ method. The eighth row presents the results of the Li method. The last row presents the results of the GF method. For all the methods, the results are presented in two ways, i.e., by contours and binary images.

  Fig. 1. Organization ACMs

Fig. 2 .

 2 Fig. 2. Spatial structure of the D2Q5 LBM lattice

1 m and 2 m 1 (

 121 average,  and  are user-controlling parameters,  and  are positive constants, are respectively the mean intensity values inside and outside the active contour and have to be computed at each iteration as follows

12 (

 12

12 (

 12

Fig. 3 .

 3 Fig. 3. Flowchart representing the process of the proposed algorithm.

P

  Positive True Positive False Negative and P True Positive True Positive False Positive are, respectively, recall and precision,  is usually set to 1 .

  Fig. 4 gives evidence of the way the proposed method improves clustering algorithms by allowing them to achieve sub-pixel accurate segmentation with smooth and closed contours. In this work, we used the mean values of the two final classes obtained by the k-means algorithm. Nevertheless, another clustering algorithm, such as the FCM or support vector machine (SVM), can be used. Fig.6(a) displays the original image with human segmentation as ground truth. As

  method for image segmentation using fuzzy clustering and lattice Boltzmann method," IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics, vol. 43, no. 3, pp. 910-920, 2013. [4] C. Benedek, X. Descombes, and J. Zerubia, "Building Development Monitoring in Multitemporal Remotely Sensed Image Pairs with Stochastic Birth-Death Dynamics," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 34, no. 1, pp. 33-50, 2012. [5] K. Zhang, Q. Liu, H. Song and X. Li, "A Variational Approach to Simultaneous Image Segmentation and Bias Correction," IEEE Transactions on Cybernetics, vol. 45, no. 8, pp. 1426-1437, 2015. [6] S. Balla-Arabé, and X-B. Gao, "Image multi-thresholding by combining the lattice Boltzmann model and a localized level set algorithm," Neurocomputing, vol. 93, pp. 106-114, 2012. . (a) Natural image with human segmentation as ground truth. (b) Segmentation result of the proposed algorithm, binary and contour representations. (c) Segmentation result of the k-means algorithm, binary and contour representations. (d) Superimposition of the k-means result with the contours (in red) extracted by the method introduced in this paper.Table II. EXECUTIVE TIMES OF THE SERIAL IMPLEMENTATION ON THE CPU (IN SECOND)

Fig. 5 .

 5 Fig. 5. Segmentation of an image corrupted by salt and pepper noise with different density values d . The first column displays the corrupted images. The second column presents the results of the proposed method. The third column presents the results of the k-means algorithm.

Fig. 7 .

 7 Segmentation of real-world images.

Fig. 9 .Fig. 10 .

 910 Fig. 9. 3D histogram showing the supervised objective evaluation. The vertical axis represents the F-measure, while the horizontal axis represents the evaluated methods.

  a spatial variable,  and  are user controlling

	positive constants, 1  and 2  are important parameters which
	can lead, if wrongly chosen, to an over or under segmentation
	result. In this work, we used 12 1  	. The values of 1 P and
	2 P are obtained as follows. The interior and the exterior
	intensity average of the initial active curve are used to initialize
	a two-phase k-means algorithm, 1 P and 2 P are the mean
	values of the two final classes obtained. This part constitutes
	the serial part of our algorithm, and it is fast since
	computational efficiency is one of the main advantages of the
	k-means algorithm.
	The derivative of	()

main E

Table I .

 I STATISTICAL RESULTS OF THE OBJECTIVE EVALUATION USING F-MEASURE (%).

	Methods	Average	Standard deviation
	A-DLSO	90.7841	0.9028
	CV	78.0608	2.7437
	HZ	76.1850	1.5957
	k-means	79.0566	1.0073
	FCM	79.9786	0.8641
	PKGC	83.1696	0.8743
	Li	81.7851	1.9599
	GF	78.1173	1.6935

Table III .

 III EXECUTIVE TIMES OF THE GPU-BASED IMPLEMENTATION (IN SECOND).

		6468	0.4394	0.4366	0.4706	0.4631	0.6523	0.5062	0.4733
	Curve evolution	0.2146	0.2107	0.2106	0.2092	0.2108	0.2160	0.2138	0.2122
	Total	0.8614	0.6501	0.6472	0.6798	0.6739	0.8683	0.7200	0.6855
	CV	85.007	84.520	85.694	84.971	84.066	85.443	137.55	84.741
	HZ	0.4497	0.4508	0.4595	0.4598	0.4603	0.4575	0.4562	04541
	K-means	0.6468	0.4394	0.4366	0.4706	0.4631	0.6523	0.5062	0.4733
	FCM	10.214	12.614	8.4572	10.355	9.3448	12.625	8.3827	12.758
	PKGC	2.6772	2.9213	2.4635	2.4909	2.4376	2.4485	2.7797	2.9501
	Li	70.2741 66.8958 74.2536 50.4556 55.4759 48.6515 79.9054 56.3074
	GF	0.9734	0.7701	0.7735	0.7898	0.7802	0.9750	0.8117	0.7993


Souleymane Balla-Arabé

Xinbo Gao

Dominique Ginhac

Fan Yang

Shape-Constrained Level Set Segmentation for Hybrid CPU-GPU Computers

Keywords: image segmentation, Massively parallel architectures, Partial differential equations, Shape prior, level set method image segmentation, Massively parallel architectures, Partial differential equations, Shape prior, level set method

Due to its intrinsic advantages such as the ability to handle complex shapes, the level set method (LSM) has been widely applied to image segmentation. Nevertheless, the LSM is computationally expensive. In order to improve the performance of the traditional LSM both in terms of efficiency and effectiveness, we propose a novel algorithm based on the lattice Boltzmann method (LBM). Using local region statistics and prior shape, we design an effective and local speed function for the LSM, from which we deduce a shape prior based body force for LBM solver. An NVIDIA graphics processing units (GPU) is used to accelerate the method. Ourintroduced algorithm has several advantages. First, it is accurate even if they are some geometric transformations (rotation angle, scaling factor and translation vector) betweenthe object to be segmented and the prior shape. Second, it is local and therefore suitable for massively parallel architectures. Third, the use of local region information allows it to deal with intensity inhomogeneities. Fourth, including shape prior allows the method to handle occlusion and noise. Fourth, the model is fast. Finally the algorithm can be used without shape prior by means of minor modification. Intensive Experiments demonstrate, objectively and subjectively, the performance of the introduced framework.

Introduction

In computer vision and pattern recognition, image segmentation [START_REF] Xu | Multiple clustered instance learning for histopathology cancer image classification, segmentation and clustering[END_REF]- [START_REF] Kim | A hierarchical image clustering cosegmentation framework[END_REF], [START_REF] Ben Salah | Mitiche: Active Curve Recovery of Region Boundary Patterns[END_REF]- [START_REF] Yang | Layered Object Models for Image Segmentation[END_REF] is a major process by which a given image is partitioned into number of meaningful and homogeneous regions, such that the union of any two neighborhood regions yields a heterogeneous segment. The task is non-trivial and more challenging in presence of noise, intensity inhomogeneities or occlusion.

Recently, optimization methods have been widely used as highly effective image segmentation strategies. They basically achieve the segmentation by minimizing a given energy function designed from the image information. Optimization methods can be roughly classified into two important classes: spatially discrete and spatially continuous representations. In spatially discrete approaches, the pixels of the image are usually considered as the nodes of a graph, and the aim of segmentation is to find cuts of this graph which have a minimal cost [START_REF] Boycov | Fast Approximate Energy Minimization Via Graph Cuts[END_REF]- [START_REF] Malik | Contour and Texture Analysis for Image Segmentation[END_REF].

As shown by Fig. 1, Active contour models (ACMs) belong to the spatially continuous approaches: the segmentation of the image plane is considered as a problem of infinite-dimensional optimization.In this class, variational methods are used in the image segmentation. The principle is to evolve the active curve in the direction of the negative energy gradient by means of an appropriate partial differential equation. The implicit active curve method or level set method (LSM) [START_REF] Ben Salah | Mitiche: Active curve recovery of region boundary patterns[END_REF]- [START_REF] Prisacariu | Nonlinear shape manifolds as shape priors in level set segmentation and tracking[END_REF] designates the class of ACM which uses the Eulerian framework, that is, the geometric representation of the evolving curve, instead of the parametric one, that is, the Lagrangian framework [START_REF] Zhu | Region Competition: Unifying Snakes, Region Growing, and Bayes/MDL for Multiband Image Segmentation[END_REF]- [START_REF] Nakhmani | Self-crossing Detection and Location for Parametric active contours[END_REF]. * Corresponding author: Tel.: +862988201838; fax: +862988201620. E-mail address: xbgao@mail.xidian.edu.cn

The LSM has several intrinsic advantages, such as the ability to easily handle complex shapes, and topological changes comparable to parametric active contours. It is, furthermore, straightforward to pass from two-dimensional (2D) to 3D space, and to add some constraints on the smoothness of the boundaries via some regularization terms.

The original idea of the LSM stems from the Hamilton Jacobiapproach, i.e., a time-dependent equation for a moving surface [START_REF] Osher | Fronts propagating and curvature dependent speed: algorithms based on Hamilton-Jacobi formulation[END_REF]- [START_REF] Balla-Arabé | Image multi-thresholding by combining the lattice Boltzmann model and a localized level set algorithm[END_REF]. In two-dimensional (2D) space, the LSM aims to evolve a given curve toward its interior or exterior normal until defining the boundary of the object of interest. The curve evolution is driven by the level set equation (LSE) which is a partial differential equation and, in its simple form can be expressed as

   V t   , (1)
where  is the level set function (LSF),   is the gradient of  and V is the speed function which drives the active contour towards the region boundaries. Two approaches are usually used to stop the evolving curve on the boundary of the desired object.The first one uses an edge indicator depending on the gradient of the image like in classical snakes and active contour models [START_REF] Kass | Snakes: Active contour models[END_REF]- [START_REF] Malladi | A topology independent shape modeling scheme[END_REF]. The second one uses some regional attributes to stop the evolving curve on the actual boundary [START_REF] Caselles | On geodesic active contours[END_REF]- [START_REF] Paragios | Geodesic active contours for supervised texture segmentation[END_REF]. The latter is more robust against noise and can detect objects with weak edges. One of the most interesting approaches was proposed in [START_REF] Chan | Active contours without edges[END_REF] where Chan and Vese introduced a level set formulation to minimize the Mumford and Shah functional [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF]. They converted the problem into a mean curvature flow problem just like the active contours. But the results were better than the classical active contours because the stopping term did not depend on the gradient of the image which reduces the dependence on strong edges and improves the robustness against noise. However the method cannot deal with intensity inhomogeneities, and is not suitable for parallel programming because at each iteration the average intensities inside and outside the contour should be computed, which increases drastically the CPU time by increasing communications between processors. Furthermore, in order to solve the LSE, most of the available methods suggest the use of some computational expensive finite difference, finite element or finite volume approximations and an explicit computation of the curvature [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF].

In [START_REF] Balla-Arabé | Level set region based image segmentation using lattice Boltzmann method[END_REF], the authors used the lattice Boltzmann method (LBM) as an alternative approach for solving the LSE for region-based image segmentation. The method can better handle the problem of time consuming because the curvature is implicitly computed and the algorithm is simple. Nevertheless, in this method the proposed stop function is not local, thus cannot deal with intensity inhomogeneities and at the same time is not suitable for GPU-based computations.

In the present paper, we propose a new supervised method which improves and overcomes some limitations of the one proposed in [START_REF] Balla-Arabé | Level set region based image segmentation using lattice Boltzmann method[END_REF]. Using local region properties of the image, we design a new localized stop function which includes shape prior information. The prior shape is used to roughly guide the evolving curve toward boundaries of the object to detect, and the local region statistics allow the detection of fine details. Moreover, since sometimes they are some small geometric transformations (rotation angle, scaling factor and translation vector) between the object to be segmented and the prior shape, we introduced a straightforward adjustment method which effectively increases the accuracy of the method. The proposed algorithm can deal with intensity inhomogeneities, occlusion and noise, and is highly parallelizable due to its local property. Furthermore, it can be easily used without shape prior by means of minor modification. The rest of this paper is organized as follows: Section II briefly introduces the LBM. In section III,we present the formulation of the proposed method. Section IV validates the performance of the proposed image segmentation method through experimental results. Section V concludes the paper.

Background

This section gives an overview of the LBM method which is used in this paper to solve the LSE.The LBM models Boltzmann particle dynamics on a 2D or 3D lattice. At first, it was designed to solve macroscopic fluid dynamics problems. The LBM is second-order accurate both in time and space, and in the limit of zero time step and lattice spacing, it yields the Navier-Stokes equations for an incompressible fluid [START_REF] He | Lattice Boltzmann model for incompressible Navier-Stokes equation[END_REF].

In this paper, we use the D2Q9 (2D with 8 links with its neighbors and one link for the cell itself) LBM lattice structure. Fig. 2 shows a typical D2Q9 model.

The evolution equation of LBM can be written as

( ) . ,              eq ii ii ii D f r t f r t f r t f r t F e c e b  (2) 2 ( ) ( 1 , 1 , [ , , ] ) ( )
where

 i e
is the velocity vector of a given link i , ( , )

and eq i f is the local equilibrium particle distribution usually given by the Bhatnager, Gross, Krook

(BGK) model 22 ( ) ( ) , ( . ()) () .           ii eq i i i i i f u A B u C u e Du e  , (4) 2 2
, fe 00 , fe ,

  i i f  (5) 1     i i i uf e  . (6)
For diffusion problems, the equilibrium function can be simplified as [START_REF] Zhao | Lattice Boltzmann based PDE solver on the GPU[END_REF] ( , ) 

 eq ii f u A  . (7)
In the case of D2Q9 model, 49

 i A
for the zero link, 19

 i A
for the axial links, 1 36

 i A
for the diagonal links and the relaxation time  is then determined using the following equation 2(2 1) 9

  , ( 8
)
where  is the diffusion coefficient. By performing the Chapman-Enskog analysis the following diffusion equation can be recovered from the LBM evolution equation (Eq. ( 2)),

()      div F t  , ( 9
)
where div is the divergence operator. Substituting  by the signed distance function  in Eq. ( 9), the LSE can be recovered. The body force F represents the link with the image data in the LBM solver.

Another approach for the LBM based level set segmentation is to add a medium between the nodes of the lattice [START_REF] Chen | Cellular automata based level set method for image segmentation[END_REF] in order to design a stop function for the evolving active contour. The particles can pass through the medium with a possibility of () i gr  , and will be punched back where they were with a possibility of 1 ( ) i gr   . They therefore modified the LBM evolution equation as

1 , 1 ( )[ , [ ( ) ( ) ( ) ( ) ( , , ] ] (1 ( )) , )                   eq i i i i i ii ii f r t g r f r t f r t f r e e t g r f r t   , ( 10
)
where σ is the convection coefficient.

TheShape-Constrained Level Set Algorithm

In this section, we first design and analyze the proposed localized shape prior algorithm, and then we give the main implementation steps.

Design and analysis of the proposed algorithm

For the level set equation defined in (Eq. ( 1)), we design the following shape prior based speed ). ( )

function 12 V ( ) ( ) () ( ) ( ( ) ) 2        shape e r e r r f r r   , (11) with ( ) ( ) ( ( ( ))
( ( )))        shape ref f r I r sign r sign r    , (12) where
               r s s s rs K I H d er K H ds   , (13) 2 (
). ( ).(1

( () )) ( ).( 1
( ))                  r s s s KI rs H KH s d er d   , ( 14
)
where H is the Heaviside function and

1 ( , 0, )          rs rs otherwise K  ( 15
)
with  the radius parameter and 𝑠 a spatial variable.K is a characteristic function used to mask local regions, it will be 1when the point𝑠 is within a ball of radius  centered at 𝑟 , and 0 otherwise.

The term ()

 r  is the non-linear curvature term used to smooth the contour. It can be written as

()       r    , ( 16
)
where  is a controlling parameter.

The proposed shape prior based level set equation is therefore

12 ( ) ( ( ( )) ( ( ))) ) 2 ( ( ) ( ) ).                 ref I r sign r sign r e r e r t         (17)
In order to make the method suitable for parallel programming, we use the local LBM to solve the proposed LSE. Since we consider  as a signed distance function, i.e., 1   , Eq. ( 17) can therefore be expressed as

)

) 2                ref I r sign r sig r e r e r t n       (18)
which is similar to Eq. ( 9) with the body force defined as

The proposed level set equation can thus be solved using the lattice Boltzmann evolution equation Eq.

(2) without the necessity of explicitly computes the computational expensive curvature term since it is implicitly handled.Furthermore, since the above defined body force is also local, it will not affect the local nature of the LBM. This will result in an adequate algorithm to massively parallel devices such as the GPU because the time consuming communications between the processors will be greatly reduced, as well as the executive time.

From Eqs. ( 11) and ( 12) we can do the followings remarks: 1-Considering an object with or without edges, ()  Vr reduces to zero when the active pixel intensity verify

( ) ( ( ( )) ( ( ))) 2 ( ) ( )         ref e r e r I r sign r sign r    , (20)
2-In the case of partial occlusion, the second term in the right hand of Eq. ( 20) makes the contour continue to evolve instead to stop on the false boundaries caused by the occlusion.

3-If the image template is not available, we just have to fix  ref  or 0   in the algorithm and to perform the segmentation without prior shape. According to the above statements, using ()  Vr as a stop function will theoretically produce interesting results. This will be confirmed by experimental results.

For the shape adjustment, we use the following method which can handle the cases of geometric transformations (translation, rotation and scaling). Let consider Fig. 3 where  is the rotation angle.

In order to determine  , the translation parameter  H and the scaling parameter S , we have the following steps. We first choose two remarkable points on the shape prior (A and B), which can be easily recover on the object to be segmented (C and D). In this paper, the selected points are vertex of some well-defined angles.O

( )( ) ( )( ).              B A D C B A D C AB CD AB CD x x x x y y y y  ( 21
)
We can then compute the rotation angle

 since 22 ( )( ) ( )( ) 1 cos( ) . ( ) ( )          B A D C B A D C B A B A x x x x y y y y S x x y y  , ( 22
)
where S is the scaling parameter defined as

2 2 2 2 [( ) ( ) ] [( ) ( ) ].           D C D C B A B A S CD AB x x y y x x y y ( 23
)
In the shape adjustment process, after the scaling and the rotation steps, the translation parameter  H can be written as Note that before the scaling and the rotation steps in the shapealignment process, we have 0   HH .

Implementation

The steps for the implementation of the proposed algorithm are outlined as follows:

a) Initialize  as a signed distance function;

b) Do the shape alignment, which consists to first put the template at the same scale with the object to segment, the scaling parameter being S . The second step is the rotation according to  , and the last step consists to translate the geometric template according to

 H ; c) Compute 1
e and 2 e according to Eqs. ( 13) and ( 14); d) Calculate ()  Fr according to Eq. ( 19); e) Resolve the proposed level set equation with the LBM evolution equation Eq. ( 2);

f) Accumulate the ( , ) i f r t 
values at each grid point using Eq. ( 5) to get the updated values of  at each point; g) Find the contour which is the zero level of  ; h) If the algorithm has not converged, i.e., , go back to step c).

Fig. 4 illustrates the flowchart of the introduced supervised level set based algorithm. It can be seen that the proposed method is simple to implement since the LBM eliminates the use of complex approximations to solve the LSE and the need to explicitly compute the curvature.

Experimental Results Analysis

In this section we firstly demonstrate the performance of the proposed image segmentation method in term of efficiency and efficacy in presence of intensity inhomogeneities and object with weak or without edges. Four level set based segmentation methods: the Chan and Vese's method (CV) [START_REF] Chan | Active contours without edges[END_REF], the Li's method [START_REF] Li | Distance regularized level set evolution and its application to image segmentation[END_REF], the Chen's method [START_REF] Chen | Cellular automata based level set method for image segmentation[END_REF] and the Aaron and Zhao's method (AZ) [START_REF] Hagan | Parallel 3D image segmentation of large data set on a GPU cluster[END_REF] are selected as benchmarks. Secondly, we demonstrate the effectiveness of the proposed method in the case of partial occlusion and noise corruption. Finally, we show that the our method is effective even without prior shape. The first, the second and the third instructions transfer respectively I ,  and ref  from the CPU to the GPU, while the fourth instruction computes the body force on the GPU using the kernel function Body_force.m programmed according to Eq. [START_REF] He | Lattice Boltzmann model for incompressible Navier-Stokes equation[END_REF]. Since all the output arguments are stored in the GPU memory, the function gather is used to transfer them back from GPU to CPU.In all the experiments the radius parameter  is set to 4 , the weighting parameter  to 3 and the controlling parameter  also to 3 .

The prior shape of the corpus callosum used in the experiments is the mean shape of human's segmentation of a training set. In the case of natural images, the more suitable prior shapes have been determined in a dataset of possible prior shapes. Then, they have been rescaled, translated and/or rotated in order to make the experimental cases more difficult.For more information about prior shapes detection or construction, one can refer to [START_REF] Charpiat | Shape Statistics for Image Segmentation with prior[END_REF]- [START_REF] Paragios | Shape Priors for Level Set Representation[END_REF].

To objectively evaluate the performance of the proposed method, we use the Levine and Nazif (LN) inter-region and intra-region contrast criterion [START_REF] Levine | Dynamic measurement of computer generated image segmentations[END_REF]. The lower is the intra-region contrast, the better is segmentation result, and the higher is the inter-region contrast, the better is the segmentation result. Fig. 5 demonstrates the proposed method on an MRI image of head and some real world images.

Row I displays the prior shapes, row II the initial contours, row III the results of the Chen's method, row IV the results of the Li's method, row V the results of AZ method, row VI the results of the CV's method and row VII the results of the proposed method. The executive times are displayed in Table 1 and the objective evaluation results in Tables 2 and3. In nearly all the presented cases the proposed method has the lowest intra-region contrast and the highest inter-region contrast. Thus, it objectively gives the best results, which can be confirmed subjectively by observing the segmented images. The Chen's and the Li's methods fail in most of the cases because they are edge based and the objects are with weak or without edges. The CV's method is based on non-local means values and is therefore not effective when they are intensity inhomogeneities like in most of the presented cases. In column (a), only the proposed method has well segmented the corpus callosum in the head magnetic resonance image. Furthermore it is the fastest one; in this case the executive time is of 0.3465 s versus 41.89 s for the CV's method. In column (b) and (c), the results of the proposed method are also pretty good, and demonstrate its effectiveness in presence of intensity inhomogeneities, weak edges and for texture image segmentation. Fig. 6 shows the performance of the proposed method in the case of partial occlusion. In all the cases, the segmentation results are promising. In Fig. 8 we demonstrate that the method we propose can be used even if the template is not available. As already specified, in this case we just need to fix ref   . The results obtained are subjectively promising.

Nevertheless, we should note that the effectiveness of our introduced model is limited in case of occlusion when the prior shape is incorrect. Furthermore, the model can lead to an inappropriate result if the recovery of the points A and B on the object to be segmented is inaccurate. This can limit the use of the method in automatic systems since it requires human interaction to be fully effective when used with prior shape.

Conclusion

In this paper, we have presented a novel shape prior level set image segmentation algorithm based on the lattice Boltzmann model. It is effective in presence of noise, intensity inhomogeneities and partial occlusion, and can be used if the prior shape is not available by means of minor modification.

Furthermore, it is local and thus highly parallelizable, which allows it to be a good candidate for parallel implementation. Experimental results using a graphics processing units on medical and real-world images have demonstrated the good performance of the proposed method in terms of efficiency, effectiveness and robustness.

Fig. 1 .

 1 Fig. 1. Organization of optimization methods.

F

  the body force, D the grid dimension, b the link at each grid point and c the length of each link which is set to 1 in this paper. The parameter  represents the relaxation time determining the kinematic viscosity  of the fluid

Fig. 2 .

 2 Fig. 2. Spatial structure of the D2Q9 LBM lattice.

r 1 e and 2 e 1 (

 121 is a spatial variable,  a controlling parameter,  a weighting positive parameter,  the level set function which is a signed distance function, ref a signed function obtain from the aligned prior shape, and I the pixel intensity. In this paper,  is positive inside the contour and negative outside, just as ref is positive inside the aligned prior shape and negative outside. The terms are local mean values of a given pixel respectively inside and outside the evolving curve. They are defined as

12

 12

1 and O 2

 2 are respectively the centers of the [AB] and[CD]. The rotation angle is therefore the angle between the straight lines (AB) and (CD). By considering the scalar product of the vectors

Fig. 3 .

 3 Fig. 3. Shape adjustment parameters.

Fig. 4 .

 4 Fig. 4.Flowchart representing the process of the introduced framework.

  The experimental environment is Matlab R2012b installed on a PC AMD Athlon (tm) 5200 processor with a clock speed of 2.31 GHz and 2 GB of RAM, and possessing the NVIDIA GPU GT 430. The optimized Matlab function arrayfun is used to execute the introduced level set algorithm on the GPU.For example, the body force is computed using the following instructions 1-Id = gpuArray(I); 2-ϕd = gpuArray(ϕ); 3-ϕ ref d = gpuArray(ϕ ref ); 4-Fd = arrayfun(@Body_force,Id, ϕd, ϕ ref d).

Fig. 7 Fig. 5 .

 75 Fig. 7 shows the performance of the proposed method in presence of noise. In row II the image is corrupted by the speckle noise with the variance 0.04  v ; in row III the image is corrupted by the Gaussian noise with the mean 0.3  m and the variance 0.3  v ; in row IV the image is corrupted by salt and pepper noise with the density 0.3  d ; finally in row V the image is corrupted by the poisson

Fig. 6 .

 6 Fig. 6. Experimental results of the proposed method in presence of occlusion. Row I: prior shapes, row II: occluded images, and row III: segmentation results of the proposed method.

Future

  work will be the development of an effective method for selecting the points A and B without human interaction. Furthermore, we will work on the extension of the proposed method to volume image segmentation and real-time objects (buildings, cars etc…) extraction in very high resolution satellite images.

Fig. 7 .

 7 Fig. 7. Experimental results of the proposed method in presence of different types of noises. Row I: prior shapes; row II: image corrupted by the speckle noise; row III: image corrupted by gaussian white noise; row IV: image corrupted by the salt and pepper noise; row V: image corrupted by the poisson noise. Column I: initial contours, and column II: segmentation results of the proposed method.

Fig. 8 .

 8 Fig. 8. Segmentation using the proposed method without prior shape. Column I: initial contours, and column II: segmentation results.

Table I .

 I COMPARISON OFEXECUTIVE TIMES (S)

	Methods	Our method	CV	AZ	Chen	Li
	Image (a)	0.3465	41.893	2.767	3.475	53.759
	Image (b)	0.9872	101.76	3.871	5.103	49.784
	Image (c)	0.6154	92.872	2.148	2.327	51.816
		Table II.	LN'S INTER-REGION CONTRAST.
	Methods	Our method	CV	AZ	Chen	Li
	Image (a)	0.2509	0.1711	0.0989	0.532	0.564
	Image (b)	0.2867	0.1254	0.0712	0.788	0.614
	Image (c)	0.2782	0.1103	0.0989	0.477	0.483

Table III .

 III LN'S INTRA-REGION CONTRAST.

	Methods	Our method	CV	AZ	Chen	Li
	Image (a)	1.96e-08	3.4e-07	0.0262	6.4e-03	9.3e-03
	Image (b)	1.72e-07	8.4e-07	0.0435	0.0643	0.0242
	Image (c)	3.44e-07	6.5e-07	0.0911	0.0891	0.0176


Frédéric Meynadier

Pacôme Delva

Christine Guerlin

Christophe Le Poncin-Lafitte

Philippe Laurent

P Wolf

Christophe Le Poncin

ACES MWL DATA ANALYSIS PREPARATION STATUS

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The ACES-PHARAO mission is an ESA space mission aiming at realizing a time scale of high stability and accuracy on board the International Space Station (ISS), and perform time and frequency transfer between the payload and ground stations during the flight. A description of its setup and main objectives may be found in reference [START_REF] Salomon | Atomic Clock Ensemble in Space:. An Update[END_REF] . The SYRTE laboratory has been heavily involved in the design and development of the PHARAO cold atom clock (Laurent et al., 2014) 2 as well as the data analysis software that will generate the mission's scientific products. This paper focusses on the current status of the data analysis.

The ACES experiment setup

In order to compare their local timescales to the ACES timescale, several metrology institutes will host a ACES Ground Terminal (GT), which will enable them to perform a two-way time transfer through microwave signals. ISS visibility will only last 300 to 500 seconds for a given station, with common-view configurations being possible if two GTs are close enough. Collected data on ground and in space is centralized at CADMOS (see fig. 1).

The MWL core data consists in 3 series of one-way time comparisons for each pass : 13.475 GHz uplink, 14.7 GHz Ku-band downlink, 2.24 GHz downlink. Combining the first two signals allows to perform a two-way time transfer, while combining the last two signals allows to determine the delay due to the ionosphere traversal [START_REF] Delva | Time and frequency transfer with a microwave link in the ACES/PHARAO mission[END_REF] . Time transfer data is determined either by determining the phase of the incoming carrier (dubbed "carrier" data), or by matching an encoded pattern at 100Mchips rate (dubbed "code" data). Code data is unambiguous, but coarse (20 ps resolution). Carrier data can theorically achieve 1ps resolution, but suffers from phase amibiguity (any measurement will only be able to give a modulo-2π phase). Combining both measurements correctly is key to optimal performance of the MWL.

Auxiliary data needed for data processing will also be contributed to the archive : instrumental calibrations, ISS orbitography, ground station positions and atmospheric parameters, etc.

Two-way microwave link

The basic principle of such a link has been presented in Delva et al. (2012) [START_REF] Delva | Time and frequency transfer with a microwave link in the ACES/PHARAO mission[END_REF] . In summary : sending a timestamp from clock A to clock B while sending a timestamp the other way, from clock B to clock A, allows (at the first order) to cancel the signal's propagation time between the two clocks, thus allowing the calculation of the desynchronisation without precise knowledge of this time of flight. However, we remain sensitive to the variation of this value between the uplink and downlink measurement : ISS orbitography and ground station positionning allow us to model this effect to the required level.

As ISS positionnal uncertainty has been found to be a major source of error, we minimize its influence by interpolating signals so that uplink arrival time matches downlink departure time, thus making the "same" error on both trips, which cancels out well (Λ configuration 3 ).

Software development

SYRTE has developped a data processing and analysis software, designed to extract scientific products from (nearly) raw ACES data. To test this software and experiment we also developped a simulation software that generates raw data in the format specified by ACES Ground Segment ICD, aiming at including all relevant effects in the calculation.

Care has been taken to separate as much as possible the development of both softwares : different developpers, different languages (Python + Numpy for the processing software, Matlab for the simulation), different approaches are used in order to avoid common interpretation mistakes.

Current status

Figure 2 is the output of the tool we have designed for validating our software during its coding : as much as possible, each effect is isolated during the process and compared to its theoretical value (i.e. the value we can calculate from simulation's input parameters). The 20 ps spread on code data is expected : it is a direct consequence of the 10 ns resolution of the raw data. Carrier data is affected by an offset of a few 10 ps, which we are currently investigating.

Tropospheric delays (2nd line of fig 2)

Tropospheric delay is a (mainly) non-dispersive effect that can reach several hundred nanosecond on each signal. Two-way measurement cancel it to a large extent, but we want to evaluate it as well as possible at higher orders. Some dispersive, second order effects have been studied by Hobiger, Piester and Baron (2013) [START_REF] Hobiger | A correction model of dispersive troposphere delays for the aces microwave link[END_REF] and will be implemented in the future.

Ionospheric delays (3rd, 4th and 5th lines of fig 2) These are due to the dispersive nature of the ionosphere, and therefore affect differently frequences f 1 and f 2 . We separated those in 2 terms, one that is proportional to 1/f 2 , the other that is proportional to 1/f 3 . Slight differences may be seen in the latter case : these have been identified as results of the difference between the model used for Earth's magnetic field in the simulation and in the processing software, and remain largely below detection threshold. STEC (Slanted Total Electron Content) is an intermediate result of this calculation and will be issued as a by-product.

Geometrical time of flight and range (6th and 7th lines of fig 2) Geometrical time of flight is the "classical, atmosphere-free" coordinate time interval between emission and reception of a given signal : As both stations move with respect to a geocentric non-rotating reference frame (GCRF), this calculation needs iterations or Taylor series developments to converge.

Range is the instantaneous distance between the two stations at a given coordinate time. This is calculated from the ISS orbitography files and the ground station coordinates, after conversion to GCRF.

Desynchronisation and associated time deviation (8th and 9th lines of fig 2) This is the end result of the processing : residual spreads are consistent with those of the pseudo ranges (first line), which is expected. The offset for carrier-phase data is the mean between the offsets on pseudo-ranges, it is a direct consequence of the offsets on pseudo-ranges which should disappear once those are removed.

Apart from this offset, no visible effects are noticed throughout the pass, showing that all other effects were correctly removed.

Conclusion

Our analysis software now implements the core functionnalities that will be needed to perform ACES data analysis. Although some important points still need to be elucidated, we are confident that this software will be operationnal on schedule, with an ACES launch foreseen in early 2017.

Figure 1 -

 1 Figure 1 -Data processing flowchart for ACES-PHARAO data. L0 to L4 data levels correspond to various levels of processing from raw data to scientifically usable measurements. The red ellipse indicates the situation of LNE-SYRTE data processing center, at the interface between partially processed data in the archive, and users.

Figure 2 -

 2 Figure2-Typical output from a comparison between the input data (i.e. theoretical values calculated by the simulation), and intermediate/final results calculated by the processing software. In this example, the simulated data includes a keplerian orbit for the ISS, earth rotation, atmospheric delays with variable parameters (troposphere + ionosphere), initial desynchronisation between the clock equal to 0.1 ms, frequency drifts due to velocity and gravitational potential for ground and space. Each point is the difference between the expected data at this date (known from the simulation's scenario), and the data actually calculated by the processing software for the same coordinate time. Ideally we should therefore get values close to zero + noise. Blue points refer to uplink (f1) signals, green points refer to downlink (f2) signals, and red points refer to desynchronisation (which combines both signals). Left column contains result for carrier-phase data, and right column for code-phase data. Sampling period is 80 ms. See text for a description of the lines.

Acknowledgements

FM acknowledges support from Centre National d' Études Spatiales for his participation to the Rencontres de Moriond 2015.