John Chaussard

Michel Couprie

Surface thinning in 3D cubical complexes

We introduce a parallel thinning algorithm with directional substeps based on the collapse operation, which is guaranteed to preserve topology and to provide a thin result. Then, we propose two variants of a surface-preserving thinning scheme, based on this parallel directional thinning algorithm. Finally, we propose a methodology to produce filtered surface skeletons, based on the above thinning methods and the recently introduced discrete λ-medial axis.

1 To be more precise, we say that a transformation Ψ "preserves topology" if X is homotopy-equivalent to Ψ (X) for any X.

Introduction

The notion of skeleton plays a major role in shape analysis and recognition. Informally, the skeleton of an object X is a subset of X which is 1) thin, 2) centered in X and 3) topologically equivalent 1 to X. To extract skeletons from geometric objects, different methods have been proposed, relying on different frameworks: discrete geometry [START_REF] Borgefors | The Euclidean distance transform: finding the local maxima and reconstructing the shape[END_REF][START_REF] Ge | On the generation of skeletons from discrete Euclidean distance maps[END_REF][START_REF] Malandain | Euclidean skeletons[END_REF][START_REF] Rémy | Exact medial axis with Euclidean distance[END_REF][START_REF] Wim | Euclidean skeletons of digital image and volume data in linear time by the integer medial axis transform[END_REF], digital topology [START_REF] Davies | Thinning algorithms: a critique and a new methodology[END_REF][START_REF] Vincent | Efficient computation of various types of skeletons[END_REF][START_REF] Talbot | Euclidean skeletons and conditional bisectors[END_REF][START_REF] Pudney | Distance-ordered homotopic thinning: a skeletonization algorithm for 3D digital images[END_REF], mathematical morphology [START_REF] Serra | Image analysis and mathematical morphology[END_REF][START_REF] Soille | Morphological image analysis[END_REF], computational geometry [START_REF] Attali | Delaunay conforming iso-surface, skeleton extraction and noise removal[END_REF][START_REF] Attali | Modelling noise for a better simplification of skeletons[END_REF][START_REF] Ogniewicz | Hierarchic Voronoi skeletons[END_REF], and partial differential equations [START_REF] Siddiqi | The Hamilton-Jacobi skeleton[END_REF].

In this paper, we focus on skeletons in the discrete 3-dimensional cubical space. In such spaces, topology preservation is usually guaranteed by the use of topology-preserving thinning methods, based e.g. on simple point deletion [START_REF] Kong | Digital topology: introduction and survey[END_REF][START_REF] Couprie | New characterizations of simple points in 2D, 3D and 4D discrete spaces[END_REF]. The centering of the skeleton may be achieved to some extent by the use of parallel thinning methods [START_REF] Bertrand | A new 3D parallel thinning scheme based on critical kernels[END_REF] or by imposing geometrical constraints such as the preservation of the centers of maximal included balls [START_REF] Davies | Thinning algorithms: a critique and a new methodology[END_REF]. However it is generally difficult, in discrete spaces, to satisfy conditions 1), 2) and 3) together. In particular, the methods evoked above fail to guarantee a thin result.

In order to overcome these limitations, we adopt in this article the framework of cubical complexes. Abstract (cubical) complexes have been promoted in particular by V. Kovalevsky [START_REF] Kovalevsky | Finite topology as applied to image analysis[END_REF] in order to provide a sound topological basis for image analysis. Intuitively, a cubical complex may be thought of as a set of elements having various dimensions (e.g. cubes, squares, edges, vertices) glued together according to certain rules (see Fig. 3). In this framework, we can say that a complex is "thin" if it does not contain any 3-dimensional element (cube). The thinning methods that will be proposed in the sequel of this paper produce thin results in this sense.

In the framework of cubical complexes, the notion of critical kernel introduced by G. Bertrand constitutes a powerful tool to study parallel homotopic thinning in any dimension, which unifies and encompasses previous works on parallel thinning [START_REF] Bertrand | On critical kernels[END_REF][START_REF] Bertrand | On parallel thinning algorithms: minimal nonsimple sets, P-simple points and critical kernels[END_REF]. Indeed, the very notion of critical kernel may be seen as thinning scheme, which consists of iteratively computing the critical kernel of the result of the previous step. Critical kernels may also be used to design new algorithms, as well as to check the topological validity of existing ones [START_REF] Bertrand | Two-dimensional parallel thinning algorithms based on critical kernels[END_REF][START_REF] Bertrand | A new 3D parallel thinning scheme based on critical kernels[END_REF]. However, thinning algorithms based on critical kernels do not guarantee that the result has a dimension strictly lower than the input. Early work on thinning in cubical complexes with the aim of reducing dimension can be found in [START_REF] Liu | 3d thinning on cell complexes for computing curve and surface skeletons[END_REF].

In our work, topology preservation will be ensured by the use of the collapse operation, which is an elementary topology-preserving transformation that has been introduced by J.H.C. Whitehead [START_REF] Whitehead | Simplicial spaces, nuclei and m-groups[END_REF] and plays an important role in combinatorial topology. It can be seen as a discrete analogue of a retraction, that is, a continuous deformation of an object onto itself.

The contributions of this article are the following. We introduce a parallel thinning algorithm with directional substeps based on collapse, which is guaranteed to preserve topology and to provide a thin result. Then, we propose two variants of a surface-preserving thinning scheme, based on this parallel directional thinning algorithm. Finally, we propose a methodology to produce filtered surface skeletons, based on the above thinning methods and the recently introduced discrete λ-medial axis [START_REF] Chaussard | A discrete lambda-medial axis[END_REF].

Basic notions

Cubical complexes

Let Z be the set of integers, we consider the family of sets F 1 0 and F 1 1 , such that

F 1 0 = {{a} | a ∈ Z} and F 1 1 = {{a, a + 1} | a ∈ Z}.
Any subset f of Z n such that f is the cartesian product of m elements of F 1 1 and (nm) elements of F 1 0 is called a face or an m-face of Z n , m is the dimension of f , we write dim(f ) = m. A 0-face is called a vertex, a 1-face is an edge, a 2-face is a square, and a 3-face is a cube (see Fig. 1). We denote by F n the set composed of all faces in

