1. Introduction
Thinning (or skeletonization) is the process of reducing an object to its skeleton. The topology preserving skeleton may be informally defined as a thinned subset of the object that retains the same topology of the original object and often many of its significant morphological properties. Thinning is a very active research area thanks to its ability of reducing the amount of information to be processed for example in medical image analysis and visualization as well as simplifying the development of pattern recognition or computeraided diagnosis algorithms. Hence, it is not surprising that thinning gained a pivotal role in a wide range of applications. An exhaustive review of the literature is beyond the scope of this work. Two dimensional skeletons have been used for digital image analysis and processing, optical character and fingerprint recognition, pattern recognition and matching, and binary image compression since a long time ago, see for example the survey paper [1]
. More recently, three dimensional skeletons have been widely used in computer vision and shape analysis
[2], in computer graphics for mesh animation [3] and in computer aided design (CAD) for model analysis and simplification [4], [5] and for topology repair [6].There is also a vast literature of applications of skeletons in medical imaging. They have been used for route planning in virtual endoscopic navigation [7], for example in virtual colonoscopy [8][12] or bronchoscopy [13]. Skeletons have also been an important part of clinical image analysis by providing centerlines of tubular structures. In particular, there is a large body of literature showing applications of skeletons to blood veins centerline extraction from angiographic images [12], [14][20], and intrathoracic airway trees classification for the evaluation of the bronchial tree structure [21]. Also protein backbone models can be produced with techniques based on skeletons [22]. Furthermore, many computeraided diagnostic tools rely on skeletons. For example, skeletons have been used to identify blood vessels stenoses [23][25], tracheal stenoses [26], polyps and cancer in colon [27] and left atrium fibrosis [28].
There are medical applications of skeletons where topology preservation is essential. Non invasively determine the threedimensional topological network of the trabecular bone [29] is a good example. Indeed, many studies demonstrate that the elastic modulus and strength of the bones is determined by the topological interconnections of the bone structure rather than the bone volume fraction [30], [31]. Therefore, topological analysis plays a fundamental role in computeraided diagnostic tools for osteoporosis [30].
Topology preserving thinning is non trivial and a vast literature, briefly surveyed in Sec. 1.1, has been dedicated to this topic. In particular, thinning by iteratively removing simple points [32] is a widely used and effective technique. It works locally and for this reason is efficient and easy to implement.
While reading the literature one may notice that thinning algorithms are claimed to be “topology preserving,” even though in most cases a precise statement of what that means is left unaddressed. This paper uses homology theory [33] to rigorously define what the virtue of being topology preserving actually consists of. This theory is less intuitive than the concepts used so far, including simple homotopy type [34], but exhibits some important theoretical and practical advantages that will be highlighted later in the paper. We remark that a homological definition of simple points has already been used in the context of skeletonization in [16], [35], but only in the case of cubical complexes. This paper generalizes this idea to cell complexes that are more general than cubical complexes. There are many applications that would benefit from an algorithm that deals with general unstructured simplicial complexes [p. 35][36]. In fact, the geometry of threedimensional objects is frequently specified by a triangulated surface, obtained for example by using an isosurface algorithm as marching cubes [p. 539][36], [12] applied on voxel data from computed tomography, magnetic resonance imaging or any other threedimensional imaging technique. Another possibility is to obtain the triangulations from the convex hull of point clouds provided for example by 3d laser scanners. Triangulated surfaces offer two potential advantages over voxel representation. They allow to adaptively simplify the surface triangulation, see for example Fig. 16.20 in [p. 549][36]. They also allow to visualize and edit the object efficiently with offtheshelf software (for example the many visualization and editing tools for stereo lithography) and without the starcase artifacts typical of voxel representation of objects with curved boundary. One may even easily print the object with additive manufacturing technology (i.e. 3d printers).
Another issue that arises reading the literature is that many different definitions of topology preserving skeleton exist. In some papers, the skeleton is obtained by removing simple pairs in the spirit of simple homotopy theory by what is well known as collapsing in algebraic topology [33]. The resulting skeleton, if no other constraints are used, has a lower dimension with respect to the input complex. On the contrary, this paper assumes that the skeleton is always a solid object of the same dimension as the initial complex. The difference is highlighted in Fig. 1.
In this paper the skeleton of a given complex is defined as a subset that is obtained from after removing a sequence of top dimensional cells. We require that the homology [33] of the initial complex is preserved during this process. In particular, a top dimensional cell can be safely removed if this does not change the homology of the complement of . Fig. 2 provides an intuitive explanation why the last requirement is desirable. This additional requirement, to the best of our knowledge, is not documented in other papers. We call this cell a simple cell, which is a generalization of the idea of simple points [32] in digital topology.
Clearly, in nontrivial cases, the skeleton is not unique.
Resorting to explicit homology computations to detect simple points as in [16], [35], [32] is quite computationally intensive, as the worstcase complexity of homology computations is cubical, see also the discussion in [32]. In this paper, we introduce a much more efficient solution by exploiting the idea of tabulated configurations, i.e. acyclicity tables, that are described in detail in Sec. 3.
Usually, a skeleton also requires to preserve the shape of the object. In this paper we show some very simple proof of concept idea how to preserve both homology and shape. Of course, this is just an example to illustrate how the idea of acyclicity tables can be used together with some additional techniques that guarantee shape preservation.
The rest of the paper is organized as follows. In Section 1.1 the prior work on thinning algorithms is surveyed. Section 1.2 analyzes the original contributions of the present paper. In Section 2 the property of being a homology preserving thinning is rigorously stated. In Section 3 the concept of acyclicity tables is introduced, whereas, in Section 4, the topology preserving thinning algorithm is presented. Section 5 discusses the results of the thinning algorithm on a number of benchmarks and, finally, in Section 6 the conclusions are drawn.
1.1. Prior work
There are hundreds of papers about thinning. Most of them fall into two categories. On one hand, there are papers using morphological operations like erosion and dilatations to obtains skeletons, see [37] and references therein. They do not guarantee topology preservation in general. The others use the idea of removing the so called simple points from the given cell complex, see [16], [35], [32]. Without pretending to be exhaustive, in the following we resume previous results.
1.1.1. 2dimensional images
1.1.2. 3dimensional cubical complexes
In case one wants to skeletonize three (or higher) dimensional images, there are much less papers available in literature. Most of them rely on case study, see [44][57]. The problem is that it is hard to prove that a rulebased algorithm is general, i.e. it removes a cell if and only if its removal does not change topology. In 3d there are more than 134 millions possible configurations for a cube neighborhood and only treating correctly all of them gives a correct thinning algorithm. References [16], [35], [32] use explicit homology computations to detect simple points.
There are a number of papers presenting thinning algorithms for 3dimensional images in which Euler characteristic is used to guarantee topology preservation, see for example [43], [54] and references therein. The problem is that Euler characteristic is a rather raw measure of topology and it is not sufficient to preserve topology in general for 3dimensional cubical complexes. For three dimensional images one needs to use both Euler characteristic and connectivity information to preserve topology, but this is not sufficient for four dimensional images.
1.1.3. 2dimensional simplicial complexes and cell complexes
All the strategies presented so far are applicable only to cubical grids (pixels, voxels, …). To our best knowledge, there are just a few papers dealing with 2d grids that are not cubical, and they are restricted to 2d binary images modeled by a quadratic, triangular, or hexagonal cell complex, see [58][61]. The main reason for the lack of results on general 2d simplicial complexes may be the absence of regularity in unstructured simplicial grids that makes casestudy algorithms very hard to devise and to implement. This gap in the literature is covered by the present paper.
1.1.4. 3dimensional simplicial complexes and cell complexes
To the best of our knowledge, we are not aware of algorithms that deal with unstructured 3d simplicial complexes or more general cell complexes. There are only some papers that find the 1dimensional skeleton by using the well known collapsing in algebraic topology [62][65]. Again, this gap in the literature is covered by the present paper.
1.2. Summary of paper contributions
In this Section the main novelties presented in this paper are summarized:

