R Mahmoudi

M Akil
email: akilm@esiee.fr

Image Segmentation Based Upon Topological Operators Real Time Implementation Case Study

Keywords: Crest Restoring, Thinning, Seg mentation, Hardware Imp lementation, Optimization Techniques, Parallelis m, Pentium, Profiling, SSE instructions set

In miscellaneous applications of image treatment, thinning and crest restoring present a lot of interests. Recommended algorith ms for these procedures are those able to act directly over grayscales images while preserving topology. But their strong consummation in term of time remains the major d isadvantage in their choice. In this paper we present an efficient hardware imp lementation on RISC processor of two powerful algorithms of thinning and crest restoring developed by our team. Proposed imp lementation enhances execution time. A chain of segmentation applied to medical imaging will serve as a concrete examp le to illustrate the improvements brought thanks to the optimization techniques in both algorith m and architectural levels. The part icular use of the SSE instruction set relative to the X86_32 processors (PIV 3.06 GHz) will allow a best performance for real time processing: a cadency of 33 images (512*512) per second is assured.

INTRODUCTION

In order to avoid important significant informat ion loss, due to the passage to binary image, nu merous operators acting directly over greyscale image were proposed [START_REF] Huang | An Improved Parallel Thinning Algorithm[END_REF] [6] [START_REF] Bertrand | Image segmentation through operators based upon topo logy[END_REF]. These kinds of operators that preserve the topology of the image are extremely useful in solving various problems such as medical image analysis, 3D model matching and character recognition. Different algorithms dealing with thinning and skeleton were studied [START_REF] Bertrand | Two-Dimensional Parallel Thinning Algorithms Based on Critical Kernels[END_REF] [4] [START_REF] Couprie | Note on fifteen 2D parallel thinning algorithms[END_REF] (in [START_REF] Couprie | Note on fifteen 2D parallel thinning algorithms[END_REF] fifteen 2D thinning algorith ms were co mpared). To grade these algorith ms, the majo rity of publications focus on the following criteria as quality criteria [START_REF] Canny | A computational approach to edge detection[END_REF] of classification i) Ho motopy: the skeleton must preserve the topology. ii) One-p ixelthickness: one-pixel-thick objects. iii) Noise immunity: Damage on skeleton caused by noise must be repaired. Other criteria relative to the imp lementation and the execution time can be introduced [START_REF] Couprie | Topological operators for greyscale image processing[END_REF] [5] [START_REF] Garcia Lorca | Efficient ASIC and FPGA imp lementation of IIR filters for real time edge detections[END_REF] to classify algorith ms as sequential (each non-significant information is detected and removed immediately), or parallel (all non -significant informat ion are firstly detected then removed). A mong all the proposed algorithms, we are interested in both operators of thinning and crest restoring introduced by Bertrand and Couprie in [START_REF] Couprie | Grayscale image processing using topological operators[END_REF]. Proposed operators [START_REF] Couprie | Grayscale image processing using topological operators[END_REF] were inserted in a powerful chain, see figure 1, which offers high segmentation quality. But the major drawback of this chain is the processing time with is very considerable. It is also impo rtant to munch that no context for specific use has been defined. We propose to imp rove this chain not only for better leg ibility but also for an execution time wh ich respects real-time image processing constraint. Medical image processing will be the new defined application context. The original chain is based on edge detection approach. It starts by extracting edges fro m the input image thanks to Deriche filter, then it calls upon the thinning operator followed by the crest restoring one and finally a thresholding process is lunched. In order to improve legib ility and also execution time, important ameliorations aiming not only algorithm's level but also architecture's level were applied: explo iting Pentiu m performance and specially the associated SEE instructions set allo w us to assure real time segmentation.

The remainder of the paper is organized as follows. Section II gives a quick review of theory bases of chooses algorith ms and exp lains the different steps of the original chain. Details about enhancement, implementation and applied optimization are g iven in section III. We show also obtained results . Finally, the conclusions recall our research work.

THEORY BASES & ORIGINAL CHAIN

In this section, we start by recalling some basic notion introduced by Bertrand and Couprie [START_REF] Couprie | Grayscale image processing using topological operators[END_REF] to modify the topology of grayscale image in order to filter out the non-significant informat ion. Then we develop different steps of the original chain.

Authors [START_REF] Couprie | Grayscale image processing using topological operators[END_REF] introduce the notion of -destructible point which allowed them to selectively alter the topology, based on a local contrast parameter ( ). Inspired by the same princip le, authors [START_REF] Couprie | Grayscale image processing using topological operators[END_REF] introduce the notion of -end and -deletable points in order to apply the same procedure on grayscale skeleton. To eliminate non significant branches and region minima, -end and -deletable points will be lowered what allows to get filtered skeleton. But the majo r disadvantage of this procedure, as any others procedures based on the dynamic, is that if two neighboring minima are separated by a high crest which contains a narrow and deep pass, the dynamic between these two minima only depends on the altitude of the pass. To cure that a crest restoring method was introduced for detecting and eliminating narrow passes, obviously on grayscale image.

As we munched in the introduction the original chain is based on edge detection approac h. The first step is edges extraction thanks to Deriche filter. Fo llo wing this detection, -destructible points will be selected under the control of a local contrast parameter for a selective mod ification of topology. We lower values of the selected up pixels until the point that they are not -destructible any more. This operation will be repeated until stability and guarantees a strict preservation of the topology. Even obtained objects at this stage are thin and lengthened; they frequently suffer fro m cuts under the effect of remaining noise. The goal is to reconnect them again. If we regard an image as a topographic relief, a threadlike object can be described like a crest line and the cracks will correspond to the narrow passes in this line. To cure this problem, authors [START_REF] Couprie | Grayscale image processing using topological operators[END_REF] propose an algorithm allowing detecting and raising these narrow passes. They introduce the notion of extensible separating point (points having under two connexit ies) and according to its local topological characteristics, separating point of the thinned image are classified in constructible point, saddle point, end point or an isolated point. Only saddle points will be raised on the same level as those of the crest. Obtained thinned skeleton with restored crest is the input of the thresholding procedure. The output image is a binary image and presents the segmentation that we are looking for. In figure 2, all steps of the original chain are illustrated. Gradient image (b) is obtained from (a) by the application of Garcia Lorca gradient operator. At first, crest has been reduced to thin lines (c) by a filtered thinning with =5, then restored (d) by executing the Crest Restoration algorithm until stability. To highlight object of interest a threshold has been applied (d) given (e). It is clear that no artifacts have been introduced by the executed chain. Compared to popular method known as hysteresis thresholding, the proposed method gives mush better results. It is also important to recall that the origin chain is a general interest chain. Its use is not framed by a specific context even if we used 2D images of human brain cut fo r testing.