Z n . Let f ∈ F n . We set f = {g ∈ F n |g ⊆ f }, and f * = f \ {f }. Any element of f
is a face of f , and any element of f * is a proper face of f . We call star of f the set f = {g ∈ F n |f ⊆ g}, and we write f * = f \ {f }: any element of f is a coface of f . It is plain that g ∈ f iff f ∈ ǧ.

A set X of faces in F n is a cell, or m-cell, if there exists an m-face f ∈ X such that X = f . The closure of a set of faces X is the set

X -= ∪{ f |f ∈ X}. The set X is F n \ X.
A finite set X of faces in F n is a cubical complex if X = X -, and we write X F n . Any subset Y of X which is also a complex is a subcomplex of X, and we write Y X. A face f ∈ X is a facet of X if f is not a proper face of any face of X. We denote by X + the set composed of all facets of X. A complex X is pure if all its facets have the same dimension. The dimension of X is dim(X) = max{dim(f ) | f ∈ X}. If dim(X) = d, then we say that X is a d-complex.

From binary images to complexes

Traditionally, a binary image (see Fig. 2a) is represented as a subset of Z 2 (or Z 3 for 3D images). Let S ⊆ Z n , the elements of S represent the pixels (or voxels) of the image, often called object pixels or black pixels. The set S = Z n \ S represents the background (set of white pixels or voxels). To obtain a complex from a subset of Z n , we consider the following bijection which associates, to each point of Z n , a facet of F n . Let x = (x 1 , . . . , x n ) ∈ Z n , define Φ(x) = {x 1 , x 1 + 1} × . . . × {x n , x n + 1}. See Fig. 2 an illustration in 2D, where the image of each pixel by Φ is a 2-face (a square). For example, the black pixel x in the first row has coordinates (3, 0); the corresponding 2-face is Φ(x) = {3, 4} × {0, 1} = {(3, 0), (3, 1), (4, 0), (4, 1)}. The map Φ is straightforwardly extended to sets: Φ(S) = {Φ(x) | x ∈ S}. Then, the complex that we associate to a set S is Φ(S) -(See Fig. 2b). Notice that the datum of an nD binary image is equivalent to the datum of a pure n-complex in F n .

Collapse

The collapse operation consists of removing two distinct elements (f, g) from a complex X under the condition that g is contained in f and is not contained in any other element of X. This operation may be repeated several times (see Fig. 3). A more precise definition follows. Let X F n , and let f, g be two faces of X. The face g is free for X, and the pair (f, g) is a free pair for X if f is the only face of X such that g is a proper face of f . In other terms, (f, g) is a free pair for X whenever ǧ * ∩ X = {f }. In this case, we say that the complex X \ {f, g} is an elementary collapse of X. It can be easily seen that if (f, g) is a free pair for X and dim(f ) = m, then f is a facet and dim(g) = (m -1).

Let X F n and Y F n be two complexes. We say that X collapses onto Y if there exists a sequence of complexes (X 0 , ..., X ℓ ) of F n such that X = X 0 , Y = X ℓ and for all i ∈ {1, . . . , ℓ}, X i is an elementary collapse of X i-1 (see Fig. 3).

Let f 0 , f ℓ be two n-faces of F n . An (n -1)-path from f 0 to f ℓ is a sequence π = (f 0 , ..., f ℓ ) of faces of F n such that for all i ∈ {0, . . . , ℓ}, either i is even and f i is an n-face, or i is odd and f i is an (n -1)-face with fi

* = {f i-1 , f i+1 }.
The following proposition will serve us to prove the thinness of our skeletons.

Proposition 1. Let X F n be an n-complex, with n > 0. Then X has at least one free (n -1)-face.

Proof. Since X is an n-complex (hence X is finite) there exists an n-face a in X and an n-face b in X. Obviously, there exists an (n -1)-path from a to b. Let h be the first n-face of π that is not in X, let k be the last n-face of π before h (thus k is in X), and let e = k ∩ h be the (n -1)-face of π between k and h.

Since k and h are the only two n-faces of F n that contain e, we see that the pair (k, e) is free for X.

Parallel directional thinning based on collapse

The most "natural" way to thin an object consists of removing some of its border elements in parallel, in a symmetrical manner. However, parallel deletion of free pairs does not, in general, guarantee topology preservation: see for example Fig. 4, where simultaneously removing all free pairs would split the object into separate components. The complex obtained by removing simultaneously from X all the pairs that are free for X.

In the framework of 2D digital topology, a popular method due to A. Rosenfeld [START_REF] Rosenfeld | A characterization of parallel thinning algorithms[END_REF] consists of dividing each thinning step into directional substeps. In each substep, only simple points that have no neighbor belonging to the object in one of the four main directions (north, south, east, west) are candidates for deletion.

However, this method cannot be straightforwardly extended to 3D. In fact, the question of knowing whether this strategy has a "natural extension" to 3D was cited among three open questions relative to digital topology by Kong, Litherland and Rosenfeld in [START_REF] Kong | Problems in the topology of binary digital images[END_REF] (question 547). An answer to this question was recently given in [START_REF] Bertrand | On parallel thinning algorithms: minimal nonsimple sets, P-simple points and critical kernels[END_REF], based on the critical kernels framework. Here, we show that the directional strategy can be directly adapted in the framework of cubical complexes, and provides interesting results.

First, we need to define the direction and the orientation of a free face. Let f ∈ F n , the center of f is the center of mass of the points in f , that is,

c f = 1 |f | a∈f a. The center of f is an element of [ Z 2 ]
n , where Z 2 denotes the set of half integers. Let X F n , let (f, g) be a free pair for X, and let c f and c g be the respective centers of the faces f and g. We denote by V (f, g) the vector (c fc g ) of [ Z 2 ] n . Let (f, g) be a free pair, the vector V (f, g) has only one non-null coordinate. We define Dir (f, g) as the index of the non-null coordinate of V (f, g). Thus, Dir () is a surjective function from F n × F n to {1, . . . , n} such that, for all free pairs (f, g) and (i, j) for X, Dir (f, g) = Dir (i, j) if and only if V (f, g) and V (i, j) are collinear. The number Dir (f, g) is called the direction of the free pair (f, g). The free pair (f, g) has a positive orientation, and we write Orient(f, g) = 1, if the non-null coordinate of V (f, g) is positive; otherwise (f, g) has a negative orientation, and we write Orient(f, g) = 0.

