Aldo Gonzalez-Lorenzo

Alexandra Bac

Jean-Luc Mari

Pedro Real

Cellular Skeletons: a New Approach to Topological Skeletons with Geometric Features

This paper introduces a new kind of skeleton for binary volumes called the cellular skeleton. This skeleton is not a subset of voxels of a volume nor a subcomplex of a cubical complex: it is a chain complex together with a reduction from the original complex. Starting from the binary volume we build a cubical complex which represents it regarding 6 or 26-connectivity. Then the complex is thinned using the proposed method based on elementary collapses, which preserves significant geometric features. The final step reduces the number of cells using Discrete Morse Theory. The resulting skeleton is a reduction which preserves the homology of the original complex and the geometrical information of the output of the previous step. The result of this method, besides its skeletonization content, can be used for computing the homology of the original complex, which usually provides well shaped homology generators.

Introduction

The notion of skeleton (or medial axis) was introduced by H. Blum in 1967 [START_REF] Blum | A Transformation for Extracting New Descriptors of Shape[END_REF]. Given a subset S ⊂ R n , its medial axis is the set of all the points in S that are centres of maximal balls included in S. It is a good descriptor of shape, it is thin, it has the same type of homotopy as S (one can be continuously deformed into the other) and it is reversible (we can reconstruct S by using the distance of each point to the boundary).

Skeletons are widely used for various applications such as video tracking [GSdA + 09], shape recognition [START_REF] Yu | Skeleton-based recognition of chinese calligraphic character image[END_REF], surface sketching [START_REF] Mari | Surface sketching with a voxel-based skeleton[END_REF] and in many other scientific domains.

In the discrete context there is no unique equivalent to the continuous definition. A discrete medial axis gives a skeleton which is not homotopically equivalent to the discrete object. Discrete skeletons are usually based on the thinning of the object, using mathematical morphology or parallel approaches: simple points [START_REF] Kong | Digital topology: Introduction and survey[END_REF], simple P-points [START_REF] Lohou | Two symmetrical thinning algorithms for 3D binary images, based on P-simple points[END_REF], critical kernels [START_REF] Couprie | Isthmus-based parallel and asymmetric 3D thinning algorithms[END_REF]. The given skeleton is a subset of voxels, much smaller, which is homotopically equivalent and should preserve in an uncertain way the geometrical features of the object. There is another class of skeletons which first computes a cubical complex associated to an object and then selects a subcomplex (see [START_REF] Chaussard | Surface thinning in 3D cubical complexes[END_REF][START_REF] Liu | A simple and robust thinning algorithm on cell complexes[END_REF][START_REF] Couprie | Hierarchic Euclidean skeletons in cubical complexes[END_REF][START_REF] Couprie | Topological maps and robust hierarchical Euclidean skeletons in cubical complexes[END_REF][START_REF] Dlotko | Topology preserving thinning of cell complexes[END_REF]).

We must note that geometrical preservation is more intuitive or heuristic in the context of discrete skeletons than in the continuous one. Hence, all previous approaches exist, with their respective strengths and weaknesses.

This paper presents a three-step method:

1. Starting from a binary volume (a set of voxels in a regular grid), we define its associated cubical complex. We propose two different constructions which encode 6 and 26-connectivity. These associated cubical complexes were separately defined in [CC09, LCLJ10] 2. We compute a skeleton of this cubical complex. We propose a method strongly based on [START_REF] Liu | A simple and robust thinning algorithm on cell complexes[END_REF] which performs elementary collapses which are encoded in a discrete gradient vector field (DGVF). It is a simple algorithm producing satisfying results. This step can also be addressed with the algorithms found in [START_REF] Liu | A simple and robust thinning algorithm on cell complexes[END_REF][START_REF] Dlotko | Topology preserving thinning of cell complexes[END_REF] or in [START_REF] Chaussard | Surface thinning in 3D cubical complexes[END_REF][START_REF] Couprie | Hierarchic Euclidean skeletons in cubical complexes[END_REF][START_REF] Couprie | Topological maps and robust hierarchical Euclidean skeletons in cubical complexes[END_REF] given a fixed parameter. As a result, the shape of the skeleton is defined 3. This step is completely new. Given the reduced cubical complex, we extend the previous DGVF in order to obtain a reduction (see [START_REF] Sergeraert | Effective homology, a survey[END_REF]) between the chain complex of the original cubical complex and a reduced one, with the property of maintaining the shape of the skeleton computed in the second step.

