Keywords: PARALLELIZATION STRATEGY, PARALLEL PROCESSING, TOPOLOGY, WATERSHED, SKELETON, SMOOTHING, SHARED MEMORY, THREADS COORDINATION

ABSTRACT

Topological features of an object are fundamental in image processing. In many applications, including medical imaging, it is important to maintain or control the topology of the image. However the design of such transformations that preserve topology and geometric characteristics of the input image is a complex task, especially in the case of parallel processing.

Parallel processing is applied to accelerate computation by sharing the workload among multiple processors. In terms of algorithm design, parallel computing strategies profits from the natural parallelism (called also partial order of algorithms) present in the algorithm which provides two main resources of parallelism: data and functional parallelism. Concerning architectural design, it is essential to link the spectacular evolution of parallel architectures and the parallel processing. In effect, if parallelization strategies become necessary, it is thanks to the considerable improvements in multiprocessing systems and the rise of multi-core processors. All these reasons make multiprocessing very practical. In the case of SMP machines, immediate sharing of data provides more flexibility in designing such strategies and exploiting data and functional parallelism, notably with the evolution of interconnection system between processors.

In this perspective, we propose a new parallelization strategy, called SD&M (Split Distribute and Merge) strategy that cover a large class of topological operators. SD&M has been developed in order to provide a parallel processing for many topological transformations.

Based on this strategy, we proposed a series of parallel topological algorithm (new or adapted version). In the following we present our main contributions:

(i)
A new approach to compute watershed transform based on MSF transform, that is parallel, preserves the topology, does not need prior minima extraction and suited for SMP machines. Proposed algorithm makes use of Jean Cousty streaming approach and it does not require any sorting step, or the use of any hierarchical queue. This contribution came after an intensive study of all existing watershed transform in the discrete case.

(ii)

A similar study on thinning transform was conducted. It concerns sixteen parallel thinning algorithms that preserve topology. In addition to performance criteria, we introduce two qualitative criteria, to compare and classify them. New classification criteria are based on the relationship between the medial axis and the obtained homotopic skeleton. After this classification, we tried to get better results through the proposal of a new adapted version of Couprie's filtered thinning algorithm by applying our strategy.

(iii) An enhanced computation method for topological smoothing through combining parallel computation of Euclidean Distance Transform using Meijster algorithm and parallel Thinning-Thickening processes using the adapted version of Couprie's algorithm already mentioned.

RÉSUMÉ

Les caractéristiques topologiques d'un objet sont fondamentales dans le traitement d'image. Dans plusieurs applications, notamment l'imagerie médicale, il est important de préserver ou de contrôler la topologie de l'image. Cependant la conception de telles transformations qui préservent à la foi la topologie et les caractéristiques géométriques de l'image est une tache complexe, en particulier dans le cas du traitement parallèle.

Le principal objectif du traitement parallèle est d'accélérer le calcul en partagent la charge de travail à réaliser entre plusieurs processeurs. Si on approche cet objectif sous l'angle de la conception algorithmique, les stratégies du calcul parallèle exploite l'ordre partiel des algorithmes, désigné également par le parallélisme naturel qui présent dans l'algorithme et qui fournit deux principales sources de parallélisme : le parallélisme de données et le parallélisme fonctionnelle.

De point de vue conception architectural, il est essentiel de lier l'évolution spectaculaire des architectures parallèles et le traitement parallèle. En effet, si les stratégies de parallèlisation sont devenues nécessaire, c'est grâce à des améliorations considérables dans les systèmes de multitraitement ainsi que la montée des architectures multi-core. Toutes ces raisons font du calcule parallèle une approche très efficace. Dans le cas des machines à mémoire partagé, il existe un autre avantage à savoir le partage immédiat des données qui offre plus de souplesse, notamment avec l'évolution du système d'interconnexion entre processeurs, dans la conception de ces stratégies et l'exploitation du parallélisme de données et le parallélisme fonctionnel.

Dans cette perspective, nous proposons une nouvelle stratégie de parallèlisation, baptisé SD&M (Split, Distribute and Merge) stratégie qui couvrent une large classe d'opérateurs topologiques. SD&M a été développée afin de fournir un traitement parallèle de tout opérateur basée sur la transformation topologique. Basé sur cette stratégie, nous avons proposé une série d'algorithmes topologiques parallèle (nouvelle version ou version adapté). Nos principales contributions sont : (i) Une nouvelle approche pour calculer la ligne de partage des eaux basée sur 'MSF transform'.

L'algorithme proposé est parallèle, préserve la topologie, n'a pas besoin d'extraction préalable de minima et adaptée pour les machines parallèle à mémoire partagée. Il utilise la même approche de calcule de flux proposé par Jean Cousty et il ne nécessite aucune étape de tri, ni l'utilisation d'une file d'attente hiérarchique. Cette contribution a été précédé par une étude intensive des algorithmes de calcule de la ligne de partage des eaux dans le cas discret.

(ii) Une étude similaire sur les algorithmes d'amincissement a été menée. Elle concerne seize algorithmes d'amincissement qui préservent la topologie. En sus des critères de performance, nous somme basé sur deux critères qualitative pour les comparer et les classés. Après cette classification, nous avons essayé d'obtenir de meilleurs résultats grâce avec une version adaptée de l'algorithme d'amincissement proposé par Michel Couprie.

(iii) Une méthode de calcul amélioré pour le lissage topologique grâce à la combinaison du calcul parallèle de la distance euclidienne (en utilisant l'algorithme Meijster) et l'amincissement/ épaississement parallèle (en utilisant la version adaptée de l'algorithme de Couprie déjà mentionné).

ACKNOWLEDGEMENTS / REMERCIEMENTS

En premier lieu, je tiens à remercier particulièrement et à témoigner toute ma reconnaissance à professeur Mohamed AKIL, qui m'a accompagner toutes ces années. Il a su dirigé mes travaux de recherche tout en me laissant suffisamment d'autonomie. Je le remercie pour son beau sens d'écoute et sa patience face à mon atonie lors des moments difficiles de rédactions et de révision des articles. Merci professeur AKIL de m'avoir enseigné la science et m'avoir appris à aimer la recherche.

Je remercie également professeur Jean-Luc GAUDIOT et professeur Olivier DEFORGES qui m'ont fait l'honneur d'être rapporteur de mon manuscrit de thèse. J'adresse mes plus sincères remerciements à professeur Antoine MANZANERAT qui a accepté de consacrer du temps à la relecture de ce manuscrit. Mes profonds remerciements s'adressent également à professeur Gilles BERTRAND qui m'a fait l'honneur d'accepter la présidence du jury.

Je remercie très vivement touts les membres du laboratoire A3SI, chercheurs, techniciens, et étudiants en thèse pour leur accueil chaleureux. Un grand merci à Thierry, Laurent, Eric, Christophe, Eva, Christine, Jean-Claude, Michel, Denis, Rostom, Jean, Dalila, John, Imran, Oussama… Sachez que je garderai un très bon souvenir de nos discussions au tour du café, du repas, ou du barbecue … Mes remerciements s'adressent plus spécialement à ceux, invisibles, qui m'ont soutenu et cru en moi : Merci à mes parents, mon frère et mes soeurs pour leur amour inconditionnel, leur soutien et tout ce qu'ils m'ont apporté mais dont ils n'ont qu'une toute petite idée. Merci à celle qui m'a enseigné que l'amour ne s'apprend pas dans les livres et qui m'a offert le plus beau cadeau du monde, mes deux filles Nersrine et Lina : Un merci très spécial à ma femme Marina qui m'a accompagné durant la préparation de cette thèse. Je lui suis très reconnaissant pour sa patience et ses encouragements dans mes moments de doute. maging applications including medical imaging (MRI 2 scan, PET 3 scan, SPECT 4 scan) 2D/3D (even 4D in some cases) implement different steps: acquisition (raw format), preprocessing, processing, analysis, interpretation and display of results. The growing needs for such applications, variety and complexity of their algorithms, new requirements in terms of performance and quality, require the development of new methods, supports and packaged software implementation of these application on new multi-processor based-architecture 5 that offer more computing power. Such software design flow must include all processing steps from high level specification (algorithm specification) to the low level specification (optimized implementation based on parallel code distributed between different processors of the platform) with respect of real-time constraints such that throughput and latency.

LIST OF TABLES

LIST OF ALGORITHMS

ALGORITHM : SCHEDULING POLICY [MAHMOUDI AND AKIL

] ---------------------------------------------------------------------31 ALGORITHM : MERGING ALGORITHM [MAHMOUDI AND AKIL] ------------------------------------------------------------------35 ALGORITHM : FLOODING WATERSHED PROCESS ------------------------------------------------------------------------------------39 ALGORITHM : FLOODING WATERSHED [VINCENT & SOILLE] ---------------------------------------------------------------------40 ALGORITHM : WATERSHED BY TOPOGRAPHIC DISTANCE PROCESS ---------------------------------------------------------------42 ALGORITHM : WATERSHED BY TOPOGRAPHIC DISTANCE [MEYER] --------------------------------------------------------------44 ALGORITHM : IFT ALGORITHM [FALCO AND AL.] ---------------------------------------------------------------------------------45 ALGORITHM : IFT -WATERSHED FROM MARKERS [LOTUFO] --------------------------------------------------------------------46 ALGORITHM : TOPOLOGICAL WATERSHED PROCESS -------------------------------------------------------------------------------47 ALGORITHM : W-DESTRUCTIBLE FUNCTION [COUPRIE AND AL.] ----------------------------------------------------------------48 ALGORITHM : HIGHESTFORK FUNCTION [COUPRIE AND AL.] -------------------------------------------------------------------48 ALGORITHM : TOPOLOGICAL WATERSHED [COUPRIE] ----------------------------------------------------------------------------49 ALGORITHM : TOBOGGAN WATERSHED PROCESS ---------------------------------------------------------------------------------51 ALGORITHM : TOBOGGAN ALGORITHM [LIN AND AL.] ---------------------------------------------------------------------------51 ALGORITHM : RESOLVE FUNCTION [LIN AND AL.] --------------------------------------------------------------------------------52 ALGORITHM : WATERSHED-CUTS ALGORITHM [COUSTY AND AL.] -------------------------------------------------------------53 ALGORITHM : STREAM FUNCTION [COUSTY AND AL.] ----------------------------------------------------------------------------54 ALGORITHM : PARALLEL WATERSHED-CUT [MAHMOUDI AND AKIL.] ----------------------------------------------------------58 ALGORITHM : FUNCTION S-LABELING [MAHMOUDI AND AKIL] -----------------------------------------------------------------59 ALGORITHM : FUNCTION F-LABELING [MAHMOUDI AND AKIL] -----------------------------------------------------------------60 ALGORITHM : THINNING ALGORITHM -SECOND VERSION [GUO AND HALL]--------------------------------------------------75 ALGORITHM : THINNING ALGORITHM [ECKHARDT AND MADERLECHNER] ----------------------------------------------------76 ALGORITHM : THINNING ALGORITHM [COUPRIE, BEZERRA AND BERTRAND ] ------------------------------------------------ 76
Discrete topology offers a range of essential tools in image processing thus applications in medical imaging calls several algorithms based on topological transformation (these algorithms have the specificity to process over input image while preserving its topology which allows keeping some important information intact). But the increasing size of processed data, due to the improved of capture devices resolution, and constraints in terms of processing time, make the development of such standard application very complex. Indeed, computing power required for these applications currently uses parallel architectures as new computer machines (for reasons of cost and availability). Only parallel processing provides a cost effective solution for this required power by increasing the number of processors and adding an efficient communication system between them which makes coding more complicated: exploring, in optimal way, intra/inter processors parallelism and work distribution among different threads. The availability, on SMP 6 machines, of multi-core processor / RISC-based cores (where the workload can be shared between different processors) makes it necessary to study and develop appropriate and effective parallelization strategies of such image processing algorithms for this type of machine. However, there is no parallelization strategy common to a set of algorithms based on topological operators to efficiently implement these algorithms on parallel machines. This strategy unifies the optimized implementation of parallel algorithms on these specific architectures via parallelization, work-load distribution and effective management of memory hierarchy. We therefore propose to study and formalize such parallelization strategy and define suitability metrics to assess performance on this type of architecture.

2 Magnetic Resonance Imaging 3 Positron Emission Tomography 4 Single Photon Emission Computed Tomography 5 Architecture interconnecting many processors: parallel processors (multi-core processor), parallel specific processors (Graphic Processor Unit). 6 Shared Memory Parallel machines I

Contributions

Research presented in this manuscript has been done in the Gaspard-Monge computer science research laboratory (LIGM) of Paris-Est University, ESIEE A3SI team, CNRS-UMLV-ESIEE (UMR 8046).

We frame our work in the field of algorithms based on topological transformation in order to study their parallelization on shared memory parallel machines (SMP machines).

Topological features of an object are fundamental in image processing. In many applications, including medical imaging, it is important to maintain or control the topology of the image. However the design of such transformations that preserve topology and geometric characteristics of the input image is a complex task, especially in the case of parallel processing. Here, the main goal of parallel processing is to accelerate computation by sharing the workload among multiple processors. In terms of algorithm design, parallel computing strategies profits from the partial order of algorithms, called also the natural parallelism present in the algorithm which provides two main resources of parallelism: data and functional parallelism. Now, from a viewpoint architectural design, it is essential to link the spectacular evolution of parallel architectures and the parallel processing. In effect, if parallelization strategies become necessary, it is thanks to the considerable improvements in multiprocessing systems and the rise of multi-core processors. And during the last decade, clock speed of processors in multi-core architectures has increased by almost two and associated cache size has increased tenfold with the addition of a third cache level L3 which ensures optimal L2 access speed while increasing the total cache. All these reasons make multiprocessing very practical. In the case of SMP machines, it adds another advantage that is the immediate sharing of data which provides more flexibility, notably with the evolution of interconnection system between processors, in designing such strategies and exploiting data and functional parallelism.

In this perspective, we propose a new parallelization strategy, called SD&M (Split Distribute and Merge) strategy that cover a large class of topological operators. SD&M has been developed in order to provide a parallel processing of any operator based on topological transformation. In practice the most effective parallel algorithm design might make use of multiple algorithm structures thus proposed strategy is a combination of the divide and conquer patterns and eventbased coordination patterns hence the name that we have assigned. Not to be confused with the mixed-parallelism approach (combining data-parallelism and task-parallelism), it is important to mention that proposed strategy (1) 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. (2) It covers only the case of recursive algorithms, while mixedparallelization strategy is effective only in the linear case. (3) It is especially designed for shared memory architecture with uniform access.

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

Based on this strategy, we proposed a series of parallel topological algorithm (new or adapted version). In the following we present our main contributions:

A new approach to compute watershed transform based on MSF 7 transform, that is parallel, preserves the topology, does not need prior minima extraction and suited for SMP machines. Proposed algorithm makes use of Jean Cousty streaming approach [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF] and it does not require any sorting step, or the use of any hierarchical queue. This contribution came after 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. This study can be seen as an update of Roerdink research [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF]. Actually, this study presents detailed description of each watershed approach, associated processing procedure followed by mathematical foundations and the algorithm of reference. Recent publications based on some approaches are also presented and discussed. Our study concludes with a classification of different algorithms studied according to solution uniqueness, topology preservation, prerequisites minima computing and linearity.

A similar study on thinning transform was conducted. It concerns five parallel thinning algorithms that preserve topology: Bernard and Manzanera [START_REF] Bernard | Improved low complexity fully parallel thinning algorithm[END_REF], Jang and Chin [START_REF] Jang | Reconstructable parallel thinning[END_REF], Eckhardt and Maderlechner [START_REF] Eckhardt | Invariant thinning[END_REF], Guo and Hall [START_REF] Guo | Fast fully parallel thinning algorithms[END_REF], and Hall [START_REF] Hall | Fast parallel thinning algorithms: Parallel speed and connectivity preservation[END_REF]. Based on the relationship between the medial axis and the obtained homotopic skeleton, we introduce two classification criteria to compare and classify them. After this classification, we tried to get better results through the proposal of a new adapted version of Couprie's filtered thinning algorithm [START_REF] Couprie | Topological operators for grayscale image processing[END_REF] by applying our strategy.

An enhanced computation method for topological smoothing through combining parallel computation of Euclidean Distance Transform using Meijster algorithm [START_REF] Meijster | A general algorithm for computing distance transforms in linear time[END_REF] and parallel Thinning-Thickening processes using the adapted version of Couprie's algorithm already mentioned.

Report organization

This thesis initially aimed to study how to parallelize topological operators on SMP machines. All along the chapters of this work, parallel version (adapted and new) of fundamental algorithms of this class will be introduced. Beyond new parallel algorithms, we will also present obtained results in terms of computation time, efficiency, cache consumption.

The second chapter will be about parallelization strategy. The global goal of the project surrounding this work was to design an efficient strategy to parallelize algorithm based on topological transform. To do so, we first start by identifying real needs for such strategy, then we define fundamental basis for any parallelization moving from finding concurrency to performance metrics of parallel programs. Then, we briefly introduce our strategy called SD&M 8 with the aim of classifying it among all existing strategy. Finally, we will conclude by presenting all details of SD&M conception: splitting, merging and merging phases.

The third chapter is about watershed transform. Our work starts with a comparative study between five different approaches: we present a review of several definitions of the watershed transform and the associated sequential algorithms, emphasizing the distinction between definition and algorithm specification. The study concludes with a classification of different algorithms according to criteria of recursion, complexity, basins computing and topology preservation. Since we identify the most suited approach to compute parallel watershed, we propose a new algorithm that is parallel, preserves the topology of the input image, does not need prior minima extraction and suited for SMP machines. Contrarily to previously published algorithms, proposed algorithm do not require any sorting step, or the use of any hierarchical queue.

The fourth chapter focuses on thinning algorithms in the framework of critical kernel. We start by resuming Couprie's study on verification methods for the topological soundness of thinning algorithms. Based on this first evaluation of 2D topological thinning algorithm, we propose to go further in this study through new quantitative and qualitative criteria. Since we have establish a new classification of thinning operators and being convinced that best performance can be reached, we propose an adapted version of Couprie's thinning algorithm that we call parallel λ -Skeleton algorithm.

In the fifth chapter, we move to another aspect of topological processing: smoothing filter. We present 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.

In the sixth chapter, we present a review of the research in the form of a critical summary of presented work restating contributions of the thesis. Future work is also presented, summarizing the next steps to follow into the research of parallelization strategy.

PARALLELIZATION STRATEGY

ultiprocessor chips make computing more efficient by exploiting parallelism which is one of the outstanding challenges of modern computer sciences. Parallel processing can be defined as the division of a fixed amount of work among different processors that run simultaneously with a common objective. The primary advantage of such processing comes exactly out of its ability to handle a large volume of tasks (or data) with reasonable latency and cadency. Obtained result should be faster with better efficiency, scalability and portability in comparison with single processor implementation. Improving these metrics doesn't depend only on the parallel programming design. The scope of algorithm design space should be expand by including finding concurrency design space, architecture design space, parallel implementation mechanisms and performance metrics of parallel programs. These steps, listed in chronological order, sets out a strategy for parallelization. However, it is generally recognized that designing such strategy is not an easy task [START_REF] Carriero | How to Write Parallel Programs. A Guide to the Perplexed[END_REF][START_REF] Chandy | An Introduction to Parallel Programming[END_REF][START_REF] Darlington | Building Parallel Applications without Programming[END_REF]. In fact, such design is significantly more complex than sequential programs design on single processor computers. As cited in [START_REF] Arjona | Architectural Patterns for Parallel Programming, Models for Performance Estimation[END_REF], some example of parallel software design [START_REF] Andrews | Concurrent Programming: Principles and Practice[END_REF][START_REF] Foster | Designing and Building Parallel Programs, Concepts and Tools for Parallel Software Engineering[END_REF][START_REF] Culler | Parallel Computer Architecture. A Hardware/Software Approach[END_REF] illustrate encountered difficulty when scientific code has been hand-crafted for specific machines and problems, at immense expense.

The study of parallelism, or how parallel process could be expressed in programming terms, start since the sixties. Dijkstra [START_REF] Dijkstra | Co-operating Sequential Processes[END_REF] was among the first to develop an initial proposal for the treatment of parallelism in a programming language by adding new concepts to sequential programming, in order to extend it into a concurrent programming. He introduced new elements such as mutual exclusion, event synchronization and critical section. The most important concept that was introduced by Dijkstra is semaphores. He didn't bring only solutions, He also raise new issues related to parallel processing such as deadlock. Ten years later, things become more formal and definitions more explicit with Hoare's work [START_REF] Hoare | Communicating Sequential Processes[END_REF]. Actually, he defines new language for the formal specification of parallel algorithms, known as Communicating Sequential Processes. His work starts by analyzing general basic structures used in programming such as assignment, sequence, repetition, and selection. Then, he introduced new structures for expressing parallelism, communication, and control of non-determinism between processes within a multiprocessor architecture. More details about both approaches can be found in [START_REF] Arjona | Architectural Patterns for Parallel Programming, Models for Performance Estimation[END_REF]. These different approaches [START_REF] Dijkstra | Co-operating Sequential Processes[END_REF][START_REF] Hoare | Communicating Sequential Processes[END_REF][START_REF] Brinch-Hansen | Distributed Processes: A Concurrent Programming Concept[END_REF] intersect at a common point: they emphasize the close link between strategy and the nature of the application to parallelize. We therefore propose to study and formalize parallelization strategy of image processing operators (topological operators) and define suitability metrics to assess performance on SMP architecture. This chapter is organized as follow: we begin, in section 2.1, by highlighting the real need for a common parallelization strategy of topological operators. After introducing the basic foundation for any successful parallelization (section 2.2), we will focus on Split Distribute and Merge (SD&M) strategy that we propose by an initial classification over all existing strategies (section 2.3), followed by detailed description of SD&M conception in section 2.4.

M

Lack of common parallelization strategy for topological operators

In 1996, Bertrand [START_REF] Bertrand | Topological approach to image segmentation[END_REF] introduced connectivity numbers for grayscale image. These numbers describe locally (in a neighborhood of 3*3) the topology of a point. According to this description any point can be characterized following its topological characteristics. He also introduced some elementary operations able to modify gray level of a point without modifying image topology. These elementary operations of point characterization present the fundamental link of large class of topological operators including, mainly, skeletonization and crest restoring algorithms [START_REF] Couprie | Topological operators for grayscale 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 2D and 3D object smoothing algorithm [START_REF] Couprie | Topology preserving alternating sequential filter for smoothing 2D and 3D objects[END_REF] and topological watershed algorithm [START_REF] Bertrand | On Topological Watersheds[END_REF]. All mentioned algorithms get also many algorithmic structure similarities. In fact associated characterizations procedures evolve until stability with induce common recursively between different algorithms. Also the grey level of any point can be lowered or enhanced more than once. Finally, all the mentioned algorithms get a pixel's array as input and output data structure. Expect in special cases where graphs are used. 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 Sienstra [START_REF] Seinstra | A software architecture for user transparent parallel image processing on MIMD computers[END_REF] for local operators and point to point operators can also be cited as example. For global operators, an adapted parallelization strategy is given in Meijster work [START_REF] Meijster | A general algorithm for computing distance transforms in linear time[END_REF]. Hence the need of a common parallelization strategy for topological operators that offers adapted algorithm structure design space. Chosen algorithm structure patterns to be used in the design must be suitable for SMP machines.

Fundamental basis for parallelization

Before defining parallelization's stages of any sequential problem, it is essential to link the spectacular evolution of parallel architectures and the parallel processing. In reality, if the parallelization strategies are so valuable, it is thanks to substantial improvements in multiprocessing systems and the rise of multi-core processors. In terms of feasibility, it will be easier to design architecture with a single fast processor (clock speed over 3 GHz) than one with many slow processors (clock speed around 1.5 GHz) with the same throughput. But during last years the clock speed of processors in multi-core architectures has increased ,see (tab. 1), by almost two and associated cache size has increased tenfold with the addition of a third cache level L3 which ensures optimal L2 access speed while increasing the total cache. These twin barriers have flipped the equation, making multiprocessing very practical and advised even for small applications. Generally five steps, see (fig. 1), are necessary to move from sequential algorithm running on single core architecture to parallel algorithm that runs with better performance on a multi-core architecture. G. Mattson and al., see [START_REF] Mattson | Patterns for parallel programming[END_REF], present the first four steps for parallel programming: Finding concurrency, algorithm structure, support structure and implementation mechanisms.

Based on Mattson logic, we define all five steps. We should like to stress the importance that we attach to the second and third phase which represent basis of our strategy as we will show later. To remain in the conceptual framework of this section, we provide only the two original organizations: by tasks and by data [START_REF] Feldmann | Dynamic scheduling on parallel machines[END_REF].

Definition Definition Definition Definition ( ( ( (2 2 2 2) ) ) ): : : : Algorithm Design Space

The set of all possible algorithm designs and algorithm design parameters that represent how the extracted concurrency can be mapped onto elementary preprocessors.

In several cases, the problem can be decomposed to a finite set of tasks. Tasks can be grouped according to several criteria: nature of the operation to achieve required operands, action-zone or returned result then groups of task can be defined. The way that the tasks within their group interact is the major feature of the concurrency. If the final solution is obtained after a single execution of all tasks and tasks dependency is null or quasi-null (temporary access to shared variables or messages exchange for synchronization), we can define the parallel task design. If processing is recursive, the problem can be solved by recursively dividing it into sub-problems, solving each sub-problem independently, and then recombining the sub-solutions into a solution to the original problem. This is the well know pattern of divide and conquer. It's important to note that the application of this principal cannot be independent from the type of the algorithm [START_REF] Wangqing | A New Divide and Conquer Algorithm for Graph-based Image and Video Segmentation[END_REF]. In other cases, global processing comes down to a continuous updating of a data structure. Thus it is better to think in terms of organizing data. G. Mattson goes further in this classification. He distinguishes between two particular cases: if the organization focuses on the distribution of data between elementary processors, then it's a simple data decomposition pattern. However, if the organization is the distribution of data between tasks groups: it is a data flow decomposition pattern. More details will be given about this pattern in SD&M strategy classification.

Definition Definition Definition Definition ( ( ( (3 3 3 3) ) ) ): : : : Architecture Design Space

It describes the set of platform that support the extension of parallel programming. Information about how instructions are executed and how memory is managed are presented in this design.

Before moving to coding, it is important to find the most appropriate architecture to support the parallel algorithm using parallel architecture design space (Def.3). This design presents standard classification of parallel computer systems [START_REF] Quinn | Parallel Computing: theory and practice[END_REF]. According to Flynn classification [START_REF] Flynn | Some Computer Organizations and Their Effectiveness[END_REF], there are four types: SISD 9 , SIMD 10 , MISD 11 and MIMD. The most significant structure encountered in the parallel application [START_REF] Steen | Overview of recent supercomputers[END_REF] is MIMD (Multiple Instruction, Multiple Data). In a MIMD machine the processors can execute different operations using their own data. Parallel processing via the application of MIMD machines offers the promise of high performance, and experience with parallel processing is accumulating rapidly. In [START_REF] Buzbee | Applications of MIMD machines[END_REF], Buzbee show, through different examples, that rapid progress is being made in the application of MIMD machines and that parallel processing can yield high performance. We distinguish between two types of MIMD computers: Shared Memory MIMD machines and Distributed Memory MIMD machines. In the case of distributed memory machines, each processor has its own memory but this does not prevent its access to the memories of other processors if necessary. Thread's synchronization and communication must also be guaranteed.

In contrast, in shared memory parallel machines, all processors share the same memory.

Although the cost of inter-processor or memory communication can be high, SMPM design still very efficient. In fact, this cost can be reduced by using the right mechanisms for parallel programming (Def.4). These mechanisms allow better exploitation of target architecture through the use of threads. The most used tools within this framework are: MPI [START_REF] Foster | Designing and Building Parallel Programs[END_REF], OpenMP [START_REF] Chandra | Parallel Programming in OpenMP[END_REF] and TBB [START_REF] Reinders | Intel Threading Building Blocks[END_REF]. After coding and running programs, it's important to evaluate efficiency, scalability and portability of the code by using performance metrics for parallel programs (Def.5). These concepts will list in detail in the last part of this chapter.

Definition Definition Definition Definition ( ( ( (5 5 5 5) ) ) ): : : : Performance metrics of parallel programs

They are a set of measurements that quantify the parallel code such as efficiency, scalability and portability.

Classification of SD&M strategy

As mentioned in last section, chosen algorithm structure patterns to be used in the design must be suitable for SMP machines. In fact, shared memory parallel machines allow access from any processor to any location through shared memory using an interconnection network (processorprocessor and memory-processor). Flexibility of these architectures does not lie in such interconnection network (which is usually predefined by manufacturer), but on shared memory. Actually, programmer perceives a shared, central, and continuous memory. Each memory location or address in unique and identical for any processor of the system. Thus, he take profits of the immediate sharing of data with is very helpful in the conception of any parallelization strategy: The global address space provides a user-friendly programming perspective to memory and data sharing between tasks is both fast and uniform due to the proximity of memory to CPUs. Communication between CPUs can be assured trough shared variables (reading/writing). Network selects appropriate memory block when reading (writing) process from a specific memory address is lunched. Data integrity is guaranteed by synchronization, and coordination mechanizes such as semaphores (Dijkstra), and monitors (Hoare) already introduced in introduction. SMP machines are also non-dedicated architecture using standard component (processor, memory...) so economically reliable. They also offer some flexibility of use in many application areas, particular image processing.

In practice the most effective parallel algorithm design might make use of multiple algorithm structures thus proposed strategy is a combination of the divide and conquer patterns and eventbased coordination patterns, see (fig. 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 [START_REF] Ramaswamy | A framework for exploiting task and data parallelism on distributed memory multicomputers[END_REF]), it is important to mention that our strategy represents the last stitch in the decomposition chain of algorithm design patterns and it provides a finegrained description of topological operators parallelization while mixed-parallelism strategy provides a coarse-grained description without specifying target algorithm. It covers only the case of recursive algorithms, while mixed-parallelization strategy is effective only in the linear case.

It is especially designed for shared memory architecture with uniform access. (Identical processors, equal access times to memory …).

SD&M strategy conception

A parallelization strategy did not aim to optimize a single metric such as speedup. Other than improved performance in terms of execution time, a good strategy has to provide a balance between efficiency, scalability, and portability to dissolve all conflict that exists between these three forces. These notions will be defined later in this section. Actually, any strategy is facing two major barriers. First, the conflict between efficiency and portability: making a program efficient almost requires that the code take into account the characteristics of the specific system on which it is intended to run, which limits portability. A design that makes use of special features of a particular programming environment (as multi-thread environment) may lead to an efficient program for that particular environment, but unusable for a different platform. Second, the conflict between scalability and portability: Improving the scalability is based on a good distribution of work over a finite number of processors for a better exploitation of the N processors' potential. This distribution limits the portability of the program since the number of processor is increased. The relative importance of these diverse metrics will vary according to the nature of the problem at hand. In our case we are dealing with a class of topological operators with common feature, as we shown in section 2.1. Shared memory parallel architectures turned out to be best suited for our needs (section 2.

Definition Definition Definition Definition ( ( ( (7 7 7 7) ) ) ): : : : Event-Based coordination

It is used when dealing with irregular, dynamic or unpredictable data flow.

The splitting phase

The Divide and Conquer pattern is applied first by recursively breaking down a 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 memory-processor), and passed between different stages of a computation.

Definition Definition Definition Definition ( ( ( (8 8 8 8) ) ) ): : : : Scalability

It is a property which exhibits performance linearly proportional to the number of processors employed.

This first stage of division will primarily affect the rate of scalability (Def. 8) of our program. To mount it, we propose the following formalization. Since speedup is the most commonly used metrics for parallel programming, it seems to be a nature choice to begin. So we assume that every program is made up of two parts, sequential and parallel, to establish the following definitions (with 2 ≥ N ):

s t
Processing time of the serial part of a program using one processor.

( )

1 p t
Processing time of the parallel part of a program using one processor.

) (N t p

Processing time of the parallel part of a program using N CPUs.

)