Let us now state an elementary property of collapse, which gives a necessary and sufficient condition under which two collapse operations may be performed in parallel while preserving topology. Proposition 2. Let X F n , and let (f, g) and (k, ℓ) be two distinct free pairs for X. The complex X collapses onto X \ {f, g, k, ℓ} if and only if f = k.

Proof. If f = k, then it is plain that (k, ℓ) is not a free pair for Y = X \ {f, g} as k = f / ∈ Y . Also, (f, g) is not free for X \ {k, ℓ}. If f = k, then we have g = ℓ,
ǧ * ∩ X = {f } (g is free for X) and ľ * ∩ X = {k} (ℓ is free for X). Thus, we have ľ * ∩ Y = {k} as ℓ = g and k = f . Therefore, (k, ℓ) is a free pair for Y .

From Prop. 2, the following corollary is immediate.

Corollary 3. Let X F n , and let (f 1 , g 1 ), . . . , (f m , g m ) be m distinct free pairs for X such that, for all a, b ∈ {1, . . . , m} (with a = b), f a = f b . The complex X collapses onto X \ {f 1 , g 1 , . . . , f m , g m }.

Considering two distinct free pairs (f, g) and (i, j) for X F n such that Dir (f, g) = Dir (i, j) and Orient (f, g) = Orient(i, j), we have f = i. From this observation and Cor. 3, we deduce the following property.

Corollary 4. Let X F n , and let (f 1 , g 1 ), . . . , (f m , g m ) be m distinct free pairs for X having all the same direction and the same orientation. The complex X collapses onto X \ {f 1 , g 1 , . . . , f m , g m }.

We say that a d-face of X is a border face if it contains a free (d -1)-face. Define Border (X) as the set of all border faces of X. We are now ready to introduce our directional thinning algorithm (Alg. 1).

Algorithm 1: ParDirCollapse(X, W, ℓ) Data: A cubical complex X F n , a subset W ⊆ X that contains faces of X which should not be removed, and ℓ ∈ N, the number of layers of free faces which should be removed from X while there exists free pairs for X in X \ W and ℓ > 0

do 1 L = Border (X) -; 2 for t = 1 → n do 3 for s = 0 → 1 do 4 for d = n → 1 do 5 E = {(f, g) free for X | g / ∈ W, f / ∈ W, 6 Dir(f, g) = t, Orient(f, g) = s, dim(f ) = d}; G = {(f, g) ∈ E | f ∈ L and g ∈ L}; 7 X = X \ G; 8 l = l -1; 9 return X; 10
Let us first comment on the notion of "layer" and the role of the sets L and G (lines 2 and 7). Intuitively, we want a single execution of the loop line 1 to only affect facets that were on the border of X at the beginning of the loop, in order to ensure that the thinning will not reduce certain parts more quickly than others due to the scanning of directions. Thus, by line 7 we ensure that the layer which is removed in one iteration of the loop line 1 is included in L, which contains all the border facets of the complex X and all the faces included in those.

Remark that different definitions of the functions Dir () and Orient () could be given, corresponding to different orders in which directions and orientations are scanned in Alg. 1, and yielding different results. However in general some arbitrary choices must be made in order to thin a 3-complex into a 2-complex. The choice of these two functions is the only arbitrary choice to be made, after which the results are uniquely defined.

For any complex X F n , any subset W ⊆ X and for any ℓ ∈ N, it follows from Cor. 4 that X collapses onto ParDirCollapse(X, W, ℓ). Furthermore, if W does not contain any n-face nor any (n -1)-face that is not a facet of X, and if ℓ = +∞, then by Prop. 1 it can easily be deduced that ParDirCollapse(X, W, ℓ) contains no n-face.

Notice that checking whether a face is free or not, is quite easy to implement and can be done in constant time, whatever the dimension of the face, if n is considered as fixed. Indeed, Alg. 1 can be used to thin a cubical complex of any dimension. When a free pair (f, g) is removed from X, it is sufficient to scan the faces of f in order to check the appearance of new free faces, avoiding to search the whole complex for free faces. Thus, Alg. 1 may be easily implemented to run in linear time complexity (proportionally to the number of faces of the complex).

As shown in Fig. 5, when the input complex has a nearly constant thickness, it is possible to obtain a surface skeleton of X by choosing a convenient value of ℓ as the last parameter of Alg. 1. However, in most cases, it is not possible to find a value of ℓ which is satisfying for the whole complex. In the next section, we explain how to use the directional thinning strategy to obtain a surface skeleton from a complex without having to tune any "thickness" parameter. In this section, two algorithms are proposed in order to perform a surfacepreserving thinning of a cubical complex. Both algorithms are based on the same principle: a single layer of faces is removed using Alg. 1, then some 2-faces are detected and kept safe from any future removal. This process is repeated until stability. For Alg. 2 (CollapseSurface), the detected 2-faces are simply those that are not contained in any 3-face of the current complex, i.e., the 2-facets of X.

Algorithm 2: CollapseSurface(X, W ) Data: A three-dimensional cubical complex X F 3 , a subset W ⊆ X that contains faces of X which should not be removed while there exists free pairs for X in X\W do 1 ParDirCollapse(X, W, 1); 6a shows the result of CollapseSurface(X, ∅), where X is the tube shown in Fig. 5a. The resulting complex is a 2-complex containing some "branches" and "surface patches" (as shown on the detailed view) which do not represent significant "surfacic features" of the original object. This motivates the introduction of a variant of this method, which imposes a more restrictive condition to preserve 2-facets from removal.

2 W = W ∪ {f ∈ X + | dim(f ) = 2}; 3 return X; 4 Fig.
We say that a 2-facet f ∈ X + is a 2d-isthmus of X if it contains no free edge, in other words, if each edge in f is included in a facet of X distinct from f . In Alg. 3 (CollapseIsthmus), the 2-facets that are detected and kept safe from further removal are the 2d-isthmuses of X.

Algorithm 3: CollapseIsthmus(X, W )

