Vis Comput DOI 10.1007/s00371-017-1407-4

ORIGINAL ARTICLE

A modified ZS thinning algorithm by a hybrid approach Lynda Ben Boudaoud1

· Basel Solaiman2 · Abdelkamel Tari1

© Springer-Verlag Berlin Heidelberg 2017

Abstract Thinning is one of the most important techniques in the field of image processing. It is applied to erode the image of an object layer-by-layer until a skeleton is left. Several thinning algorithms allowing to get a skeleton of a binary image are already proposed in the literature. This paper investigates several well-known parallel thinning algorithms and proposes a modified version of the most widely used thinning algorithm, called the ZS algorithm. The proposed modified ZS (MZS) algorithm is implemented and compared against seven existing algorithms. Experimental results and performances evaluation, using different image databases, confirm the proposed MZS algorithm improvement over the seven examined algorithms both in terms of the obtained results quality and the computational speed. Moreover, for an efficient implementation (on Graphics Processing Units), a parallel model of the MZS algorithm is proposed (using the Compute Unified Device Architecture, CUDA, as a parallel programming model). Evaluation results have shown that the parallel version of the proposed algorithm is, on average, more than 21 times faster than the traditional CPU sequential version. Keywords Thinning algorithm · Binary image · Thinning rate · Thinning speed · Connectivity · Noise sensitivity · GPU

B

Lynda Ben Boudaoud [email protected]

1

Laboratoire d’Informatique Médicale, LIMED, Faculté des Sciences Exactes, Université de Bejaia, 06000 Bejaia, Algeria

2

Image and Information Processing Department, IMT Atlantique Institute, 29238 Brest, France

1 Introduction In digital image processing, thinning is a fundamental early preprocessing step usually applied in many computer vision applications [18]. Examples include pattern recognition, fingerprint classification, biomedical diagnosis. The thinning process intends to reduce the original binary object containing different thicknesses to a thin representation of a unit width called a skeleton. Thinning is a very useful step which reduces the amount of data to be processed (image complexity) reducing, thus, the time needed for objects recognition and making the shape analysis process easier and more intuitive. In the past few decades, many thinning algorithms have been proposed in the literature [3,6,7,10,13,15–17,20]; most of them are iterative parallel algorithms [12]. Iterative parallel thinning algorithms are based on deleting successive layers on the edges of the considered pattern until a suitable skeleton is obtained in a parallel way. Therefore, all pixels can be independently examined in each iteration [10]. These algorithms can be mainly divided into three categories according to the used approach: the directional approach, the subfield approach and the fully parallel approach [2,9]: – The directional approach: This approach consists on “breaking” the thinning step (iteration) into several thinning sub-steps (sub-iterations) each of which is based on directions, or, on a combination of directions (north, west, east, south). The deletion conditions are changed from one sub-iteration to another, and all border pixels belonging to the same direction are removed in parallel; – The subfield approach: In this approach, the image is partitioned into subfields. At each iteration, a parallel operator is only applied over one member (subfield) of the partition. At a given iteration step, only border pixels (in

123

L. Ben Boudaoud et al.

the active subfield) are potential candidates to be deleted in parallel; – The fully parallel approach: In this approach, no subiteration, or subfield, takes place. The same thinning operator (i.e., deletion criterion) is used at each iteration (using an extended k × k neighborhood, with k ≥ 3). In general, a good thinning algorithm should assume some important properties. Some of these properties are given as follows: – Connectivity (or topology) preservation: the algorithm should preserve the connectivity of the pattern (this notion has been extensively discussed in [2]); – One-pixel thickness: the resulting skeleton should be made of curves, i.e., one-pixel-thick objects; – Shape preservation: the shape of a considered pattern should be preserved. For example, an object having the shape of the capital letter “E” should not be transformed into another object like “F”; – Medial location: the resulting skeleton should lie in the middle of the original shape; – Prevention of excessive erosion: the length of lines and curves in the resulting skeleton should be preserved; – Noise immunity: the slight noise appearing near the boundary of an object should not intensively affect the resulting skeleton. This means that no spurious branches should be contained in the resulting skeleton or, at least, they should be minimized; – High computational speed: the algorithm should be as fast as possible. In fact, increasing the speed of the thinning process will speed up the application using the thinning algorithm. The above-enumerated fundamental properties are the most important features for which thinning methods should consider. Some of these features are, somehow, contradictory so it is hard to fulfill all these requirements simultaneously. For this reason, different algorithms may not satisfy all of these properties. Nevertheless, a good algorithm is assumed to satisfy most of them. In this paper, a modified ZS algorithm, called MZS, allowing to compute the skeleton of binary images and satisfying the entire above good algorithm properties is proposed. A preliminary version of the proposed algorithm has been given in [2]. Indeed, in this previous work, a partial version of the proposed algorithm is presented and compared with the ZS algorithm using only two performance measures (which are the thinning speed and the thinning rate) on a small scope of images. In the present paper, the complete version of the MZS algorithm is given with a demonstration of its validity regarding to connectivity (topology) preservation. The scope of experiments is extended to cover different

123

widely known databases (e.g., KIMIA, DRIVE and FVC). Moreover, a comparison is conducted with seven known algorithms in the literature (based on visual and numerical performance measures). Five performance measures have been used, namely the thinning time, the connectivity measure, the noise sensitivity measure, the thinning speed and the thinning rate. Moreover, a parallel model and implementation of the MZS algorithm (on Graphics Processing Units, GPU, using CUDA) is proposed. This implementation leads to substantial execution time improvements compared to the traditional CPU version (where the thinning process has been accelerated, on average, 21 times). The remainder of this paper is organized as follows. Section 2 gives some basic notions about binary images. In Sect. 3, seven recent and/or most popular parallel thinning algorithms, namely ZS algorithm [20], LW [11], WT [19], AW [1], JK [7], Tarabek [17] and Davit [8], are studied. Section 4 presents the proposed MZS algorithm, and the proof of its topology (connectivity) preservation is shown. Section 5 defines the used performance measurements and, then, presents the results of visual analysis and performances evaluation for the MZS and the considered algorithms on a broad range of test images along (a discussion is also given for each image’s type taken from different databases). Section 6 exhibits the proposed multi-threaded implementation on GPU of the MZS algorithm. Finally, in Sect. 7, we round off the paper by conclusions.

