Ramzi Mahmoudi
email: mahmoudr@esiee.fr

Mohamed Akil
email: akilm@esiee.fr

Enhanced Computation Method of Topological Smoothing On Shared Memory Parallel Machines

To prepare images for better segmentation, we need preprocessing applications, such as smoothing, to reduce noise. In this paper, we present an enhanced computation method for smoothing 2D object in binary case. Unlike existing approaches, proposed method provides a parallel computation and better memory management, while preserving the topology (number of connected components) of the original image by using homotopic transformations defined in the framework of digital topology. We introduce an adapted parallelization strategy called split, distribute and merge (SDM) strategy which allows efficient parallelization of a large class of topological operators. To achieve a good speedup and better memory allocation, we cared about task scheduling and managing. Distributed work during smoothing process is done by a variable number of threads. Tests on 2D grayscale image (512*512), using shared memory parallel machine (SMPM) with 8 CPU cores (2× Xeon E5405 running at frequency of 2 GHz), showed an enhancement of 5.2 with cache success rate of 70%.

Introduction

Smoothing filter is the method of choice for image preprocessing and pattern recognition. For example, the analysis or recognition of a shape is often perturbed by noise, thus the smoothing of object boundaries is a necessary preprocessing step. Also, when warping binary digital images, we obtain a crenellated result that must be smoothed for better visualization. The smoothing procedure can also be used to extract some shape characteristics: by making the difference between the original and the smoothed object, salient or carved parts can be detected and measured.

Smoothing shape has been extensively studied and many approaches have been proposed. The most popular one is the linear filtering by Laplacien smoothing for 2D-vector [START_REF] Taubin | Curve and surface smoothing without shrinkage[END_REF] and 3D mesh [START_REF] Liu | A novel volume constrained smoothing method for meshes[END_REF]. Other approach by morphological filtering can be applied directly to the shape [START_REF] Asano | Active contour model based on mathematical morphology[END_REF] or to curvature plot of the object's contour [START_REF] Leymarie | Curvature morphology[END_REF]. Unfortunately none of these operators preserve the topology (number of connected components) of the original image. In 2004, our team introduced a new method for smoothing 2D and 3D objects in binary images while preserving topology [START_REF] Couprie | Topology preserving alternating sequential filter for smoothing 2D and 3D objects[END_REF]. Objects are defined as sets of grid points, and topology preservation is ensured by the exclusive use of homotopic transformations defined in the framework of digital topology [START_REF] Kong | Digital topology: introduction and survey[END_REF]. Smoothness is obtained by the use of morphological openings and closings by metric discs or balls of increasing radius, in the manner of alternating sequential filters [START_REF] Sternberg | Grayscale Morphology[END_REF]. The authors' efforts have brought about two major issues such as preserving the topology and the multitude of objects in the scene to smooth out without worrying about memory management, latency or cadency of their filter. This paper describes an enhanced computation method of topological smoothing filter that assure better performance. We present also a new parallelization strategy, called Split Distribute and Merge (SD&M). Our strategy is designed specifically for topological operator's parallelization on shared memory architectures. The new strategy is based upon the exclusive combination of two patterns: divide and conquer and event-based coordination. This paper is organized as follows: in section 2, some basic notions of topological operators are summarized; the original smoothing filter is introduced. In section 3, parallelization strategy, that has been adopted, is introduced. We define the class of operators that our strategy may cover. Motivations for using SMP machines are also cited. Threads coordination and tasks scheduling are discussed. In section 4, the new parallel smoothing method is introduced and evaluations of acceleration, efficiency and success rate of cache memory access are also presented and discussed. Finally, we conclude with summary and future work in section 5.

Theoretical background