Data: A three-dimensional cubical complex X F 3 , a subset W ⊆ X that contains faces of X which should not be removed while there exists free pairs for X in X\W do 1 ParDirCollapse(X, W, 1);

2 W = W ∪ {f ∈ X | f is a 2d-isthmus of X}; 3 return X; 4
Notice that, following Prop. 1, the results of Alg. 2 and Alg. 3 are thin (they do not contain any 3-face), provided that the input parameter W is set to ∅. However, thanks to parameter W , it is possible to constrain these algorithms to preserve selected parts of the original object, in addition to the surface parts that are automatically detected. Fig. 6b shows the result of CollapseIsthmus(X, ∅), where X is the skewed tube shown in Fig. 5a. The resulting complex contains less branches and spurious surfaces (as shown on the detailed view) than the complex obtained with CollapseSurface(X, ∅).

The results obtained show that both algorithms presented above allows one to obtain a 2-dimensional skeleton from a 3-complex, containing important "shape information" from the initial object. However, even if Alg. 3 produces better results than Alg. 2, it fails in obtaining a completely satisfactory skeleton. In the next section, we introduce a methodology which includes a filtering step based on the Euclidean distance, providing a better robustness to contour irregularities and a lower sensitivity to the orientation of shape features in the space.

Skeleton filtering

A major difficulty when using the skeleton in applications (e.g., shape recognition, shape analysis), is its sensitivity to small contour perturbations, in other words, its lack of stability. A recent survey [START_REF] Attali | Stability and computation of the medial axis -a state-of-the-art report[END_REF] summarizes selected relevant studies dealing with this topic. This difficulty can be expressed mathematically: the transformation which associates a shape to its skeleton (or medial axis) is only semi-continuous. This fact, among others, explains why it is usually necessary to add a filtering step (or pruning step) to any method that aims at computing a skeleton. Hence, there is a rich literature devoted to skeleton pruning, in which different criteria were proposed in order to discard "spurious" skeleton points or branches.

In 2005, Chazal and Lieutier [START_REF] Chazal | The lambda medial axis[END_REF] introduced the λ-medial axis and studied its properties, in particular those related to stability. A major outcome of [START_REF] Chazal | The lambda medial axis[END_REF] is the following property: informally, for "regular" values of λ, the λ-medial axis of X remains stable under perturbations of X that are small with regard to the Hausdorff distance. Typical non-regular values are radii of locally largest maximal balls.

5.1 Intuitive presentation of the λ-medial axis Consider a bounded subset X of R n , as for example, for n = 2, the region enclosed by the solid curve depicted in Fig. 7 (left). The medial axis of X consists of the points of X that have several closest points on the boundary of X. For example in Fig. 7, the boundary points closest to x are a and b, the boundary points closest to x ′ are a ′ and b ′ , and the only boundary point closest to x ′′ is a ′′ .

Let λ be a non-negative real number, the λ-medial axis of X is the set of points x of X such that the smallest ball including all boundary points that are closest to x, has a radius greater than or equal to λ. Notice that the 0-medial axis of X is equal to X, and that any λ-medial axis with λ > 0 is included in the medial axis. We show in Fig. 7 Intuitively, in 2D, every "bump" on the boundary of the object generates a medial axis branch. The goal of the filtering is to eliminate those branches which are due to unsignificant contour irregularities (noise). Roughly speaking, in the λ-medial axis, the branches that remain correspond to contour features having a "width" greater than the parameter value λ.

Definition of the discrete λ-medial axis

The original definition of the λ-medial axis holds and make sense in the continuous Euclidean n-dimensional space. In [START_REF] Chaussard | A discrete lambda-medial axis[END_REF], we introduced the definition of a discrete λ-medial axis (DLMA) in Z n . We evaluated experimentally its stability and rotation invariance. Furthermore, we introduced a variant of the DLMA which may be computed in linear time, for which the results are very close to those of the DLMA.

Notice that the DLMA applies on a binary image (i.e., a set of voxels or a subset of Z 3 ), not on a complex. Let x, y ∈ R n , we denote by d(x, y) the Euclidean distance between x and y, in other words, d(x, y) = (

n k=1 (y k -x k ) 2 ) 1 2 . Let S ⊆ R n , we set d(y, S) = min x∈S {d(y, x)}.
Let x ∈ R n , r ∈ R + , we denote by B r (x) the ball of radius r centered on x, defined by B r (x) = {y ∈ R n | d(x, y) ≤ r}.

For each point x ∈ Z n , we define the direct neighborhood of x as N (x) = {y ∈ Z n | d(x, y) ≤ 1}. The direct neighborhood comprises 2n + 1 points.

Let S be a nonempty subset of Z n , and let x ∈ Z n . The projection of x on S, denoted by Π S (x), is the set of points y of S which are at minimal distance from x ; more precisely, Π S (x) = {y ∈ S | ∀z ∈ S, d(y, x) ≤ d(z, x)}. Let S be a finite subset of Z n , we denote by R(S) the radius of the smallest ball enclosing S, that is,

R(S) = min{r ∈ R | ∃y ∈ R n , B r (y) ⊇ S}.
Transposing directly the definition of the λ-medial axis to the discrete grid Z n would yield unsatisfactory results (see [START_REF] Chaussard | A discrete lambda-medial axis[END_REF]), this is why we need the following notion. Let S ⊆ Z n , and let x ∈ S. The extended projection of x on S, denoted by Π e S (x), is the union of the sets Π S (y), for all y in N (x) such that d(y, S) ≤ d(x, S).

Let S be a finite subset of Z n , and let λ ∈ R + . We define the function F S which associates, to each point x of S, the value F S (x) = R(Π e S (x)). The discrete λ-medial axis of S, denoted by DLM A(S, λ), is the set of points x in S such that F S (x) ≥ λ.

Fig. 6c shows the DLMA of the object presented in Fig. 5a, with λ = 18. For more details, illustrations and performance analysis, see [START_REF] Chaussard | A discrete lambda-medial axis[END_REF].

Surface skeleton based on the DLMA

Let us now describe our methodology to obtain a thin filtered surface skeleton from a binary image (a set of voxels). Let S denote our original set, a finite subset of Z 3 .

