A parallel algorithm for skeletonizing images by using spiking neural P systems

A parallel algorithm for skeletonizing images by using spiking neural P systems Daniel Dı ´az-Pernil a, n , Francisco Pen ˜a-Cantillana b , Miguel A. Gutie ´ rrez-Naranjo b a CATAM Research Group – Department of Applied Mathematics I University of Seville, Spain b Research Group on Natural Computing – Department of Computer Science and AI, University of Seville, Spain Keywords: Spiking neural P systems Skeletonizing Membrane computing Guo and Hall algorithm abstract Skeletonization is a common type of transformation within image analysis. In general, the image Bis a skeleton of the black and white image A, if the image Bis made of fewer black pixels than the image A, it does preserve its topological properties and, in some sense, keeps its meaning. In this paper, we aim to use spiking neural P systems (a computational model in the framework of membrane computing) to solve the skeletonization problem. Based on such devices, a parallel software has been implemented within the Graphics Processors Units (GPU) architecture. Some of the possible real-world applications and new lines for future research will be also dealt with in this paper. 1. Introduction Computer vision [1] is probably one of the most promising challenges for computer scientists in the near future. This growing research area needs contributions from many other scientific fields such as artificial intelligence, pattern recognition, signal processing, neurobiology, psychology or image processing. Computer vision deals with the automated processing of images from the real world in order to extract and interpret pieces of information. Under a computational point of view, a digital image is a function from a two-dimensional surface; each point in the surface is mapped into a set of features such as brightness or colour. The computational treatment of such mappings (or digital images) is the basis of many current applications in computer vision such as optical character recognition (OCR), biometrics, automotive safety or surveillance. In this paper, we focus on the problem of skeletonizing an image. Skeletonization is one of the approaches for representing a shape with a small amount of information by converting the initial image into a more compact representation and keeping the meaning features. The conversion should remove redundant information, but it should also keep the basic structure. Skeletonization is usually considered as a pre-process in pattern recognition algorithms, but its study is also interesting by itself for the analysis of line-based images as texts, line drawings, human fingerprints or cartography. Many problems in the processing of digital images have features which make it suitable for techniques inspired by nature. One of them is that the treatment of the image can be parallelized and locally solved. Regardless how large is the image, the process can be performed in parallel in different local areas of it. Another interesting feature is that the local information needed for a pixel transformation can also be easily encoded in the data structures used in natural computing. In the literature, we can find many examples of the use of natural computing techniques for dealing with problems associated to the treatment of digital images. One of the classic examples is the use of cellular automata [2]. Other efforts are related to artificial neural networks as in [3]. In this paper, we use spiking neural P systems, a computational model in the framework of membrane computing. Spiking neural P systems (SN P systems, for short) were introduced in [4] as a new class of distributed and parallel computing devices, inspired by the neurophysiological behavior of neurons sending electrical impulses (spikes) along axons to other neurons. SN P systems are the third model of computation in the framework of membrane computing, 1 together with the cell-like model inspired by the compartmental structure and functioning of a living cell and the tissue-like model, based on n Corresponding author. Tel.: þ34 95 455 69 21. E-mail addresses: [email protected] (D. Dı ´az-Pernil), [email protected] (F. Pen ˜a-Cantillana), [email protected] (M.A. Gutie ´rrez-Naranjo). 1 We refer to [5] for a comprehensive presentation in this area and the P system web page http://ppage.psystems.eu, for the up-to-date information. intercellular communication and cooperation between cells in a tissue. Recently, membrane computing techniques have been used for solving problems from digital image. Different P system models have been used for dealing with images, as in [6] where cell-like P systems are used [7–10], where tissue-like P systems are used [11], where the symmetric dynamic programming stereo (SDPS) algorithm [12] for stereo matching was implemented by using simple P modules with duplex channels or even [13], where membrane computing and quantum-inspired evolutionary algorithms are combined. The choice of SN P systems for implementing our bio-inspired solution have some advantages and drawbacks. On the one hand, the use of the flow of spikes among neurons for encoding information allows to use an alphabet with a unique symbol and then, the theoretical model can be translated into a hardware programming language in a natural way. On the other hand, this neurophysiological inspiration has a computational cost, and other bio-inspired implementations obtain better experimental times. The main contribution of this paper comes from the theoretical side, by presenting the design of a classical skeletonizing algorithm in a novel bio-inspired computational model. Nonetheless, the theoretical advantages of the membrane computing techniques for computer vision need a powerful software and hardware for an effective implementation. In this paper, we present a parallel software developed by using a device architecture called Compute Unified Device Architecture (CUDA TM ). It is a general purpose parallel computing architecture that allows the parallel NVIDIA 2 Graphics Processors Units (GPUs) to solve many complex computational problems in a more efficient way than on a Central Processing Unit (CPU). GPUs constitute nowadays a solid alternative for high performance computing, and the advent of CUDA allows programmers a friendly model to accelerate a broad range of applications. The way GPUs exploit parallelism differs from multi-core CPUs, which raises new challenges to take advantage of its computing power. GPU is especially well-suited to address problems that can be expressed as data-parallel computations. The paper is organized as follows: First, we present the restricted model of SN P systems used in the paper and recall the Guo and Hall algorithm for skeletonizing images. In Section 4, the design of the SN P system for skeletonizing images is presented. Next, we show illustrative examples of the use of our implementation. Finally, Section 6 is dedicated to conclusions and future work. 2. Spiking neural P systems SN P systems can be viewed as an evolution in membrane computing corresponding to a shift from cell- to neural-like architectures. In SN P systems, the processing elements are called neurons and are placed in the nodes of a directed graph, called the synapse graph. The computation is performed by sending electrical impulses among the neurons through the synapses. Such electrical impulses are encoded via a single object type, namely the spike, which is placed in the neurons. The number of copies of such object determines the electrical charge of the neuron. Each neuron may also contain rules which allow to send spikes (possibly with a delay) to other neurons, or to remove a given number of spikes from it (firing and forgetting rules). Firing rules allow a neuron to send information to other neurons in the form of electrical impulses which are accumulated at the target cell. Forgetting rules remove from the neuron a predefined number of spikes. The application of every rule is determined by checking the contents of the neuron against a regular set associated with the rule. In each time unit, if a neuron can use a rule, then it must be used. If two or more rules could be applied, i.e., there are two or more rules associated to the neuron and the current number of spikes satisfies the conditions of such rules, then only one of them is applied. The rule to be applied is nondeterministically chosen among all the candidates. Thus, the rules are used in the sequential manner in each neuron, but neurons work in parallel with each other. A global clock is assumed, marking the time for the whole system, and hence the functioning of the system is synchronized. Other biological features have been explored in the framework of SN P systems. One such extension (with mathematical motivation) was introduced in [14], where a neuron can emit more than one spike, if the number of emitted ones is not greater than the consumed ones. Other variants including astrocytes [15], weights, which modify the number of spikes that arrives to a neuron according to the quality of the link between neurons [16–18], anti-spikes [19], or neuron division [20] have also been considered. In this paper, we will consider SN P systems with weights in the synapses and rules without delay. 3 In this way, the restricted model of SN P systems used in this paper can be formally described as follows. A spiking neural P system of degree mZ1 is a construct of the form: P ¼ðO, s 1 , s 2 ,..., s m ,syn,inÞ, where 1. O¼fagis the singleton alphabet (ais called spike); 2. s 1 , s 2 ,..., s m are neurons, of the form s i ¼ðn i ,R i Þ,1rirm, where (a) n i Z0 is the initial number of spikes contained in s i ; (b) R i is a finite set of rules of the following two forms: (1) Firing rules E=a c -a d , where Eis a regular expression 4 over a, and cZdZ1 are integer numbers; if E¼a b (with ban integer number, bZc), then the rule is usually written in the following simplified form: a b =a c -a d ; (2) Forgetting rules a b =a c - l , for band cinteger numbers with bZcZ1. 3. synDf1,2,...,mgf1,2,...,mg W is the set of synapses between neurons, where the set of weights is W ¼ Q \½0,1, i.e., the set of rational numbers between 0 and 1. The synapses also verifies that ði,i,kÞ=2syn for 1rirm,kA W . 4. in is the label of the input neuron of P . A firing rule E=a c -a d AR i can be applied in neuron s i if it contains bspikes, bZc, and a b belongs to the language associated to E. For applying rules of type a b =a c -a d , the neuron must contain exactly bspikes. The execution of these rules removes c spikes from s i (thus leaving the remaining spikes in s i ), and sends dspikes to all the neurons s j such that ði,j,kÞAsyn. The number of spikes that arrives to the neuron jthrough the synapse ði,j,kÞ depends on the number dof emitted spikes and the quality of the synapse, encoded by the weight k. In each neuron, the number of spikes after a computation step is the number of 2 http://www.nvidia.com. 3 For a more general description, see [4]. 4 Along this paper, we use regular expressions of type a n a n , with nAN. In this case, the language associated to a n a n is the set fa nþk 9kANg. For more details about regular expressions, see, for example, [21]. non-consumed spikes plus the contribution of other neurons via the corresponding synapses. Let us consider that dspikes are emitted through a synapse ði,j,kÞ. The contribution of s i to s j is an increase of bdkcspikes, where bvcdenotes the largest integer not greater than v.Aforgetting rule a b =a c - l can be applied in neuron s i if it contains exactly b spikes. The execution of this rule simply removes all the cspikes from s i (thus leaving bcspikes). Aconfiguration of the system is described by the numbers /n 1 ,n 2 ,...,n m Sof spikes present in each neuron. At the beginning of the computation, the number of spikes in each neuron s i is n i , but in the input neuron (with label in): If the input of the computation is N, then, in the initial configuration, the number of spikes in the input neuron is n in þN. Example 1. Let us consider the SN P system P ¼ðfag, s 1 , s 2 , s 3 ,syn,1Þshown in Fig. 1, with s 1 ¼ð7,R 1 Þ, s 2 ¼ð3,R 2 Þ, s 3 ¼ ð8,R 3 Þand the set of rules and synapses: R 1 ¼fR 11 a 5 a n =a 5 -a 3 g R 2 ¼fR 21 a 5 =a 3 -a 2 g R 3 ¼fR 31 a 8 =a 3 - l ,R 32 a 7 =a 2 -a 1 g syn ¼fð1,2,0:7Þ,ð2,1,1Þ,ð2,3,0:8Þ,ð3,2,0:5Þ,ð1,3,0:4Þg We will consider the input N¼2, so, in order to start the computation, 9 spikes (7þ2) are placed in the neuron 1. Notice that R 11 will be applied if the neuron 1 contains at least 5 spikes. R 31 is the unique forgetting rule in the SN P system. Since the input is N¼2, the initial configuration is C 0 ¼/9,3,8S. From this initial configuration, rules R 11 and R 31 can be applied. The rule R 11 consumes 5 spikes from the neuron 1 and sends three spikes to the neurons 2 and 3. These three sent spikes are multiplied by the corresponding weights before arriving to the target neurons. The weight of the synapses between the neurons 1 and 2 is w 12 ¼0:7, so the application of the rule R 11 produces an increase of b30:7c¼2 spikes in the neuron 2. Bearing in mind that w 13 ¼0:4, the application of the rule R 11 increases in b3 0:4c¼1 the number of spikes in the neuron 3. The forgetting rule R 31 deletes three spikes from neuron 3 and hence, the new obtained configuration is C 1 ¼/4,5,6S. Now, the unique applicable rule is R 21 . This rule consumes three spikes from neuron two and sends two spikes to the neurons 1 and 3. According with the corresponding weights, the number of the spikes in the neuron 1 is increased in b21c¼2 and the number of the spikes in the neuron 3 is increased in b20:8c¼1. By applying these modifications, the obtained configuration is C 2 ¼/6,2,7S. From this configuration, rules R 11 and R 32 are applicable. The effects of R 11 have been described above. The rule R 32 removes two spikes from the neuron 3 and sends one spike to neuron 2. This spike is multiplied by the corresponding weight, w 32 ¼0:5, so it does not produce any increase in the number of spikes in neuron 3, since b0:51c¼0. With these changes, the obtained configuration is C 3 ¼/1,4,4S. No more rules can be applied and C 3 is the halting configuration. 3. Guo and Hall algorithm Skeletonization is a common transformation in image analysis. The concept of skeleton was introduced by Blum in [22], under the name of medial axis transform. There are many different definitions of the skeleton of a black and white image and many skeletonizing algorithms, 5 but in general, the image Bis a skeleton of the image A, if it has fewer black pixels than A, preserves its topological properties and, in some sense, keeps its meaning.Fig. 2 illustrates this idea. In this paper, we focus on an iterative procedure of thinning: roughly speaking, the border black pixels are removed as long as they are not considered significant. The remaining set of black pixels is called the skeleton. Among the parallel algorithms, special attention deserves the socalled 1-subcycle parallel algorithms or fully parallel algorithms [24]. Our bio-inspired design is based on a classical skeletonizing algorithm, the Guo and Hall algorithm [24,25]. In this algorithm, the pixels are examined for deletion in an iterative process. First of all, given an pqimage, it is divided into two subsections. One of the sections is composed of the pixels a ij such that iþjis even. Alternatively, the second sub-section corresponds to the pixels a ij such that iþjis odd. The algorithm consists of two sub-iterations where the removal of redundant pixels from both sub-sections are alternated, i.e., in each step only the pixels of one of the subsections are evaluated for its deletion. The decision is based on a 3 3 neighborhood. Given a pixel P0, a clockwise enumeration P1,...,P8 of its eight neighbor pixels is considered, as shown in Fig. 3(a). As usual, for each iAf1,...,8g, Pi is considered as a Boolean variable, with the truth value 1 if Pi is black and 0 if Pi is white. In order to decide if a pixel P0 is deleted in the corresponding iteration sub-cycle, two parameters are evaluated: BðP0Þ¼ X i¼8 i¼1 Pi Fig. 1. Example: the input neuron 1 has 9¼7þ2 spikes at the starting configuration. Fig. 2. A hand-written word and its skeletonization. This skeletonization and the ones shown in Figs. 7 and 8 have been obtained by using the CUDA implementation of the corresponding SN P system. 5 A detailed description is out of the scope of this paper. For a survey, see e.g., [23]. CðP0Þ¼ð:P24ðP33P4ÞÞþð:P44ðP53P6ÞÞ þð:P64ðP73P8ÞÞþð:P84ðP13P2ÞÞ BðP0Þcounts how many pixels in the neighborhood of P0 are black. CðP0Þevaluates the connectivity of the pixel P0. Notice that for isolated black pixels, the connectivity is 0, and for pixels surrounded by eight black pixels, the connectivity is 4. According to the Guo and Hall algorithm, in each iteration, an evaluated black pixel P0 is deleted (changed to white) if and only if all of the following conditions are satisfied. Guo and Hall conditions: 1. BðP0Þ41; 2. CðP0Þ¼1; This condition is necessary for preserving local connectivity when Pis deleted. 3. ðP14P34P54P7Þ3ðP24P44P64P8Þ¼FALSE;Intuitively, this condition is satisfied if P0 is not the central pixel of a cross. For example, let us consider as P0 the central pixel in the image of Fig. 3(b). In this case, BðP0Þ¼341, CðP0Þ¼1, and the third condition is also satisfied. Hence, P0 will be deleted in the corresponding sub-cycle iteration. 4. SN P systems for skeletonizing In this paper, we will show how to use SN P systems for skeletonizing images. In particular, we use SN P systems for implementing the Guo and Hall algorithm. Without losing generality, we will consider each image as a mapping I:f1,...,pg f1,...,qg-fblack,whiteg, with Ið1,kÞ¼Iðp,kÞ¼Iðj,1Þ¼Iðj,qÞ¼white for all kAf1,...,qgand jAf1,...,pg, i.e., the image has pqpixels and all the pixels on the border are white. Given a pixel (i,j), we can use the enumeration of the pixels used in the previous section to represent the neighborhood of the pixel P0in(i,j). Such a neighborhood will be represented as a list ½H0,...,H8, where for rAf0,...,8g,Hr¼1ifPr is a white pixel and Hr¼0ifPr is a black one. 6 This representation of the neighborhood can be done in a more compact way, by encoding the neighborhood as a number 7 in f0,...,511g: codði,jÞ¼ X 8 r¼0 Hr 2 r For example, in Fig. 3(b), the 3 3 neighborhood can be encoded as [0,0,0,0,1,1,1,1,1], or, shortly, 2 4 þ2 5 þ2 6 þ2 7 þ2 8 ¼ 496. Since the decision of removing a black pixel (changing to white) depends on its 3 3 neighborhood and there is a bijective correspondence among the sets of all the possible neighborhoods and the possible encodings f0,...,511g, it is easy to check that the pixel in (i,j) must be removed if it belongs to the set DEL ¼ 6 1214182426283036384446 48 50 56 58 60 62 66 72 74 96 98 104 106 112 114 120 122 124 126 132 134 140 142 144 146 152 154 156 158 164 166 172 174 176 178 184 186 188 190 192 194 200 202 224 226 232 234 240 242 248 250 252 258 262 264 266 270 286 288 290 294 296 298 302 318 384 386 390 392 394 398 414 416 418 422 424 426 430 448 450 454 456 458 462 480 482 486 488 490 496 498 504 8 > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > : 9 > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > ; The complementary set of encodings will be denoted by DEL: DEL ¼fsAf0,...,511g9s=2DELg For each image of size pq, a SN P system will be provided. The input of the SN P system is a non-negative integer which represents the number of iterations in the skeletonizing process. Formally, given an pqimage, we associate to it the following SN P system is degree ðpqÞþ2: P ¼ðO, s 11 , s 12 ,..., s pq , s odd , s even ,syn,oddÞ, i.e., a SN P system with a neuron for each pixel in the image plus two extra neurons, s odd and s even . The input neuron is s odd .  O¼fagis the singleton alphabet.  s odd ¼ð512,a 513 a n =a 513 -a 512 Þ and s even ¼ð0,a 512 =a 512 -a 512 Þ. s ij ¼ðn ij ,R ij Þ,iAf1,...,pg,jAf1,...,qg, where n i,j ¼0ifi¼13i¼p3j¼13j¼q codði,jÞotherwise ( R 1k ¼R pk ¼R j1 ¼R jq ¼|for all kAf1,...,qgand jAf1,...,pg. For the remaining (i,j) R ij ¼fa b =a 511 -a 256 9b¼rþ5124rADELg[ fa b =a 512 - l 9b¼rþ5124rADELg  syn ¼ð S p1 i¼2 S q1 j¼2 syn ij Þ[syn odd [syn even [f/odd,even,1S, /even,odd,1Sg, where syn odd ¼/odd,ði,jÞ,1S9iAf2,...,p1g,jAf2,...,q1g, iþj odd () syn even ¼/even,ði,jÞ,1S9iAf2,...,p1g,jAf2,...,q1g, iþj even () and for all iAf2,...,q1g,jAf2,...,q1g, syn ij ¼ /ði,jÞ,ðiþ1,jÞ,1=2 0 S,/ði,jÞ,ði1,jÞ,1=2 4 S, /ði,jÞ,ðiþ1,jþ1Þ,1=2 1 S,/ði,jÞ,ði1,j1Þ,1=2 5 S, /ði,jÞ,ði,jþ1Þ,1=2 2 S,/ði,jÞ,ði,j1Þ,1=2 6 S, /ði,jÞ,ði1,jþ1Þ,1=2 3 S,/ði,jÞ,ðiþ1,j1Þ,1=2 7 Sg 8 > > > > > < > > > > > : 9 > > > > > = > > > > > ; The SN P system has one neuron s ij for each pixel of the image plus two extra neurons s odd and s even . The neurons corresponding to the border of the image has zero spikes in the initial configuration, the remaining neurons s ij corresponding to the pixels Fig. 3. (a) Enumeration of the pixels in a 3 3 neighborhood; (b) 3 3 neighborhood with encoding [0,0,0,0,1,1,1,1,1], or, shortly, 2 4 þ2 5 þ2 6 þ2 7 þ2 8 ¼496 and (c) scheme of the weights of the synapses. 6 Notice that we encode black pixels as 0 and white pixels as 1 for an easier implementation with SN P systems. In the previous section, we keep the opposite encoding in order to keep continuity with the literature. 7 Similar ideas are also used in [23]. of the image (called hereafter, the regular neurons) have codði,jÞ spikes in the initial configuration. The neurons s odd and s even have 512 and zero spikes, respectively. We will add to the 512 spikes in s odd as many spikes as indicated as input for starting the computation. The neuron s odd has only one rule a 513 a n =a 513 -a 512 which is applied if the neuron has at least 513 spikes. The neuron s even has also one rule a 512 =a 512 -a 512 . The regular neurons have two types of rules: Firing and forgetting ones, which will be applied if the number of spikes is exactly b¼rþ512 with rADEL or rADEL, respectively. With respect to the synapses, a regular neuron corresponding to a pixel Pis linked to the eight neurons corresponding to the eight neighbor pixels of P, with the weights 1=2 i where iAf0,...,7gfollows an anti-clockwise enumeration of the pixels starting in the east pixel, as shown in Fig. 3(c). The neurons s odd and s even are linked each other. The neuron s odd is also linked to all the regular neurons s ij with iþjodd and, analogously, s even is linked to all the regular neurons s ij with iþjeven. 4.1. How does it work? In order to understand how the SN P system works, first we observe that at the initial configuration, the set of regular neurons s ij encodes the image which will be skeletonized. We will show that, at any time, these neurons encode the successive images obtained in the iterative process of deleting black pixels according to the Guo and Hall algorithm. Let us remark that if the number of spikes in s ij is even, then the corresponding pixel is black; otherwise, if the number of spikes is odd, then the corresponding pixel is white. This is easily derived from the definition of codði,jÞ. Another observation to be considered is that the parity of the number of spikes in a neuron never changes, since the number of spikes received or removed is always an even amount, except by the application of the rule a b =a 511 -a 256 . As we will see in Lemma 2, the application of this rule is interpreted as the deletion of the corresponding black pixel in the Guo and Hall algorithm and it is applied once at most in each neuron. Before explaining the different steps of the process, let us consider a pixel (i,j) in the image and a black pixel adjacent to (i,j). We identify this black pixel to vAfP1,...,P8gaccording to the clockwise enumeration described above. Let us suppose that the black pixel in vbelongs to the selected subsection in the current step of the Guo and Hall algorithm and it satisfies the conditions to be deleted. Let us consider now the neurons s ij and s v corresponding to the pixels in (i,j) and v. As we will show in Lemma 1, the three conditions for v(it is black, it belongs to a selected subsection and it satisfies the Guo and Hall conditions to be deleted) indicate that the number of spikes in s v is b¼rþ512 with rADEL. In this case the rule a b =a 511 -a 256 is applied. As pointed out above, the application of this rule changes the parity of the number of spikes s r (from even, since the pixel is black, to odd) and this change is interpreted as a deletion of the pixel in v. We focus on the influence of the deletion of the black pixel in v (or, equivalently, the application of the rule a b =a 511 -a 256 in s v ) on the neuron s ij . Since the number of spikes in s ij at the initial configuration is codði,jÞ¼ P 8 i¼0 Hi 2 i and, in this configuration, the pixel vis black, according to the encoding, this means that Hv is zero, or, in other words, 2 v does not appear in the encoding of the environment as an addition of powers of 2. The deletion of the pixel in vchanges the environment of (i,j) and then, since the number of spikes in s ij represents such environment, the number of spikes must change. In particular, the change corresponds to turn Hv to 1, or, in other words, to add 2 v to the number of spikes in s ij . In order to check that this happens, it is suffices to see that the rule a b =a 511 -a 256 sends 256 ¼2 8 spikes from neuron s v to s ij , but these 2 8 spikes must be multiplied by the corresponding weight in f1=2 0 ,...,1=2 7 g, so only 2, 2 2 ,2 3 ,y,2 7 or 2 8 spikes arrive to s ij , depending on the value of v. A simple inspection shows that the number of spikes that arrives to s ij is exactly 2 v when the black pixel vis deleted. Bearing in mind these considerations, we show that for an input NAf1,...,513g, the computation steps of the SN P system correspond to the iterative process of the Guo and Hall algorithm where the first selected subsection corresponds to pixels with iþj odd. In such way, we will show the following statements:  Statement 1: The set of regular neurons is split into two subsections. One of the sections is composed of the neurons s ij such that iþjis even. Alternatively, the second sub-section corresponds to the neurons s ij such that iþjis odd. Both subsections are alternatively selected, starting with the odd subsection.  Statement 2: In each computation step, only neurons corresponding to the selected subsection are evaluated. The evaluation consists of determining if the neighborhood of the pixel associated to the neuron satisfies the conditions of the Guo and Hall algorithm to be deleted.  Statement 3: If an evaluated black pixel satisfies the Guo and Hall conditions to be deleted (see Section 3), then the number of spikes in the corresponding neuron changes from even to odd.  Statement 4: In each configuration, the number of spikes in the regular neurons is the codification of an image, according to the encoding described above. The first key point of the algorithm is that the image is split into two subsections which will be explored alternatively. One black pixel will be considered for its deletion only if it belongs to the subsection selected in the current step. We consider that a regular neuron s ij is selected at the step rif its number of spikes in the configuration C r is greater than or equal to 512. Otherwise, if its number of spikes is lower than 512, then the neuron is not selected. In the algorithm, the regular neurons with iþjodd and even are alternatively selected. The selection of subsections is performed by the neurons s odd and s even which send, alternatively, 512 spikes to the regular neurons s ij with iþjodd and even, respectively. Lemma 1. Let s ij be a regular neuron and NAf1,...,513gthe input of the SN P system. For rAf0,...,N1g  If iþjis odd,then the number of spikes in s ij is greater than or equal to 512 in the configuration C 2rþ1 and it is lower than 512 in the configuration C 2r .  If iþjis even,then the number of spikes in s ij is greater than or equal to belongs to512 in the configuration C 2rþ2 and it is lower than511 in the configuration C 2rþ1 . Proof. Let us observe that in the initial configuration, the number of spikes in a regular neuron s ij is codði,jÞo512; the number of spikes in s odd is 512þNand there is zero spikes in the neuron s even . From this initial configuration, the unique applicable rule is a 513 a n =a 513 -a 512 in the neuron s odd (since NZ1Þ. After applying this rule, in the configuration C 1 , the number of spikes in s odd is N1; the number of spikes in s even is 512; and the number of spikes in s ij is codði,jÞþ512 if iþjis odd and codði,jÞif iþjis even. Let us focus now on s odd and s even . The unique neuron that sends spikes to s odd is s even and, analogously, the unique neuron that sends spikes to s even is s odd . In the configuration C 1 , the spikes in s odd and s even are N1 and 512, respectively. Since Nr513, the rule in s odd cannot be applied in this configuration, but the rule in s even can be applied, so the spikes in s odd and s even in the configuration C 2 are 512þN1 and 0, which is similar to the situation in the initial configuration, so we have that for rAf1,...,Ng, the number of spikes in s odd and s even and in the configuration C 2r are 512þNrand 0. Notice that at the configuration C 2N , the number of spikes in s odd and s even are 512 and 0, respectively, and no more rules are applied in these neurons. According to the number of spikes in s odd and s even in the odd and even configurations, and taken into account their synapses, then we have that, for rAf0,...,N1g, at the configuration C 2rþ1 , the regular neurons with iþjodd has at least 512 spikes; and at C 2rþ2 , the regular neurons s ij with iþjeven has at least 512 spikes. In order to complete the proof, it is necessary to prove that for rAf0,...,N1g, at the configuration C 2rþ1 , the regular neurons with iþjeven has at most 511 spikes; and at C 2rþ2 , the regular neurons s ij with iþjodd has at most 511 spikes. Let us start by considering a regular neuron s ij with iþjodd at the configuration C 1 . As we show above, its number of spikes is b¼512þrwith r¼codði,jÞr511. Depending on rADEL or rADEL, one of the rules a b =a 511 -a 256 or a b =a 512 - l is applied.  Let us suppose that rADEL. In particular, this means that the corresponding pixel is black and the applied rule is a b =a 511 -a 256 . The number of spikes in the configuration C 2 is equal to the spikes in the configuration C 1 ð512þrÞ, minus the consumed ones 511 plus the contribution of other neurons. Since s odd does not send any spike in this step, the unique contribution to the number of spikes comes from other regular neurons. Each contribution is the addition of 2 i spikes to the spikes in s ij , but, bearing in mind that the pixel is black, then 2 0 does not appears in the decomposition of codði,jÞas sum of powers of 2, and it cannot be added as a contribution of other neuron, so rplus the contribution of other neurons is at most 510 and then, the number of spikes in s ij in C 2 is lower than 512.  If r=2DEL, then the rule a b =a 512 - l is applied. Since 512 spikes are consumed and rplus the contributions of the other regular neurons is at most 511, then the number of spikes in s ij is lower than 512, also in this case. This reasoning can be also applied to show that the number of spikes of the neurons s ij with iþjeven is lower than 512 in the configuration C 3 and in general we have that for rAf0,...,N1g, at the configuration C 2rþ1 , the regular neurons with iþjeven has at most 511 spikes; and at C 2rþ2 , the regular neurons s ij with iþj odd has at most 511 spikes. & Lemma 1 shows that the property to have at least 512 spikes changes alternatively from neurons s ij with iþjodd and even. From this result, it is easy to check the second statement, since the rules in the regular neurons can be only applied if the number of spikes is at least 512. This means that only in such cases the neuron is considered for evaluation. Evaluating a neuron consists on deciding if the rule a b =a 511 -a 256 or a b =a 512 - l is applied, but such decision depends on the set DEL which are the set of encodings of the neighborhood such that the central pixel must be deleted. The next key point of the algorithm is the deletion of pixels. By definition of codði,jÞ, a regular neuron s ij has an odd number of spikes if and only if it represents a white pixel. Analogously, a regular neuron s ij has an odd number of spikes if and only if it represents a white pixel. Bearing in mind this coding of black pixels, deleting a black pixel in a computation step consists on removing an odd amount of spikes from a neuron with an even amount of spikes. Lemma 2. Let us consider a black pixel and a step of the Guo and Hall algorithm,such that the pixel belongs to the selected subsection and it satisfies the conditions of the algorithm to be deleted. Then,in the corresponding step of the SN P system computation,the corresponding regular neuron s ij will pass from an odd amount of spikes to an even amount. Proof. According to the previous construction, a black pixel which belongs to the selected subsection in the Guo and Hall algorithm and verifies the conditions to be deleted has associated a regular neuron with 512þrspikes, rADEL. In this case, the rule a b =a 511 -a 256 is applied. Bearing in mind that ris even, the contributions of other neurons is even and an odd number of spikes is consumed, in the next configuration, the number spikes in the neuron is odd. & Since all the rADEL are even and the contribution of other neurons is always even, then the rule a b =a 511 -a 256 with b¼512þris applied at most once in each neuron. This means that if a pixel is deleted (changed from black to white) it never becomes black to white, and the iterative process of thinning the image is also carried out in the SN P system. Finally, to sum up these statements. We claim the following result. Theorem. The set of regular neurons of the SN P system encodes in each configuration the successive images obtained in the iterative process of thinning of the Guo and Hall algorithm, by taking as black the pixels with an even amount of spikes and white the neurons with an odd amount of spikes. 5. Experimental simulation Simulation of different variants of P systems has been widely studied in the last years. Nonetheless, up to now, no implementation in vivo or in vitro has been carried out and the natural way to explore the behavior of designed P systems is to simulate it in conventional computers. A short description of some of these simulators can be found in [26,27]. Currently, a big effort is being developed in the P-lingua project [28,32], by combining an efficient simulation engine with an ad-hoc programming language. In this paper, a software tool based on the design of the SN P system has been implemented by using Compute Unified Device Architecture (CUDA TM )[29]. The experiments have been performed on a computer with a CPU AMD Athlon II 4 645, which allows to work with four cores of 64 bits to 3.1 GHz. The computer has four blocks of 512 KB of L2 cache memory and 4 GB DDR3 to 1600 MHz of main memory. The used graphical card (GPU) is an NVIDIA Geforce GT240 composed of 12 Stream Processors with a total of 96 cores to 1340 MHz. It has 1 GB DDR3 main memory in a 128 bits bus to 700 MHz. So, the transfer rate obtained is by 54.4 Gbps. The used Constant Memory is 64 KB and the Shared Memory is 16 KB. Its Compute Capability level is 1.2 (from 1.0 to 2.1). In the following, we comment the workflow of our application, where every square in Fig. 4 represents a parallel step using CUDA. Step1: First of all, the Initial State is created in parallel, by encoding the pixels from the input image. The matrix is created and stored in the device memory and it is represented as the box STATES MATRIX. Step2: The initial STATES MATRIX is consumed via the execution of the SN P system algorithm. In parallel, all the neurons are checked. If a rule can be applied, then, according to the SN P system model, it is triggered. In order to spread in parallel the number of spikes, four new matrices with the new spikes are created in device memory, called SPIKES MATRIX. The synchronization of the emitted spikes is made via these four matrices. Step3: Next, the state of the SN P system is updated in parallel from the current state and the four SPIKES MATRIX (which are consumed). The next state is stored in STATES MATRIX. Step4: The neurons are checked again. If the number of the spikes in a neuron has changed (i.e., if a rule has been triggered), then the implementation comes back to Step 2. Otherwise, it goes to Step 5. Step5: The algorithm finishes by creating the output image from the halting configuration of the P system. The flow of pieces of information between the HOST and the DEVICE is the following:  The process starts in HOST. The needed information for the parallel calculus is copied into DEVICE. Then, the process that triggers the threads from DEVICE is called.  The threads are executed in parallel from DEVICE. The threads perform the calculus and the control comes back to the HOST.  The HOST copies the results from the DEVICE in order to get the information obtained in parallel. The execution of the threads in parallel depends on the available resources in the GPU: CORES and memory. The execution threads are grouped in 2-dimensional blocks (Block Threads). Each block has 16 16 threads. Such Block Threads are grouped in a 2-dimensional grid. For an image of pqpixels, at least a grid with ðp=16Þðq=16ÞBlock Threads of 16 16 threads are needed. The implementation deals with Nblocks of threads for the complete image in our GPU of 96 cores. Let us notice if the height and width of the image are not multiples of 16, then more threads Fig. 4. Workflow of the implementation. Fig. 5. Scheme of the threads. Fig. 6. Example of image with traffic signals. than pixels are needed. In this case, in each block some threads are not used. Fig. 5 shows the sequential steps of the process. In each step the 96 GPU cores are working in parallel. If the number of blocks is not multiple of 96, in the last step there are unused blocks. 5.1. Examples A first example is shown in Fig. 2. Skeletonizing hand-written texts is one of the challenges of skeletonizing, since the skeleton keeps the topological structure and meaning of the original and the text can be easily stored. 8 Next, we provide several examples of a realistic recognizing problem with applications in the automotive industry. In Fig. 6, we can see a photograph taken in a road. It has been binarized by using a threshold method by using a threshold 100 on a gray scale 0,...,255. We can see that the skeletonized images in Fig. 7 keep Fig. 7. (a) The binarization of the image from Fig. 6, (b) inverse binarization of the image, (c) its skeletonizing and (d) the skeletonizing of the inverse thresholding. Fig. 8. Original images and their skeletons. 8 For a description of the use of skeletonization in image processing, see, for example, [30]. the information of the traffic signals and they can be used in a further pattern recognition problem. In Fig. 8, two more examples off skeletonizing real images are shown. We finish this section by showing the results of some experiments performed with our implementation. We have taken 36 totally black images of nnpixels, 9 from n¼512 to n¼5120 with a regular increment of 512 pixels of side. By taking this initial number of pixel and this increment we are sure that the number of threads for each block is a divisor of the size of the image. Fig. 9 shows the time in milliseconds of our software tool inspired in the designed SN P system for implementing the Guo and Hall algorithm for 1, 30, 60 and 90 steps in the skeletonizing process. Fig. 10 shows a comparison between the implementation of the Guo and Hall algorithm via Cellular Automata (CA) presented in [31] and the implementation via SN P system shown in this paper. The time of both implementations have been obtained in the same machine from totally black images. The obtained time shows that the current implementation of the SN P system on a GPU takes more time that the CA implementation. One of the reasons for this drawback is that the use of an alphabet with only one object, the spike a, does not fit in the GPU architecture. In this way, further research is needed. In the one hand, new features of the theoretical model can be explored in order to adapt the SN P system to the new hardware architecture. On the other hand, the designed SN P system can be implemented in different hardware architectures by looking for improvements in the efficiency. Fig. 11 shows the speed-up of the implementation of the SN P system in the GPU vs. the sequential one, i.e., the ratio of the time taken by both implementations. The study has been made on the same set of black images used in the comparison shown in Fig. 9. As expected, ratio increases when the size of the image or the number of steps increases. 6. Conclusions The development of new bio-inspired parallel techniques provides a chance for revisiting classical algorithms. In this paper, we have consider a classical algorithm for skeletonizing images, but many other algorithms can be considered. In particular, the bio-inspired computing techniques have features as the encapsulation of the information, a simple representation of the knowledge and parallelism, which are appropriate with dealing with digital images. The main contribution of this paper is to show that classical parallel algorithms, as the Guo and Hall algorithm, can be implemented in a new computational paradigm inspired in the Fig. 9. Experimental time obtained for our SN P system implementation of the Guo and Hall algorithm 36 totally black images of nnpixels, from n¼512 to n¼5120 with a regular increment of 512 pixels of side. Fig. 10. Experimental comparison between the implementation of the Guo and Hall algorithm via Cellular Automata [31] and the implementation via SN P systems presented in this paper. 9 Theoretically, this is the worst case, since the time inverted by the algorithm depends on the size of the biggest black connected component of the original image.