In this section, we recall some basic notions of digital topology [START_REF] Kong | Digital topology: introduction and survey[END_REF] and mathematical morphology for binary images [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. We define also the homotopic alternating sequential filters [START_REF] Couprie | Topology preserving alternating sequential filter for smoothing 2D and 3D objects[END_REF]. For the sake of simplicity, we restrict ourselves to the minimal set of notions that will be useful for our purpose. We start by introducing morphological operators based on structuring elements which are balls in the sense of Euclidean distance, in order to obtain the desired smoothing effect.

We denote by ℤ the set of relative integers, and by Ε the discrete plane 2 ℤ . A point x ∈ Ε is defined by

{ } 4 1 1 2 2 ( ) ; 1 x y E y x y x Γ = ∈ -+ -≤ , { } 8 1 1 2 2 ( ) ; max , 1 x y E y x y x Γ = ∈ - -≤ .
For general case, we define

{ } * ( ) ( ) \ n n x x x Γ = Γ with { } 4,8 n ∈ . Thus y is said n-adjacent to x if * ( ) n y x ∈ Γ
. We say also that two points x and y of X are n-connected in X if there is a n-path between these two points. The equivalence classes for this relation are n-connected components of X . A subset X of E is said to be n-connected if it consists of exactly one n-connected component. The set of all n-connected components of X which are n-adjacent to a point x is denoted by [ ] , n C x X . In order to have a correspondence between the topology of X and the topology of X , we use n-adjacency for X and nadjacency for X , with ( , ) n n equal to (8; 4) or (4; 8).

Informally, a simple point p of a discrete object X is a point which is inessential to its topology. In other words, we can remove p from X without changing its topology

. A point x X ∈ is said simple if each n- component of X contains exactly one n-component of { } \ X x and if each n -component of { } X x ∪ contains exactly one n -component of X . Let X E
⊂ and x E ∈ , two connectivity numbers defined as follows ( # X = cardinality of X ):

( )

* 8 ( , ) # , n T x X C x x X   = Γ ∩   ; ( ) * 8 ( , ) # , n T x X C x x X   = Γ ∩   .
The following properties allows us to locally characterize simple points [START_REF] Kong | Digital topology: introduction and survey[END_REF][9] hence to implement efficiently topology preserving operators:

x E ∈ is simple for ( , ) 1 X E T x X ⊆ ↔ = and ( , ) 1 T x X = .
The homotopic alternating sequential filter is a composition of homotopic cuttings and fillings by balls of increasing radius. It takes an original image X and a control image C as input, and smoothes X while respecting its topology and geometrical constraints implicitly represented by C . A simple illustration is given by figure 1. Smoothed image (b) is obtained using HAS filter with a radius equal to five and four connectedness ( 4 Γ ).More example can be found in [5]. ...

C D D C D C n n n

HASF HF HC HF HC =

In the previous formula, C n HC (i) refers to homotopic cutting of X by n B with constraint set C and D n HF (ii) refers to homotopic filling of X by n B with constraint set D . These two homotopic operators can be defined as follows:

( ) ( ) , C n HC X H Y V = * With ( ) ( ) ( ) ( ) { , n n Y H X X C V Y X ε δ = ∪ = ∩ (i) ( ) ( ) , D n HF X H Z W = With ( ) ( ) ( ) { , ( ) n n Z H X X D W Y X δ ε = * ∩ = ∪ (ii)
We recall that ( ) , H Z W is an homotopic constrained thinning operator. It gives the ultimate skeleton of Z constrained byW . The ultimate skeleton is obtained by selecting simple point in increasing order of their distance to the background thanks to a pre-computed Euclidian distance map [START_REF] Danielsson | Euclidean distance mapping[END_REF]. We recall also that ( )

, H Y V *
is an homotopic constrained thickening operator. It thickens the set of Y by iterative addition of points which are simple for Y and belong to the set V until stability.

We have provided, in this section, the theoretical underpinnings for studying topological transforms. We introduced also the homotopic alternating sequential filter which constitutes base for topological smoothing.

Parallelization Strategy

In this section, we start by defining the class of topological operators. We also present our motivation to parallelize these algorithms on parallel shared memory machines. Then, we will introduce different steps of our approach after making a brief classification over existing strategies. We will focus especially on distribution phase and tasks scheduling over different processors. Scheduling and merging algorithms are presented and discussed. To illustrate both algorithms, scenarios are also introduced and discussed.

Class of topological algorithms

In 1996, Bertrand and Couprie [START_REF] Bertrand | Topological approach to image segmentation[END_REF] introduced connectivity numbers for grayscale image. These numbers describe locally (in a neighborhood of 3x3) the topology of a point. According to this description any point can be characterized following its topology. They also introduced some elementary operations able to modify gray level of a point without modifying image topology. These elementary operations of point characterization present a fundamental link between large class of topological operators including, mainly, skeletonization and crest restoring algorithms [START_REF] Couprie | Topological operators for greyscale image processing[END_REF]. This class can also be extended, under condition, to homotopic kernel and leveling kernel transformation [START_REF] Bertrand | Image segmentation through operators based on topology[END_REF], topological watershed algorithm [START_REF] Bertrand | On Topological Watersheds[END_REF] and topological smoothing algorithm [START_REF] Couprie | Topology preserving alternating sequential filter for smoothing 2D and 3D objects[END_REF] which is the subject of this article. All mentioned algorithms get also many algorithmic structure similarities. In fact associated characterizations procedures evolve until stability which induce common recursion between different algorithms. The grey level of any point can also be lowered or enhanced more than once. Finally, all mentioned algorithms get a pixel's array as input and output data structure. It is important to mention that, to date, this class has not been efficiently parallelized like other classes as connected filter of morphological operator which recently has been parallelized in Wilkinson's work [START_REF] Wilkinson | Concurrent Computation of Attribute Filters on Shared Memory Parallel Machines[END_REF]. Parallelization strategy proposed by Seinstra [START_REF] Seinstra | A software architecture for user transparent parallel image processing[END_REF] for local operators and point to point operators can also be cited as example. For global operators, Meijster strategy [START_REF] Meijster | A general algorithm for computing distance transforms in linear time[END_REF] shows also consistence. Hence the need of a common parallelization strategy for topological operators that offers an adapted algorithm structure design space. Chosen algorithm structure patterns that will be used in the design must be suitable for SMP machines.

In reality, although the cost of communication (Memory-processor and inter-processors) is high enough, shared memory architectures meet our needs for different reasons: (i) These architectures have the advantage of allowing immediate sharing of data with is very helpful in the conception of any parallelization strategy (ii) They are non-dedicated architecture using standard component (processor, memory ...) so economically reliable (iii) They also offer some flexibility of use in many application areas, particular image processing.

Split Distribute and Merge Strategy

In practice, most effective parallel algorithm design might make use of multiple algorithm structures thus proposed strategy is a combination of the divide and conquer pattern and event-based coordination pattern, see figure 2. Hence the name that we have assigned: SD&M (Split Distribute and Merge) strategy. Not to be confused with mixed-parallelism approach (combining data-parallelism and task-parallelism), it is important to mention that our strategy (i) represents the last stitch in the decomposition chain of algorithm design patterns and it provides a fine-grained description of topological operators parallelization while mixed-parallelism strategy provides a coarse-grained description without specifying target algorithm. (ii) It covers only the case of recursive algorithms, while mixed-parallelization strategy is effective only in the linear case. (iii) It is especially designed for shared memory architecture with uniform access.

Split phase

The Divide and Conquer pattern is applied first by recursively breaking down the problem into two or more sub-problems of the same type, until these become simple enough to be solved directly. Splitting the original problem take into account, in addition to the original algorithm's characteristics (mainly topology preservation), the mechanisms by which data are generated, stored, transmitted over networks (processor-processor or memoryprocessor), and passed between different stages of computation.

Distribute phase

Work distribution is a fundamental step to assure a perfect exploitation of multi-cores architecture's potential. We'll start by recalling briefly some basic notion of distribution techniques then we introduce our minimal distribution approach that is particularly suitable for topological recursive algorithms where simple point characterization is necessary. Our approach is general and applicable to shared memory parallel machines. Critical cases are also introduced and discussed.

Indeed there are two main types of scheduler. There are those designed for real-time systems (RTS). In this case, the most commonly approaches used to schedule real-time task system are: Clock-Driven, Processor-Sharing and Priority-Driven. Further description of different scheduling approaches can be found in [START_REF]Imprecise and Approximate Computation[END_REF][START_REF]Foundations of Real-Time Computing: Scheduling and Resources Management[END_REF][START_REF] Leung | REAL-TIME SCHEDULING ANALYSIS report[END_REF]. According to [START_REF] Leung | REAL-TIME SCHEDULING ANALYSIS report[END_REF] the Priority-Driven is far superior the other approaches. These schedulers must provide an operational RTS: completed work and delivered results on a timely basis. Other schedulers are designed for Non Real-time system. In this case, schedulers are not subject to the same constraints. Thus, "Symmetric Multiprocessing" scheduler distributes tasks to minimize total execution time without load balancing between processors, see figure 3 (a). On multi-core architectures, this can lead to high occupancy rate of one processor while the others are free. We propose a novel tasks scheduling approach to prevent improper load distribution while improving total execution time, see figure 3 (b). In literature, there are several schedulers that provide a balanced distribution of tasks such as RSDL "Rotating Staircase Deadline" [START_REF] Kolivas | RSDL completely fair starvation free 64 interactive cpu scheduler[END_REF] which incorporates a foreground-background descending priority system (the staircase) with run-queue managing minor and major epochs (rotation and deadline). Other scheduler, as CFS "Completely Fair Scheduler" [START_REF] Molnar | Modular Scheduler Core and Completely Fair Scheduler[END_REF], shows consistence. It handles resource allocation for executing processes, and aims to maximize overall CPU utilization while maximizing interactive performance. These schedulers are based on tasks uniformity principle. Through the tasks homogeneity, better distribution can be achieved and total execution time reduced. Unfortunately, these schedulers are not available in all operating system versions especially for small system. Based on the same principle of tasks uniformity, we propose a new scheduling algorithm, simpler to implement and more adapted to topological algorithm implementation.

I7

Let be a basic non-preemptive scheduler 'Basic-NPS', { }

a T P T ≠ ∅ ∧ ≠ ∅ then : x y T P ⇒ ; x T T T ∈ y a P P ∈ .
In this scheduler, each processor will treat at maximum max / max( )

i i j T m t t p P = → ≤ tasks with { } 1, 2,..., j n ∈
. Then, the worst case to process T is { }

1 ... ( ) max max ,..., max T T k i j i K T T p T p < < = → → . As proof, let suppose that it exist a set ( ) L T as ( ) ( ) L T K T ≥ ∑ ∑
. As 'Basic-NPS' manage ( ) L T and ( ) K T , so we can introduce the following: ( )

L T m ≤ and ( ) K T m ≤ . Thus, if ( ( ) ( )) L T K T ≥ ∑ ∑ then it exists at least one task { } l , with ( ) k K T ∈ , such as: ( ) A B C ∧ ∧ with ( ( )), ( ( )), ( ) A l L T B l K T C l k = ∈ = ∉ = >
. This is impossible according to the definition of ( ) K T

which was defined as the worst case.

Algorithm 1: Scheduling policy EndIf 20. EndWhile Algorithm 1 describes 'Basic-NPS' policy. The first step consists on asking operating system to determine the number of available processor. Depending on this number, algorithm will generate process. One active process will be assigned for each available processor. These new processes will belong to the SHED_FIFO class in order to ensure preemption and especially to avoid context switching. Process will only stop running if work is complete or less frequently when another process, belonging to the same class, with higher priority requesting processor. The global execution will stop if there no more task to process.

Merging phase

The key problem of each parallelization is merging obtained results. Normally this phase is done at the end of the process when all results are returned by all threads what usually means that only one output variable is declared and shared between all fighting threads. But as we mentioned in section 3.1, we are dealing with a dynamic evolution and if we take into account different steps of simple point detection then pixel characterizations, we can plan the following: The original shared data structure, containing all pixels, is divided

into n research zones { } 1 2
, ..., n z z z . We associate one thread from the following list { } 1 2

, ..., n T T T to each zone. Each thread can browse freely its zone and if it detects target pixel types, it lowers characterized pixel and it pushes its eight neighbors in one of the available FIFO queues. A queue is said available if only one thread (owner) is using it. One queue cannot be shared by more than two threads so if no queue is available, threads can create a new one and become owners.

Since two threads finished, they directly merge and a new thread is created and then same process is lunched again. New created thread will inherit queue shared between his parents. Thus it can restart research. It is also important to mention that there is no hierarchical order in thread merging, only criteria is finishing time. We mention also that one neighbor cannot be inserted twice. It is a precaution in order to minimize consumed cache. More formal description of merging techniques is given in by algorithm 2.

It is important to highlight similarity and difference that may exist between our merging algorithm and KPN [START_REF] Kahn | The Semantics of a Simple Language for Parallel Programming[END_REF]. In effect, both are deterministic and do not depend on execution order. But KPN algorithm may be executed in sequentially or in parallel with the same outcome while our merging algorithm is designed only for parallel execution. KPN support recurrence and recursion while our merging algorithm support only recursion.

In large scale application, KPN showed consistence. Examples include Daedalus project [START_REF] Nikolov | Daedalus: Toward Composable Multimedia MP-SoC Design[END_REF] where generated KPN models are used to map process into FPGA architecture. Ambric architectures [START_REF] Halfhill | Ambric's New Parallel Processor[END_REF] implement also a KPN model using bounded buffers to create massively DMP Machines based on structural object programming model.

In a narrower framework limited to simple point characterization, the implementation of such a model will be very expensive and it would be better to find an easier and more specific algorithm. We have introduced, in this section, three necessary steps to implement our parallelization strategy (SDM). It is important to mention that some similarity may exist between our split/merging phases and alpha-extension/betareduction phases from structural perspective. Actually both approaches intended to put in place more guarantees that the parallelism will actually be met. But uses contexts are different. In effect, Jean Paul Sansonnet [START_REF] Giavitto | Introduction à 8 1/2 Rapport interne[END_REF][START_REF] Giavitto | 8 1/2 : data-parallélisme et data-flow[END_REF][START_REF] Mahiout | Distribution and scheduling data-parallel dataflow programs on massively parallel architectures, SMS-TPE '94: Software for Multiprocessors and Supercomputers[END_REF] team introduced alpha-extension (diffusion) and beta-reduction (merging) notions for stream manipulation in the framework of Declarative Data Parallel language definition and there techniques cannot be applied without a scalar function. While our proposal is restricted to topological characterization in the framework of topological operator's parallelization and no scalar function is required during the application of these two phases.

Parallel smoothing filter

In this section we start by analyzing overall structure of original algorithm. Then we continue with the parallelization of Euclidean distance, thinning and thickening algorithm. We conclude by a performance analysis of the entire smoothing topological operator. Obtained execution time, efficiency, speedup and cache misses will be introduced and discussed. As we have shown in Section 2, smoothing algorithm receives as input a binary image and maximum radius. It uses two procedures for homotopic opening and closing, see figure 5 (a) (b). The call is looped to ensure an ongoing relationship between input and output. The opening process is a consecutive execution of erosion, thinning, dilatation and thickening. While closure procedure ensures the same performance of the four consecutive functions with single difference: the erosion instead of dilatation. Thinning and thickening ensure the topological control of erosion and dilatation. This control is based on researching and removing of all destructible points. When destructible point is deleted, its neighbors are reviewed to ensure that they are not destructible either.

Euclidean distance computing 4.1.1 Study on Euclidian Distance algorithms

During previous evaluation, 4SED [START_REF] Danielsson | Euclidean distance mapping[END_REF] algorithm was used for Euclidean distance computation. So we are looking for another algorithm that is faster, and parallelizable. New algorithm must have an Euclidean distance computation error less than, or equal to, that produced by 4SED in order to maintain homotopic characteristics of the image. In literature, several algorithms for Euclidean distance computing exist. Lemire [START_REF] Lemire | Streaming Maximum-Minimum Filter Using No More than Three Comparisons per Element[END_REF] and Shih [START_REF] Shih | Fast Euclidean distance transformation in two scans using a 3x3 neighborhood[END_REF] algorithms are bad candidates because Lemire's algorithm does not use Euclidean circle as structuring element. Then homotopic property will not be preserved. Shih's algorithm has a strong data dependency which penalizes parallelization. In [START_REF] Cuisenaire | Fast Euclidean Distance Transformation by Propagation Using Multiple Neighborhoods[END_REF], Cuissenaire propose a first algorithm for Euclidian distance computing, called PSN "Propagation Using a Single Neighborhood" that uses the following element structure:

2 2 4 ( ) 1 x x y x d p q q p q p   = ∨ - + - <     (i)
He also proposes a second algorithm, called PMN "Propagation Using Multiple Neighborhood" that uses eight neighbors. In [START_REF] Cuisenaire | Fast and exact signed Euclidean distance transformation with linear complexity[END_REF], he also proposes a third algorithm with 3/2 ( ) n ο complexity, which offers an accurate computation of the Euclidean distance. Only drawback of this third algorithm is computation time which is very important and goes beyond the two algorithms mentioned above. Even if computing error produced by PSN is greater than computing error produced by PMN, it is comparable to that produced by 4SED. Low data dependence and ability to operate on 3D images, makes PSN algorithm a potential candidate to replace 4SED.

Meijster [START_REF] Meijster | A general algorithm for computing distance transforms in linear time[END_REF] proposes an algorithm to compute exact Euclidean distance. Algorithm complexity is ( ) n ο and it operates in two independent, but successive, steps. First step is based on looking over columns then computing distance between each point and existing objects. Second step includes same treatment looking over lines. It is important to note that strong independence between different processing steps and computing error equal to zero makes Meijster algorithm another potential candidate to replace 4SED. Algorithm is also able to operate on 3D images. Theory analysis of Meijster and Cuissenaire algorithms can be found in Fabbri's work [START_REF] Fabbri | 2D Euclidean distance transform algorithms: A comparative survey[END_REF].

In the following, we propose first analysis based on different algorithms implementation in order to compare between them. We have implemented 4SED algorithm using a fixed size stack. This stack uses a FIFO queue and it has small size while 4SED algorithm does not need to store temporal image. Results are directly stored into the output image, we will retain this implementation because 4SED assessment serve only as reference for comparison. For PSN implementation, we used stacks with dynamic sizes. Memory is allocated using small blocks defined at stack creation. When an object is added to queue, algorithm will use available memory of last block. If no space is available, a new block is allocated automatically. Block size is proportional to image size (N x M / 100). Finally we used a simple memory structure to implement Meijster algorithm. A simple matrix was used to compute distance between points and object of each column and three vectors were used to compute distance in each line. We recall that this comparison is done in order to select the best algorithm among three candidates. If we assume that minimum distance of an empty group K is ∞ and z K ∀ ∈ , we have ( )

y z + ∞ = ∞ then
( ) EDT p formula can be written as follow: Thus we can split the Euclidian distance transform procedure into two steps. The first step is to scan columns and compute EDT for each column y . Second step consists on repeating the same procedure for each line. In the following we start by detailing these two steps: In the first step ( , )

x G p y can be computed through the two following sub functions:

( , ) min : ( , )

T x x x x G p y p b b b y = - = , ( , ) min : ( , ) B x x x x G p y b p b b y = - = with 0 x b n ∀ ≤ ≤ . To compute ( , ) T x
G p y and ( , )

B x
G p y , we scan each column y from top to bottom using the two following formula:

( , ) ( , 1) 1

T x T x G p y G y p = -+ ( , ) ( , 1) 1 B x B x
G p y G y p = + + . Thus sequential algorithm of the first step can be written as follows. Complexity order is ( ) n m ο × .

Algorithm 3: Meijster O.V [1 st Step]

Data : m:colums, n:lines, b:image

1. Forall [ ] 1 .. 0 - ∈ m y do 2. If ( ) B y ∈ , 0 then [ ] 0 .. 0 = y g 3.
else [ ] ∞ = y g .. 0

endif

5.

/* G T */ 6. for ( ) For each row u , we note that there is, for the same point p , the same value of ( , ) f p y for different values of y , so we can introduce the concept of "region of column ".

1 = x to ( )
Let S be the set of y points such that ( , ) f p y is minimal and unique. The formula of S , 0 y u

∀ ≤ ≤ , is ( ) m i n : ( , ) ( , ) p S u y f p y f p i = ≤ . 0 i u u m ∀ ≤ ≤ ∧ ≤ . Let T be

Let

( , ) Sep i u be the separation between regions of i and u , defined by: ( , ) ( , )

f p i f p u ≤ 2 2 2 2 ( ) ( , ) ( ) ( , ) y x y x p i G p i p u G p u ⇔ -+ ≤ - + 2 2
( , ) ( ) / 2( 1)

x p y Sep i u u i Dif u p ⇔ = -+ -= With 2 2
( ( , ) ( , ) )

x x Dif G p u G p i = -
. Thus lines will be processed, from left to right then from right to left. During the first term, from left to right, two vectors S and T will be created. These two vectors will contain respectively all regions and all intersections. During the second treatment, from right to left, we compute f for each value of S . f is also computed for each respective values of T . Algorithm 4 is associated to second step. For the first term, complexity order is 2( ) q m u + -whereas complexity order of the second term is only m .

The independence of data processing between rows and columns is the key to apply of SDM parallelization strategy. In the first stage, column processing, we can define data interdependence by the following equation:

{ }

( , ) min ( , ), ( , )

x T x B x G p y G p y G p y = { } ( , ) 0 ( , ) ( , ) x T x if p y B T x G p y else G p y ∈ ⇔ = { } ( , ) min ( 1, ), ( , ) B x B x T x G p y G p y G p y ⇔ = +
It follows that values of each column y of G, depends only on lines: x p , 1

x p + and 1

x p -. Similarly, at the second stage, we can introduce the following interrelationship:

( ) ( , ( ))

p Edt p f p S q = . Then (0 ), (0 ) ( ) y u i u u m ∀ ≤ ≤ ≤ ≤ Λ < , ( ) p S u = min : ( , ) ( , ) y f p y f p i ≤ . Thus, if ( ( ))
p u T q = so ( 1) q q =which imply the following: ( ) ( ( ), ) 1

x p p p T u Sep S q u = + .

Algorithm 4: Meijster O.V [2 nd Step]

Data : b:image, g: G_Table, m: columns, n:lines Result :

1. Forall [ ] 1 .. 0 - ∈ n x do 2. 0 = q 3. [ ] 0 0 = s 4.
[ ] 0 0 = t

5.

/* First part */ 6. for ( )

1 = u to ( ) 1 - m do 7. [ ] [ ] ( 0) (( , ), ) A q f x t q s q   = ≥ Λ   8.
[ ]

(( , ), ) B f x t q u = 9. while ( ) A B > then ( 1) q q ← + 10.
end while 11. 12.

if ( 0) q < then ( 0) q ← 13. [ ] ( 0 ) s u ← 14. else [ ] ( , , ) 1 w Sep s q u x ← + 15. if ( ) w m < then ( 1) q q ← + 16.
[ ] for ( )

s q u ← 17. [ ] t q w ←
1 u m = -to ( ) 0 do 23. [ ] [ ] , (( , ), ) Edt x u f x u s q = 24. if [ ] ( ) u t q = then ( 1) q q ← -

endif 26. endfor 27. end forall

According to this formalization, values of ( , )

f p i and ( , )

x Sep i u are independent of modified data. So using two vectors S and T , a private variable q for each line ensures complete independence in writing. We start applying the splitting step by sharing the columns and lines processing between multiple processors. A thread can process one or more columns and the number of threads used will depend on the number of processors. The results returned by all threads in this first stage will be merged in order to start lines processing. In the following we introduce the parallel version of Meisjter algorithm for both steps. Associated algorithm complexity is ( )

( ) / n m N ο × . (
) n m × refers to image size and N refers to the number of processors.

Algorithm 5: Meijster P.V [1

y m y y t

= < = + do 2. If ( ) B y ∈ , 0 then [ ] 0, 0 g y ← 3. else [ ] 0, g y ← ∞ 4. endif 5.
/* G T */ 6. for ( ) Proposed parallel version of Meijster algorithm was implemented in C using OpenMP directives. Speedup for numbers of threads equal to 1, 2, 4, 8, and 16 were determined. The efficiency measure Ψ ( n ) is given by the following formula with n the number of processors: Ψ ( n ) = seq. time /( n *para. time) (ii) Times were performed on eight-core (2× Xeon E5405) shared memory parallel computer, on Intel Quad-core Xeon E5335, on Intel Core 2 Duo E8400 and Intel mono-processor Pentium 4 660. The minimum value of 5 timings was taken as most indicative of algorithm speed. More information about architectures characteristics are given in Section 4.

1 = x to ( ) 1 - n do 7. if [ ] B y x ∈ , then [ ] , 0 g x y ← 8. else [ ] [ ] , 1

Algorithm 6: Meijster P.V [2 nd

Step]

1. For max ( , , ) x t x n x x t = < = + do 2. 0 = q ; [ ] 0 0 = s ; 3. [ ] 0 0 = t ; 4.
/* First part */ 5. for ( )

1 = u to ( ) 1 - m do 6. [ ] [ ] ( 0) (( , ), ) A q f x t q s q   ← ≥ Λ   7.
[ ]

(( , ), ) B f x t q u ← 8. while ( ) A B > do ( 1) q q ← + 9.
end while 10. 11.

if ( 0) q < then ( 0) q ← 12.
[ ]

( 0 ) s u ← 13. else [ ] ( , , ) 1 w Sep s q u x ← + 14. if ( ) w m < then ( 1) q q ← + 15.
[ ] for ( ) The measurements were done on 2D binary image (512*512). If we can get a satisfactory outcome for this standard, it will be the same for smaller size images. View cache size limits, larger image will not be tested.

s q u ← 16. [ ] t q w ←
1 u m = -downto ( ) 0 do 22. [ ] [ ] , (( , ), ) Edt x u f x u s q ← 23. if [ ] ( ) u t q = then ( 1) q q ← -
Figure 7 shows that number of instructions to compute Euclidian distance drops from an average of 9.5x10 8 using 4SED algorithm down to 7.6x10 8 ms with Meijster algorithm. Despite the passage from a sequential version running on single core to a parallel version running on 8 processors, acceleration is only multiplied by 1.6 as shown in figure 8 (a).This can be explained by the choke point between columns processing and lines processing. Waiting time between these two treatments significantly penalizes acceleration. Figure 8 (b) shows that efficiency variation depends on the number of threads. It is also proportional to the number of processors. Moving to 3, 5 or 7 threads (odd number) decreases significantly the efficiency which reaches its maximum each time that the number of threads is equal the number of processors.

Thinning and thickening computing

Algorithms of thinning and thickening are almost the same. The only difference between them is the following: in thinning algorithm, destructible points are detected then their values are lowered. In thickening algorithm, constructible points, are detected then their values are increased. For parallelization, we will apply the same techniques introduced in [START_REF] Mahmoudi | Parallel image thinning through topological operators on shared memory parallel machines[END_REF]. We propose a similar version using two loops. Target points are initially detected then their value lowered or enhanced according to appropriate treatment. The set of their eight (or four) neighbors are copied into a "buffer" and rechecked. This treatment is repeated until stability. In the following, we present an adapted version of Couprie's thinning algorithm.

Algorithm 7: Adapted Version Thinning Algo.

← -

Endwhile

Unfortunately direct application of introduced parallel processing is not possible with the set of all points. Some points, called critical points, cannot be eliminated in parallel because initial topology of the image may be broken. Figure 9 illustrates this case: Critical points of an input image (a) are identified in (b). If these points are deleted in one iteration (c) topology necessary is broken (d).

To resolve this problem, we propose that research areas assigned to each thread must be composed of at least six lines (of the image). Each thread will use two buffers to treat each three lines thus four buffers are used to treat six lines as shown in 9 (e). In some cases, a neighbor of a destructible point is detected on the border of a contiguous area. To prevent that such neighbor escape to recheck, it must be injected to buffer of the right thread. Let's suppose that a point p ∈Z 2 is considered as destructible by T2, so its value will be lowered and its four neighbors { }

1 2 3 4 , , , v v v v should be rechecked. Neighbors { } 1 2 4
, , v v v belong to Z 2 so they will be push in T2 buffers. The neighbor { } 3 v belongs to Z 3 so it will stack T3 buffers.

Performance evaluation of introduced adapted version of Couprie's algorithm is shown in figure 10. On eight cores architecture, acceleration does not exceed 3.4. Such moderate result can be explained by critical borders processing. Regarding efficiency, the best performance is achieved when the number of thread is equal to the number of processors. If this equality is not ensured, the efficiency decreases. The problem threads' add number still persists. The next step is to combine the parallel version of Meijster algorithm and the adapted version of Couprie's algorithm to build the parallel processing of topological smoothing.

Global analyses

In this section, we present a global evaluation of the parallel smoothing operator. We start by presenting performance evaluations in terms of acceleration and efficiency. Then we evaluate cache memory consumption.

Execution time

We implemented two versions of the proposed parallel topological smoothing algorithm, the first one using 'Symmetric Multiprocessing' scheduler and the second one using 'basic-NPS' scheduler. Wall-clock execution times for numbers of threads equal to 1, 2, 4, 8, and 16 were determined. The minimum value of 2 timings was taken as most indicative of algorithm speed. The measurements were done on 2D binary image (512*512). Results of the second implementation on the eight-core are shown in the following figure. We note that number of instructions drops from an average of 1879 x10 8 FI with a single thread down to 1652 x10 8 ms with 8 threads. As expected, the speed-up for the second implementation using 'basic-NPS' scheduler is higher than for the one using "Symmetric Multiprocessing" scheduler, thanks to balanced distribution of tasks. A remarkable result about speedup is also shown in figure 12 (a). In fact, speed-up increases as we increase the number of threads beyond the number of processors in our machine (eight cores). In the first implementation, using "Symmetric Multiprocessing" scheduler, the speedup at 8 threads is 1.9 ± 0.01. However, for the second implementation, using our scheduler, the speedup has increased to 5.2 ± 0.01. Another common result between different architecture is stability of execution time on each n-core machine since the code uses n or more threads.

Cache Memory Evaluation

As memory access is a principal bottleneck in current-day computer architectures, a key enabler for high performance is masking the memory overhead. If we starts from basic theory that two classic cache design parameters dramatically influence the cache performance: the block size and the cache associativity. So the simplest way to reduce the miss rate is to increase the block size even it increases the miss penalty. The second solution is to decrease associatively in order to decrease hit time thus to retrieve a block in an associative cache, the block must be searched inside of an entire set since there is more than one place where the block can be stored.

Unfortunately, we are dealing with non-reconfigurable architectures with caches whose associativity and block size are predefined by the manufacturer. Nowadays, new approaches to reduce cache miss are developed such as taking advantage of locality of references to memory or using aggressive multithreading so that whenever a thread is stalled, waiting for data, the system can efficiently switch to execute another thread. Despite their power, the application of both approaches remains limited. In fact, applications of locality approach still experimental even with Larrabee technology introduced by Intel. And the aggressive multithreading approach has been specially designed for graphics processing engines, which manage thousands of in-flight threads concurrently. So it is not recommended for general SMP machines with limited number of processors and threads. With all these limitations, the most intuitive solution is to rely on the scheduling. Thanks to our basic-NPS scheduler, we have balanced the charges then prevent context switching thus we minimize caches misses. In the following we present our experimental analysis. We consider a commonly used Intel processor configuration (More details are given by table 2). Number of processor varies from one to eight. The frequency varies between 1,73 GHz and 3,4 GHz. The L1 caches have at least a 32-byte block size, while capacity vary between 16 Kbytes and 32 Kbytes, and for the associativity, only eight ways is considered. The L2 caches have at least a 64-byte block size, while capacities vary between 512 Kbytes and 6 Mbytes, and the associativity varies between two and twenty four ways.

(A) (B) The scheduler relies on our basic-NPS scheduling policy. As a result of this experiment, see figure 13 (A-1), we found that three performance regions are clearly evident: In the leftmost region, as long as the cache capacity can effectively serve the growing number of threads, increasing the number of threads improves performance, as more processors are utilized. This area is generally identified as cache-efficiency zone. At some point, the cache becomes too small for the growing stream of access requests, so memory latency is no longer masked by the cache and instruction cache misses reduce more moderately. As the number of available threads again increases, the multithread efficiency zone (on the right) is reached, where adding more threads improves performance up to the maximal performance of the machine, or up to the bandwidth wall. Balanced workloads offer higher locality and better exploit the cache and hence expand the cache efficiency zone to the right and up. An outstanding example is given by the following table which summarizes number of L1 instruction misses on Intel Dual Core T1400 architecture using SMP scheduling policy and Basic-NPS scheduling policy. We note that number of instruction misses drops from an average of 18844 L1 Instr. misses (using SMP) with two threads down to 6030 L1 Instr. misses (using Basic-NPS) usually with two threads. Here success rate is largely above the average of 50%. The same rate will be practically maintained when increasing the number of threads.

Figures 14 (A) and 16 (A-1) show so much load balancing and implicitly context switching between processes can affect performance in terms of reading data from caches. However, improvement in writing data, see Fig. 14 (B) and Fig. 16 (B-1), in two caches remains modest. When there are more computation instructions per memory access, performance climbs more steeply with additional threads. This is because as more instructions are available for each memory access, fewer threads are needed to fill the stall time resulting from waiting for memory.

Conclusion

Topological characteristics are fundamental attributes of an object. In many applications, it is mandatory to preserve or control the topology of an image. Nevertheless, the design of transformations which preserve both topological and geometrical features of images is not an obvious task, especially for parallel processing.

In this paper, we have presented a new parallel computation method for topological smoothing through combining parallel computation of Euclidean Distance Transform using Meijster algorithm and parallel Thinning-Thickening processes using an adapted version of Couprie's algorithm.

We have also presented a new parallelization strategy called SDM (Split Distribute and Merge). Proposed strategy is partially based on divide and conquers principle associated to event-based coordination techniques. Further than smoothing operator, SDM Strategy can be applied for a large class of topological operators as we shown in section 3.1. In addition to identified conditions during splitting step, we introduced an adapted scheduler called basic-NPS (Basic -Non Preemptive Scheduler) able to distribute in balanced way a set of active tasks on available processors. Finally we introduced an adapted merging policy designed especially for dynamic system evolving until stability.

Parallel topological operator computation poses many challenges, ranging from parallelization strategies to implementation techniques. We tackle these challenges using successive refinement, starting with highly local operators, which process only by characterizing points and then deleting target pixels, and gradually moving to more complex topological operators with non-local behavior. In future work, we will study parallel computation of the topological watershed [START_REF] Bertrand | On Topological Watersheds[END_REF].

Fig. 1 :

 1 Fig.1: (a) Input image (b) Smoothed imageBased on this filter, Authors[START_REF] Couprie | Topology preserving alternating sequential filter for smoothing 2D and 3D objects[END_REF] introduce a general smoothing procedure with a single parameter to control smoothing degree. Let C X ⊆ , r ∈ ℕ and D X ⊆ with X any finite subset of E . The homotopic alternating sequential filter ( ) HASF of order n , with constraint sets C and D , is defined as follows:

Fig. 2 :

 2 Fig.2: SDM strategy classification

Fig. 3 :

 3 Fig.3: (a) Non-real time symmetric task distribution (b) Task distribution based upon uniformity principle

Fig. 4 :

 4 Fig.4: Illustration of merging phase with four threads In Figure 4, we give an illustration of the merging algorithm with four threads. The original shared data structure is divided into 4 research areas{

Fig. 5 :

 5 Fig.5: Overall structure of the original smoothing algorithm

Fig. 6 :

 6 Fig.6: Comparison between different EDT algorithms

Figure 6

 6 Figure 6 describes obtained results by different implementations on single processor architecture P4. During this evaluation we used binary test image (200x200). We have also varied ball radius. We used Valgrind software to evaluate different designs. Callgrind tool returns the cost of implementing of each program by detecting IF (Instruction Fetch). Results show that PSN algorithm is the most expensive in all cases (for any radius). Meijster algorithm is moderately faster than 4SED. The output images returned by Meijster algorithm hold the best visual quality while Euclidean distance computation error is almost zero thus our efforts will be brought on Meijster algorithm parallelization.

Fig. 7 :Fig. 8 :

 78 Fig.7: Instruction distribution (Meijster Alg)

Fig. 9 :

 9 Fig.9: (a) (b) (c) (d) Critical point illustration (e) Research Area Assignment Through this organization threads can start running in parallel on Z 11 , Z 21 and Z 31 .Once processing is completed threads can restart running on Z 12 , Z 22 and Z 32 . In some cases, a neighbor of a destructible point is detected on the border of a contiguous area. To prevent that such neighbor escape to recheck, it must be injected to buffer of

Fig. 10 :

 10 Fig.10: (a) Performance Evaluation of Couprie's Algorithm (b)Efficiency Evaluation of Couprie's Algorithm.

Fig. 11 :

 11 Fig. 11: Number of Instructions and tasks distribution using 'Basic-NPS'

Fig. 12 :

 12 Fig. 12: (a) Global Performance improvement (b) Global efficiency improvement For better readability of our results, we tested also efficiency of our algorithm on various architectures (see figure 12(b)) using the ( ) n ψ formula introduced earlier. For parallel time ratio we used best obtained time with 8 threads ('basic-NPS' scheduler).

Fig. 13 :

 13 Fig. 13: (A) Instruction -L1 misses

Fig. 14 :

 14 Fig. 14: (A) Data Read -L1 misses (B) Data Write -L1 misses

  the ball of radius r centred on x , defined by , where d is a distance on Ε . We denote by r B the map which associates to each x in Ε the ball ( ) . The Euclidean distance d on Ε is defined by: E denotes the set of all subsets of E . Let r be an integer, the dilation by r B is the operator r

	1 ( , ) 2 x x with { ( ) r B x y E d x y r , ( , ) = ∈ ≤ 1/ 2 } 2 2 ( , ) A B d x y   = -  with
	A ( =	1 x y 1 -	)	and	B ( =	2 x y -	2	)
									δ defined by ( ) r X δ	= ∪	( ) B x x X r ∈	( ) X P E ∀ ∈	. The
	ball r B is termed as the structuring element of the dilation. The erosion by r B is the operator r ε defined by
	duality: r ε	r = * . δ

i x ∈ ℤ . Let x ∈ Ε , r ∈ ℕ , we denote by ( ) r B x r B x . An operator on E is a mapping from ( ) P E into ( ) P E , where ( ) P Now, we introduce notion of simple point which is fundamental for the definition of topological operators in discrete spaces. We give a definition of local characterization of simple points in 2 E = ℤ . Let consider two neighbourhoods relations 4 Γ and 8 Γ defined for each point x E ∈ by:

  will detect new target point so it will lower characterized pixel (in 3 z ) then push neighbors in 1 F before continue browsing. 3 T don't need to create new FIFO queue since 1 F is available. 1 T and 3 T will repeat this procedure twice. Since they finish browsing, they merge and new thread 5

				T is born. 5 T will start browsing only
	1 F . Since it detect new target point so it will lower characterized pixel (in 5 z = 1 z + 3 z ) then push neighbors in
	3 F that it has created before continue browsing. Similarly 2 T and 4 T will generate the creation of 2 F and 6 T .
	Here 6 T will eventually merge with 5 T to give birth to 7 T . Finally there will be a single thread 7 T which will
	brows 3 F without detection any target points.
	Algorithm 2: Merging technique
	1. Z : Set of research zones
	2. T : Set of threads
	3.	_ FIFO Q : Set of available FIFO queues
	4.	T P : Target pixel type ; D P : Detected pixel
	5. For all zones (	) ∈ do: Z Z i
	6.	Parallel_browsing ( i T , i Z ) ;
	7. EndFor
	8. For each thread (	) ∈ do: T T i
	9.	If (pixel_caract( i T , T P )==True) then
	10.	modify_value( D P );
	11.	If ( (	_ FIFO Q ≠ ∅ then )
	12.	usedstatus(	_ j FIFO Q ,true) ;
	13.	insert_neighbors( i T , D P ,	_ j FIFO Q );
	14.	Else : add_new_fifo (	_ FIFO Q )
	15.		usedstatus(	1 FIFO Q + ,False) ; _ j
	16.		insert_neighbors( i T , D P ,	1 FIFO Q + ); _ j
	17.	EndIf ;
	18.	EndIf ;
	19. EndFor;

T is the first to detect target point (constructible, destructible …) so it lowers characterized pixel (in 1 z )

8 | P a g e and it pushes its eight (or four) neighbors in FIFO queue 1

F that it has created before continue browsing. Later, 3 T

  It can reach 64.44% of total time with a radius equal to 5. However time for topological characterization is only 8.89% for 1 million points. These finding remain the same if we increase image size. Beyond (512*512), computing time of point characterization becomes considerable.

		r=5	200x200 r=10	r= ∞	r=5	168x288 r=10	r= ∞
	EucDis	64.44	54.93	46.67	59.25	49.79	35.25
	(%)
	TopCar	8.89	13.89	18.15	11.58	16.50	24.03
	(%)

Tab.1: Time execution rate of E.D and topological characterization functions

A preliminary assessment of first implementation code, see Table

1

, shows that Euclidean distance computing (EucDis) takes more time than topological point characterization (Topcar). For an image of (200*200), computation time of E.D with an infinite radius is 46.67% while point characterization of 2.4 million points occupies only 18.15%. If we limit radius between 5 and 10, computation time of E.D. continues to increase.

  the set of points with coordinate greater than, or equal to, horizontal coordinate of the intersection with a region:

	( ) T u p	=	S e p	p	x	(	S	p	(	u	-	1 ) , ) 1 u + .

1. while (

	15.			x	←	( pop stack	2)
	16.			( ) v neighbors x ←
	17.			i ←	0
	18.			While ( 8) i < do
	19.						if [ ] ( v i stack 1) ∉	then
	20.						[ ] push v i stack ( , 1)
	21.						endif
	22.			endwhile
	23.	endwhile
	24.	max	iter		max	iter	1
		[ ] input x is destructible) do
	2.	( , push x stack	1)
	3.	x x ← +	1
	4. endwhile
	5. output input ←
	6. While (	1 stack ≠ ∅ ∧ ) (max	iter	0) > do
	7.	While (	1 stack ≠ ∅ do )
	8.			x	←	( pop stack 1)
	9.			if (	[ ] output x is destructible) then
	10.					[ ] output x	←	_ ( ) reduce pt x
	11.						( , push x stack	2)
	12.			endif
	13.	end while
	14.	While (	2 stack ≠ ∅ do )


Ramzi Mahmoudi
email: mahmoudr@esiee.fr

Mohamed Akil
email: akilm@esiee.fr

The Watershed Transform

Keywords: Watershed transform, Flooding, Path-cost minimization, Topology preservation, Local condition, Minimum spanning forest

In the framework of mathematical morphology, watershed transform (WT) represents a key step in image segmentation procedure. In this paper, we present a thorough analysis of some existing watershed approaches in the discrete case: WT based on flooding, WT based on path-cost minimization, watershed based on topology preservation, WT based on local condition and WT based on minimum spanning forest. For each approach, we present detailed description of processing procedure followed by mathematical foundations and algorithm of reference. Recent publications based on some approaches are also presented and discussed. Our study concludes with a classification of different watershed transform algorithms according to solution uniqueness, topology preservation, prerequisites minima computing and linearity.

INTRODUCTION

The watershed concept began with Maxwell [START_REF] Maxwell | On hills and dales[END_REF] who introduces the theory behind representing physical characteristics of a land by means of lines drawn on a map. He highlights relationships between the numbers of hills, dales and passes which can co-exist on a surface. Subsequently, through the work of Beucher and al. [START_REF] Beucher | Use of watersheds in contour detection[END_REF], watershed transform was introduced to image segmentation and nowadays it represents one of the basic foundations of image processing [START_REF] Beucher | The morphological approach to segmentation: the watershed transformation[END_REF].

In this framework, the most simplified description of the watershed approach is to consider a grayscale image as a topographic surface: the gray level of a pixel becomes the elevation of a point, the basins and valleys correspond to dark areas, whereas the mountains and crest lines correspond to the light areas. If topographic relief is flooded by water, watersheds will be the divide lines of the attraction's domains of rain falling over the region [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF] or sources of water springing from reliefs' peaks. Another synopsis has shown consistency is that topographic surface is immersed in a lake with holes pierced in local minima. Catchment basins will fill up with water starting at these local minima, and, at points where water coming from different basins would meet, dams are built. As a result, the topographic surface is partitioned into different basins separated by dams, called watershed lines. Figure 1 gives a very symbolic description of the mentioned approach. In fact, it shows trends in the use of watershed transform for image processing. Despite its simplicity, this concept has been formalized in different ways giving rise to several definitions of watershed transform. In the discrete case, which is our main interest in this paper, this problem is amplified since there is no unique definition of the path that the drop of water would follow. This led to a multitude of algorithm to compute watershed transform. Some of these algorithms don't even meet associated watershed definition. We also note that some definitions take the form of algorithm specification which makes the distinction between algorithm specification and implementation very complicated. This problem in literature has been partially resolved in Roerdink and al. [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF] ten years ago. Actually authors presented a critical review of several definitions of the watershed transform and the associated sequential algorithms. Even they discuss various issues which often cause confusion in the literature; they don't go further in the classification or comparison of different approaches. They instead focus on parallelization aspect. In other more recent publications, authors tentatively drawn a comparison chart of some watershed transform definition to serve their end goals: showing the relationships that may exist between some discrete definition of watershed [START_REF] Audigier | Watershed by image foresting transform, tie-zone, and theoretical relationship with other watershed definitions[END_REF] or showing that most classical watershed algorithm do not allow the retrieval of some important topological features of the image [START_REF] Najman | Watershed algorithms and contrast preservation[END_REF].

The purpose of this paper is to introduce an intensive study of all existing watershed transform in the discrete case: WT based on flooding, WT based on path-cost minimization, watershed based on topology preservation, WT based on local condition and WT based on minimum spanning forest. Indeed, for each approach, we start by giving informal definition, then we present processing procedure followed by mathematical foundations and the algorithm of reference.

Recent publications based on some approach are also presented and discussed. Our study concludes with a classification of different algorithms according to criteria of recursion, complexity, basins computing and topology preservation. This paper is organized as follows: in section 2, different approaches to compute watershed are presented. In section 3, we draw a comparison between the various presented algorithms. Finally, we conclude with summary and future work in section 4.

WATERSHED TRANSFORMATION

In this section, we propose a comparative study of different approaches to compute the watershed transform in the discrete case. The goal is to identify most suited watershed transform for parallel processing. This study can also be seen as an update of Roerdink work [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF] presented ten years ago. Indeed, for each approach, we present processing procedure, mathematical foundations and sequential algorithm. Recent publications based on some approach are also presented and discussed.

Watershed based on flooding

Based on flooding paradigm [START_REF] Soille | Morphological Image Analysis[END_REF][START_REF] Dougherty | Hands-on Morphological Image Processing[END_REF][START_REF] Najman | Watersheds, mosaics and the emergence paradigm[END_REF], the intuitive idea underlying this method comes from geography. Since grayscale image can be seen as topographic surface, the intensity of a pixel can be considered as the altitude of a point. Now, let immerge this surface in still water, with holes created in local minima. Water fills up basins starting at these local minima. As described in algorithm 1, the filling of basins is an iterative process that involves gradually raising the water level from Alt min to Alt max . Algorithm must, for each iteration, fill existing basins (extension regions) and possibly create new basins (new regions). We denote by LR the region list. Dams will be built where waters coming from different basins meet.

{ } Alt p f D p T Alt ≤ ∈ = ) ( / min
It define a recursion with the gray level Alt increasing from Alt min to Alt max , the basin associated with the minima of f are successively expanded. Let X Alt denote the union of the set of basins computed at level Alt. A connected component of the threshold set T Alt+1 at level Alt+1 can be either a new minimum or an extension of the basin in X Alt . The geodesic influence zone (IZ) of X Alt within T Alt+1 can be computed resulting in an update X Alt+1 . Let MIN h denote the union of all regional minima at altitude Alt. Then we can introduce the following definition:

Definition 2.1 (Flooding watershed) { } min min min ) ( / Alt Alt T Alt p f D p X = = ∈ = ) ( 1 1 1 Alt T Alt Alt X IZ MIN X ALT+ ∪ = + + The watershed ) ( f Wshed of f is the complement of max Alt X in D : max / ) ( Alt X D f Wshed =
Vincent and Soille [START_REF] Vincent | Watersheds in digital spaces -An efficient algorithm based on immersion simulations[END_REF] presented an original and efficient implementation (algorithm 2) of the flooding watershed. This implementation uses FIFO queue and it needs two steps:

(1) Sorting pixels in increasing order of grey values (Alt min ; Alt max )

(2) Flooding process: All nodes with grey level h are first given the initial label. Then those nodes that have labeled neighbors from the previous iteration are inserted in the queue, and from these pixels geodesic influence zones are propagated inside the set of initial pixels.

In their study [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF] They also propose some modification to implement the recursion (2.3) exactly. In line 13, all pixels with [ ] h p im ≤ have to be masked, the queue has to be initialized with basin pixels only (drop the disjunct lab[q] = wshed in line 15), the resetting of distances (line 50) has to be done in line 14, and the propagation rules in lines 32-47 have to be slightly changed. [ ]

; WSHED p lab ← 42. end_if 43.
else if ((lab[q]=MASK) and (dist[q]= 0)) then //q is plateau pixel 44.

[ ]

; 1 + ← curdis q dist 45.
) [ ]

; curlab r lab ←
From the introduction about immersion simulation above, we can see that the level-by-level method during the flood procedure is uniform. Unfortunately, this can cause over-segmentation in several cases. Based on this original simulating immersion, Shengcai and Lixu [START_REF] Shengcai | A Novel Implementation of Watershed Transform Using Multi-Degree Immersion Simulation[END_REF] propose, in 2005, a novel implementation using multi-degree immersion. To our knowledge, it is the last update that can be found in the literature, the proposed implementation resists to oversegmentation problem effectively. It starts by redefining the threshold set of f at level Alt.

Instead of the original formula (2.1), they propose the following one:

{ } Alt p Diff p f D p T Alt ≤ - ∈ = ) ( ) ( / With ) ( p Diff
refers to the immersion level when the flooding reaches pixel p . Segmentation results are sensitive to this function. In fact, if 0 ) ( = p Diff , (2.4) can be seen as a special case of (2.1). Other extreme case, if

) ( p Diff

reaches its maximum values, all pixels susceptible to be dumped, will be. According to the user requirement, ∑

) ( p Diff

∈ - = ) ( ) ( ) ( ) ( p Neighbor q conn q f p f p Diff (2.4) (2.5)
Obtained results through two implementations of the original and the multi-degree watershed shows that multi-degree immersion method resists the over-segmentation problem effectively. Indeed, the number of detected region in an image brain (181*217*181 voxel volume), decreases from 10991, using the old method, to only 35 using the new method. Computation time and consumed memory size are practically the same. More information about implementation can be found in [START_REF] Shengcai | A Novel Implementation of Watershed Transform Using Multi-Degree Immersion Simulation[END_REF].

Watershed based on path-cost minimization

In this class, there are two possible approaches. The first one associates a pixel to a catchment basin when the topographic distance is strictly minimal to the respective regional minimum. While the second one builds a forest of minimum-path trees, each tree representing a basin. In the following we start by introducing watershed by topographic distance [START_REF] Meyer | Topographic distance and watershed lines[END_REF] before moving to the watershed by image foresting transform [START_REF] Falcão | The Image Foresting Transform: Theory, Algorithms, and Applications[END_REF][START_REF] Lotufo | The Ordered Queue and the Optimality of the Watershed Approaches[END_REF]. Q = all non-marked neighbors 9. end_if

end_while

Based on the drop of water principle, the intuitive idea behind topographic watershed approach in the steepest descent path principal [START_REF] Meyer | Topographic distance and watershed lines[END_REF][START_REF] Najman | Watershed of a continuous function[END_REF]. A drop of water falling on a topographic relief flows down, as "quickly" as possible, until it reaches a regional minimum. Let Ν → D f : be a digital grey value image. Let M S be the set of markers, pixels where the flooding shall start, are chosen. Each is given a different a label Lab . Topographic watershed process can be described by algorithm 3. Let us note that priority level when inserting neighbors (line 2) corresponds to the gray level of the pixel. In line 6, only neighbours that have already been labelled are compared. Finally, only neighbors (in line 8) that

are not yet in the priority queue are pushed into the priority queue. The watershed lines set are the complement of the set of labeled points.

For mathematical formulation of the mentioned process, we follow here the presentation in [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF] which is based on [START_REF] Meyer | Topographic distance and watershed lines[END_REF]. For the sake of simplicity, we restrict our self to the minimal set of notion that will be useful for our propos. We start by introducing the topographic distance. Let us consider

) ( p N G
as the set of neighbors of pixel p , and ) , ( q p d as the distance associated to edge ) , ( q p . Then the lower slope ) ( p LS of f at a pixel p can be defined as follow:

{ }         - = ∪ ∈ ) , ( ) ( ) ( ) ( ) ( q p d q f p f MAX p LS p p N q G (2.6)
And the cost for walking from pixel a p to a neighboring pixel q can be defined as:

( )        = + < > = ) ( ) ( / ) , ( . ) ( ) ( 2 1 ) ( ) ( / ) , ( ). ( ) ( ) ( / ) , ( ). ( ) , cos( q f p f q p d q LS p LS q f p f q p d q LS q f p f q p d p LS q p
The topographical distance between p and q is the minimum of the topographical distances ) , ( q p T f π along all paths between p and q :

[ ] ) , ( ) , ( q p T q p T f q p f MIN π π → ∈ =
We recall that the topographical distance along a general path ) ,..., ( 0

l q p = π is defined as ) , ( cos ) , ( ) , ( 1
1 1 + - + ∑ = i i l i i i f p p t p p d q p T π
Finally we can define the topographic watershed for a grey value image f , with * f the lower completion of f . Each pixel which is not in a minimum has a neighbor of lower grey value (

LC f f = * ). Let I i i m ∈ ) (
be the collection of minima of f . The basin

) ( i m CB of f corresponding to a minimum I i i m ∈ ) (
is defined as a basin of the lower completion of f :

{} { } ) , ( ) ( ) , ( ) ( : / , ) ( * * * * j f j i f i i m p T m f m p T m f i I j D p m CB + < + ∈ ∈ =
And the watershed is the set of points which do not belong to any catchment basin:

c i I i m CB D f Wshed ) ( ( ) ( ∈ ∪ ∩ =
Several shortest paths algorithms for the watershed transform with respect to topographical distance can be found in the literature but the reference algorithm is that of Fernand Meyer. In the following we present variant of Meyer algorithm with integrate the lower slope of the input image as introduced Roerdink and al. [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF]. [ ] ; [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF]. for all (v∈NG(u)) in the future (w.r.t scan order) for u do [START_REF] Lin | Comparison Between Immersion-Based and Toboggan-Based Watershed Image Segmentation[END_REF].

0 ← v lab 5. [ ] ; ∞ ← v dist
if

[ ] [ ] [ ] ( ) u dist v u t u dist < + , cos then 25. [ ] [ ] ( ) v u t u dist v dist , cos + ← 26. [ ] [ ] u lab v lab ← 27.
false stable ←

28.

else if

[ ] WSHED v lab ≠ and [ ] [ ] [ ] ( ) v dist v u t u dist = + , cos then 29. if [ ] [ ] ( ) u lab v lab ≠ then 30. [ ] WSHED v lab = 31.

false stable ←

The second approach to compute a watershed based on path-cost minimization, as we introduced in the beginning, consists on building a forest of minimum-path trees where each tree represent a basin. This approach is described in the framework of image foresting transform [START_REF] Falcão | The Image Foresting Transform: Theory, Algorithms, and Applications[END_REF]. The IFT defines a minimum-cost path forest in a graph, whose nodes are the image pixels and whose arcs are defined by an adjacency relation between pixels. The cost of a path in this graph is determined by a specific path-cost function, which usually depends on local image properties along the path, such as color or gradient. The roots of the forest are drawn from a given set of seed pixels. For suitable path-cost functions, the IFT assigns one minimum-cost path from the seed set to each pixel, in such a way that the union of those paths is an oriented forest, spanning the whole image. The IFT outputs three attributes for each pixel: its predecessor in the optimum path, the cost of that path, and the corresponding root. Returned solution is usually obtained in linear time and requires a variant of the Dijkstra [START_REF] Dijkstra | A Note on Two Problems in Connection with Graphs[END_REF] , Moore [START_REF] Pape | Implementation and Efficiency of Moore Algorithms for the Shortest Root Problem[END_REF] or Dial's shortest-path algorithm [START_REF] Dial | A Computational Analysis of Alternative Algorithms and Labeling Techniques for Finding Shortest Path Trees[END_REF].

For mathematical formulation of the IFT-watershed, we start by defining some basic notions of image foresting transform as introduced in [START_REF] Falcão | The Image Foresting Transform: Theory, Algorithms, and Applications[END_REF]. Actually, an image ImgIn can be seen as a pair ( , )

J I where J refers to a finite set of pixels and I refers to a mapping that assigns to each pixel ( ) p J ∈ , a pixel value ( ) I p in some arbitrary value space. Distinct binary relation between pixels of J will define an adjacency relation A . Once the adjacency has been fixed, ImgIn can be interpreted as a directed graph, whose nodes are image pixels and whose arcs are pixel pairs in A .

Before moving to path cost definition, let's remember that a sequence of pixels The IFT-watershed assumes that seeds pixel correspond to regional minima of the image or to markers that can be considered as imposed minima. The max-arc path-cost function max f is the same as (2.12). We remember that ( ) h t is a fixed but arbitrary handicap cost for any paths (2.12)

(2.13)

(2.14)
starting at pixel t . We remember also that ( , ) w s t is the weight of arc , s t A ∈ , ideally, higher on the object boundaries and lower inside the objects.

There are two usual arc weight functions : a)

( ) q N p ∈ 8. If { } ( ( ) max ( ), ( ) ) C q C p I q > then 9.

{ }

( ) max ( ), ( ) C q C p I q ← ; 10.

Insert q in the FIFO with cost ( ) C q ; 11.

( ) ( ) wshed p wshed q ←

Lotufo and al. [START_REF] Lotufo | IFT-Watershed from Gray-Scale Marker[END_REF] introduce the IFT-watershed from markers, algorithm 6, which can be computed by a single IFT where the labeled markers are root pixels. Proposed algorithm use hierarchical FIFO queue (HFQ).

In this case, the root map can be replaced by the label map which corresponds to the catchment basins and the used path-cost function is given by: { } {

Topological watershed

The original concept behind topological watershed [START_REF] Couprie | Topological grayscale watersheds transform[END_REF] is to define a "topological thinning" that transforms the image while preserving some topological properties, namely the number of connected components of each lower cross-section as we will explain in the following.

Before introducing the topological watershed process, algorithm 7, we define some basic notions. Some of these notions will be resumed in next paragraph for mathematical formulation. Let F be If p 1 is a marker pixel. Otherwise (2.15) International Journal of Image Processing (IJIP), Volume (5): Issue (4), 2011 11 grayscale image and λ be a grey level, the lower cross-section λ F is the set composed of all the points having an altitude strictly lower than λ . A point x is said to be W-destructible for F if its altitude can be lowered by one without changing the number of connected components of λ

F , with ) (x F k =
. A map G is called a W-thinning of F if it may be obtained from F by iteratively selecting a W-destructible point and lowering it by one. A topological watershed of F is a Wthinning of F which contains no W-destructible point. The major feature of this transform is to produce a grayscale image.

The following algorithms give a global description of the computing process. Note that this process is repeated on loop until no W-destructible point remains.

: Topological watershed process

1. For all p in E , check the number of connected components of the lower cross-section at the level of p which are adjacent to p .

Lower the value of p by one if this number is exactly one

For mathematical formulation, we follow description provided in [START_REF] Najman | Watersheds, mosaics and the emergence paradigm[END_REF]. We start by defining a simple point in a graph, in a sense which is adapted to the watershed, and then we extend this notion to weighted graphs through the use of lower sections [START_REF] Couprie | Topological grayscale watersheds transform[END_REF].

Couprie and al. define a transform that acts directly on the grayscale image, by lowering some points in such a manner that the connectivity of each lower cross-section λ F is preserved. The regional minima of the result, which have been spread by this transform, can be interpreted as the catchment basins. The formal definition relies on a the following particular notion of simple point:

Definition 2.3 (Simple point) Let ( ) Γ = , E G be a graph and let E X ⊂ . The point X x ∈ is simple (for X ) if the number of connected components of

{ }

x X ∪ equal to the number of connected components of X . In other words, x is simple (for X ) if x is adjacent to exactly one connected component of X . Now, we can define more formerly destructible point, and the topological watershed:

Definition 2.4 (Topological watershed) Let ( ) E F F ∈ , E x ∈ and ) (x F k = . The point x is destructible (for F ) if x is simple (for k F ).
We say that ( )

E F W ∈
is a topological watershed of F if W may be derived from F by iteratively lowering destructible points by one until stability (that is, until all points of E being nondestructible forW ). Actually checking whether a point is w-destructible, or not, cannot be done locally if the only available information is the graph ( ) , E Γ and the function F and a point may also be lowered several times until it is no more w-destructible. Couprie and al. [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF] propose a new algorithms making possible to perform this test on all the vertices of a weighted graph in linear time, and also to check directly how low the W-destructible point may be lowered until it is no more w- destructible, thanks to the component tree which may be built in quasi-linear time. In the following, we introduce Couprie's functions to identify W-destructible point.

(2.16) (2.

Algorithm 8 : Function W-destructible

Input : Ψ ), ( , F C F ; 1. ← V Set of element of ) (F C pointed by ) (q Ψ for all q in ) ( 1 p - Γ ; 2. If ) ( φ = V then return [ ] φ , ∞ ; 3. [ ] ← m m c k , HighestFORK ) ), ( ( V F C ; 4. If [ ] [ ] φ , , ∞ = m m c k then return ) min(V ; 5. If )) ( ( p F k m ≤ then return [ ] m m c k , else return [ ] φ , ∞ ;
Previous algorithm gives correct results with regard to the definition (2.17) and is linear in time complexity with respect to the number of neighbors of p .

Checking whether a point is w-destructible or not, involves the computation of the highest fork of different elements of the set V(p),see algorithm 9. This may require a number of calls to BLCA (Binary lowest common Ancestor) which is quadratic with respect to the cardinality of V (p): every pair of elements of V (p) has to be considered.

Algorithm 9 : Function HighestFORK

Input : C : a component tree, V : a set of components of C 1. [ ]← 1 1 , c k ) min(V ; // let [k2,c2]...[kn,cn]be the other elements of V 2. 1 k k m ← ; 3. 1 c c m ← ; 4. for i from 2 to n do 5. [ ]← c k, BLCA [ ] [ ]) , , , , (
m m i i c k c k C ; 6. If [ ] [ ] i i c k c k , , ≠ then 1 k k m ← ; 7. 1 c c m ← ; 8. 9. If ) ( 1 k k m = then return [ ] φ , ∞ else return [ ] m m c k , ;
The HighestFork algorithm returns the highest fork of the setV , or the indicator [ ] φ , ∞ if there is no highest fork. This algorithm makes (n-1) calls of the BLCA operator, where n is the number of elements inV .

Let C be a component tree, let V be a set of components of C, we denote by min(V ) an element of V which has the minimal altitude. For following algorithm, we assume that C is represented in a convenient manner for BLCA.

Thus, we must propose a criterion for the selection of the remaining W-destructible points, in order to avoid multiple selections of the same point. Couprie and al. introduce the idea to give the greatest priority to a W-destructible point which may be lowered down to the lowest possible value. They prove that an algorithm which uses this strategy never selects the same point twice. A priority queue could be used to select W-destructible points in the appropriate order. Here, we present their specific linear watershed algorithm which may be used when the grayscale range is small.

Algorithm 10 : Topological watershed [Couprie]

Data : Ψ ), ( , F C F ; Result : F ; 1. For k from min k to ) 1 ( max - k do φ ← k L 2.
For all

) ( E p ∈ do 3. [ ] ) ), ( , , ( , Ψ - ← F C p F le Destructib W c i 4. If ) ( ∞ ≠ i then 5. { }; 1 1 p L L i i ∪ ← - - 6. { } ; 1 - ← i p K 7. { }← p H pointer to [ ] c i, ; 8. end if 9.
end for 10.

For k from min k to ) 1 ( max - k do 11. While ) ( k L p ∈ ∃ do 12. { } p L L k k / = ; 13. If ) ) ( ( k p K = then 14. ; ) ( k p F ← 15.
)

; ( ) ( p H p ← Ψ 16.
For all

) ( ), ( ( q F k p q < Γ ∈ do 17. [ ] ) ), ( , , ( , Ψ - ← F C q F le Destructib W c i ; 18. If ) ( ∞ = i then ∞ ← ) ( p k ; 19. Else if )) 1 ( ) ( ( - ≠ i p k then 20. { } q L L i i ∪ ← - - 1 1 ; 21.

{ }← q H

pointer to [ ] c i, ;

Watershed transform based on local condition

There is a big similarity between this approach and the drop of water one. Actually basin surface increases in a progressive manner. The local condition of label continuity is iteratively applied along the steepest descent path that reaches the basin minimum. The downhill algorithm, the hill climbing algorithm and the toboggan algorithm are based on this approach. More details of the first two algorithms are given by [START_REF] Lin | Comparison Between Immersion-Based and Toboggan-Based Watershed Image Segmentation[END_REF][25] and the toboggan algorithm will be detailed later in this section. Differences between these three algorithms lie in the processing strategy and data structure as shown in [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF].

For mathematical formulation, we follow description provided in [START_REF] Audigier | Watershed by image foresting transform, tie-zone, and theoretical relationship with other watershed definitions[END_REF] As we mentioned earlier, we will introduce the toboggan algorithm [START_REF] Lin | Comparison Between Immersion-Based and Toboggan-Based Watershed Image Segmentation[END_REF][START_REF] Mortensen | Toboggan-based intelligent scissors with a four-parameter edge model[END_REF] as a reference of the local condition watershed approach. Actually this algorithm is referred as a drainage analogy. It seeks to identify the steepest descent from each pixel of the gradient magnitude of the input image to a local minimum of the surface. Then pixels that belong to the same minima are merged by assigning them a unique label. Sets of pixels having the same label will define catchment basins. The resulting watershed regions are divided by a boundary path which will build the watershed lines.

Let consider

: G D R +
→ as a gradient magnitude image, where D is the indexing domain of the image. D can be decomposed into a finite number of disjoint level sets since pixels are sorted in the increasing order. Sets can be denoted by:

{ } | ( ) h D p D G p h = ∈ = .

( ) h G p =

3.

{ } min ( ), ( )

MIN h G q q Neighbor p = ∈ 4.

If ( )

MIN h h > then 5.

{ } ( | ( ) ) & &( ( ))

MIN S q G q h q Neighbor p = = ∈ For mathematical foundations, we will follow the notations in [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF] to present some basic definitions to handle with minimum spanning forest cuts and watershed-cuts.

Let G be graph with ( ( ), ( ))

G V G E G = . ( ) V G is a finite set of vertex of G . Unordered pairs of ( )
V G , called also edge of G , constitute the element of ( ) E G set. Let denote the set of all maps from E to ℤ by F and we consider that any maps of F weights the edges of G . Let F F ⊆ and ( ) u E G ∈ , F( ) u will refers to the altitude of u and M(F) will refers the graph whose vertex set and edge set are, respectively, the union of the vertex sets and edge sets of all minima of F . Let X and Y be two sub-graphs of G . We say that Y is a forest relative to X if Y is an extension of X and for any extension Z X ⊂ of X , we have Y=Z whenever V(Z)=V(Y) (i) Y is said a spanning forest relative to X (for G ) if Y is a forest relative to X and if V(Y)=V . In this case, there exists a unique cut S for Y . It is composed by all edges of G whose extremities are in two distinct components of Y . Since Y is an extension of X, it can be seen that this unique cut S (induced by Y ) is also a cut for X .

(ii) Y is said a minimum spanning forest relative to X (for F, in G) if Y is a spanning forest relative to X and if the weight of Y is less than or equal to the weight of any other spanning forest relative to X . In this case, S is considered as a minimum spanning forest cut for X .

Trough these equivalences, Cousty demonstrate that the set S E ⊆ is a minimum spanning forest cut for M(F) if and only if S is a watershed cut of F , that can be computed by any minimum spanning tree algorithm. And he proposes a linear algorithm to compute it using a new 'stream' notion that we will not detail in this section. Only the stream algorithm will be introduced. Now, before presenting the watershed-cuts algorithm we just recall the definition of the minimal altitude of an edge. Let denote by ( ) F x the map from V to ℤ such that for any x V ∈ , ( ) F x is the minimal altitude of an edge which contains x . Then a path 0 ,..., l x x π =

, is considered as a path of steepest decent for F (in G ) if for any [ ]

1, i l ∈ , { } 1 1 ( , ) ( ) i i i F x x F x - - - = .
In the following, we introduce the watershed-cuts computing algorithm and stream function. [ ] L,Lab Stream(V,E,F,ψ,x) ← L ← L ∪ {z}; // z is the only well of L 11.

L' ← {z}; // switch the in-depth exploration first 12.

breadth_first ← FALSE 13. Else 14.

L ← L ∪ {z}; // therefore z is a well of L 15.

L' ← L' ∪ {z}; // continue exploration in width first 16. return [L,-1]

Classification of watersheds transforms

In this section we will learn from different syntheses present in Roerdink [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF] and Audigier [START_REF] Audigier | Watershed by image foresting transform, tie-zone, and theoretical relationship with other watershed definitions[END_REF] works. The following table summarizes some characteristic of introduced watershed transform. Selected criteria are justified by our objective to identify the most suitable algorithm for parallel implementation.

runs also in linear time with O(n + m) complexity. Lotufo and Cousty's algorithm run also in linear time. Cousty's algorithm is executed at most O(|E|) times.

Conclusion

In this paper, we have presented an intensive study of all existing watershed transform in the discrete case: WT based on flooding, WT based on path-cost minimization, watershed based on topology preservation, WT based on local condition and WT based on minimum spanning forest.

The first major contribution in this paper is the global nature of the proposed study. In fact, for each approach, we start by giving informal definition, then we present processing procedure followed by mathematical foundations and the algorithm of reference. Recent publications based on some approach are also presented and discussed. The contribution concerns algorithms' classification according to criteria of recursion, complexity, basins computing and topology preservation.

In our future work we will present a parallel MSF-Watershed algorithm suitable for shared memory parallel machines based on Jean Cousty stream function already presented.

FIGURE 1 :

 1 FIGURE 1: (a) Cleavage fractures in steel, (b) contour of (a) obtained truth watershed definition introduced by Beusher and al. [2] in 1979, (c) Maximum intensity projection of original human lower limb (d) Bone tissue removed using mask extended with 3D watershed transform introduced by Straka and al.[5] in2003.

  can be even a constant function or a function computed according to the local characteristic of p . In their paper, Shengcai and Lixu define it as introduced in (2.5) with) ( p Neighboris the set of all p neighbors. And conn refers to the predefined connectivity:

Algorithm 3 :

 3 Watershed by topographic distance process 1. Foreach (marked area ∈ M

Algorithm 4 :

 4 Watershed by topographical distance[Meyer] Data : Lower complete image im on a dgital grey scale image G=(D,E) with cost Result : Labelled watershed image lab on D 1. #define WSHED 0 //label of the watershed pixels 2. //Uses distance image dist. On output, div[v]=im[v], for all v ∈ D 3. for all (v∈D) do //Initialize 4.

17

 17

  )

Algorithm 14 :

 14 Watershed-cuts algorithm[Cousty and al.] Data : (V,E,F) : Edge-weighted graphs; Result: Ψ: a flow mapping of F1. Foreach ( V) x ∈ do ψ(x) NO_LABEL ← ∈ such that (ψ(x)=NO_Label) do 4.

Algorithm 15 : 3 .

 153 Stream function [Cousty and al.] Data : (V,E,F) : Edge-weighted graphs; Ψ : a label of V; x : point of V; Result : [L, lab] : L is a flow obtained from x (source of L) ; lab is the associated label to an Θflux included in L or (-1). 1. L ←{x} 2. L'←{x} // the set of sources not yet explored of L While there exists (y ∈ L') do 4.L' ← L' \ {y}; 5.breadth_first ← TRUE ; 6.While (breadth_first) and (∃ {y,z} ∈ E / z ∉L and F({y,z})= F(y)) do 7.If (Ψ(z) # No_label) then 8.return [L,Ψ(z)] // exist an Θflow L already labelled 9.

  Alt min and Alt max the minimum and the maximum value of. The threshold set of at level

	Algorithm 1: Flooding watershed process
	1. for level from Altmin to Altmax
	2.	// Action 1 : Extend existing region
	3.	foreach (R∈ LR) do Growing [R] until level Alt;
	4.	end_for
	5.	// Action 2 : Create new region
	6.	foreach (Pixel P ∈ level)
	7.	if (Pixel P is not associated to any region R) then
	8.	Create new region [R] in LR ;
	9.	Add Pixel P to region [R];
	10.	Growing [R] until level ;
	For mathematical formulation of the mentioned process, let	f :	D	→	Ν	be a digital grey value
	image, with Alt is:

  We introduce the max-arc path-cost function max f that will be used later. Note that ( ) h t and ( , ) w s t are fixed.

	1 2 , ,..., k t t t ( ) org π the f π can be defined since the search to paths start in a given π = 1) < < -constitute a path. In the following we will denote by i k } max + ∈ for (1 1 ( , ) i i t t A ( ) t h t = { , ) max ( ), ( , ) s t f w s t π where f max ( f π ⋅ = For IFT use, a specific function ( ) S set ( ) S J ⊂ of seed pixels. { ( ) ( ( ) ) ( ) S f if org S otherwise f π π π ∈ +∞ Now, we can introduce the spanning forest concept. We remember that a predecessor map is a function P that assigns to each pixel t J ∈ either some other pixel J ∈ or a distinctive marker J ∉ Μ . Thus, a spanning forest (SF) can be seen as a predecessor map which contains no cycles. Definition 2.2 (Spanning forest) For any pixel t J ∈ , a spanning forest P defines a path * ( ) P t recursively as t if ( ) P t = M and * ( ) , P s s t ⋅ if ( ) P t s = ≠ M , we denote by 0 ( ) P t the initial pixel of * ( ) P t . Algorithm 5: IFT Algorithm [Falco and Al.] Input: Img = ( , ) J I : Image, A JxJ ⊂ : Adjacency relation, f :path-cost function Output: P : optimum path forest, , α β : two sets of pixels with J α β ∪ = . 1. Set { } α ← , J β ← //Initialize 2. for all pixels t J ∈ do //Initialize 3. Set ( ) P t ← M 4. end for 5. while β φ ≠ do //Compute 6. remove from β a pixel s such that * ( ( )) f P s is minimum, 7. add s to α 8. for each pixel t such that ( , ) s t A ∈ 9. If * * ( ( ) , ) ( ( )) f P s s t f P t ⋅ < then 10. Set ( ) P t s ← ; Falcan and al. algorithm describes IFT computing. Its algorithm is based on Dijkstra's procedure for computing minimum-cost path from a single source in a graph and returns an optimum-path forest P for a seed-restricted cost function s f or any pixel t with finite cost ( *( )) s f P t . Pixels will origin max ( ) belong to a tree whose root is a seed pixel.

i t of π and by ( ) dst π the destination k t of π . Now let assume given a function f that assigns to each path π a path cost ( ) f π , in some totally ordered setυ of cost values.

  Data : I : input image, wshed : labeled marker image Result : wshed : watershed catchment basins Aux : C : cost map, initialized to infinity; FIFO :hierarchical FIFO queue

		1 ( , ) w s t	=	( ) J s J t ( ) -	, where ( ) J s refers to the intensity of pixel of s . In that case IFT-
		Watershed is said by dissimilarity.
	b)	2 ( , ) w s t G t ( ) =	, where	( ) G t is the morphological gradient of ImgIn at pixel t . In that
		case IFT-Watershed is said on gradient.
	Algorithm 6 : IFT-watershed from markers [Lotufo]
	1. for all pixels	( ) 0 wshed p ≠ do //Initialize
	2.	( ) C p	←	( ) I p
	3.	Insert p in FIFO with cost ( ) C p
	4. end for
	5. while FIFO do //Propagation
	6.	p ← remove from FIFO
	7.	for each

  Lin and al[START_REF] Lin | Comparison Between Immersion-Based and Toboggan-Based Watershed Image Segmentation[END_REF] define the following pixels classes: Class 1 C refers to all pixels p in h D with an altitude strictly greater than the altitude of its lowest neighbor. Class 2Based on this process, authors introduce the following algorithm to compute watershed. Toboggan Algorithm[Lin and al.]

								C refers to all pixels p in h D belonging to a
	connected component with one or more catchment basin and	1 p C ∉ . Finally, class 3 C refers to all
	pixels p in h D belonging to a connected component without any catchment basin. Thus we can
	give a global description of the computing process.
	Algorithm 11 : Toboggan watershed process
	1. Records the sliding directions for all	(	1 p C ∈	) ( ∪ ∈ p C	2	)	in D
		a. Records the lowest neighbours of all	(	1 p C ∈	)	in D .
		b. Grown region from all	(	1 p C ∈	)
	2. Assign label for all	(	3 p C ∈	)
	3. Assign label to each unlabeled image by first tobogganing then backtracking using best first
		search.
	Algorithm 12: Data : Img : a gradient magnitude image;
	Result : L : a label image, Q : empty FIFO queue ;
	1.	For all (	) ∈ do //Simulation of sliding for all C1 pixels p D
	2.

Watershed transform based on minimum spanning forest

  The original idea is very close to the second case of the path cost minimization based watershed that consist on building a spanning forest from a graph as we introduced in section 2.2. Actually, the beginning was with Meyer[START_REF] Meyer | Minimum Spanning Forests for Morphological Segmentation[END_REF] who proposes compute watershed transform from a weighted neighborhood graph whose nodes are the catchment basins corresponding to the minima of the image. Arcs of the graph, that separate neighbor catchment basins, are weighted by the altitude of the pass between these basins. Extracted minimum spanning forests define partitions that are considered solution of watersheds. It's important to mention that returned solutions are multiple. Authors established also the links between the minimum spanning forest and flooding from marker algorithms. Trough Meyer's bases, Cousty and al.[START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF] introduce the watershed-cuts and establish the optimality of this approach by showing the equivalence between the watershed-cuts and the separations induced by minimum spanning forest relative to the minima.

	5.	End for	(2.18)
	6.	If S has a unique label α then
	7.	( ) L p α =
	8.	Else
	9.	( ) L p RIDGE label _ =
	10.	End if
	11. End if
	2.5
	6.	( ) SlidingList p S =

International Journal of Image Processing (IJIP), Volume (5): Issue (4), 2011

( )

For all ( ( )) q Neighbor p ∈ and ( ( )

If ( ( ) )

SlidingList q φ = then 19.

Append ( )

( )

Append ( ) p to ( ) SlidingList q

24.

End if 25.

End while ( ) S SlidingList p =

3.

For all ( ) q S ∈ do

4.

Resolve(q) The starting point is the definition space; we note that IFT-Watershed and MSF-Watershed definitions are limited to the discrete space while the other watersheds definitions are spread into continue space. IFT-Watershed, MSF-watershed and LC-Watershed form the region based watershed transform family since pixels are assigned to basins. Flooding-Watershed, TD-Watershed and Topological-Watershed form the line based watershed family since some pixels are labeled as watershed. Only Topological-Watershed defines lines that consistently separate basins while Flooding-Watershed and TD-Watershed merely swing between thick and disconnected watershed lines. Through definitions, only Flooding-Watershed and TD-Watershed return unique solution while all other definitions return multiple solutions. Note that set of solutions returned by the IFT-Watershed can be unified by creating litigious zones when solutions differ [START_REF] Audigier | Watershed by image foresting transform, tie-zone, and theoretical relationship with other watershed definitions[END_REF].

Watershed based on

All six algorithms, that don't exactly include their definitions, return unique solution but don't preserve the number of connected components of the original input image. Actually, Vincent-Soille, Meyer and Lin's algorithm don't preserve important topological features. Only Lotufo, Couprie and Cousty's algorithm are correct from this point of view.

Regarding computing process, only Flooding-Watershed needs pixel's sorting while others transforms will pass this costly step. But this does not preclude associated algorithms to use hierarchical structures when implementing. Except Cousty's algorithm that doesn't need any hierarchical queue. Vincent-Soille and Meyer's algorithms impose also a prior minima computation, which is not the case of the others. For complexity, observe that Vincent and Soille algorithm runs in linear with respect to the number N of pixels in the image which is processed. In most current situations of image analysis, where the number of possible values for the priority function is limited and the number of neighbors of a point is small constant, Couprie's algorithm