First, we compute the discrete λ-medial axis of S for a chosen value of λ. From here, we will consider the complex X = Φ(S) -and the set W = Φ(DLM A(S, λ)).

Notice that we have no guarantee that the complex W -be topologically equivalent to the complex X (see counter-examples in [START_REF] Chaussard | A discrete lambda-medial axis[END_REF]). Thus our next step consists of computing, using Alg. 1, the complex Y = ParDirCollapse(X, W, +∞) that contains W and that is topologically equivalent to X. Now, as W is made of 3-faces (see the close-up of Fig. 6c), Y is a 3-complex and cannot be seen as thin in a strict sense. Then to achieve our initial goal we compute Z = CollapseIsthmus(Y, ∅) which is both thin (following Prop. 1 there is no 3-face in this complex, see the close-up of Fig. 6d for an illustration) and topologically equivalent to X (by Cor. 4). The centering of Z in X is achieved thanks to the use of the DLMA, based on the Euclidean distance. The parameter λ can be tuned in order to adjust the filtering to the characteristics (size, smoothness of contours . . .) of the input shape, and to the requirements of the user.

Figures 6d and8 show various results using this methodology. It can be seen that the resulting 2-complexes indeed capture the main surfacic features of the original objects, without spurious branches or surface patches. A significant advantage of the 2D nature of the obtained skeletons, is to enable an easy analysis of important shape features such as intersections of surface parts.

Conclusion

We introduced in this work a methodology for obtaining filtered, thin surface skeletons of 3D objects. This methodology is based on a robust medial axis extraction method (DLMA), a parallel directional thinning algorithm based on collapse, and a surface detection strategy. The algorithms involved in this methodology can be implemented to run in linear time complexity.

Aknowledgements

This work has been partially supported by the "ANR BLAN07-2 184378 Mi-croFiss" project.

Fig. 1 .

 1 Fig. 1. Graphical representations of: (a) a 0-face, (b) a 1-face, (c) a 2-face, (d) a 3-face.

Fig. 2 .

 2 Fig. 2. (a): A binary image, represented by a set S ⊆ Z 2 (black pixels). (b): The complex X = Φ(S) -.

Fig. 3 .

 3 Fig. 3. (a-e): A series of complexes that illustrate a sequence of collapse operations.

Fig. 4 .

 4 Fig. 4. (a): A 2-complex X, each free pair (fi, gi) is symbolized by an arrow from gi to fi. (b):The complex obtained by removing simultaneously from X all the pairs that are free for X.

Fig. 5 .

 5 Fig. 5. (a): A 3-complex X having the shape of a thick tube with a skewed axis. (b): The result of ParDirCollapse(X, ∅, 10) is not thin, some 3-faces remain in the object. (c): The result of ParDirCollapse(X, ∅, 15) is a 2-complex. (d): The result of ParDirCollapse(X, ∅, ∞) does not contain any cube nor any square, it is composed of edges and points.

Fig. 6 .

 6 Fig. 6. Results of different algorithms performed on the complex X presented in Fig 5a. (a) : CollapseSurface(X, ∅). (b) : CollapseIsthmus(X, ∅). (c) : The discrete λ-medial axis (Sec. 5) of X, with λ = 18. (d) : Result of the method described in Sec. 5.3 performed on X, with λ = 18.

Fig. 7 .

 7 Fig. 7. Illustration of the λ-medial axis. Left: Points x, x ′ and x ′′ and their respective closest boundary points. Top right: λ-medial axis with λ = ǫ, a very small positive real number. Bottom right: λ-medial axis with λ = d(a ′ , b ′ ) + ǫ, where d(a ′ , b ′ ) is the distance between a ′ and b ′ .

Fig. 8 .

 8 Fig. 8. Results of the methodology described in Sec. 5.3. On the left, the original object; on the right, the result obtained by our method. First row: A bumped and skewed tube, λ = 20. Second row: A rocker arm, λ = 8. Third row: A pelvis bone, λ = 5.


Nicolas Ngan
email: ngann@esiee.fr

Geoffroy Marpeaux
email: marpeaug@esiee.fr

Eva Dokladalova
email: e.dokladalova@esiee.fr

Mohamed Akil
email: akilm@esiee.fr

Franc ¸ois Contou-Carrère
email: francois.contou-carrere@sagem.com

Memory system for a dynamically adaptable pixel stream architecture

Keywords: memory, FPGA, adaptable, embedded, image processing, ring

Nowadays, embedded vision systems have to face new hard requirements involved by modern applications: realtime processing of high resolution images issued by multiple image sensors. Recently, a new adaptable ring-based interconnection network on chip has been proposed. Based on adaptive datapath, it allows handling of multiple parallel pixel streams.

In this paper, we present a new hierarchical memory system proposed for this adaptable ring-based architecture. The design of its different levels is discussed and we show how the memory system adapts dynamically with respect to the datapath and data access management in the interconnection network. We also present the timing performance and area occupation measured on an FPGA prototype.

I. INTRODUCTION

Modern embedded vision systems have to face growing demands in image sensor multiplication [START_REF] Van Der Horst | A real-time stereo smartcam, using fpga, simd and vliw[END_REF]. Moreover, they have to implement a wide variety of applications (i.e. zoom, image filtering, fusion, face detection). Those applications introduce constraints in the amount of data, the processing time and latency (i.e. the delay allowed from input to output). Also, the system has to support variable real-time constraints from 25 to 50 frames per second with an image size up to the 1080p HD resolution (1920x1080). Such requirements represent the challenge of how to design the system to run efficiently multiple-sensor applications.

If it has been demonstrated in [START_REF] Chen | Coarse-grained reconfigurable image stream processor for digital still cameras and camcorders[END_REF] and [START_REF] Hauck | Reconfigurable Computing: The Theory and Practice of FPGA-Based Computation[END_REF] that the computing power does not represent a major problem, the bottleneck remains on the management of multiple streams and the efficient parallel data access. To illustrate, we can cite some critical applications such as image fusion [START_REF] Stathaki | Image Fusion: Algorithms and Applications[END_REF], [START_REF] Blum | Multi-Sensor Image Fusion and Its Applications[END_REF], multiscale approaches [START_REF] Cohen | Wavelets and Multiscale Signal Processing[END_REF] and especially temporal processing with multiple delayed frames.