2 Some basic notions A binary image S is represented by a matrix P, of size m × n, where P(i, j) represents the binary value of the pixel at the position (i, j). If P(i, j) = 1 (resp. P(i, j) = 0) then, the pixel (i, j) is black (resp. white). Black pixels form the foreground of S or the object O, and the white ones form the background which is the complement of the foreground O. As it is shown in Fig. 1, a pixel “ p” at a location (i, j) has two horizontal and two vertical neighbors. This set of four pixels, noted N4 ( p), is called the 4-neighborhood of the pixel “ p”. Four diagonal pixels, noted N D ( p), represent the diagonal neighborhood of the pixel “ p”. The points of N4 ( p) and N D ( p) gathered together are called the 8-neighborhood of the pixel “ p” and referred to as N8 ( p)/N8 ( p) = N4 ( p) ∪ N D ( p). Two pixels are said to be 4−, 8− or D− connected if they are adjacent in some sense, i.e., – they are neighbors (N4 ( p), N D ( p) or N8 ( p)), respectively, and – their intensity values are the same (“o” or “1” for binary images).

A modified ZS thinning algorithm by a hybrid approach Fig. 3 The ZS 3 × 3 neighborhood

Fig. 1 Pixel “ p” different neighborhood configuration

Fig. 4 Templates to retain the pixel p1

Fig. 2 Human body (in the left) and its skeleton (in the right)

The application of a thinning process on any input binary image produces another binary image as an output which is the skeleton. An example of a result of thinning the drawn human body is given in Fig. 2.

3 Literature review In this section, a selection of recent and/or widely known and used parallel thinning algorithms is presented. Widely known thinning algorithms are simple and have always guided comparisons of new algorithms, whereas the most recent ones try to tackle as much as possible of the previously expressed properties of a good thinning algorithm. 3.1 The ZS algorithm A very popular and well-proved thinning algorithm is the ZS algorithm proposed by Zhang and Suen [20]. It is an iterative parallel thinning algorithm operating on a 3×3 neighborhood as shown in Fig. 3. The ZS algorithm is a directional algorithm which is broken up into two sub-iterations. The first sub-iteration aims to delete the southeast boundary pixels and the northwest corner pixels, while the second one aims to delete the northwest boundary pixels and the southeast corner pixels (i.e., the opposite orientations).

Fig. 5 MZS algorithm flowchart

123

L. Ben Boudaoud et al. Fig. 6 Number of triangles of P(i, j)

In the first sub-iteration, the contour point p1 is deleted from the pattern, if it satisfies the following conditions: (a) (b) (c) (d)

Fig. 7 Thinning results of a 2 × 2 square given by the eight evaluated algorithms

2 ≤ B( p1 ) ≤ 6 A( p1 ) = 1 p2 × p4 × p6 = 0 p4 × p6 × p8 = 0

In the second sub-iteration, the contour point p1 is deleted from the pattern, if it satisfies the following conditions: (a) (b) (c’) (d’)

2 ≤ B( p1 ) ≤ 6 A( p1 ) = 1 p2 × p4 × p8 = 0 p2 × p6 × p8 = 0

where A( p1 ) or A(P(i, j)) can be computed easily by calculating the number of 0 −→ 1 or 1 −→ 0 in the clockwise direction of the current pixel P(i, j) in the order p2 , p3 , . . . , p2 . B( p1 ), or B(P(i, j)), represents the number of foreground pixels in the 8-neighborhood of P(i, j), that is:

B( p1 ) =

9 

Pi .

Fig. 8 Thinning results of 135◦ line given by the eight evaluated algorithms

(1)

i=2

If any condition is not satisfied then the pixel p1 will not be deleted from the foreground. Although ZS is a relatively simple, fast and an efficient algorithm, it has some drawbacks, mainly [2]: 1. The loss of connectivity due to the complete disappearance of the 2 × 2 square patterns and all the patterns thinned to 2 × 2 square in the thinning process. This makes the algorithm invalid from the topological (connectivity) point of view (this will be shown in Sect. 5.1); 2. The resulting skeleton does not respect the one-pixel thickness owing to the presence of redundant pixels in the skeleton; 3. The ZS algorithm suffers from the excessive erosion of diagonal lines (to be shown ZS result in Sect. 5.1).

123

Fig. 9 Thinning results of 45◦ line given by the eight evaluated algorithms

3.2 The LW algorithm By changing the neighborhood condition, Lu and Wang have proposed an improved ZS algorithm, called LW algorithm, to solve the problem of excessive erosion in slanting line [11], by changing the neighborhood condition. In fact, the LW algorithm retains the same diagonal lines, but they are of two-pixel wide. Moreover, this algorithm also suffers from the problem of the complete disappearance of the 2×2 square (as in the case of the ZS algorithm).

A modified ZS thinning algorithm by a hybrid approach

