Membrane parallelism for discrete Morse theory applied to digital images

Membrane parallelism for discrete Morse theory applied to digital images Raúl Reina-Molina · Daniel Díaz-Pernil · Pedro Real · Ainhoa Berciano Abstract In this paper, we propose a bio-inspired membrane computational framework for constructing discrete Morse complexes for binary digital images. Our approach is based on the discrete Morse theory and we work with cubical complexes. As example, a parallel algorithm for computing homology groups of binary 3D digital images is designed. Keywords Discrete Morse theory · Digital Imagery · Homology groups · Membrane Computing 1 Introduction Nowadays, to determine in a fast and accurate way topology-related information of technologically relevant mathematical structures has become a very important question. Within the context of Digital Imagery, computational topology has acquired a significant role in fields of applications like computer vision, digital image processing or medical imaging. Working with cubical cell complex representations of digital R. Reina-Molina (B)·D. Díaz-Pernil ·P. Real CATAM Research Group, Department of Applied Mathematics I, University of Seville, Seville, Spain e-mail: [email protected] D. Díaz-Pernil e-mail: [email protected] P. Real e-mail: [email protected] A. Berciano Department of Didactic of Mathematics and Experimental Sciences, University of the Basque Country, Bilbao, Spain e-mail: [email protected] 50 R. Reina-Molina et al. images, we rewrite here a non-negligent part of the computational machinery underlying in discrete Morse theory [13], using bio-inspired models of Membrane Computingarea[25]. Discrete Morse theory (DMT) deals with the combinatorial homotopy analysis of cell complexes using discrete Morse functions or, equivalently, discrete vector fields. Membrane Computing deals with distributed and (massively) parallel computing models inspired from cellular membranes and locally processing multisets of symbol objects. This paper is a foundational effort to properly combine the computational power of bio-inspired computing with the detailed topological analysis of cubical complexes that DMT provides us. In this way, we implement a membrane computing framework for constructing “simplified” digital images with regard to some special topologically consistent properties mainly related with homology groups. As a result of the inherent parallelism involved in the definition of digital images and due to the current improvement in parallel computational architectures, several studies are devoted to parallel implementation of sequential algorithms in order to take advantage of parallel hardware. The development of parallel algorithms is not a trivial question and involves an adequate management of task’s concurrency and a performance, memory efficiency and scalability trade-off. Natural Computing has been presented as a productive research field that can develop no ordinary computational models inspired by Nature that can solve several difficult problems in an elegant and simple way. Here, to maintain a strong link between this computational framework and algebraic topology tools is a priority. In this way, we obtain a topological description of the digital objects which allows us to understand, efficiently compute and handle advanced topological information of kD digital objects and images [2,8]. For this purpose, we present a new methodology to develop algorithms in topological analysis of digital images. Within this bioinspired setting, our approach mainly focus on performance more than on memory efficiency or scalability. This paper is the first one using membrane computing for DMT techniques. Various DMT-based approaches to compute discrete functions defined on cells of a cell complex have been already used. Reininghaus et al. [29,30] apply DMT methods to analyze vector fields. Bauer et al. [1] computed simplified two-dimensional scalar functions while ensuring that the input function is modified by no more than a threshold d and all surviving critical point pairs have persistence greater than 2d. Cazals et al. [7] and Lewiner et al. [17] successfully employed Forman’s discrete Morse theory to compute Morse complexes of piecewise-linear functions and show applications to segmentation, visualization, and mesh compression. Gyulassy et al. [15] also used a DMT-based formulation to develop an efficient algorithm for computing Morse (more concretely, Morse-Smale) complexes of large 3D data. Robins et al. [30] proposed an algorithm to compute the Morse complex of 2D and 3D gray-scale digital images modeled as discrete functions on cubical complexes. The algorithm computes the Morse complex with provable guarantees on its correctness with respect to the critical cells. In [20], an algorithm for computing homology using a recursive tree-based technique for generating Morse complexes is established. Peterka et al. [27] introduced a set of building blocks for implementing parallel algorithms, which leverage high performance computing clusters. In particular, they discuss a parallel implementation of the discrete Morse theory based on the algorithm proposed by Gyulassy et al. [15]. Membrane parallelism for DMT in digital images 51 Günther et al. [14] design a memory efficient algorithm to compute the Morse–Smale complex for 3Ddata and use the complex to compute persistent homology groups. The discrete gradient field is computed using a parallel variant of the method proposed by Robins et al. followed by an efficient computation of the boundary map representing the Morse complex. Finally, the parallel algorithm for computing the discrete gradient field given in [22] is based on a novel description of the discrete Morse function followed by a two-step algorithm computing the cells of the Morse–Smale complex. The algorithms are implemented using a hybrid multi-core implementation. Thepapers[6,9,10],dealwiththe computationofsomecombinatorial andalgebraic aspects of the topology of the image. This paper is structured as follows. First, some aspects related with cubical homology and DMT are reviewed. Then, basic concepts and techniques from Membrane computing are given. Next, a theoretical framework for designing algorithms in this “hibrid” context and some examples of its use are presented. Finally, some conclusions are discussed,where the most important fact is the theoretical viability of framework shown before, that allows us to implement algorithms computing the homology groups of 3D images using parallel computing tools. 2 Discrete Morse theory In this section we introduce the topological background required by our framework. First of all, we review the required combinatorial structure for the topological spaces that are used. We mainly follow the process introduced by Kaczy´nski et al. [16] with some minor changes. Specifically, Kaczy´nski et al. uses closed cubes as main combinatorial objects while we use open cubes. Acubical cell σis a finite product of intervals: σ=I1×···×Ik⊂Rk where Ijis an interval with integer extremes (mj,mj+1)and length one or the singleton {mj}, denoted as (mj), for each j∈{1,...,k}. The interval Ijis referred as the jth component of σand denoted by Ij(σ). The set of all cubical cells in Rkis denoted by Kk. The set of all cubical cells is K= ∞  k=1 Kk(1) We usually require the cubical cells to be bounded. Hence, we define Kk ∗,n={σ∈Kk:0≤inf Ip(σ), sup Ip(σ) < n,1≤p≤k}(2) where kand nare nonnegative integer numbers. Given a cubical cell σin Rk, its embedding number k is denoted by emb σ.The dimension of σis defined as the number of unitary intervals in its expression as product of intervals and is denoted by dim σ. The set of all elementary cubes with dimension pis denoted by Kp. The set of all elementary cubes in Rkwith dimension pis denoted by Kk p. Whenever the dimension of a cubical cell require to be made 52 R. Reina-Molina et al. explicit, it is denoted as a superscript between parenthesis. Therefore, dim σ(p)=p. We also explicitly indicate the dimension of a cell σ(p)writing that σis a p-cell. Theclosure ofacubicalcellcanbedecomposedinto theunionoflower-dimensional cubical cells. If δand σare two cubical cells in Rkof any dimension and δ⊂σ, then δis a face of σand is denoted as δ≤σ.Ifδis a face of σand δ= σ, then δis a proper face of σ, denoted as δ<σ.Ifδis a face of σand dim δ=dim σ−1 then δ is a primary face of σ, denoted by δ∈∂σ. Acubical complex is a collection Kof cubical cells with the same embedding number and such that, for every cubical cell σ∈K, all of its primary faces are in the complex. We denote the set of p-cells in Kas Kp. As usual, we define the chain complex C(K)as the graded Z-module {Cp(K)}p∈Z where Cp(K)is the free abelian Z-module generated by the cubical p-cells in K.The boundary map ∂in the chain complex is defined in [16]. Concretely, for any cubical cell σ,∂σ is an alternating sum over the faces of σand such that ∂◦∂=0. As usual, let Zp=ker ∂pthe subgroup of cycles and Bp=Im ∂p+1the subgroup of boundaries. The pth cubical homology group of Kis Hp(K)=ZpK Bp(K) The cubical homology of Kis the collection H∗(K)={Hp(K)}p∈Z. Recall that, for p<0orp>dim K,isHp(K)=0. 2.1 Discrete Morse theory So far, we have recalled all the required notions from cubical complexes. Now, we introduce the needed elements from DMT. A real-valued smooth map defined over a compact d-manifold is a Morse function if all its critical points have non-singular Hessian matrix and no two critical points have the same function value [15]. Morse functions allow to endow the manifold with a cellular structure. A CW complex is, roughly speaking, a topological space endowed with a decomposition into smaller pieces called cells, being homeomorphic to Euclidean spheres. A formal definition can be found in [13]. Discrete Morse theory introduced by Forman [13] adapts Morse Theory to CW complexes instead of smooth manifolds. Although DMT has been developed for finite regular CW complexes, we restrict our research to finite cubical complexes, which are a particular type of CW-complexes. Let Kbe a cubical complex and f:K→Ra function that assigns scalar values to every cell of K.fis a discrete Morse function if, for every cubical cell σ(p)∈K the following conditions hold: –#{τ(p+1)∈K|τ>σ∧f(τ) ≤f(σ )}≤1 –#{μ(p−1)∈K|μ<σ∧f(τ) ≥f(σ )}≤1 A discrete Morse function fcan be thought as a discrete function that is increasing with respect to cell dimension except for, at most, two cells for each cell. Concretely, Membrane parallelism for DMT in digital images 53 for each cell σthere is at most one face μ(μ∈∂σ) with f(μ) ≥f(σ) and there is at most one facet τ(σ∈∂τ) 0 with f(τ) ≤f(σ). A cubical cell σ(p)is critical if one of the following conditions hold: –#{τ(p+1)∈K|τ>σ∧f(τ) ≤f(σ )}=0 –#{μ(p−1)∈K|μ<σ∧f(τ) ≥f(σ )}=0 Given a cubical cell σits barycenter is the point b(σ) =inf I1(σ ) +sup I1(σ) 2,...,inf Ik(σ) +sup Ik(σ ) 2. As the closure of a cubical cell is convex, its barycenter is an interior point. We define a discrete vector as a pair of cells {σ(p)<τ (p+1)}. Vectors are represented as arrows from the barycenter of the lower dimension cell to barycenter of the higher dimension cell. A discrete vector field V on Kis a collection of vectors in K such that each cubical cell in the vector belongs to, at most, one pair of V. Formally, a discrete vector field is a map V:K→K∪{0}such that: 1. for each σ∈K,ifV(σ ) = 0 then dim V(σ ) =dim σ+1 2. for each σ∈Kp, either V(σ ) =0orσ∈∂V(σ) 3. if σ∈Im(V)then V(σ ) =0 4. for each σ∈Kp #{μ(p−1)∈K|V(μ) =σ}≤1 Given a discrete vector field Von K,aV-path of dimension p,γ, is a sequence of cubical p-cells σ0,σ 1,σ 2,...,σ rsuch that 1. If V(σi)=0, then σi+1=σi. 2. If V(σi)= 0, then σi+1<V(σi)with σi+1= σi. The set of V-paths is denoted as Γ(V). AV-path γ=σ0,σ 1,...,σ ris called closed is σr=σ0and is called nonstationary if σ1= σ0. Given two V-paths γ1=σ0,...,σ rand γ2=σr,...,σ r+s,theconcatenation γ2◦γ1is defined by γ2◦γ1=σ0,...,σ r,σ r+1,...,σ r+s(3) Two V-paths γ1and γ2aresaidtobeconnectible if there is a V-path γ3such that γ3=γ2◦γ1 AV-path γ1is said to be extendable by another V-path γ2whenever there is a V-path γsuch that γ2=γ◦γ1, and is denoted by γ1γ2. It can be easily verified that (Γ (V), )is a partially ordered set. Whenever a V-path γ1is extendable by a V-path γ2so that γ2=γ◦γ1as well as γis non-stationary, it is denoted as γ1γ2. 54 R. Reina-Molina et al. Given a V-path γ=σ0,σ 1,...,σ r,themultiplicity of γis defined as m(γ ) = r−1  i=0 V(σi)=0 σi,∂V(σi)σi+1,∂V(σi)(4) Given two connectible paths γ1and γ2, then m(γ2◦γ1)=m(γ2)·m(γ2)(5) AV-path γis said to be maximal if there is no V-path γsuch that γγ. A maximal path is called critical if its last cell is critical. Notice that, if γ=σ0,...,σ ris a maximal critical path, then there is a critical cell μsuch that σ0,∂μ = 0. Given a discrete Morse function, a special discrete vector field called discrete gradient vector field is defined so that f(V(σ)) ≤f(σ ).Asshownin[13], a cubical complex can be transformed into another homotopically equivalent following a series of simplicial collapses,1where each of them collapses both cubical cells in each vector in the corresponding discrete gradient vector field. We recall here one of the main results of DMT. Theorem 1 ([13], 9.3) A discrete vector field V is the gradient vector field of some discrete Morse function if and only if there are no non-stationary closed V-paths. A vector field Vis extended to a graded group homomorphism V:Cp(K)→ Cp+1(K)such that V(σ(p))=±τ(p+1)if {σ<τ}is a vector in V 0 otherwise (6) where the sign is chosen so that σ, ∂V(σ )=−1. The reduced (discrete time) flow map, denoted as ˜ φis defined by ˜ φ=id +∂◦V(7) The reduced flow map associates a (p+1)-chain ˜ φ(σ) to any p-cell σ, such that ∂V(σ) =˜ φ(σ)−σ. Hence, σ, ˜ φ(σ)=0. The map ˜ φestablishes a way of following the flow determined by the vector field. Concretely, if there is a V-path of length r from a p-cell σto another p-cell σ, then ˜ φr(σ), σ  = 0. This fact can be proved by induction in the length of the V-path. Definition 1 Let Kbe a cubical complex and fbe a Morse discrete function on K. Let Cp(K)denote the p-chains on Kand Mp⊆Cp(K)the span of the critical p-cells. The graded group M={Mp}p∈Zis called the space of Morse chains. 1 In [13] Forman works with simplicial complexes, however the mathematical scaffolding provided by DMT can be settled with no change to finite cubical complexes. Membrane parallelism for DMT in digital images 55 In [13] Forman makes use of the set of V-paths to build the boundary map of the Morse complex associated to a given complex and an acyclic vector field. It is shown in the following result. Theorem 2 ([13],7.1) Thereareboundarymaps ˜ ∂p:Mp→Mp−1,for each p ∈Z, so that ˜ ∂p−1◦˜ ∂p=0 and such that the differential complex 0˜ ∂k+1 −−→ Mk ˜ ∂k −→ Mk−1 ˜ ∂k−1 −−→··· ˜ ∂2 −→ M1 ˜ ∂1 −→ M0 ˜ ∂0 −→ 0 calculates de homology of K. That is, if we define Hp(M,˜ ∂) =ker ˜ ∂p Im ˜ ∂p+1 then, for each p Hp(M,˜ ∂) ≡Hp(K). In [13] an explicit expression for the boundary operator in the Morse complex (M,˜ ∂) is also given in terms of V-paths. Theorem 3 ([13], 7.3) For any critical cubical (p+1)-cell τ,let˜ ∂τ be defined as ˜ ∂τ = σ∈Mp cσ,τ σ where cσ,τ = γ∈Γ(τ,σ) m(γ ) and Γ(τ,σ)is the set of maximal critical V-paths (for a given gradient vector field V) which go from a maximal face of τto σ. The complex (M,˜ ∂) and K have isomorphic homology groups. 2.2 Examples In Fig. 1a 2D cubical complex Kis shown (Fig. 1a). In this example, cubical cells are referred as ordered tuples of vertices. A vector field Vis shown in Fig. 1b. As the previous vector field is acyclic, it is in fact the gradient vector field for some discrete Morse function f, whose critical cells are showed in Fig. 1c. 56 R. Reina-Molina et al. (a) (b) (c) (d) (e) Fig. 1 a Example of cubical complex, bgradient vector field, ccritical cells where two kinds of critical cells are shown: ordinary critical cells representing a hole at any dimension and other critical cells showing a “step” configuration which will be canceled later (see [13], section 11, for details), dcritical complex and efinal homological complex with null differential The reduced flow map ˜ φmaps a cell σ(p)into a p-chain c, such that all the cells in care those in ∂V(σ) except σ. For example, ˜ φapplied to the 1-cell spanned by P4 and P5, denoted by <P4,P5>(see Fig. 1a) returns φ ˜ φ ˜(P4, P5) =P4, P9+P9, P10−P5, P10 Informally speaking, the application of φ˜ means the annihilation of some cells involved in the vector field, as well as the “extension” of critical cells in the sense of creating critical “super-cells” in terms of chains. For example, the application of to P4, P5 not only removes P4, P5 and P4, P9, P10, P5, but it extends the cell P0, P4, P5, P1 to the chain P0, P4, P5, P1+P4, P9, P10, P5. This previous example provides us the key idea in the development of our framework. The homology of a 3D cubical complex can be calculated in only two stages. First the complex is reduced by the construction of a gradient vector field and the calculation of the associated Morse complex (see Fig. 1d). Second, the process is repeated in the Morse complex and, finally, the result is a chain complex with null differential whose pth homology group is easily calculable and coincides with the pth homology group of the cubical complex (see Fig. 1e). 2.3 Encoding images as cubical complexes and cubes as tuples In this section, we work with a cubical cell version of a digital image [16]. A common problem appearing in the development and implementation of algorithms with cubical Membrane parallelism for DMT in digital images 57 complexes is the requirement of a huge amount of memory resources. Below, an efficient encoding of cubical cells is described. Given two non-negative integers kand n≥2, the following auxiliary sets are defined Tk n:= {0,1,2,...,n−1}k,(8) Tk n:= Tk 2n−1,(9) In,k:= Tnk={0,1,2,...,nk−1},(10) In,k:= T(2n−1)k={0,1,2,...,(2n−1)k−1}.(11) The set of points2for the source images is the set Tk n⊂Zkequipped with a cubic neighborhood function, described as follows: two points i=(i1,...,ik)and j=(j1,..., jk)are said (2k)-adjacent if k p=1|ip−jp|=1. The neighborhood function is given by N(i1,...,ik)=⎧ ⎨ ⎩ (j1,..., jk)∈Tk n: k  p=1 |ip−jp|=1⎫ ⎬ ⎭ This function restricted to k=2 defines 4-adjacency and to k=3 defines 6-adjacency. Let I:Tk n→{0,1}be a k-D binary image of size nk, where the set of points in the object (or black points) is given by I−1(1).LetK=K(I)be the cubical cell complex associated to I.InK, the 0-cells are points in the object, the 1-cells represent pairs of (2k)-adjacent points, the 2-cells unit squares where its edges are pairs of (2k)-adjacent points, and so on. In general, each p-cell is a p-dimensional unit hypercube whose edges are determined by pairs of (2k)-adjacent points. A cubical cell is encoded as a tuple through the function T:Kk→Zkdefined as follows T(σ) =(sup I1(σ ) +inf I1(σ),...,sup Ik(σ ) +inf Ik(σ)) (12) The main property of Tis given in the following lemma. Lemma 1 The function Tis a bijection and T(Kk ∗,n)=Tk n. Proof Firstatall,letusprovethatTis injective. Let σand μbe two distinct cells. Then, let j0∈{1,2,...k}be the first index such that Ij0(σ ) = Ij0(μ). The following cases may come up: –Ij0(σ) =(a,a+1)and Ij0(μ) =(b,b+1)with a= b. Then inf Ij0(σ) =a= b=inf Ij0(μ) 2The reader is supposed to be familiar with concepts of Image Algebra. For a detailed text, see [31]. 64 R. Reina-Molina et al. Theorem 4 Given a k-D binary image I, the corresponding P system ΠCub builds the cubical complex for Iin two computation steps. Proof (Remark)LetI:Tk n→{0,1}be a binary k-D image. Let Kbe the cubical complex built from I. In the initial configuration, the first membrane consists on objects Bifor i∈I−1(1). In this situation, only some rules in R1can be selected, because the other rules require the presence of objects σi.TherulesinR1that can be selected are those for l∈In,ksuch that all the points in gen(σ (l)) are in I−1(1). Therefore, the first computation step builds the objects σlthat represents cubical cells in the complex and, also, introduce the mark for critical cells (Cl). In the second configuration, no rule in R1can be selected since the proper objects representing each cubical cell act as inhibitors for this rules. However, the rules in R2 and R3are susceptible to be selected. In fact, the selected rules in R2introduce into the first membrane the objects that symbolize the boundary map (δ± ij). The selected rules in R3introduce the objects that mark the (positive) direction of the boundary barycentric vector vσμ foreverypairofcellsσ, μ ∈Ksuch that σ∈∂μ and vσμ ·ud=1 2. The application of rules in R2and R3prevent the selection of any other rule in R2or R3any time, because the inserted objects act as inhibitors. Hence, the next computation step (the third) makes the system reach the halting condition, inasmuch as no rule can be selected.  Remark 1 If σis in the boundary of μ, then both cells share all its intervals except one, degenerated in σand non-degenerated in μ. We suppose here that this interval is the d-th. Then, all the coordinates of the vector vσμ are 0, except the d-th one, whose value is 1 2and, hence, vσμ =1 2ud. Definition 5 Let Kbe a cubical complex and ∂the associated boundary map. Let V be a gradient vector field in K.ThePsystemΠFlow is defined as ΠFlow =(Γ, Σ, E,w 1,RFlow,Pri,iin,iout) where – The alphabet of objects is given by Γ=Γ0∪{Eij :i,j∈In,k,i= j}. – The input alphabet is given by Σ={σi:i∈In,k(T(K))}∪{δ± ij :i,j∈In,k(T(K)) ∧σ(i), ∂σ ( j)=±1} ∪{V± ij ,Vi,Vj:i,j∈In,k(T(K)) ∧V(σ(i)) =±σ(j)} ∪{Ci:i∈In,k(T(K)) ∧V(σ (i)) =0∧V−1(σ(i)) =∅} ∪{Eij :i,j∈In,k(T(K)) ∧σ(j), ∂V(σ(i))=0∧i= j} – The environment alphabet is given by E=Γ0. – The set of rules Ris given by –R1≡{σiσjσlδs1 il δs2 jlV−s1 il }¬{φiφj}|1,λ/φ−s1s2 ij φiφj,0for i,j,l∈In,k, i= j= l= iand s1,s2∈{+,−} –R2≡{σiCi}¬{φ+ ii }|1,λ/φ+ ii ,0for i∈In,k Membrane parallelism for DMT in digital images 65 –R3≡{φs1 i1i2Vs2 i2j1φs3 j1j2}|1,Ei2j1/φ−s1s2s3 ij ,0for i,j,l∈In,k,s1,s2,s3∈ {+,−} and i= j= l= i. –R4≡{δs1 li CiCj}¬{˜ δs1s2 ji }|1,φs2 lj /˜ δs1s2 ji ˜ Ci˜ Cj˜σi˜σj,0for i,j,l∈In,kwith i= j= l= i. –R5≡1,˜ δ+ ij ˜ δ− ij/λ, 0for i,j∈In,k,i= j. –Pri =∅ –iin =iout =1. Theorem 5 Let K ⊂Kk ∗,nbe a cubical complex and V a gradient vector field in K. The module ΠFlow computes the Morse complex for K associated to V in, at most, 4+log2νcomputation steps, where ν≤nkis the length of the greatest V-path on K. Proof (Remark) The initial configuration for ΠFlow consists in all the elements that codify the cubical complex K(objects σifor cubical cells and Cimarking critical cells), its boundary map (objects δ± ij) and the gradient vector field (objects V± ij ). There also are objects Eij depicting that σiand σjare distinct maximal faces in V(σi).In this situation, only rules in R1and R2can be selected, as the other rules require the existenceofobjectsφ± ij introducedinthefirstmembranebytheserules.Theapplication of rules R1creates one object φ± ij for each pair of cubical cells σi,σj∈Ksuch that ˜ φ(σi), σj=±1. Recall that ˜ φ:C→Cgiven by ˜ φ=id +∂◦Vis the reduced (discrete time) flow of the gradient vector field [13]. Theapplicationofrules R2previouslyselectedintroduceobjectsφ+ ii foreachcritical cubical cell σi, representing that, for each critical cell μis ˜ φ(μ) =μ. In the current situation, the second configuration, only rules in R3,R4or R5are available to be selected. Recall that objects φ± ij encodes the presence of a V-path from σito σjwith multiplicity equals to ±1. Rules in R3concatenate two V-paths γand γif there is a cell σ such that the ending cell of γ, named σ(i), the starting cell of γlay on its boundary and V(σ(i)) =σ. Hence, rules in R3are applied until a maximal path is constructed. This process needs, at most, log2νsteps, where νis the length of the longest maximal V-path. Rules in R4get some objects from the environment, encoding the cells of the Morse complex. Indeed, for each critical path from σito σjwith multiplicity ±1, there is an object φ± ij, hence for each critical cubical cell σlwith σi,∂σ l = 0 there will be as many objects φ± ij as the number of critical V-paths from σito σjwith multiplicity ±1, respectively. Finally, rules in R5cancel pairs of objects ˜ δ+ ij,˜ δ− ij. The halting configuration of ΠFlow has as many objects ˜ δ± ij as the incidence of σiin the boundary of σj. Notice that the entire process needs, at most, 4 +log2νsteps, where νis the length of the longest critical V-path.  In the following section, we show an example of the use of our framework by designing a dP scheme for solving one Computational Algebraic Topology related problem. 66 R. Reina-Molina et al. 5 dP systems for the calculation of the homology groups for a 3D binary image We deal here with a membrane algorithm for computing an optimal gradient vector field for a 3D digital object I. Working with a cubical “cellularization” K(I)of I,we show that the critical cells of dimension p(p=0,1,2,3) of this particular gradient vector field coincide with the pth Betti number of K(I). Discrete Morse theory establishes a tool to simplify a cell complex6while the homology groups are kept. This can be used to calculate the homology groups of appropriate cell complexes. Concretely, the technique explained below require the complex to be torsion free. An acyclic vector field over a cubical complex Kis built, removing as many cubical cells as possible. In case that the resulting Morse complex has non-null differential, another acyclic vector field is constructed to accomplish the so called critical cell cancellation (see section 11 in [13]). In this way, a sequence of Morse reductions is built until a complex with null-differential is “reached”. To ensure the finiteness of this process, the cubical complex has to be torsion free. Figure1presentsanexample of the calculation of homologygroupsfora2Dcubical complex.In this examplethe strategypresented aboveis shown.First, an acyclicvector field is built such that almost every cubical p-cell σis paired to one of its co-faces μ(cubical (p+1)-cell with σ∈∂μ). The algorithm used to determine the acyclic vector field returns two kind of critical cubical cells: on one hand, cubical p-cells that represent a p-dimensional hole and, on the other hand, 1-cells that represents an “step” configuration (as P11,P12or P15,P16in Fig. 1). Therefore, another acyclic vector field is required tocancel these spurious critical cells. Notice that the first Morse complex is not a cubical complex but a CW complex. However, it is usually much more simpler than the original cubical complex. For binary 3D images, it suffices to build up to two acyclic vector fields to calculate an homotopically equivalent cell complex with null differential, which makes homology group calculation trivial. In the paragraphs below, two P modules are defined to perform this task: ΠVector and ΠVector2. Definition 6 Let I:T3 n→{0,1}be a k-D binary image. The module ΠVector is defined as ΠVector ={Γ,Σ,E,w 1,R,Pri,iin,iout} where –Γ=Γ0 –Σ={σi:i∈In,k(T(K(I)))} –E=Γ0 –w1=∅ –Ris the following set of rules: –R1≡¬{ViVj}|1,σ iσjδ± ijdipU+ ij,dCiCj/σiσjδ± ijU+ ij,dV∓ ij ViVj,0for i,j∈ In,3,i= j,0≤p≤3 and 1 ≤d≤3 6CW complex in [13] and cubical complex in this paper. Membrane parallelism for DMT in digital images 67 These rules take a critical cubical p-cell (σi), a critical cubical (p+1)-cell (σj) with σi∈∂σjand such that the barycentric vector vσiσjis 1 2ud. Then, objects V± ij are brought from the environment (this represents that V(σi)=±σj)as long as objects Viand Vj(used to prevent that a cubical cell is used in two different vectors in the vector field). – The priorities are given by Pri ={R1[i,j,p,d]>R1[i,j,p,d]:i,j,j∈In,3∧0≤p≤3∧1 ≤d<d≤3}∪{R1[i,j,p,d]>R1[i,j,p,d]:i,i,j∈In,3∧0 ≤p≤3∧1≤d<d≤3}∪{R1[i,j,p,d]>R1[j,j,p+1,d]: i,j,j∈In,3∧0≤p≤2∧1≤d≤3} where R1[i,j,p,d]represents the corresponding instance of rule R1for the given values of the parameters. The priorities in the first set are used to assign to a cubical cell σia cubical cell σj with the corresponding barycentric vector parallel to the vector udwith the lowest das possible. The priorities of the second kind are used to choose between two cubical cells σi and σiwhen another cubical cell σjcan be assigned to be the image by V.The criteria is again to select the cubical cell with the lowest barycentric vector. The last priority type is used to choose whether a cubical cell is chosen as the source cell of a vector or the end cell of another. In this case it is used as the end cell. –iin =iout =1 Informally speaking, the vector field tries to “go” in the direction of u1. When it is not possible to go in that direction, it is used u2and, if it cannot be performed, u3is used. Hence, the flow of Vcan be thought as a waterfall going downside until it goes from left to right and, finally, rear to back. Definition 7 Let I:T3 n→{0,1}be a k-D binary image. The module ΠVector2is defined as ΠVector2={Γ,Σ,E,w 1,R,Pri,iin,iout} where –Γ=Γ0 –Σ={σi:i∈In,k(T(M))}where Mis the Morse complex built from the acyclic vector field created by ΠVector. –E=Γ0 –w1=∅ –Ris the following set of rules: –R1≡¬{ViVj}|1,σ iσjδ± ijCiCj/σiσjδ± ijV∓ ij ViVjV,0for i,j∈In,3,i= j. These rules take two cells that can be used to build a vector and create that vector. Notice that, once a cell is used in a vector, it will not be used in another. 68 R. Reina-Molina et al. –Pri =∅. –iin =iout =1 This P module cancels pairs of critical cells of ΠVector due to the previously mentioned “step” configuration. The dP scheme ΔBin3DHom =(Γ0,Σ,Π Cub,Π Vector,Π Flow,Π Vector2,R,iin,iout) with four modules defined below, computes the homology groups for the cubical complex generated from I. –Σ={Bi:i∈I−1(1)} – The modules are given by ΠCub This module builds a cubical complex from the 3D binary image I. ΠVector This module creates a maximal acyclic vector field. ΠFlow This module determines the Morse complex associated to the vector field created above. ΠVector2This module creates a maximal acyclic vector field in the Morse complex. – The distribution rules are given by –R 1≡(Cub,σ iCi→σiCi,Vector)for i∈In,3 –R 2≡Cub,δ± ij →δ± ij,Vectorfor i,j∈In,3 –R 3≡Cub,U+ ij,d→U+ ij,d,Vectorfor i,j∈In,3and 1 ≤d≤3 –R 4≡Cub,dip →dip,Vectorfor i∈In,3and 0 ≤p≤3 –R 5≡(Vector,σ i→σi,Flow)for i∈In,3 –R 6≡(Vector,Ci→Ci,Flow)for i∈In,3 –R 7≡Vector,δ± ij →δ± ij,Flowfor i,j∈In,3,i= j –R 8≡Vector,dip →dip,Flowfor i∈In,3and 0 ≤p≤3 –R 9≡Vector,V± ij →V± ij ,Flowfor i,j∈In,3,i= j –R 10 ≡Flow,˜σi˜ Ci→σiCi,Vector2for i∈In,3 –R 11 ≡Flow,˜ δ± ij →δ± ij,Vector2for i,j∈In,3,i= j –R 12 ≡({V}|Vector2,σ i→σi,Flow)for i∈In,3 –R 13 ≡({V}|Vector2,Ci→Ci,Flow)for i∈In,3 –R 14 ≡{V}|Vector2,δ± ij →δ± ij,Flowfor i,j∈In,3,i= j –R 15 ≡{V}|Vector2,V± ij →V± ij ,Flowfor i,j∈In,3,i= j –iin =Cub, iout =Vector2 The modules ΠVector and ΠVector2are the only modules that do not rely on the framework. In other words, the user of the framework only has to deal, most of the times, with writing modules for calculating acyclic vector fields, as the tasks for building the cubical complex associated to the binary image and the building of the Morse complex associated to it, are left to the framework. Membrane parallelism for DMT in digital images 69 Theorem 6 The dP scheme ΔBin3DHom defined above calculates the homology groups for a given binary 3D image I:T3 n→{0,1}in O(log2n)computation steps. Proof (Remark) The result is proved by following the behavior of the dP scheme ΔBin3DHom for an input image I. First of all, the input alphabet for ΠCub is set up from Ias ΠCub is the input module. No distribution rules can be selected until one module reaches its halting stage. By Theorem 4,ΠCub halts in two steps. Then, the only distribution rules that can be selected are R1,R2,R3and R4. After it application, the module ΠVector starts calculating a vector field. Rules R1inΠVector buildavectorfromacellσitoaco-faceσj.Theprioritiesensures that all the p-cells are mapped by V, the vector field, to a p+1-cell, prioritizing that the vector joining its barycenter has as lowest index as possible. Also, the rules ensure that the vector field is acyclic, as all the vectors are “parallel” to the vectors in the canonicalbase inR3withpositivedirection.Thisworkisdone inonlyone computation step. This way of building a vector field creates two kind of critical cells. First of all, there are critical cells due to homology generators. Second, there are critical cells due to “steps” in the cubical complex. In Fig. 1a, the cells P11,P16and P11,P12 are in such position. It looks somewhat an “step” in a stair in a bottom view. This configuration of critical cell is due to the priority relation among rules in ΠVector. Hence, we have two kinds of critical cells, those representing an homology generator and those cells that must have to be canceled in later computation steps. Once ΠVector halts, the only distribution rules that can be selected are R5to R9, sending objects to ΠFlow. Concretely, those rules send cells (σi), boundary map (δ± ij), critical cells (Ci) and vector field V± ij . Module ΠFlow calculates the Morse complex associated to the previous vector field in, at most, 4+3log2(n−1). To prove this note that the longest V-path in a 3D cubical complex with, at most, (n−1)31-cells have length, at most, (n−1)3and, then, use Theorem 5. Once the Morse complex is built, it is sent to the module ΠVector2by distribution rules R10 and R11. This module builds an acyclic vector field that “removes” all the cells that do not represent an homology generator. This computation is completed in only one step. If some vector field is built, which means that there are cells that can be canceled, at least one object Vcomes from the environment that ensures the selection and application of distribution rules R12 to R15, sending the Morse complex and the gradient vector field to the module ΠFlow, where another Morse complex is built and sent to ΠVector2. Once the second Morse complex has came to ΠVector2the dP scheme ΔBin3DHom reaches the halting state, as no other distribution rule can be selected. This takes place as the only critical cells that are not homology generators are 1-cells that are in a “step” configuration, and those cells are removed by the second vector field. Formally, we have reduced the initial cubical complex to another complex with null differential, which means that all the cells are generators of the homology groups. Recall that the sum of the computation steps required by each module is, in the worst case, at most 12 +6log2(n−1), which is O(log2n). The full calculation process of homology groups can be followed in Fig. 1.The module ΠCub builds the cubical complex in Fig. 1a in two computation steps. Then, 70 R. Reina-Molina et al. the module ΠVector builds the acyclic vector field in Fig. 1b in one computation step. Following, the module ΠFlow constructs the (first) Morse complex in 7 steps (because the longest V-path in the complex has length 6). Figure 1c shows the two kinds of critical cells. On one hand, P19 is a critical 0-cell due to an homology generator. On the other hand, P12 and P16 are 0-cells due to “step” configuration. P5,P6is a critical 1-cell due to “hole” while P11,P12and P15,P16are critical 1-cells originated by a “step” configuration. Next, the module ΠVector2specify, in one computation step, another acyclic vector field in the complex in Fig. 1d which leads to the computation in 5 computation steps of the homological complex in Fig. 1e, that has null differential, so that all the cells are homology generators. 6 Conclusions and future work In this paper a theoretical massively parallel framework to develop algorithms in homological analysis of 3D imagery has been presented. We focus our study in 3D digital images, due to the fact that “linear topology” or homology information can incorporate torsion algebraic data in higher dimension beyond classical DMT context. In fact, the computatibity problem of torsion homology subgroups for kD(k≥4) digital objects is a priority at long-term for this incipient and worthwhile research line. At short term, we intend to implement this set of algorithms in massively parallel hardware as graphics cards supporting CUDA. We also work on applications of our parallel theoretical framework to other image processing tasks, as segmentation or skeletonization. We are also interested in the representational power of our framework, as a parallel implementation of the works of Real and Molina-Abril [21,28]. References 1. Bauer, U., Lange, C., Wardetzky, M.: Optimal topological simplification of discrete functions on surfaces. Discrete Comput. Geom. 47(2), 347–377 (2012) 2. Berciano, A., Molina-Abril, H., Real, P.: Searching high order invariants in computer imagery. Appl. Algebra Eng. Commun. Comput. 23(1–2), 17–28 (2012) 3. Carnero, J., Díaz-Pernil, D., Gutiérrez-Naranjo, M.A.: Designing tissue-like P systems for image segmentation on parallel architectures. In: del Amor, M.A.M., P˘aun, Gh., de Mendoza, I.P.H., Romero- Campero, F.J., Cabrera, L.V. (eds.) Ninth Brainstorming Week on Membrane Computing, pp. 43–62. Fénix Editora, Sevilla, Spain (2011) 4. Christinal,H.A., Díaz-Pernil,D., Real, P.: Segmentationin 2Dand 3Dimage usingtissue-likeP system. In: Bayro-Corrochano, E., Eklundh, J.O. (eds.) Progress in Pattern Recognition, Image Analysis, Computer Vision, and Applications, Lecture Notes in Computer Science, vol. 5856, pp. 169–176 (2009) 5. Christinal, H.A., Díaz-Pernil, D., Real, P.: Region-based segmentation of 2D and 3D images with tissue-like P systems. Pattern Recogn. Lett. 32(16), 2206–2212 (2011) 6. Christinal, H.A., Díaz-Pernil, D., Real, P.: P systems and computational algebraic topology. Math. Comput. Model. 52(11–12), 1982–1996 (2010) 7. Cazals, F., Chazal, F., Lewiner, T.: In: Proceedings of the Nineteenth Annual Symposium on Computational Geometry, pp. 351–360. ACM (2003) 8. Díaz-Pernil, D., Christinal, H.A., Gutiérrez-Naranjo, M.A., Real, P.: Using membrane computing for effective homology. Appl. Algebra Eng. Commun. Comput. 23(5–6), 233–249 (2012) Membrane parallelism for DMT in digital images 71 9. Díaz-Pernil, D., Gutiérrez-Naranjo, M.A., Molina-Abril, H., Real, P.: Designing a new software tool for digital imagery based on P systems. Nat. Comput. 11(3), 381–386 (2012) 10. Díaz-Pernil, D., Gutiérrez-Naranjo, M.A., Pérez-Jiménez, M.J., Riscos-Núñez, A.: A linear-time tissue P system based solution for the 3-coloring problem. Electr. Notes Theor. Comput. Sci. 171(2), 81–93 (2007) 11. Díaz-Pernil, D., Gutiérrez-Naranjo, M.A., Pérez-Jiménez, M.J., Riscos-Núñez, A.: Solving subset sum in linear time by using tissue P systems with cell division. In: Mira, J., Álvarez, J.R. (eds.) IWINAC (1). Lecture Notes in Computer Science, vol. 4527, pp. 170–179 (2007) 12. Díaz-Pernil, D., Gutiérrez-Naranjo, M.A., Real, P., Sánchez-Canales, V.: Computing homology groups in binary 2D imagery by tissue-like P systems. Roman. J. Inf. Sci. Technol. 13(2), 141–152 (2010) 13. Forman, R.: Morse theory for cell complexes. Adv. Math. 134, 90–145 (1998) 14. Gunther, D., Reininghaus, J., Hotz, I., Wagner, H.: Memory-efficient computation of persistent homology for 3D images using discrete morse theory. In: 2011 24th SIBGRAPI Conference on Graphics, Patterns and Images (Sibgrapi), pp. 25–32 (2011) 15. Gyulassy, A., Bremer, P.T., Hamann, B., Pascucci, V.: A practical approach to Morse–Smale complex computation: scalability and generality. IEEE Trans. Vis. Comput. Graph. 14(6), 1619–1626 (2008) 16. Kaczy´nski, T., Mischaikow, K., Mrozek, M.: Computational homology. Applied Mathematical Sciences. Springer, Berlin (2004) 17. Lewiner, T., Lopes, H., Tavares, G.: Applications of Forman’s discrete Morse theory to topology visualization and mesh compression. IEEE Trans. Vis. Comput. Graph. 10(5), 499–508 (2004) 18. Martín-Vide, C., Pazos, J., P˘aun, Gh., Rodríguez-Patón, A.: A new class of symbolic abstract neural nets:tissue Psystems.In: Ibarra,O.H., Zhang,L. (eds.) COCOON.Lecture Notesin Computer Science, vol. 2387, pp. 290–299 (2002) 19. Martín-Vide, C., P˘aun, Gh., Pazos, J., Rodríguez-Patón, A.: Tissue P systems. Theoret. Comput. Sci. 296(2), 295–326 (2003) 20. Molina-Abril, H., Real, P.: Homological optimality in discrete Morse theory through chain homotopies. Pattern Recogn. Lett. 33(11), 1501–1506 (2012) 21. Molina-Abril, H., Real, P.: Homological spanning forest framework for 2D image analysis. Ann. Math. Artif. Intell. 64(4), 385–409 (2012) 22. Nithin, S., Vijay, N.: Parallel computation of 3D Morse–Smale complexes. Comput. Graph. Forum 31(3), 965–974 (2012) 23. Peña-Cantillana,F., Díaz-Pernil,D., Berciano, A., Gutiérrez-Naranjo, M.A.: A parallel implementation of the thresholding problem by using tissue-like P systems. In: Real, P., Díaz-Pernil, D., Molina-Abril, H., Berciano, A., Kropatsch, W.G. (eds.) CAIP (2). Lecture Notes in Computer Science, vol. 6855, pp. 277–284 (2011) 24. Peña-Cantillana, F., Díaz-Pernil, D., Christinal, H.A., Gutiérrez-Naranjo, M.A.: Implementation on CUDA of the smoothing problem with tissue-like P systems. Int. J. Nat. Comput. Res. 2(3), 25–34 (2011) 25. P˘aun, A., P˘aun, Gh.: The power of communication: P systems with symport/antiport. New Gener. Comput. 20(3), 295–306 (2002) 26. P˘aun, G., Pérez-Jiménez, M.J.: Solving problems in a distributed way in membrane computing: dP systems. Int. J. Comput. Commun. Control 5(2), 238–252 (2010) 27. Peterka, T., Ross, R., Gyulassy, A., Pascucci, V., Kendall, W., Shen, H.-W., Lee, T.-Y., Chaudhuri, A.: Scalable parallel building blocks for custom data analysis. 2011 IEEE Symposium on Large Data Analysis and Visualization (LDAV), pp. 105–112 (2011) 28. Real, P., Molina-Abril, H., Kropatsch, W.: Homological tree-based strategies for image analysis. Comput. Anal. Images Patterns (LNCS) 5702, 326–333 (2009) 29. Reininghaus, J., Hotz, I.: Combinatorial 2d vector field topology extraction and simplification. In: Topological Methods in Data Analysis and Visualization, pp. 103–114. Springer, Berlin, Heidelberg (2011) 30. Reininghaus,J.,Lowen, C.,Hotz,I.:Fastcombinatorialvectorfieldtopology.IEEETrans. Vis.Comput. Graph. 17(10), 1433–1443 (2011) 31. Ritter, G.X., Wilson, J.N., Davidson, J.L.: Image algebra: an overview. Comput. Vis. Graph. Image Process. 49(3), 297–331 (1990)