Our approach presents several advantages: the topology preservation throughout the thinning is guaranteed by the Discrete Morse Theory; our method dissociates connectivity and skeleton extraction, thus the thinning algorithm becomes independent of the connectivity relation; the obtained reduction accelerates the computation of its homology since the reduced complex contains fewer cells. Moreover, starting the homology computation from this cellular skeleton should produce well shaped homology generators, as they are included in the skeleton. This reveals an advantage of computing homology using Discrete Morse Theory: we can control the shape of the homology generators.

In Sect. 2, we introduce all the necessary definitions for understanding our method. In Sect. 3, we describe our approach. Section 4 shows some results of our framework on some binary volumes. We finish this paper by presenting our conclusion and our future perspectives.

Preliminaries

Binary Volumes and Cubical Complexes

A 3D binary volume is a set of voxels centred on integer coordinates. We will describe it by the set of the coordinates of its elements.

The rest of this section is derived from [START_REF] Kaczynski | [END_REF]. For a deeper understanding of these concepts, the reader can refer to it. An elementary interval is an interval of the form [k, k+1] (nondegenerate) or a set {k} (degenerate), also denoted as [k, k], where k ∈ Z. An elementary cube in R n is the Cartesian product of n elementary intervals, and the number of nondegenerate intervals in this product is its dimension. An elementary cube of dimension q will be called q-cube for short, or even q-cell, since cubical complexes are a special kind of cell complexes.

Given two elementary cubes P and Q, we say that P is a face of Q if P ⊂ Q and we note it P < Q. It is a primary face if the difference of their dimensions is 1. The definitions of the dual concepts coface and primary coface are immediate.

A cubical complex is a set of elementary cubes with all of their faces. The boundary of an elementary cube is the collection of its primary faces.

A cubical complex can be completely described by its Hasse diagram. It is a directed graph whose vertices are all the elementary cubes, and whose arrows go from each cube to its primary faces. In this paper we will usually not make the distinction between the vertices and the elementary cubes they represent, so we will mix the terms vertex, cube and cell.

Chain Complexes

A chain complex (C * , d * ) is a sequence of R-modules C 0 , C 1 , . . . (called chain groups) and homomorphisms d 1 : C 1 → C 0 , d 2 : C 2 → C 1 , . . . (called differential
or boundary operators) such that d q-1 d q = 0, ∀q > 0, where R is some ring, called the ground ring or ring of coefficients.

Given a cubical complex, we define its chain complex (with coefficients in Z 2 ) as follows:

-C q is the free group generated by the q-cubes of the complex. Their elements are called q-chains. d q gives the "algebraic" boundary, which is the sum of the primary faces of the q-cubes.

Discrete Morse Theory

Discrete Morse Theory was introduced by Robin Forman as a discretization of the Morse Theory [START_REF] Forman | A user's guide to discrete Morse theory[END_REF]. The main idea is to obtain some homological information by means of a function defined on the complex. This function is equivalent to a discrete gradient vector field and we will rather use this notion.

A discrete vector field on a cubical complex is a matching on its Hasse diagram, that is a collection of edges such that no two of them have a vertex in common. From a Hasse diagram and a discrete vector field we can define a Morse graph: it is a graph similar to the Hasse diagram except for the arrows contained in the matching, which are reverted. These arrows will be called integral arrows, and the others, differential arrows.

A V-path is a path on the Morse graph which alternates between integral and differential arrows. A discrete gradient vector field is a discrete vector field which does not contain any closed V-path. A critical vertex (or critical cell) is a vertex which is not paired by the matching. Figure 1 shows a cubical complex with a DGVF.

A DGVF can be given by a set of elementary collapses. An elementary collapse [Whi50] consist of removing a free pair from a cell complex, that is a cell Fig. 1: A DGVF over a cubical complex. Red segments represent integral arrows going from one cell to some of its primary cofaces. The critical cells are represented in blue.