3.3 The WT algorithm Wu and Tsai have proposed a new parallel thinning algorithm, called WT algorithm, for binary images [19] which can be considered as a fully parallel thinning algorithm. In fact, the authors proposed 14 templates for thinning: twelve 3 × 3 templates, one 3 × 4 and one 4 × 3 templates in the overall set. If the neighborhood of a current black pixel matches any of the templates, then, this pixel will be deleted. The WT algorithm is known as one of the fastest parallel thinning algorithms. It allows obtaining skeletons with the properties of perfectly 8-connected and noise insensitive. However, it does not preserve the structures of original objects in some patterns producing, thus, biased skeletons and cutting corners [6] . 3.4 The AW algorithm Ahmed and Ward have proposed a rotation invariant rulebased thinning algorithm (called AW algorithm) for character recognition [1]. This algorithm can also be considered as a fully parallel one since there is a unique pass applying the same conditions to each pixel at any iteration. The AW algorithm uses 20 rules, or templates, during the thinning process. To avoid discontinuities and extraneous pixels, due to the parallel application of 20 rules, the AW algorithm adds conditions deriving templates. Although, the AW algorithm produces skeletons in their central lines which are rotation invariant, preserving the shape of the original objects, this algorithm, as the ZS algorithm, fails to produce a one-pixel width skeleton. Furthermore, it is restricted to thin symbol and character digital patterns. 3.5 The JK algorithm Jagna and Kamakshiprasad have proposed a new thinning algorithm (called JK algorithm) for binary images. It is a modification of the ZS algorithm aiming to improve it and to solve the excessive erosion problem. The JK algorithm deletes, in both sub-iterations, the southeast boundary pixels and the northwest corner pixels. This algorithm preserves the connectivity of the input image and shows better performances than the ZS algorithm [7]. 3.6 The Tarabek algorithm Tarabek has proposed a robust parallel thinning algorithm for patterns recognition, (called Tarabek algorithm) [17], allowing to overcome the disadvantages of the ZS algorithm, by incorporating (into the ZS algorithm) some additional conditions and, then, applying a post-processing step using an extended 4 × 4 neighborhood. Tarabek algorithm improves

the ZS algorithm, preserves connectivity and is insensitive to boundary noise [17]. 3.7 The Davit algorithm Davit has proposed a slight modification to the ZS thinning algorithm (called Davit algorithm) and reported that Davit algorithm is efficient and better preserves the connectivity for fingerprint images. This algorithm uses the same deletion approach of the ZS algorithm and some conditions have been added in order to improve and to preserve the structure of the image and to stop “unwanted” removal of lines [8].

4 The MZS algorithm In this section, a modified ZS algorithm is proposed. It aims to settle the ZS algorithm’s drawbacks previously expressed, to improve it and to produce better skeletons. The proposed MZS algorithm uses a 3×3 neighborhood as shown in Fig. 3, where – p1 is a black pixel selected for deletion; – B( p1 ) is the number of the black pixels in the 8neighborhood of p1 ; – C( p1 ) is the number of 8-connected black components in the 8-neighborhood of p1 , defined as follows [13]: C( p1 ) = ¬ p2 ∧ ( p3 ∨ p4 ) + ¬ p4 ∧ ( p5 ∨ p6 ) + ¬ p6 ∧ ( p7 ∨ p8 ) + ¬ p8 ∧ ( p9 ∨ p2 )

(2)

The MZS algorithm is based on the directional approach used by the ZS algorithm and is combined with the subfield approach. Two subfields are defined according to the parity of pixels. The first subfield is the set of even pixels (i.e., those for which the sum of their coordinates i + j is even) and the second subfield is similarity composed of the set of odd pixels. The MZS algorithm is described as follows: in the first sub-iteration, the pixel p1 is considered as a candidate for deletion if it satisfies the following conditions: (a) (b) (c) (d) (e)

(i + j)mod2 = 0 C( p1 ) = 1 2 ≤ B( p1 ) ≤ 7 p2 × p4 × p6 = 0 p4 × p6 × p8 = 0

In the second sub-iteration, the pixel p1 is considered as a candidate for deletion if it satisfies the following conditions: (a) (i + j)mod2 = 0

123

L. Ben Boudaoud et al.

(b) (c) (d) (e)

C( p1 ) = 1 1 ≤ B( p1 ) ≤ 7 p2 × p4 × p8 = 0 p2 × p6 × p8 = 0

where C( p1 ) = 1 means that there is only one 8-connected component in the neighborhood of p1 . This implies that it is deletable when p1 is a boundary pixel and the deletion will not disconnect the 3 × 3 neighborhood. With the conditions (c) and (d), it ensures that no hole is created in the 3 × 3 window. 2 ≤ B( p1 ) ≤ 7 means that we are examining a larger set of border points than those considered by the ZS algorithm. This implies the deletion of more boundary pixels. Indeed, the ZS algorithm could have removed more pixels but it was limited by its conditions [2]. To preserve the 2 × 2 square patterns (or in general, the digital patterns that can be reduced to a square), to avoid excessive erosion of diagonal lines and to minimize the generation of spurious branches, an additional condition to the second sub-iteration is integrated. Let p ∗ ∈ N D ( p1 ) be one of the diagonal pixels surrounding p1 which are also odd pixels. During the examination p1 ’s neighborhood, if B( p1 ) = 1, and the black neighbor is a p ∗ pixels, we examine again B( p ∗ ). Then, the pixel p1 is retained if it satisfies B( p ∗ ) ≤ 2. In other words, the pixel p1 is retained if it matches any of the conditioned templates (see Fig. 4), where x can take either 0 or 1 value and B( p ∗ ) ≤ 2. To show the correctness of the proposed MZS algorithm in terms of topology preserving, the following enumerated properties of connectivity preservation (cited in [2]) are demonstrated: 1. Proposition: No object in S can be completely deleted by the MZS algorithm in any iteration. Proof: Let us assume the contrary hypothesis, thus, suppose that an object O is deleted in one iteration, more precisely in only one sub-iteration/subfield. If we consider the neighborhood of any pixel p1 in O, this neighborhood will contain from 0 − 4 ones in the subfield of p1 . The neighborhood with only 0 or 1 ones fails the conditions (c) of the two subfields and the neighborhood with 2, 3 or 4 ones cannot be 8-connected, thus nor deleted, because of failing conditions (b) of the two subfields. Therefore, p1 cannot be deleted, which contradicts the contrary hypothesis. 2. Proposition: The MZS algorithm does not connect any originally disjoint object in O, nor does it disconnect an existing hole or creates any new hole in O. Proof: It is obvious , since the MZS algorithm converts only ones to zeros and only creates new zeros which are

123

