Yukiko Kenmochi
email: y.kenmochi@esiee.fr

Atsushi Imiya
email: imiya@faculty.chiba-u.ac.jp

Combinatorial Boundary of a 3D Lattice Point Set

Keywords: boundary extraction, combinatorial surface, polyhedral complex

come

Introduction

Several algorithms have been presented for boundary extraction [START_REF] Artzy | The theory, design, implementation, and evaluation of a three-dimensional surface detection algorithm[END_REF][START_REF] Kong | Digital topology: introduction and survey[END_REF][START_REF] Kovalevsky | Finite topology as applied to image analyses[END_REF][START_REF] Voss | Discrete Images, Objects, and Functions in Z ![END_REF] and surface generation [START_REF] Artzy | The theory, design, implementation, and evaluation of a three-dimensional surface detection algorithm[END_REF][START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF][START_REF] Franc ¸on | Discrete combinatorial surfaces[END_REF][START_REF] Kong | Continuous Analogs of Axiomatized Digital Surfaces[END_REF][START_REF] Lorensen | Marching cubes: a high-resolution 3d surface construction algorithm[END_REF][START_REF] Morgenthaler | Surfaces in three-dimensional images[END_REF][START_REF] Wyvill | Data structure for soft objects[END_REF] of 3-dimensional digital images for purposes of visualization, calculation of geometric features such as surface areas, calculation of topological features such as Euler characteristics, numerical analysis for deformable objects, etc. Even if both topological objects, borders and surfaces, are required simultaneously (sometimes implicitly) for many applications as listed above, it is not easy to nd a useful theory allowing to discuss both topological concepts for each point in a 3-dimensional lattice space, i.e. for each voxel in a 3-dimensional digital image. Note that there are many approximation techniques, but we are interested in completely discrete techniques because our input are digital images and our computations for image analysis are also digital. Furthermore, such discrete techniques may bring us special geometric and topological properties which will be seen only in discrete spaces. The eld of those studies are called discrete/digital geometry and topology [START_REF] Kong | Digital topology: introduction and survey[END_REF] and useful properties may provide us new ef cient algorithms for the applications listed above.

Even in the Euclidean space, it is not easy to draw relations between borders in the sense of general topology and surfaces in the sense of combinatorial topology [START_REF] Moise | Geometric Topology in Dimensions 2 and 3[END_REF]. More discussions on the historical backgrounds may be found in Section 2.

In this paper, we tackle a problem for clarifying relations between border points and surface points, i.e. points which are vertices of polyhedral surfaces, in a 3dimensional lattice space. To solve the problem, we use polyhedral complexes such that all vertices are lattice points and the adjacent vertices are neighboring each other in the sense of 3-dimensional digital topology [START_REF] Kong | Digital topology: introduction and survey[END_REF]. Such polyhedral complexes are called discrete polyhedral complexes and they enable us to give a topology or a polyhedral surface to a set of border points. By using discrete polyhedral complexes, we also present a new algorithm for extracting border points which constitute a polyhedral surface. We, therefore, succeed to extract border points and their surface structures simultaneously.

The de nition of border points is based on general topology [START_REF] Hausdorff | Set Theory[END_REF][START_REF] Moise | Geometric Topology in Dimensions 2 and 3[END_REF] and it has been shown that we can obtain border points by a set operation using neighborhoods [START_REF] Kong | Digital topology: introduction and survey[END_REF][START_REF] Soille | Morphological Image Analysis: Principles and Applications[END_REF][START_REF] Voss | Discrete Images, Objects, and Functions in Z ![END_REF] (see [START_REF] Franc ¸on | Discrete combinatorial surfaces[END_REF] in Section 2). Because we need to carry out the set operation for each point in a 3-dimensional lattice space, i.e. each voxel in a 3-dimensional digital image, the computational time is linear in the size of a digital image.

In two dimensions, some ef cient border tracking algorithms have been proposed by using curve structures of border points [START_REF] Kong | Digital topology: introduction and survey[END_REF][START_REF] Rosenfeld | Connectivity in digital pictures[END_REF] such that a set of border points is given as a sequence of points (or pixels) and each point (or pixel) has exactly two neighboring points (or pixels). Each border point is tracked by a left-hand-onwall border following algorithm from the previous point in a sequence; therefore, we do not have to scan all points in a whole digital image; the computational time becomes linear in the number of border points.

In three dimensions, a completely different approach from that in two dimensions is commonly used because of the dif culty of nding surface structures of border points. An algebraic-topology based approach is taken so that unit cubes (or voxels) whose centroid are lattice points are rst considered and then for border tracking the common faces between two voxels centered at points F in an object region and G in a complement of the region are considered [START_REF] Artzy | The theory, design, implementation, and evaluation of a three-dimensional surface detection algorithm[END_REF]. Such faces are represented by an ordered pair (F, G). Therefore, boundaries are represented by surfaces which are sets of square faces and whose topological structures are given as cellular complexes as shown in [START_REF] Kovalevsky | Finite topology as applied to image analyses[END_REF]. We can also consider the set of all points F of such pairs (p, q) as the border [START_REF] Voss | Discrete Images, Objects, and Functions in Z ![END_REF]. However, for such a set of points, we can not obtain any topological surface structures.

There are some axiomatic de nitions of discrete surfaces such that all points of discrete surfaces are lattice points and not voxel faces [START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF][START_REF] Franc ¸on | Discrete combinatorial surfaces[END_REF][START_REF] Kong | Continuous Analogs of Axiomatized Digital Surfaces[END_REF][START_REF] Morgenthaler | Surfaces in three-dimensional images[END_REF]. However, the relations between border points and those surface points are not yet clari ed; for example, we can nd easily some border points which cannot be points of discrete surfaces de ned in [START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF], called simplicity surfaces, as shown in Figure 1. Conversely, the connectedness of border points are shown in [START_REF] Kong | A Theory of Binary Digital Pictures[END_REF][START_REF] Kong | Concepts of digital topology[END_REF], but the concept of connectedness is clearly not suf cient for providing surface structures.

In order to give the relations between border points and surface points, we need a different approach. In the sense of combinatorial topology [START_REF] Alexandrov | Combinatorial topology I[END_REF], this is a special formulation of a triangulation problem for border points. After presenting historical backgrounds in the next section, we de ne 3-dimensional discrete polyhedral complexes and give the combinatorial boundary which contains 2-dimensional surface structures in Section 3. In Section 4, we present an algorithm to provide a combinatorial boundary from any given 3-dimensional lattice point set. Because our algorithm is similar to a marching cubes algorithm [START_REF] Lorensen | Marching cubes: a high-resolution 3d surface construction algorithm[END_REF][START_REF] Wyvill | Data structure for soft objects[END_REF] using a look-up table, our computational time is linear in the size of a 3-dimensional digital image. We then derive the relations between borders in the sense of general topology and our combinatorial boundaries. From these relations, we nally conclude that our combinatorial boundary extraction algorithm gives a triangulation for border points, simultaneously with border points, with respect to a given 3-dimensional lattice point set.

U A (x) = {y ∈ 4 n : x -y < A} (1)
In [START_REF] Hausdorff | Set Theory[END_REF], the term boundary is used instead of frontier. In this paper, we keep the term boundary for combinatorial boundary in the sense of combinatorial topology and follow the terminology in [START_REF] Moise | Geometric Topology in Dimensions 2 and 3[END_REF] to distinguish between boundaries in general topology and combinatorial topology.

of radius e > 0 de ne a basis of open sets for this Euclidean space.

If a point x in A ⊂ R n is such that there exists a neighborhood U A (x) ⊆ A, then it is called an interior point of A. Otherwise, a point x ∈ A is called a border point of A. Let I nt (A) and Br(A) be the sets of all interior and border points such that

I nt (A) = {x ∈ A : U A (x) ⊆ A}, (2)
Br(A) = A\Int(A), (3)
called the interior and border of A, respectively. Then we have A = I nt (A)∪Br(A).

Let Ā be the complement of A such that R n = A ∪ Ā. Then, the interior points of Ā are also the exterior points of A. The union of the borders of A and Ā yields the frontiers Fr(A) and Fr( Ā) such that

Fr(A) = Br(A) ∪ Br( Ā) = Fr( Ā). ( 4
)
Figure 2 shows examples of the border and frontier of a point set A in R .

In this paper we also consider the combinatorial boundary of an n-dimensional polyhedral complex K [START_REF] Alexandrov | Combinatorial topology I[END_REF][START_REF] Moise | Geometric Topology in Dimensions 2 and 3[END_REF] so that we treat boundaries as (n -1)-dimensional quasi-manifolds [START_REF] Kovalevsky | A topological method of surface representation[END_REF] as shown in Figure 1. If K is a 3-dimensional polyhedral complex, then the combinatorial boundary ¶K is the set of all 2-polyhedra s of K such that s lies only in one 3-polyhedron of K, together with all of its faces.

The precise de nitions of polyhedral complexes and combinatorial boundaries will be given in Section 3. Let K be a 3-dimensional polyhedral complex which is a triangulated 3-manifold with boundary and |K| be the union of the elements of K, with the subspace topology induced by the topology of R n . Then, the relation between the frontier and the combinatorial boundary is derived such that

| ¶K| = Fr(|K|) (5)
if |K| is closed; see [START_REF] Moise | Geometric Topology in Dimensions 2 and 3[END_REF] for the proof.

2.1.2

Borders and Boundaries in Z and Z ! Let us consider the set Z n of all lattice points in R n such that their coordinates are all integers. For any point set V ∈ Z n , which is given as an object component in an n-dimensional binary image, borders are also de ned similarly to borders in R n . In this paper, we consider the cases n = 2,3.

Traditionally, the following m-neighborhoods

N m (x) = {y ∈ Z n : x -y ≤ t } with t = 1, √ 2 (resp. t = 1, √ 2,
√ 3) are in common used for x ∈ Z (resp. x ∈ Z ! ), and m = 4,8 (resp. m = 6,18,26) stands for the cardinality of these neighborhood systems [START_REF] Kong | Digital topology: introduction and survey[END_REF]. In distinction to e-neighborhoods of (1), the radius t is only one of the three numbers 1, √ 2 or √ 3. It follows that these m-neighborhoods do not establish a basis of open sets of a topology on Z n , and that image analysis normally only assumes adjacency graphs in Z n for de ning concepts of connectedness [START_REF] Kong | Digital topology: introduction and survey[END_REF].

Let m ∈ {4, 8} for n = 2 and m ∈ {6, 18,26} for n = 3. If a point x in V ⊂ Z n is such that N m (x) ⊆ V, then x is called an interior point (with respect to mneighborhoods) [START_REF] Kong | Digital topology: introduction and survey[END_REF][START_REF] Voss | Discrete Images, Objects, and Functions in Z ![END_REF]. The set of interior points of V is called the interior of V and denoted by

Int m (V) = {x ∈ V : N m (x) ⊆ V} (6)
similarly to (2) for A ⊂ R n . If a point x ∈ V is not an interior point of V, then x is called a border point of V with respect to m-neighborhoods [START_REF] Kong | Digital topology: introduction and survey[END_REF][START_REF] Voss | Discrete Images, Objects, and Functions in Z ![END_REF]. The set of all border points of V is called the m-border of V, denoted by