The claim “topology preserving thinning” is rigorously defined, for any cell complex, by means of homology theory.

A novel topology preserving thinning algorithm that removes simple cells is introduced. Conceptually this algorithm falls into the category of thinning algorithms based on simple points and generalizes all previous papers. In fact, the acyclicity tables introduced in this paper give a classification of all possible simple points that can occur in a given cell complex. Therefore, no rules are needed since all of them are encoded into the acyclicity tables.

The most important advantage of the novel approach is that acyclicity tables are automatically filled in advance, for any cellular decomposition, with homology computations performed by a computer. Therefore, once the tables are available, the implementation of a thinning algorithm is straightforward since identifying simple cells requires just queering the acyclicity table. No other topological processing is needed.

The fact that acyclicity tables are filled automatically and correctly, for all possible configurations, provides a rigorous computerassisted mathematical proof that the homologybased thinning algorithm preserves topology. It is also verified, simply by checking all acyclic configurations, that using Euler characteristic is not enough to ensure preservation of topology in 3dimensional or higher dimensional cubical and simplicial complexes. However, when one checks both Euler characteristic and that the number of connected component before and after cell removal remains one, then topology is preserved. Checking Euler characteristic together with connectivity does not suffice to preserve topology in 4d.

The acyclicity tables for simplicial complexes of dimension 2, 3 and 4 and for cubical complexes of dimension 2 and 3, that can be freely used in any implementation of the proposed algorithm, are provided as supplemental material at [66]. This way, we dispense readers to implement homology computations to produce the acyclicity tables.

The thinning algorithm, unlike the standard collapsing of algebraic topology [33], does not require the whole cell complex data structure but it uses only the top dimensional elements of the complex, with obvious memory saving.