four-connected to existing zeros from conditions (d) and (e) of the two sub-iteration/subfield. 3. Proposition: The MZS algorithm does not disconnect any object in O. Proof: It can be clearly seen that the condition (b) ensures that no object can be disconnected in any subfield. Therefore, for a pixel p1 to be a deletion candidate, it must have one 8-connected component in its neighborhood, ensuring, thus, the existence of a path linking the 8 neighboring pixels before deletion. Thus, no path can be disconnected and the MZS algorithm verifies this property. 4. Proposition: The MZS algorithm does not connect any hole of O to another distinct hole or to the background of O. Proof: Let us assume that the contrary hypothesis holds and that two originally unconnected components of the background O become four-connected for the first time at iteration k (via a path W containing a set of one or more ones, T , which were deleted at iteration k). Since no two members of either subfield are four-connected through pixels in the same subfield, none of the members of T are four-adjacent. Thus, in the path W, any member x of T is surrounded by two pixels y and z which were 0 before. For at least one member of T , x, y and z were not four-connected before iteration k, but are four-connected by x in iteration k. Otherwise, the two distinct components O were already connected before iteration k, which contradicts the hypothesis. The flowchart describing the MZS algorithm steps is plotted in Fig. 5.

5 Sequential implementation: tests and results All the algorithms, previously detailed, as well as the MZS algorithm were implemented in C language using PBM images on core i5 intel CPU and, subsequently, tested them on a large number of binary images. This involved several patterns ranging from alphabet of various international languages like Chinese, Arabic, English to different shapes of objects like animals, symbols, handwritten characters, numbers, fingerprints, retinal taken from the DRIVE public retinal images database, the famous Kimia database, the FVC (Fingerprint Verification Competition) database, and our database composed of different shapes. In order to evaluate the MZS algorithm, comparisons among these algorithms are conducted on the basis of the following performance measures: – Thinning Rate (TR): or the degree of thinness of the image. This measure was proposed in [21] by the following formula:

A modified ZS thinning algorithm by a hybrid approach

TR = 1 − TM1 ÷ TM2

(3)

where T M1 stands for the total triangle count of the thinned image given by: TM1 =

n  m 

TC(P(i, j))

S(P(i, j))

(9)

(4) where

and P(i, j) = 1 refers to a black pixel. T C is a function which counts the number of black triangles which can be created from P(i, j) and its neighboring pixels given by the following formula and illustrated in Fig. 6: TC( p1 ) = p1 ∗ [( p8 ∗ p9 ) + ( p9 ∗ p2 ) + ( p2 ∗ p3 ) +( p3 ∗ p4 )]

(5)

TM2 represents the largest number of black triangles that an image could have computed. This number is computed as follows: TM2 = 4 ∗ [max(m, n) − 1]2

 S(P(i, j)) =

1, if A(P(i, j)) > 2, 0, otherwise.

(10)

– Execution Time (ET) or thinning time is the real time taken by the thinning process to thin the input image on a given computer. ET is calculated as the mean of 100 different experiments. – Thinning Speed (TS) This measure intends to evaluate the number of pixels thinned per time unit (second) [2]. It is computed by: DP ET

(11)

DP = OP − SP

(12)

TS = (6)

where m, n stand for the dimensions of the matrix P representing the image. When TR = 1, the image is perfectly thinned. – Connectivity Measure (CM) This measure is used to evaluate the connectivity of the output skeleton. It is defined as the number of end points and discrete points in the thinned images. The closer the number of these points in the final skeleton to that of the original image, the higher degree of connectivity is in the skeleton (with respect to the original image) [21]. This measure is given by the following formula: CM =

SM =

n  m  i=0 j=0

i=0 j=0

n  m 

Hence, the fewer the cross-points, the higher the immunity of the algorithm to noise. The SM is defined as follows [21]:

where – DP or (deleted pixels) is the number of black pixels deleted during thinning; – OP or (object pixels) is the number of black pixels of the image before thinning; – SP or (skeletal pixels) is the number of black pixels; remaining after thinning or the number of pixels in the skeleton. 5.1 Tests and results

S(P(i, j))

(7)

(8)

In this section, both visual and numerical results of the MZS and the seven thinning algorithms on various kinds of shapes (from different databases) are reported. Different algorithms are also evaluated in terms of connectivity, erosion, thinness, noise sensitivity, time execution and speed.

– Sensitivity measure (SM) This measure targets the estimation of the sensitivity to noise. Images usually contain noise as a perturbation in the outline of the object which causes deformation, spurious branches (offshoots and tabs) in the final skeletons. It is obvious that a crosspoint exists wherever there is a spurious branch. For a given image, the total number of its cross-points is fixed. The excess cross-points are caused by the vulnerability to noise and the inadequacy of the thinning algorithm.

Test 1: 2 × 2 square and diagonal lines First, different algorithms are evaluated on the 2 × 2 square, 135◦ and 45◦ diagonal lines. The application results of the different thinning algorithms on these patterns are shown in Figs. 7, 8 and 9 successively. As it can be witnessed from Fig. 7, there are two skeletal pixels in the MZS algorithm results, that are not redundant, and one skeletal pixel in Tarabek algorithm, whereas all the other algorithms (i.e., ZS, LW, WT, AW, JK and Davit) make the structure completely vanishing. Hence, all algorithms fail

i=0 j=0

where  S(P(i, j)) =

1, if B(P(i, j)) < 2, 0, otherwise.

123

L. Ben Boudaoud et al.

Fig. 11 White skeletons superposed on the original turtle produced by the eighth algorithms

Fig. 10 White skeletons superposed on the original duck produced by the eight algorithms

in maintaining the connectivity which is a basic requirement for all thinning algorithms. In Figs. 8 and 9, the MZS algorithm deals successfully with each diagonal pattern. There are neither excessive erosion (the length of the diagonal lines is remain unchanged) as in ZS and Davit algorithms nor redundant pixels as in LW, AW and JK algorithms. Furthermore, it produces better diagonal skeleton shape than those of Tarabek and WT algorithms (which also produce reasonably good result). Test 2: Bundle Images For this experimentation test, different algorithms are evaluated on binary real-life bundle patterns which are: duck, turtle and camel. After the application of different thinning processes, we obtain for each image its corresponding superimposed white skeleton over the black original pattern as shown in Figs. 10, 11 and 12. For the duck image (see Fig. 10), except the algorithms AW, WT and MZS, all other algorithms (i.e., ZS, LW, JK, Tarabek and Davit) do not preserve the shape of the original duck. In fact: – The ZS, Tarabek and Davit algorithms erase one finger and one part of the tail; – The LW algorithm produces the fingers correctly but one part of the tail is missing;