with a primary face which does not have any other coface. A collapse is a sequence of elementary collapses. The homotopy type of a complex is invariant under collapses. The free pairs of a collapse define a DGVF.

Reduction

The Effective Homology theory [START_REF] Sergeraert | Effective homology, a survey[END_REF] provides a tool that establishes a strong relation between two chain complexes, called reduction. Formally, a reduction between two chain complexes (C * , d * ) and (C * , d * ) is a triple of homomorphisms (h * , f * , g * ) such that:

-h q : C q → C q+1 for every q ≥ 0 -f q : C q → C q is a chain map (f d = d f ) -g q : C q → C q is also a chain map (gd = dg) -gf = 1 -dh -hd -f g = 1 C -hh = hf = hg = 0
This notion is a special case of chain contraction [START_REF] Eilenberg | On the groups h(π, n), i[END_REF] or strong deformation retraction data [START_REF] Lambe | Applications of perturbation theory to iterated fibrations[END_REF]. It is a usual tool for reducing chain complexes in order to compute their homology. We will use it to define the Morse complex of a cubical complex endowed with a DGVF. The exact definition is given in [START_REF] Gonzalez-Lorenzo | Computing homological information based on directed graphs within discrete objects[END_REF]. Roughly speaking, the Morse complex is a cell complex composed of the critical cells, which is homotopically equivalent to the original cubical complex.

Our Framework

Overview

In this section we present the structure of our framework. It should be noted that, although the algorithms are designed for 3D volumes, they can be generalized to any dimension: one only needs to replace parameter 3 by n.

Let us briefly recall the structure of our approach, previously described in the introduction. We start with a binary 3D volume. Depending of the connectivity relation which we want to use (6 or 26), we build its associated cubical complex with one of the methods explained in Sect. 3.2. Next, we make a homotopic thinning of the cubical complex based on elementary collapses encoded in a DGVF. Finally we minimize the number of cells in the Morse complex by preserving its shape in Algorithm 4.

The final result is a cellular skeleton: a reduction from the original complex to a reduced one. Then we give a representation of this skeleton by showing the cells of the chain f (σ) (see Sect. 2.4) for every critical cell.

Construction of the Cubical Complex: Choosing the Connectivity

As explained earlier, the first step of our approach consists in building the cubical complex associated to the digital volume and the connectivity relation chosen.

Hence, there are two encodings: one for 6-connectivity (that is 2n-connectivity in dimension n) and other for 26-connectivity ((3 n -1)-connectivity). Figure 2 illustrates these two cubical complexes associated to the same binary volume.

Fig. 2: Left: a binary volume. Center: its primal associated cubical complex. Right: its dual cubical complex

The primal associated cubical complex. We encode a binary volume equipped with 26-connectivity into a cubical complex (called primal associated cubical complex ). In this case, the construction is quite elementary as every voxel x = (x 1 , x 2 , x 3 ) generates the 3-cube [

x 1 , x 1 + 1] × [x 2 , x 2 + 1] × [x 3 , x 3 + 1]
and all its faces. This method was already presented in [START_REF] Chaussard | Surface thinning in 3D cubical complexes[END_REF].

The dual associated cubical complex. Another approach consists in encoding a binary volume equipped with 6-connectivity into a cubical complex (called dual associated cubical complex ). Let us first adapt the notion of clique to our context: a d-clique is a maximal (in the sense of inclusion) set of voxels such that their intersection is a d-cube. First, for every voxel (in fact 3-clique) x = (x 1 , x 2 , x 3 ) of the volume, we add the 0-cube σ

= [x 1 ] × [x 2 ] × [x 3 ].
Then, for every d-clique (d < 3) in the volume, we add to the cubical complex a (3 -d)cube such that its vertices are the voxels of the d-clique. This approach was used in [START_REF] Liu | A simple and robust thinning algorithm on cell complexes[END_REF].

Homotopic Thinning Algorithm

This step performs a homotopic thinning of the cubical complex. This is done by establishing a DGVF, which can be seen as a set of elementary collapses (deletion of free pairs). Actually, this DGVF describes the relation between the original complex and the thinned one (the Morse complex) in terms of a reduction. Let us point out that our approach does not have to deal with simple points, critical kernels, etc.