OPTIMIZATION & RESULTS

To obtain the best execution time, it is necessary to have an overall view of the intervening agents and to maintain in spirit reached objective: In this wo rk, we deal with a segmentation process executed on Pentiu m microprocessor (P4 with 128 b its). We will start by giving a first evaluation of the original chain's execution t ime to be able to fix the first targeted functions for optimization. Then we exp lain the algorithmic improvements brought before finishing with the architectural ameliorations.

Within the framewo rk of algorith mic improvements, a new segmen tation's chain will be proposed. The new chain is inspired by the old one proposed by Bertrand and Couprie [START_REF] Couprie | Grayscale image processing using topological operators[END_REF] and adapted for med ical images processing. Also basing itself on the contour approach, it uses the two principal algorith ms of thinning and crests restoration while including three new algorithms to ensure a dynamic redistribution phase, a double thresholding phase and a labeling phase. The new chain is a co mplete chain which offers a high segmentation quality and it is much faster than the origina l one. Concerning architectural side, applied optimization's techniques have allowed a spectacular profit in segmentation's times. Let take the case of an image o f 768 Kbytes size (512* 512), its segmentation's time passed fro m 3,6s to 0.03s that allo ws a possible cadency of 33 image/s.

First evaluati on

In our analyses, we called upon Gprof as being the standard profiler of GNU world. After first evaluation, it proved that segmentation time of the original chain is extremely important. Indeed, only image s having sizes lower than 48kbytes can be segmented in less than one second, the reason is the quasi reduced number of the treated pixels of these small sizes images. But since the mo ment that we process images having size higher then 750kbytes, the segmen tation time exceeds one second. For example images (512* 512) take 3.6s to be segmented. According to this profiling, we note that thinning phase consumes more than 50% of seg mentation's total time. The iterat ive criterion (repetition until stability) of the two algorithms (thinning and crest restoring) can explain this percentage. We should not also forget that the pixels resultants from contours detection phase are numerous and cover the 255 gray's levels, this big number of level slows down the algorithm of thinning. The algorithm of crest restoring is slowed down by the extensive damage on the crest's lines. These damages are due to the presence of the noise in the acquired images. We will try to solve these problems and to remove obstacles which penalize the segmentation's time.

Algorithmic enhancement

The mod ifications associated to this level aim the thinning and crest restoring algorith ms in first. The goals of these modifications can be classified as follows: the first priority to minimize the number of p rocessed pixels by the thinning algorith m, then to reduce damage in crest lines by reducing noise area, to improve in term of quality the legibility of the final resulting image, and obviously to min imize throttling nodes in targeted algorithms.

Ne w chain organization

The first factor which slows down the thinning phase is the number of treated pixels. The direct passage from the contour detection towards thinning penalizes enormously the segmentation time. For this reason, we propose to advance the thresholding phase before thinning and that will allow to reduce the number of treated pixels and consequently to accelerate thinning. Thanks to this reorganizat ion, we obtain an average profit of about 72%. The seg mentation's interval of the whole tested images passed from [0:14]s to [0:2]s. The best profits were relating to big size images. Results of the new reorganized chain are pro mising and make it possible to look for real -times segmentation follo wing other improvements. At this stage, only nine images (512* 512) are seg mented per second.

New dynamic redistributi on phase

When we introduced crest restoring phase. We mentioned that cuts on thinned objects are due to the remaining noise in the image. As we explain it before if two neighboring min i ma are separated by a high crest which contains a narrow and deep pass, the dynamic between these two min ima only depends on the altitude of the pass .Thus to reduce the number of processed pixels by the algorithm of crest restoring we propose to add a new dynamic red istribution phase at the beginning of the chain. To achieve this aim, we chose the dynamic's redistribution as pretreatment's method. This technique consists in transforming image's histogram so that this one occupies all dynamics. It is very i mportant mention that the topology will be p reserved even the associated dynamic of each pixel will change. It is also impo rtant to mention that thanks to the thresholding process, enhanced noise by the dynamic's redistribution will be eliminated. The addition of this new dynamic's redistribution phase has improved the quality of objects resulting from the thinning phase. The objects are less damaged by noise. Crests lines are more continuous with fewer disruptions. Thus it will be fewer narrow passes to be detected and eliminated by crest restoring algorith m. Fig. 3. Adding new labeling phase

New Labeling phase

The threshold, in the original chain, aimed to highlight the skeletal objects. When we advanced this phase before the thinning, we changed its algorith m and we have defined a new role (minimizing the number of processed pixels).

Thus its original purpose is no longer assured. For this we propose to add a new phase of labeling to address this problem. The new algorith m has to increase the readability of the final image to ensure a better quality of segmentation and it should not penalize the segmentation time that we seeks to improve since the beginning of this work. Once the restoration of crest is completed, we can associate to the s cene the list of areas delimited by closed contours , see figure 3. These areas will be located by a label. This label is characterized a color. Operator of area labeling associates the same color for all p ixels of the same area.

In figure 4, all steps of the new chain are illustrated. After the dynamic redistribution of the original image (a) we obtain (b). Grad ient image (c) is obtained from (b) by the application of Garcia Lorca gradient operator then a double thresholding process is lunched to get (d). Crest line has firstly been reduced to thin lines (e) by a filtered thinning with =5, then restored (f) by executing the Crest Restoration algorithm until stability. To highlight object of interest a labelling algorithm has been applied to (f) given (g). The new chain provides great segmentation's quality. The various objects of interest are present in the final image. They are presented in more exp licit way than ones offered by the original chain.

Architectural level opti mization

Algorith mic optimization can be necessary in very specific cases, but it remains always independent of target architecture. However on RISC machines which have a large registers and where the effectivene ss requires the reunification of the identical instructions to profit fro m the effect p ipeline, we can imp lement a code even more powerful than the one obtained after algorithmic optimizat ion. This aim can be achieved by explo iting SIMD technology which rests on a revolutionary technique, namely to treat several data by loading only one instruction.

Compiler opti mization opti ons

We begin by activating the optimization options of the used compiler GCC. There is a very slight improvement only for large sizes images. Pro fit in segmentation's time is very limited, even with the activation of -O3 option.

Memory and register management