123

– The JK algorithm erases one finger and tail. It deviates from the medial line and thus produces a very bad skeleton which does not represent the original shape of the duck. – The Davit algorithm produces a broken pieces skeleton, it is a disconnected skeleton. It does not preserve the connectivity of the pattern. The AW and WT algorithms reflect the shape of the duck; their skeletons contain some spurious branches (at the level of abdomen and ankle) which is undesirable. It is visually clear that the skeleton produced by the MZS algorithm is better than all the others. It can be clearly observed that the geometric characteristics are completely preserved. Indeed, the paw and the tail are preserved; it takes the central line and there are no spurious branches. For the turtle image (see Fig. 11), the MZS and Tarabek algorithms produce a “clean” skeleton, preserving the medial, topology and geometry of the pattern, whereas other skeletons produced by the other algorithms appear either with spurious branches and redundant pixels like (in the AW, JK, LW, WT and Davit algorithms) or deviating from the center (like the JK algorithm). Moreover, Davit algorithm produces a disconnected skeleton structure. For the camel image (see Fig. 12), the MZS and Tarabek algorithms produce skeletons preserving all the desirable characteristics, whereas other skeletons produced by the other algorithms (like in LW, WT, AW and JK algorithms), they result in thicker curves and spurious branches. Moreover, skeleton produced by the JK algorithm is pulled to the

A modified ZS thinning algorithm by a hybrid approach

and the seven other algorithms and this on the different tests performed previously. It is worthwhile to notice that the best values are shown as black boldfaced values and the worst ones are shown as red boldfaced values. Initial values before thinning are given in column X 0 : Object Points (OP), Thinness Measure (TR0), Connectivity Measure (CM0) and Sensitivity Measure (SM0). After thinning, the performance measures, shown in the third column (PM), are calculated for the various algorithms. From the resulting values, comparisons (see Table 1) show that the MZS algorithm behaves better than most of the other algorithms in terms of:

Fig. 12 White skeletons superposed on the original camel produced by the eight algorithms

left of the original camel shape moving away from the middle. In the case of Davit algorithm, a disconnected camel skeleton is generated. Test 3: Sparse images Fingerprints and retinal blood vessels binary images are also considered for the evaluation of different algorithms. After the application of different thinning processes, the corresponding black skeleton is obtained. Based on these skeletons, it is possible to extract the bifurcation points, the end points and to identify each segment and extract its characteristics. It is worthwhile to notice that the MZS algorithm radiates a clear fingerprint skeletons at one-pixel width despite of the boundary noise contained in the original images, whereas other algorithms suffer from two-pixels width, spurious branches and broken skeleton (surrounded with red, blue and green dashed circles, respectively, as shown in Figs. 13 and 14). The blood vessels final skeletons produced by the evaluated algorithms preserve the connectivity and the shape except Davit algorithm which generates a disconnected segments. Moreover, some algorithms (like the ZS, LW, AW and JK algorithms) produce thicker segments (see Fig. 15). Table 1 shows the performance results according to the five performance measures, previously detailed, for the MZS

– Thinness: For the five images, the Thinness Measure (TR) value produced by all the algorithms is greater than its corresponding TR0, which means that thinning is done by each one. However, it should be noticed that the MZS TR value is the closest to 1 among others. This means that the MZS does not contain any redundant pixels, or, it contains a very small number of redundant pixels leading, thus, to produce thinner and better skeletons of width 1. It is considered as the best thinner; – Execution time (ET): The MZS algorithm is faster than all other algorithms in thinning bundle patterns, but for sparse ones, the WT and LW algorithms seem to be slightly faster than the MZS algorithm (while other algorithms come after). The execution time ET and the speed TS measures confirm this fact; From Table 1, we can clearly notice that the MZS execution time (ET) is the best (lesser value) among other values in the cases of turtle, duck and camel images, whereas in the cases of fingerprint and retinal images, the MZS is slower than the LW and WT algorithms by few milliseconds. On the other hand, the worst execution time is produced by the AW, Tarabek and Davit algorithms (see Table 1). In Fig. 16 and 17, our target is to analyze the behavior of these algorithms (MZS, LW, WT) in terms of execution time versus fragmentation of an object to various parts while maintaining the same number of black points. For this, eight images are considered, where, in the first image, the object is focused in one part, and, in the remaining images, the object is fragmented increasingly. For each image, the execution time, in milliseconds, is computed for the MZS and its competitive algorithms (i.e., WT and LW algorithms) in sparse images. It is worthwhile to notice that the three algorithms behave in the same way. In fact, the more fragmented the object, the less time the algorithms have for thinning. Nevertheless, the MZS algorithm is better for the first three images (Image1, Image2 and Image3) and the more the object is fragmented, the three algorithms show almost similar behavior, and WT and LW become slightly better.

123

L. Ben Boudaoud et al.

Fig. 13 Black fingerprint (1) skeletons given by the eight considered algorithms

123

A modified ZS thinning algorithm by a hybrid approach Fig. 14 Black fingerprint (2) skeletons given by the eight considered algorithms

From this analyze, we can conclude that the more boundary pixels are present in the image, the more the LW, WT and MZS algorithms become faster. The WT and LW algorithms gain time when the number of boundary pixels

increases (because processing boundary pixels implies limited rules and conditions, which in turns, leads to less processing and thus, less consuming time). This is unlike the MZS algorithm which deletes an extended number of

123

L. Ben Boudaoud et al.

Fig. 15 Black retinal skeletons given by the eight evaluated algorithms