Recently, we have proposed a new adaptable ring-based interconnection network on chip [START_REF] Ngan | Dynamically adaptable architecture for real-time video processing[END_REF]. The main principle at the system level is shown in Figure 1. An application loader precharges the application context into the processing system from the application memory. The system then dynamically adapt its datapath according to the running imaging application. This datapath adaptation allows handling parallel computation for multiple streams. It involves the need for a performant adaptable memory system. In particular, that memory system has to be adapted dynamically in datapath and data access. In this paper, we propose a new memory system applied on the published adaptable ring-based architecture. That memory system is dynamically adaptable and it is organized in two levels: a frame level and a line level. The frame memory level allows data accesses by blocks of data and the line level allows random data accesses.

The paper is organized as follows. The Section II is a presentation of the Multi Data Flow Ring (MDFR) architecture. MDFR is an original adaptable ring-based interconnection network on chip. The dynamically adaptable memory system organized in two levels is described in Section III and IV. Results in space and timing performance are then presented for its first hardware prototyping in Section V.

II. ADAPTABLE RING-BASED PIXEL STREAM ARCHITECTURE

An imaging application is defined as one or a set of associated algorithms. We can cite for example a 2D convolution. An algorithm describes a sequence of operations to carry out by a specified task T i . Thus, the imaging application is described as a predefined sequence of a set of tasks {T i }. From data flow point of view, that sequence defines an association of different execution modes: pipeline, parallel and pipelineparallel as shown in Figure 2. We call a pipeline execution mode a highly optimized datapath used to connect in a serie computing ressources in order to process a given pixel stream. A parallel execution mode is considered when there are datapath divisions to have computing ressources processed tasks in parallel. Finally, a pipeline-parallel execution mode is a combination of pipeline and parallel execution mode.

Let consider in the following of this paper that those tasks T i are implemented in a processing element P E i and that the execution modes are realized by a defined datapath. Depending on the application, those execution modes might also have an impact on the memory requirements. For instance in parallel execution mode, processed parallel streams might need to be synchronized and a memory is then required. The size of the memory can vary from a few lines to a complete frame. Also, each T i can require different manners to access data.

We distinguish two different data access modes in memory: data access by lines and data access by pattern. We define a data access by lines as a linear access and data access by pattern as a random access for accessing a specific block of data.

A. Architecture description

The proposed architecture called Multi Data Flow Ring (MDFR) is described in Figure 3. It is a dynamically adaptable architecture for real-time imaging applications in embedded vision systems. It can be adapted in both datapath and memory accesses. The architecture contains two data router types: Data Flow Routers (DFR) and I/O Stream Gate Manager (SGM) nodes linked together forming a ring topology. We define a ring channel as a unidirectionnal communication link between routers. DFRs and SGMs are connected through four data ring channels : two rings oriented in clockwise direction and the two others in the inverse direction as shown in Figure 3. Those two directions allow the pixel stream to reach a PE through DFR from both directions in the shortest possible datapath.

Main datapath adaptations are made with the DFRs to realize application execution modes. A DFR is a router in charge of redirecting the pixel streams to its connected processing element P E i executing a task T i . Its main actions are switching input pixel streams for its connected PE or forwarding it to the next DF R in the ring. It contains buffers to synchronize input and output streams, mutiplexers, configuration registers and pixel counters to track the status of PE data output.

SGMs are the master units of the architecture. They control the pixel streams in the ring channels. A SGM is a particular data router that contains dedicated ports for external input and Fig. 3. Multi DataFlow Ring Architecture output streams. Each SGM has an access to the Frame Buffer System (FBS) which is in charge to store and manage the frames from pixel streams.

A SGM has two levels of memory (Figure 4). The first level is the Line Memory System (LMS) inside the SGM and the second level are the Frame Memories inside the Frame Buffer System (FBS). The main controller inside the SGM can make a read or write request to the controller called Frame Server inside the FBS. Pixels stored in LMS are then transferred from the FBS or to the FBS. Each memory level has specific controllers ensuring memory accesses and datapath adaptations. Table I summarizes the memories characteristics (memory type, master, data access mode). For the SGMs, the LMS is considered as a local memory and the FBS as a shared global memory. The LMS can be read in a predefined way to extract different patterns from input frames to be sent in ring channels (Figure 3). In the following of this paper, we focus on the memory system adaptations (datapath and access mode) with the architecture description of the SGM and the FBS.

III. STREAM GATE MANAGER ARCHITECTURE

A SGM interfaces three types of pixel streams sources from ring channels ports, the external port or the FBS port. Let K c be the number of ring channels connected to the SGM, K ext the number of external input ports and K fbs the number of input ports from FBS. Consequently, the maximum number of input pixel streams K to manage by the SGM would be :

K = K c + K ext + K fbs (1)
An overview of SGM architecture is given in Figure 5. We can see that Channel 1 and 2 are the two channel rings oriented in clockwise direction Channel 3 and 4 those oriented in the inverse direction. In this architecture, each SGM has one input port from FBS (K fbs = 1) and one external input port(K ext = 1). Two ring channels are connected for each direction (K c = 4). Consequently, a SGM has to manage at maximum K = 6 input pixel streams. Those values are chosen in order to minimize the SGM size and to have the possibility to realize any execution mode (pipeline, parallel, pipeline-parallel).

The SGM contains a Line Memory System (LMS) which is a set of local memories for buffering incoming pixel streams from ring channels, external and FBS. The SGM main controller called the Pixel Stream Controller interfaces each ring channel inputs (Figure 5). Its role is to analyze incoming pixel streams and to check the next datapath and memory access adaptation instruction in the Instruction Memory. From those instructions, it then sends requests to the Datapath and Memory Access Configurator. This unit contains registers that inform the datapath configuration in the LMS. It decides dynamically datapath and memory access modifications in LMS according to occupied channels. The Pixel Stream Controller can also request a reading of a new pixel stream or a writing of an incoming pixel stream to FBS. Reading and writing on FBS are controlled by specific controller units called Frame Reader and Frame Writer. Architecture descriptions of the Pixel Stream Controller and the Datapath and Memory Access Configurator units are out of the scope of this article focusing only on the memory system.

A. Line Memory System Architecture