Some SSE2 instructions require that its arguments be aligned on 16-byte boundaries when moving to (/from) memory (otherwise the instruction throws exception). For that we have defines the follo wing procedure allowing the redirection of _aligned_malloc (ms RCT) and the memo rizing of allocated size. Others those variables of transitions, the memo ry align ment is very important at res ervation's time o f large tables: matrices of the pixels, table of the topological nu mbers.

Even if we imp rove memory access time, it is certain that if operands are placed in the registers the execution time will be faster. But the number of available reg isters is limited, so it is almost impossible to eliminate all memo ry access for intermediates results. Let take the example of the associated values for each point in thinning algorithm. For each pixel, we need to compute its index (according to the eight neighbors) and we need also to compute associated alpha. To avoid double memory access, we coded both parameters in the same integer (by shift): index will be code on 24 bits and alpha on the 8 remaining bits. Thus for each point we will use only 32 b its what enables to charge values of four points at the same register (see figure 5).

Intra-Register vectorization

Pro mote the pipeline is one of the key points to improve the performance of programs. To explo it the potential parallelism offered by the Pentiu m 4, reg ister vectorizat ion can be applied by unrolling serial loops into a form that takes advantage of instructions provided by SSE2.

Anyone who looks carefully to thinning algorithm realizes that the two phases of detection and removing of destructible point are independent. During the detection, the eight neighbors of each pixel are charged according to a precise chronological order. Thus, neighbors of any pixel (x) and its neighbor (x+1) can be exanimate at t he same time without ambiguity. We also point out that each pixel indicated as a destructible point is coded with its alpha in the same integer (32 bits) then piled up. To benefit fro m the 128 availab le bits, we can launch up to four instructions of stacking in parallel. We draw attention to the fact that the vectorization is effect ive only in the case of large sizes loops. If the number of iterat ion is small, no contribution will be v isible. It is also very important to check variab les sharing within th e loop. Even in the case of a conditional dependency, the probability to increase latency, instead of decrease it, is great.

External loop unrolling

The whole of applied operations to characterize point's nature uses piles. Let us take the case of restoring crest algorith m, which uses in a first stage a pile to store the pre-divergent points with their alpha. The order in which points are stored has no importance thus we can launch this loop several times in parallel. For efficient loops unrolling, extra branches must be eliminated fro m the loops, unrolling also enhances other optimizations, such as prefetch, scalar replacement, and instruction scheduling. Finding the unroll factor is important to determine the amount of loop expansion. In the present case, the management of the pile is the only obstacle. If we launch more than two loops at the same time, it will be d ifficult to manage the pile's sharing resulting in auto matic increase in latency. Thus, we are limited to a double external loop unrolling for the better manages of the pile.

CONCLUSION

In this work, we proposed a new image processing chain, which offers high quality segmentation respecting real time expiries. The new chain is perfectly suited to med ical applicat ion. Acting directly over grayscale images while preserving their topologies, all topological operators and in particular thinning and crest restoration operators have given birth to a powerfu l image segmentation of 2D b rain's image cuts. The association of algorithmic optimization techniques to architectural optimization techniques relative to Pentium processor (P4 3.06Gh z), allows to expect a cadency of 33 images (512 * 512) per second, see figure 7. These real-time performance confirms the interest of an "Algorithm Architecture Adequation" approach, see figure 6. For better adaptation between the chain and the target architecture, we begin by carving algorithms to ensure better stability and readability of code. Then, we make calls to techniques that exploit the best performance of the processor.

Fig. 1 .

 1 Fig. 1. Original segmentation chain

Fig. 2 .

 2 Fig. 2. Obtained results by the original chain (d) (e) (b) (c) (a)

Fig. 4 .

 4 Fig. 4. Obtained results given by the new chain

  void *benchmalloc (size_t size) 01: void *ptr = (void *) _aligned_malloc(4*size); 02: if(ptr) alloc_mem += size; 03: return ptr;

Fig. 5 .

 5 Fig. 5. Register management

Fig. 6 .

 6 Fig. 6. Summary of performance -(a) the original chain's performance, (b) after reorganizing, (c) after dynamic redistribution, (d) after algorithmic optimization, (e) after compiler optimization and (f) after architectural optimization.

Fig. 7 .

 7 Fig. 7. Summary of performance in terms of caden cy


Ramzi Mahmoudi
email: mahmoudr@esiee.fr

Mohamed Akil
email: akilm@esiee.fr

Petr Matas
email: matasp@esiee.fr

Parallel Image Thinning Through Topological Operators On Shared Memory Parallel Machines

In this paper, we present a concurrent implementation of a powerful topological thinning operator. This operator is able to act directly over grayscale images without modifying their topology. We introduce an adapted parallelization methodology which combines split, distribute and merge (SDM) strategy and mixed parallelism techniques (data and thread parallelism). The introduced strategy allows efficient parallelization of a large class of topological operators including, mainly, ߣ-leveling, skeletonization and crest restoring algorithms. To achieve a good speedup, we cared about coordination of threads. Distributed work during thinning process is done by a variable number of threads. Tests on 2D grayscale image (512*512), using shared memory parallel machine (SMPM) with 8 CPU cores (2× Xeon E5405 running at frequency of 2 GHz), showed an enhancement of 6.2 with a maximum achieved cadency of 125 images/s using 8 threads.

Introduction

In 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] Anil | Fundamentals of Digital Image Processing[END_REF]. Topological thinning and skeletonization are ones of the most cardinal operators for this kind of preprocessing, especially since the development, by our team, of an efficient thinning algorithm able to act directly over grayscale image [START_REF] Couprie | Topological operators for grayscale image processing[END_REF]. Using topological operators allows topology preservation which results in conservation of important significant information [START_REF] Everat | New topological operators for segmentation[END_REF]. This conservation was impossible in the case of binary image processing [START_REF] Ojala | Generalized local binary pattern operators for multi-resolution grayscale and rotation invariant texture classification[END_REF]. Early thinning algorithms were designed for serial implementations, but since parallel computers are available several approaches have been developed with parallel processing [START_REF] Heydorn | Optimization and Performance Analysis of Thinning Algorithms on Parallel Computers[END_REF][START_REF] Bourbakis | Design of an array processor for parallel skeletonization of images[END_REF]. In [START_REF] Heydorn | Optimization and Performance Analysis of Thinning Algorithms on Parallel Computers[END_REF], Heydorn presents a concept for an implementation of different parallel thinning algorithms on parallel processors. The emphasis is put on a good parallelization using fine granularity and the simultaneous usage of vectorization.