123

A modified ZS thinning algorithm by a hybrid approach

Fig. 16 The fragmented images Fig. 18 Thinning time behavior of the competitive algorithms (LW, WT and MZS) according to the number of object pixels in different bundle images

Fig. 17 Thinning time behavior of the competitive algorithms (LW, WT and MZS) according to different fragments

boundary pixels in order to produce thinner skeleton. Figure 18 shows the behavior of the competitive algorithms in terms of execution time (according to the number of object pixels contained in different bundle images). It is important to notice that the thinning time (for the three algorithms) increases as the number of object points contained in the bundle patterns increases. The MZS algorithm needs less time, slightly better than the LW and WT algorithms. Nevertheless, when the number of object points increases in the bundle images, the MZS algorithm becomes the fastest and the time gain achieves thousands of milliseconds as in the case of rabbit image (21643). – Connectivity: It is obvious that all the thinning algorithms produce a connectivity value CM greater than the initial C M0 (because the thinned skeleton would have more end and discrete pixels than those of the original image). From Table 1, the worst connectivity values are produced by Davit algorithm, where its CM value is the greatest one compared to others except the case of fingerprint2 (where

it is surpassed by the WT algorithm values). These values coincide with the visual analysis. Davit thinning leads to patterns disconnections and thus, the generation of more end pixels. In fingerprint2, the WT algorithm skeleton results in more spurious branches, producing thus, more end pixels. Best values of CM (which represents the closest values to its CM0) are produced, in most cases, by the MZS and LW algorithms, representing the higher degree of connectivity. According to the CM values of the overall tests, the proposed MZS algorithm provides a very good connectivity; – Sensitivity to noise: It is obvious that SM values will be greater than SM0 when operating the thinning process (because of the creation of cross-points in the skeleton). The proposed algorithm shows a good insensitivity to noise. The noise sensitivity measure SM of the MZS algorithm (over different experiments) shows better immunity values. From a visual analysis point of view, this corresponds to few or no spurious branches. However, other algorithms like Davits, produce better SM values in the case of duck and camel at the expense of connectivity performance. AW, LW algorithms are the most sensitive to noise producing the worst SM values in most experiments. Table 2 summarizes the visual results given in Figs. 7, 8, 9, 10, 11, 12, 13, 14 and 15 and the numerical results given in Table 1. Three symbols are used to classify each algorithm according to the six characteristics of a good thinning algorithm, where (−): means that the algorithm does not give a satisfactory result for the corresponding characteristic.

123

L. Ben Boudaoud et al. Table 1 Comparing thinning performances of the eight evaluated algorithms Images Turtle

Duck

Camel

Finger1

Finger2

X0

ZS

LW

WT

AW

JK

Tarabek

Davit

MZS

O P = 4568

SP

174

190

178

188

192

164

159

171

T R0 = 0.651248

TR

0.999612

0.999664

0.999920

0.999472

0.999728

1.00

0.999696

1.00

C M0 = 0

E T (ms)

7.595

6.947

8.405

13.679

6.261

8.043

11.286

6.179

S M0 = 0

T S ( p/s)

578530.501

630201.476

522327.453

320187.317

698907.651

547574.792

390643.871

711590.761

CM

4

5

5

4

5

4

7

4

SM

2

3

3

3

3

2

1

2

O P = 21643

SP

570

596

544

569

564

525

518

538

T R0 = 0.594783

TR

0.999612

0.999566

0.999958

0.999723

0.999754

0.999962

0.999608

0.999981

C M0 = 0

E T (ms)

99.125

97.663

120.460

177.711

96.123

105.916

94.199

81.710

S M0 = 0

T S ( p/s)

212589.096

215505.508

175153.296

118585.738

219292.127

199384.247

224258.936

258290.515

CM

4

4

7

8

6

4

29

3

SM

4

5

6

8

6

4

3

6

O P = 6900

SP

417

521

453

506

559

369

389

381

T R0 = 0.495902

TR

0.998332

0.997762

0.999639

0.997807

0.998873

0.99991

0.998332

0.99999

C M0 = 0

E T (ms)

13.2

12.637

15.507

24.53

17.019

14.081

20.02

10.286

S M0 = 0

T S ( p/s)

491118.558

504803.584

415749.739

260655.205

372587.139

457139.738

327413.129

549413.071

CM

9

8

13

12

11

9

22

9

SM

7

10

11

11

10

7

6

7

O P = 37538

SP

9115

9302

8740

9079

9181

8610

9058

8487

T R0 = 0.769854

TR

0.997671

0.998084

0.999840

0.998060

0.998234

0.999908

0.997738

0.999963

C M0 = 172

E T (ms)

34.461

28.385

24.827

39.013

29.538

41.876

31.544

27.548

S M0 = 28

T S ( p/s)

824775.890

1112302.613

1159923.955

729471.910

960002.586

690804.484

902859.869

1054544.379

CM

303

269

320

300

302

304

348

275

SM

28

53

52

34

29

28

27

26

O P = 46171

SP

10125

10401

9114

10065

10223

8933

10088

8916

T R0 = 0.892029

TR

0.998055

0.998402

0.999908

0.99818

0.998478

0.999962

0.998075

0.999977

C M0 = 10

E T (ms)

57.939

35.805

28.24

41.969

37.209

49.955

80.699

38.375 970803.813

S M0 = 2

Retinal1

PM

T S ( p/s)

622142.296

999010.407

1312218.209

860305.974

966111.51

745425.218

447295.788

CM

124

104

155

119

120

124

143

120

SM

44

80

80

53

44

44

43

42

O P = 24658

SP

8976

9061

7803

8945

9030

7757

8863

7736

T R0 = 0.940945

TR

0.997645

0.998033

0.999712

0.997785

0.998055

0.999857

0.997718

0.999878

C M0 = 82

E T (ms)

49.387

33.299

34.017

43.919

47.095

48.432

63.369

36.499

S M0 = 15

T S ( p/s)

316895.598

471970.384

483844.565

359080.845