The LMS architecture is presented in Figure 6. It allows both datapath and memory adaptations for six input pixel streams. Datapath adaptations are realized with four 2-input multiplexers, one 4×4 crossbar, and 4-input multiplexers. The LMS contains four memory units called Line Memory Channel (LMC). As shown in Figure 6, each LMC is connected to one ring channel at its output. Consequently, the Datapath and Memory Access Configurator (Figure 5) has to select the ring channel output with the crossbar (Figure 6). Note that pixel streams from ring channels can bypass direcly the LMS without any stream modification from input registers to output registers. That possibility is important when a pixel stream has to cross several times SGMs in the ring for reusing or finding a particular PE. Each LMC can accept up to one source type pixel stream from ring channel, external and FBS. A LMC can output pixel stream to a ring channel, to external output or to the FBS for writing. A LMC contains two main memories: Memory Up and Memory Down. Memory Up is used to store streams to be sent to a ring channel and Memory Down is used either to ouput the stream to an external port or to write the buffered stream to the FBS. Thus, a pixel stream from ring channels can be oriented to any LMC to access an external port or to be written in FBS. The memory access adaptation is decoded by the Pixel Stream Extractor. As shown in Figure 7, it contains a Command Decoder unit which receives datapath and memory access requests from the Datapath and Memory Access Configurator in the SGM (Figure 5). Datapath is decoded to configure the input mux a and b to select the right pixel streams to process. The adress generator units for Memory Up and Memory Down are then configured to write input pixel streams linearly and to realize an access mode. An access mode by lines configures the adress generators unit to read the pixels linearly. Patterns are predefined for each address generators and memory addresses are computed according to the line width value of the frame. For instance, it can be configured to jump pixels for a fast frame downscaling or access particular block of pixels. Obviously, the pixel block size access is limited by the size of Memory Up or Down for a giving line width value and pixel granularity. A defined access mode configures also the Pixel Format units according to the pixel granularity in the stream such as 8, 16 or 32 bits only. The Pixel Format unit adapts the output pixels data of the Memory Up and Down to the ring channel data size by using byte masks and succesive required concatenations. For instance, four 8-bits pixels to get the size of a 32-bits data in the ring channels.

B. Datapath and Memory Access adaptation control

LMS has two levels of datapath control. The first level is selecting the right LMC for pixel streams coming from ring channels and redirecting the LMC streams output either to a ring channel, to external or to FBS (Figure 6). The second level interests the selection of Memory Up and Memory Down inside a LMC (Figure 7). For an input pixel stream, there can be five datapath instruction requests from the PixelStream controller:

• FW(D) : Forward the input pixel stream to the direct associated ring channel • FW(C) : Change the direction of the input pixel streams to another ring channel • READ : Read another pixel streams in FBS to be processed • STORE : Store the input pixel stream to FBS • OUT : Exit the input pixel stream to the external port Table II summarizes the cost in memory (Memory Up and Down) and datapath utilization (Ring channels) for an input pixel stream from a ring channel. The Datapath and Memory Access Configurator then takes a decision on the datapaths according to the maximum input streams capacity K and the availability of Memory Up and Down. After selecting the datapath, if a pixel stream has to be stored in Memory Up or Down, the Datapath and Memory Access Configurator has the possibility to define a specific access mode by lines or patterns by configuring the address generators in each LMC.

(D) 0 0 1 0 0 0 FW(C) 1 0 1 0 0 0 READ 1 0 0 0 1 0 STORE 0 1 0 0 0 1 OUT 0 1 0 1 0 0

IV. FRAME BUFFER SYSTEM ARCHITECTURE

A. Frame Slot Memory definition

We define a frame slot as a frame at a specific time index ts. The time index ts represents the number of frames delayed from the current processed frame. Frames are stored a frame memory at a physical memory address. A frame memory can then contains several frame slots.

Let M and N be the width and height of a frame slot. Let r be the pixel data size in bits. The memory block size S F would be :

S F = M * N * r (2)
We call frame slot memory the memory block of size S F storing a frame slot. For example, the frame slot memory size for a complete HD greyscale image of 1920x1080 in 8 bits would be about 2 Mbits. Thus, let S LM the size of the frame memory. The number P of available frame slot memory in the frame memory would be :

P = S LM S F ( 3
)
The Frame Buffer System (FBS) is based on the proposition to divide dynamically a frame memory in several frame slot memories. Each frame slot memory is identified by attributes such as the time index ts, the pixel granularity id and the processing element n which has processed the frame slot.

B. Frame Buffer System architecture

The FBS architecture is presented in Figure 8. It contains a Frame Server which receives commands (read or write) from SGMs and controls the adaptation units. The adaptation units are the The Memory Access Controller and the Data Ouput Controller. The Memory Access Controller generates the adresses to access the frame memories and the Data Ouput Controller adapt dynamically the frame memories ouput datapaths with a crossbar to reach any SGMs.

The FBS contains as many frame memory as the number of SGMs in the ring (See Figure 8 with four SGMs). Each SGM is connected to its private frame memory in writing. The set of frame memories are designed to get a private access in writing but a shared access in reading for all SGMs.

Each frame memory is divided in frame slot memories and those slots are managed by the Frame Slot Manager. All frame slots memory attributes (time index ts, the pixel granularity id and the processing element n) are registered The FBS architecture has several advantages. First, the physical memory adressing is completely transparent for the SGMs. Frame read or write requests to the Frame Server are done indirectly with the frame slot attributes and not with the physical memory address. Thus, it allows reusability and evolution to any frame memory size in the FBS. The frame management stays completely independent. The Frame Register Tables size only depends on the defined number of frame slot memory to manage and not the memory size. Secondly, it allows processing a full frame by blocks sequentially or in parallel. If the incoming frame size is higher than the defined frame slot memory size, then it is automatically divided and considered as several lower frame slots to be processed. It is particularly interesting for processing elements limited in frame size. Thus, imaging applications can also be accelerated by propagating the pixel stream to each SGM and duplicating the frame to each frame memory. Each SGMs can then access in parallel to separate parts of the frame considered as lower frame slots. For instance, an HD 1920 × 1080 image can be considered as sixteen 480 × 270 frame slots to be processed.