This paper describes an adapted parallelization methodology combining split, distribute and merge (SDM) strategy based upon the well-known principle of divide-and-conquer and thread coordination which allows an efficient parallelism for introduced thinning operator on shared memory machines. Proposed strategy can also be applied for any topological operator having the same characteristics based on elementary operations of point characterization and similar algorithmic structure as we will demonstrate later.

This paper is organized as follows: in section 2, some basic notions of topological operators are summarized; the original algorithm of thinning is introduced. We define also the class of operators that our parallelization strategy may cover. In section 3, parallelization strategy, that has been adopted, is introduced. In section 4, Thread coordination and synchronization is discussed. In section 5, as a concrete example of introduced strategy and threads synchronization techniques, a parallel version of thinning algorithm is presented. Experimental analyzes results of different implementations are also presented and discussed. Finally, we conclude with summary and future work in section 6.

Topological thinning operator

Skeletonization and thinning are major applications of topology in image processing. A great number of thinning algorithms for binary images have been developped [START_REF] Lam | Thinning methodologies: a comprehensive survey[END_REF]. The use of this kind of images assumes a prior segmentation which implies a loss of information. Some attention has been given to the development of thinning algorithms acting directly over grayscale images. Dyer and Rosenfeld [START_REF] Dyer | Thinning Algorithms for Gray-Scale Pictures[END_REF] proposed an algorithm based on a weighted notion of connectedness. The thinning is done directly over the graylevels of the points but, as the authors showed, the connectivity of objects is not always preserved. Other works [START_REF] Shiaw-Shian | A new thinning algorithm for grayscale images by the relaxation technique[END_REF] use an implicit image binarization into a background and a graylevel foreground. The graylevel information guides the removal of points of the foreground that are simple, in the binary sense. This technique makes it possible to obtain certain desired geometric properties. Inspired by this technique, M. Couprie and G. Bertrand [START_REF] Couprie | Topological operators for grayscale image processing[END_REF] propose a filtered thinning method that allows to selectively simplify the topology, based on a local contrast parameter λ .To achieve this simplification; they introduce the notion of λ-destructible point which is more flexible then the notion of destructible point. This algorithm is the one we are going to present en details in the following to illustrate our parallelization strategy.

Theoretical 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 n-adjacency for the section ܨ of Ϝ , we must use ݊-adjacency for the section ܨ with (݊, ݊) = [START_REF] Dyer | Thinning Algorithms for Gray-Scale Pictures[END_REF][START_REF] Ojala | Generalized local binary pattern operators for multi-resolution grayscale and rotation invariant texture classification[END_REF] or [START_REF] Ojala | Generalized local binary pattern operators for multi-resolution grayscale and rotation invariant texture classification[END_REF][START_REF] Dyer | Thinning Algorithms for Gray-Scale Pictures[END_REF]. 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 with T(x, X) and T(x, X) the two connectivity numbers defined as follows (# ܺstands for the cardinal of ܺ):

,ݔ(ܶ ܺ) = ܥ# ,ݔ[ Γ ଼ * )ݔ( ∩ ܺ]; ,ݔ(ܶ ܺ) = ܥ# ,ݔ[ Γ ଼ * )ݔ( ∩ ܺ];
So we can define the four neighborhoods:

߁ ାା ,ݔ( )ܨ = ሼݕ ∈ ߁ ଼ * ;)ݔ( )ݕ(ܨ > )ሽݔ(ܨ ߁ ା ,ݔ( )ܨ = ሼݕ ∈ ߁ ଼ * ;)ݔ( )ݕ(ܨ ≥ )ሽݔ(ܨ ߁ ିି ,ݔ( )ܨ = ሼݕ ∈ ߁ ଼ * ;)ݔ( )ݕ(ܨ < )ሽݔ(ܨ ߙ ି ,ݔ( )ܨ = ൜ ,)ݕ(ܨሼݔܽ݉ ݕ ∈ ߁ ିି ,ݔ( ,)ܨ ݂݅߁ ିି ,ݔ( )ܨ ≠ ∅ሽ )ݔ(ܨ ݁ݏ݅ݓݎ݁ℎݐ
We define also some associated connectivity numbers:

ܶ ା ,ݔ( )ܨ = ܥ# ,ݔ[ ߁ ା ,ݔ( ])ܨ ܶ ାା ,ݔ( )ܨ = ܥ# ,ݔ[ ߁ ାା ,ݔ( ])ܨ ܶ ିି ,ݔ( )ܨ = ܥ# ,ݔ[ ߁ ିି ,ݔ( ])ܨ
Furthermore, the connectivity numbers allow the classification of the topological characteristics of a point:

ݔ is a peak point if T ା ,ݔ( )ܨ = 0. ݔ is a k-divergent if T ିି ,ݔ( )ܨ = ݇ ℎݐ݅ݓ ݇ > 1.
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 is destructible or x is k-divergent and 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) > ߣ.

Original algorithm

∀ F ∈ φ, we say that G ∈ φ is a skeleton of F if G is
obtained from F by iteratively selecting a destructible and non-end point in F and lowering it down to α ି (x, F), 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, 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 λ = 10.

Algorithm : λ -Skeleton (input : Ϝ ∈ ߮, ߣ ∈ Ν; output : Ϝ) 1. Repeat until stability 2. Among all the points which are λ-deletable and not λ-end 3.

Select a point x of minimal value ; 4.