329166.493

350670.634

249253.964

456355.791

CM

124

96

124

109

121

124

206

120

SM

192

192

186

191

177

200

182

194

Bold black values correspond to the best values for the current measure. Bold red values correspond to the worst values for the current measure

(+): means that the algorithm gives a satisfactory result for the corresponding characteristic. (++): means that the algorithm gives more satisfactory result for the corresponding characteristic.

6 GPU parallel implementation of the MZS algorithm Today’s world everything goes around how to get fast solutions. The thinning techniques are fundamental operations

123

and important building blocks of many image processing algorithms. If these basic operations could be implemented faster, then all the applications using these operations will get automatically faster. One of the solutions allowing the acceleration of the image thinning process consists on the use of GPUs which are inexpensive parallel architectures. It is worthwhile to recall the fact that not all image processing algorithms are suitable for parallelization on GPU. However, the analysis of the MZS algorithm shows its intrinsic characteristic of data parallelism. It performs the same computations on black pixels at each sub-iteration, which

A modified ZS thinning algorithm by a hybrid approach Table 2 Comparison of results of different algorithms

Algorithms

ZS

LW

WT

AW

JK

Tarabek

Davit

MZS

1Pixel-width

−

−

−

−

−

+

−

++

Connectivity Preservation

−

−

−

−

−

++

−

++

Excessive erosion

−

−

+

+

−

−

−

++

Medial location

+

+

+

+

−

+

+

++

Noise Sensitivity

+

−

−

−

+

+

+

++

Speed

−

++

++

−

−

−

−

++

makes it a good candidate to fit the GPU architecture. In this section, the Compute Unified Device Architecture (CUDA) technology and the GPU architecture are briefly described. The results obtained by the parallel implementation of the MZS algorithm compared to its serial implementation are, then, detailed. 6.1 CUDA (Unified Device Architecture) CUDA is a parallel computing architecture and programming model initially developed by NVIDIA [14]. It enables important increase in terms of computing performances by exploiting the advantage of the Graphics Processing Unit (GPU) power using high-level languages (like C language). The CUDA software architecture consists of a host program that runs on the CPU and one or more parallel kernels running on the GPU device.

Fig. 19 General architecture of a GPU and its memory hierarchy (Note that the actual layout is much more complex and differ for each GPU)

6.2 GPU architecture A GPU architecture is generally composed of several Graphic Multiprocessors (MPs) ranging from 2 to 32 in some powerful graphic cards; each comprising a set of Streaming Processors (SPs). It has also different types of memory (i.e., global, constant, texture and local). Operations on the global memory can be read/write from all threads but each thread has only access to its private local memory, while the access to the constant and texture memory is on read only [4,5]. Each MP has a shared memory that allows a read/write access by all the SP it contains while all the MPs have access to the global/constant/texture memory of the device (see Fig. 19). The GPU architecture belongs to the family of SIMD parallel architectures (Single Instruction Multiple Data). In other words, in a GPUs architecture, all threads execute the same instructions. 6.3 Tests and results The parallel version of the MZS algorithm using CUDA (Compute Unified Device Architecture) was implemented on NVIDIA GeForce GTX 950M graphic card.

Fig. 20 MZS thinning execution time over six different images (with different sizes) on CPU and GPU

Several experiments were conducted to compare the MZS implementations performances of GPU and CPU. In the first experiment, as it is shown in Fig. 20, the parallel MZS algorithm was evaluated on the same images that have been previously been used in the sequential context and a comparison of the thinning times was conducted. As it can be clearly witnessed, the parallel version shows an interesting thinning execution time (i.e., much faster than the CPU implementation) over different images of different

123

L. Ben Boudaoud et al.

Fig. 21 The MZS thinning execution time over different sizes of the duck image on CPU and GPU

Fig. 22 The MZS thinning execution time over different sizes of the letter “L” image on CPU and GPU

Fig. 23 MZS thinning speedup evolution on different duck image sizes

Fig. 24 MZS thinning speedup evolution on different alphabet letter “L” image sizes

7 Conclusions sizes. In the second experiment, we have opted to use the same image by scaling its size to achieve 16 million pixels and then compare the thinning time and speedup given by the two implementations. Figures 21 and 22 show the thinning execution time spent by the two CPU and GPU implementations over the same image with a scaling by factor 2 on the duck image and the alphabet letter “L”. Notice that the serial MZS thinning time grows with the size of the image achieving hundreds of seconds. The process becomes very slow when the image size exceeds a million of pixels. In the counterpart, the parallel MZS on the GPU presents a substantial thinning time improvement. Figures 23 and 24 show the acceleration results of the parallel MZS algorithm on GPU corresponding to the previous experiment (scaling duck and alphabet letter) achieving 25 times faster than the CPU-MZS.

123

In this paper, a modified version of the ZS image thinning algorithm based on combining its directional approach with a subfield approach is proposed. The proposed algorithm was compared against seven recent and/or popular image thinning algorithms on a variety of images of different shapes including characters, symbols, animals, fingerprints, etc., extracted from several databases (KIMIA, DRIVE, FVC), independently of the application area. Considered algorithms fulfill some properties of a good image thinning algorithm and fail fulfilling some other properties. However, the performance evaluation results and experiments show that the proposed MZS algorithm strongly satisfies most requirements leading, thus, to believe that the proposed algorithm can be considered as an effective and a robust parallel image thinning algorithm. Indeed, the proposed algorithm:

A modified ZS thinning algorithm by a hybrid approach

– Produces perfectly connected thinned skeletons having single pixel width (this accurately reflects the shape of the original pattern); – Is not affected by boundary noise in a short time. Furthermore, the good performances and the visual quality of the obtained results lead to consider the proposed algorithm as being outperforming existing image thinning algorithms. Nevertheless, the MZS thinning on the CPU becomes time-consuming whenever the image size is increased. This makes the performances of the MZS algorithm strongly dependent on the considered application. The parallel implementation of the MZS algorithm on a GPU using CUDA framework presents a substantial time improvement compared with the sequential implementation with a speedup factor that may reach the value 21. As a perspective, the application of the proposed MZSbased GPU algorithm into fingerprint and retinal identification systems is actually under study in order to compare the performance accuracy of matching with other existing methods. Acknowledgements This work was carried out in the framework of research activities of the laboratory LIMED which is affiliated to the Faculty of Exact Sciences of the University of Bejaia and the Image and Information processing department of IMT Atlantique Institute. The authors would like to thank the referees for their comments.