Br m (V) = V\Int m (V). (7)
Equation ( 7) corresponds to [START_REF] Alexandrov | Combinatorial topology I[END_REF].

In terms of mathematical morphology [START_REF] Soille | Morphological Image Analysis: Principles and Applications[END_REF] it follows that an interior set Int m (V) of ( 6) coincides with the erosion of V with the structure element N m (o) where o is the origin of Z n [START_REF] Voss | Discrete Images, Objects, and Functions in Z ![END_REF]. We see that (7) also de nes Br m (V) via a set operation such as

Br m (V) = {x ∈ V : N m (x) ∩ V = / 0}. ( 8
)
This formulation is possible because the radii t of N m (x) is constant in Z n . In consequence, no set operation corresponding to (8) exists for Br(A) in R n of (3). Figure 3 shows examples of the 4-borders of V ∈ Z and of the complement V.

Let us consider the boundary points of V ⊂ Z n , corresponding to the frontier points of A ⊂ R n , in the sense of general topology. From (4), a point set A ⊂ R n and the complement Ā has the frontier which is the common boundary as shown in Figure 2 (d). Similarly, we can de ne the m-boundary of V as the union of the m-borders of V and of V. Such m-boundaries are used for the composition of boundaries by contributions from both participating sets [START_REF] Lorensen | Marching cubes: a high-resolution 3d surface construction algorithm[END_REF][START_REF] Sloboda | On the topology of grid continua[END_REF]. In digital image analysis, however, Br m (V) and Br m ( V) are considered separately [START_REF] Kong | Digital topology: introduction and survey[END_REF][START_REF] Rosenfeld | Connectivity in digital pictures[END_REF][START_REF] Voss | Discrete Images, Objects, and Functions in Z ![END_REF] not only for

We follow the terminology in [START_REF] Kong | Digital topology: introduction and survey[END_REF], even if the term inner point is used instead of interior point in [START_REF] Voss | Discrete Images, Objects, and Functions in Z ![END_REF], to make a correspondence between interior points in R n and Z n . boundary tracking but also for thinning ! . They are called internal and external m-boundaries, respectively [START_REF] Soille | Morphological Image Analysis: Principles and Applications[END_REF].

Such existence of different boundaries in a discrete space has been already pointed out by W. K. Clifford. In [START_REF] Clifford | The postulates of the science of space[END_REF], he explained it using an example of a heap of white marbles on the top of which black marbles are put. The boundary of the white part would be a layer of white marbles and the boundary of the black part would be a layer of the black marbles, that is, the two adjacent parts have different boundaries when they are divided into two parts. He also referred to the Aristotelian de nitions of continuity and discontinuity: the continuity as that of which two adjacent parts have the same boundary; and the discontinuity or discreteness as that of which two adjacent parts have different boundaries. This paper adopts the boundary approach as follows: the internal boundary, simply called m-border, Br m (V) de nes the initial point set, we consider the combinatorial boundary of an n-dimensional polyhedral complex K such that all vertices are in Br m (V). We call such n-dimensional polyhedral complexes discrete polyhedral complexes to distinguish them with other (general) polyhedral complexes. This is a special formulation of a triangulation problem for Br m (V) in the sense of combinatorial topology [START_REF] Alexandrov | Combinatorial topology I[END_REF]. In the following sections, we give de nitions which are necessary for construction of discrete polyhedral complexes and present a solution for the triangulation problem. We afterwards derive the relations between Br m (V) and the combinatorial boundary of a discrete polyhedral complex ¶K with respect to the relation of ( 5) for R n (see Theorem 6 and Corollary 7 in Section 5). A set B ⊂ Z n is said to be connected or m-connected if any pair of x, y ∈ B has a point sequence x = x, x , . . . , x k = y such that all x i ∈ B and x i+ ∈ N m (x i ) [START_REF] Kong | Digital topology: introduction and survey[END_REF].

In Z , it is known that the m-border Br [START_REF] Franc ¸on | Discrete combinatorial surfaces[END_REF][START_REF] Artzy | The theory, design, implementation, and evaluation of a three-dimensional surface detection algorithm[END_REF] [START_REF] Kong | Digital topology: introduction and survey[END_REF][START_REF] Rosenfeld | Connectivity in digital pictures[END_REF]. All m-border points are then tracked as a sequence of points such as x , x , x ! , . . . and every point x i in the sequence is found as an element of m ′ -neighborhood of the previous point x i- [START_REF] Kong | Digital topology: introduction and survey[END_REF]. We see in this approach that the de nition of a curve is implicitly given as a sequence of points. In other words, border points are tracked by using the curve structures in Z .

m (V) of V ⊂ Z is m ′ -connected if V is m ′ - connected without hole where (m, m ′ ) = (4, 8),
In Z ! , it has been shown in [START_REF] Kong | A Theory of Binary Digital Pictures[END_REF][START_REF] Kong | Concepts of digital topology[END_REF] that the m-border Br m (8) of 8 ⊂ Z ! is m ′connected if 8 and 8 are m ′ -and m-connected respectively for any pair (m, m ′ ) ∈ {6, 18,26} × {6, 18,26} \ (6, 6). " Similarly to the case of two dimensions, for border tracking of 8 in Z ! , we need a de nition of a surface instead of that of a curve in Z . Clearly, the connectivity is not enough for representing the structures of surfaces such as triangulated surfaces.

2.2.2

Surface Representation in Z ! The de nition of surfaces in Z ! is more complicated than that of curves in Z . There exist various de nitions of two-dimensional surfaces in Z ! . The approaches are mainly classi ed into the following four types:

(1) the graph-theory based approach: a surface is de ned as a set of lattice points which satis es some conditions based on the neighborhood relations or the connectedness [START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF][START_REF] Morgenthaler | Surfaces in three-dimensional images[END_REF][START_REF] Voss | Discrete Images, Objects, and Functions in Z ![END_REF]. Every point on surfaces is considered to have a characteristic of spatial separation according to the Jordan surface theorem in a local sense. (2) the algebraic-topology based approach: surfaces are de ned as the combinatorial boundaries of 3-dimensional cellular complexes. In [START_REF] Artzy | The theory, design, implementation, and evaluation of a three-dimensional surface detection algorithm[END_REF][START_REF] Herman | Discrete multidimensional Jordan surfaces[END_REF][START_REF] Kovalevsky | Finite topology as applied to image analyses[END_REF], cells are considered to be unit cubes (or voxels) whose centroids correspond to lattice points and surfaces are represented by sets of faces of unit cubes. In [START_REF] Kenmochi | Boundary extraction of discrete objects[END_REF], simplicial complexes are used instead of cellular complexes so that the vertices of simplexes are all lattice points. (3) the combinatorial-manifold based approach: surfaces are triangulated and any point on a surface is topologically equivalent to the central point of an open disc [START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF][START_REF] Franc ¸on | Discrete combinatorial surfaces[END_REF][START_REF] Kong | Continuous Analogs of Axiomatized Digital Surfaces[END_REF]. (4) the analytical approach: geometric surfaces such as planes and spheres are de ned by using inequalities in Z ! instead of using equations in R ! [START_REF] Andres | Discrete Analytical Hyperplanes[END_REF][START_REF] Andres | The Discrete Analytical Hyperspheres[END_REF][START_REF] Debled-Rennesson | A new approach to digital planes[END_REF].

The analytical approach can be applied if only geometric objects such as planes and spheres are considered. In this paper, we would like to treat any free-form objects. Thus, we cannot take the analytical approach.

The graph-theory based approach is the most classic, but is also axiomatic. Since it contains only neighborhood relations and not topological structures, the combinatorialmanifold based approach has been taken in [START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF][START_REF] Kong | Continuous Analogs of Axiomatized Digital Surfaces[END_REF] for making comparison between the graph-theory based approach and the combinatorial-manifold based approach. Clearly, the combinatorial-manifold based approach has the strong power for inves-tigating topological structures, but it is not evident that a set of border points can become a combinatorial manifold. For example, a set of border points may not construct a manifold as shown in Figure 1 (b) and (c); they are called quasi-manifolds [START_REF] Kovalevsky | A topological method of surface representation[END_REF] and (b) is also called a pseudo-manifold [START_REF] Stillwell | Classical Topology and Combinatorial Group Theory[END_REF].

For border tracking in Z ! , therefore, the algebraic-topology based approach based on voxels is commonly used [START_REF] Artzy | The theory, design, implementation, and evaluation of a three-dimensional surface detection algorithm[END_REF] such as tracking the common faces between two voxels centered at the points p in V ⊂ Z ! and q in V = Z ! \ V. Such faces are represented by the ordered pair (p, q). Since q ∈ N $ (p) ∩ V, the set of all p of such pairs (p, q) becomes equal to Br $ (V) of ( 8). In this approach, the surface is represented by a set of square faces of voxels and the topological structures of cellular complexes, i.e. voxels, voxel faces, etc., are shown in [START_REF] Kovalevsky | Finite topology as applied to image analyses[END_REF].

Because we would like to consider triangulated surfaces on the points of Br m (V), we need another notion based on algebraic topology. In this paper, we extend our notions of discrete simplexes in [START_REF] Kenmochi | Boundary extraction of discrete objects[END_REF] to discrete convex polyhedra and give the de nition of discrete polyhedral complexes instead of discrete simplicial complexes in [START_REF] Kenmochi | Boundary extraction of discrete objects[END_REF]. The following sections are devoted for presenting triangulation of Br m (V).

Discrete Polyhedral Complexes and Combinatorial Boundaries

In this section, we give de nitions of a polyhedral complex which consists of a nite set of convex polyhedra such that the vertices are all points in Z ! and any adjacent vertices are m-neighboring. Given a lattice point set V ⊂ Z ! , such a polyhedral complex is introduced for giving a complicial representation, i.e. an object representation by a complex, of V. An algorithm for obtaining a polyhedral complex from V will be presented in the next section. Similar complicial representations for V are also found, for examples, in [START_REF] Klette | m-Dimensional cellular spaces[END_REF][START_REF] Kong | Concepts of digital topology[END_REF][START_REF] Tourlakis | Some results in computational topology[END_REF]. The differences between our complicial representation and them will be discussed in Section 6.

Convex Polyhedra and Polyhedral Complexes in R n

For the de nitions of convex polyhedra and polyhedral complexes in R n , we follow the notions in [START_REF] Ziegler | Lectures on polytopes[END_REF]. # Similar notations are also seen in [START_REF] Alexandrov | Combinatorial topology I[END_REF][START_REF] Stillwell | Classical Topology and Combinatorial Group Theory[END_REF].

De nition 1 A convex polyhedron I is the convex hull of a nite set of points in some R d .

The dimension of a convex polyhedron I is the dimension of its af ne hull. An n-dimensional convex polyhedron I is abbreviated to an n-polyhedron. For in-# Instead of the term convex polyhedra, polytopes is used in [START_REF] Ziegler | Lectures on polytopes[END_REF] stance, a point is a 0-polyhedron, a line segment is a 1-polyhedron, a triangle is a 2-polyhedron, and a tetrahedron is a 3-polyhedron. A linear inequality a • x ≤ z is valid for I if it is satis ed for all points x ∈ I. A face of I is then de ned by any set of the form @ = I ∩ {x ∈ R @ : a • x = z} where a • x ≤ z is valid for I. For instance, a 3-polyhedron which is a tetrahedron has four 0-polyhedra, six 1-polyhedra and four 2-polyhedra for its faces. The point of a 0-polyhedron, the endpoints of a 1-polyhedron and the vertices of 2-and 3polyhedra are called the vertices of each convex polyhedron. De nition 2 A polyhedral complex K is a nite collection of convex polyhedra such that (1) the empty polyhedron is in K, (2) if I ∈ K, then all faces of I are also in K, (3) the intersection I ∩ J of two convex polyhedra I, J ∈ K is a face both of I and of J.