As a proof of concept, an open source C++ implementation that works for 3dimensional simplicial complexes is provided to the reader as supplemental material at [66]. We remark that the code is optimized for readability and memory usage and not for speed.
2. Topology preserving thinning by preserving homology
When one claims that an algorithm ”preserves topology,” in order to give a precise meaning to this statement, one needs to specify which topological invariant is preserved. In the literature, the invariant is assumed to be, in most cases implicitly, the so called homotopy type [33]. The problem of this choice is that this strong topological invariant in general is not computable according to Markov [67]. This is the reason why in this paper we propose to use homology theory which is computable in place of homotopy theory, even if it is weaker than the former. Indeed, homology seems to be the strongest topological invariant that can be rigorously and efficiently computed. Therefore, every time we claim that topology is not changed, implicitly we mean that the homology is not changed.
Homology groups may be used to measure and locate holes in a given space. Zero dimensional holes are the connected components. One dimensional holes are handles of a given space, whereas two dimensional holes are voids totally surrounded by the considered space (i.e. cavities). One can look at a dimensional hole as something bounded by a deformed sphere. A space is homologically trivial (or acyclic) if it has one connected component and no holes of higher dimensions. A rigorous definition of homology groups is not presented in this paper due to the availability of rigorous mathematical introductions in any textbook of algebraic topology as [33] and the lack of space. For a more intuitive presentation for non mathematicians one one may consult [68], [69].
In this paper, we consider in particular two standard ways of representing spaces, namely the simplicial and cubical complexes. A nsimplex is the convex hull of points in general position (point, edge, triangle, tetrahedron, 4dimensional tetrahedron). A simplex spanned with vertices is denoted by . By a face of a simplex we mean the simplices spanned by a proper subset of vertices spanning . A simplicial complex is a set of simplices such that for every simplex and every face of simplex , . Pixels (2cubes) and voxels (3cubes) are widely used in image analysis. They form a Cartesian grid, that is a special case of grid where cells are unit squares or unit cubes and the vertices have integer coordinates. Even though we assume to deal with a Cartesian grid, the results presented in this paper hold also for more general grids such as a rectilinear grid, that is a tessellation of the space by rectangles or parallelepipeds that are not, in general, all congruent to each other. Therefore, we define a cubical complex as a set of cubes such that for every cube and for every being face of we have , . We want to stress here that we assume every cell to be closed, i.e. if a cell is present in a complex, so do are its faces.
In the iterative thinning algorithm presented in this paper, the top dimensional cells (simplices, voxels) are iteratively removed from the object. Homology theory is used to ensure that removing of a given cell (simplex/cube) does not change the topology of the object. If removal of a cell does not change the topology, the cell is said to be simple. Due to efficiency reasons, the homology cannot be recomputed after removing every single element. In fact, one may compute the homology of a cell complex for instance with [70] software, but the worst case computational complexity is cubical. Therefore, the main idea is to rely on the so called Mayer–Vietoris sequence [33]. Let us express the considered space as , where is a single top dimensional simplex or voxel. The sequence states that once the intersection is homologically trivial, then the homology of is the same as the homology of . This important result is the key of the method presented in this paper. In fact, it implies that, in order to preserve the homology of with respect to , one should check the homology of the intersection . In practice, this may be easily performed with the [70] software.
The main novelty of this paper is to present a different idea to speed up the computations. Let be a simplex or voxel. By we denote the boundary of , i.e. all lower dimensional cells which are entirely contained in the closure of . The idea is based on the observation that in there are not too many elements, namely

in case of triangle (dimensional simplex);

in case of tetrahedron (dimensional simplex);

in case of 4dimensional simplex (i.e. the convex hull of points in in general position);

in case of pixel (2dimensional cube);

in case of voxel (3dimensional cube).
By configuration we mean any subset of . When looking for simple cells, the configuration characterizes the way intersects the complement of the set that we aim to thin. A configuration is acyclic if its homology—computed as the homology of the corresponding chain complex—is trivial. Since the number of all possible configurations in is , where is the number of boundary elements of , one may precompute the homology of all the configurations and store them in a lookup table. In this case, homology computations are done only in a preprocessing stage and once and for all. After creating such an acyclicity table, one may instantly (i.e. in time) get the answer whether the intersection is homologically trivial or not. This is the strategy that we aim to use in the thinning algorithm. Next section shows how to use the acyclicity tables and how to obtain them automatically.
3. Use of acyclicity tables and their generation
Let us consider a generic cell of a cell complex. Let us fix an order of all boundary elements of . We consider all subsets of the set and enumerate them in the following way. For the number of a subset is . The acyclicity table of is an array of size having, at position , the value true if the configuration is acyclic and false otherwise.
Let us describe how the acyclicity table is constructed and used starting by considering a tetrahedron (3dimensional simplex), see Fig. 3a. Let us enumerate vertices, edges and faces of tetrahedron as in Fig. 3a.
We now introduce an ordering on boundary elements of the 3dimensional simplex (bold numbers are the indexes of elements in the given order):
(3.1) 
Let be the indexes of elements in the considered configuration (i.e. bold numbers corresponding to elements that are present in the configuration). The index of the configuration in the acyclicity table is computed with
The acyclicity table is automatically generated in advance as follows. All possible configurations of the elements are automatically generated and the homology group of each of these configurations is computed with the [70] software. If a configuration turns out to be acyclic, then a true is set to the place in the array corresponding to the examined configuration, false otherwise.
Example 3.1.
Suppose the 3dimensional simplex is given as input. Let and be the maximal elements in the configuration (i.e. the configuration consists of those elements and the vertices , edges , , that are the faces of ). This configuration needs to be mapped into the 3dimensional simplex model presented in Fig. 3a. Hence, we have the following mapping between vertices: 0, 1, 2, 3. It is naturally extended to the mapping on simplices. Namely, the triangle is mapped to 0,1,2 in the 3dimensional simplex model, whereas vertex is mapped to vertex 3. Therefore, the elements in this configuration are:

Vertices: 0, 1, 2, 3 (indices 1, 2, 3, 4);

Edges: 0,1, 0,2, 1,2 (indices 5, 6, 8);