References 1. Ahmed, M., Ward, R.: A rotation invariant rule-based thinning algorithm for character recognition. IEEE Trans. Pattern Anal Mach Intell 24(12), 1672–1678 (2002) 2. Boudaoud, L.B., Sider, A., Tari, A.: A new thinning algorithm for binary images. In: Control, Engineering and Information Technology (CEIT), 2015 3rd International Conference on, pp. 1–6. IEEE (2015) 3. Chen, W., Sui, L., Xu, Z., Lang, Y.: Improved Zhang-Suen thinning algorithm in binary line drawing applications. In: Systems and Informatics (ICSAI), 2012 International Conference on, pp. 1947–1950. IEEE (2012) 4. Cheng, J., Grossman, M., McKercher, T.: Professional Cuda C Programming. John Wiley & Sons, London (2014) 5. Couturier, R.: Designing Scientific Applications on GPUs. Chapman and Hall/CRC, Boca Raton (2013) 6. Deng, W., Iyengar, S.S., Brener, N.E.: A fast parallel thinning algorithm for the binary image skeletonization. Int. J. High Perform. Comput. Appli. 14(1), 65–81 (2000) 7. Jagna, A., Kamakshiprasad, V.: New parallel binary image thinning algorithm. ARPN J. Eng. Appl. Sci. 5(4), 64–67 (2010) 8. Kocharyan, D.: An efficient fingerprint image thinning algorithm. Am. J. Softw. Eng. Appl. 2(1), 1–6 (2013) 9. Kong, T.Y., Rosenfeld, A.: Topological Algorithms for Digital Image Processing, vol. 19. Elsevier, Amsterdam (1996) 10. Lam, L., Lee, S.W., Suen, C.Y.: Thinning methodologies-a comprehensive survey. IEEE Trans. Pattern Anal. Mach. intell. 14(9), 869–885 (1992)

11. Lü, H., Wang, P.S.P.: A comment on a fast parallel algorithm for thinning digital patterns. Commun. ACM 29(3), 239–242 (1986) 12. Palágyi, K.: A 3-subiteration 3D thinning algorithm for extracting medial surfaces. Pattern Recognit. Lett. 23(6), 663–675 (2002) 13. Rosenfeld, A.: Digital Picture Processing. Academic Press, New York (1976) 14. Sanders, J., Kandrot, E.: CUDA by Example: An Introduction to General-Purpose GPU Programming. Addison-Wesley, Upper Saddle River (2010) 15. Song, C., Pang, Z., Jing, X., Xiao, C.: Distance field guided L 1 median skeleton extraction. Vis. Comput (2016). doi:10.1007/ s00371-016-1331-z 16. Tang, Y.Y., You, X.: Skeletonization of ribbon-like shapes based on a new wavelet function. IEEE Trans. Pattern Anal. Mach. Intell. 25(9), 1118–1133 (2003) 17. Tarabek, P.: A robust parallel thinning algorithm for pattern recognition. In: Applied Computational Intelligence and Informatics (SACI), 2012 7th IEEE International Symposium on, pp. 75–79. IEEE (2012) 18. Wolberg, G.: Skeleton-based image warping. Vis. Comput. 5(1), 95–108 (1989) 19. Wu, R.Y., Tsai, W.H.: A new one-pass parallel thinning algorithm for binary images. Pattern Recognit. Lett. 13(10), 715–723 (1992) 20. Zhang, T., Suen, C.Y.: A fast parallel algorithm for thinning digital patterns. Commun. ACM 27(3), 236–239 (1984) 21. Zhou, R., Quek, C., Ng, G.S.: A novel single-pass thinning algorithm and an effective set of performance criteria. Pattern Recognit. Lett. 16(12), 1267–1275 (1995)

Lynda Ben Boudaoud is a PhD student at the university of Bejaia. She got her Master and Licence diploma from the University of Bejaia (Algeria) in 2011 and 2013, respectively. Her research activities revolve around high-performance computing paradigms, programming and applications in imaging, applied mathematics and information security. She is a member of the Medical Image Processing group (TIM) at the LIMED laboratory (Laboratoire d’Informatique Médicale).

Basel Solaiman received the Telecommunication Engineering degree from the Ecole Nationale Supérieure des Télécommunications (ENST) de Bretagne, in 1983, and the D.E.A, Ph.D. and H.D.R degrees from the Université de Rennes-I, France, in 1983, 1988 and 1997, respectively. In 1983, he joined the Laboratoire d’Electronique de Philips, Paris, France, where he worked on image compression standards. In 1992, he joined the ENST de Bretagne (TELECOM-

123

L. Ben Boudaoud et al. Bretagne), Brest, where he is actually a Full Professor and the Head of the Image and Information Processing Department. He has published several academic books, over 180 journal papers, several book chapters and over 230 conference proceedings and technical reports. He chaired several international conferences focusing on his research activities. Application domains of his research activities range from civilian into military ones: medical, remote sensing, underwater imaging and knowledge mining.

Abdelkamel Tari is a Full Professor and Head of LIMED (laboratory of Medical Computing) since February 2013 and team leader of data mining. He held different positions at the university of Bjaia (Head of Computing department, Operational Research department, Head of Doctoral School) and at the National School of Computing ESI (ex.INI), Algiers between 1986–1994. After his Diploma in Mathematics at USTHB (Algiers), he held his

123

Master by Research in Lancaster University (England). He defended his PhD Thesis in Computer Sciences in 2008 at the university of Bjaia and his accreditation to supervise (HDR) at ESI (ex-INI) in 2010.