Ϝ()ݔ ≔ ߙ ି ,ݔ( ;)ܨ (

Class of operators based upon point characterization in the grayscale image case

Bertrand [START_REF] Anil | Fundamentals of Digital Image Processing[END_REF][START_REF] Bertrand | Topological approach to image segmentation[END_REF] introduced connectivity numbers for grayscale image as showed in section 2.1. 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 based on w-thinning 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 recursivety 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. 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 strategies

Multiprocessor chips make computing more efficient by exploiting parallelism which is one of the outstanding challenges of modern computer sciences. Exploiting such parallelism depends on the way of scheduling tasks to different processors such that the tasks can be computed simultaneously in parallel [START_REF] Mattson | Patterns for Parallel Programming[END_REF]. Computing each individual task in parallel using all the processors and computing tasks one after the other is Data Parallelism [START_REF] Hatcher | Data-Parallel Programming on MIMD Computers[END_REF]. For both strategies, programming challenges arise at all scales of multiprocessor systems: at the small scale, processors within a single chip need to coordinate access to shared memory locations; at the large scale, processors in a super computer need to coordinate routing of data. It is also possible to combine the mentioned strategies for better scheduling [START_REF] Feldmann | Dynamic scheduling on parallel machines[END_REF]; such strategy is called Mixed Parallelism. In this case, challenges are also related to the asynchronous criteria of modern computers: activities can be halted without warning by interrupts, preemption or frequently by cache misses.

In a more global frame, better strategies taking advantage of such parallelism to improve computational speed are based on the well known principle of divide and conquer. The application of this principle cannot be independent from the type of algorithm [START_REF] Wangqing | A New Divide and Conquer Algorithm for Graph-based Image and Video Segmentation[END_REF]. Indeed, application of this principle to divide the initial problem and then application of Mixed Parallelism strategy during parallel subproblems processing seems to be sufficient. But it is only true for static parallel algorithms in which each thread can achieve its work "independently" from the other. Low-level image processing algorithms are a good example of this class because they have a high degree of locality allowing different segments of the image to be treated independently by different processors [START_REF] Seinstra | A software architecture for user transparent parallel image processings[END_REF]. Other global operators like Fourier transform and Euclidean distance transform are separable, allowing parallelization by defining a direction for computing pixels [START_REF] Meijster | A general algorithm for computing distance transforms in linear time[END_REF].

For target algorithms, as we shown in section 2.3, get some iterative criteria and evolutes until stability. Intermediate results need also to be stored. Each time, that a pixel is lowered, a new process for inserting its neighborhoods is launched. So threads need imperatively to communicate and to share the same queue; this is why we return to dynamically parallel algorithms in which threads can interact with one another. Through parallelization strategies presented in the beginning of this section, we see that for an interprocessor parallelism based on divide and conquer principle; better performance can be achieved by the use of mixed parallelism, since it allows us to combine SDM-Strategy and coordination of threads. And as our processing continues until stability, we primarily focus on an approach where data parallelism is used at upper levels. At lower levels of the processing, we will switch to threads parallelism and coordination to compute parallel read/write for managing cacheresident data. If we observe carefully the studied class, we see that there are two fundamental stages: the first one is to characterize a point. Then, according to the nature of this point, we decide to eliminate it (modify its value) or not. If one pixel is lowered, it becomes necessary to re-examine its eight neighbors. So we can follow these steps to apply divide and conquer principle for our class of algorithms.

How to Split

In upper level, search space is subdivided into smaller regions, and bounds are found on all solutions contained in each sub-region under consideration. Usually, dividing original image is not advised, when dealing with topological operators, warning topology is not preserved. But pixel characterization procedure can be split into sub-procedures. So we can characterize in parallel more than one point during a single iteration.

How to distribute

All algorithms associated to our class are executed in a loop until stability for example: no more λdeletable and not λ-end for thinning. Thus we can specify translation states. The initial system will undergo an evolution until reaching a stable state. During the evolution, sub-procedures defined in section 3.1 are distributed among used threads. Usually thread evolution is uniform, but, due to data dependency, thread evolution must be dynamic. The number of threads is changing during the whole processing procedure. So the second stage of the algorithm can be realized. If a point is characterized, its value is lowered and its eight neighbors will be inserted in a FIFO queue. Since one thread terminates, it will generate a new thread to repeat the same work with new inserted neighbors. This is how we can plan distribution.

How to merge

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 threads. But as we mentioned in section 3.2, we are dealing with a dynamic evolution so we can plan the following: since two threads finished, they directly merge and a new thread is created. This implies the creation of some shared FIFO queue containing all inserted neighbors by both two parent threads. Only one shared data structure will contain pixels lowered by all threads. In threads merging, there is no hierarchical order, the only criteria is finish time. It is also important to mention that only newly created threads can modify the created FIFO queue and one neighbor cannot be inserted twice. It is a precaution in order to minimize consumed cache. In Fig. 2, we illustrate the introduced SDMstrategy. The original shared data structure, containing all pixels, is divided in n research zones ሼz ଵ , z ଶ , . . z ୬ିଵ , z ୬ ሽ. We associate one thread from the following list ሼT ଵ , T ଶ , . . T ୬ିଵ , T ୬ ሽ to each zone. Each thread can browse freely its zone and if it detects target types, it lowers the characterized pixel and it pushes its eight neighbors in one of the shared concurrent FIFO queues. One queue cannot be shared by more than two threads. There is no hierarchical order in merging of threads. Queues are attributed for the first two threads which have arrived (first-come, first-served). Since two threads finish their work, a new thread is created to browse their FIFO queue. For storing new value of possibly found target point, the new thread has full access to both original zones. Example, thread T ୟ merging from T ଵ and T ଶ got access to Z ୟ with Z ୟ = Z ଵ + Z ଶ .

Coordination of threads

Here is the second major challenge in multi-core multithread architecture programming. In an ideal case, moving from one-core to multi-core should provide nfold 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 interprocessor coordination. Let's come back to our algorithm, consider eight threads which cross eight search spaces in order to characterize pixels then push its eight neighbors in a FIFO queue. This kind of analysis was evoked in many researches. Let's focus on Amdahl's Law [START_REF] Amdahl | Validity of the single-processor approach to achieving large scale computing capabilities[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. Definition of the speedup S of a work is the ratio between the time it takes one processor to complete the work versus the time it takes ݊ concurrent processors to complete the same work. Amdahl's Law defines the maximum speedup ܵ that can be achieved by ݊ processors collaborating on an application, where  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 ( ݊)

⁄ seconds and the sequential part takes (1 -) seconds. Overall, the parallelized computation takes (1 - + ) seconds.

So the speedup is:

(݊) = ଵ ଵିା .
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 in section 3.3, 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 lock alternatives for sharedmoney multiprocessors[END_REF], for illustration we propose figure 4, 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.

Performance testing

Parallel ૃ -Skeleton algorithm

Now we present a parallel version of the thinning according to the concepts 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 sub-space 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 ୱ୬ will be the newly defined set to explore since the 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 ୩୧ାଵ ).

Algorithm :Dynamically Parallel λ -Skeleton 1. ݎܨ ݈݈ܽ  ∈ ݁݉ܫ ݀ 2.
݂݅ ൫݉ < ݁݉ܫ )( < ݉ ௦௨ ൯ݐℎ݁݊ ܧ ← ܧ ∪ ሼሽ;

Experimental analyses

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. 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 measure Ψ(݊) is given by the following formula With n the number of processors:

Ψ(݊) = ݈ܽ݅ݐ݊݁ݑݍ݁ݏ ݁݉݅ݐ (݊ * ⁄ ݈݈݈݁ܽݎܽ )݁݉݅ݐ
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 1). 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 Figure 5 On the eight-core machine, wall-clock execution time for the first implementation using a lock-based 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 spinwait 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 figure 6 and figure 8 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 the speedup has increased to 6.2 ± 0.01. Another common result between figure 5 and figure 7 is stability of execution time on each n-core machine since the code uses n or more threads. 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 Figure 9, 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 ܵ(݊) = ்(ଵ)