+ =

Total processing time of the serial and parallel part of the program using ) (

N CPUs. ) 1 ( ) 1 (
) ( )) ( 1 ( ) ( N t t N t N N p s p + = - = α β
Scaled percentage of the parallel part of the program using N CPUs

Table 2 : Basic definitions of processing time

Now we can formalize the fixed-size speedup, which fixes the problem size and emphasizes how fast a problem can be solved. By first theoretical approach speedup can be seen as the ratio of a quantity of works by a period of time:

Definition Definition Definition Definition ( ( ( (9 9 9 9) ) )

): : :

: Speedup-(1 st equation) Speedup = Time Work = T W
A second formal definition can be given by applying Amdahl's law [START_REF] Amdahl | Validity of the Single Processor Approach to Achieving Large Scale Computing Capabilities[END_REF] so the speedup can be defined by the ratio of total processing time of the serial and parallel part of the program using one processor by the total processing time of the same parts using N processors.

Definition Definition Definition Definition ( ( ( (10 10 10 10) ) )

): : :

: Speedup-(Amdahl approach) Speedup= ) ( ) 1 ( N t t T T =             + + N t t t t p s p s ) 1 ( ) 1 (
This formula can be written differently using non-scaled percentage (1)

β
Speedup=             - + N ) 1 ( 1 ) 1 ( 1 β β
An alternative formulation referred to as Gustafson's law [START_REF] Gustafson | Reevaluating Amdahl's Law[END_REF] exists. This formulation calibrates the serial percentage according to the total parallel time using N processors.

)

) ( 1 ( ) ( ) ( N N N t T α α - + = )) ( 1 ( ) ( ) 1 ( N N N t T α α - + =
Thus we can define the speedup as follows:

Definition Definition Definition Definition ( ( ( (12 12 12 12) ) )

): : : W . We recall that the efficiency is considered as the ratio of speedup by the number of processor (More details about efficiency will be given in the second section). Ideally, an algorithm should be effective on wide range of numbers of processing elements, from two up to decade. So if N η and ' N η represent the optimal efficiency rate, we can draw the following equation using (Def. 9):

: Speedup-(3 d equation) Speedup=       ) ( ) 1 ( N t t T T = [ ] )) ( 1 ( ) ( N N N α α - + = [ ] ) ( ) 1 ( N N N
If ' N N η η = ' ) ' ( ) ( N N Speedup N N Speedup = ⇔ ' ) ' ( ) ( N N Speedup N N Speedup = ⇔ ) ' ( '* ) ( * ' N t N W N t N W T N T N = ⇔ ) ' , ( * ' * ) ' ( ) ( ' N N N W N W N t N t N N T T ψ = = ⇔
Thus it follows that the only parameter that provides a linear performance in proportion to the number of processors (Def. 8) is the ratio

      ' N N W W .
Hence the importance of splitting step.

Unfortunately such impact can't be shown by applying simply Gustafson approach (Def. 12). Scalability will be express only in term of the number of processor.

If ' N N η η = ' ) ' ( ) 1 ' ( ' ) ( ) 1 ( N N N N N N N N α α - - = - - ⇔ [ ] [ ] ) ' ( ) 1 ' ( ' ) ( ) 1 ( ' N N N N N N N N α α - - = - - ⇔ [ ] [ ] ) ' ( ) 1 ' ( ' * ) ( ) 1 ( ' '* N N N N N N N N N N α α - - = - - ⇔ [ ] [ ] ) ' ( ) 1 ' ( ) ( ) 1 ( ' N N N N N N α α - = - ⇔         + - =         + - ⇔ ) ' ( ) 1 ' ( ) ( ) 1 ( ' N t t t N N N t t t N N p s s p s s         + - =         + - ⇔ ) ' ( * ) 1 ' ( * ) ( * ) 1 '*( N t t t N N N t t t N N p s s p s s         + - =         + - ⇔ ) ' ( * ) ' * ( ) ( * ) ' '* ( N t t t N N N N t t t N N N p s s p s s         + - =         + - ⇔ ) ' ( ) 1 ' ( * ) ( ) 1 '*( N t t N N N t t N N p s p s       - =       - ⇔ ) ' ( ) 1 ' ( * ) ( ) 1 '*( N t N N N t N N T T ) ' , ( ) 1 ' ( * ) 1 '*( ) ' ( ) ( N N N N N N N t N t T T ψ ==       - - =       ⇔

Distribution phase

We attach great importance to work distribution because it 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 synchronization 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.

The main challenge when performing parallel operations on simple point characterization is the dynamic nature of work distribution. Since the workload is not known a priori, assigning work units to different cores in advance is impossible. In the literature, there are two main approaches for multi-core work distribution: the first one, called work queues approach, consist on using a shared work queue in main memory and control access to it via synchronization primitive. The second approach is work stealing. In this case, every core has a separate work queue which is still accessible to other processors. Cores can steal work units from others' queues whenever their own queue is empty. However, all these techniques do not currently work well on new architectures as Xeon for many reasons. Primarily, work-stealing has also been known to be cache-unfriendly for some applications due to randomized stealing [START_REF] Acar | The data locality of work stealing[END_REF]. For tasks that share the same memory footprints, randomized locality oblivious work-stealing schedulers do nothing to ensure scheduling of these tasks on workers that share a cache. This significantly limits, not only scalability, but also efficiency (Def. 13) for some memory-bandwidth bounded applications on machines that have separate caches.

Definition Definition Definition Definition ( ( ( (13 13 13 13) ) ) ): : : : Efficiency

It is the cost of what is actually produced or performed with what can be achieved with the same consumption of resources (processor frequency, memory size, surface, etc.). It is an important factor in determination of productivity.