A simple thinning algorithm with satisfying results was given in [START_REF] Liu | A simple and robust thinning algorithm on cell complexes[END_REF]. The algorithm is described in three steps:

Step 1: Thinning Perform an iterative thinning: at each iteration, all free pairs are identified and then collapsed while it is possible. For every cell σ, I(σ) is the first iteration after which σ has no cofaces and R(σ) is the iteration in which σ is removed. If I(σ) is defined,

M abs (σ) = R(σ) -I(σ) and M rel (σ) = 1 - I(σ) R(σ)
Step 2: Clustering Given some thresholds ε q abs , ε q rel and τ q (k = 1, 2), consider the set B of the cells scoring higher than ε q abs and ε q rel . Remove from this set those cells whose connected component size is fewer than τ q .

Step 3: Thinning Repeat the first step while maintaining the cells in B.

Note that when we identify all the free pairs, we can find several pairs (τ > σ 1 ), (τ > σ 2 ), . . . and we must choose one of them. This choice was defined as arbitrary in [START_REF] Liu | A simple and robust thinning algorithm on cell complexes[END_REF] and it was pointed out that this should be studied. Our contribution in this step is an alternative to the simple iterative collapse which makes some of these collapses order-independent. It is described in detail in Algorithm 2, which calls Algorithm 1. Indeed, for each elementary collapse between a maximal cell τ and one of its primary faces σ 1 , . . . , σ n , we must choose one of them. We partially solve this problem by performing multiple elementary collapses (see Algorithm 1) in the same iteration in order to remove some of these cells. Sometimes the choice of the cell σ 1 is irrelevant, as shown in Fig. 3. This seems to be related to the notion of simple cell [DS14, Def. 3.3].

Algorithm 1 AdvancedCollapse

Require: K a cubical complex; [τ, σ1, . . . , σn] cells such that (τ > σi) are free pairs; it the iteration; R, I two maps over K; V a DGVF over K. Ensure: R , I , V extensions of R, I, V .

1: V ← V, R ← R, I ← I, Q an empty queue 2: V ← (σ1, τ ); R (σ1), R (τ ) ← it 3: for all σi, i = 2, . . . , n do 4: I (σi) ← it 5: for all ρ < σ critical primary face do 6:

Q ← ρ 7: while Q not empty do 8:

ρ ← Q, C ← set of its critical primary cofaces 9:

if C is empty then 10:

Make an elementary collapse if there exists π < ρ critical primary face with only one critical coface; update V , R and I 11:

if C = {σi} then 12: V ← (ρ, σi); R (ρ), R (σi) ← it 13:
for all ρ < σi critical primary face do 14:

Q ← ρ return (R , I , V )
Fig. 3: Up: the choice of the first elementary collapse does not affect the collapse. Down: only the two first (left and center) collapses are independent of the first elementary collapse. In such a situation, a choice is necessary.

Cell Clustering: Minimizing the Number of Cells

After the previous step, the shape of the skeleton is already defined and we have reduced the number of cubes. Nevertheless, this can be improved. Algorithm 4, which calls Algorithm 3, describes this step.

This step, which is the main novelty of this article, has the following property. Let V be the DGVF computed at step 2 and V its extension returned at the end Algorithm 2 Alternative to Step 1: Thinning Require: K a cubical complex. Ensure: R, I two maps K; V a DGVF over K.

1: it ← 0 2: M ← set of the maximal critical cells (without critical cofaces) 3: for all σ ∈ M do 4: I(σ) ← it, if it is not defined 5: repeat 6:

it ← it + 1 7:

for all τ ∈ M of dimension > 0 do 8:

(τ, σ1), . . . , (τ, σn) ← the free pairs containing τ 9:

(R, I, V) ← AdvancedCollapse(K, [τ, σ1, . . . , σn] , it, R, I, V) 10:

M ← set of the maximal critical cells 11:

for all σ ∈ M do 12: I(σ) ← it, if it is not defined 13: until idempotency 14: return (R, I, V) Algorithm 3 BlockedCollapse Require: K a cubical complex, V a DGVF over K, B a set of (blocked) cells Ensure:

V an extension of V 1: V ← V 2: FreePairs ← ∅ 3: for all critical cell σ ∈ K not in B do 4:
if σ has only one critical primary coface τ not in B then 5:

FreePairs ← (σ, τ ) 6: for all (σ, τ ) in FreePairs do 7:

if τ is critical then 8:

V ← (σ, τ ) return V of the step 3: for each critical cell σ of V, there exists one and only one critical cell σ of V such that σ appears in the chain f (σ).

A direct consequence of this property is that when we visualize the chains f (Cr) (Cr the set of critical cells in V ), we obtain the skeleton computed at step 2. Hence, this DGVF preserves the geometric structure of V .

Algorithm 4 has two motivations:

-Following the previous property, by displaying the chains f (Cr) we obtain a cell complex homotopically equivalent to the initial cubical complex, in which the q-cells are unions of q-cubes. This can be considered as a classification of the skeleton in terms of manifolds. -This step accelerates a later homology computation based on the reduction induced by the resulting DGVF. Moreover, as homology generators are contained in the skeleton, they are supposed to be well shaped, according to the geometric properties of the skeleton.

Algorithm 4 CellClustering

Require: K a cubical complex endowed with a DGVF V Ensure: V an extension of V. 1: V ← V 2: FinalCells ← ∅ 3: for q = 2, 1 do 4: BlockedCells ← ∅ 5:

for all critical (q -1)-cell σ ∈ K do 6:

if σ has = 2 critical primary cofaces then 7:

BlockedCells ← σ 8: repeat 9:

Take any critical q-cell γ not in FinalCells 10:

FinalCells ← γ 11: repeat 12:

V ← BlockedCollapse(K, V , BlockedCells ∪ FinalCells) 13:

until idempotency 14:

until idempotency return V

Validation and Discussion

Our algorithms have been implemented in C++ using the library DGtal [DGt].

In the following we only consider the dual associated cubical complex. Also, as in [START_REF] Liu | A simple and robust thinning algorithm on cell complexes[END_REF] we will consider the thresholds ε q abs = 0.05 • L, ε q rel = 0.05 and τ q = (0.05 • L) k , where L is the width of the bounding box.

Seeing points or edges belonging to big cellular skeletons can be difficult. Hence, we propose a voxelized version of the cellular skeleton. Its construction is quite simple: every

cubical cell [a 1 , a 2 ] × [b 1 , b 2 ] × [c 1 , c 2 ]
produces all the possible voxels (a 1 , b 1 , c 1 ), (a 2 , b 1 , c 1 ), (a 1 , b 2 , c 1 ), (a 2 , b 2 , c 1 ), (a 1 , b 1 , c 2 ), . . .. These voxels are coloured following the dimension of the cubical cell which created them: red for 2-cubes, green for 1-cubes and blue for 0-cubes.

Figure 4 compares some cellular skeletons obtained with our algorithm and [LCLJ10] over some voxels sets. Let us remark that both give similar results, even on the third example, which does not keep all the fingers. There is the need to study further how different fixed orders in the elementary collapses affect both algorithms. Also, we observe that upper skeletons are thinner. Hence, we intend to study whether we can obtain even more similar skeletons by adjusting the thresholds.

Conclusion and Future Perspectives

The present paper introduces a new kind of skeleton for binary volumes which is a chain complex together with a reduction, which is obtained by a three-step method. It works for different connectivity relations and it does not make use of look-up tables. Our main contribution is the third step, which is completely new and provides an alternative skeleton representation, synthesizing further homology computations and geometric representation (each cell in the cellular skeleton actually stands for a piece of manifold in the geometric skeleton).

Our future goals are:

to better understand the difference between choosing the primal or the dual associated cubical complex. The primal associated cubical complex contains more cells, and this can affect the complexity of our algorithms; to study different thresholds for the second step. Since it differs from [START_REF] Liu | A simple and robust thinning algorithm on cell complexes[END_REF],

we need to study further appropriate thresholds; to estimate the complexity of our method.

Fig. 4 :

 4 Fig. 4: Top: cellular skeletons computed using the proposed method. Bottom: cellular skeletons computed with the algorithm in [LCLJ10].

Acknowledgement

The hand and the statue datasets are courtesy of Michel Couprie, and the Happy Buddha is courtesy of the Stanford University Computer Graphics Laboratory.