The dimension of K is the largest dimension of a convex polyhedron in K.

Note that any K is a partially ordered set which can be identi ed with a topological space called a discrete space; the proof is nd in Section 6.1 of [START_REF] Alexandrov | Combinatorial topology I[END_REF].

Discrete Convex Polyhedra and Discrete Polyhedral Complexes

Now we consider polyhedral complexes such that the vertices of convex polyhedra are all lattice points in Z ! and the adjacent vertices are m-neighboring for m = 6,18,26. For constructing such polyhedral complexes, we rst consider all possible convex polyhedra such that all vertices are lattice points and any adjacent vertices of a convex polyhedron are m-neighboring each other for m = 6,18,26. Such convex polyhedra and polyhedral complexes are called discrete convex polyhedra and discrete polyhedral complexes hereafter. The constraints allow us to look for a discrete convex polyhedron which is not larger than the unit cubic region as follows.

Let us consider all possible convex polyhedra in a unit cubic region such that the vertices of each convex polyhedron are the vertices of a unit cube. A unit cube has eight lattice points for the vertices. For each lattice point we assign the value of either 1 or 0 and call the point a 1-or 0-point, respectively. There are 256 con gurations of 1-and 0-points for the eight lattice points in a unit cubic region. In fact, we can reduce the number of the con gurations from 256 to 23 with considering the congruent con gurations by rotations as shown in Table 1. $ For each con guration, we obtain a convex polyhedron such that the vertices of the polyhedron are 1-points. We then classify each convex polyhedron into a set of discrete convex polyhedra with the dimension of n = 0,1,2,3 and with the mneighborhood relations between the adjacent vertices for m = 6,18,26 as shown in Table 2. From Table 2, we see that there are a nite number of discrete convex polyhedra for each neighborhood system and for each dimension from 0 to 3. For the abbreviation, we call the n-dimensional discrete convex polyhedra in Table 2 discrete n-polyhedra hereafter. For any neighborhood system, an isolated point of con guration P1 in Table 2 is regarded as a discrete 0-polyhedron. Similarly, the line segment for con guration P2a is regarded as a discrete 1-polyhedron for any neighborhood system because the adjacency between two points are m-neighboring for any m = 6,18,26. However, the line segment of con guration P2b is not considered to be a discrete 1-polyhedron for the 6-neighborhood system, but considered to be a discrete 1-polyhedron for the 18-and 26-neighborhood systems. The line segment of con guration P2c is considered to be a discrete 1-polyhedron only for the 26-neighborhood system. Table 2 illustrates that we have one, two and three of discrete 1-polyhedra for the 6-, 18-and 26-neighborhood systems, respectively. A discrete 2-polyhedron is always bounded by discrete 1-polyhedra which are the faces of the discrete 2-polyhedron. Therefore, all discrete 2-polyhedra for the 6-neighborhood system have the point con guration of P4a. For the 18-and 26-neighborhood systems, four and ve different discrete 2-polyhedra are considered, respectively. In a similar way, a discrete 3-polyhedron is bounded by discrete 2-polyhedra which are the faces of the discrete 3-polyhedron. The discrete 3-polyhedra for each neighborhood system are illustrated in the last line of Table 2.

In Table 2, we see that for each m-neighborhood system, m = 6,18,26, every n ′dimensional face of any discrete n-polyhedron for n ′ < n is also a discrete n ′polyhedron. This is important because it enables us to construct a discrete polyhedral complex which is a nite collection of discrete convex polyhedra satisfying the three conditions in De nition 2 for each m-neighborhood system. Hereafter, we call an n-dimensional discrete polyhedral complex, for short, a discrete n-complex.

If we cannot decompose a discrete n-polyhedron into other discrete n-polyhedra in one of the neighborhood systems, such a discrete n-polyhedron also called a discrete n-simplex [START_REF] Kenmochi | Boundary extraction of discrete objects[END_REF]. In R ! , any n-dimensional simplex has n + 1 vertices [START_REF] Alexandrov | Combinatorial topology I[END_REF] while there exist, in Z ! , discrete n-simplexes which has more than n + 1 vertices such as the discrete simplexes of P4a and P8 for the 6-neighborhood system in Table 2. In mathematics such as combinatorial topology, simplexes are sometimes more $ For the proof that the 23 con gurations are complete, see in the appendix B of [START_REF] Kong | A Theory of Binary Digital Pictures[END_REF]. Note that there are 22 con gurations in [START_REF] Kong | A Theory of Binary Digital Pictures[END_REF] because symmetry is also considered. focused on than cells or convex polyhedra. It is because polygonal 2-polyhedra are too general compared with triangular 2-simplexes. In our case, however, if we only use discrete simplexes for triangulation of a subset V of Z ! , the simplicial decomposition of V may not be accomplished for 18-neighborhood system even if it is accomplished for 6-and 26-neighborhood systems [START_REF] Kenmochi | Boundary extraction of discrete objects[END_REF]. In this paper, therefore, we show that triangulation of V is succeeded for any neighborhood system by using not only discrete simplexes but also discrete convex polyhedra in Table 2.

Combinatorial Boundaries as Discrete Polyhedral Complexes

Before de ning combinatorial boundaries, we give some topological notions for discrete polyhedral complexes [START_REF] Alexandrov | Combinatorial topology I[END_REF]. A discrete n-complex K is said to be pure if every discrete n ′ -polyhedron of K where n ′ < n is a face of some discrete npolyhedron. Figure 4 illustrates examples of pure and non-pure discrete 3-polyhedra for the 26-neighborhood system. If K is any subset of K, the complex consisting of all the elements of K and of all the elements of K each of which is a face of at least one element of K is called the combinatorial closure Cl(K ) of K in K.

We consider a discrete polyhedral complex C as a topological representation of V ⊂ Z ! , i.e., as a topological space by topologizing V; note that we topologize V but not the whole space of Z ! . Because we require our boundary representation to contain the surface structures such as triangulated surfaces, we consider a pure discrete 3-subcomplex O ⊆ C and de ne the boundary ¶O of O for the combinatorial boundary of V; a procedure for obtaining ¶O from V will be presented in the next section. The notion of such combinatorial boundary ¶O is based on algebraic topology [START_REF] Stillwell | Classical Topology and Combinatorial Group Theory[END_REF].

De nition 3 Let O be a pure discrete 3-complex and H be the set of all discrete 2-polyhedra in O each of which is a face of exactly one discrete 3-polyhedron in O. The boundary of O is de ned as ¶O = Cl(H).

From De nition 3, we obtain the following proposition. Proposition 4 The boundary ¶O of a pure discrete 3-complex O is a pure discrete 2-subcomplex of O.

Note that the union of all discrete convex polyhedra in ¶O may not form a manifold but form a non-manifold such as a pseudo-manifold [START_REF] Stillwell | Classical Topology and Combinatorial Group Theory[END_REF] and a quasi-manifold [START_REF] Kovalevsky | A topological method of surface representation[END_REF] as shown in Figure 1 according to the de nition.

Because discrete convex polyhedra are de ned for each m-neighborhood system where m = 6,18,26, a discrete polyhedral complex C, a discrete pure 3-polyhedron O and the combinatorial boundary ¶O are also de ned for each m-neighborhood system. When we insist a m-neighborhood system considering for them, they are denoted by C m , O m and ¶O m instead.

Combinatorial Boundary Extraction

This section presents a procedure for obtaining the boundary ¶O m of a pure discrete 3-complex O m from any nite set V ⊂ Z ! for each m = 6,18,26. The procedure is divided into three steps as shown in Figure 5: (I) decompose V into discrete n-polyhedra where n = 0,1,2,3 such that those discrete n-polyhedra constitutes a discrete polyhedral complex C m (from (a) to (b) in Figure 5), (II) from C m , obtain a pure discrete 3-subcomplex O m ⊆ C m (from (b) to (c) in Figure 5), and (III) extract the boundary ¶O m of O m for each m-neighborhood system, m = 6,18,26 (from (c) to (d) in Figure 5). After explaining such a procedure in subsection 4.1, we also presents a practical algorithm for obtaining ¶O m directly from V in subsection 4.2.

4.1

Theory of the Procedure 4.1.1

Step 1: Polyhedral Decomposition of V

The decomposition of V into discrete convex polyhedra is achieved in two steps.

For each x = (i, j, k) in Z ! , let D(x) = {(i + e , j + e , k + e ! ) | e i = 0 or 1}.

We say that the points of V are 1-points and the points of Z ! \ V are 0-points. For each x ∈ Z ! , we locally consider a discrete polyhedral complex C m (x) as follows.

If a discrete n-polyhedron s for an m-neighborhood system exists with respect to an con guration of 1-points in D(x) in Table 2, we set C m (x) to be a collection of s and its faces where n = 0,1,2,3. Otherwise, we consider discrete n-polyhedra s such that n is as large as possible where n ≤ 3 and the vertices of s are all 1-points in D(x) and set C m (x) to be a collection of such ss and their faces. For each 1-point con guration of D(x), we then obtain a discrete polyhedral complex C m (x) for each m = 6,18,26 as shown in Tables 3, 4 and 5, respectively.

Now let

C m = ∪ x∈Z ! C m (x) (9)
and we verify that + m is mostly a discrete polyhedral complex satisfying the conditions in De nition 2; there are few exceptional cases that we need to replace + m (x)s in ( 9) to obtain a discrete polyhedral complex + m only for m = 18.

Say that + m (x) and + m (y) are adjacent if ,(x) ∩ ,(y) = / 0. Their adjacency types are classi ed into the following three #(,(x) ∩ ,(y)) = 1,2 or 4 (and never 3) where #(A) represents the number of elements of the set A. The adjacency types and the conceivable polyhedral decomposition at the joint are illustrated in Table 6.

For each adjacent pair of + m (x) and + m (y), let + m (x, y) = + m (x) ∪ + m (y).

(

) 10
We then verify, from Tables 3, 4 and 5, that + m (x, y) is mostly a discrete polyhedral complex satisfying the conditions of De nition 2; only for m = 18, there are two exceptional cases that we need to replace both + m (x) and + m (y) in ( 10) to obtain a discrete polyhedral complex + m (x, y).

First, let us consider the case of #(,(x) ∩ ,(y)) = 1. As shown in the rst line of Table 6, the common point z is either 1-or 0-point. If z is a 0-point (Case 1), + m (x) and + m (y) include no common discrete convex polyhedron. Thus, we simply obtain + m (x, y) by [START_REF] Herman | Discrete multidimensional Jordan surfaces[END_REF] as empty. If z is a 1-point (Case 2), both + m (x) and + m (y) include a common discrete 0-polyhedron s . Let us introduce the notion of the skeleton Sk(s) of a discrete convex polyhedron s such as the set of all vertices of s [START_REF] Alexandrov | Combinatorial topology I[END_REF]. Then, we have Sk(s ) = {z}. Thus, we obtain a discrete 0-complex + m (x, y) = {s } by [START_REF] Herman | Discrete multidimensional Jordan surfaces[END_REF].