Face: 0,1,2 (index 11).
Consequently, the index of this configuration is
One may check, at this position of the provided acyclicity table for 3dimensional simplices, that this configuration is not acyclic^{1}^{1}1Clearly it cannot be, since it has two connected components..
In the same spirit, we introduce an ordering for the 2dimensional simplex
(3.2) 
and for the 4dimensional simplex
(3.3) 
In the case of cubes, unlike the case of simplices, the model cube is expressly needed to specify the location of vertices in the cube^{2}^{2}2This happens because in a cube not all vertices are connected with edges as in case of simplices. Therefore, the model cube is needed to point out the incidences of the vertices.. The model for  and dimensional cubes is represented in Fig.s 3b and 3c. The ordering for the 2dimensional cube is
(3.4) 
whereas, for the 3dimensional cube (voxel) is
(3.5) 
Of course, in order to compute the index in the acyclicity table, exactly the same procedure as the one described for the 3dimensional simplex is used.
Historically, the acyclicitiy tables for cubes [71] and simplices [72] were introduced in order to speed up homology computations. In this paper we provide an even stronger result. Not only the homology of the initial set and its skeleton is the same, but one can construct a retraction from the initial set to its skeleton. The existence of retraction implies the isomorphism in homology, but the existence of retraction is a stronger property than homology preservation. We demonstrated the existence of a retraction by a bruteforce computer assisted proof, i.e. checking all acyclic configurations. Thus, the following lemma holds.
Lemma 3.2.
For every acyclic configuration in the boundary of 2, 3 or 4dimensional simplices and 2 or 3dimensional cubes (denoted as ) there exist a simple homotopy retraction from to .
At the end of this section, let us define more rigorously a simple cell.
Definition 3.3.
A cell in a complex is simple if is acyclic.
In the supplemental material, we already provide the acyclicity tables for dimensional simplices and dimensional cubes (pixels, voxels), in such a way that the reader can safely bypass the step of constructing them. We note that we do not provide tables for higher dimensional simplices or cubes, since the memory required to store them is huge^{3}^{3}3All configurations for dimensional cube require almost PB. All configurations for dimensional simplex require PB. On the contrary, the acyclicity table for the 3dimensional simplices provided as supplemental material requires no more than kB..
4. Topology preserving thinning algorithm
In this section we propose a simple thinning technique that iteratively removes simple cells. The algorithm is valid both for cubes and simplices provided that the corresponding acyclicity table is used. We want to point out that the algorithm works on top dimensional cells (cubes, simplices). Therefore—unlike the case of homological algorithms or collapsing—there is no need to generate the whole lower dimensional cell complex data structure^{4}^{4}4This structure have to be generated only locally for the boundary of a cell when checking if is simple.. The input of the algorithm consists of a list of top dimensional cells in the considered set. The output is a subset of being its skeleton.
At the beginning, we present a first version of the algorithm that preserves only the topology of . At the beginning, one searches the list to find all the cells that are simple and store them into a queue . Then, the queue is processed as long as it is not void. In each iteration, an element is removed from the queue . Then, with the acyclicity table, one has to check if is simple in the set . We want to point out that elements already removed form the considered set in previous iterations are treated as the exterior of at a given iteration. If is simple in , then it is removed from the set . In this case, all neighbors of ^{5}^{5}5A neighbor of cell/simplex is any cell/simplex such that . that are still in are added to the queue . The details of the presented procedure are formalized in Alg. 1.
We want to stress that Alg. 1 is just an illustration. It may be turned into an efficient implementation by using more efficient data structures (for instance removing from the list can be replaced by a suitable marking the considered element.) Also searching for intersection of with current should be performed by using hash tables that, for the sake of clarity, are not used explicitly in Alg. 1. Let us now discuss the complexity of the algorithm. Clearly the for loop requires operations. We assume that one can set and check a flag of every cell in a constant time. This flag indicates if a cell is removed from or not. Every cell appears in the while loop only times, where is maximal number of neighbors of a top dimensional cell in the complex. Therefore, the while loop performs at most iterations before its termination. The time complexity of every iteration is , which means that the overall complexity of the procedure is . Typically the number is a dimension dependent constant and, in this case, the complexity of the algorithm is . The same complexity analysis is valid for Alg. 2.
We now present in Alg. 2 a simple idea that enables to preserve the shape of the object in addition to its topology. We stress that the aim of this second algorithm is just to show how to couple topology and shape preservation.
In Alg. 2 there is one basic difference with respect to Alg. 1. In Alg. 2, after removing a single external layer of cells, a check is made at line 16 to determine whether all cells that remain in are already in the boundary of . Once they are, the thinning process terminates. The topology is still preserved due to line 10. The additional constraint used at line 16 of Alg. 2 is very simple and it gives acceptable results in practice. It may be easily coupled with other techniques to preserve shape already described in literature.
Finally, we discuss the situation when one wants to keep the skeleton attached to some pieces of the external boundary of the mesh. In this case, when testing whether a top dimensional cell is simple, one should consider as elements in . In other words, elements from are not considered as an interface between the object to skeletonize and its exterior.
4.1. Proofs
Now we are ready to give a formal definition of skeleton.
Definition 4.1.
Let us have a simplicial or cubical complex . A skeleton of , denoted by , is a set of top dimensional simplices or cubes such that:

is obtained from by iteratively removing top dimensional elements , provided that the intersection of with complement is acyclic. Consequently, homology groups of and are isomorphic;