Efficiency = N η =       N Speedup =       - - N N N N ) ( ) 1 ( α according to (Def. 11) =       T Nt W according to (Def. 9)
It is also important to mention that using memory fence operations, consistency can be enforced, but with relatively high overhead. Even if memory consistency were not a problem, busy waiting such as by spinning on a lock variable is relatively inefficient on an architecture with high memory latency and hardware multi-threaded execution can also lead to priority inversion and prevent other threads on the same core from performing useful work.

The main idea of our approach is to imagine an auto-supplying task system, see (fig. 4). Keeping the local queues filled will be our major goal. Local threads should never have to broadcast over processors because of an empty queue. They should always find something in their local queues due to an auto-supplying task system that allows an automatic check of different queues then permanent redistribution of tasks to maintain a certain balance between all processors. Let's consider a distributed job list where no jobs are duplicated anywhere so each processor local list is unique and exclusive and jobs can be moved between processor only before they go into "executing" status. In spite of each processor balances itself by requesting or stealing work units from others' queues whenever its own queue is empty, we design an auto-supplying system using a shared work queue located in principal memory that supply the processor shortage of work. System maintains a minimum level of work in its queue by importing extra work from others' processors queues. This load balancing feature will keeps all processors busy. Before detailing our approach, we 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] Natarajan | Imprecise and Approximate Computation[END_REF][START_REF] Tilborg | Foundations of Real-Time Computing: Scheduling and Resources Management[END_REF][START_REF] Leung | REAL-TIME SCHEDULING ANALYSIS[END_REF]. According to [START_REF] Leung | REAL-TIME SCHEDULING ANALYSIS[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 (fig. 5. 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 (fig. 5. 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 an adapted scheduling algorithm, simpler to implement and more adapted to topological algorithm implementation.

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

{ } 1 2 , ,..., k T t t t = is the set of all tasks, { } 1 2
, ,..., In this scheduler, each processor will treat at maximum max / max( )

T i T t t t = is the set of tasks to process with T T T ⊂ , { } 1
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 there is 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 [Mahmoudi and Akil]

Aux: T : Set of all tasks, P : Set of all processors 1. While ( )

T ≠ ∅ repeat : 2. T N = Nbr_active_tasks() ;
3.

P N = Nbr_ available_processors();

4.

If ( 0) P N ≠ then 5. If ( T N < P N ) then 6.
For each processor 'Basic-NPS' policy is described by (Alg. 1). 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. This approach, and despite the centralized aspect, does not depend on the number of processors (number of processor is not predefined) or the minimum load of processors witch make this approach more generic and the parallelization strategy more portable. In fact portability is increasingly cited as a desirable goal in parallelization strategy conception.

Definition Definition Definition Definition ( ( ( (14 14 14 14) ) ) ): : : : Portability

It is a property which assures that parallel programs are both code portable and performance portable to various parallel machines.

Portability:

% 100 * ) ( B n T n n S S t b = → Ρ Ρ = [ ] [ ] % 100 * ) ( ) 1 ( ) ( ) 1 ( BS TS N N N N N N α α - - - -
Despite disagreement about the exact meaning of "Portability", we can consider (Def.14) as working definition. According to James D. Mooney [START_REF] Mooney | Strategies for Supporting Application Portability[END_REF] the primary goal of portability is to facilitate the activity of porting an application from an environment in which it currently operates to a new or target environment. This activity has two major aspects:

(i) Transportation: physical movement of the program's instructions and data to the new environment. (ii) Adaptation: modification of the information as necessary to work satisfactorily in the new environment.

We skip "Adaptation" which involves higher level modifications that might be necessary to adjust the program to work with aspects of the new environment that are intentionally or unavoidably different from the old one. We focus on physical transportation which includes use of compatible media or communication channels between processor, and interpreting and translating file formats, character codes, data representations, processor design. Standard languages and portable compilers bridge the gap between programs and the variety of CPU interfaces that exist in target environments. However, many of these mechanisms still define only part of the environment interface that many applications need. Elements such as file structures, memory management, or especially asynchronous event handling are not adequately defined by most language standards or library specifications. When requirements for communication, concurrency, or timing constraints exist, conventional languages are clearly insufficient.

Here is another major challenge in multi-core multithread architecture programming. In an ideal case, moving from one-core to multi-core should provide n fold increase in computational power. But practically, it is something that never happened. In fact, all existing computational problems cannot be efficiently parallelized without incurring the costs of inter-processor coordination.

This kind of analysis was evoked in many researches. Let's come back to (Def.10) and we focus on Amdahl's law [START_REF] Ramaswamy | A framework for exploiting task and data parallelism on distributed memory multicomputers[END_REF]. It captures the notion that the extent to which we can speed up any complex work is limited by percentage of the sequential part in the executed work. Amdahl Through this formula, for the given problem and an eight-core machine, Amdahl's law says that even if we parallelize 90% of the solution, but not the remaining 10%, then we end up with only four-fold speedup, and not the expected eight-fold speedup. In fact, these additional parallel parts involve substantial communication and coordination.

In our dynamic parallelization strategy, as we explained later in nest section, each two threads will share only one FIFO queue in order to push neighbors of lowered pixels. Intuitively we are going to opt towards a solution with a simple lock-based shared FIFO queue. Associated push and pop methods will be synchronized by a mutual exclusion lock. Even if this implementation is a correct concurrent FIFO queue, because each method accesses and updates fields while holding an exclusive lock, the method calls take effect sequentially. And according to Amdahl's law, this sequential communication can substantially affect the performance of our program as a whole. In multi-core architecture, such synchronization technique can also be the origin of costly overheads. Even if we opt to second method based on lock-free solution [START_REF] Anderson | Real Time computing with lock-free shared Object[END_REF] in order to minimize the overheads, it is demanded that at least one thread (of all the threads that are executing the push or pop function at one moment) is progressing (inserting or extracting pixels from or to the FIFO queue). Unfortunately, we do not know in advance how many parallel threads will call push or pop functions. And method calls still take effect sequentially. Other solution is wait-free technique [START_REF] Herlihy | Wait-Free Synchronization[END_REF], it is required that a process finishes within a finite number of execution steps. Something that we cannot also guarantee because we cannot predict how many points will be characterized and then how many pixels will be inserted in the FIFO queue.

Finally we decide to move to spin-wait mechanism [START_REF] Anderson | The performance of spin locks alternatives for shared money multiprocessors[END_REF], illustration is given by (fig. 6), and a thread waiting to push an item might spin for a brief duration without being added to the queue of waiting threads. As a result, the thread is effectively put to sleep without relinquishing the remainder of its CPU time slot. It is potentially more efficient to spin and wait, instead of using either lock-free or wait-free mechanisms, because those force a thread context switch, which is one of the most expensive operations performed by the operating system.

The 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 (Algo. 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.

An illustration of the merging algorithm with four threads is given by (fig. 7). The original shared data structure is divided into 4 research areas{

} 1 2 3 4 , , , z z z z . Threads { } 1 2 3 4
, , , T T T T will start browsing different zones in parallel.

If ( ( _ ) FIFO Q ≠ ∅ then 8.
Usedstatus(

_ j FIFO Q ,true) ; 9. Insert_neighbors( i T , D P , _ j FIFO Q ); 10. Else : add_new_fifo ( _ FIFO Q ); 11.
Usedstatus

( 1 _ j FIFO Q + ,False) ; 12. Insert_neighbors ( i T , D P , 1 _ j FIFO Q + );

13.

EndIf ; EndIf ; 14. EndFor;

Conclusion

In this chapter we have identified needs for a common parallelization strategy of topological operators on shared memory parallel architectures. After studying theoretical basis of parallelization strategies, we presented our approach called SD&M. 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 ½[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[END_REF] team introduced alpha-extension (diffusion) and betareduction (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.

TOPOLOGICAL WATERSHED

he watershed is now used as a fundamental step in many powerful segmentation procedures. Judged by the great diversity of applications in which this transform appears relevant, watershed can be seen as one of the most popular segmentation tools coming from the field of mathematical morphology. In general, five large classes of algorithms to compute watershed transform can be figured out. The first one is based on flooding approach, the second is based on a path-cost minimization, third class is based on topological approach, fourth class is based on local condition and finally a fifth class based on minimum spanning approach. These classes will be studied in details in the second section of this chapter.

In addition to this study, we will focus also in the parallel implementation of watershed. In fact, the difficulty of such implementation resides in the non-local criteria of the watershed transform. The decision whether a pixel belongs to a basin cannot be based on purely local considerations as Roerdink and al. [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF] introduced. He also explains that results given by some algorithms depend on the order in which pixels are treated during execution. In the sequential case, fixing the scanning order can resolve this problem, so that a deterministic result is obtained. In a parallel implementation this is no longer true since the outcome depends on the relative time instants at which different processors treat the pixels, and this is unpredictable in the case of asynchronous processors. Throughout this section, we will leave in search of more adapted parallel algorithm to compute watershed transform suited for shared memory parallel machines.

We start by presenting an intensive study of all existing watershed transform in the discrete case in section (3.1). This study concludes with a classification of different algorithms according to criteria of recursion, complexity, basins computing and topology preservation in section (3.2). Once best approach identified, we present, in section (3.3), a new algorithm to compute watershed that is parallel, preserves the topology of the input image, does not need prior minima extraction and suited for SMP machines. Based on stream notion introduced by Cousty [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF], this algorithm does not require any sorting step, or the use of any hierarchical queue. Experimental analyzes such as execution time, performance enhancement and cache consumption are also presented and discussed.

Watershed transformations

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 T 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 8 gives a very symbolic description of the mentioned approach. In fact, it shows trends watershed transform use 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 chapter, 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 by 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]. In the following we present our analysis of watershed transform (WT) 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 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.

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 (Alg.

) ( / 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 (Alg. 4) of the flooding watershed. This implementation uses FIFO queue and it needs two steps:

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

(2) Flooding process: All nodes with grey level h are first given an initial label. Then those nodes that have labeled neighbors from the previous iteration are inserted in the queue, and then 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] Roerdink and Meijster have removed two points of inconsistency in Vincent algorithm's recursion. (i) Only pixels with grey value h are masked for flooding (line 13), instead of all non-basin pixels of (level h ≤ ), as (Def.16) would require. This explains why labels of 'wshed-pixels' (line 15) are also propagated with labels of catchment basins. (ii) If a pixel is adjacent to two different basins, it is initially labeled 'wshed'. But it is allowed to be overwritten at the current grey level by another neighbor's label, if that neighbor is part of a basin (lines [START_REF] Reinders | Intel Threading Building Blocks[END_REF][START_REF] Ramaswamy | A framework for exploiting task and data parallelism on distributed memory multicomputers[END_REF].

They also propose some modification to implement the recursion (Def.16) 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 [START_REF] Buzbee | Applications of MIMD machines[END_REF][START_REF] Foster | Designing and Building Parallel Programs[END_REF][START_REF] Chandra | Parallel Programming in OpenMP[END_REF][START_REF] Reinders | Intel Threading Building Blocks[END_REF][START_REF] Ramaswamy | A framework for exploiting task and data parallelism on distributed memory multicomputers[END_REF][START_REF] Amdahl | Validity of the Single Processor Approach to Achieving Large Scale Computing Capabilities[END_REF][START_REF] Gustafson | Reevaluating Amdahl's Law[END_REF][START_REF] Acar | The data locality of work stealing[END_REF][START_REF] Natarajan | Imprecise and Approximate Computation[END_REF][START_REF] Tilborg | Foundations of Real-Time Computing: Scheduling and Resources Management[END_REF][START_REF] Leung | REAL-TIME SCHEDULING ANALYSIS[END_REF][START_REF] Kolivas | RSDL completely fair starvation free 64 interactive cpu scheduler[END_REF][START_REF] Molnar | Modular Scheduler Core and Completely Fair Scheduler[END_REF][START_REF] Mooney | Strategies for Supporting Application Portability[END_REF][START_REF] Anderson | Real Time computing with lock-free shared Object[END_REF][START_REF] Herlihy | Wait-Free Synchronization[END_REF] have to be slightly changed. [ ]

; INIT p lab ← 9.
[ ] ; 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 (Def.15), they propose the following one: Definition (17): Threshold set of (f) at level (Alt) [Shengcai and [ ] ; curlab r lab ← pixels susceptible to be dumped, will be. According to the user requirement, ) ( p Diff can be even a constant function or a function computed according to the local characteristic of pixel p . In their paper, Shengcai and Lixu define as the following.

0 ← p dist //dist
) ( p Neighbor will denote the set of all p neighbors. And conn will refer to the predefined connectivity:{ } 8 , 4 .

Definition (18): Set of all connected neighbor [Shengcai and Lixu]

∑ ∈ - = ) ( ) ( ) ( ) ( p Neighbor q conn q f p f p Diff
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 brain crop frame (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 is 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 label Lab .

Topographic watershed process can be described by (Alg. 5). Let us note that priority level when inserting neighbors (line 2) corresponds to the gray level of the pixel. In line 6, only neighbors that have already been labeled 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 labeled points set.

For mathematical formulation of the mentioned process, we follow 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. And 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 the edge ) , ( q p . Then the lower slope ) ( p LS of f at a pixel p and the cost for walking from a pixel p to a neighboring pixel q can be defined as following: Definition ( 19): Lower slope of (f) at a pixel p [Roerdink and al.]

{ }         - = ∪ ∈ ) , ( ) ( ) ( ) ( ) ( q p d q f p f MAX p LS p p N q G

Definition (20): Cost functions [Roerdink and al.]

( )

       = + < > = ) ( ) ( / ) , ( . ) ( ) ( 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 also 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 . Note that each pixel, which is not in a minimum, has a neighbor of lower grey value (

LC f f = * ).

Definition (21): Topographic watershed [Roerdink and al.]

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 Fernand Meyer one. In the following we present a variant of Meyer algorithm with integrate the lower slope (Def. [START_REF] Brinch-Hansen | Distributed Processes: A Concurrent Programming Concept[END_REF]) of the input image as introduced Roerdink and al. [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF]. [ ] ;

0 ← v lab 5.
[ ] [START_REF] Bertrand | On Topological Watersheds[END_REF]. for all (v∈NG(u)) in the future (w.r.t scan order) for u do [START_REF] Wilkinson | Concurrent Computation of Attribute Filters on Shared Memory Parallel Machines[END_REF].

; ∞ ← 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 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 11) i k < <constitute a path. In the following we will denote by ( ) org π the origin i t of π and by ( ) dst π the destination k t of π . Now let assume a given function f that assigns to each path π a path cost ( ) f π , in some totally ordered setυ of cost values. 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.

Definition (22):

Max-arc path-cost function [Falcão and al.] max ( ) ( )

f t h t = { } max max ( , ) max ( ), ( , ) f s t f w s t π π ⋅ =
For IFT use, a specific function ( ) S f π can be defined since the search to paths start in a given 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 (23): Spanning forest [Falcão and al.]

For any pixel t J ∈ , a spanning forest P defines a path *( ) Falcão and al. algorithm describes IFT computing. Its algorithm is based on Dijkstra's procedure [START_REF] Dijkstra | A Note on Two Problems in Connection with Graphs[END_REF] 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 twith finite cost ( *( )) s f P t . Pixels will belong to a tree whose root is a seed pixel.

P
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 (Def. 22). We remember that ( ) h t is fixed but arbitrary handicap cost for any paths 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.

Definition (24):

Arc weight functions

1 ( , ) ( ) ( ) w s t J s J t = -

( , ) ( ) w s t G t =

There are two usual arc weight functions (Def 24): (w 1 ) where ( ) J s refers to the intensity of pixel of s . In that case IFT-Watershed is said by dissimilarity and (w 2 ) where ( ) G t is the morphological gradient of ImgIn at pixel t. In that case IFT-Watershed is said on gradient.

In [START_REF] Lotufo | IFT-Watershed from Gray-Scale Marker[END_REF] Lotufo and al. introduced IFT-watershed from markers, (Algo. 8), 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 the following formula. Note that, the maximum of the set { } p ← remove from FIFO 7.

for each

( ) 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 ←

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 global process to compute topological watershed, (Algo. 9), we define some basic notions that will be resumed in next paragraph for mathematical formulation : Let F be 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. Thus a basic definition of the topological watershed and a global description of its computing process can be given by (Def. 26) and (Alg. 9). Note that this process is repeated on loop until no more 'w-destructible' point remains.

Definition (26): Topological watershed [Couprie and Bertrand]

A topological watershed of F is a 'w-thinning' of F which contains no W-destructible point. The major feature of this transform is to produce a grayscale image.

Algorithm 9 : Topological watershed process

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

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

For deep 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 of topological watershed (Def. 28) relies on the following particular notion of simple point.

Definition (27): Simple point [Couprie and Bertrand]

Let

( ) Γ = , E G be a graph and let E X ⊂ . A point X x ∈ is simple (for X ) if the number of connected components of { } x X ∪ is 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 .

Definition (28): Topological watershed [Couprie and Bertrand]

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 non-destructible forW ).

Actually checking whether a point is 'w-destructible', or not, cannot be done locally if only available information are graph ( ) , E Γ and a function F . Couprie and al. [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF] propose a new algorithms (Algo. 10) 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.

Algorithm 10 : w-destructible function [Couprie and al.]

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 (Def. 28) 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 (Algo. 11). 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 11 : HighestFork Function [Couprie and al.]

Input : C : a component tree, V : a set of components of C 1. [ ] ← 1 1 , c k ) min(V ; // let [k 2 ,c 2 ]...[k n ,c n ]be the other elements of V 2. 1 k k m ← ; 3. 1 c c m ← ; 4. for ifrom 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. If ) ( 1 k k m = then return [ ] φ , ∞ else return [ ] m m c k , ;
The HighestFork function returns the highest fork of the set V , 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 the 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 'wdestructible' points in the appropriate order. Here, we present their specific linear watershed algorithm which may be used when the grayscale range is small.

Algorithm 12 : Topological watershed [Couprie]

Input: Ψ ), ( , F C F ; Output: 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] Meyer | Topographic distance and watershed lines[END_REF][START_REF] Najman | Watershed of a continuous function[END_REF] 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]. In witch Audigier and al.

start by presenting the following catchment basin formulation

{ } ( ) V, L( ) L( ) LC i i CB m v v m = ∈ =
, since local condition watershed assigns to each pixel the label of some minimum i m . Thus watershed can be defined as follow. We recall that the condition of '

{ } ( ) steepest v ≠ P

' means that ' v ' has at least one lower neighbor.

Definition (29): Watershed based on local condition [Audigier and Lotufo]

For any lower complete image CB L , a function L assigning a label to each pixel is called watershed segmentation if: a) ( )

( ) i j L m L m i j ≠ ∀ ≠ /{ } k m is the set of minima of LC L . b) For each pixel v with { } ( ) steepest v ≠ P , ( ), L( ) L( ) steepest p v v p ∃ ∈ = P .
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 topographic 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 = ∈ = .

For all (

) p D ∈ do //Simulation of sliding for all C 1 pixels 2.

( )

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 = = ∈ 6.

( )

SlidingList p S = 7. Q p ← 8. ( ) 0 GrowingDist p = 9.
Else if : 10.

( )

SlidingList p φ = 11.
End if 12. End for 13. ( )

While

h G p = 17.
For all ( ( )) q Neighbor p ∈ and ( ( )

) G q h = do 18.

If ( ( ) )

SlidingList q φ = then 19.

Append ( ) p to ( ) SlidingList q 20.

( )

GrowingDist q d = 21. Q ; q ← 22. Else If ( ( ) ) GrowingDist q d = then 23.
Append ( ) p to ( ) SlidingList q ( ) S SlidingList p =

3.

For all ( ) q S ∈ do

4.

Resolve(q) 5.

End for 6.

If S has a unique label α then 7.

( )

L p α = 8. Else 9. ( ) _ L p RIDGE label = 10. End if

End if

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. Actually, the beginning was with Meyer [START_REF] Meyer | Minimum Spanning Forests for Morphological Segmentation[END_REF] who proposes to 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.

For mathematical foundations, we will follow 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 later. 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.

Definition (30): Minimal altitude of an edge [Cousty and al.]

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 - - - = .

Algorithm 16 : Watershed-cuts algorithm [Cousty and al.]

Input : (V,E,F) : Edge-weighted graphs; Output: Ψ: a flow mapping of F

1. Foreach ( V) x ∈ do ψ(x) NO_LABEL ← ; 2. _ 0 nb labs ← 3. Foreach ( V)
x ∈ such that (ψ(x)=NO_Label) do

4.

[ ] L,Lab Stream(V,E,F,ψ,x) ← Else Foreach ( L) y ∈ do ψ(y) labs ←

;

Previous algorithm computes a watershed-cut using a stream function that is described as the following:

Algorithm 17 : Stream function [Cousty and al.]

Input : (V,E,F) : Edge-weighted graphs; Ψ : a label of V; x : point of V; Output: [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).

Classification of watershed algorithms

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 transforms. Selected criteria are justified by our objective to identify the most suitable algorithm for parallel implementation.

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]. 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.

Basic notions and definitions

Based on Cousty approach [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF], we define and illustrate stream notion. For the sake of simplicity, we restrict ourselves to the minimal set of notions that will be useful for our purpose.

We denote by V an edge-weighted graph. Let L V ⊂ . We say that L is a stream if, for any two points x and y of L , there exists, in L , either a path from x to y or from y to x , with steepest descent for F . We recall that steepest descent principal is already defined on section 3. L , and there is from y to x a path L , with steepest decent for F . If