்()

. Then the defined efficiency Ψ(݊) = ܶ ௦ (݊ * ⁄ ܶ ) can be written Fig. 8: Performance improvement for the second implementation using a spin-wait FIFO queue.

as 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. Fig. 9: Efficiency improvement

Ψ(݊) = ܶ ௦ (݊ * ⁄ ܶ ) = S (

Conclusion

In this paper, we have presented a new parallel version of the λ-Skeleton algorithm. We have also presented a adapted parallelization strategy combining Split Distribute and Merge (SDM) strategy and mixed parallelism techniques. SDM-strategy was a conditional application of the well known principle of divide and conquer. Associated mixed parallelization techniques were data parallelism at upper levels and thread parallelism at lower levels of the processing. The first major contribution in this paper is the non-specific nature of the proposed parallelization strategy. In fact, the introduced strategy can be applied to a large class of topological operators introduced in section 2.3.

The second contribution concerned threads parallelism and more specifically threads coordination and communication during computing dependently parallel read/write for managing cache-resident data which present a substantial problem. The problem addressed by this paper is how to deal with shared FIFO queue which requires inter-process coordination and communication in an essential way. And thanks to combination of spinning and waiting techniques, the proposed algorithm shows a good degree of speed-up using eight threads (about 6.2 on eight cores of the 2× Xeon E5405, about 3.1 on the Quad-cores of Xeon E5335 and 1.8 on Core 2Duo E8400).

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

  a) (b) Fig. 1. (a): after Deriche gradient operator; (b) filtered skeleton with λ = 10 .

Fig. 2 :

 2 Fig. 2: Merging of threads and associated area of activities.

Fig. 4 :

 4 Fig. 4: Spin-wait Synchronization

Fig. 5 :

 5 Fig. 5: wall-clock execution time for the first implementation using a lock-based shared FIFO queue.

Fig. 6 :

 6 Fig.6: Performance improvement for the first implementation using a lock-based shared FIFO queue.

Fig. 7 :

 7 Fig. 7: wall-clock execution time for the second implementation using a spin-wait shared FIFO queue.

  be written as follows: Ψ(݊) = ଵ * (ଵି)ା

  based shared FIFO queue Using private-shared FIFO queue Using lock-based synchronization Using spin-wait synchronization

Table 1 .

 1 and Figure 7. Characteristics of processors

		P4 660	E8400	E5335	E5405
	CPU	3.6	3	2	2
	Speed	GHz	GHz	GHz	GHz
	Bus	800	1333	1333	1333
	Speed	MHz	MHz	MHz	MHz
	L2 Size	4 MB	6 MB	8 MB	12 MB
	L2 Speed 3.6	3	2	2
		GHz	GHz	GHz	GHz


Ramzi Mahmoudi
email: mahmoudr@esiee.fr

Mohamed Akil
email: akilm@esiee.fr

Real Time Topological Image Smoothing On Shared Memory Parallel Machines

Smoothing filter is the method of choice for image preprocessing and pattern recognition. We present a new concurrent method for smoothing 2D object in binary case. Proposed method provides a parallel computation while preserving the topology by using homotopic transformations. We introduce an adapted parallelization strategy called split, distribute and merge (SDM) strategy which allows efficient parallelization of a large class of topological operators including, mainly, smoothing, skeletonization, and watershed algorithms. To achieve a good speedup, we cared about task scheduling. Distributed work during smoothing process is done by a variable number of threads. Tests on 2D binary image (512*512), using shared memory parallel machine (SMPM) with 8 CPU cores (2× Xeon E5405 running at frequency of 2 GHz), showed an enhancement of 5.2 thus a cadency of 32 images per second is achieved.

Introduction

Smoothing is a fundamental step to reduce noise and prepare the image for subsequent processing such as segmentation. 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. 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 the curvature plot of the object's contour [START_REF] Leymarie | Curvature morphology[END_REF]. Unfortunately none of these operators preserve the topology (number of connected components) of the original image. In 2004, our team introduces a new method for smoothing 2D and 3D objects in binary images while preserving topology [START_REF] Couprie | Topology preserving alternating sequential filter for smoothing 2D and 3D objects[END_REF].

Objects are defined as sets of grid points and topology preservation is ensured by the exclusive use of homotopic transformations defined in the framework of digital topology [START_REF] Kong | Digital topology: introduction and survey[END_REF]. Smoothness is obtained by the use of morphological openings and closings by metric discs or balls of increasing radius, in the manner of alternating sequential filters [START_REF] Sternberg | Grayscale Morphology[END_REF]. The authors' efforts have brought about two major issues such as preservation of the topology and the multitude of objects in the scene to smooth out without worrying about the latency of their filter. This paper describes a new parallel method for topological smoothing which provides real time image processing. We present also an adapted parallelization strategy, called Split Distribute and Merge (SD&M). Our strategy is designed specifically for topological operators parallelizing on shared memory architectures. The new strategy is based upon the excusive combination of two patterns: divide and conquer and event-based coordination. This paper is organized as follows: in section 2, some basic notions of topological operators are summarized; the original smoothing filter is introduced. In section 3, parallelization strategy, that has been adopted, is introduced. We define the class of operators that our parallelization strategy may cover. In section 4, concurrent computation method of topological smoothing is presented. Experimental analyzes results of different implementations are also presented and discussed. Finally, we conclude with summary and future work in section 5.

Theoretical background

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

1/2 2 2 1 1 2 2 ( , ) ( ) ( ) d x y x y x y   = - - -   .
An operator in E is a mapping from ( ) P E into ( ) P E , where ( ) { }

P E denotes the set of all subsets ( ) E ⊂ . Let r be an integer, the dilation ( ) r δ by r B is defined by ( ) ( ) r x X r X B x δ ∈ = ∪ , ( )
X P E ∀ ∈
4 1 1 2 2 ( ) ; 1 x y E y x y x Γ = ∈ -+ -≤ , { } 8 1 1 2 2 ( ) ; max , 1 x y E y x y x Γ = ∈ - -≤ .
In the following we will denote by n a number such that 4 n = or 8 n = thus we define:

{ } * ( ) ( ) \ n n x x x Γ = Γ . We say that a point y E ∈ is n adjacent - to x E ∈ if * ( ) n y x ∈ Γ
. Thus we can generalize for two points x , y belonging to X . They are n-connected in X if there is a n-path in X between these two points.

The equivalence classes for this relation are the nconnected components of X . A subset X of E is said to be n-connected if it consists of exactly one nconnected component. The set composed of all nconnected components of X which are n-adjacent to a point x is denoted by [ ] , n C x X . In order to have a correspondence between the topology of X and the topology of X , we use 'n-adjacency' for X and ' 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 the topology of X . In other words, we can remove the point p from X without changing the topology of X .

The point x X

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

( )

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

x E ∈ is simple for ( , ) 1 X E T x X ⊆ ↔ = and ( , ) 1 T x X =
The homotopic alternating sequential filter is a composition of homotopic cuttings and fillings by balls of increasing radius. It takes an original image X and a control image C as input, and smoothes X while respecting the topology of X and the geometrical constraints implicitly represented by C , simple illustration is given by figure 1. Based on this filter, authors [START_REF] Couprie | Topology preserving alternating sequential filter for smoothing 2D and 3D objects[END_REF] introduce a general smoothing procedure with a single parameter which allows controlling the degree of smoothing.

Let X be any finite subset of E , let C X ⊆ , r ∈ ℕ and D X ⊆ . The homotopic alternating sequential filter ( ) HASF of order n with constraint sets C ; D is defined as follows: ( ) ( )

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

( )

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

Parallelization Strategy

In this section, we start by defining the class of topological algorithms. Then we present our motivation to parallelize it on parallel shared memory machines. Finally we will introduce different steps of our parallelization strategy. We will focus especially on distribution phase and tasks scheduling over different processors.

Class of topological algorithms

In 1996, Bertrand and Couprie [START_REF] Bertrand | Topological approach to image segmentation[END_REF] introduced connectivity numbers for grayscale image. These numbers describe locally (in a neighborhood of 3x3) the topology of a point. According to this description any point can be characterized following its topological characteristics. They also introduced some elementary operations able to modify gray level of a point without modifying image topology. These elementary operations of point characterization present the fundamental link of large class of topological operators including, mainly, skeletonization and crest restoring algorithms [START_REF] Couprie | Topological operators for greyscale image processing[END_REF]. This class can also be extended, under condition, to homotopic kernel and leveling kernel transformation [START_REF] Bertrand | Image segmentation through operators based on topology[END_REF], topological watershed algorithm [START_REF] Bertrand | On Topological Watersheds[END_REF] and topological 2D object smoothing algorithm [START_REF] Couprie | Topology preserving alternating sequential filter for smoothing 2D and 3D objects[END_REF] which is the subject of this article. All mentioned algorithms get also many algorithmic structure similarities. In fact associated characterizations procedures evolve until stability 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. 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] Meijster | A general algorithm for computing distance transforms in linear time[END_REF] for local operators and point to point operators can also be cited as example. Hence the need of a common parallelization strategy for topological operators that offers an adapted algorithm structure design space. Chosen algorithm structure patterns that will be used in the design must be suitable for SMP machines.

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

Split Distribute and Merge Strategy

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

Split phase

The Divide and Conquer pattern is applied first by recursively breaking down 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 computation.

Distribute phase

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

In effect, non-real-time system scheduler doesn't know, in advance, necessary time to perform each task. Thus, "Symmetric Multiprocessing" scheduler distributes tasks to minimize total execution time without load balancing between processors. We propose a novel tasks scheduling approach to prevent improper load distribution while improving total execution time. 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] and CFS "Completely Fair Scheduler" [START_REF] Molnar | Modular Scheduler Core and Completely Fair Scheduler[END_REF]. These schedulers are based on tasks uniformity principle. Then, the worst case to process T will be:

{ } 1 2
( ) max ( ), max ( ),..., max ( )

m K T T T T = .
To demonstrate that let suppose that exist a set ( )

L T as ( ) ( ) L T K T ≥ ∑ ∑ . As 'NPS' manage ( ) L T and
( ) K T , so we can introduce the following: ( )

L T m ≤ and ( ) K T m ≤ . Now, If ( ( ) ( )) L T K T ≥ ∑ ∑ then it exist at least one task { } l , with ( ) k K T ∈ , such as : ( ( )) ( ( )) ( ) l L T l K T l k ∈ ∧ ∉ ∧ > .
This is impossible according to the definition of ( ) K T . We remember that ( ) K T was defined as the worst case.

3.2.3

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 threads. In the case of topological operators, we are dealing with a dynamic evolution process so we can plan the following: since two threads finished, they directly merge and a new thread is created. In thread's merging, there is no hierarchical order, the only criteria is finish time.

Parallel smoothing filter

In this section we start by analyzing overall structure of the original algorithm. Then we move to the parallelization of the Euclidean distance algorithm, thinning algorithm and thickening algorithm. We conclude by a performance analysis of the entire smoothing topological operator.

As we have shown in Section 2, the algorithm receives as input a binary image and the maximum radius. It uses two procedures for homotopic opening and closing. 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. The closure procedure ensures the same performance of the four consecutive functions with a single difference: the erosion instead of dilatation.

Euclidean distance computing

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. The new algorithm must have an Euclidean distance computation's 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 and al. algorithm has 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, called PSN "Propagation Using a Single Neighborhood" that uses the following structure element:

( ) ( )

2 2 4 ( ) 1 x x y x d p q q p q p   = ∨ - + - <     .
He also proposes a second algorithm, called PMN "Propagation Using Multiple Neighborhood", which uses the eight neighbors. In [START_REF] Cuisenaire | Fast and exact signed Euclidean distance transformation with linear complexity[END_REF], Cuissenaire proposes a third algorithm with 3/2 ( ) n ο complexity, which offers an accurate computation of the Euclidean distance. The only drawback of this algorithm is computation time witch 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. The low data dependence and the ability to operate on 3D images, makes PSN algorithm a potential candidate for replacing 4SED. Meijster [START_REF] Meijster | A general algorithm for computing distance transforms in linear time[END_REF] proposes an algorithm to compute exact Euclidean distance. The 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 for lines.

It is important to note that strong independence between different processing steps and computing error equal to zero makes Meijster algorithm a potential candidate to replace 4SED. It is also able to operate on 3D images. Theory analysis of Meijster and Cuissenaire algorithms can be found in Fabbri's work [START_REF] Fabbri | 2D Euclidean distance transform algorithms: A comparative survey[END_REF].