C. FBS Datapath and Memory Access control

The FBS uses adaptation algorithms for reading and writing on the different memories. It contains a reading queue in case of simultaneous reading requests on the same frame memory.

For writing a frame, the Frame Server has to decode the attributes of the requested frame slot (the time index ts, the pixel granularity id and the processing element n) and the pixel stream size q (number of pixels). It asks the Frame Slot Manager for a free frame slot memory. From the frame slot attributes, the Frame Slot Manager can address the Frame Register Tables and returns to the Frame Server the physical memory address of the last free slot memory (Figure 9). The physical memory address is then transmitted to the Memory Access Controller with q to initialize the writing adress generator unit of the frame memory. During the writing, the Frame Slot Manager updates the table by linking the frame slot memory adress to the attributes.

For reading a frame slot, the Frame Server has to decode also first the attributes to get the physical memory address. If the frame memory is still busy, it puts the read command in the reading queue. If the frame memory is available, then the Frame Server configures the Data Output Controller to adapt the frame memory output datapath. The reading adress generator unit of the frame memory is configured with the physical memory address and the parameter q.

V. HARDWARE PROTOTYPING

The presented adaptable memory system has been prototyped on an Altera FPGA Stratix III EP3SL150. For evaluation of this implementation, we focus on the memory system adaptation cost in area (logic elements) and time (latency). We consider working with frame slots of size 10 × 100 Bytes (B). The whole system is working at the same frequency (maximum frequency of 210 MHz) and we use the on-chip memory to implement the memory system. We define the ring channels data size to 32 bits and the instruction memory to contain 8 × 1 B instructions.

A. SGM implementation

All the SGM control units (Pixel Stream controller, Frame Reader, Frame Writer and Datapath and Memory access configurator) are implemented as Finite-State Machines (FSMs). The Pixel Stream controller contains 16 B FIFO memories to store pixels during adaptations. The Memory Up and Down in LMCs are implemented as 32-bits dual-port memories to fit 2 × 100 B (two lines). Table III shows the implementation results for one SGM with the parameters K c = 4, K ext = 1 and K fbs=1 . Note that the results for the control units depend only on the parameters K c , K ext and K fbs . Table IV shows the SGM latency in clock cycles for one pixel stream in a pipeline execution mode and for each adaptation requests described in Table II. To cross a SGM, each stream has to wait the datapath and memory adaptations inside the SGM. The fastest adaptation is forwarding directly (FW(D)) the stream because the datapath in LMS is just modified at the first level in 2 clock cycles. If the stream has to be stored an LMC (i.e. with instructions FW(C), READ, STORE, OUT) then it adds 4 clocks cycles to communicate with the Pixel Extractor and wait for its handshake (signal ack in Figure 6 and7). Reading or writing to the FBS adds 2 clock cycles to initialize the Frame Reader and Writer (Figure 5). β is the time for the address generators to be initialized for the access pattern. In our case, β = 1 for a linear access and β = 3 for accessing a block. The time α is the latency for buffering the pixels in LMC according to the access mode. For instance, for accessing the frame each two lines of 100 pixels, α = 50 clock cycles. Note that using parallel or pipeline-parallel execution mode only impact the datapath latency for FW(D) to 6 because the stream has to be stored in LMC in order to use a second ring channel.

B. FBS implementation

For the FBS implementation, we use four frame memory implemented as 32-bits dual-port memories of size 16384 B to fit 10 frame slots of 10 × 100 B. Frame Server and Frame Slot Manager are implemented as FSMs. The Memory Access controller is a FSM with several counters for each address generators units. Table V presents the FBS implementation results for managing 10 Frame slots. Note that the size of frame slot does not have any impact on the size of control units, datapath adaptation unit and memory access controller.

Table VI shows the FBS latency for read and write requests. Checking the Frame Register table requires 3 clock cycles to get the translation between the attributes to the physical memory address. A register store the last free frame slot to From those results, we can conclude that the memory system latency cost for an input stream to cross a SGM is from 4 to 14 clock cycles. The latency cost for the FBS is 8 clock cycles and it can be overlapped by the SGM adaptation latency. Obviously, the more combined instructions or number of input streams are used and the more efficient the adaptation is. For long pixel stream pipeline, thoses latencies become negligible and the ring structure remains performant. Besides, those small latencies allows to have reduced sized FIFO buffers for each SGM in the ring.

VI. CONCLUSION

In this paper, we have presented a new adaptable hierarchical memory system to address streaming imaging applications with several pixel streams. It can then be adapted dynamically to manage different frames from image sensors or delayed frames stored in memory. The memory adaptation is both based on datapath and memory accesses to satisfy each execution mode. A first FPGA prototype has been presented and evaluated. Future works will focus on power consumption.

Fig. 1 .

 1 Fig. 1. Adaptable Ring-Based Architecture Context

Fig. 2 .

 2 Fig. 2. Illustration of different execution modes

Fig. 4 .

 4 Fig. 4. Memory Hierarchy Levels

Fig. 5 .

 5 Fig. 5. I/O Stream Gate Manager (SGM) architecture

Fig. 6 .

 6 Fig. 6. LMS architecture : Kc = 4, Kext = 1 and K fbs = 1 1) Line Memory Channel architecture: The LMC architecture is presented in Figure 7.

Fig. 7 .

 7 Fig. 7. Line Memory Channel architecture

Fig. 8 .

 8 Fig. 8. Frame Buffer architecture overview

Fig. 9 .

 9 Fig. 9. Frame Register Tables use for frame slot 1

  . δ 1 is the waiting time for reading if the frame memory is busy (from 0 to 250 clock cycles equivalent of reading one 10 × 100B frame). δ 2 is the waiting time for writing if the Frame Buffer is still updating the Frame Register Tables (from 0 to 30 clock cycles for 10 slots).

		TABLE V
	FBS: AREA ESTIMATION (ALTERA EP3SL150)
		ALUTs	Registers Memory (bits)
	Control units	157	308	197632
	Datapath Adaptation unit	264	0	0
	Memory access controller	160	256	0
	Frame Memories (4)	0	0	524288
	FBS (Total)	581	564	721920
	FBS Occupation (%)	0.5	0.4	12.8
	accelerate the writing