A Parallel Implementation for Computing the Region-Adjacency-Tree of a Segmentation of a 2D Digital Image

A Parallel Implementation for Computing the Region-Adjacency-Tree of a Segmentation of a 2D Digital Image Fernando D´ıaz-del-R´ıo1, Pedro Real1(B), and Darian Onchis2 1H.T.S. Informatics’ Engineering, University of Seville, Seville, Spain [email protected],[email protected] 2Faculty of Mathematics, University of Vienna, Vienna, Austria [email protected] Abstract. A design and implementation of a parallel algorithm for computing the Region-Adjacency Tree of a given segmentation of a 2D digital image is given. The technique is based on a suitable distributed use of the algorithm for computing a Homological Spanning Forest (HSF) structure for each connected region of the segmentation and a classical geometric algorithm for determining inclusion between regions. The results show that this technique scales very well when executed in a multicore processor. Keywords: Digital image ·Segmentation ·RAG ·Parallel algorithm 1 Introduction An important high level processing task in image understanding is to find a topological and structured description of the segmentation previously performed at low level processing, mostly independent of size and contrast characteristics of the extracted regions. The most usual representations are graph-based and the reasoning used in the region-node connectivity calculus involves two topological properties: adjacency and inclusion. The Region Adjacency Graph (RAG) [18–20] can be considered as the germ notion of all these models and it is composed of a set of nodes, one by region of the image, and there is an edge between two nodes if and only if the two corresponding regions are neighbors. Within the traditional context of square pixel-based 2-dimensional digital images, the most usual adjacency relationships between regions employed for the RAG are 4 or 8-neighborhoods. In order to exclusively highlight the ambient isotopic property “to be surrounded by” as adjacency relationship between regions or boundaries of regions, the notion of Region-Adjacency Tree (or RAG tree, for short) (also called homotopy tree, inclusion tree or topological tree) is created [6,20,21]. Restricted to binary 2D digital images, the RAG tree contains all the homotopy type information of the foreground object (black object) but the converse is, in general, not true [21]. Aside from image understanding applications [2], RAG trees have encountered exploitation niches in geoinformatics, rendering, dermatoscopics image, biometrics,... [3,4]. Common operations for modifying RAGs are splitting and merging of regions. Most of the methods using these split-and-merge techniques are employed in “improving the quality of initial oversegmentations”. Some split-and-merge RAG-based methods consider low-level information of the image for the region merging task [10,12,13,22]. Others take into account not only local knowledge of the image but also the extraction of “global” information about image models [1,8,9,23]. The contribution of this paper consists of the design and implementation of an parallel algorithm for computing the RAG of a given segmentation, based on a suitable distributed use of the algorithm for computing an HSF structure for each 4-connected region. The paper has the following sections. Section 2is devoted to recall the machinery for computing the structures needed for HSF construction already given in [7]. Next, a parallel algorithm for constructing RAG starting from a HSF is showed in Sects. 3and 4. Section 5is devoted to show the advantages of our parallel implementation and to present the scalability results, and finally the conclusions are summarized. 2 About HSF Framework for RAG Tree Computation In this paper, we specify an algebraic-topological technique for parallel computing a RAG of a given segmentation. This modus operandi have already used in [7,16] for developing a topologically consistent framework for parallel computation in 2D digital image context. Succinctly, starting from a 2-dimensional abstract cell (cubical) complex C2analogous of a digital image Iin which the pixels are the 0-cells, we compute in parallel another cell complex D2topologically equivalent to the first one, having the same 0-cells but in which the number of face and coface relationships in the canonically associated partial-order set of cells is reduced to a “minimum”. Taking into account that C2is contractible, D2will have only one cell (critical 0-cell) with no coface neighbors. In [7], it is employed the name of Morse spanning forest (MrSF)forD2. Let us note that we use the nickname of MrSF instead of that of MSF, due to the fact that MSF is an usual abbreviation for Minimum Spanning Forest in graph theory. If the concrete strategy for “cutting” coface connexions used in any processing unit associated Pto each 0-cell mainly depends on the localization in Pof the pixels of the “foreground” and the “background” of a given segmentation S={Ri}k i=0 of I,theMrSF D2is, in principle, a good candidate structure for topologically analyzing Sand, in particular, to compute its associated RAG. A much better and natural strategy for the RAG-problem is to compute a sort of optimal MrSF (called HSF in [7]) for each region Riof Sand to find an efficient way for “minimally” connecting those HSF structures (interpreted as connected components) between them using adjacency relationships. By “minimal” connectivity, we mean to use only one edge for connecting two adjacent connected components (nodes of the RAG). Since we use an object-based digital image analysis, the design of an HSF framework for 8-connected digital objects can straightforwardly be done simply by slightly changing the external aspect of the processing units employed. An important topological property in the 2D digital setting is that, given a black region-of-interest Rof an image, a white maximal connected region surrounded by R(“no there” or hole of Rin homology theory) is always 4-connected independently whether Ris 4 or 8-connected. We suppose here that all the regions Riof Sare 4-connected. Although this is, in principle, a strong topology restriction that almost never is accomplished in real image segmentation, a suitable mathematical morphological preprocessing of the image can solve this obstacle. Algorithms of the next section part from the information structures that contain the HSF and the contours of the regions of interest. They can be computed (see [7]) on a time complexity order (for a parallel execution) of O(log c), being cthe number of corners detected on a 2D digital image. To begin with, a brief summary of the extraction of the HSF and the contours follows. From now on, in order to favor the computational methods and techniques allowing the practical implementation of a parallel architecture, we reduce the amount of mathematical results to a minimum. In [7] the processing was focused on detecting in an exhaustive manner the outer and inner contours of the different curves that are present at a presegmented image. Thus, global information could be extracted. The implementation in [7] was essentially a sort of CCL algorithm (Connected-Component Labeling, see [11,15]) only applied to those pixels that form corners. This process was divided in three sections. The first part built the MrSF forest and determine some matrices containing the image corners and their characteristics. The second part was devoted to discover the outer and inner contour curves by following consecutive corners in vertical and horizontal directions. The result was two matrices Hcycle,V cycle, containing the corner tags of the cycles that stood for the detected curves, and other structures with the information about the sink and sources cells (for the present work they have been comprised into the vector crit idx). This part seeks for the global information relating the different corners that form the objects inside an image. Finally, a third section finished with the construction of a HSF for the foreground ROI of a binary image, by finding the appropriate pairs of sources and sinks, and by doing the arrow reversing process in parallel for each curve and each line. With previous processing, the contour information of the different image regions is comprised into two matrices. Matrix Crow,col,type contains one row per corner and three columns. Column 1 indicates the row index of the corner (or alternatively the Y coordinate). Column 2 holds the column index (or alternatively the X coordinate), and the last column contains the type (number of quadrant, with negative sign if it is external to the current curve, that is, a number among 1, 2, 3, 4, -1, -2, -3, -4). The matrix with the information of closed curves is Hcycle. This matrix can be reduced to have one row per curve and as many columns as corners have each curve. The numbers contained in the elements of this matrix are indexes to the rows of Crow,col,type. Finally vector crit idx indicates the row index of each critical cell (according to the order given by Matrix Crow,col,type). Two counters are obtained by previous algorithms: nof curves, which is the number of curves detected, and the vector nof corn curve, whose elements contain the number of corners of each curve. In order to clarify this notation, a basic image with two closed curves (each one with a hole) and five curves (one of them is the image border) is depicted in Fig. 1. In this case, nof curves is 5 and the five elements of nof corn curve contains 4. The elements of these matrices are given in Fig. 2. One of the rows of Hcycle includes its corner sequence (16,11,10,5); in particular the first row. Its critical 0-cell has been represented by corner 5. Similarly for the other curves. Fig. 1. A simple image with three curves (one of them is the image border). The crosses of the different segments with a dotted vertical line that parts from critical cell (tagged as 3) serve to detect inclusions. An important fact of Hcycle is that the corner tags are stored in a sequential order following a “contour” along the curve. This will permit us to find the number of times that a set of segments move in one direction or in the opposite. Moreover, the first hop of this contour is always a horizontal segment, which is important to extract a set of only horizontal pairs of corners (segments). 3 Parallel Computation of the Inclusion Relations for the MrSF Structures Using these matrices, the parallelization of an inclusion searching algorithm can be done in an efficient way, as depicted by the pseudo code of Algorithm 1 (see Fig. 3). Our implementation is inspired by the classical ray-crossing algorithm for the point in polygon problem [14]. However in our case, this problem can Fig. 2. Input matrices for a simple image with three curves. be computed more efficiently because instead of arbitrary polygons, only horizontal and vertical concatenated segments exist. Thus only simple computer operations like comparisons are sufficient to detect crossings. Furthermore, the parallel processing becomes more evident. Taking into account previous considerations, the codification style and the notation of Algorithm 1, is different from the usual iterative/procedural conventions to describe a pseudo-code. Most studies use an OpenMP-like notation or the for each paradigm to express more clearly the thread parallelism that can be exploited (see references along this work). As our aim is to describe the inherent parallelism that can be exploited for extracting the RAG Tree, our notation follows that of OCTAVE/MATLAB, so time complexity orders can be easily computed. The last stage to compute the RAG Tree consists of the conversion of matrix Tinc into another simplified matrix TTM that holds the topological tree (also in a matrix form). This is done by Algorithm 2 (Fig. 4)inafewsteps. The first line computes the number of columns of the TTM, which is given by those curves that have no other inside, by counting the number of columns of Tinc that are empty. If the level of “depth” of a curve is the number of curves that it encloses, then the maximum level of “depth” of all the curves plus one gives the number of rows (tree levels) of the TTM (see step 1). Step 2 contains the key point that will serve to order the rows of Tinc, in order to compute the RAG tree: each row of matrix Cidx contains the right sequences. This is due that Fig. 3. Algorithm 1: Pseudo code for the algorithm that computes the inclusion relation of a set of curves. the minimum row coordinates of a group of nested curves identify the order of inclusions (when sorting them). Next line reorders Tinc using these sequences Cidx, meanwhile step 4 adds the own curve numbers in a new column. Finally the last step 5 extracts those lines that are the leaves of the RAG tree and transpose the result to obtain TTM. Thus, those curves that do contain any other curve are forgotten in the final TTM matrix. Fig. 4. Algorithm 2: Pseudo code for the algorithm that computes the inclusion relation of a set of curves. 4 An Example of RAG Tree Parallel Computation In order to explain those Algorithms, let us examine the case of curve 5, whose critical cell is tagged with a 3 (Fig. 1). The searching of the curves that surround any other is based on detecting the crosses of a vertical line with the critical cell of each curve. In Fig. 1, the vertical line is depicted on the column 4, where critical cell of curve 5 was placed by the HSF construction [7]. First step extracts the tags of both extremes of any possible segment, which give us: Fextreme =(1610;212;113;414;315) Sextreme =(115;717;618;919;820) After that, those tags are converted into matrices Cleft,Crigth, which are filled with the column indexes of all the pair of corners that represent horizontal segments. In our case: Cleft =(113;711;810;37;46) Crigth = (13 1; 11 7; 10 8; 7 3; 6 4) Step 4 extract the tag of the critical cell of this curve, that is crit cell (5)=3. Coordinates rcrit cell =10 and ccrit cell =4 are obtained by accessing to Crow,col,type. Step 5 gives us the indexes over Cleft and Crigth (for each curve) of those segments that cross the vertical line on clockwise (CW) or anticlockwise (AW) directions. In our case, the find function returns: CWidx = (1; Empty matrix; Empty matrix; 1; Empty matrix) AWidx = (2; Empty matrix; Empty matrix; 2; Empty matrix) Matrices ECWidx,EAWidx are provided for those cases where the column coordinate is the same for the critical cell and a segment extreme. For those cases, the proper corners that determine this crossing are selected (step 6). At next step 7, the tags of the crosses are found, by indexing on Fextreme and Sextreme. For our example, CWidx tag =(10; Empty matrix; Empty matrix; 14; Empty matrix); AWidx tag =(6; Empty matrix; Empty matrix; 4; Empty matrix) The row coordinates of previous tags CWidx tag and AWidx tag are found by accessing the first column of Crow,col,type in the operations of step 8. This gives us: CWrow =(13; Empty matrix; Empty matrix; 7; Empty matrix) AWrow =(1; Empty matrix; Empty matrix; 11; Empty matrix) The first two operations of step 9 count the number of crosses (non zero elements). Function nnz returns the values: sup crosses =(1;0;0;1;0) inf crosses =(1;0;0;1;0) The other operation saves the minimum row coordinates of all the segments, which will be used to find the order in which a curve is surrounded by others (see next algorithm). Finally step 10 computes the vector inclusion, that is, the absolute of the element-by-element product (sup crosses .* inf crosses), because sup crosses and inf crosses can be 1 or -1 when a segment crosses the vertical line. Note that a curve include the current cell, if and only if there is one cross with this vertical line above this cell and only one below it. Four our case, the fifth row of inclusion tells us that curves 1 and 4 surround the current curve 5: inclusion (5, :) = (1 ; 0 ; 0 ; 1 ; 0) This vector is multiplied element-by-element by the numbers (1: nof curves), so at step 11 the matrix Tinc is marked with these numbers. In our example, the fifth row is set to (1 0 0 4 0). In addition, the matrix Tinc min row is filled with the minimum row coordinates of all the segments of surrounding curves. At the same time, these values serve to fill a vector rows to be deleted with the number of times that a curve include any other. This vector will help to delete those rows that must not appear at the final RAG tree. In our example: Tinc = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 00000 10000 12000 10000 10040 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ This means that curve 2 is surrounded by 1 (the image borders, in row 1), curve 3 by 1 and 2 (in row 3), curve 4 by 1, and curve 5 by 1 and 4 (row 7). Using Algorithm 2, it is obvious that the resultant matrix for this example TTM results: ⎛ ⎝ 11 24 35 ⎞ ⎠ 5 Implementation and Testing Results A complete implementation in OCTAVE/MATLAB has been built following previous pseudo codes. In order to check its scalability, a non-fully functional implementation of Algorithm 1 (because its time computing is very much bigger than that of Algorithm 2) has been done in C++ using OpenMP directives. Using previous notation, the next advantages are done evident: (a) It can be easily demonstrated that Algorithm 1 scales well whit the number of Processing Elements of the target parallel architecture (like multi or manycore, GPUs, SIMD kernels, ...). This also eases the codification for a parallel computer. For similar reasons, those sentences that can be executed in parallel are grouped in the same step (using lines preceded with letters). (b) OpenMP codes can be written directly through this notation, in most of the cases by converting the matrix processing into nested loops, and preceding the outer loop (usually devoted to the rows) with the directive #pragma omp parallel for. (c) The memory access pattern (which is in most occasions [5] a critical point for the performance that can be achieved in GPUs) can be clearly observed. (d) Matrix operations that cannot be done in an element-by-element manner (like matrix inversions, matrix multiply, etc.) are avoided. (e) Only plain matrices have been used, preventing complex data structures (trees, chained lists, stacks, etc.), very inefficient when running on massive data parallel computer architectures. (f) Finally, Algorithms have no conditional sentences. The avoidance of conditional sentences in the hot spot zones prevents the so-called thread divergence for GPUs, which is one the main reasons why the performance on this platforms diminishes [5].