In the case of #(,(x) ∩ ,(y)) = 2, there are two common points z and z as shown in the second line of Table 6. Since each of z and z is either 1-or 0-point, there are three possible con gurations of 1-and 0-points for the pair of z and z . If both z and z are 0-points (Case 1), both + m (x) and + m (y) include no common discrete polyhedron. Thus, we obtain + m (x, y) by [START_REF] Herman | Discrete multidimensional Jordan surfaces[END_REF] as empty. If either of z and z is 1-point (Case 2), the 1-point becomes the common 0-polyhedron s in + m (x) and + m (y). Thus, we obtain a discrete 0-complex + m (x, y) = {s} by [START_REF] Herman | Discrete multidimensional Jordan surfaces[END_REF]. If both z and z are 1-points (Case 3), + m (x) and + m (y) have a discrete 1-polyhedron s and its 0-dimensional faces as the common discrete polyhedra such that Sk(s) = {z , z }.

Thus, we obtain a discrete 1-complex + m (x, y) = Cl({s}) by [START_REF] Herman | Discrete multidimensional Jordan surfaces[END_REF].

In the case of #(,(x) ∩ ,(y)) = 4, let z i for i = 1,2,3,4 be the four common points.

Since each point is either 1-or 0-point, there are six con gurations of 1-and 0points for the four points as shown in the last line of Table 6. It also shows the possible common discrete polyhedra of + m (x) and + m (y) for each con guration;

• the empty set in Case 1;

• a discrete 0-polyhedron I such that Sk(I ) = {z } in Case 2;

• a discrete 1-polyhedron I and its 0-dimensional faces such that Sk(I ) = {z , z ! } in Case 3;

• two discrete 0-polyhedra I and I such that Sk(I ) = {z } and Sk(I ) = {z " } for the 6-neighborhood system and a discrete 1-polyhedron I with its 0dimensional faces such that Sk(I ) = {z , z " } for the 18-and 26-neighborhood systems in Case 4;

• two discrete 1-polyhedra I and I with their faces such that Sk(I ) = {z , z ! } and Sk(I ) = {z ! , z " } for the 6-neighborhood system, and a discrete 2-polyhedron I with its faces such that Sk(I ) = {z , z ! , z " } for the 18-and 26-neighborhood systems in Case 5;