There is no top dimensional element that has an acyclic intersection with complement (i.e. the process of removing such elements has been run as long as possible.)
We want to point out that sometimes, due to some deep phenomena arising in simple homotopy theory, some skeleton may be redundant. For instance it is possible to have a skeleton of a 3dimensional ball that is a Bing’s house [73] instead being a single top dimensional element. In general it is impossible to avoid this issue due to some intractable problems in topology.
In the follwing, we formally show that the skeleton obtained from Alg. 1 satisfies Def. 4.1. This fact is shown with a sequence of two simple lemmas.
Lemma 4.2.
The homology of and are isomorphic.
Proof.
The proof of this lemma is a direct consequence of the Mayer–Vietoris sequence [33]. Let be the elements removed during the course of the algorithm (enumeration is given by the order they were removed by the algorithm.) Let us show that, for every , homology of and homology of ^{6}^{6}6The difference used in the formulas in this proof is not a set theoretic difference. All the objects are assumed to contain all their faces. are isomorphic. Let us write the Mayer–Vietoris sequence in reduced homology for :
The intersection is acyclic. This is because the intersection of with the set complement is checked in the acyclicity tables to be acyclic. Once it is, also is acyclic. Therefore, is trivial. Also, since is a simplex or cube, it is acyclic. This provides being trivial (we are considering reduced homology.) Consequently from the exactness of the presented sequence we have the desired isomorphism between and . The conclusion follows from a simple induction. ∎
Lemma 4.3.
After termination of the algorithm there is no element that has an acyclic intersection with the complement.
Proof.
Let be the elements removed during the course of the algorithm (enumeration is given by the order they were removed by the algorithm.) Suppose, by contrary, that a exists such that it has an acyclic intersection with the complement. Let denotes the index of last element among that has nonempty intersection with . If , then would be put to the queue in the line 5 of Alg. 1 and removed from in the line 11 of the algorithm, since no change to its intersection with complements is made by removing . If , then after removing the intersection of with complement does not change. Therefore, it is acyclic after removing for . When Alg. 1 removes in the line 12, is added to the list and it is going to be removed in the line 11, since removing for does not affect the acyclicity of the intersection of with the complement. In both cases we showed that is removed from by Alg. 1. Therefore, a contradiction is obtained. ∎
5. Experimental results
5.1. Assessing aortic coarctation and aneurism
Skeletons can be used in computeraided diagnostic tools for coarctation and aneurism, by evaluating the transverse areas of any vessel structure, see for example [23].
5.1.1. Aortic coarctation
Aorta coarctation is a congenital heart defect consisting of a narrowing of a section of the aorta. Surgical or catheterbased treatments seek to alleviate the blood pressure gradient through the coarctation in order to reduce the workload on the heart. The pressure gradient is dependent on the anatomic severity of the coarctation, which can be determined from patient data. Gadoliniumenhanced magnetic resonance angiography (MRA) has been used in a 8 year old female patient to image a moderate thoracic aortic coarctation, see Fig. 9a. Fig. 9b shows a rendering of the 3D triangulated surface, obtained by segmenting the MRA data, which models the ascending aorta, arch, descending aorta, and upper branch vessels. The interior of the surface has been covered with 94756 tetrahedra. The skeleton of this vessel structure, obtained with Alg. 2, is shown in Fig. 9c.
5.1.2. Cerebrovascular aneurism
Cerebrovascular aneurysms are abnormal dilatations of an artery that supplies blood to the brain. Magnetic resonance imaging (MRI) has been used to image the cerebral circulation in a 47 year old female patient, see Fig. 5a. Fig. 5b shows a rendering of the 3D triangulated surface, obtained from the segmentation of the MRI data. The interior of the surface has been covered with 390,081 tetrahedra. The skeleton of this vessel structure, obtained with Alg. 2, is shown in Fig. 5c.
5.2. Analysis of pulmonary airway trees
Pulmonary arteries connect blood flow from the heart to the lungs in order to oxygenate blood before being pumped through the body. Skeletons have been used for quantitative analysis of intrathoracic airway trees in [21]. A 3D triangulated surface, shown in Fig. 6b, represents the 3D model of pulmonary airway trees of a 16 year old male patient obtained by segmenting data from computed tomography (CT) images, see Fig. 6a. The interior of the surface is covered with 236,433 tetrahedra. The topology preserving skeleton obtained by Alg. 2 is shown in Fig. 6c.
5.3. Extracting centerline for virtual colonoscopy
A 3D triangulated surface that represents the 3D model of a colon is obtained by segmenting data from computed tomography (CT) images, see Fig. 7. The interior of the surface is covered with 2,108,424 tetrahedra. The topology preserving skeleton obtained by Alg. 2, which may be used as a colon centerline to guide a virtual colonoscopy, is shown in black in Fig. 7.
5.4. Computing topological interconnections of bone structure
A 3D model of a human bone belonging to a 61 year old male patient has been obtained from a stack of thresholded 2D images acquired by Xray MicroCT scanning [74]. In particular, a region of interest (ROI) of size mmmm ( pixels, pixel size m) is selected in the trabecular region. A stack of 195 2D images has been considered, resulting in a volume of interest (VOI) of approximately mmmmmm. From this 3D model, consisting of about 2.4 millions voxels, a 3D triangulated surface has been obtained, see Fig. 8a. The interior of this surface is covered with 688,773 tetrahedra. The topology preserving skeleton obtained by Alg. 1 is shown in Fig. 8b.
5.5. Some other non medical examples.
6. Conclusion
This paper introduces a topology preserving thinning algorithm for cell complexes based on iteratively culling simple cells. Simple cells, that may be seen as a generalization of simple points in digital topology, are characterized with homology theory. Despite homotopy, homology theory has the virtue of being computable. It means that, instead of resorting to complicated rulebased approaches, one can detect simple cells with homology computations. The main idea of this paper is to give a classification of all possible simple cells that can occur in a cell complex with acyclicity tables. These tables are filled in advance automatically by means of homology computations for all possible configurations. Once the acyclicity tables are available, implementing a thinning algorithm does not require any prior knowledge of homology theory or being able to compute homology. The fact that acyclicity tables are filled automatically and correctly for all possible configurations provides a rigorous computerassisted mathematical proof that the homologybased thinning algorithm preserves topology. We believe that such rigorous topological tools simplify the study of thinning algorithms and provide a clear and safe way of obtaining skeletons.
References
 [1] L. Lam, S.W. Lee, C.Y. Suen, “Thinning MethodologiesA Comprehensive Survey,” IEEE Trans. Pattern Anal. Mach. Intell., Vol. 14, pp. 869–885, 1992.
 [2] H. Sundar, D. Silver, N. Gagvani, S. Dickinson, “Skeleton based shape matching and retrieval,” Shape Modeling International 2003, pp. 130–39, 2003.
 [3] N.D. Cornea, D. Silver, P. Min, “CurveSkeleton Properties, Applications, and Algorithms,” IEEE Trans. Vis. Comput. Graph., Vol. 13, No. 3, pp. 530–548, 2007.
 [4] T. Dey, K. Li, J. Sun, D. Cohen–Steiner, ”Computing Geometry–aware Handle and Tunnel Loops in 3D Models,” ACM Trans. Graph., Vol. 27, article 45, 2008.
 [5] T. Dey, K. Li, J. Sun, “Computing handle and tunnel loops with knot linking,” Comput. Aided. Des., Vol. 41, pp. 730738, 2009.
 [6] Q.Y. Zhou, T. Ju, S.M. Hu, “Topology repair of solid models using skeletons,” IEEE Trans. Vis. Comput. Graph., Vol. 13, pp. 675685, 2007.
 [7] D.S. Paik, C.F. Beaulieu, R.B. Jeffrey, G.D. Rubin, S. Napel, “Automated flight path planning for virtual endoscopy,” Med. Phys., Vol. 25, pp. 629–637, 1998.
 [8] Y. Ge, D.R. Stelts, J. Wang, D.J. Vining, “Computing the Centerline of a Colon: A Robust and Efficient Method Based on 3D Skeletons,” J. Comput. Assist. Tomogr., Vol 23, No. 5, pp 786–794, 1999.
 [9] M. Wan, Z. Liang, Q. Ke, L. Hong, I. Bitter, A. Kaufman, “Automatic Centerline Extraction for Virtual Colonoscopy,” IEEE Trans Med. Imaging, Vol. 21, pp. 1451–1460, 2002.
 [10] R.L. Van Uitert, R,M. Summers, “Automatic correction of level set based subvoxel precise centerlines for virtual colonoscopy using the colon outer wall,” IEEE Trans Med. Imaging, Vol. 26, No. 8, pp. 1069–1078, 2007.
 [11] M. Ding, R. Tong, S.H. Liao, J. Dong, “An extension to 3D topological thinning method based on LUT for colon centerline extraction,” Comput. Methods Programs Biomed., Vol. 94, Np. 1, pp. 39–47, 2009.
 [12] A. Huang, H.M. Liu, C.W. Lee, C.Y. Yang, Y.M. Tsang, “On Concise 3D Simple Point Characterizations: A Marching Cubes Paradigm,” IEEE Trans. Med. Imaging, Vol. 28, No. 1, pp. 43–51, 2009.
 [13] A.P. Kiraly, J.P. Helferty, E.A. Hoffman, G. McLennan, W.E. Higgins, “Threedimensional path planning for virtual bronchoscopy,” IEEE Trans. Med. Imaging, Vol. 23, No. 11, pp. 1365–1379, 2004.
 [14] K. Haris, S.N. Efstratiadis, N. Maglaveras, C. Pappas, J. Gourassas, G. Louridas, “Modelbased morphological segmentation and labeling of coronary angiograms,” IEEE Trans. Med. Imaging, Vol. 18, No. 10, pp. 1003–1015, 1999.
 [15] A.F. Frangi, W.J. Niessen, R.M. Hoogeveen, T. van Walsum, M.A Viergever, “Modelbased quantitation of 3D magnetic resonance angiographic images,” IEEE Trans. Med. Imaging, Vol. 18, No. 10, pp. 946–956, 1999.
 [16] M. Niethammer, A.N. Stein, W.D. Kalies, P. Pilarczyk, K. Mischaikow, A. Tannenbaum, “Analysis of blood vessel topology by cubical homology”, Proc. of the International Conference on Image Processing 2002, Vol. 2, 969972, 2002.
 [17] S.R. Aylward, E. Bullitt, “Initialization, Noise, Singularities, and Scale in Height Ridge Traversal for Tubular Object Centerline Extraction, IEEE Trans. Med. Imaging, Vol. 21, No. 2, pp. 61–75, 2002.
 [18] S.R. Aylward, J. Jomier, S. Weeks, E. Bullitt, “Registration and Analysis of Vascular Images,” Int. J. Comput. Vision, Vol. 55, No. 23, pp. 123–138, 2003.
 [19] M. Straka, M. Cervenansky, A. La Cruz, A. Kochl, M. Sramek, E. Groller, D. Fleischmann, “The VesselGlyph: focus & context visualization in CTangiography,” Proc. of IEEE Visualization 2004, pp. 385–392, 2004.
 [20] M. Schaap et al., “Standardized evaluation methodology and reference database for evaluating coronary artery centerline extraction algorithms,” Med. Image Anal., Vol. 13, pp. 701714, 2009.
 [21] J. Tschirren, G. McLennan, K. Palágyi, E.A. Hoffman, M. Sonka, “Matching and anatomical labeling of human airway tree,” IEEE Trans. Med. Imaging, Vol. 24, pp. 1540–1547, 2005.
 [22] M.L. Baker, S.S. Abeysinghe, S. Schuh, R.A. Coleman, A. Abrams, M.P. Marsh, C.F. Hryc, T. Ruths, W. Chiu, T. Ju, “Modeling protein structure at near atomic resolutions with Gorgon,” J. Struct. Biol., Vol. 174, No. 2, pp. 360–373, 2011.