1 2 L L ≺ the L is a stream. If there is no stream under L , L is considered as an stream - ≺ . Now any stream L which contains streams - ≺ is itself an streams - ≺ .
Basic illustration of stream notion is given by (fig. 9

Parallel stream computing

For that propose, following algorithm will assigns, in parallel way, a label to each point of the graph. Actually, from each non-labeled point x , a stream L composed of non-labeled points and whose top is x is computed. It is important to mention that streams computing at this level are completely independent than streams can be completely computed in parallel, see (fig. 11). For N point (x 1 , x 2 … x n ), their associated flows are simultaneously extracted: The stream function (Alg. 17) halts at line 16 when all bottoms of L have been explored or, at line 8, if a point z already labeled is found. In the former case, the returned set L is an stream -≺ . In the latter case, the label lab of z is also returned and there exists a bottom y of L such that , y z is a path with steepest descent. Thus, there is an stream -≺ 1 L , under L , included in the set of all vertices labeled lab.

(L 1 , L 2 … L n ).
Remark that, in stream function, the use of breadth-first iterations is required to ensure that produced set L is always an stream -≺

. Otherwise, if only depth-first iterations were used, Stream could be stuck on plateaus (connected sub-graphs of G with constant altitude) since some bottoms of L would never be explored. [L i ,lab i ] ← Stream (V, E, F, Ψ,x i ) ; // to get associated stream for each x i 6. nb-fusion = i ; 7. while ( nb_fusion != 1 ) 8.

for (j =0 ; j <= nb_fusion ; j+=2) do // lunch (nb_fusion/2) process at once 9.

if (L j ∩ L j+1 ) = ∅ then s-labeling ([L j ,lab j ] , nb_labs) ; 10.

s-labeling ([L j+1 ,lab j+1 ] , nb_labs) ; 11.

else f-labeling ([L j ,lab j ] , [L j+1 ,lab j+1 ] , nb_labs) ;

12. nb-fusion = nb-fusion / 2 ;

The procedure, in the parallel case, is based on the idea of labeling and merging two obtained flows at once. If two flows (to merge) 'L i ' and 'L i+1 ' contain no common summit, (L i ∩ L i+1 )=∅, meaning there are no common wells between the two sources 'x i ' and 'x i+1 ', in foreach (y ∈ L) do Ψ(y) ← lab ; 6. return NULL If the two flows (to merge) 'L i ' and 'L i +1 ' contain common summit, (L i ∩ L i+1 ) ≠ ∅, meaning there are common wells between the two sources 'x i ' and 'x i +1 ', The major problem in concurrent merging of multiple flows is summed up in labels assignment.

If two streams share the same well, which label should be given to involved pixels? Proposed solution is inspired from the flooding paradigm. Indeed, we start by studying all possible cases of merging two water streams gushing from different sources, see (fig. 13). Our goal is to identify which stream will be the first to reach the well. This latest will mark the well by its own label. The starting point is the steepest decent approach with the following conditions: (i) Water flow rate is identical for all sources, (ii) Flow surface is perfectly smooth and (iii) Runoff velocity is uniform for each flow If these conditions are fully met, three factors come into play to determine the flow velocity: the source altitude, distance between source and sink, and finally the slope. In fact, topographic slope particularly influence the runoff. The inclination of the slope is surely the most important

X1 X2 Xk-1 Xk Xi Xj Xm Xn B A C
topographical aspect. Normally, its impact is limited on short slope. It is more visible on longer slope even if runoff needs a certain distance to reach its maximum velocity. Mathematic formulation of flow medium speed can be is given by the Chezy formula:

( )

1/ 2 * c V c h s =
introduced in 1769. 'C' refers to roughness coefficient of Chezy. 'S' refers to the slope, and 'h' refers to the altitude of the source. If we draw the truth table with these three factors (d: distance, s: slope, s: altitude), by varying one parameter each time, we can identify only five possible cases: The two streams have the same altitude, slope and distances that separates sources from well. In the 2nd case, both flows traverse the same distance but slopes and sources altitudes are different. In the 3rd case, the two streams run down the same slope but they travel different distances since sources' altitudes are different. In the 4th case, the altitude is the same for both sources, but traveled distances and slopes are different. Finally, the altitude of the sources, the distances separating them from the well and the slopes of followed paths are different for the two streams.

The question now is: does one of these five situations necessarily appear when merging. If we are dealing with two stream -≺ , this problem does not arise because we are forced to generate a new label for identified wells (line 1, Algo. 20). Also, if one of the two streams includes an stream -≺

, it means there exist a label already generated that we can assign to the common wells (line 5, Algo. 20). Finally, if both streams include stream -≺ then two labels already exist. In order to decide which one to assign, we compute approximately the flow's average speed using Chezy formula. It is important to mention that gray level of a pixel represent its altitude. Slope and distance between sources and wells can be computed trough pixels coordinate. According to fixed conditions, roughness coefficient is equal to one. 19. Return Null

Performance testing

In this section we present an overall assessment of the parallel watershed operator. We begin by presenting test conditions. Then, obtained results in terms of execution time and cache use are presented and discussed. Based on these results, we compute efficiency and scalability of our implementation. We enhance discussion on scalability by computing the amount of work required to reach the average speed. Unfortunately, portability of our application will not be assessed for purely technical reasons. We note that execution time drops from an average of 4636 ms with a single thread on one CPU down to 713 ms with 8 threads on 8 CPUs. The speed up was computed using formula Ts/Tp with Ts for 1 CPU = 4360 ms (definition 11 chapter 2). A remarkable result about speedup is also shown in table 4. In fact, speed-up increases as we increase the number of threads beyond the number of processors in our machines. In the first implementation, using two CPUs, the speedup at 2 threads is 1.37 ± 0.01. However, for the second implementation, using 8 CPUs, the speedup has increased to 6.11 ± 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. For better illustration we establish execution time and speedup curve (see fig. 14). In the following we present our experimental analysis. We consider a commonly used Intel processor configuration. Number of processor varies from one to eight. The frequency varies between 1,73 GHz and 3,4 GHz as shown in table 5. 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 8 Mbytes, and the associativity varies between two and sixteen ways.

As a result of this experiment, fig. 15 (A), we found that two 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. 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 instruction, L1 and L2 data misses on four architectures using SMP scheduling policy. We note that number of instruction increase from an average of (34x10 6 ) instr. on 1 CPU to (790x10 6 ) instr. on 8 CPUs.

To highlight cache performance, we compute wait status which refers to the delay experienced by processor when accessing external L2 caches each time that information is missing in L1.

Since L1miss is followed either by an L2hit (success) or L2miss, wait status can be computed by following formula: the sum of L2hit and L2miss.

We suppose that L2 access time is estimated at 10 cycles (in hit case) and 100 cycles (in miss case). In the sequel, we turn to efficiency evaluation, using def. 13 (with ts = 360ms), in order to describe exploitation degree of each processor in used SMP machines. As introduced in [START_REF] Kumm | Parallel Computing Efficiency: Climbing the Learning Curve[END_REF], this profiling will highlight limitations introduced by parallel watershed implementation on SMP machines. Indeed, the efficiency decreases of 30% when switching from mono core architecture to dual cores architecture. Despite a slight increase on quad cores architecture, the efficiency is 20% lower than that measured with 1Cpus. More details are shown in Causes for losses of efficiency can be explained by the following reasons, partially introduced in [START_REF] Jain | Fundamentals of Digital Image Processing[END_REF] as parallel computing delays: (i) I/0 delays due to the need to distribute parallel data across local PE data stores. (ii) Communication delays, due to the need for PES to access data which is not located in their local data stores. (iii) Set-up delays due to the set-up of control and processing logic and the network for inter-PE communication. In further evaluation, we extend speedup profiling of parallel watershed computing into scalability analysis. According to Intel theoric study [START_REF] Fruehe | Planning Considerations for Multicore Processor Technology[END_REF], very high scalability can be achieved on multicore architecture. By way of example, dual-core architecture offers a scalability of roughly 80% for the second processor, depending on the OS, application, compiler, and other factors. That means the first processor may deliver 100% of its processing power, but the second processor typically suffers some overhead from multiprocessing activities. As a result, the two processors do not scale linearly. Thus, a dual-processor system does not achieve a 200% performance increase over mono-core architecture, but instead provides approximately 180% of the performance that a single-processor system provides.

In our evaluation framework, we first introduce the average unit speed. This parameter, seen as the ration between achieved speedup and the number of processor, will be very useful to determinate scalability. We can also extend this definition into the maximum average speed which is defined as ratio of maximum achieved speedup by number of processor. max( ) ( ) ( )

speed
N W W N × > , thus ( , ') 1 N N ψ < .
To calculate different values of efficiency foreach architecture, we must first determine the necessary amount of work ' W , as shown in As the number of thread increase, a linear speedup has been observed (see fig. 17). Also, the speedup improves as the problem sizes increases. Note that when number of thread exceed number of cores, total execution time dramatically reduces. The difference between each efficiency curve with the ideal curve (constant efficiency equal to 1) decreases as the number of thread increases.

Conclusion

In this chapter, 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.

First contribution is the global nature of the proposed study. In fact, for each approach, we give informal definition, then we presented processing procedure followed by mathematical foundations and the algorithm of reference. Recent publications based on some approach are also presented and discussed.

Second contribution concerns classification of watershed algorithms according to criteria of recursion, complexity, basins computing and topology preservation.

Third contribution concern a new algorithm to compute watershed that is parallel, preserves the topology of the input image, does not need prior minima extraction and suited for SMP machines. This algorithm does not require any sorting step, or the use of any hierarchical queue.

A global description of the computing process is given by fig. 18. Through this illustration, we show links between parallel watershed-cut and the SD&M strategy application. In fact, splitting step is applied directly on input graph when selecting sources. Unlike conventional technique of division such as pixel division, or block division, the source selection is completely random. Associated steam computing is fully parallel (read mode data accesses). Then distribution depends only on the available processors. This flexibility in data manipulation allowed us to obtain very good results especially in terms of efficiency fig. 14 (b) without using the 'Basic-NPS' scheduler. Finally, the merging step contains procedures of s-labeling and f-labeling. Through these two functions, we have remained confident in our approach for merging streams two by two (algorithm already introduced in section 2.4.3. Experimental analyzes such as execution time, performance enhancement, cache consumption, efficiency and scalability are also presented and discussed.

Note that our algorithm can't be applied directly over grayscale image. Actually three major steps are needed: the passage from grayscale image to edge-weighted graphs, then the application of the parallel watersheds-cut algorithm on the plot and finally the visualization of the graph in the Khalimsky space [START_REF] Khalimsky | Topological graph theory foundations of design and control in multidimensional discrete systems[END_REF][START_REF] Khalimsky | Topological structures in computer science[END_REF].

TOPOLOGICAL THINNING

n many computer vision applications, standard techniques of pattern recognition are thinning algorithms. As a preprocessing stage, these algorithms have been used for the recognition of handwriting or printed characters, fingerprints, chromosomes and biological cell structures, etc. [START_REF] Jain | Fundamentals of Digital Image Processing[END_REF]. Topological thinning and skeletonization are ones of the most cardinal operators for this kind of preprocessing. In literature, several 2D parallel thinning methods can be found, see [START_REF] Arcelli | Parallel thinning of binary pictures[END_REF][START_REF] Bernard | Improved low complexity fully parallel thinning algorithm[END_REF][START_REF] Chin | A one-pass thinning algorithm and its parallel implementation[END_REF][START_REF] Hall | Fast parallel thinning algorithms: Parallel speed and connectivity preservation[END_REF][START_REF] Pavlidis | An asynchronous thinning algorithm[END_REF]. Proving that such an algorithm always preserve topology is not an easy task, even in 2D. We say that an algorithm preserves topology if obtained skeleton through thinning method has the necessary information to reconstruct the original image. The proofs found in the literature are often combinatorial and hardly extendable to 3D, a fortiori to higher dimensions.

Couprie [START_REF] Couprie | Note on fifteen 2d parallel thinning algorithms[END_REF] present a study of fifteen parallel thinning algorithms, see table 13, based on the framework of critical kernels. He proves that ten among these fifteen algorithms indeed guarantee topology preservation, and give counter-examples for the five other ones. He also investigates, for some of these algorithms, the relation between the medial axis and the obtained homotopic skeleton. Z. Guo and R.W. Hall [START_REF] Pavlidis | An asynchronous thinning algorithm[END_REF] (Second version) 1992 10 Z. Guo and R.W. Hall [START_REF] Pavlidis | An asynchronous thinning algorithm[END_REF] (Third version) 1992 11 B.K. Jang and R.T. Chin [START_REF] Jang | One-pass parallel thinning: Analysis, properties, and quantitative evaluation[END_REF] 1992 12 B.K. Jang and R.T. Chin [START_REF] Rutovitz | Pattern recognition[END_REF] 1993 13 U. Eckhardt and G. Maderlechner [START_REF] Pavlidis | A flexible parallel thinning algorithm[END_REF] 1993 14 S.S.O. Choy, C.S.T. Choy and W.C. Siu [START_REF] Choy | New single-pass algorithm for parallel thinning[END_REF] 1995 15 T. Bernard and A. Manzanera [START_REF] Khalimsky | Topological structures in computer science[END_REF] 1999 16 M. Couprie, F. N. Bezerra and G Bertrand [START_REF] Couprie | Topological operators for grayscale image processing[END_REF] 2001

Table 13 : Classification of thinning algorithm according to topology preservation

According to the above study, algorithms proposed by Pavlidis in 1981 [START_REF] Pavlidis | An asynchronous thinning algorithm[END_REF][START_REF] Pavlidis | A flexible parallel thinning algorithm[END_REF], by Chin and al. in 1987 [87], by Holt and al. in 1987 [93], by Hall in 1989 [START_REF] Hall | Fast parallel thinning algorithms: Parallel speed and connectivity preservation[END_REF], by Guo and Hall in 1992 [START_REF] Guo | Fast fully parallel thinning algorithms[END_REF] (3 variants), by Jang and Chin in 1993 [START_REF] Jang | Reconstructable parallel thinning[END_REF], by Eckhardt and Maderlechner in 1993 [START_REF] Eckhardt | Invariant thinning[END_REF], and by Bernard and Manzanera in 1999 [START_REF] Bernard | Improved low complexity fully parallel thinning algorithm[END_REF] preserve topology. However, algorithms proposed by Rutovitz [START_REF] Rutovitz | Pattern recognition[END_REF], by Zhang and Wang [START_REF] Zhang | A modified parallel thinning algorithm[END_REF], by Wu and Tsai [START_REF] Wu | A new one-pass parallel thinning algorithm for binary images[END_REF], by Jang and Chin [START_REF] Jang | One-pass parallel thinning: Analysis, properties, and quantitative evaluation[END_REF] and by Choy and al. [START_REF] Choy | New single-pass algorithm for parallel thinning[END_REF] produce a Skeleton that don't allows to reconstruct the original image. We propose to extend this study by adding a new algorithm proposed by Couprie and al. [START_REF] Couprie | Topological operators for grayscale image processing[END_REF]. This algorithm is also able to act directly over grayscale image without modifying topology.

I

Based on this first 'classification' of 2D thinning algorithm, we go further in this study and we start by proposing new classification, in section 4.2, that takes into account not only the criterion of topology preservation but also the proximity of the obtained skeleton from the medial axis, necessary execution time and cache consumption. Since suited algorithm for parallelization is identified and convinced that best performance can be achieved, we propose a concurrent implementation of a powerful topological thinning algorithm that is originally introduced by Couprie [START_REF] Couprie | Topological operators for grayscale image processing[END_REF]. In section 4.3, we start by introducing theoric background, we give also some illustration of original algorithm before presenting parallel lambda skeleton procedure. Based on SD&M strategy, Distributed work during thinning process is done by a variable number of threads. Tests on 2D grayscale image (512x512), using shared memory parallel machine (SMPM) with 8 CPUs cores (2 × Xeon E5405 running at frequency of 2 GHz), showed an enhancement of 6.2 with a maximum achieved cadency of 125 images/s using 8 threads.

(1) (

Classification of thinning algorithms

In order to make a quantitative evaluation of selected algorithms, we start by presenting two selection criteria based on Jang and Chin's study [START_REF] Jang | One-pass parallel thinning: Analysis, properties and quantitative evaluation[END_REF]. Authors propose to compute ratio between the number of pixels contained in the maximal disks (obtained from the skeleton) and the effective number of pixels in the original image. This measure, see (Def. 31), is used to determine the proximity of the skeleton from the medial axis as proof of skeletal connectivity and convergence. () Area function is a pixel's counter, ' S refers to the number of pixels contained in the maximal disks and S refers to the effective number of pixels (in the original image).

Definition (31): Proximity of the skeleton from medial axis (

M m ) ' ( ) ( ) m Area S M Area S =
Inspired by Jang and al. study, we adopt this criterion of skeleton proximity from medial axis with some modification, see (Def.32). Actually, we privilege the ratio between number of pixels of the skeleton which belong to the medial axis (Ai) and number of pixels in the skeleton (Ni). Our goal is to identify the algorithm that returns the most centered skeleton. Thus, most faithful algorithm is the one whose C1 value is closest to 1. An initial assessment, see table 14, using three different images shown in fig. 19, has allowed us to draw first curve (see fig. 20). Unfortunately this quantitative measure is not sufficient to evaluate algorithms. Actually, even if (C1) tends to 1, the gap between the number of skeleton's pixels belonging to medial axis and the number of reference pixels of medial axis can be great. For example, let's consider two algorithms. The first one provides 50 pixels belonging to medial axis, the total number of skeleton's pixels is also 50 and the number of reference pixels is 100. The second algorithm provides respectively 100, 150 and 100. Thus, (C1) for the first algorithm is equal to 1 even if medial axis doesn't include many points from the medial reference axis; therefore, the deviation is wide 0.5. For the second algorithm, (C1) is over 1 even if all pixels of the median reference axis are included in reference median axis. To resolve this problem, we propose a second criterion (C2) to identify how mush generated Skeleton is symmetric and close to reference medial axis, see Def. 33. Thus we compute the ratio between number of pixels of the skeleton which belong to the medial axis (Ai) and number of pixels of reference medial axis (R i ).

Definition (32): Skeleton connectivity criteria (C

1 ) 1 ( ) ( ) i i Area A C Area N = Shape (1) Shape (2) Shape (3) N 1 A 1 N 2 A 2 N 3 A 3 1 T.
Definition [START_REF] Foster | Designing and Building Parallel Programs[END_REF] ) except that it has a bad factor (C 1 ) which doesn't exceed 0.7. Proposed algorithm by Chin [START_REF] Jain | Fundamentals of Digital Image Processing[END_REF] has the worst classification according to two criteria so it can be automatically rejected. Although, Holt [START_REF] Couprie | Note on fifteen 2d parallel thinning algorithms[END_REF] and Hall [START_REF] Arcelli | Parallel thinning of binary pictures[END_REF] algorithms have similar factors, Couprie [START_REF] Bernard | Improved low complexity fully parallel thinning algorithm[END_REF] and Jang [START_REF] Rutovitz | Pattern recognition[END_REF] shows that Hall algorithm is an enhanced version of Holt algorithm, therefore we chose it for assessment. The three versions of Guo [START_REF] Pavlidis | An asynchronous thinning algorithm[END_REF] algorithm have similar results but second version (Guo, Hall 1992 -b) has a better classification in the two graphs thus this version will be also selected. Jang [START_REF] Rutovitz | Pattern recognition[END_REF], Eckhardt [START_REF] Pavlidis | A flexible parallel thinning algorithm[END_REF] and Bernard [START_REF] Khalimsky | Topological structures in computer science[END_REF] have also a very good classification in both graphs and can be kept. Finally we can conclude by the upper comparison chart that ranks algorithm by qualitative descending factor [(C1) AND (C2)].

In the following, we conduct an initial evaluation of selected algorithms. Tested images are those shown in fig. 19 ( We note good data management when reading or writing. Algorithms that suffer from some problems of data management are rather those of Couprie, Guo and Eckhardt. These results are most visible especially when computing cache consumption for image 2. If we combine implementation features of these algorithms (execution time and cache use) with classification criteria previously presented (C 1 && C 2 ), we can draw following classification on table 18. Even if Guo and Hall algorithm has best performance, first three algorithms can be classified in the same family. In fact, they have many similarities and close performances. Second family is composed of two algorithms, Jang and Bernard. Despite their low cost of cache consumption, execution times are quite high. Based on Couprie study, we will develop more all of these five algorithms to identify the most suitable one for parallel implementation on SMPM.

For the sake of simplicity, we restrict ourselves to the minimal set of notions that will be useful for our purpose. Y ← set of pixels X which are GH-deletable 3.

X ← \ X Y Y ← set of pixels in X witch are both simple and perfect 3. To achieve this simplification; they introduce the notion of -destructible point which is more flexible then the notion of destructible point. In fact, a point is said to be a -deletable point (for F), λ being a positive integer, if it is either a -destructible point, or a peak point such that Fሺxሻ െ α ି ሺx, Fሻ λ. We remind that a point x is said -destructible if it satisfies one of the two following conditions:

X ← \ X Y
(i) x is destructible or x is -divergent. (ii) At least k-1 connected components c ୧ of Γ -ሺx, Fሻ are such that Fሺxሻ െ F ି ሺc ୧ ሻ λ, with i ൌ ሼ1, … , k െ 1ሽ.
Let X ؿ Z ଶ and x א X, x is an end point (for X) if #ሺΓ ୬ כ ሺxሻ ת Xሻ=1. Let F ؿ φ and x Ԗ Z ଶ , x is an end point (for F) if it is an end point for the set F ୩ with k ൌ Fሺxሻ. A point is said to be -end point (for F) if it is an end point for F and if: Fሺxሻ െ α ି ሺx, Fሻ ߣ. A ← medial axis of X Y ← set of pixels in X which match Bernard and Manzanera masks 3.

X ← \ X Y 4. Until Y φ =
Algorithmic structures, previously presented, are very similar. An iterative process is always launched until stability. In each iteration, some pixels are selected according to a set of criteria then output image is updated. The only algorithm that requires pre-treatment is Jang and Chin one (pixels of the medial axis are calculated). Unlike other algorithms, Couprie procedure doesn't exempt selected pixels from the output image but it changes their values. This will allow greater independence in the data processing if we ever plan to launch a parallel processing. For this reason, we selected the latter algorithm to propose a parallel version suitable for (SMPM).

Parallel lambda-skeleton algorithms

In this section, some basic notions of topological operators are first summarized; Illustration of the original algorithm is also introduced. Then, new adapted version of thinning is introduced. Experimental analyzes results of different implementations are also presented and discussed.

Theorical background

First, we recall some basic notions of grayscale images. A 2D grayscale image may be seen as a map ܨ from ܼ ଶ to ܼ. For each point ݔ א ܼ ଶ , ܨሺݔሻ is the graylevel value of .ݔ We denote by ߮ the set composed by all maps from ܼ ଶ to ܼ. Let ԋ א ߮, the section of ԋ at the level ݇ is the set ܨ composed of all point ݔ א ܼ ଶ such that ܨ ݇. As for the binary case, if we use the nadjacency for the section ܨ of ԋ , we must use ݊-adjacency for the section ܨ with ሺ݊, ݊) = (8,4) or (4,8).

We remind that for two points ݔሺݔ ଵ , ݔ ଶ ሻ, ݕሺݕ ଵ , ݕ ଶ ሻ ߳ ܼ ଶ , we consider that ݕ is 4-adjacent to

ݔ if ݕ| ଵ െ ݔ ଵ | ݕ| ଶ െ ݔ ଶ | 1, and ݕ is 8-adjacent to ݔ if max ሺ|ݕ ଵ െ ݔ ଵ |, ݕ| ଶ െ ݔ ଶ |ሻ 1.
In the following, we consider the two neighborhoods relations Γ ସ and Γ ଼ defined by, for each point x א Z ଶ ,

߁ ସ ሺݔሻ ൌ ሼݕ א ܼ ଶ | ݕ ݏ݅ 4 adjacent ݐ ݔሽ, ߁ ଼ ሺݔሻ ൌ ሼݕ א ܼ ଶ | ݕ ݏ݅ 8 adjacent ݐ ݔሽ.
For more general presentation, we will define Γ ୬ כ ሺxሻ ൌ Γ ୬ ሺxሻ\ሼxሽ. We will also denote by ԋ the complementary map of ԋ. We note that the complementary sets of the section of ԋ are section of ԋ. In all the rest of this paragraph, we will note n=8 for the section of ԋ , thus we must use ݊=4 for ԋ. It is also important to mention that a non-empty connected component ܺ of a section ܨ of ܨ is a (regional) maximum for ܨ if ܺ ת ܨ ାଵ ൌ  and a set ܺ ؿ ܼ ଶ is a regional minimum for ܨ if it is a regional maximum for ԋ. Let F א φ, the point x א Z ଶ is destructible (for F) if x is a simple for F ୩ , with k ൌ Fሺxሻ. We remind that a point x is said simple for X ؿ Z ଶ if Tሺx, Xሻ ൌ 1 and Tሺx, Xሻ ൌ 0 . We define also some associated connectivity numbers:

(i) (ii) (iii)
Furthermore, the connectivity numbers allow the classification of the topological characteristics of a point:

(i) is a peak point if . (ii) is a -divergent if .
A point is said to be a -deletable point (for ), being a positive integer, if it is either adestructible point, or a peak point such that . We remind that a point is said -destructible if it satisfies one of the two following conditions: is destructible or isdivergent and at least k-1 connected components of - are such that , with .

Let and , is an end point (for ) if =1. Let and , is an end point (for ) if it is an end point for the set with . A point is said to be -end point (for ) if it is an end point for and if: .

Illustration of original algorithm

, we say that is a skeleton of if is obtained from by iteratively selecting a destructible and non-end point in and lowering it down to , until stability. In order to get a filtered skeleton, that is to eliminate non significant branches and regional minima, Bertrand and Couprie allow -deletable and not -end to be lowered. It is important to mention that each time that a pixel is lowered, its eight neighbors must be reexamined to be sure that topology is still preserved. In Figure 1

Parallel thinning algorithm

Now, we present a parallel version of the thinning according to strategy previously discussed. Let the map F from Z ଶ to Z represent the input grayscale image. For each point x א Z ଶ , Fሺxሻ is the graylevel value of x. We denote by φ the set composed by all maps from Z ଶ to Z.

Let ԋ א φ, the section of ԋ at the level k is the set F ୩ composed of all point x א Z ଶ such that F ୩ k.

Let T be the set of type sought in the characterization of pixels. For thinning algorithm: T ൌ ሼλ-deletable and not λ-end pointsሽ. It is important to mention that points from T can also be end-point and isolate-point for crest restoring.

We will refer to global search space by Ime, and associated map (from Z ଶ to Z ) to each subspace Ime ୧ is F ୧ . For each point x א Z ଶ , F ୧ ሺxሻ is the graylevel value of x in the search space Ime ୧ .

The following dynamically parallel λ-Skeleton algorithm (it is adapted for two concurrent threads, but it can be easily extended to N threads) starts by dividing the search space. m ୧୬ and m ୱ୳୮ define sub-region bounds. Since the distributed work starts, each thread will lower each characterized pixel and then push its eight neighbors in E ୱ୬ . E ୱ୬ is the set of all selected neighbors and it is shared between only two threads. E ୱ୬ is the new defined set to explore since threads finished. Newly characterized pixels are pushed in a private set called E ୩୧ . The pixel set assigned to the newly generated thread is nothing else than E ୱ୬ and the associated search space is ሺሺIme ୧  Ime ୧ାଵ ሻ  E ୩୧  E ୩୧ାଵ ሻ.

Experimental analysis

The proposed parallel λ -Skeleton algorithm was implemented in C in two variants: the first implementation, based on a simple lock-based shared FIFO queue, using OpenMP critical directive. The second is based on a spin-wait FIFO queue, already introduced in section 2.4.2. Wall-clock execution times for numbers of threads equal to 1, 2, 4, 8, and 16, for each one of these implementations, were determined. The efficiency Ψሺ݊ሻ, already defined in section 2.4.1., is also measured.

Times were performed on eight-core (2× Xeon E5405) shared memory parallel computer of the Faculty of Electrical Engineering and Communication of Brno University, on Intel Quad-core Xeon E5335, on Intel Core 2 Duo E8400 and Intel mono-processor Pentium 4 660. Each processor of the Xeon E5405 and E5335 runs at 2 GHz and both of the two machines have 4 GB of RAM. The E8400 processor runs at 3 GHz. The Pentium processor runs at 3.6 GHz (see Table 17). The last two machines have 2 GB of RAM. The minimum value of 5 timings was taken as most indicative of the speed of the algorithm. The measurements were done on 2D grayscale image (512*512) of real brain MRI. Results of the two implementations are shown in fig. 28 (a)(b).

On the eight-core machine, wall-clock execution time for the first implementation using a lockbased shared FIFO queue drops from an average of 40.211 ms for a single thread down to 28.458 ms at 8 threads. For the second implementation using spin-wait FIFO queue, wall-clock execution time drops from an average of 41.889 ms for a single thread down to 8.282 ms at 8 threads. As expected, the speed-up for the second implementation using Private-Shared FIFO queue is higher than for the one using lock-based shared FIFO queue, because context changing were nearly eliminated. A remarkable result shown in (fig. 29) is the fact that the speed-up increases as we increase the number of threads beyond the number of processors in our machine (eight cores). For the first implementation, the speedup at 8 threads is 1.7 ± 0.05. However, for the second implementation For better readability of our results, we tested the efficiency of our algorithm on various architectures using the Ψሺ݊ሻ formula introduced earlier with fixed serial time equal to 48.247 ms. For parallel time we use best parallel time obtained using 8 threads. As can be seen in fig. 30, second implementation is more efficient that the first one in all architectures. It is also suitable to return to Amdahl's law, introduced in section 4, in order to explain obtained results.

In fact the global speed up formula is ܵሺ݊ሻ ൌ . Thus if the number of cores increases, the speedup also increases (more work can be done simultaneously with more threads). On the other hand the efficiency will decrease. that is parallel, preserves makes link between research area into will lower each characterized contains the set of Unfortunately, conventional given good results modest. Therefore we have applied our approach for better performance. faster than original version proposed by ), using shared memory parallel of 2 GHz), showed using 8 threads.

TOPOLOGICAL SMOOTHING

moothing 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, Couprie and Bertrand [START_REF] Couprie | Topology preserving alternating sequential filter for smoothing 2D and 3D objects[END_REF] introduced a new method for smoothing 2D and 3D objects in binary images while preserving topology. 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. Inspired by their approach, we propose a new algorithm for topological smoothing that is parallel and preserves topology.

This chapter is organized as follows: in section 5.1, some basic notions of topological operators are summarized; the original smoothing filter is introduced. In section 5.2, the new parallel smoothing method is introduced. Evaluations of acceleration, efficiency and success rate of cache memory access are presented and discussed in section 5.3. Finally, we conclude with summary in section 5.4.

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 More examples can be found in [START_REF] Couprie | Topology preserving alternating sequential filter for smoothing 2D and 3D objects[END_REF]. ( ) ( )

, 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 a 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.

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 each operator. Obtained execution time, efficiency, speedup and cache misses will be introduced and discussed.

As we have shown previous section, smoothing algorithm receives as input a binary image and maximum radius. It uses two procedures for homotopic opening and closing, see fig. 33 (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 (already defined in section 4.2.1). When destructible point is deleted, its neighbors are reviewed to ensure that they are not destructible either. In the following, we propose first analysis of selected algorithms (Danieslson-4SED version [START_REF] Danielsson | Euclidean distance mapping[END_REF], Cuisenaire-PSN version [START_REF] Cuisenaire | Fast Euclidean Distance Transformation by Propagation Using Multiple Neighborhoods[END_REF] and Meijster [START_REF] Meijster | A general algorithm for computing distance transforms in linear time[END_REF]) based on their 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. Results show that PSN algorithm is the most expensive in all cases (for any radius). Meijster algorithm is moderately faster than 4SED. Their curves are practically parallel and their returned values are proportional. However, difference between 4SED and PSN curves is more visible, it become larger when we increase radius. 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. . 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 with 0 x b n ∀ ≤ ≤ : ( , ) min : ( , ) /* G T */ 6. for ( ) Input : b:image, g: G_Table, m: columns, n:lines

T x x x x G p y p b b b y = - = , ( , ) min : ( , )
1 = x to ( ) 1 - n do 7. if [ ] B y x ∈ , then [ ] 0 .. = y x g 8. else [ ] [ ] 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.
[ ] for ( )

(( , ), ) B f x t q u = 9. while ( ) A B > then ( 1) q q ← + 10. end while 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. [ ] s q u ← 16. [ ] t q w ←
1 u m = -to ( ) 0 do 22. [ ] [ ] , (( , ), ) Edt x u f x u s q = 23. if [ ] ( ) u t q = then ( 1) q q ← -

endif

25.

Endfor

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.

= 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   ← ≥ Λ   [ ]

end Forall

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. Figure 35 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 fig. 36 (a).This can be explained by the choke point between columns processing and lines processing. Waiting time between these two treatments significantly penalizes acceleration.

Thinning and thickening computation

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 previous chapter. Inspired from Couprie approach, 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. ( , 1) push x stack

Global analysis

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 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 fig. 40 (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 'basic-NPS' scheduler, the speedup has increased to 5.2±0.01. Another common result between different architecture is stability of execution time on each ncore machine since the code uses n or more threads. For better readability of our results, we tested also efficiency of our algorithm on various architectures (see fig.

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 inflight 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 the '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 20). 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. The scheduler relies on our basic-NPS scheduling policy. As a result of this experiment, see fig. 41 (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. Moreover, the shape of the performance curve depends on how fast the cache hit rate degrades as a function of the number of threads. Any success access to L1 will eliminate an attempt to access to L2 thus performance curve, fig. 43 (A-1)(A-2), will evaluate in the same way. By reducing the number of cache miss from instruction cache, processor or thread of execution has not to wait (stall) until the instruction is fetched from main memory which immediately impact execution time. 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.

reliable and they also offer some flexibility of use in many application areas, particular image processing.

Based on this strategy, we proposed a series of parallel topological algorithm (new or adapted version). In the following we present our main contributions:

Major contribution concern a new algorithm to compute watershed that is parallel, preserves the topology of the input image, does not need prior minima extraction and suited for SMP machines. This algorithm does not require any sorting step, or the use of any hierarchical queue. Links between parallel watershed-cut and the SD&M strategy application can be described as follow: The splitting step is applied directly on input graph when selecting sources. Unlike conventional technique of division such as pixel division, or block division, the source selection is completely random. Associated steam computing is fully parallel (read mode data accesses).

Then distribution depends only on the available processors. This flexibility in data manipulation allowed us to obtain very good results especially in terms of efficiency without using the 'Basic-NPS' scheduler. Finally, the merging step contains procedures of s-labeling and f-labeling. Through these two functions, we have remained confident in our approach for merging streams two by two. Experimental analyzes such as execution time, performance enhancement, cache consumption, efficiency and scalability are also presented and discussed. Note that proposed parallel watershed-cut algorithm was preceded by an intensive study of different 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. For each approach, we give informal definition, then we presented processing procedure followed by mathematical foundations and the algorithm of reference.

Recent publications based on some approach are also presented and discussed. This study led us to classification of watershed algorithms according to criteria of recursion, complexity, basins computing and topology preservation.

Second contribution concern an adapted algorithm to compute skeleton that is parallel, preserves the topology of the input image and suited for SMP machines. Links between parallel watershed-cut and the SD&M strategy application can be described as follow: first step is dividing research area into different sub-region bounds. Since distribution start, each thread will lower each characterized pixel and then push its eight neighbors in available FIFO queue. Each queue contains the set of all selected neighbors and it is shared between only two threads.

Conventional synchronization techniques such that lock-based shared FIFO queue have not given good results and performance of our algorithm remained modest. Therefore we have applied our approach based on spin-wait FIFO queue for better performance. Dynamic lambda skeleton algorithm becomes five times faster than original version proposed by Couprie and al. [START_REF] Couprie | Topological operators for grayscale image processing[END_REF]. Tests on 2D grayscale image (512x512), using shared memory parallel machine (SMPM) with 8 CPUs cores (2 × Xeon E5405 running at frequency of 2 GHz), showed an enhancement of 6.2 with a maximum achieved cadency of 125 images/s using 8 threads. Note that proposed parallel lambda-skeleton algorithm was preceded by an intensive study of sixteen thinning algorithms in the frame work of critical kernels. We conclude this study by the classification of these algorithms according to five selection criteria: (i) preservation of topology, (ii) skeleton watershed transform, followed by classification of 3D watershed regions using the previously acquired bone mask. Third step, masks are deleted to cope with partial volume effect.

Even if this method is designed specifically for bone extraction (robustness and efficiency of the application is shown in fig. 45), Authors [START_REF] Straka | 3D Watershed Transform Combined with a Probabilistic Atlas for Medical Image Segmentation[END_REF] believe that it is applicable for a large class of software where object that might have variable densities throughout the dataset need to be identified. Other application on modern 3D medical imaging turns around cardiology software. In many clinical applications, it is essential to obtain precise information on the size and the function of the left ventricle (LV). Trough 3D images acquisition at different times of the heart cycle, Magnetic Resonance (MR) imagery offers a complete and perfect morphological left ventricle characterization. High precision in extracted measures make Magnetic Resonance imagery the method of reference in left ventricle study and analysis. But like any application with high precision in medical imaging, generated volume of data is very large. Thus, segmentation of such datasets is not only complex but also very expensive.

Based on this segmentation framework, Jean Cousty [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF][START_REF] Cousty | Automated, Accurate and Fast Segmentation of 4D cardiac MR images[END_REF] provides a new automated method to segment the left ventricular myocardium in 4D (3D+t) cine-MR images, see fig. 45. For object recognition, he used Exact Euclidean distance transforms to take into account prior geometric properties and Homotopic transforms to guarantee topological soundness of the segmentations. In order to assure time continuity of the successive 3D objects, he used a watershed transform that can be applied directly on the 4D sequence considered as a whole. To this end, he considered the watershed cuts in the watershed transform, and a temporal component of a 3D+t image gradient.

In the future, we plan to study this application and revisit used algorithms: replacing the watershed algorithm by the parallel watershed-cut to assess, in practice, the contribution of this alternative.

FIGURE

  FIGURE 1 : OVERVIEW OF PARALLEL COMPUTING STEPS...................................................................................................... FIGURE 2 : CIRCLE OF DECISION FOR THE ALGORITHM STRUCTURE DESIGN SPACE ............................................................. FIGURE 3 : SD&M STRATEGY DESIGN ................................................................................................................................ FIGURE 4 : AUTO-SUPPLYING TASK SYSTEM ....................................................................................................................... FIGURE 5: ILLUSTRATION OF TASK DISTRIBUTION ON MULTI-CORE ARCHITECTURE ............................................................ FIGURE 6: SPIN-WAIT SYNCHRONIZATION ........................................................................................................................... FIGURE 7: ILLUSTRATION OF MERGING PHASE WITH FOUR THREADS ................................................................................... FIGURE 8 : WATERSHED APPLICATION 1979-2003 ............................................................................................................... FIGURE 9 : STREAM NOTION ILLUSTRATION FOLLOWING COUSTY APPROACH [1]. ............................................................... FIGURE 10 : WATERSHED COMPUTING PRINCIPAL ................................................................................................................ FIGURE 11 : FLOW COMPUTE ILLUSTRATION ....................................................................................................................... FIGURE 12 : MERGING ILLUSTRATION ................................................................................................................................. FIGURE 13 : MERGING TECHNIQUES .................................................................................................................................... FIGURE 14 : (A) EXECUTION TIME (B) PERFORMANCE IMPROVEMENT [PARALLEL WATERSHED ALGO.] ............................... FIGURE 15 : CACHE PROFILING [PARALLEL WATERSHED] .................................................................................................... FIGURE 16 : EFFICIENCY IMPROVEMENT [PARALLEL WATERSHED ALGO.] ........................................................................... FIGURE 17 : PERFORMANCE IMPROVEMENT FOR PARALLEL WATERSHED ............................................................................ FIGURE 18 : SEGMENTATION CHAIN BASED ON PARALLEL WATERSHEDS-CUT ..................................................................... FIGURE 19: USED SHAPES FOR THINNING ALGORITHMS COMPARISON [79] .......................................................................... FIGURE 20: ALGORITHM CLASSIFICATION ACCORDING TO SKELETON CONNECTIVITY CRITERIA (C 1 ) .................................. FIGURE 21: ALGORITHM CLASSIFICATION ACCORDING TO SKELETON SYMMETRY CRITERIA (C 2 ) ....................................... FIGURE 22: EXECUTION TIME -SERIAL VERSIONS ON MONOCORE MACHINE [THINNING ALGO.] ......................................... FIGURE 23: CACHE PROFILING OF SERIAL VERSIONS ON MONOCORE MACHINE [THINNING ALGO.] ..................................... FIGURE 24: GUO AND HALL MASKS ..................................................................................................................................... FIGURE 25: JANG AND CHIN MASKS ..................................................................................................................................... FIGURE 26: BERNARD AND MANZANERA MASKS ................................................................................................................. FIGURE 27 : FILTERED SKELETON ILLUSTRATION [THINNING ALG.] .................................................................................... FIGURE 28 : EXECUTION TIME [PARALLEL THINNING] .......................................................................................................... FIGURE 29 : PERFORMANCE IMPROVEMENT [PARALLEL THINNING] ..................................................................................... FIGURE 30 : EFFICIENCY IMPROVEMENT [PARALLEL LAMBDA SKELETON] .......................................................................... FIGURE 31 : ILLUSTRATION ON DYNAMIC LAMBDA SKELETON PROCESS ............................................................................. FIGURE 32 : SMOOTHING ILLUSTRATION.............................................................................................................................. FIGURE 33 : OVERALL STRUCTURE [ORIGINAL SMOOTHING ALGORITHM] ........................................................................... FIGURE 34 : EXECUTION TIME [DANIESLSON, CUISENAIRE AND MEIJSTER ALGO.] ............................................................. FIGURE 35 : EVALUATION OF INSTRUCTION DISTRIBUTION (MEIJSTER ALG.) ...................................................................... FIGURE 36: (A) PERFORMANCE EVALUATION (B) EFFICIENCY EVALUATION [MEISJTER ALGO.] .......................................... FIGURE 37 : (A) (B) (C) (D) CRITICAL POINT ILLUSTRATION (E) RESEARCH AREA ASSIGNMENT .......................................... FIGURE 38 : (A) PERFORMANCE EVALUATION (B) EFFICIENCY EVALUATION [COUPRIE'S ALGORITHM]. ............................. FIGURE 39: TASKS DISTRIBUTION USING 'BASIC-NPS' [PARALLEL TOPOLOGICAL SMOOTHING] ......................................... FIGURE 40 : (A) GLOBAL PERFORMANCE IMPROVEMENT (B) GLOBAL EFFICIENCY IMPROVEMENT ...................................... FIGURE 41 : (A-1) INSTRUCTION -L1 MISSES; (A-2) ZOOM ON (A-1) [PARALLEL TOPOLOGICAL SMOOTHING] ..................... FIGURE 42 : (A) DATA READ (B) DATA WRITE -L1 MISSES [PARALLEL TOPOLOGICAL SMOOTHING] ................................ FIGURE 43 : (A-1) INSTRUCTION (B-1) DATA READ (B-1) DATA WRITE -L2 MISSES ......................................................... FIGURE 44 : ILLUSTRATION OF 3D WATERSHED TRANSFORM APPLICATION FOR MEDICAL IMAGE SEGMENTATION ......... FIGURE 45 : ILLUSTRATION OF 4D WATERSHED TRANSFORM APPLICATION FOR MEDICAL IMAGE SEGMENTATION .........

FIGURE 24 :

 24 FIGURE 1 : OVERVIEW OF PARALLEL COMPUTING STEPS...................................................................................................... FIGURE 2 : CIRCLE OF DECISION FOR THE ALGORITHM STRUCTURE DESIGN SPACE ............................................................. FIGURE 3 : SD&M STRATEGY DESIGN ................................................................................................................................ FIGURE 4 : AUTO-SUPPLYING TASK SYSTEM ....................................................................................................................... FIGURE 5: ILLUSTRATION OF TASK DISTRIBUTION ON MULTI-CORE ARCHITECTURE ............................................................ FIGURE 6: SPIN-WAIT SYNCHRONIZATION ........................................................................................................................... FIGURE 7: ILLUSTRATION OF MERGING PHASE WITH FOUR THREADS ................................................................................... FIGURE 8 : WATERSHED APPLICATION 1979-2003 ............................................................................................................... FIGURE 9 : STREAM NOTION ILLUSTRATION FOLLOWING COUSTY APPROACH [1]. ............................................................... FIGURE 10 : WATERSHED COMPUTING PRINCIPAL ................................................................................................................ FIGURE 11 : FLOW COMPUTE ILLUSTRATION ....................................................................................................................... FIGURE 12 : MERGING ILLUSTRATION ................................................................................................................................. FIGURE 13 : MERGING TECHNIQUES .................................................................................................................................... FIGURE 14 : (A) EXECUTION TIME (B) PERFORMANCE IMPROVEMENT [PARALLEL WATERSHED ALGO.] ............................... FIGURE 15 : CACHE PROFILING [PARALLEL WATERSHED] .................................................................................................... FIGURE 16 : EFFICIENCY IMPROVEMENT [PARALLEL WATERSHED ALGO.] ........................................................................... FIGURE 17 : PERFORMANCE IMPROVEMENT FOR PARALLEL WATERSHED ............................................................................ FIGURE 18 : SEGMENTATION CHAIN BASED ON PARALLEL WATERSHEDS-CUT ..................................................................... FIGURE 19: USED SHAPES FOR THINNING ALGORITHMS COMPARISON [79] .......................................................................... FIGURE 20: ALGORITHM CLASSIFICATION ACCORDING TO SKELETON CONNECTIVITY CRITERIA (C 1 ) .................................. FIGURE 21: ALGORITHM CLASSIFICATION ACCORDING TO SKELETON SYMMETRY CRITERIA (C 2 ) ....................................... FIGURE 22: EXECUTION TIME -SERIAL VERSIONS ON MONOCORE MACHINE [THINNING ALGO.] ......................................... FIGURE 23: CACHE PROFILING OF SERIAL VERSIONS ON MONOCORE MACHINE [THINNING ALGO.] ..................................... FIGURE 24: GUO AND HALL MASKS ..................................................................................................................................... FIGURE 25: JANG AND CHIN MASKS ..................................................................................................................................... FIGURE 26: BERNARD AND MANZANERA MASKS ................................................................................................................. FIGURE 27 : FILTERED SKELETON ILLUSTRATION [THINNING ALG.] .................................................................................... FIGURE 28 : EXECUTION TIME [PARALLEL THINNING] .......................................................................................................... FIGURE 29 : PERFORMANCE IMPROVEMENT [PARALLEL THINNING] ..................................................................................... FIGURE 30 : EFFICIENCY IMPROVEMENT [PARALLEL LAMBDA SKELETON] .......................................................................... FIGURE 31 : ILLUSTRATION ON DYNAMIC LAMBDA SKELETON PROCESS ............................................................................. FIGURE 32 : SMOOTHING ILLUSTRATION.............................................................................................................................. FIGURE 33 : OVERALL STRUCTURE [ORIGINAL SMOOTHING ALGORITHM] ........................................................................... FIGURE 34 : EXECUTION TIME [DANIESLSON, CUISENAIRE AND MEIJSTER ALGO.] ............................................................. FIGURE 35 : EVALUATION OF INSTRUCTION DISTRIBUTION (MEIJSTER ALG.) ...................................................................... FIGURE 36: (A) PERFORMANCE EVALUATION (B) EFFICIENCY EVALUATION [MEISJTER ALGO.] .......................................... FIGURE 37 : (A) (B) (C) (D) CRITICAL POINT ILLUSTRATION (E) RESEARCH AREA ASSIGNMENT .......................................... FIGURE 38 : (A) PERFORMANCE EVALUATION (B) EFFICIENCY EVALUATION [COUPRIE'S ALGORITHM]. ............................. FIGURE 39: TASKS DISTRIBUTION USING 'BASIC-NPS' [PARALLEL TOPOLOGICAL SMOOTHING] ......................................... FIGURE 40 : (A) GLOBAL PERFORMANCE IMPROVEMENT (B) GLOBAL EFFICIENCY IMPROVEMENT ...................................... FIGURE 41 : (A-1) INSTRUCTION -L1 MISSES; (A-2) ZOOM ON (A-1) [PARALLEL TOPOLOGICAL SMOOTHING] ..................... FIGURE 42 : (A) DATA READ (B) DATA WRITE -L1 MISSES [PARALLEL TOPOLOGICAL SMOOTHING] ................................ FIGURE 43 : (A-1) INSTRUCTION (B-1) DATA READ (B-1) DATA WRITE -L2 MISSES ......................................................... FIGURE 44 : ILLUSTRATION OF 3D WATERSHED TRANSFORM APPLICATION FOR MEDICAL IMAGE SEGMENTATION ......... FIGURE 45 : ILLUSTRATION OF 4D WATERSHED TRANSFORM APPLICATION FOR MEDICAL IMAGE SEGMENTATION .........

95 INTRODUCTION 1 . 1

 9511 ALGORITHM : THINNING ALGORITHM [JAN AND CHIN]---------------------------------------------------------------------------77 ALGORITHM : THINNING ALGORITHM [BERNARD AND MANZANERA] ----------------------------------------------------------77 ALGORITHM : DYNAMICALLY PARALLEL Λ -SKELETON [MAHMOUDI AND AKIL] --------------------------------------------80 ALGORITHM : E.D.T ALGORITHM -1 ST STEP -ORIGINAL VERSION [MEIJSTER] -----------------------------------------------90 ALGORITHM : E.D.T ALGORITHM -2ND STEP -ORIGINAL VERSION [MEIJSTER] ----------------------------------------------92 ALGORITHM : E.D.T ALGORITHM -1ST STEP -PARALLEL VERSION [MAHMOUDI AND AKIL] -------------------------------93 ALGORITHM : E.D.T ALGORITHM -2ND STEP -PARALLEL VERSION [MAHMOUDI AND AKIL] ------------------------------93 ALGORITHM : ADAPTED VERSION THINNING ALGO. [MAHMOUDI AND AKIL] ------------------------------------------------Context and motivations

Figure 1 :

 1 Figure 1 : Overview of parallel computing steps

Figure 2 :

 2 Figure 2 : Circle of decision for the algorithm structure design space

Figure 3 :

 3 Figure 3 : SD&M Strategy Design

  3). Therefore, Split Distribute and Merge strategy, that we propose, combines two patterns: Divide and Conquer pattern (Def.6) and Event-Based Coordination (Def.7). In the following we detail all three phases of SDM strategy. Divide and Conquer pattern It is based on multi-branched recursion. It solve problem by recursively dividing it into subproblems. After solving each sub-problem independently, it recombines the sub-solutions into a solution to the original problem.

α

  percentage of the serial part of the program using one processor.Scaled percentage of the serial part of the program using N CPUs

  Scaled percentage of the parallel part of the program using one processor.

  an algorithm runs on a first architecture using N processors with N η efficiency. Shared amount of work is N W . The same program runs on a second architecture using ' N processor with ' N η efficiency. Shared amount of work is ' N

Figure 4 :

 4 Figure 4 : Auto-Supplying task system

Figure 5 :

 5 Figure 5: Illustration of task distribution on multi-core architecture

Figure 6 :

 6 Figure 6: Spin-wait synchronization

Figure 7 :

 7 Figure 7: Illustration of merging phase with four threads

( a )

 a Cleavage fractures in steel, (b) contour of (a) obtained truth watershed definition introduced by Beusher and al. [56] 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.[59] in 2003.

Figure 8 :

 8 Figure 8 : Watershed application 1979-2003

Algorithm 4 :

 4 Flooding watershed [Vincent & Soille] input : Digital grey scale image G=(D,E,im) Output: Labeled watershed image lab on D 1. #define INIT -1 //initial value of lab image 2. #define MASK -2 //initial value of each level 3. #define WSHED 0 //label of the watershed pixels 4. #define FICTITIOUS (-1,1) //fictitious pixel ∉ D is the current label 6. fifo_init(queue) 7. for all (P∈D) do 8.

Algorithm 5 :

 5 Watershed by topographic distance process1. Foreach (marked area ∈ M

Algorithm 6 :

 6 Watershed by topographic distance [Meyer]Input : Lower complete image im on a digital grey scale image G=(D,E) with cost Output : Labeled 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.

Algorithm 8 : 5 .

 85 if (p 1) is already a marker pixel. And ( ) i I p refer to the value of the pixel i p in the image I .Definition (25): Path-cost functions [Lotufo and al.] IFT -Watershed from markers [Lotufo]Input : I : input image, wshed : labeled marker image Output : wshed : watershed catchment basins Aux : C : cost map, initialized to infinity; FIFO :hierarchical FIFO queue 1. while FIFO do //Propagation 6.

Algorithm 14 :

 14 Toboggan Algorithm [Lin and al.]Input : Img : a gradient magnitude image;Output: L : a label image, Q : empty FIFO queue ; 1.

  Q φ ≠ do //Simulation of keep-sliding for all C 2 pixels

  y ∈ do ψ(y) nb_labs ← ; 8.

  ): (a) the red graphs superimposed are the minima of corresponding functions. Let us consider G and F as associated graph and depicted function, (b) the sets { } , , , L a b e i = and { } , , j m n are two examples of streams, (c) the set stream since there is no path in ' L , between iand k , with steepest descent for F . Note that the sets { } , a b and { } b are respectively the set of bottoms and tops of L . Here the sets L is under the stream { } , , j m n and thus { } , , , , , , a b e i j m n is also a stream. There is no stream under{ } , , , a b e i and { } , , , , , , a b e i j m n . They are considered as two streams -≺ and they contain the set{ } , a b which is the vertex set of minimum of F .

Figure 9 :

 9 Figure 9 : Stream notion illustration following Cousty approach [1].

( a )Figure 10 :

 a10 Figure 10 : watershed computing principal

  Each flow 'L i ' is composed of point not yet labeled and whose source is x i . Stream function proposed by Cousty, pleaded in line 5 (Alg. 18), is launched N times. It allows the extraction of L i∈{1,2,3…n} . Intuitively, it explores the path of greatest slope, by mixing iterations first in-depth and width of the approaches. The main invariants of this function are: (i) The set 'L' is, for each iteration, a stream (flow). (ii) The set L' (line 2 -stream function) includes all wells of L not yet explored.

( a ) 11 :

 a11 Partition of input image (b) Parallel stream computation Figure Flow compute illustration At the end of flow function executing, a family , of N streams (L 1 , L 2 … L n ) whose elements must be labeled is generated. The initial procedure [in the iterative case] is to assign a new label (nb_labs) to each 'L i ' element if the latter is a stream -≺ . If it is not the case, the old returned label lab, of the stream -≺ 'H i ', included in 'L i ', is assigned to the different elements of 'L i '. Now if we want to launch this procedure in parallel manner, N/2 flows can be treated at ones. Algorithm 18 : Parallel watershed-cut [Mahmoudi and Akil.] Input : (V, E, F) : Edge-weighted graphs; Output : Ψ : Flow partition of F 1. foreach x ∈ V do Ψ(x) ← No-Label ; // No data dependency -FULL PARALALISM 2. nb_labs ← 0 ; // Global shared attributed label 3. i ←0 ; treated stream 4. foreach (x ∈ V) such as (Ψ(x) = = No-Label) do // lunch N process in parallel 5.

  see (fig 12.b). In this case, merging is more complicated. We developed f-labeling procedure (Algo. 19) able to make fusion in the following special cases: (i) 'L i ' and 'L i +1 'are two streams -≺ , (ii) 'L i ' and 'L i +1 'are two streams including two streams -≺ , (iii) 'L i ' is an stream -≺ and 'L i +1 ' is a stream including an stream -≺ . (a) Merging streams with common wells (b) Merging streams with common wells Figure 12 : Merging illustration

Figure 13 :

 13 Figure 13 : Merging techniques

Algorithm 20 : 6 .

 206 Function f-labeling [Mahmoudi and Akil] Input: (L a , lab a , L b , lab b , nb_labs) 1. // L a AND L b ARE two stream -≺ 2. if (laba = -1) && (labb = -1) then 3. nb_labs ++ ; 4. Attrib_lab(L a ,L b ,nb_lab) ; 5. // L a OR L b INCLUDES an stream -else if ((lab a ≠ -1) && (lab b = -1)) 7. Attrib_lab(L a ,L b ,lab a ) ; 8. else if ((lab a = -1) && (lab b ≠ -1AVspeed(L a )> AVspeed(L b )) 12. Attrib_lab(L a ,L b ,lab a ) ; 13. else Attrib_lab(L a ,L b ,lab b ) ; 14. Return NULL 15. Function Attrib_lab(L 1 ,L 2 ,lab) : 16. foreach (z ∈ {L 1 ∩ L 2 }) do Ψ(z) ← lab ; 17. foreach (x ∈ L 1 ) such as (Ψ(x) = = No-Label) do Ψ(x) ←lab; 18. foreach (y ∈ L 2 ) such as (Ψ(y) = = No-Label) do Ψ(y) ←lab;

Figure 14 :

 14 Figure 14 : (a) execution time (b) performance improvement [parallel watershed Algo.]

  lost time during memory access, we simply multiply the wait status by P4 660 frequency (3.6 GHz) and E5405 frequency (2 GHz). Thus we realize that estimated lost time on 8CPU is insignificant compared to lost time on 1CPU. This result is very visible on the E8400 and E5335 architectures. For the E5405 architecture, result is less visible due to cache structure: While E5405 is considered as eight CPUs architectural, but physically they are two Quads on the same chip (L2 = 2x4Mb).

Figure 15 :

 15 Figure 15 : Cache profiling [parallel watershed]

Figure 16 :

 16 Figure 16 : Efficiency improvement [parallel watershed algo.]

Figure 17 :

 17 Figure 17 : Performance improvement for parallel watershed

Figure 18 :

 18 Figure 18 : Segmentation chain based on parallel watersheds-cut

Figure 19 :

 19 Figure 19: Used shapes for thinning algorithms comparison [89]

Figure 20 :

 20 Figure 20: Algorithm classification according to Skeleton connectivity criteria (C 1 )

Figure 21 :

 21 Figure 21: Algorithm classification according to Skeleton symmetry criteria (C 2 )

Figure 22 :

 22 Figure 22: Execution time -serial versions on monocore machine [Thinning Algo.]

( 1 )Figure 23 :

 123 Figure 23: Cache profiling of serial versions on monocore machine [Thinning Algo.]

  G H x ) is said delatable if and only if:(i) D( G H x )=1, (ii) G( G H x )=0, (iii) B( G H x )>2 and(iv) Neighborhood of ( G H x ) does not match any the following masks.

Figure 24 :

 24 Figure 24: Guo and Hall masks Algorithm 21 : Thinning algorithm -Second version [Guo and Hall]

4 .

 4 Until Y φ = Couprie and al., see (Algo. 23), introduced new notions and operators in the frame work of cross-section topology. In particular, the notion of destructible λpixels which allows to selectively simplify the topology, based on a local contrast parameter λ .

Algorithm 23 :=

 23 Thinning algorithm [Couprie, Bezerra and Bertrand ]Jan and chin, see (Algo. 24), use the medial axis (of the input object) when computing Skeleton to enhance connectivity, unit-width convergence, medial axis approximation, noise immunity, and efficiency.

Figure 25 :

 25 Figure 25: Jang and Chin masks Algorithm 24 : Thinning algorithm [Jan and Chin]

Y=Figure 26 :

 26 Figure 26: Bernard and Manzanera masks Algorithm 25 : Thinning algorithm [Bernard and Manzanera]

  , we illustrate this method on a gradient image (a) obtained from a 2D grayscale image of an MRI brain section by Deriche gradient operator. (b) is obtained by a filtered thinning with Full algorithm is already introduced in last section (Algo. 23). (a) (b) (a) After Deriche gradient operator; (b) filtered skeleton with Figure 27 : Filtered Skeleton illustration [Thinning alg.]

( a )Figure 28 :

 a28 Figure 28 : Execution time [parallel thinning]

Figure 29 :

 29 Figure 29 : Performance improvement [parallel thinning]

்ሺଵሻ்ሺሻ.

  Then the defined efficiency Ψሺ݊ሻ ൌ ܶ ௦ ሺ݊ כ ⁄ ܶ ሻ can be written as Ψሺ݊ሻ ൌ ܶ ௦ ሺ݊ כ ⁄ be written as follows: Ψሺ݊ሻ ൌ ଵ כሺଵିሻା

Figure 30 :

 30 Figure 30 : Efficiency improvement [parallel lambda Skeleton]

  1 2 ( , ) x x with i x ∈ ℤ . Let x ∈Ε , r ∈ ℕ , we denote by ( ) r B x the ball of radius r S 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 fig. 32. Smoothed image (b) is obtained using Homotopic Alternative Sequential filter (HAS) with a radius equaled to five and four connectedness ( 4 Γ ).

  Input image (b) Smoothed image

Figure 32 :

 32 Figure 32 : Smoothing illustration

Figure 33 :

 33 Figure 33 : Overall structure [Original smoothing algorithm]

Figure 34

 34 Figure[START_REF] Chandra | Parallel Programming in OpenMP[END_REF] 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).

Figure 34 : 5 . 2 . 2

 34522 Figure 34 : Execution time [Danieslson, Cuisenaire and Meijster Algo.]

  , we scan each column y from top to bottom using the two first step can be written as follows. Complexity order is ( ) n m ο × .

Algorithm 27 :

 27 E.D.T algorithm -1 st Step -Original version [Meijster]Input : m:colums, n:lines, b:image

Figure 36 (

 36 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.

Figure 35 :Figure 36 :

 3536 Figure 35 : Evaluation of instruction distribution (Meijster Alg.)

Algorithm 31 :

 31 Adapted Version Thinning Algo. [Mahmoudi and Akil]

Figure 38 :

 38 Figure 38 : (a) Performance Evaluation (b) Efficiency Evaluation [Couprie's Algorithm].

  '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.

Figure 39 :

 39 Figure 39: Tasks distribution using 'Basic-NPS' [parallel topological smoothing]

  40 (b)) using the ( ) n ψ formula introduced earlier. For parallel time ratio we used best obtained time with 8 threads ('basic-NPS' scheduler).

Figure 40 :

 40 Figure 40 : (a) Global Performance improvement (b) Global efficiency improvement

Figure 41 :

 41 (A-1) Instruction -L1 misses; (A-2) zoom on (A-1) [parallel topological smoothing]

Figure 42 :

 42 (A) Data Read (B) Data Write -L1 misses [parallel topological smoothing]

Figure 43 : (A- 1 )

 431 Figure 43 : (A-1) Instruction (B-1) Data Read (B-1) Data Write -L2 misses

  Maximum Intensity Projection (MIP) of the original dataset, (b) MIP, bone tissue partially removed using mask obtained with probabilistic atlas. (c) MIP, bone tissue removed using mask extended with 3D watershed transform, (d) MIP, bone tissue removed using the dilated mask. (Patella and sacro-coccygeal bone not modeled).[59]

Figure 44 :

 44 Figure 44 : Illustration of 3D Watershed Transform application for Medical Image Segmentation

  Schematic view of three orthogonal sections of the objects of interest in left ventricle images (b) Three orthogonal section of 3D MR images of the left ventricle. (c) A three dimensional rendering of these objects of interest[START_REF] Cousty | Automated, Accurate and Fast Segmentation of 4D cardiac MR images[END_REF]

Figure 45 :

 45 Figure 45 : Illustration of 4D Watershed Transform application for Medical Image Segmentation

  1 : OVERVIEW OF PARALLEL COMPUTING STEPS...................................................................................................... FIGURE 2 : CIRCLE OF DECISION FOR THE ALGORITHM STRUCTURE DESIGN SPACE ............................................................. FIGURE 3 : SD&M STRATEGY DESIGN ................................................................................................................................ FIGURE 4 : AUTO-SUPPLYING TASK SYSTEM ....................................................................................................................... FIGURE 5: ILLUSTRATION OF TASK DISTRIBUTION ON MULTI-CORE ARCHITECTURE ............................................................ FIGURE 6: SPIN-WAIT SYNCHRONIZATION ........................................................................................................................... FIGURE 7: ILLUSTRATION OF MERGING PHASE WITH FOUR THREADS ................................................................................... FIGURE 8 : WATERSHED APPLICATION 1979-2003 ............................................................................................................... FIGURE 9 : STREAM NOTION ILLUSTRATION FOLLOWING COUSTY APPROACH [1]. ............................................................... FIGURE 10 : WATERSHED COMPUTING PRINCIPAL ................................................................................................................ FIGURE 11 : FLOW COMPUTE ILLUSTRATION ....................................................................................................................... FIGURE 12 : MERGING ILLUSTRATION ................................................................................................................................. FIGURE 13 : MERGING TECHNIQUES ....................................................................................................................................

FIGURE 14 : (A) EXECUTION TIME (B) PERFORMANCE IMPROVEMENT [PARALLEL WATERSHED ALGO.] ............................... FIGURE 15 : CACHE PROFILING [PARALLEL WATERSHED] .................................................................................................... FIGURE 16 : EFFICIENCY IMPROVEMENT [PARALLEL WATERSHED ALGO.] ........................................................................... FIGURE 17 : PERFORMANCE IMPROVEMENT FOR PARALLEL WATERSHED ............................................................................ FIGURE 18 : SEGMENTATION CHAIN BASED ON PARALLEL WATERSHEDS-CUT ..................................................................... FIGURE 19: USED SHAPES FOR THINNING ALGORITHMS COMPARISON [79] .......................................................................... FIGURE 20: ALGORITHM CLASSIFICATION ACCORDING TO SKELETON CONNECTIVITY CRITERIA (C 1 ) .................................. FIGURE 21: ALGORITHM CLASSIFICATION ACCORDING TO SKELETON SYMMETRY CRITERIA (C 2 ) ....................................... FIGURE 22: EXECUTION TIME -SERIAL VERSIONS ON MONOCORE MACHINE [THINNING ALGO.] ......................................... FIGURE 23: CACHE PROFILING OF SERIAL VERSIONS ON MONOCORE MACHINE [THINNING ALGO.] .....................................

  GLOBAL PERFORMANCE IMPROVEMENT (B) GLOBAL EFFICIENCY IMPROVEMENT ...................................... FIGURE 41 : (A-1) INSTRUCTION -L1 MISSES; (A-2) ZOOM ON (A-1) [PARALLEL TOPOLOGICAL SMOOTHING] ..................... FIGURE 42 : (A) DATA READ (B) DATA WRITE -L1 MISSES [PARALLEL TOPOLOGICAL SMOOTHING] ................................ FIGURE 43 : (A-1) INSTRUCTION (B-1) DATA READ (B-1) DATA WRITE -L2 MISSES ......................................................... FIGURE 44 : ILLUSTRATION OF 3D WATERSHED TRANSFORM APPLICATION FOR MEDICAL IMAGE SEGMENTATION .........

..... FIGURE 30 : EFFICIENCY IMPROVEMENT [PARALLEL LAMBDA SKELETON] .......................................................................... FIGURE 31 : ILLUSTRATION ON DYNAMIC LAMBDA SKELETON PROCESS ............................................................................. FIGURE 32 : SMOOTHING ILLUSTRATION.............................................................................................................................. FIGURE 33 : OVERALL STRUCTURE [ORIGINAL SMOOTHING ALGORITHM] ........................................................................... FIGURE 34 : EXECUTION TIME [DANIESLSON, CUISENAIRE AND MEIJSTER ALGO.] ............................................................. FIGURE 35 : EVALUATION OF INSTRUCTION DISTRIBUTION (MEIJSTER ALG.) ...................................................................... FIGURE 36: (A) PERFORMANCE EVALUATION (B) EFFICIENCY EVALUATION [MEISJTER ALGO.] .......................................... FIGURE 37 : (A) (B) (C) (D) CRITICAL POINT ILLUSTRATION (E) RESEARCH AREA ASSIGNMENT .......................................... FIGURE 38 : (A) PERFORMANCE EVALUATION (B) EFFICIENCY EVALUATION [COUPRIE'S ALGORITHM]. ............................. FIGURE 39: TASKS DISTRIBUTION USING 'BASIC-NPS' [PARALLEL TOPOLOGICAL SMOOTHING] ......................................... FIGURE 40 : (A) FIGURE 45 : ILLUSTRATION OF 4D WATERSHED TRANSFORM APPLICATION FOR MEDICAL IMAGE SEGMENTATION .........

TABLE 1 :

 1 PROCESSORS PERFORMANCE EVOLUTION ............................................................................................................

TABLE 2 :

 2 BASIC DEFINITIONS OF PROCESSING TIME ............................................................................................................

TABLE 3 :

 3 COMPARISON BETWEEN MAIN WATERSHED TRANSFORM .....................................................................................

TABLE 4 :

 4 TESTED IMAGE [PARALLEL WATERSHED] ............................................................................................................

TABLE 5 :

 5 USED PROCESSORS FEATURES [PARALLEL WATERSHED ALG.] .............................................................................

TABLE 6 :

 6 WALL CLOCK (MS) -[PARALLEL WATERSHED ALG.] ...........................................................................................

TABLE 7 :

 7 PERFORMANCE IMPROVEMENT [PARALLEL WATERSHED ALG.] ...........................................................................

TABLE 8 :

 8 CACHE PROFILING [PARALLEL WATERSHED] ......................................................................................................

TABLE 9 :

 9 EFFICIENCY IMPROVEMENT [PARALLEL WATERSHED ALGO.] .............................................................................

TABLE 10 :

 10 MAXIMUM AVERAGE UNIT SPEED [PARALLEL WATERSHED] ............................................................................

TABLE 11 :

 11 AVERAGE SPEED UNIT [PARALLEL WATERSHED ALGO.]...................................................................................

TABLE 12 :

 12 SCALABILITY PROFILING [PARALLEL WATERSHED ALGO.] ................................................................................

TABLE 13 :

 13 CLASSIFICATION OF THINNING ALGORITHM ACCORDING TO TOPOLOGY PRESERVATION ...................................

TABLE 14 :

 14 EVALUATION OF PIXELS' SKELETON (A I )

AND (N I ) [79] .....................................................................................

TABLE 15 :

 15 EVALUATION OF MEDIAL AXIS (REFERENCE) PIXELS (R FINDING CONCURRENCY DESIGN SPACE ..................................................................................................... DEFINITION (2): ALGORITHM DESIGN SPACE ........................................................................................................................ DEFINITION (3): ARCHITECTURE DESIGN SPACE ................................................................................................................... DEFINITION (4): PARALLEL IMPLEMENTATION MECHANISMS ................................................................................................. DEFINITION (5): PERFORMANCE METRICS OF PARALLEL PROGRAMS ...................................................................................... DEFINITION (6): DIVIDE AND CONQUER PATTERN ................................................................................................................. DEFINITION (7): EVENT-BASED COORDINATION .................................................................................................................... DEFINITION (8): SCALABILITY .............................................................................................................................................. DEFINITION (9): SPEEDUP-(1 ST EQUATION) ........................................................................................................................... DEFINITION (10): SPEEDUP-(AMDAHL APPROACH) ............................................................................................................... DEFINITION (11): SPEEDUP-(2 ND EQUATION) ........................................................................................................................ DEFINITION (12): SPEEDUP-(3 D EQUATION) ......................................................................................................................... DEFINITION (13): EFFICIENCY .............................................................................................................................................

I ) [79] ..........................................................................

TABLE 16 : TOP FIVE THINNING ALGORITHM ACCORDING TO SKELETON CONNECTIVITY AND SYMMETRY .......................... TABLE 17 : USED PROCESSORS FEATURES [THINNING ALG.]................................................................................................ TABLE 18 : CLASSIFICATION OF THINNING ALGORITHM ...................................................................................................... TABLE 19 : TIME EXECUTION RATE [SMOOTHING ALGORITHM] ........................................................................................... TABLE 20 : HARDWARE CONFIGURATION [PARALLEL SMOOTHING ALG.] ............................................................................ TABLE 21 : L2 -INSTRUCTIONS MISSES [TOPOLOGICAL SMOOTHING] ..............................................................................

LIST OF DEFINITIONS

DEFINITION (1): DEFINITION (14): PORTABILITY ............................................................................................................................................ DEFINITION (15): THRESHOLD SET OF (F) AT LEVEL (ALT) [VINCENT AND SOILLE] ................................................................ DEFINITION (16): FLOODING WATERSHED ............................................................................................................................ DEFINITION (17): THRESHOLD SET OF (F) AT LEVEL (ALT) [SHENGCAI AND LIXU] ................................................................. DEFINITION (18): SET OF ALL CONNECTED NEIGHBOR [SHENGCAI AND LIXU] ....................................................................... DEFINITION (19): LOWER SLOPE OF (F) AT A PIXEL P [ROERDINK AND AL.] ............................................................................ DEFINITION (20): COST FUNCTIONS [ROERDINK AND AL.] .................................................................................................... DEFINITION (21): TOPOGRAPHIC WATERSHED [ROERDINK AND AL.] ..................................................................................... DEFINITION (22): MAX-ARC PATH-COST FUNCTION [FALCÃO AND AL.] .................................................................................. DEFINITION (23): SPANNING FOREST [FALCÃO AND AL.] ....................................................................................................... DEFINITION (24): ARC WEIGHT FUNCTIONS .......................................................................................................................... DEFINITION (25): PATH-COST FUNCTIONS [LOTUFO AND AL.] ............................................................................................... DEFINITION (26): TOPOLOGICAL WATERSHED [COUPRIE AND BERTRAND] ............................................................................ DEFINITION (27): SIMPLE POINT [COUPRIE AND BERTRAND] ................................................................................................ DEFINITION (28): TOPOLOGICAL WATERSHED [COUPRIE AND BERTRAND] ............................................................................ DEFINITION (29): WATERSHED BASED ON LOCAL CONDITION [AUDIGIER AND LOTUFO] ........................................................ DEFINITION (30): MINIMAL ALTITUDE OF AN EDGE [COUSTY AND AL.] .................................................................................. DEFINITION (31): PROXIMITY OF THE SKELETON FROM MEDIAL AXIS (M M ) ............................................................................. DEFINITION (32): SKELETON CONNECTIVITY CRITERIA (C 1 ) ................................................................................................... DEFINITION (33): SKELETON SYMMETRY CRITERIA (C 2 ) .........................................................................................................

Pentium 4 Processor Extreme Edition Intel Core 2 Duo Processor E4300 Intel Xeon Processor X7560 Code name

		Prescott	Conroe	Nehalem-EX
	# of Cores	1	2	8
	# of Threads	2	2	16
	Clock Speed	3.73 GHz	1.8 GHz	2.266 GHz
	Cache	2 MB L2 Cache	2 MB L2 Cache	24 MB L3 Cache
	Bus Type	FSB	FSB	QPI
	System Bus	1066 MHz	800 MHz	6.4 GT/s

Table 1 : Processors performance evolution

 1

  Total processing time of the serial and parallel part of the program using one processor.

	t	T	) 1 (	=	t	s	+	t	p	(
	t	T	(	N	)		t	s		t	p	(	N	)

  's Law defines the maximum speed up S that can be achieved by n processors collaborating on an application, where k is the fraction of the work that can be executed in parallel. Assume, for simplicity, that a single processor completes the work in one second. With n concurrent

	processors, the parallel part takes ( ) p n	seconds and the sequential part takes (1 ) p -	seconds.
	Overall, the parallelized computation takes ((1 ) p -+	p n	)	seconds. Thus, ( ) 1/ ((1 ) S n p = -+ . ) p n

  3), the filling of basins is an iterative process that involves gradually raising the water level from Alt min to Alt max . Algorithm must, foreach iteration, fill existing basins (extension regions) and possibly create new basins (new regions). L R will denote region list. Dams will be built where waters coming from different basins meet.This threshold set define a recursion with a 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 . Thus we can introduce the following definition. MIN h will denote the union of all regional minima at altitude Alt.

	For mathematical formulation of the mentioned process, let	f :	D	→	Ν	be a digital grey value
	Alt T	p D f p Alt = ∈ ≤
	Algorithm 3 : Flooding watershed process
	1. for level from Alt min to Alt max
	2.		// Action 1 : Extend existing region
	3.		foreach (R∈ L R ) 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 L R ;
	9.		Add Pixel P to region [R];
	10.		Growing [R] until level ;
	Definition (16): Flooding watershed
		min	{	min	}	min

image, with Alt min and Alt max the minimum and the maximum value of. The threshold set of f at level Alt can be defined as: Definition (

15

): Threshold set of (f) at level (Alt)

[Vincent and Soille]

{ } / ( )

  is a work image of distances 10. end_for 11. SORT pixels in increasing order of grey values ( h min , h max ) // starting flooding process 12. for h = h min to h max do //geodesic SKIZ of level h-1 inside level h 13.

			lab	p		←	MASK	;
	15.		if (p has a neighbour q) with ((lab[p] > 0 or lab[q]=WSHED)) then
				//initialize queue with neighbours at level h of current basins or watersheds
	16.			dist	[ ] ; 1 ← p
	17.				fifo	_	add	(	q	,	queue	);
	18.		end_if
	19.	end_for
	20.	curdist	←	; 1
	21.	fifo	_	add	(	FICTITOUS	,	); queue
	22.	loop //extend basins
	23.	p ←	fifo	_	( remove	); queue
	24.	if (p = FICTITIOUS) then
	25.		if	(	fifo	_	empty	(	queue	))	then

for all (P∈D) with im[p]=h do // mask all pixels at level h //these are directely accessible because of the sorting step 14.

[ ]

  Assign label to each unlabeled image by first tobogganing then backtracking using best first search.

	Lin and al [77] define the D with an altitude strictly greater than C refers to all pixels p in h following pixels classes: Class 1 the altitude of its lowest neighbor. Class 2 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 1. Records the sliding directions for all 1 2 ( ) ( ) p C p C ∈ ∪ ∈ 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 ∈ in Algorithm 13 : Toboggan watershed process 3.

h D belonging to a connected component without any catchment basin. Thus we can give a global description of the computing process (Algo. 13) followed by the toboggan algorithm (Algo. 14).

End For Algorithm 15 : Resolve function [Lin and al.]

	33.					If ( ) L q is not assigned then
	34.						0 ( )) L q L p ( ( ) =
	35.						Q	←	p
	36.					End if
	37.				End for
	38.		End while
	39. End if
	40. End for
	41. For all (	) ∈ do // Tobogganing -Depth first search p D
	42. Resolve( ) p
	43. Input : Pixel site p
	1.	If ( ) L p is not assigned then
	2.
	24.			End if
	25.	End while
	26. For all	(	0 p D ∈ and )	(	0 ( ) SlidingList p	) = do // labeling C 3 pixels φ
	27.	If	0 ( ) L p is not assigned then
	28.		0 ( ) L p =new_label
	29.		0 ( ) h G p =
	30.		While Q φ ≠ do
	31.			p ←	Q
	32.			For all (	( )) p Neighbor p ∈	and ( ( ) G q h = do )

  1.2. Now, Let consider a stream L , we say that x L ∈ is a top of L if the altitude of x is greater than, or equal to the altitude of any y L ∈ . If the altitude of x is less than the altitude of any y , then x is considered as a bottom of L .

	Let consider two disjoint streams 1 L and 2 L , with 1 L L φ ∩ = . Let L be the union of both 2 streams with 1 2 L L L = ∪ .We say that 1 L is under 2 L , written 1 L L ≺ , if there exist a top x of 2 1 L and a bottom y of 2

Algorithm 19 : Function s-labeling [Mahmoudi and Akil]

	Input: (L, lab, nb_labs)
	1. if (lab = -1) then // L is	≺	-	stream
	2.	nb_labs ++ ;
	3.	foreach (y ∈ L) do Ψ(y) ← nb_labs ;
	4. else
	5.
					A	B	X1 X2	X8
					L1
					L2
					L8
					A	B

this case the merging is simple, for each flow 'L i ' and 'L i+1 ', see (fig 12.a). Note that s-labeling function (Algo. 19) launches only the initial procedure [used previously in the iterative case].

Table 4 : Tested Image [parallel watershed]

 4

	TESTED IMAGES

Table 5 : Used processors features [parallel watershed alg.]

 5 Wall-clock execution times for numbers of threads equal to 1, 2, 4, 8, 16 and 32 were determined. Different characteristics of used architectures are presented in table 5. The minimum value of 2 timings was taken as most indicative of algorithm speed. Results of implementation on the different architecture are shown in the following table.

		1 CPU	2 CPUs	4 CPUs	8 CPUs
	1 Thread	4638	4448	4898	5190
	2 Threads	5321	3182	3114	3092
	4 Threads	4898	3303	1384	1709
	8 Threads	5253	3253	1639	713
	16 Threads	5129	3278	1744	990
	32 Threads	5190	3303	1794	1235

Table 6 : wall clock (ms) -[parallel watershed Alg.]

 6

Table 7 : Performance improvement [parallel watershed Alg.]

 7

Table 8 : Cache profiling [parallel watershed]

 8

			Nbr. Instr.	D1 miss	L2d miss	WScm
		1 Thread	34.598.772	194.274	141.738	14 699
	1CPU	2 Threads 4 Threads	34.204.145 34.340.721	235.764 268.779	162.775 168.539	17 007 17 856
		8 Threads	34.441.168	288.664	161.337	17 406
		1 Thread	90.783.715	221.199	198.110	20 041
	2CPUs	2 Threads 4 Threads	90.875.188 116.984.996	229.857 284.697	198.282 207.448	20 143 21 517
		8 Threads	152.704.881	284.753	207.448	21 517
		1 Thread	334.816.008	251.241	215.443	21 902
	4CPUs	2 Threads 4 Threads	339.998.650 334.020.732	251.982 204.315	215.582 215.860	21 922 21 470
		8 Threads	474.895.119	295.774	228.187	23 494
		1 Thread	784.648.688	265.884	216.760	22 167
	8CPUs	2 Threads 4 Threads	784.745.461 789.432.158	271.859 279.951	219.487 224.625	22 472 23 015
		8 Threads	790.804.849	282.122	231.142	23 624

table 9

 9

		Intel P4-660	Intel Dual C. E8400	Intel C2 Quad E5335	Intel Xeon E5405
	EFFICIENCY RATE	0,94	0,69	0,79	0,76

, see also fig.

16.

Table 9 : Efficiency improvement [parallel watershed Algo.]

 9

Table 10 : Maximum Average Unit Speed [parallel watershed]

 10

	1
	0,8
	0,6
	0,4
	0,2
	0
	1 CPU	2 CPUs	4 CPUs	8 CPUs

Table 11 ,

 11 to reach Average Speed Unit us Av . Note that chosen Average Unit Speed is 0,787 on 4 CPUs using 4 Threads (Associated W = 334.020.732).

			1 CPU	2 CPUs
	1 Thread	us Av	0,940	0,787	-
		Work (W')	34.598.772	28.967.270
	2 Threads	us Av Work (W')	-	0,685 90.875.188	0,787 104.406.968
	4 Threads	us Av Work (W')	-		-
	8 Threads	us Av Work (W')	-		-
			4 CPUs	8 CPUs
	1 Thread	us Av	-		-
		Work (W')
	2 Threads	us Av Work (W')	-		-
	4 Threads	us Av Work (W')	0,787 334.020.732	-
	8 Threads	Average Unit Speed us Av	-	0,763 790.804.849	0,787 815.679.444

Table 11 : Average speed Unit [parallel watershed Algo.]

 11 The scalability results for parallel watershed processing are shown in table 12. Our experiments demonstrate a very good scalability across all tested architectures.

		1 CPU	2 CPUs	4 CPUs	8 CPUs
	1 CPU	1	0,554	0,338	0,284
	2 CPUs		1	0,626	0,512
	4 CPUs			1	0,819
	8 CPUs				1

Table 12 : Scalability profiling [parallel watershed Algo.]

 12

Table 14 : Evaluation of pixels' Skeleton (A i ) and (N i ) [89]

 14

:

  Skeleton symmetry criteria (C 2 )A second assessment using same images showed in fig.19, has allowed us to draw second curve (see fig.21). Measurements are presented in table 14 and table15.

	C	2	=	( ) ( ) i i Area A Area R
					Shape (1)	Shape (2)	Shape (3)
	Medial axis (reference)	564	1359	2178

Table 15 : Evaluation of Medial axis (reference) pixels (R i ) [89]

 15

Table 16 : Top five thinning algorithm according to Skeleton connectivity and symmetry

 16

	According to fig. 20 and fig. 21, first algorithm proposed by Pavlidis [87] has a very good factor
	(C 2

Table 17 : Used processors features [Thinning alg.]

 17

Table 18 : Classification of thinning algorithm

 18 Guo and Hall, see (Algo. 21), introduced the notion of delatable pixel G Hx . Based on Boolean expressions defined bellow, they classify a pixel as delatable if some conditions are hold: Let

	1	Guo, Hall 1992 (b)
	2	Eckhardt, Maderlechner 1993
	3	Couprie and al. 2000
	4	Jang, Chin 1993
	5	Bernard, Manzanera 1999

pixel. (iii) Inner boundary pixel: it

  is a boundary pixel which has an interior pixel as strong neighbor.Thus, a pixel ( ) p in X can be classified as simple, if it is a boundary pixel and if there exist exactly one strong connected component of pixels of X in the neighborhood of ( )

								p which is
	strongly connected to ( ) p .
	Finally, inner boundary pixel ( ) p is classified as perfect, if there exists a strong neighbor
	Γ	( ) i p	of ( ) p which is interior and such that ( ) j p Γ	∉	X	with	( 4) mod 8 j i = +	.
	Algorithm 22 : Thinning algorithm [Eckhardt and Maderlechner]
	Input: set X
	Output: set X
		1. repeat
		2.
		4. Until Y φ

=

Eckhardt and Maderlechner, see (Algo. 22), introduced the notion of simple and perfect pixels. First, they establish the following definitions: (i) Interior pixel: it is a pixel in X having all its four strong neighbors in X (ii) Boundary pixel: it is a pixel in X which is not interior

  Tሺx, Xሻ and Tሺx, Xሻ are the two connectivity numbers defined as follows (# ܺstands for the cardinal of ܺ): ܶሺ,ݔ ܺሻ ൌ ܥ# ሾ,ݔ Γ ଼ כ ሺݔሻ ת ܺሿ; ܶሺ,ݔ ܺሻ ൌ ܥ# ሾ,ݔ Γ ଼ כ ሺݔሻ ת ܺሿ;

	So we can define the four neighborhoods:
	߁ ାା ሺ,ݔ ܨሻ ൌ ሼݕ א ߁ ଼ כ ሺݔሻ; ܨሺݕሻ ܨሺݔሻሽ ߁ ା ሺ,ݔ ܨሻ ൌ ሼݕ א ߁ ଼ כ ሺݔሻ; ܨሺݕሻ ܨሺݔሻሽ ߁ ିି ሺ,ݔ ܨሻ ൌ ሼݕ א ߁ ଼ כ ሺݔሻ; ܨሺݕሻ ൏ ܨሺݔሻሽ ߙ ି ሺ,ݔ ܨሻ ൌ ൜ ݉ܽݔሼܨሺݕሻ, ݕ א ߁ ିି ሺ,ݔ ܨሻ, ݂݅߁ ିି ሺ,ݔ ܨሻ ് ሽ ܨሺݔሻ ݁ݏ݅ݓݎ݄݁ݐ

Algorithm 26 : Dynamically Parallel ૃ -Skeleton [Mahmoudi and Akil] Input : m:colums, n:lines, b:image

	70
	60 1. For all p א Ime do 2. if ൫m ୧୬ ൏ Ime ୧ ሺpሻ ൏ m ୱ୳୮ ൯then E ୧ ՚ E ୧  ሼpሽ; 3. Repeat until stability 4. E ୱ୬ ՚ ; 5. While ሺk ് 0ሻthen 6. For all p א E ୧ do 7. if ሺp א Tሻ then F ୧ ሺxሻ ՚ α ି ൫x, F ୧ ሺxሻ൯; 8. 20 30 40 50 1 cores Wall-clock time [ms] 2 cores E ୱ୬ ՚ E ୱ୬  ሼeight p neighborsሽ; 9. 10. endif 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. end while clean ሼE ୧ , E ୩୧ , E ୩୧ାଵ ሽ; if ሺE ୧ ൌ ሻ then k ՚ 0; Ime ୧ ՚ Ime ୧  E ୩୧  E ୩୧ାଵ ; Ime ୧ ՚ Ime ୧ Ime ୧ାଵ ; E ୧ ՚ E ୱ୬ ; endif else E ୩୧ାଵ ՚ E ୩୧ାଵ  ሼpሽ; E ୱ୬ ՚ E ୱ୬  ሼeight p neighborsሽ; if ሺp א Tሻ then Ime ୧ାଵ ሺxሻ ՚ α ି ൫x, Ime ୧ାଵ ሺxሻ൯; Number of threads For all p א E ୧ାଵ do 0 1 2 4 8 16 20 0 4 cores 10 else E ୩୧ ՚ E ୩୧  ሼpሽ; 8 cores

algorithm -2nd Step -Original version [Meijster]

	Algorithm 28 : E.D.T
			g	x	,	y	=	g	x	+	, 1	y	+
	9.	endif
	10.	endfor
	11. 12.	/* G B */ for ( = n x	-		2	)	downto ( ) 0 do
	13.	if [ x g	+	] [ ] y x g , 1 < y ,	then
	14.	g	[ ] [ , = g y x	x	+	, 1	y	] 1 +
	15.	endif
	16.	Endfor; endforall

Algorithm 29 : E.D.T algorithm -1st Step -parallel version [Mahmoudi and Akil] 1. For

		(	, y t y m y y t , = < = +	max	)	do
	2.	If ( ) B y ∈ , 0	then [ ] 0, g y ←	0
	3.					else [ ] 0, g y ← ∞
	4.	endif
	5. 6.	/* G T */ for ( ) 1 = x	to ( ) 1 -n	do
	7.		if [ ] B y x ∈ ,	then [ ] , g x y ←	0
	8.		else [ ] [ , g x y g x y 1, ← +	]	+	1
	9.		endif
	10.	Endfor
	11. 12.	/* G B */ for ( = n x	-	2	)	downto ( ) 0 do
	13.		if	(	[ g x y g x y ] [ ] 1, , ) + <	then
	14.			[ ] [ , g x y g x y 1, ← +	]	+	1
	15.		endif
	16.	endfor
	17.

Endforall

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 *p. time) (ii), already presented in section 2.4.1. 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 Table

20. Algorithm 30 : E.D.T algorithm -2nd Step -parallel version [Mahmoudi and Akil] 1. For

	(	, x t x n x x t , = < = +	max	)	do
	2.

  combine the parallel version of Meijster algorithm and the adapted version of Couprie's algorithm to build the parallel processing of topological smoothing.

	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 )
	15.			x	←	( pop stack	2)
	16.			( ) v neighbors x ←
	17.			i ←	0
	18.			While ( 8) i < do
	19.				if [ ] ( v i stack 1) ∉

then is to

Table 20 : Hardware configuration [parallel smoothing alg.]

 20

Table 21 : L2 -Instructions Misses [Topological smoothing]

 21

Last version : October 10,

Split Distribute & Merge

LIST OF FIGURES

Watershed based on

Vincent & Soille [START_REF] Vincent | Watersheds in digital spaces : An efficient algorithm based on immersion simulations[END_REF] Meyer [START_REF] Meyer | Topographic distance and watershed lines[END_REF] Lotufo [START_REF] Lotufo | The Ordered Queue and the Optimality of the Watershed Approaches[END_REF] Couprie [START_REF] Couprie | Topological grayscale watersheds transform[END_REF] Lin [START_REF] Lin | Comparison Between Immersion-Based and Toboggan-Based Watershed Image Segmentation[END_REF] Cousty [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF] Defined 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 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.

Trough this analysis, (Algo. 16) holds best characteristics. The fact that sorting step is not required, hierarchical queue is not used and minima are not computed, make it an excellent candidate for parallelization on shared memory architecture.

Construction of parallel topological watershed

In this section, we start by introducing some basic definitions of stream notion which is crucial to the flooding paradigm. Then, we introduce in detail our parallel watershed-cut. Illustration of parallel compution process is given. Execution time and cache consumption are performed and analyzed. Efficiency and scalability are also presented and discussed.

Conclusion

In this chapter, we have presented an intensive study of work of critical kernels: Pavlidis [ and Hall [START_REF] Guo | Fast fully parallel thinning algorithms[END_REF] (3 variants), Jang and Chin [ Manzanera [START_REF] Bernard | Improved low complexity fully parallel thinning algorithm[END_REF], Rutovitz [START_REF] Rutovitz | Pattern recognition[END_REF], Zhang and Wang [ Couprie and al. [START_REF] Couprie | Topological operators for grayscale image processing[END_REF].

First contribution is limited to the criteria: (i) preservation of topology execution time and (v) cache consumption algorithm as the most suitable algorithm Second contribution concern an adapted algorithm to compute the topology of the input image and suited for SD&M strategy and dynamic lambda different sub-region bounds. Since distribution pixel and then push its eight neighbors in all selected neighbors and it is sh synchronization techniques such that and performance of our algorithm remain based on spin-wait FIFO queue Dynamic lambda skeleton algorithm Couprie and al. [START_REF] Couprie | Topological operators for grayscale image processing[END_REF]. Tests on 2D machine (SMPM) with 8 CPUs an enhancement of 6.2 with a maximum achieved cadency of an adapted algorithm to compute skeleton that is parallel, preserves the topology of the input image and suited for SMP machines. Fig. 31 makes link between lambda skeleton algorithm: first step is dividing research area into Since distribution start, each thread will lower each characterized pixel and then push its eight neighbors in available FIFO queue. Each queue contains all selected neighbors and it is shared between only two threads. Unfortunately synchronization techniques such that lock-based shared FIFO queue have not given good results and performance of our algorithm remained modest. Therefore we have applied our approach wait FIFO queue, already introduced in section 2.4.

For general case, we define

. We say also that two points x and y of X are n-connected in X if there is an 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 nconnected 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 . To guarantee correspondence between X topology and X topology, we use n-adjacency for X and n -adjacency 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

∈ , two connectivity numbers defined as follows ( # X = cardinality of X ):

( )

The following properties allow us to locally characterize simple points [START_REF] Kong | Digital topology: introduction and survey[END_REF][START_REF] Bertrand | Simple points, topological numbers and geodesic neighborhoods in cubic grids[END_REF] hence to implement efficiently topology preserving operators:

A preliminary assessment of first implementation code, see

Study on Euclidean 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 only four neighbors (on element structure). 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 [9] 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].

Let's move to the second step. We start by defining 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".

Let S be the set of y points such that ( , ) f p y is minimal and unique. The formula of S ,

. 0 i u u m ∀ ≤ ≤ ∧ ≤ . Let T be the set of points with coordinate greater than, or equal to, horizontal coordinate of the intersection with a region:

Let ( , )

Sep i u be the separation between regions of i and u , defined by the following with 2 2

( ( , ) ( , ) )

( , ) ( ) / 2( 1)

.

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 . For each respective values of T , f is computed. Algorithm 28 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 SD&M parallelization strategy. In the first stage, column processing, we can define data interdependence by the following equation:

{ }

( , ) min ( , ), ( , )

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 ) ( )

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 = + .

20.

[ ] ( , 1) push v i stack

← -

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. 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 { }

, , 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 fig. 38 (a) (b). 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

Conclusion

In this chapter, 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.

Introduced smoothing filter is parallel, preserve the topology and suited for SMP implementation. SD&M strategy was applied twice. First time, when computation E.D.T, the splitting step starts by sharing columns and rows to scan between different processes. A thread can process one or more columns (rows). Number of threads will depend only on number of processors. Second time, when computing thinning (thickening), splitting step starts by dividing research area into different sub-region bounds. Since distribution start, each thread will lower each characterized pixel and then push its eight neighbors in available FIFO queue. Each queue contains the set of all selected neighbors and it is shared between only two threads. Unfortunately, obtained results using SMP scheduling policy, are not sufficient especially for cache consumption. For this reason, we move to PSN scheduler when distributing work. Finally we apply the same approach of fusion threads in pairs when computing E.D.T, thinning and thickening.

CONLUSION

his chapter presents the conclusions of the thesis, an extended summary of the research work and restatements of contributions. Future work section summarizes the next steps to follow into the research of parallelization strategy.

Contribution

The main idea behind the present thesis is basically to identify the best approach to parallelize image processing operators based on topological transform. This approach, formally described as a strategy, was designed for shared memory parallel machines.

Parallelization strategy is the science of planning and marshalling architecture resources for their most efficient and effective use. Its main objective is to provide a mechanism to design parallel algorithm by identifying existing concurrency, modifying their structure before any further implementation (with right techniques) takes place. The idea is to take into consideration algorithm characteristics and performance issues during parallel algorithm design from the beginning. So, introduced strategy (SD&M) aims to give a full description of parallel topological algorithm design approach, helping parallel designers to obtain powerful and enhanced computation method of topological transforms on shared memory parallel machines. If we consider a given topological operator, the parallelization strategy describes how such algorithm can be modified with respect of the following algorithmic characteristics: (i) Topology : structure allowing full topology preservation (ii) Processing: structure providing parallel processing functionalities and best work distribution. During parallel processing, structure assuring coordination and communication are also used.

Thus, using SD&M strategy, a powerful parallel operator based on topological transform will be the result of restructuring and splitting original amount of work between different processor. Work distribution and processing activities take into account coordinating and communication.

In more formal way, SD&M strategy can be described as the combination of divide and conquer patterns and event-based coordination patterns hence the name that we have assigned. Note that introduced strategy represents the last stitch in the decomposition chain of algorithm design patterns and it provides a fine-grained description of topological operators' parallelization. It covers recursive algorithms and it is especially designed for shared memory architecture with uniform access. Although the cost of communication (Memory-processor and inter-processors) is high enough, shared memory architectures appear to be among adapted platforms for this type of processing. Actually, these architectures have the advantage of allowing immediate sharing of data with is very helpful in the conception of any parallelization strategy. They are nondedicated architecture using standard component (processor, memory, buses ...) so economically T connectivity, (iii) skeleton symmetry, (iv) execution time and (v) cache consumption. Through this classification, we identified Couprie's algorithm as the most suitable algorithm for parallelization on shared-memory architectures.

Third contribution concern 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. Introduced smoothing filter is parallel, preserve the topology and suited for SMP implementation. SD&M strategy was applied twice. First time, when computing E.D.T, the splitting step starts by sharing columns and rows to scan between different processes. A thread can process one or more columns (rows). Number of threads will depend only on number of processors. Second time, when computing thinning (thickening), splitting step starts by dividing research area into different sub-region bounds. Since distribution start, each thread will lower each characterized pixel and then push its eight neighbors in available FIFO queue. Each queue contains the set of all selected neighbors and it is shared between only two threads. Unfortunately, obtained results using SMP scheduling policy, are not sufficient especially for cache consumption. For this reason, we move to PSN scheduler when distributing work. Finally we apply the same approach of fusion threads in pairs when computing E.D.T, thinning and thickening.

Perspectives

A possible extension of present work is the use of produced parallel algorithms and powerful multicore/multithread architectures to improve modern 3D medical imaging software in which segmentation procedure plays a crucial role.

Applications based on multiple detector-row computed tomography (CT) provide huge highresolution volumetric datasets which is extremely hard to inspect without computer aided preprocessing. Straka and al. [START_REF] Straka | 3D Watershed Transform Combined with a Probabilistic Atlas for Medical Image Segmentation[END_REF] work aimed at visualization and treatment planning of peripheral arterial occlusive disease by means of CT-angiography. In fact, for the visualization and the diagnostic assessment of vascular disease it is useful if tissues in images are segmented. To highlight data volume and complexity of such application, note that CT-angiography of the peripheral arteries is performed using multiple-detector row computed tomography. Thus obtained series contain about 2000 transverse images (512x512-12bit/pixels). The analysis of such data volume using currently medical workstation takes about four hours.

A first contribution to solve this problem is based on algorithmic solution proposed by Straka [START_REF] Straka | 3D Watershed Transform Combined with a Probabilistic Atlas for Medical Image Segmentation[END_REF]. The basic idea of his method, for bone tissue labeling in computed tomography data, is reflected in the combined application of two techniques: a-prior knowledge derived from a density based probabilistic atlas is used to locate characteristic parts of bones and a watershed transform to identify spatially coherent sub volumes, regardless of their density. First step consist a pre-computed atlas of bone density information to assign a bone-probability to each voxel. Second step consist on partition the whole volume in homogenous regions using the