• in Case 6, a discrete 2-polyhedron I and its faces such that Sk(I) = {z , z , z ! , z " } for any neighborhood system except for the cases illustrated in Figures 6 (a) and (c) which appear only when m = 18.

We now discuss the exceptional cases which appear only for the 18-neighborhood system as illustrated in Figures 6 (a) and (c). In Figure 6 (a), both adjacent unit cubes ,(x) and ,(y) have the con gurations P5a. In such case, + & (x) includes two discrete 2-polyhedra I and I such that Sk(I ) = {z , z , z ! } and Sk(I ) = {z , z " , z ! }, and + & (y) includes two discrete 2-polyhedra I and I ! such that Sk(I ) = {z , z , z " } and Sk(I ! ) = {z , z ! , z " }. Thus, if + & (x, y) is obtained by [START_REF] Herman | Discrete multidimensional Jordan surfaces[END_REF], such + & (x, y) does not become a discrete polyhedral complex. In order to construct + & (x, y) which is a discrete polyhedral complex, we therefore replace discrete polyhedral complexes + & (x) and + & (y) from Figure 6 (a) to (b) so that

+ & (x) = Cl({J }) ∪Cl({J }), + & (y) = Cl({J }) ∪Cl({J ! })
where Sk(J ) = {z , z , z # }, Sk(J ) = {z , z ! , z # }, Sk(J ) = {z ! , z , z $ }, and Sk(J ! ) = {z " , z ! , z $ }. We then obtain

+ & (x, y) = + & (x) ∪ + & (y) = Cl(J ) ∪Cl(J ) ∪Cl(J )
where Sk(J ) = {z , z }, Sk(J ) = {z , z ! } and Sk(J ) = {z ! , z " }.

In Figure 6 (c), adjacent unit cubes ,(x) and ,(y) have the pair of con gurations P5a and P4a. In such case, + & (x) includes two discrete 2-polyhedra I and I such that Sk(I ) = {z , z , z ! } and Sk(I ) = {z , z " , z ! }, and + & (y) includes a discrete 2-polyhedron I such that Sk(I ) = {z , z , z ! , z " }. Thus, if + & (x, y) is obtained by [START_REF] Herman | Discrete multidimensional Jordan surfaces[END_REF], such + & (x, y) does not become a discrete polyhedral complex. In order to construct + & (x, y) which is a discrete polyhedral complex, we therefore replace discrete polyhedral complexes + & (x) and + & (y) from Figure 6 (c) to (d)

so that C & (x) = Cl({J ! }) ∪Cl({I }), C & (y) = Cl({I }) ∪Cl({I })
where Sk(J ! ) = {z , z , z ! , z # }. We then obtain

C & (x, y) = C & (x) ∪ C & (y) = C & (y) = Cl(I ) ∪Cl(I ).
Consequently, setting C m (x) for each x ∈ Z ! referring to Tables 3 We can also obtain O m directly from V without considering C m such that

O m = ∪ x∈Z ! O m (x)
where O m (x) is a pure discrete 3-complex at each unit cubic region D(x). Each O m (x) is easily obtained by referring to one of Tables 7, 8 and 9 instead of one of Tables 3, 4 and 5 for C m (x). We easily create Tables 7, 8 and 9 by making C m (x) in Tables 3, 4 and 5 to be pure. Note that O & (x) will be replaced as an empty set if the 1-point con gurations at D(x) and its adjacent D(y) are as illustrated in Figure 6 (a). From the above procedure, we consequently obtain the following proposition.

Proposition 5 Given a nite subset V ⊂ Z ! , the combinatorial boundary ¶O m is uniquely obtained for any m-neighborhood system, m = 6,18,26.

Algorithm of Combinatorial Boundary Extraction

For practical use, we present an effective algorithm of generating ¶O m directly from V by referring to Table 10, which is a similar table used for the marching cubes method [START_REF] Lorensen | Marching cubes: a high-resolution 3d surface construction algorithm[END_REF][START_REF] Wyvill | Data structure for soft objects[END_REF], for each neighborhood system. The comparison between the marching cubes method and our method is discussed in [START_REF] Kenmochi | Marching cubes method with connectivity[END_REF].

We obtain Table 10 from Tables 3, 4 and 5 as follows. First we look only for discrete 2-polyhedra of C m (x) at each unit cubic region D m (x) because ¶O m is a pure discrete 2-complex; ¶O m does not contain more than three-dimensional discrete convex polyhedra and less than two-dimensional discrete convex polyhedra which are not faces of any discrete 2-polyhedra. We then classify each discrete 2-polyhedron s of C m (x) in Tables 3, 4 and 5 into the following four types, A-1, A-2, B-1, B-2:

A. s is a face of a discrete 3-polyhedron d ∈ C m (x) so that s = d ∩ {x ∈ R ! : a • x = z} ( 12
)
where a • x ≤ z is valid for d, and A-1. s is located at a face of a unit cube D(x) so that

{x ∈ R ! : a • x ≤ z} ∩ D(x) = D(x); (13)
A-2. s is located inside a unit cube D(x) so that (13) does not hold; B. s is not a face of any discrete 3-polyhedron d ∈ C m (x), i.e., there is no discrete 3-polyhedron d which satis es [START_REF] Kenmochi | Marching cubes method with connectivity[END_REF], and B-1. s is located at a face of a unit cube D(x), so that either of the equations [START_REF] Klette | Digital Topology for Image Analysis, Part I; Basics and Planar Image Carriers[END_REF] holds where s = s ∩ {x ∈ R ! : a • x = z}; B-2. s is located inside a unit cube D(x) so that neither ( 14) nor (15) holds.

{x ∈ R ! : a • x ≤ z} ∩ D(x) = D(x), ( 14
) {x ∈ R ! : a • x ≥ z} ∩ D(x) = D(x)
For each C m (x), we obtain the set of discrete 2-polyhedra of each type, A-1, A-2, B-1, B-2, which is denoted by TP (x), TP (x), TP (x) or TP (x), respectively. In the following, we show that only two types TP (x) and TP (x) can be found in ¶O m .

For s ∈ TP (x), s does not belong to ¶O m if there exits a discrete 3-polyhedron d at a unit cube D(y) adjacent to D(x) such that s is a face of d. If there is no such d at D(y), then s ∈ TP (y). For s ∈ TP (x), s ∈ C m \ O m and thus s /

∈ ¶O m . Setting I m (x) = ∪ I∈62 (x) Cl(s), (16)
J m (x) = ∪ I∈62 (x)
Cl(s), [START_REF] Kong | Continuous Analogs of Axiomatized Digital Surfaces[END_REF] in Table 10 we obtain a pure discrete 2-complex T m (x) = J m (x) ∪ I m (x). [START_REF] Kong | Concepts of digital topology[END_REF] for every x ∈ Z ! . The arrow of every s in Table 10 indicates the side where the half space {x ∈ R ! : a • x > z} exists; roughly speaking, it is oriented to the exterior of ¶O m and is useful for visualization as a normal vector of each s. Note that either J m (x) or I m (x) is empty for any T m (x) in Table 10 except for the con guration P5a of the 18-neighborhood system.

We then see that any t ∈ I m (x) constitutes ¶O m , thus,

I m (x) ⊂ ¶O m . ( 19
)
For a discrete 2-polyhedron s ∈ J m (x), if s ∈ J m (y) at an adjacent unit cube D(y) to D(x) as shown in Figure 7, then

s ∈ C m \ O m , (20)
and otherwise s ∈ ¶O m . [START_REF] Kovalevsky | A topological method of surface representation[END_REF] Therefore, we need to verify [START_REF] Kovalevsky | A topological method of surface representation[END_REF] (or ( 20)) for each s ∈ J m (x) for constructing ¶O m , while every s ∈ I m (x) is always in ¶O m from [START_REF] Kong | Digital topology: introduction and survey[END_REF]. Such veri cation is achieved in step 2.3 in Algorithm 1. The special treatment for the cases illustrated in Figures 6 (a) and (c) which occur only for the 18-neighborhood system is also considered in step 2.2 in Algorithm 1. For the algorithm, we set an input lattice space to be nite such as 10; 2.2 if m = 18, check if each pair of T m (x) and T m (y) for y = (i-1, j, k), (i, j -1,k),(i, j, k -1) where x = (i, j, k) is in the case as illustrated in Figure 6 (a) or (c); if so, replace T m (x) and T m (y) from Figure 6 (a) to (b), or (c) to (d);

W = {(i, j, k) ∈ Z ! : 1 ≤ i ≤ L, 1 ≤ j ≤ M, 1 ≤ k ≤ N}. Algorithm 1 input: A subset V of W. output: A combinatorial boundary ¶O m .
2.3 if J m (x) of T m (x) is not empty, then check for each s ∈ J m (x) if there exists y such that s ∈ J m (y) where y = (i -1, j, k), (i, j -

as shown in Figure 7; if so, replace T m (x) and T m (y) with Cl(T m (x) \ Cl({s})) and Cl(T m (y) \Cl({s}));

3 obtain ¶O m = ∪ x∈9 T m (x).
end From Algorithm 1, it is obvious that we obtain the combinatorial boundary ¶O m for each m = 6,18,26 from any nite set V ⊂ Z ! . Some experimental results of combinatorial boundaries ¶O m for the various inputs V, such as digitized sphere, cube, torus and catenoid, with respect to m = 6,18,26 are shown in Figures 8,9, 10 and 11. Those inputs of volume data are made by Volgen [30].

Relations between Borders and Combinatorial Boundaries

We already introduced the notion of the skeleton Sk(s) of a discrete convex polyhedron s such as the set of the vertices of s in the previous section [START_REF] Alexandrov | Combinatorial topology I[END_REF]. Let ¶O m be the combinatorial boundary obtained by Algorithm 1 from a given V ⊂ Z ! for m = 6, 18, 26. We call the union of the skeletons of all discrete convex polyhedra of ¶O m the skeleton of ¶O m and it is denoted by Sk( ¶O m ). We then obtain the following relations between the skeleton Sk( ¶O m ) and the border Br m ′ (V) of ( 8). Those relations are considered to be the discrete version of the relation ( 5) in R ! . Theorem 6 The border Br m ′ (V) and the skeleton Sk( ¶O m ) of the combinatorial boundary ¶O m obtained from a nite subset V ⊂ Z ! have the relations such that

Br $ (V) = Sk( ¶O $ ) ∪ (Sk(C $ ) \ Sk(O $ )) (22) = Sk( ¶O & ) ∪ (Sk(C & ) \ Sk(O & )) \ A ($,&) , (23)
Br & (V) = Sk( ¶O $ ) ∪ (Sk(C $ ) \ Sk(O $ )) \ A (&,$) , (24)
Br $ (V) = Sk( ¶O $ ) ∪ (Sk(C $ ) \ Sk(O $ )).

(

) 25
where

A (m ′ ,m) = ∪ x∈Z ! A (m ′ ,m) (x)
so that A (m ′ ,m) (x) is shown in the right column of Table 11 as the set of black points at a unit cube D(x) only when D(x) has a 1-point con guration P5a or P7 only for (m ′ , m) = [START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF][START_REF] Kong | Concepts of digital topology[END_REF] or [START_REF] Kong | Concepts of digital topology[END_REF][START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF], respectively. Note that A ($,&) (x) for the con guration P5a is empty if it has no adjacent unit cube whose con guration is also P5a as shown in Figure 6 (b).

A pair (m ′ , m) of neighborhood systems which is considered in Theorem 6 is [START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF][START_REF] Sloboda | On the topology of grid continua[END_REF], [START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF][START_REF] Kong | Concepts of digital topology[END_REF], [START_REF] Kong | Concepts of digital topology[END_REF][START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF] or [START_REF] Sloboda | On the topology of grid continua[END_REF][START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF] and similar pairs (m ′ , m) are also seen in the relations between Br ′ m (V) and its m-connectivity mentioned in section 2.2.1 [16,18].

From Theorem 6, we can derive the following corollary which has more similar formulas to (5) for R n shown in [START_REF] Moise | Geometric Topology in Dimensions 2 and 3[END_REF].

Corollary 7 For a pure discrete 3-complex O m where m = 6, 18, 26, we have the relations such that

Br $ (Sk(O $ )) = Sk( ¶O $ ), Br $ (Sk(O & )) = Sk( ¶O & ) \ A ($,&) , Br & (Sk(O $ )) = Sk( ¶O $ ) \ A (&,$) , Br $ (Sk(O $ )) = Sk( ¶O $ ).
To emphasize the difference between the continuous and discrete cases, however, we refer to the relations ( 22), ( 23), ( 24) and ( 25) of Theorem 6 rather than those of Corollary 7. The difference is that there is an additional term which is the second term for the union in the right side of each equation of ( 22), ( 23), ( 24) and ( 25) while there is no such additional term in [START_REF] Clifford | The postulates of the science of space[END_REF]. The second term Sk(C m ) \ Sk(O m ) is a set of vertices which are not included in any discrete 3-polyhedra but included in less than three-dimensional discrete convex polyhedra of C m . In Figure 5, we can see Sk(C m ) \ Sk(O m ) such as the right point which is included in a discrete 1polyhedron in (b) but does not appear in (c). Because we have a step for removing less than three-dimensional parts from three-dimensional parts, namely, obtaining O m from C m such as a procedure from (b) to (c) in Figure 5, we need to have the second term Sk(C m ) \ Sk(O m ) to compare with the border points Br m ′ (V) based on general topology. Note that not only discrete 1-polyhedra but also discrete 2and 0-polyhedra may exists in C m \ O m . Therefore, for the 3-dimensional border extraction, we take account of three kinds of dimensions for removing the parts C m \ O m whose dimension are reduced, i.e. zero, one and two dimensions, while we have to take account of two kinds of dimensions, i.e. zero and one dimension, for the 2-dimensional border extraction. This difference causes the dif culty of 3dimensional border tracking problem as we already mentioned in Subsection 2.2.1.

The third terms A ($,&) and A (&,$) which only appear in ( 23) and ( 24) respectively show the difference between Sk( ¶O & ) and Sk( ¶O $ ) of ( 23) and ( 22) and the difference between Br & (V) and Br $ (V) of ( 24) and ( 25), respectively. Note that A ($,&) rarely becomes non-empty; A ($,&) is not empty only if we have a pair of adjacent unit cubes whose 1-point con gurations are both P5a as shown in Figure 6 (b).

For proving Theorem 6, we need the following two lemmas.

Lemma 8 For a unit cubic region D(x), setting

CubeBr m ′ (V; x) = {y ∈ V ∩ D(x) : D(x) ∩ N m ′ (y) ∩ V = / 0} (26)
for each m ′ = 6,18,26, we have

Br m ′ (V) = ∪ x∈Z ! CubeBr m ′ (V; x). ( 27
) Proof. Since N m ′ (y) ∩ V = ∪ x∈Z ! (D(x) ∩ N m ′ (y) ∩ V),
we obtain ( 27) from ( 8) and ( 26). 2

The points in CubeBr m ′ (V; x) are illustrated for every possible con guration of 1points in D(x) in Table 11.

Lemma 9 At each unit cubic region D(x) for x ∈ Z ! , setting T m (x) to be a discrete 2-complex given by Table 10, C m (x) to be a discrete polyhedral complex given by Tables 3, 4 and5, and O m (x) to be a pure discrete 3-complex of C m (x) by Tables 7,8 and 9, we have

Sk( ¶O m ) ∪ (Sk(C m ) \ Sk(O m )) = ∪ x∈Z ! (Sk(T m (x)) ∪ (Sk(C m (x)) \ Sk(O m (x)) \ Sk(T m (x)))). (28)
Proof.

(I) First we show the following inclusion:

∪ x∈Z ! (Sk(T m (x)) ∪ (Sk(C m (x)) \ Sk(O m (x)) \ Sk(T m (x)))) ⊂ Sk( ¶O m ) ∪ (Sk(C m ) \ Sk(O m )). (29)
Let us consider the two pure discrete 2-complexes J m (x) and I m (x) of ( 16) and ( 17) for each T m (x) of [START_REF] Kong | Concepts of digital topology[END_REF]. From [START_REF] Kong | Concepts of digital topology[END_REF], we then derive the relation Sk(T m (x)) = Sk(I m (x)) ∪ Sk(J m (x)), (30) and from [START_REF] Kong | Digital topology: introduction and survey[END_REF],

Sk(I m (x)) ⊂ Sk( ¶O m ). (31)
The point con gurations of Sk(I m (x)) are illustrated in Table 11. Let us consider a vertex z ∈ Sk(J m (x)). Because J m (x) is a pure discrete 2-complex, any z is included in a discrete 2-polyhedron s ∈ J m (x) and such a s satis es either [START_REF] Kovalevsky | Finite topology as applied to image analyses[END_REF] or [START_REF] Kovalevsky | A topological method of surface representation[END_REF]. If z is a vertex of s of (20),

z ∈ Sk(C m ) \ Sk(O m ) (32)
and if z is a vertex of s of ( 21),

z ∈ Sk( ¶O m ). (33)
Thus,

Sk(J m (x)) ⊂ Sk( ¶O m ) ∪ (Sk(C m ) \ Sk(O m )) (34) Therefore, if z ∈ Sk( ¶O m ) ∪ (Sk(C m ) \ Sk(O m )), then z ∈ Sk(T m (x) ∪ (Sk(C m (x)) \ Sk(O m (x)) \ Sk(T m (x)))
and it contradicts (37).

From (I) and (II), we thus obtain [START_REF] Stillwell | Classical Topology and Combinatorial Group Theory[END_REF]. 2

Proof of Theorem 6. For (m ′ , m) = (6, 26), [START_REF] Sloboda | On the topology of grid continua[END_REF][START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF], we have

CubeBr m ′ (V; x) = Sk(T m (x)) ∪ (Sk(C m (x)) \ Sk(O m (x)) \ Sk(T m (x))) (38)
for any x ∈ Z ! from (30) and Table 11. Thus, from Lemmas 8 and 9, we obtain [START_REF] Lorensen | Marching cubes: a high-resolution 3d surface construction algorithm[END_REF] and [START_REF] Rosenfeld | Connectivity in digital pictures[END_REF].

For (m ′ , m) = (6, 18), if we have the case as shown in Figure 6 (b) for the con guration P5a of D(x), we see that

CubeBr m ′ (V; x) = Sk(T m (x)) ∪ (Sk(C m (x)) \ Sk(O m (x)) \ Sk(T m (x))) \ A ($,&) (x)
from Table 11, and otherwise we have (38). Thus, we obtain [START_REF] Moise | Geometric Topology in Dimensions 2 and 3[END_REF].

For (m ′ , m) = (18, 6), for the con guration P7 of D(x), we see that

CubeBr m ′ (V; x) = Sk(T m (x)) ∪ (Sk(C m (x)) \ Sk(O m (x)) \ Sk(T m (x))) \ A (&,$) (x)
from Table 11, and otherwise we have (38). Thus, we obtain (24). 2

Conclusions and Discussions

In this paper, we gave a solution to one of the important problems in three-dimensional image analysis; is it possible to give a triangulation of border points Br m (V) such that all vertices of triangulated surfaces are border points and adjacent vertices are m-neighboring for m = 6,18,26? Our answer is yes. We also succeed to present Algorithm 1 which gives such a triangulated surface ¶O m from any nite subset V ⊂ Z ! . We insists that the calculation time is linear to the size of V, i.e. the size of a 3-dimensional digital image, and it is the same as that of the set operation [START_REF] Franc ¸on | Discrete combinatorial surfaces[END_REF] for obtaining Br m (V) from V even if our algorithm provides a pure discrete 2polyhedron ¶O m which contains not only a point set Sk( ¶O m ) but also the combinatorial topological structures of ¶O m . Theorem 6 which indicates discrete versions of the relation [START_REF] Clifford | The postulates of the science of space[END_REF] shows that ¶O m ′ becomes a triangulation of Br m (V) if we choose a good pair such as (m, m ′ ) = (6, 18), [START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF][START_REF] Sloboda | On the topology of grid continua[END_REF], [START_REF] Kong | Concepts of digital topology[END_REF][START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF], [START_REF] Sloboda | On the topology of grid continua[END_REF][START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF]. Note that there may be extra points of Sk(C m ′ ) \ Sk(O m ′ ) if Br m (V) contains some lattice points where we cannot put any discrete 3-polyhedron because of their con gurations such as the con guration around the right point of Figure 5 (b). Our discrete polyhedral complex is useful to analyse the reasons why we have to ignore such points, i.e. points of Sk(C m ′ ) \ Sk(O m ′ ) for triangulation of Br m (V). It is also interesting that the possible pairs for (m, m ′ ) are similar to the pairs (a, b) for b-connectedness of a-borders [START_REF] Kong | Concepts of digital topology[END_REF] as we mentioned in Section 2.2.1.

Improvement of the Combinatorial Boundary Tracking Algorithm

It may be also possible to present more effective combinatorial boundary tracking algorithms whose calculation time is linear to the number of border points if we succeed to investigate every possible local con gurations of combinatorial boundaries. In fact, such an effective border tracking algorithm for three-dimensional digital image is already presented by using an algebraic-topology based approach by using voxel faces [START_REF] Kovalevsky | A topological method of surface representation[END_REF] We only need to extend the algorithm for discrete polyhedral complexes instead of their cellular complexes.

Comparison with Other Polyhedral Complexes in Z n

We took the combinatorial/algebraic-topology based approach by using discrete polyhedral complexes for giving a solution to the triangulation problem. Due to the strong powers for topological problems in discrete spaces, similar complicial representations for a nite subset V ⊂ Z ! are also seen in different literatures [START_REF] Klette | m-Dimensional cellular spaces[END_REF][START_REF] Kong | Concepts of digital topology[END_REF][START_REF] Tourlakis | Some results in computational topology[END_REF], for example. For our term of discrete polyhedral complexes C m for V, they use the different terms: cellular complexes [START_REF] Klette | m-Dimensional cellular spaces[END_REF], continuous analogs [START_REF] Kong | Concepts of digital topology[END_REF] and polyhedra [START_REF] Tourlakis | Some results in computational topology[END_REF]. Because their aims are different, the ways of obtaining C m from V are also different.

Continuous analogs are presented for de ning a digital fundamental group whose concept is used for three-dimensional thinning. During three-dimensional thinning, they need to preserve a digital topology whose criteria are given by using the concepts of connectedness and of a digital fundamental group. For a digital fundamental group, they need to consider a region of interest and also its complement, and therefore consider topologies for the whole Z ! , not only for V ⊂ Z ! as we do in this paper. In [START_REF] Kong | Concepts of digital topology[END_REF], one example for a set of continuous analogs is presented. They are different from our discrete polyhedral complexes in the geometric sense;

for example, some continuous analogs may have augmented points which are not lattice points but centroids of lattice cubes as their vertices. On the other hand, if we consider discrete polyhedral complexes C m (V) and C m ′ ( V) choosing some pairs for (m, m ′ ) for V and V, then we do not know if they satisfy the conditions of continuous analogs or not. Because such discussion is beyond the subjects of this paper, we leave it for our future work.

Even if the aims in [START_REF] Klette | m-Dimensional cellular spaces[END_REF][START_REF] Tourlakis | Some results in computational topology[END_REF] are different from ours such as calculation of topological equivalence between two different subsets of Z ! [START_REF] Tourlakis | Some results in computational topology[END_REF], we see that cellular complexes [START_REF] Klette | m-Dimensional cellular spaces[END_REF] and polyhedra [START_REF] Tourlakis | Some results in computational topology[END_REF] are the same as our discrete polyhedral complexes C $ (V) for the 6-neighborhood system. This is because the shapes of discrete convex polyhedra for m = 6 such as cubes, squares, unit line segments, etc. can be seen in lattice grids and they are straightforward to topologize Z ! . In fact, if we topologize Z ! instead of V ⊂ Z ! in the same way of C $ (V), i.e. C $ (Z ! ), we see Khalimsky space [START_REF] Khalimsky | Pattern analysis of N-dimensional digital images[END_REF] which is well known in digital image analysis. In [START_REF] Klette | Digital Topology for Image Analysis, Part I; Basics and Planar Image Carriers[END_REF] it is also shown that Khalimsky space is homeomorphic to Kovalevsky's nite topology [START_REF] Kovalevsky | Finite topology as applied to image analyses[END_REF] for the case Z .

.ECKHAI !&

Figure 3 (

 3 b) and (c) show examples of the internal 4-boundary Br " (V) and the external 4-boundary Br " ( V), respectively.

, 4 and 5

 5 with taking account of the additional replacements of Figure6for m = 18, we uniquely obtain C m by[START_REF] Hausdorff | Set Theory[END_REF] for any m = 6, 18, 26 from any V ⊂ Z ! .4.1.2 Step 2: Construction of a Pure Discrete 3-ComplexAssume that the dimension of C m is three. Let G to be the set of all discrete 3polyhedra in C m . In order to obtain a pure discrete 3-complex O m from C m , we remove all discrete n-polyhedra which are not included in any discrete 3-polyhedrain C m for every n < 3, such that O m = Cl(G).(11)If C m is less than three dimensions, G is empty and thus O m is also empty. This occurs when C m contains only discrete 0-, 1-and 2-polyhedra and have no discrete 3-polyhedron. We consider that C m \ O m each of whose element has less than three dimensions is caused by the limited resolution of a digital image. If we would like to change the dimensions of elements in C m \ O m into three, we may need to increase the resolution of a digital image at C m \ O m . This is natural because our aim is to obtain surface structures from border points to calculate the shape information such as surface areas and curvatures. In order to obtain surface structures, apparently an isolated point is not adequate and we need to increase the image resolution to have more points around the isolated point. From (11), it is clear that O m is uniquely obtained from C m . Examples of the procedure for obtaining O $ from C $ are seen in Figures 4 (from (b) to (a)) and 5 (from (b) to (c)).

4. 1 . 3

 13 Step 3: Boundary Extraction of a 3D Pure Discrete Complex From De nition 3, the boundary ¶O m of O m is derived from the set H of discrete 2-polyhedra in O m each of which is a face of exactly one discrete 3-polyhedron in O m . Because H is uniquely obtained from O m , ¶O m is also uniquely obtained from O m .

begin 1

 1 set m = 6,18 or 26; 2 for each x ∈ W do 2.1 obtain T m (x) by referring to Table

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .Fig. 11 .

 1234567891011 Fig.1. Examples of manifolds and non-manifolds: (a) a manifold, (b) a pseudo-manifold[START_REF] Stillwell | Classical Topology and Combinatorial Group Theory[END_REF], and (c) a quasi-manifold[START_REF] Kovalevsky | A topological method of surface representation[END_REF]. Each central black point is a border point which cannot be a point of simplicity surfaces[START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF].

! For thinning of 8 ⊂ Z n , we consider simple points which we can remove without collapsing the criteria of digital topology[START_REF] Kong | Digital topology: introduction and survey[END_REF]. Obviously, simple points are related to border/boundary points of 8 ⊂ Z n .

" This holds if a digital picture space (Z ! , m ′ , m) is weakly normal, i.e. it has one of such pairs (m, m ′ ) from Proposition 7.4.1 in[START_REF] Kong | Concepts of digital topology[END_REF].

Acknowledgements A part of this work was done when the rst author stayed at Laboratoire A2SI, ES-IEE, Cité Descartes, B.P. 99, 93162 Noisy-Le-Grand Cedex, France, with the support of the JSPS Postdoctoral Fellowships for Research Abroad from 2000 to 2002. She thanks Gilles Bertrand and Michel Couprie at ESIEE for useful discussion on the subject. Parts of the experiments were made by Yusuke Oda at Department of Information Technology, Okayama University.

for any x ∈ Z ! . The point con gurations of Sk(J m (x)) are also shown in Table 11. From (30), [START_REF] Voss | Discrete Images, Objects, and Functions in Z ![END_REF] and (34), we then obtain Sk(T m (x)) ⊂ Sk( ¶O m ) ∪ (Sk(C m ) \ Sk(O m )).

(35)

For each discrete convex polyhedron s ∈ C m (x) \ O m (x) \ T m (x), if s ∈ T m (y) at other unit cube D(y) adjacent to D(x), we have [START_REF] Kovalevsky | Finite topology as applied to image analyses[END_REF] or [START_REF] Kovalevsky | A topological method of surface representation[END_REF], and otherwise we have [START_REF] Kovalevsky | Finite topology as applied to image analyses[END_REF]. Consequently, if a vertex z ∈ Sk(C m (x))\Sk(O m (x))\ Sk(T m (x)) is a vertex of s of [START_REF] Kovalevsky | A topological method of surface representation[END_REF], we have [START_REF] Ziegler | Lectures on polytopes[END_REF], and if z is a vertex of s of [START_REF] Kovalevsky | Finite topology as applied to image analyses[END_REF], we have [START_REF] Wyvill | Data structure for soft objects[END_REF]. Thus,

are also shown in Table 11.

From ( 35) and (36), we have

for every x ∈ Z ! , and thus we obtain the inclusion [START_REF] Tourlakis | Some results in computational topology[END_REF].

(II) Now we verify if there exists a point

for any x ∈ Z ! . Considering a point z ∈ V ∩ D(x) which satis es (37) for a point x ∈ Z ! , we see that

from Tables 3,4, 5, and Tables 10 and11, namely,

Tables Table 1 All possible 23 con gurations of 1-and 0-points for the eight lattice points in a unit cubic region. With considering the congruent con gurations by rotations, we obtain all 256 con gurations from them. All discrete n-polyhedra for n = 0, 1, 2, 3 such that all vertices are lattice points in Z ! and the adjacent vertices are m-neighboring for m = 6, 18, 26. Note that discrete n-polyhedra with asterisks are called discrete n-simplexes in the reference [START_REF] Kenmochi | Boundary extraction of discrete objects[END_REF].

discrete convex polyhedra Discrete convex polyhedral decomposition + $ (N) with respect to every 1-point con guration of a unit cubic region ,(N). Discrete convex polyhedral decomposition + & (N) with respect to every 1-point con guration of a unit cubic region ,(N). Discrete convex polyhedral decomposition + $ (N) with respect to every 1-point con guration of a unit cubic region ,(N). Three adjacency types of two unit cubic regions ,(N) and ,(O) such that #(,(N) ∩ ,(O)) = 1, 2 and 4. For each adjacency type, all possible con gurations of 1-and 0-points and a discrete convex polyhedral decomposition are shown. "" Table 7 Three-dimensional polyhedral decomposition O $ (N) corresponding to the con guration of 1-points in a unit cubic region D(N). Three-dimensional polyhedral decomposition O & (N) corresponding to the con guration of 1-points in a unit cubic region D(N). Three-dimensional polyhedral decomposition O $ (N) corresponding to the con guration of 1-points in a unit cubic region D(N). The look-up table which provides a one-to-one correspondence between an con guration of 1-points in a unit cubic region D(x) and a pure discrete 2-complex T m (x) for the combinatorial boundary ¶O m of the set V of all 1-points with respect to each m = 6,18,26. The arrows are oriented to the exterior of ¶O m .

For each 1-point con guration of a unit cube D(x), the con gurations of points of CubeBr m ′ (V; x) for m ′ = 6,18,26, Sk(I m (x)), Sk(J m (x)) and Sk(C m (x)) \ Sk(O m (x)) \ Sk(T m (x)) for m = 6,18,26 are shown with A (m ′ ,m) (x) for the adjustment in the cases of (m ′ , m) = (6, 18), [START_REF] Kong | Concepts of digital topology[END_REF][START_REF] Couprie | Simplicity surfaces: a new de nition of surfaces in Z ![END_REF]. "'


L Najman
email: l.najman@esiee.fr

Using Mathematical Morphology for Document Skew Estimation

We propose a concise definition of the skew angle of document, based on mathematical morphology. This definition has the advantages to be applicable both for binary and grey-scale images. We then discuss various possible implementations of this definition, and show that results we obtain are comparable to those of other existing algorithms.

INTRODUCTION 1.Foreword

We are going to apply the mathematical morphology formalism to document skew estimation. The use of mathematical morphology in the context of document image analysis allows us to give a clean definition of the skew angle of a document, through the definition of a new operator based on the closing by a segment and on the rose of directions. [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF] Our definition, relying on a well-established formalism, is straightforwardly applicable both to binary and grey-scale images. We can then propose several different possible efficient implementations of this operator in different contexts.

Skew estimation and document

The automatic processing of document images is now widespread and is performed for a variety of reasons including, for example optical character recognition. Often there are problems in the automatic processing because the document image is skewed, and it is advisable to detect or estimate the skew angle, and correct the skew, before applying the further image processing.

Various techniques have been proposed for automatic skew detection in document images. [START_REF] O'gorman | Executive briefing: document image analysis[END_REF] These are usually methods based on clustering of nearest neighbors, methods based on Hough transforms, or method involving projection profiles. This latest method is of interest for the comprehension of the one presented in this paper, and is thus described in the next section. Let us note that it is not straightforward to extend existing methods to grey-scale images. Numerous reviews can be found in the abundant existing litterature [START_REF] Amin | Comparative study of skew detection algorithms[END_REF][START_REF] Okun | Document skew estimation without angle restriction[END_REF] on skew estimation.

A WELL-KNOWN SKEW ESTIMATION ALGORITHM

A popular method for skew detection employs projection profiles. A projection profile is a histogram of the number of ON-pixel accumulated along parallel sample lines taken through the document. The profile may be at any angle, but often is taken horizontally along rows or vertically along columns, and these are called the horizontal and vertical projection profiles, respectively. For a document whose text lines span horizontally, the horizontal projection profile will have peaks whose widths are equal to the character height and mimimu height valleys whose widths are equal to the between-line spacing. For multicolumn documents, the vertical projection profile will have a plateau for each column, separated by valleys for the between-column and margin spacing.

The most straightforward use of the projection profile for skew detection is to compute it at a number of angles close to the expected orientation. [START_REF] Postl | Detection of linear oblique structures and skew scan in digitized documents[END_REF] For each angle, a measure is made of the variation in the bin heights along the profile, and the one with the maximum variation gives the skew angle. At the correct skew angle, since scan lines are aligned to text lines, the projection profile has maximum height peaks for text and valleys for between-line spacing. Modifications and improvements can be made to this general technique to more quickly iterate to correct the skew angle and more accurately determine it. [START_REF] Baird | The skew angle of printed documents[END_REF][START_REF] Bloomberg | Measuring document image skew and orientation[END_REF] All implementations of this kind of techniques rely on the maximization of a non-linear cost function. For instance, the cost function can be the sum of the squares of the successive differences of the projection profile bin heights in a given direction. No linear cost function can be chosen: if it were the case, we would compute a result proportional to the total number of all the black pixels, whatever the direction.

MATHEMATICAL MORPHOLOGY AND DOCUMENT

Founded in the sixties by G. Matheron and J. Serra, mathematical morphology is a theory for the analysis of spatial structures. It is called morphology because it aims at analysing the shape of unknown objects through the use of known shapes called structuring elements. A modern reference on this subject is the book of P. Soille. [START_REF] Soille | Morphological Image Analysis[END_REF] A structuring element is a small set used to probe the image under study. The shape and the size of the structuring element must be adapted to the geometric properties of the image objects to be processed. For instance, linear structuring elements are suited to the extraction of linear objects.

Basic morphological operators

The two basic morphological operators are the dilation and the erosion.

The dilation answers the question: does the structuring element hit the set? For a binary image X and a structuring element B, it is given by

X ⊕ B = {x|B x ∩ X = ∅} (1)
where B x = {b + x|b ∈ B} is the translated of B at point x, and for a grey-scale image f , by

(f ⊕ B)(x) = max b∈B f (x + b) (2)
The erosion is the dual of the dilation, and answers the question: does the structuring element fit the set? For a binary image X and a structuring element B, it is given by

X B = {x|B x ⊆ X} (3)
and for a grey-scale image f , by

(f B)(x) = min b∈B f (x + b) (4)
Advanced operators can be derived by combining the dilation and the erosion. For instance, a closing is a dilation followed by an erosion. Under the name of RLSA, and for a linear structuring element, the closing enjoys a wide popularity in the document analysis community.

Run-Length Smoothing closing algorithm

Let L 0 be a horizontal segment of unit length. The Run-Length Smoothing closing algorithm RLSA 9 is a closing of the binary image X by γL 0 :

RLSA(X) = X ⊕ γL 0 γL 0 ( 5
)
where γ is a size parameter.

The RLSA is a popular method for smoothing characters into smal, unrecognizable blobs. Assuming that the skew is corrected, the method merges characters into words, and words into text lines.

Using the morphological description of the RLSA, we can compute RLSA on lines oriented at an angle α different from 0, or we can apply it to a grey-scale image f . For an image A (binary or grey-scale), we thus define the operator

RLSA α (A) = A ⊕ γL α γL α (6)
Note that computing RLSA α is an idea quite similar to computing a projection profile at an angle α: we are looking for runs of pixels oriented in the α direction.

Rose of directions

Another useful operator is a kind of granulometry, in the sense that we are filtering our image according to a given family of structuring elements, and that we are computing a measure (namely the surface) on the resulting filtered image. In our case, we are not going to filter according to the size of the structuring element, but according to its direction.

Let L α be a segment of unit length and of angle α with respect to the horizontal direction, and let β ∈ a size parameter. We define the following operator ρ(α):

ρ(α) = Area(A βL α ) (7)
Ideally, the function ρ should be maximum at the predominant angle for the image. For a document image, this function should be maximum for the skew angle.

In the real world, however, it is difficult to extract from the function ρ applied directly on the original image useful information: except for very small β, the erosion will remove nearly all characters present on the document.

An operation related to equation 7 is the erosion not with L α but with a pair of points P α where

P α = {(0, 0)(cos α, sin α)} . (8)
Area of erosions by pair of points separated by a fixed distance but with varying orientations are sometimes represented in a polar diagram which is called a rose of directions. 1 This is the curve (ρ(α), α) α∈[0,360] .

Let us note that the morphological covariance K of an image F consists in measuring the volume (or the surface) of the image eroded by a pair of points P α .

K(f ; P α ) = V ol(f P α (f )) (9)
For binary images, it reduced to

K(A; P α ) = Area(A ∩ A (cos(α),sin(α)) ) (10)
which is the same as the rose of directions.

See [START_REF] Tuzikov | Extraction of grid lines on stamped metal pieces using mathematical morphology[END_REF][START_REF] Peyard | A dedicated hardware system for the extraction of grid patterns on stamped metal sheets[END_REF] for a recent industrial application on metal sheet. For binary images, the morphological covariance is the autocorrelation function. The correlation has been used for skew estimation, for instance using autocorrelation of projections. Correlation methods [START_REF] Yan | Skew correction of document images using interline cross-correlation[END_REF][START_REF] Chaudhuri | Robust detection of skew in document images[END_REF] require large text regions on a page and takes more time to compute than methods like profile projection.

Morphological skew estimator

Equation 7 cannot be applied directly on the original document image, and we have to filter the image before applying the operator. According to what has been presented in section 3.2, the logical idea is to filter the image by RLSA α . Thus, we are going to compute

ν(α) = Area(RLSA α (A) βL α ) (11)
and the skew angle of image A is provided by

skew(A) = arg max α ν(α) (12) = arg max α Area(A ⊕ γL α (γ + β)L α ) (13)
Let us note that for a precision of one degree, due to the grid structure of the image, the length γ + β of the line structuring element should be greater than 1 tan 1 = 58 pixels. Typically, the lenght γ + β can be of the order of the length of a line of the document. The value of γ should be of the order of the size of a letter of the document or the size of the space between words, whichever is the greater.

The comparison of this formula with the profile projection method is enlightening. Indeed, the projection profile is replaced by a closing and the non-linear cost function necessary for the projection profile method is replaced by the erosion morphological operator, which is non-linear. So, what have we gained? We now have an elegant definition for the skew angle of document, relying on a well-founded theory; moreover, operators of this theory have been intensively studied from the point of view of algorithmics. Indeed a straightforward implementation of equation 11 is prohibitive. In the sequel of the paper we are going to discuss various possible efficient implementations of this formula. Another benefit is that this formula is also directly applicable on grey-scale images (replacing the area by the sum), which is not the case for most skew estimation algorithms.

We illustrate equation 11 on the image 1. Figure 1 shows the erosion of RLSA α computed for various angles. We can see that the closer we are to the correct angle, the darker the image is.

Morphological algorithm for skew estimation

The naive method is to evaluate ν(α) at many α in a global search. This is expensive to compute for typical document images, containing more than 10 6 pixels. Thus, a traditional optimization algorithm that requires fewer function evaluations is applied. Given a initial bracketing interval, Brent's method 14 returns the position of the maximum to within a specified tolerance. It is advantageous to subsample the image prior to the optimization search. Downsampling and other implementation issues are discussed in section 5.

Brent's method is a kind of parabolic interpolation in which the values of six parameters a, b, u, v, w and x, are monitored. The parameters a and b are the limits of a bounding interval in which the minimum is located, x is the point with the lowest function value found so far, w is the point with the second lowest function value found so far, u is the point at which the function was evaluated most recently, and v is the previous value of w. The method is iterative.

According to Brent's method, parabolic interpolation is attempted fitting through the points x, v and w. In order to be considered acceptable, the current parabolic-interpolation step must (i) produce a new minimum which falls within the bounding interval (a, b), and (ii) imply a movement (amount of change) from the best current value, x, that is less than half the movement of the step before last. This second criterion ensures that the successive steps of the method will lead to convergence. In the worst case, where successive steps approximately alternate between parabolic steps and golden sections, there will ultimately be convergence thanks to the golden sections.

IMPLEMENTATIONS

The raw algorithm for computing erosion or dilation consists in calculating the minimum or maximum of pixels in the structuring element for each pixel of the image. There are (n -1) min/max comparisons per image pixel for a structuring element of n pixels. This number can be drastically reduced by appropriate algorithms and data structures. This section describes the best algorithms adapted to our problem.

Binary implementations 4.1.1. Constant time implementation

There exists several possibilities for binary implementation of dilation and erosion. For instance, one can use the method described by Nadadur and Haralick, [START_REF] Nadadur | Recursive binary dilation and erosion using digital line structuring elements in aritrary orientations[END_REF] which uses bresenham lines and requires only one maximum operation for erosion and one minimum operation for dilation transform, and one thresholding step and one translation step per result pixel. Such algorithms run at constant time with all possible length line structuring elements. However, we can find implementations requiring less than one operation per result pixel. This is the object of the next section.

Logarithmic decomposition

A logarithmic decomposition of convex sets has been proposed by Boomgaard and van Balen, [START_REF] Van Den Boomgaard | Methods for fast morphological image transform using bitmapped binary images[END_REF] based on a definition of extreme sets of a convex set given by Pecht. 17 Moreover, since w pixels of a binary image can be stored in a w-bit data type word, a logical operator can be simultaneously applied to w pixels using a bitwise operator. For the completeness of this paper, we are going to provide the main idea.

It is a well known fact in mathematical morphology that convex sets are divisible with respect to dilation. Let n be an integer, and S a convex set. We define nS = S ⊕ . . . ⊕ S n terms nS has the same shape as S but n times larger. Using the associativity of the dilation, X ⊕ nS can be computed by s succesive dilation by S. This is called linear decomposition. We can go further if we suppose that S is a two-dimensional convex symmetric set. In this case the linear decomposition can be replaced [START_REF] Pecht | Speeding up successive minkowski operations[END_REF] with

nS = S⊕ E(S) ⊕ . . . ⊕ E(S) (n -1) terms
where E(S) is the extreme set of S. For discrete sets, the extreme set E(S) contains just the vertices of the continuous convex hull of S. The linear decomposition can furthermore be replaced by the more efficient logarithmic decomposition. Using the extreme set, we have 2S = S ⊕ E(S), and by associativity 4S = 2S ⊕ E(2S) = S ⊕ E(S) ⊕ E(2S), and so on. In the general case for arbitrary n, we have

nS = S ⊕ E(S) ⊕ E(2S) ⊕ . . . ⊕ E(2 m-1 S) ⊕ E((n -2 m-1 )S) ( 14
)
where

m = [log 2 (n)].
If we apply that formula to our case, as an example we have

8L 0 = {• • •} ⊕ {• • •} ⊕ {• • • • •} ⊕ {• • • • • • • • •} (15)
which has a total of 9 pixels instead of the 19 pixels of 8L 0 .

For an efficient algorithmic implementation, we are using a bit-mapped representation of binary images. Each addressable word of memory/CPU contains 32 (or 64) horizontally consecutive pixels. Furthermore, morphological operations such as dilation or closing can be computed in parallel, leading to fast algorithms. Indeed, a dilation is the logical OR operation, and an erosion is the logical AND operation. The precise description of this implementation can be found in the paper by Boomgaard and van Balen. [START_REF] Van Den Boomgaard | Methods for fast morphological image transform using bitmapped binary images[END_REF] Thus, knowing that the surface can be efficiently computed by a hash table in the last step, one step of the skew estimation is computed in O((log(γ) + log(γ) log(β)) nm w ) where γ and β are the parameters of the rose of directions, nm is the number of pixels of the image, and w is the size of a machine word.

Fourier implementation

On dedicated hardware, it is possible to implement morphological dilation and erosion using a Fourier transform. [START_REF] Mazille | Mathematical morphology and convolution[END_REF][START_REF] Killinger | Morphological filtering using a fourier transform hologram[END_REF] Let us give the basic idea. The indicator function k A of a set A is given by k

A (x) = 1 if x ∈ A and k A (x) = 0 otherwise. The convolution of a function f by a function c is defined by c f (x) = u c(u)f (x -u).
The relationship between mathematical morphology and convolution lies in the following equation

k A k Ǩ (x) = Area(A ∩ K x ) ( 16
)
where Ǩ is the reflected set of K with respect to the origin. We have

A ⊕ K = {x|k A k Ǩ (x) > 0} (17)
A K = {x|k A k Ǩ (x) = Area(K)} (18)
As computing a convolution is easier using Fourier transform, using special processors which are currently commercially available based on the Fast Fourier Transform, the skew estimation algorithm described in this paper can be straightforwardly implemented on such processors.

Grey-scale implementation

Our formula for skew estimation can directly be used on grey-scale images, prior to binarization. We thus have to look for implementations of min/max filter for grey-scale images, and several of them can been found that run in constant time. For linear structuring elements, van Herk 20 proposed a recursive algorithm requiring 3 min/max comparisons per image pixel whatever the length of the line segment. This algorithm has been extended [START_REF] Soille | Recursive implementation of erosions and dilations along discrete lines at arbitrary angles[END_REF] to bresenham lines at arbitrary angles. However, this algorithm does not give the exact result, and Soille and Talbot [START_REF] Soille | Directional morphological filtering[END_REF] showed that an exact translation invariant implementation of min/max filters along a line segment of slope in the form of an irreducible fraction dy/dx can be achieved at the cost of 2 + k min/max comparisons per image pixel, where k = max(|dx|, |dy|).

TESTING THE ALGORITHM AND COMPARISON WITH OTHERS SKEW DETECTION METHODS

The tests we made (using the logarithmic decomposition implementation) validate the skew estimation, for a global computation time less than 0, 75 seconds on a Pentium III 733 Mhz for a 1214x1151 pixels image, including call to the wrapping code necessary to our development framework. This timing is comparable to timing of commercial packages or existing algorithms, albeit slower than some. Let us note that the test implementation is quite generic and portable, and thus does not take advantage of any processor-dedicated optimization.

We have to estimate the quality of the result of the morphological estimator, and to compare its accuracy with respect to other estimators. As document database, we have chosen the UW English Document Image Database -I. [START_REF] Chen | Uw english document image database -(i)[END_REF] This database contains 979 scanned images. For each of the scanned images, an estimate of the actual image skew is provided in the database, with an expected error of from 50 to 100 milli-degrees.

Choice of estimator for comparison

To our knowledge, there exist two methods using morphological tools as an important brick of the estimation.

• Chen and Haralick 24 used a more involved method that started with threshold reduction, [START_REF] Bloomberg | Image analysis using threshold reduction[END_REF] applied recursive morphological closings and openings to close up textlines and remove ascenders and descenders, determined connected components, fit the best line to the points in each set of connected components, and estimated a global skew by discarding outlier lines and averaging the remaining ones. The skew error reported was greater than 0.3 degrees on about 10 percent of the images in the UW English Document Image Database. [START_REF] Chen | Uw english document image database -(i)[END_REF] • Das and Chanda [START_REF] Das | A fast algorithm for skew detection of document images using morphology[END_REF] propose to use RLSA 0 , i.e. a closing of the image using line structuring element. That forms solid black bands corresponding to text lines. They open the closed image with small square structuring element to remove bumps. Then they register all 0 to 1 transitions, which provide lines. The skew is computed from the biggest lines. The method is reported to work for small skew angles.

One can note that contrary to our method, none of these methods is fully morphological We will not compare with those methods. D. Bloomberg, both in a paper [START_REF] Bloomberg | Measuring document image skew and orientation[END_REF] and on his web site, [START_REF] Bloomberg | Leptonica web pages[END_REF] claims that notwithstanding the size of skew estimation bibliography, there is only one method that can be recommended for both speed and accuracy. It is based on the maximization of the differential skew function, and is a variant of the projection profile method.

The sum for the i th sloped line at angle α is denoted by s(i, α), with s(i, α) = pixels in(i,α)scan,j∈w

p(i, j) (19)
S(α) = scanlines

s 2 (i, α) (20)
The differential skew function S δ (α) is defined as the sum of squares of differential counts on adjacent lines:

S δ (α) = scanlines (s(i, α) -s(i -1, α)) 2 (21)
As noted by Bloomberg, 7 in contrast to the sum skew function S(α), the width of the differential function can be much less than the textline asperity ratio. D. Bloomberg provides an excellent open source C implementation of this method. [START_REF] Bloomberg | Leptonica web pages[END_REF] The differential skew estimator is reported to work well on the UW English Document Image Database -I and we are going to compare this method with the morphological skew estimator.

Experiment and measurement on scanned images

It is well known that the assessment of skew measurement performance on scanned images is problematic. The first problem is that the actual image skew is not known. The approach taken with the UW database was to measure the skew of each image manually. Unfortunately, as noted in, 7 observation of many of the scanned images indicates that there is no single global skew. The most common scan artifact is that the skew at the top and bottom of the image differ, often by two hundred milli-degrees or more. There are also a few "outliers" where the measured angle is far from correct. Significant measured bias is occasionally induced by text fragments with non-conforming skew imaged from adjacent book pages.

For each of the 979 images of the UW database, we compute both the result of the morphological skew estimator, and the result of the differential skew estimator, and we compare those results with the "groundtruth" from the UW database. For the morphological skew estimator, we obtain the best results with parameters 64 for the dilation and 512 for the erosion, which corresponds respectively to the space between word and to the size of a line of text. The morphological skew estimator has a mean absolute error of 0.2 and a root-mean-square of 0.25 with respect to the ground-truth, while the differential skew estimator has a mean absolute error of 0.15 and a root-mean-square of 0.25. By downsizing the images by a factor 2, we obtain for the morphological skew estimator a mean absolute error of 0.25 and a root-mean-square of 0.31. Figure 3 gives the graph of the histogram and the cumulative histogram of absolute errors with respect to the estimated ground-truth on the UW database, both for the differential skew estimator and for the morphological skew estimator. One can see that the differential skew estimator behaves slightly better than the morphological skew estimator. We believe that part of the difference between the two methods lies in the fact that the two implementations have different optimization methods.

CONCLUSION

This paper demonstrates how to use mathematical morphology to accurately and rapidly estimate the skew angle of a document image. Testing shows that this method is competitive with the ones existing on the market.

The benefits of this paper are the following:

• We now have a clean mathematical definition of the skew angle of document as a function of morphological operators, that works both for binary and grey-scale images.

• We can find algorithmic implementations that are competitive, both in term of speed and precision.

Further testing is needed to assess if the accuracy of the skew estimation is enhanced if we apply our formula directly on grey-scale images, without any binarization of the document image.

Part of this work has been done while the author was with Océ Print Logic Technologies, 1, rue Jean Lemoine 94105 Créteil France. Patent Pending.

3 Figure 1 .

 31 Figure 1. Run-Length Smoothing Closing eroded at several angles for a document with skew of -3 degreeFigure 2 displays an example of the curve we maximize to obtain the correct skew. The maximum of this curve is at the correct skew angle.

Figure 2 .

 2 Figure 2. Curve of ν(α) for α ∈ [-7, 3] for a document with skew of -3 degree. The maximum of the curve is at -3 degree.

Figure 3 .

 3 Figure 3. Histogram (a) and cumulative histogram (b) of absolute errors with respect to the estimated "ground-truth" on the UW database both for the differential skew estimator and for the morphological estimator. 20 bins were used for the histogram and 50 bins for the cumulative histogram.