[23]
K. Kitamura, J.M. Tobis, J. Sklansky, “Estimating the 3D skeletons and transverse areas of coronary arteries from biplane angiograms,”
IEEE Trans. Med. Imaging, Vol. 7, No. 3, pp. 173–187, 1988.  [24] C. Kirbas, F. Quek, “A review of vessel extraction techniques and algorithms,” ACM Comput. Surv., Vol. 36, No. 2, pp. 81–121, 2004.
 [25] Y. Yang, L. Zhu, S. Haker, A.R. Tannenbaum, D.P. Giddens, “Harmonic skeleton guided evaluation of stenoses in human coronary arteries,” Med. Image Comput. Assist. Interv., Vol. 8, pp. 490–497, 2005.
 [26] E. Sorantin, C.S. Halmai, B. Erdôhelyi, K. Palágyi, L.G. Nyúl, L.K. Ollé, B. Geiger, F. Lindbichler, G. Friedrich, K. Kiesler, “SpiralCTbased assessment of tracheal stenoses using 3Dskeletonization,” IEEE Trans. Med. Imaging, Vol. 21, pp. 263–273, 2002.
 [27] J.T. Ferrucci, “Colon cancer screening with virtual colonoscopy: promise, polyps, politics,” Am. J. Roentgenol., Vol. 177, No. 5, pp. 975–988, 2001.
 [28] D. Ravanelli, E. dal Piaz, M. Centonze, G. Casagranda, M Marini, M. Del Greco, R. Karim, K. Rhode, A. Valentini, “A novel skeleton based quantification and 3D volumetric visualization of left atrium fibrosis using Late Gadolinium Enhancement Magnetic Resonance Imaging,” IEEE Trans. Med. Imaging, 2014, in press.
 [29] B.R. Gomberg, P.K. Saha, H.K. Song, S.N. Hwang, F.W. Wehrli, “Topological analysis of trabecular bone MR images,” IEEE Trans. Med. Imaging, Vol. 19, No. 3, pp. 166–174, 2000.
 [30] F.W. Wehrli, B.R. Gomberg, P.K. Saha, H.K. Song, S.N. Hwang, P.J. Snyder, “Digital Topological Analysis of In Vivo Magnetic Resonance Microimages of Trabecular Bone Reveals Structural Implications of Osteoporosis,” J. Bone Miner. Res., Vol. 16, No. 8, pp. 1520–1531, 2001.
 [31] P.K. Saha, Y. Xu, H. Duan, A. Heiner, G. Liang, “Volumetric Topological Analysis: A Novel Approach for Trabecular Bone Classification on the Continuum Between Plates and Rods,” IEEE Trans. Med. Imaging, Vol. 29, No. 11, pp. 1821–1838, 2010.
 [32] M. Couprie, G. Bertrand, “New Characterizations of Simple Points in 2D, 3D, and 4D Discrete Spaces”, IEEE Trans. Pattern Anal. Mach. Intell., Vol. 31, No. 4, pp. 637–648, 2009.
 [33] A. Hatcher, “Algebraic Topology”, Cambridge University Press (Available Online), 2002.
 [34] M.M. Cohen, “A course in simple homotopy theory,” Graduate texts in mathematics vol. 10, SpringerVerlag, New York, 1973.
 [35] M. Niethammer, W.D. Kalies, K. Mischaikow, A. Tannenbaum, “On the detection of simple points in higher dimensions using cubical homology,” IEEE Trans. Image Process., Vol. 15, No. 8, pp. 2462–2469, 2006.
 [36] P.J. Frey, P.L. George, “Mesh generation: application to finite elements,” Hermes Science Publishing, Oxford, UK, 2000.
 [37] T.Y. Kong, A. Rosenfeld, “Digital topology, Introduction and survey”, Comput. Vision Graphics Image Process, Vol. 46, pp. 357–393, 1989.
 [38] T.Y. Zhang, C.Y. Suen, “A fast parallel algorithm for thinning digital patterns,” Commun. ACM, Vol. 27, pp. 236239, 1984.
 [39] C.M. Holt, A. Stewart, M. Clint, R.H. Perrott, “An improved parallel thinning algorithm,” Commun. ACM, Vol. 30, No. 2, pp. 156–160, 1987.
 [40] B.K. Jang, R.T. Chin, “Analysis of Thinning Algorithms Using Mathematical Morphology,” IEEE Trans. Pattern Anal. Mach. Intell., Vol. 12, No. 6, pp. 54–551, 1990.
 [41] M. Ahmed, R. Ward, “A rotation invariant rulebased thinning algorithm for character recognition,” IEEE Trans. Pattern Anal. Mach. Intell., Vol. 24, No. 12, pp. 1672–1678, 2002.
 [42] P.I. Rockett, “An improved rotationinvariant thinning algorithm,” IEEE Trans. Pattern Anal. Mach. Intell., Vol. 27, No. 10, pp. 16711674, 2005.
 [43] M. Ashwin, G. Dinesh, “A new sequential thinning algorithm to preserve topology and geometry of the image,” International Journal of Mathematics Trends and Technology, Volume 2. Issue 2, pp. 1–5, 2011.
 [44] S. Lobregt, P.W. Verbeek, F.C.A. Groen, “ThreeDimensional Skeletonization: Principle and Algorithm,” IEEE Trans. Pattern Anal. Mach. Intell., Vol. 2, No. 1, pp. 75–77, 1980.
 [45] T.C. Lee, R.L. Kashyap, “Building skeleton models via 3D medial surface/axis thinning algorithm,” CVGIP: Graphical models and image processing, Vol. 56, No. 6, pp. 462–478, 1994.
 [46] P.P. Jonker, O. Vermeij, ”On skeletonization in 4d images,” Proceeding SSPR ’96 Proceedings of the 6th International Workshop on Advances in Structural and Syntactical Pattern Recognition, SpringerVerlag London, UK 1996.
 [47] C. Pudney, “DistanceOrdered Homotopic Thinning: A Skeletonization Algorithm for 3D Digital Images,” Computer Vision and Image Understanding, Vol. 72, No. 3, pp. 404–413, 1998.
 [48] G. Borgefors, I. Nyström, G. Sanniti Di Baja, “Computing skeletons in three dimensions,” Pattern Recognition, Vol. 32, No. 7, pp. 1225–1236, 1999.
 [49] P.K Saha, B.B. Chaudhuri, D.D. Majumder, “A new shape preserving parallel thinning algorithm for 3D digital images,” Pattern Recognition, Vol. 30, No. 12, pp. 1939–1955, 1997.
 [50] K. Palágyi, A. Kuba, “A Parallel 3D 12Subiteration Thinning Algorithm,” Graphical Models and Image Processing, Vol. 61, No. 4, pp. 199–221, 1999.
 [51] K. Palágyi, E. Balogh, A. Kuba, C. Halmai, B. Erdöhelyi, E. Sorantin, K. Hausegger, “A Sequential 3D Thinning Algorithm and Its Medical Applications,” Information Processing in Medical Imaging, Lecture Notes in Computer Science, Vol. 2082, pp. 409–415, 2001.
 [52] K. Palágyi, “A 3subiteration 3D thinning algorithm for extracting medial surfaces,” Pattern Recognition Letters, Vol. 23, No. 6, pp. 663–675, 2002.
 [53] W. Xie, R.P. Thompson, R. Perucchio, “A topologypreserving parallel 3D thinning algorithm for extracting the curve skeleton,” Pattern Recognition, Vol. 36, No. 7, pp. 1529–1544, 2003.
 [54] T.A. Aberra, “Topology Preserving Skeletonization of 2D and 3D Binary Images,” Master Thesis, Department of Mathematics, Technische Universitat Kaiserlautern, 2004.
 [55] T. Ju, M.L. Baker, W. Chiu, “Computing a family of skeletons of volumetric models for shape description,” ComputerAided Design, Vol. 39, pp. 352360, 2007.
 [56] C. Lohou, J. Dehos, “Automatic Correction of Ma and Sonka’s Thinning Algorithm Using PSimple Points,” IEEE Trans. Pattern Anal. Mach. Intell., Vol. 32, No. 6, pp. 1148–1152, 2010.
 [57] G. Németh , P. Kardos , K. Palágyi, “Thinning combined with iterationbyiteration smoothing for 3D binary images,” Graphical Models, Vol. 73, No. 6, pp. 335–345, 2011.
 [58] P. Wiederhold, S. Morales, “Thinning on Quadratic, Triangular, and Hexagonal Cell Complexes,” Combinatorial Image Analysis, Vol. 4958, Lecture Notes in Computer Science, pp. 13–25, 2008.
 [59] P. Wiederhold, S. Morales, “Thinning on cell complexes from polygonal tilings, Discrete Applied Mathematics,” Vol. 157, No. 16, pp. 3424–3434, 2009.
 [60] P. Kardos, K. Palgyi, “On Topology Preservation for Triangular Thinning Algorithms,” Combinatorial Image Analaysis, Lecture Notes in Computer Science, Vol. 7655, pp. 128–142, 2012.
 [61] P. Kardos, K. Palgyi, “Topologypreserving hexagonal thinning,” International Journal of Computer Mathematics, Vol. 90, No. 8, pp. 1607–1617, 2014.
 [62] T. Dey, J. Sun, “Defining and computing curveskeletons with medial geodesic function”, Proc. of the fourth Eurographics symposium on Geometry, pp. 143–152, 2006.
 [63] L. Tcherniavski, P. Stelldinger, “A Thinning Algorithm for Topologically Correct 3D Surface Reconstruction,” VIIP 2008, Proc. 8th IASTED International Conference on Visualization, Imaging, and Image Processing, J.J. Villanueva (Ed.), pp. 119–124, ACTA Press, 2008.
 [64] J. Cao, A. Tagliasacchi, M. Olson, H. Zhang, Z. Su, “Point Cloud Skeletons via Laplacian Based Contraction,” Shape Modeling International Conference (SMI), pp.187–197, 2123 June 2010.
 [65] L. Liu, E.W. Chambers, D. Letscher, T. Ju, “A simple and robust thinning algorithm on cell complexes,” Comp. Graph. Forum, Vol 29, No. 7, pp. 2253–2260, 2010.
 [66] P. Dłotko, R. Specogna, “Thin It code,” http://www.thinit.org, http://www.diegm.uniud.it/elettrotecnica/web/thinit/.
 [67] A.A. Markov, “Insolubility of the Problem of Homeomorphy,” (English translation available online on A. Zomorodian webpage), 1958.
 [68] P. Dłotko, R. Specogna, “Physics inspired algorithms for (co)homology computations of threedimensional combinatorial manifolds with boundary,” Comput. Phys. Commun., Vol. 184, No. 10, pp. 2257–2266, 2013.
 [69] P. Dłotko, R. Specogna, “Cohomology in 3d magnetoquasistatics modeling,” Commun. Comput. Phys., Vol. 14, No. 1, pp. 48–76, 2013.
 [70] Computer Assisted Proofs in Dynamic, www.redhom.ii.uj.edu.pl.
 [71] M. Mrozek, P. Pilarczyk, N. elazna, “Homology Algorithm Based on Acyclic Subspace,” Computers and Mathematics with Applications, Vol. 55, pp. 2395–2412, 2008.
 [72] P. Brendel, P Dłotko, M Mrozek, N. elazna, “Homology Computations via Acyclic Subspace,” Computational Topology in Image Context, LNCS 7309, pp. 117–127, 2012.
 [73] R.H. Bing, “Some aspects of the topology of 3manifolds related to the Poincaré Conjecture,” Lectures on Modern Mathematics II, T.L. Saaty ed., Wiley, pp. 93–128, 1964.
 [74] E. Perilli, F. Baruffaldi, “Proposal for shared collections of Xray microCT datasets of bone specimens,” ICCB03, 2426 September 2003, Zaragoza, Spain. The data are available at \(http://www.biomedtown.org\).
Comments
There are no comments yet.