In the following, we propose a first analysis based on different algorithms implementation in order to compare between them. We have implemented 4SED algorithm using a fixed size stack. It uses a FIFO queue and it has a 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 serves only as reference for comparison. For PSN implementation, we used stacks with dynamic sizes. The memory is allocated using small blocks defined at stack creation.

When an object is added to the queue, the algorithm will use the available memory of the last block. If no space is available, a new block is allocated automatically. Block size is proportional to the image size (N x M / 100). Finally we used a simple memory structure for the implementation of Meijster algorithm.

A simple matrix was used to compute distance between points and object of each column. Three vectors were used to compute distance in each line.

Fig.2. Overall structure of the original algorithm

Results show that PSN algorithm is the most expensive in all cases (for any radius). Meijster algorithm is moderately 5x faster than 4SED. The output images returned by Meijster algorithm hold the best visual quality while Euclidean distance computation error is almost zero. So our efforts will be brought on Meijster algorithm parallelization.

We denote by I the input image with m columns and n rows. We denote by B an object included in I . The idea is to compute, for each point p I p B ∈ ∧ ∉ , the separating distance between p and the closest point b with b B ∈ . This amounts to compute the following matrix:

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

x y dt p p EDT p   =   With ( ) EDT p =
( , ) ( , 1) 1 For each row u , we note that there is, for the same point p , the same value of ( , ) f p y even if y change its values, so we can introduce the concept of "region of column ".

T x T x G p y G y p = -+ ( , ) ( , 1) 1 B x B x G p y G y p = + + .
Let S be the set of all y points such that ( , ) f p y is minimal and unique. The formula of S , 0 y u ∀ ≤ ≤ , will be ( ) min : ( , ) ( , )

p S u y f p y f p i = ≤ .
Let T be the set of all points having coordinates greater than, or equal to, horizontal coordinate of the intersection with a region. The formula of

T , ( (0 )
( )) i u u m ∀ ≤ ≤ ∧ ≤
, will be the following :

( ) ( ( 1), ) 1 x p p p T u Sep S u u = - + Let ( , )
Sep i u be the separation between regions of i and u , defined by: ( , ) ( , )

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

x x Dif G p u G p i = - .
Thus lines processing will be, at the beginning, from left to right then from right to left. During the first term, from left to right, two vectors S and T will be created. These two vectors will contain respectively all regions and all intersections. During the second treatment, from right to left, we compute f for each value of S . f is computed also for each respective values of T .

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

{ }

( , ) min ( , ), ( , )

x T x B x G p y G p y G p y = ; { } ( , ) 0 ( , ) ( , )
x

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

x p + and 1

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

( ) ( , ( ))

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

( ) ( ( ), ) 1

x p p p T u Sep S q u = + .
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 the algorithm Meisjter for both steps.

Associated algorithm complexity is ( ) [ ] 0 0 = s

( ) / n m N ο × . ( ) n

4.

[ ] 0 0 = t

5.

/* First part */ 6. for ( )

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

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

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

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

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

Thinning and thickening computing

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

Experimental analyses

Proposed parallel version of topological smoothing algorithm was implemented in C using OpenMP directives. We implemented two versions, the first one using 'Symmetric Multiprocessing' scheduler and the second one using 'non-preemptive' scheduler. Wallclock execution times 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: On the eight-core machine, Number of instructions in the second implementation drops from an average of 1879 x10 6 FI with a single thread down to 202x10 6 ms with 8 threads. As expected, the speed-up for the second implementation using 'non preemptive' scheduler is higher than for the one using "Symmetric Multiprocessing" scheduler, thanks to balanced distribution of tasks. A remarkable result, shown in figure 4, is the fact that speed-up increases as we increase the number of threads beyond the number of processors in our machine (eight cores). In a first implementation, using "Symmetric Multiprocessing" scheduler, the speedup at 8 threads is 1.9 ± 0.01. However, for the second implementation, using our scheduler, the speedup has increased to 5.2 ± 0.01. Thus execution time is decreased from 0.16± 0.005s to 0.03± 0.001s. This performance has allowed a cadency of 32 images per second (see Figure 5). This real-time performance confirms the interest of proposed parallelization strategy. Another common result between different architecture is stability of execution time on each n-core machine since the code uses n or more threads (see figure 4). For better readability of our results, we tested the efficiency of our algorithm on various architectures (see figure 5) using the Ψ ( n ) formula introduced earlier with s t sequential time on mono-processor Pentium 4 660. For parallel time we use best parallel time obtained using 8 threads with Non Preemptive scheduler.

Ψ ( n ) = s t / ( n * p t ),

FIG. 6. Efficiency improvement

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

In this paper, we have presented a new parallel computation method for topological smoothing. We have also presented an adapted parallelization strategy. SDM-strategy was a conditional application of the well known principle of divide and conquers associated to event-based coordination techniques. First major contribution in this paper is the parallel computation method for image smoothing allowing real-time processing with topology preservation. Second contribution is the non-specific nature of proposed parallelization strategy. In fact it can be applied for a large class of topological operators as we shown in section 3.1. Third contribution concern tasks distributions. We presented a non-preemptive scheduler 'NPS', simpler to implement and more adapted to particular topological algorithm processing.

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

. 2 E 4 Γ and 8 Γ

 248 The ball r B is termed as the structuring element of the dilation. The erosion ( ) introduce the notion of simple point witch is fundamental for the definition of topology preserving transformations in discrete spaces. We start by giving a definition of local characterization of simple points in = ℤ . Let consider two neighbourhoods relations defined by, for each point x E ∈ :

Fig. 1 .

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

  Figure 2 describes obtained results by the different three implementations. During this evaluation we used a binary test image (200x200). We have also modified 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).

.

  If we assume that minimum distance of an empty group K is ∞ and z K ∀ ∈ , we have ( ) y z + ∞ = ∞ then ( ) EDT p formula can be written as follow: split the Euclidian distance transform procedure into two stages. The first step is to scan columns and compute EDT for each column y . The 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:

Let

stFIG. 3 .

 3 FIG. 3. Number of Instructions and tasks distribution using 'NPS'

FIG. 4 .FIG. 5 .

 45 FIG. 4. Performance improvement using 'NPS'.

Algorithm 2: Parallel Version Meijster [2 nd Step] 1. For

  m × refers to image size and N refers to the number of processors.

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

25. endif 26. endfor 27. end forall