Nicolas Normand
email: nicolas.normand@polytech.univ-nantes.fr

Pierre Evenou
email: pierre.evenou@polytech.univ-nantes.fr

Medial Axis Lookup Table and Test Neighborhood Computation for 3D Chamfer Norms

Chamfer distances are discrete distances based on the propagation of local distances, or weights, defined in a mask. The medial axis, i.e. the centers of maximal balls (balls which are not contained in any other ball), is a powerful tool for shape representation and analysis. The extraction of maximal disks is performed in the general case by testing the inclusion of a ball in a local neighborhood with covering relations usually represented by lookup tables.

The proposed method determines if a mask induces a norm and in this case, computes the lookup tables and the test neighborhood based on geometric properties of the balls of chamfer norms, represented as H-polytopes. The method does not need to repeatedly scan the image space, and improves the computation time of both the test neighborhood detection and the lookup table computation.

Introduction

The distance transform DT X of a binary image X is a function that maps each point p to its distance to the closest background pixel i.e. with the radius of the largest open disk centered in p included in the image. Such a disk is said to be maximal if it is not contained in another disk also included in X. The set of centers of maximal disks, also called the medial axis, is a convenient description of binary images for many applications ranging from image coding to shape recognition. Its attractive properties are reversibility and (relative) compactness.

Algorithms for computing the distance transform are known for various discrete distances [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF][START_REF] Montanari | A method for obtaining skeletons using a quasi-euclidean distance[END_REF][START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF][START_REF] Coeurjolly | Optimal separable algorithms to compute the reverse Euclidean distance transformation and discrete medial axis in arbitrary dimension[END_REF][START_REF] Strand | Weighted distances based on neighborhood sequences[END_REF]. In this paper, we will focus on chamfer (or weighted) distances which are defined by a set of weighted vectors described by a mask, called the chamfer mask. The classical medial axis extraction method is based on the removal of non maximal disks in the distance transform. It is thus mandatory to describe the covering relation of disks, or at least the transitive reduction of this relation. For simple distances, this knowledge is summarized in a local maximum criterion [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF] but a most general method uses lookup tables for that purpose [START_REF] Borgefors | Centres of maximal discs in the 5-7-11 distance transforms[END_REF].

In this paper, we propose a two-phase method that determines if a mask induces a norm and in this case, computes the lookup tables and the test neighborhood based on geometric properties of the balls of chamfer norms.

The first phase starts from the 3D chamfer mask and produces a triangulation of the chamfer neighbors inversely pondered by the chamfer weights. The algorithm conceptually operates on the convex hull of these weighted vectors and creates a triangulation of this convex hull. It produces two results: a norm condition check and a description of the geometry of the chamfer balls for the given chamfer mask. The second phase is performed only for norms; from the general geometry obtained during the previous phase, each chamfer ball is described as the intersection of a set of half-spaces (H-description). The norm condition is sufficient for the balls to be convex, ensuring the validity of their H-description. This H-description is then used to compute the test neighborhood T and the LUT values.

Basic notions, definitions and known results about chamfer disks and medial axis lookup tables are recalled in section 2. Then section 3 justifies the use of polytope formalism in our context and presents the principles of the method. In section 4, a triangulation algorithm is given for the 3D case. The computation of the test neighborhood T and of the LUT values is explained in section 5. Finally, section 6 gives speed-up figures for the overall algorithm compared to a reference implementation [START_REF] Rémy | Medial axis for chamfer distances: computing look-up tables and neighbourhoods in 2D or 3D[END_REF]. Notice that we do not address the actual computation of the medial axis which remains unchanged.

This paper is a more detailed version of [8] with an added method for the triangulation of 3D balls and the detection of norm conditions. The algorithms for LUT and T computation were adapted to the 3D case from [8].

Chamfer Distances and Medial Axis

Discrete Distances

Definition 1 (Discrete distance, metric and norm). Consider a function d : Z n × Z n → N and the following properties ∀x, y, z ∈ Z n , ∀λ ∈ Z:

1. positive definiteness d(x, y) ≥ 0 and d(x, y) = 0 ⇔ x = y , 2. symmetry d(x, y) = d(y, x) , 3. triangle inequality d(x, z) ≤ d(x, y) + d(y, z) , 4. translation invariance d(x + z, y + z) = d(x, y) , 5. positive homogeneity d(λx, λy) = |λ| • d(x, y) .
d is called a distance if it verifies conditions 1 and 2, a metric with conditions 1 to 3 and induces a norm if it also satisfies conditions 4 and 5.

2

P R E P R I N T Most discrete distances are built from a definition of neighborhood and connected paths (path-based distances), the distance between two points being equal to the length of the shortest path that joins them [START_REF] Rosenfeld | Distances functions on digital pictures[END_REF]. Distance functions differ by the neighborhoods used to build paths and by the way path lengths are measured. For the simple distance d 4 (denoted by d in [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF]), defined in the square grid Z 2 , each pixel has four neighbors located at its top, left, bottom and right edges. Similarly, for distance d 8 (d * in [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF]), each pixel has four extra diagonally located neighbors. In both cases, d 4 and d 8 , the length of a path is defined as its number of displacements, whereas it is measured as a weighted sum of displacements for chamfer distances [START_REF] Montanari | A method for obtaining skeletons using a quasi-euclidean distance[END_REF][START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF] or by the displacements allowed at each step for neighborhood sequence distances [START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF][START_REF] Rosenfeld | Distances functions on digital pictures[END_REF], or even by a mixed approach of weighted neighborhood sequence paths [START_REF] Strand | Weighted distances based on neighborhood sequences[END_REF].

For a given distance d, the closed ball B ≤ and open ball B < of center c and radius r are the sets of points of Z n :

B < (c, r) = {p : d(c, p) < r} ; B ≤ (c, r) = {p : d(c, p) ≤ r} . (1) Since the codomain of d is N, we have ∀r ∈ N, d(c, p) ≤ r ⇔ d(c, p) < r + 1. So: ∀r ∈ N, B ≤ (c, r) = B < (c, r + 1) . ( 2
)
Definition 2 (Distance transform). The distance transform DT X of the binary image X is a function that maps each point p to its distance to the closest background pixel:

DT X : Z n → N DT X (p) = min d(p, q) : q ∈ Z n \ X . ( 3
)
Alternatively, since all points at a distance less than DT X (p) from p belongs to X (B < (p, DT X (p)) ⊂ X) and at least one point at a distance equal to DT X (p) is not in X (B ≤ (p, DT X (p)) = B < (p, DT X + 1) ⊂ X) then DT X (p) is the radius of the largest open disk centered in p included in the image:

DT X (p) = max r : B < (p, r) ⊂ X . ( 4
)
Definition 3 (Medial axis). The medial axis MA X of the binary image X is the set of centers of maximal open balls in X valued with their radii:

MA X : Z n → N MA X (p) = 0 if ∃ q, r s.t. B < (p, DT X (p)) B < (q, r ) ⊂ X , DT X (p) otherwise .
(5)

Chamfer Distances

Definition 4 (Chamfer mask [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF]). A weighting M = ( -→ v ; w) is a vector -→ v of Z n associated with a weight w (or local distance). A chamfer mask M is a central-symmetric set of weightings having positive weights and non-null displacements, and containing at least one basis of

Z n : M = {M i ∈ Z n × N * } 1≤i≤m P R E P R I N T
The grid Z n is symmetric with respect to the hyperplanes normal to the axes and to the bisectors (G-symmetry). This divides Z n in 2 n .n! subspaces (there are 8 of them for Z 2 and 48 for Z 3 ). The particular subspace x n ≥ . . . ≥ x 1 ≥ 0 is called the generator cone or simply generator and is denoted G. From every point in Z n , we can determine its unique G-symmetrical point in G by ordering the absolute values of its components in decreasing order. Conversely, from any point p in G, we can derive all its G-symmetrical points by generating all n 1 ! permutations of the components of p and all 2 n2 combinations of their signs where n 1 and n 2 are respectively the number of different absolute values of the components of p and the number of non null components.

Chamfer masks are usually restricted to G for simplicity. Weightings are then only given in G and we denote M| G = {M i ∈ G × N * } 1≤i≤m , the chamfer mask restricted to G. A usual ordering of the points in G is the lexicographical order; p is before q if for some i: p 0..i-1 = q 0..i-1 and p i < q i . According to this ordering, a common naming scheme assigns alphabetic letters to visible points (such that gcd(p) = 1). For instance, in the 2D grid:

-→ a = (1 0), -→ b = (1 1), -→ c = (2 1), -→ d = (3 1), -→ e = ( 3 2
)... The corresponding chamfer weights are respectively named a, b, c, d, e... Definition 5 (Chamfer distance [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF]). Consider the chamfer mask

M = {( -→ v i ; w i ) ∈ Z n × N * } 1≤i≤m .
The chamfer (or weighted) distance between two points p and q is:

d(p, q) = min λ i w i : p + λ i -→ v i = q, λ i ∈ N, 1 ≤ i ≤ m . ( 6
)
Paths between two points p and q can be produced by chaining displacements. The length of a path is the sum of the weights associated with the displacements and the distance between p and q is the length of the shortest path. Any chamfer masks defines a metric [START_REF] Verwer | Local distances for distance transformations in two and three dimensions[END_REF]. However a chamfer mask only generates a norm when proper conditions on the mask neighbors and on the corresponding weights permits a triangulation of the ball in influence cones [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF][START_REF] Rémy | Normes de chanfrein et axe médian dans le volume discret[END_REF]. When a mask induces a norm then all its balls are convex and therefore can be represented as polytopes.

Chamfer Medial Axis

For simple distances d 4 and d 8 , the medial axis extraction can be performed by the detection of local maxima in the distance map [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF]. Chamfer distances raise a first complication even for small masks as soon as the weights are not unitary. Since all possible values of distance are not achievable, two different radii r and r may correspond to the same set of discrete points. The radii r and r are said to be equivalent. Since the distance transform labels pixels with the greatest equivalent radius, criteria based on radius difference fail to recognize equivalent disks as being covered by other disks. In the case of 3 × 3 2D masks or 3 × 3 × 3 3D masks, a simple relabeling of distance map values with the smallest equivalent radius is sufficient [START_REF] Arcelli | Finding local maxima in a pseudo-Euclidian distance transform[END_REF][START_REF] Svensson | Digital distance transforms in 3D images using information from neighbourhoods up to 5×5×5[END_REF]. However this method fails for greater masks and the most general method for medial axis extraction from

P R E P R I N T T = -→ a = (1 0); -→ b = (1 1); -→ c = ( 2
B < (O, 14) ⊆ B < (O + -→ c , 23) as shown in (c) but B < (O, 14) ⊆ B < (O + -→ c , 22) because of the point (-2, 1) (d). In terms of closed balls, B ≤ (O, 13) ⊆ B ≤ (O + -→ c , 22) but B ≤ (O, 13) ⊆ B ≤ (O + -→ c , 21).
the distance map involves lookup tables (LUT) that represent for each neighbor -→ v i in a set called the test neighborhood T and for each radius r 1 , the minimal open ball covering

B < (O, r 1 ,) in direction -→ v i [6]: Lut-→ v i (r 1 ) = min r 2 : B < (O + -→ v i , r 1 ) ⊆ B < (O, r 2 ) .
Equivalently, using closed balls (considering (2)):

Lut-→ v i (r 1 ) = 1 + min r 2 : B ≤ (O + -→ v i , r 1 -1) ⊆ B ≤ (O, r 2 ) . ( 7
)
As noticed by Thiel, the test neighborhood T is not necessarily equal to the chamfer mask M [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF].

Medial Axis LUT Coefficients

A general method for LUT coefficient computation was given by Rémy and Thiel [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF][START_REF] Thiel | Les distances de chanfrein en analyse d'images : fondements et applications[END_REF][START_REF] Rémy | Computing 3D medial axis for chamfer distances[END_REF]. The idea is that the disk covering relation can be extracted directly from values of distance to the origin. If d(O, p) = r 1 and d(O, p + -→ v i ) = r 2 , we can deduce the following:

p ∈ B ≤ (O, r 1 ) = B < (O, r 1 + 1) , p + -→ v i / ∈ B ≤ (O, r 2 -1) = B < (O, r 2 ) , P R E P R I N T hence B < (O + -→ v i , r 1 + 1) ⊂ B < (O, r 2 ) and Lut-→ v i (r 1 + 1) > r 2 . If ∀p, d(O, p) ≤ r 1 ⇒ d(O, p + -→ v i ) ≤ r 2 then Lut-→ v i (r 1 + 1) = r 2 + 1. Finally, Lut-→ v i (r) = 1 + max d(O, p + -→ v i ) : d(O, p) < r .
This method only requires one scan of the distance function for each displacement -→ v i . Moreover, the visited area may be restricted according to the symmetries of the chamfer mask. The order of complexity is about O(mL n ) for m neighbors if we limit the computation of the distance function to an image of size L n .

Medial Axis Test Neighborhood

Thiel observed that the chamfer mask is not adequate to compute the medial axis [15, p. 81 (O, 351). This inclusion relation is neither detected with the vector -→ c = (2 1) nor with the other vectors of the chamfer mask. Remy and Thiel then introduced a LUT Mask (called test neighborhood and denoted by T (R) here) for that purpose [START_REF] Rémy | Normes de chanfrein et axe médian dans le volume discret[END_REF]. T (R) is the minimal set of vectors sufficient to detect the medial axis for shapes whose inner radius (the radius of a greatest ball) is less than or equal to R. In the previous d 14,20,31,44 example, the point (4, 2) is not in the chamfer mask but should be in T (R) for R greater than 350.

A test neighborhood incompleteness produces extra points in the medial axis (undetected ball coverings). A general method for both detecting and validating T is based on the computation of the medial axis of all disks [START_REF] Rémy | Medial axis for chamfer distances: computing look-up tables and neighbourhoods in 2D or 3D[END_REF]. When T is complete, the medial axis is restricted to the center of the disk, when extra points remains, they are added to T . This neighborhood determination was proven to work in any dimension n ≥ 2. However it is time consuming even when taking advantage of the mask symmetries. [START_REF] Ziegler | Lectures on Polytopes[END_REF] Definition 6 (Polyhedron). A convex polyhedron is the intersection of a finite set of half-hyperplanes.

H-Polytopes and Chamfer Balls

General H-Polytopes

Definition 7 (Polytope).

A polytope is the convex hull of a finite set of points.

Theorem 1 (Weyl-Minkowski). A subset of Euclidean space is a polytope if and only if it is a bounded convex polyhedron.

As a result, a polytope in R n can be represented either as the convex hull of its k vertices (V-representation) or by a set of l half-planes (H-representation):

P = conv({p i } 1≤i≤k ) = p = k i=1 α i p i : α i ∈ R + and k i=1 α i = 1 , (8) P R E P R I N T P = x : Ax ≤ y , ( 9
)
where A is a l × n matrix, y a vector of n values that we name H-coefficients of P . Given two vectors -→ u and -→ v , we denote

-→ u ≤ -→ v if and only if ∀i, -→ u i ≤ -→ v i .
Definition 8 (Simplicial cone). A simplicial cone C o,U from a point o is a cone of dimension m defined by a set U of m independant vectors.

In a simplicial cone, each point is representable by a unique (up to a permutation) non negative combination of the vectors of the cone, i.e.:

p ∈ C o,U ⇒ -→ op = m i=1 α i -→ u i , where (α i ) ∈ R m + is unique.
Definition 9 (Discrete polytope). A discrete polytope Q is the intersection of a polytope P in R n with Z n (Gauss dicretization of P).

Definition 10 (Unimodular cone). A unimodular cone C o,U is a simplicial cone defined by a set of integral vectors U that generate all the integral points of C o,U .

In a unimodular cone, all discrete points are representable as a unique non negative integral combination of the vectors of the cone.

A cone C o,U of dimension n is unimodular if and only if | det(U )| = 1.
Then U is the basis of a unimodular point lattice, equivalent to Z n , so each integral point can be reached.

Minimal Representation

Many operations on R n polytopes in either V or H representation often require a minimal representation. The redundancy removal is the elimination of unnecessary points in the V-representation or unnecessary inequalities in the H-representation of polytopes. Since our purpose is mainly to compare Hpolytopes defined with the same matrix A, no inequality removal is needed. However, for some operations, H-representations of discrete polytopes must be minimal in terms of H-coefficients.

Definition 11 (Minimal parameter representation).

A minimal parameter H-representation of a discrete polytope P , denoted H-representation, is a H-representation of P = x : Ax ≤ y such that y is minimal:

P = {x ∈ Z n : Ax ≤ y} and ∀i ∈ [1..l], ∃x ∈ P : A i x = y i , ( 10
)
where A i stands for the i th line of the matrix A.

The H function, introduced for convenience, gives the minimal parameter vector for a given polytope P : H(P ) = max Ax : x ∈ P . As a consequence, In the H case, the three equalities are verified for the same point [START_REF] Borgefors | Centres of maximal discs in the 5-7-11 distance transforms[END_REF][START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF]. Notice that although coefficient values are minimal, this representation is still redundant: the second inequality could be removed.

{x : Ax ≤ H(P )} is the H-representation of P = {x : Ax ≤ y}.

H-Polytope Translation

Let P = {x : Ax ≤ y} be a H-polytope. The translated of P by -→ v which is also the Minkowski sum of P and { -→ v } is:

(P )-→ v = P ⊕ { -→ v } = {x + -→ v : Ax ≤ y} = {x : Ax ≤ y + A -→ v } . ( 11
)
The translation of a minimal representation gives a minimal representation.

Covering Test

Let P = {x : Ax ≤ y} and Q = {x : Ax ≤ z} be two polyhedra represented by the same matrix A but having different sets of H-coefficients y and z. We have:

y ≤ z ⇒ P ⊆ Q . ( 12
)
Furthermore, if the H-description of the enclosed polyhedron has minimal coefficients, the condition is also necessary:

y = H(P ) ≤ z ⇔ P ⊆ Q . ( 13
)
3.2. Geometry of the Chamfer Ball [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF][START_REF] Rémy | Normes de chanfrein et axe médian dans le volume discret[END_REF] Describing balls of chamfer norms as polygons in 2D and polyhedra in higher dimensions is not new [START_REF] Verwer | Local distances for distance transformations in two and three dimensions[END_REF]. Thiel and others have extensively studied chamfer ball geometry from this point of view [START_REF] Rémy | Normes de chanfrein et axe médian dans le volume discret[END_REF][START_REF] Thiel | Les distances de chanfrein en analyse d'images : fondements et applications[END_REF][START_REF] Borgefors | Weighted digital distance transforms in four dimensions[END_REF]. Our purpose is to introduce properties specific to the H-representation of these convex balls.

We can deduce from (1) and ( 6) a recursive construction of chamfer balls:

B ≤ (O, r) = B ≤ (O, r -1) ∪ 0≤i≤m B ≤ (O + -→ v i , r -w i ) . ( 14
)
This construction is sufficient to obtain an exact description of the balls in the case of chamfer norms as we will see later (section 5.1).

P R E P R I N T Definition 12 (Rational ball). Consider a chamfer mask M, the rational unitary ball or simply rational ball B R is the convex hull of the rational points

-→ v i /w i , ( -→ v i , w i ) ∈ M: B R = conv -→ v i w i : ( -→ v i ; w i ) ∈ M = m i=1 α i -→ v i w i : α i ∈ R + , m i=1 α i = 1 . ( 15
)
Note that the rational ball, as it is introduced here, is convex by definition and is different, in this regard, from the equivalent rational ball described in [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF][START_REF] Rémy | Normes de chanfrein et axe médian dans le volume discret[END_REF][START_REF] Fouard | 3-D chamfer distances and norms in anisotropic grids[END_REF].

Proposition 1. The homothetic of B R , rB R contains B ≤ (O, r): B ≤ (O, r) ⊆ rB R .
Proof. Let p be a point of B ≤ (O, r) and

m i=1 λ i -→ v i , λ i ∈ N be a minimal path between O and p. Then d(O, p) ≤ r and d(O, p) = m i=1 λ i w i . We can describe p as a combination of O, r -→ v 1 /w 1 , . . . , r -→ v n /w n : p = r-d(O,p) r •O+ m i=1 λiwi r • r wi -→ v i , with r-d(O,p) r ≥ 0, λiwi r ≥ 0 ∀i and r-d(O,p) r + m i=1 λiwi r = 1.
Then, by definition of conv:

p ∈ conv(O, r w1 -→ v 1 , . . . , r wn -→ v n ) ⊂ rB R and every point of B ≤ (O, r) is in rB R .
Definition 13 (Normal vector [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF]). We call normal vector to a facet F of B R , the unique vector -→ F orthogonal to the facet such that ∀p ∈ F, -→ F • p = 1.

In the simplicial cone spanned by vectors

-→ v 1 w1 . . . -→ v n wn , -→ F • -→ v 1 w1 . . . -→ v n wn = (1 . . . 1
). This implies:

-→ F = (w 1 , . . . , w n ) • ( -→ v 1 | . . . | -→ v n ) -1 . ( 16
)
The

vector α = ( -→ v 1 | . . . | -→ v n ) -1 p give the unique representation of p as a combina- tion of -→ v 1 , . . . , -→ v n . If α has integer coefficients then there is a path α( -→ v 1 | . . . | -→ v n ) to p and its length is (w 1 , . . . , w n ) • α = -→ F • p.
The normal vector is equivalent to the elementary displacement and discrete gradient of the cone intercepted by F as defined for chamfer norms in [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF] and [START_REF] Rémy | Normes de chanfrein et axe médian dans le volume discret[END_REF]. Note that the normal vector is defined here for general chamfer masks, whether they induce a norm or not.

For instance, with the chamfer norm d 5,7,11 , the point (3, 1) is in the cone spanned by the vectors -→ a = (1 0) and -→ c = (2 1) and the weights involved are 5 and 11. The distance between the origin and the point (3, 1) is then [10, (4.32)]:

d(O, (3, 1)) = 5 11 • 1 2 0 1 -1 • 3 1 = 5 1 • 3 1 = 16 .
P R E P R I N T

Mask Conditions for Chamfer Norms

Proposition 2 (Distance lower bound). If F is a facet of the rational ball and -→ F its normal vector, then:

d(O, p) ≥ -→ F • p . ( 17
)
Proof. Let p = i λ i -→ v i be a minimal path to p, so d(O, p) = i λ i w i . The facet F is supported by the hyperplane {x ∈ R n : -→ F • x = 1}.
Due to convexity, the unitary ball is included in the halfspace {x ∈ R n :

-→ F • x ≤ 1}. This applies to the vertices of the unitary ball -→ v i /w i :

-

→ F • -→ v i /w i ≤ 1. By linearity of the dot product, -→ F • -→ v i ≤ w i and -→ F • ( i λ i -→ v i ) = i λ i -→ F • -→ v i ≤ i λ i w i . Hence, -→ F • i λ i -→ v i is always less than or equal to the length of the path i λ i -→ v i .
This lower bound applies for any point p, whether or not it is located in a cone that intercepts F; d does not need to induce a norm.

Proposition 3 (Homegeneity in a unimodular cone). Let C O,U be a unimodular cone that intercepts a unique facet F of B R . The distance d is homogenous in the cone and equal to -→ F • p, ∀p ∈ C O,U if and only if:

d(O, O + -→ u ) = -→ F • -→ u , ∀ -→ u ∈ U . ( 18
) Proof. According to proposition 2, d(O, p) ≥ -→ F • p, ∀p ∈ C O,U . Since C O,U is unimodular, p is reachable by a non negative integral combination of vectors of U , -→ Op = λ i -→ u i .

By triangle inequality and proposition hypotheses, d(O, p) ≤ λ i d(O, O + -→ u

i ) = λ i -→ F • -→ u i = -→ F • p. So -→ F • p is
→ u = α • ( -→ v 1 | . . . | -→ v n ). The components of α = ( -→ v 1 | . . . | -→ v n ) -1 • -→ u are not necessarily integers (they are if ∆ = det( -→ v 1 | . . . | -→ v n ) = 1). However, ∆( -→ v 1 | . . . | -→ v n ) -1 has

integral coefficients hence ∆ -→ u can be uniquely represented as an integral combination of the vectors

-→ v i : ∆ -→ u = ∆( -→ v 1 | . . . | -→ v n ) • α, so there is a path to O + ∆ -→ u whose length is: ∆(w 1 . . . w n ) • α = ∆(w 1 . . . w n ) • ( -→ v 1 | . . . | -→ v n ) -1 • -→ u = ∆ -→ F • p ≥ d(O, ∆ -→ u ). Thus d(O, O + ∆ -→ u ) < ∆d(O, O + -→ u
) and the distance is not homogenous in the cone.

According to this proposition, the homogeneity in a nD unimodular cone can be determined by only testing the vectors of the cone (which are not necessarily in the chamfer mask). We derive a necessary and sufficient condition for norms (a sufficient condition, based on a unimodular triangulation, is given in [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF]):

Corollary 1 (Norm condition). Let U i be a set of vectors sets U i with the following properties:

P R E P R I N T • each unimodular cone C i = C O,Ui intersects a unique facet of B R ;
• the union of the cones covers the entire space Z n .

The chamfer mask M induces a norm if and only if:

d(O, O + -→ u ) = -→ F • -→ u , ∀ -→ u ∈ i U i . ( 19
)
Proposition 4 (Direct distance formulation). If d induces a norm, the distance from O to any point p is:

d(O, p) = max 1≤i≤l -→ F i • p ( 20
)
where l is the number of facets of B R and -→ F i is the normal vector to the i th facet of B R . This formula does not require to determine in which cone lies p.

Proof. The inequality d(O, p) ≥ -→ F i • p holds for all facets due to proposition 2 and the equality d(O, p) =

-→ F i • p holds, according to proposition 3, for all the facets of B R (there is at least one of them) that intersect [O, p).

H-Representation of Chamfer Norm Balls

The H-representation of chamfer balls is directly derived from [START_REF] Chand | An algorithm for convex polytopes[END_REF]:

p ∈ B ≤ (O, r) ⇔ max 1≤i≤l { -→ F i • p} ≤ r ⇔ A M • p ≤ y (21)
where A M is a H-representation matrix depending only on the chamfer mask M. The number of rows in A M is equal to the number l of facets of B R , each line of the matrix A M is computed with [START_REF] Rémy | Computing 3D medial axis for chamfer distances[END_REF] and y is a column vector whose values are r. For instance, the H-representation matrix of d 5,7,11 balls restricted to G is A M = 5 1 4 3 where (5 1) and ( 43) are the normal vectors of the two facets of B R and

B ≤ (O, r) = p ∈ Z n : 5 1 4 3 • p ≤ r r .
Note that if a polytope is not centered in O, the simplification due to symmetries do not hold and the full set of H-coefficients is needed, unless we ensure that the H-coefficents for the hyperplanes in the generator G are greater than H-coefficients for the corresponding symmetric cones. This is the case when a G-symmetric polytope is translated by a vector in G. Proof. The smallest ball that covers the polytope P must contain its furthest point from the origin: min r ∈ N :

Proposition 5 (Furthest point). Let A M be the matrix defined by a chamfer mask M generating a norm. The furthest point from the origin in the H-polytope P = {p : A M • p ≤ H(P )} is at a distance equal to the greatest component of H(P ).

Proof. By construction of A M , (20) is equivalent to d(O, p) = max

P ⊆ B ≤ (O, r) = max p∈P d(O, p) = max 1≤i≤l H i (P ) .
Definition 14 (Covering function). We call covering function of a set X of points of Z n the function C X which assigns to each point p of Z n , the radius of the minimal ball centered in p covering X:

C X : Z n → N C X (p) = min r : X ⊆ B ≤ (p, r) .
The covering function of the chamfer ball B ≤ (O, r) at point p gives the radius of the minimal ball centered in p that contains B ≤ (O, r) or, by central symmetry of the chamfer balls, to the radius of the minimal ball centered in O covering B ≤ (p, r) and therefore it is the maximal component of the H-representation of B ≤ (p, r):

C B ≤ (O,r) (p) = max H(B ≤ (p, r) = max H(B ≤ (O, r)) + A M • p . ( 22
)
One can notice that the covering function of the zero radius disk is equal to the distance function, as is the distance transform of the complement of this disk:

C B ≤ (O,0) (p) = DT Z n \{O} (p) = d(p, O) = d(O, p) .

Definition 15 (Covering cone).

A covering cone C o,U in the covering function C X is a unimodular cone defined by a vertex o and a set of vectors U , such that each ball centered in p is included in the ball centered in p + -→ u , -→ u ∈ U :

∀p ∈ C o,U , ∀ -→ u ∈ U, B ≤ (p, C X (p)) B ≤ (p + -→ u , C X (p + -→ u )) . Proposition 7. If C o,U is a covering cone in C B ≤ (O,r
) then for any point q in C o,U \ U \ {o} there exists p such that the ball centered in q includes the ball centered in p:

B ≤ (O, r) B ≤ (p, C B ≤ (O,r) (p)) B ≤ (q, C B ≤ (O,r) (q)) . Proof. If U is not empty, there is -→ u ∈ U such that p = q --→ u is in C o,U \ {o}. Then, by definition of C o,U , B ≤ (p, C B ≤ (O,r) (p)) B ≤ (q, C B ≤ (O,r) (q)). In addition, B ≤ (O, r) ⊆ B ≤ (p, C B ≤ (O,

. l] and a point o such that

A Mj • -→ u i , ∀ -→ u i ∈ U and H j (B ≤ (o, r)) are maximal then C o,U is a covering cone in C B ≤ (O,r) .
Proof. Let j be the row number of a maximal component of H(B ≤ (o, r)) and A M • -→ u i , ∀ -→ u i ∈ U , then j is a maximal component of any positive linear combination of these vectors. Let p be any point in

C o,U , p = o + i λ i -→ u i , λ i ∈ N. From (22), we deduce that, in C o,U , C B ≤ (O,r)
takes the form of an affine function of the components λ i :

C B ≤ (O,r) (p) = max H (B ≤ (o, r)) + λ i A M -→ u i = H j (B ≤ (o, r)) + λ i A Mj -→ u i .
In the same way, the radius of the minimal ball centered in p

+ -→ u , -→ u ∈ U covering B ≤ (O, r) is: C B ≤ (O,r) (p + -→ u ) = H j (B ≤ (o, r)) + λ i A Mj -→ u i + A Mj -→ u = C B ≤ (p,C B ≤ (O,r) (p)) (O + -→ u ) .
In other words,

B ≤ (p, C B ≤ (O,r) (p)) is contained in B ≤ (p + -→ u , C B ≤ (O,r) (p + -→ u ))
and C o,U is a covering cone according to definition 15.

Triangulation of 3D Chamfer Balls

The first stage of the method consists in decomposing the generator cone G in unimodular cones that each intercepts a unique facet of the rational ball. This triangulation has three goals:

1. compute the distance gradients which are vectors normal to the facets of the rational ball. These distance gradients are the rows of the matrix A M ; 2. provide a set of unimodular cones suitable for the recursive partition of the generator G used in the next stage; 3. verify that the chamfer mask generates a norm by checking homogeneity in each cone.

Convex Hull Vertices in the Generator Cone

The expected result is a triangulation of the part of the rational ball B R that lies in the generator G. The input data consists in vectors from the chamfer mask M| G divided by their corresponding weights (to form the equivalent unitary ball). In the algorithm, these points in Q n are represented as weightings, and all computations are based on integer arithmetics. The weightings of the chamfer mask are not sufficient alone because:

G conv -→ v w , ( -→ v ; w) ∈ M = conv -→ v w , ( -→ v ; w) ∈ M| G .
P R E P R I N T

The intersection of B R with G introduces new vertices that are points in edges of B R that intersect the faces of G. This is illustrated in fig. 3a; vectors -→ e = (2 1 1) and -→ o = (4 1 1) are not in the chamfer mask and thus can not be vertices of B R , but they appear in the intersection of B R with (O, { -→ a , -→ c }). Since the edges of B R are not known yet, all the intersections of the faces of G with couples of rational vertices of B R must be added to the input of the algorithm. These points are rational combinations of the vertices of the intersected edges, but are also represented as weightings with integer coefficients. Among those, the majority of the points that lie in the interior of B R will be discarded during the computation of the convex hull.

Gift-Wrapping Convex Hull and Triangulation of the Rational Ball

Many different algorithms are available for convex hull computation. We chose to implement an modified version of the algorithm known as gift wrapping [START_REF] Chand | An algorithm for convex polytopes[END_REF] for two reasons: first it naturally produces a triangulation of the convex hull i.e. a partition in simplicial cones and second, a subpartition of these simplicial cones in unimodular cones and the face lattice of the cones can be easily integrated as by-products of the algorithm.

Input: A set of points S find a starting edge of the convex hull (p, q) 1 push (p, q) in , the list of active edges Algorithm 1: (Simplified) gift wrapping algorithm [START_REF] Chand | An algorithm for convex polytopes[END_REF] Starting with a first edge, the gift wrapping algorithm (1) iteratively extends the convex hull by appending a new triangle and its edges. Since the convex hull is computed in G, one vertex of the starting edge can be chosen as the point in the first axis of G that is closest to the origin: p = (x 1 , 0, 0) (with maximal x 1 if several points exist on this axis) and the second vertex is the point q that maximizes the angle (O, p, q) in a facet of G containing p. The algorithm proceeds until exhaustion of active edges. { e, h} The intersection of the rational ball for mask 7,8,11,14 and G is pictured in fig. 3b. A unimodular triangulation without any new point of the face (

P R E

P R I N T -→ c -→ e -→ o -→ b -→ d -→ a z x y (a) Rational ball -→ a -→ b -→ c -→ d -→ e -→ o -→ h z y (b) Unimodular triangulation
{ b, h, o} { b, e} G { c, e} { a, d} { b, d, o} { a, d, o} { b, d} { d, o} { b, c} { b, o} { b, h} { b, e, h} { a, o} { a} { c} { b} { e} { d} { h} { o} { b, c, e} ∅ { h, o}
-→ b = (1 1 0), -→ d = (2 1 0), -→ e = (2 1 1), -→ o = (4 1 1)) can not be made because det( -→ d , -→ e , -→ o ) = det( -→ b , -→ e , -→ o ) = 2
. This is solved by introducing the vector -→ h . The cone lattice corresponding to this unimodular triangulation is pictured in fig. 4.

Homogeneity check

The unimodular triangulation of B R produces a set of cones that cover the generator G and their corresponding face normal vectors. Using corollary 1, by comparing d(O, O + -→ u ), obtained by propagating the chamfer weights, with -→ F i • -→ u , where -→ F is the normal vector of one of the faces that intersects [O, -→ u ), we determine if the distance is homogenous in G, so by symmetry, in Z 3 .

LUT and Test Neighborhood Computation for 3D chamfer norms

H-Representation of Chamfer Balls

The computation of the LUT is based on a H-representation of the chamfer norm balls. All share the same matrix A M which depends only on the chamfer mask (21). The H-coefficients of the balls are computed iteratively from the ball of radius 0, B ≤ (O, 0) = x : A M x = 0 using ( 14) and [START_REF] Verwer | Local distances for distance transformations in two and three dimensions[END_REF]:

H(B ≤ (O, r)) = max H(B ≤ (O, r -1)); H(B ≤ (O, r -w i )) + A M -→ v i , i ∈ [1..m] .

LUT values

LUT values are directly obtained from the covering function (22):

Lut-→ v i [r] = 1 + C B ≤ (O,r-1) (O + -→ v i ) = 1 + max H(B ≤ (O, r -1)) + A M • -→ v i .

Test neighborhood

The algorithm 2 begins with an empty test neighborhood T and seek direct covering relations in balls of increasing radii like the method in [START_REF] Rémy | Medial axis for chamfer distances: computing look-up tables and neighbourhoods in 2D or 3D[END_REF]. However, in [START_REF] Rémy | Medial axis for chamfer distances: computing look-up tables and neighbourhoods in 2D or 3D[END_REF], the covering relations are seen from the perspective of the covering balls (in the distance map, requiring propagation of distances), whereas in our case, they are considered from the point of view of the covered balls, i.e. in the covering function that can be computed directly without propagation using the direct formula (22). The algorithm uses the lattice of cones L (fig. 4), resulting from the triangulation phase, as a partition of G. Each cone that is not a covering cone is recursively partitioned by procedure visitCone as depicted in fig. 5: a nD cone C o,U is divided in 2 n subcones corresponding to all the subsets of U . Direct covering relations at the vertices of all the visited cones are checked by procedure visitPoint using a direct computation of the covering function in all know neighbors in T .

Input: L lattice of simplicial unimodular cones, maximal radius R T ← ∅;

Results

An implementation of these algorithms was developed in the C language. It produces output in the same format as the reference algorithm [START_REF] Rémy | Medial axis for chamfer distances: computing look-up tables and neighbourhoods in 2D or 3D[END_REF] so that outputs can be compared character-to-character. Tests were done on various chamfer masks and different maximal radii and helped discover a propagation issue in the reference algorithm. This problem corrected, the results are almost always identical except for insignificant cases close to the maximal radius for which covering radii exceed the maximum. These radii are handled differently by both algorithms but since they exceed the maximum radius, there is no impact on the medial axis computation. The run times of both reference and proposed algorithms are given in table 1 for various sizes and a few masks. These figures show that the overall mean computational complexity is linear with the maximal radius for the proposed method whereas it is of order r 4 in 3D for the reference algorithm.

Conclusion

In this paper, methods to compute both the chamfer LUT and the test neighborhood were presented. Speed gains from the reference algorithm [START_REF] Rémy | Medial axis for chamfer distances: computing look-up tables and neighbourhoods in 2D or 3D[END_REF] are attributable to the representation of chamfer balls as H-polytopes. With this description, we avoid the use of weight propagation in the image domain and

Figure 1 :

 1 Figure 1: (a) Test neighborhood and (b) lookup table for distance d 5,7,11 extracted from [7, fig. 14]. The value Lut-→ c (14) = 23 (in boldface) means that B < (O, 14) ⊆ B < (O + -→ c , 23) as shown in (c) but B < (O, 14) ⊆ B < (O + -→ c , 22) because of the point (-2, 1) (d). In terms of closed balls, B ≤ (O, 13) ⊆ B ≤ (O + -→ c , 22) but B ≤ (O, 13) ⊆ B ≤ (O + -→ c , 21).

Figure 2 :

 2 Figure2: H-representations of a discrete G-symmetrical polytope P (restricted to the first octant). Dashed lines: a H-representation of P . Thick lines: Hrepresentation of P . In the H case, the three equalities are verified for the same point[START_REF] Borgefors | Centres of maximal discs in the 5-7-11 distance transforms[END_REF][START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF]. Notice that although coefficient values are minimal, this representation is still redundant: the second inequality could be removed.

  simultaneously a lower and upper bound to d(O, p) then the equality holds and d is homogenous in C O,U . Conversely, consider a vector -→ u of U such that d(O, p) = p • -→ F , so d(O, p) > p • -→ F , by proposition 2. In the facet F, we always can find n = dim(F) + 1 vertices -→ v 1 /w 1 . . . -→ v n /w n such that -→ u lies in the simplicial cone spanned by the vectors -→ v 1 . . . -→ v n . As a consequence -→ u can be represented by a unique non negative combination of n linearly independent vectors: -

AProposition 6 (

 6 i A Mi • p , thus: max p∈P d(O, p) = max p∈P max 1≤i≤l Mi • p = max 1≤i≤l H i (P ) .P R E P R I N T Minimal covering ball). The radius of the minimal ball centered in O that contains all points of a discrete H-polytope P represented by the matrix A M and the vector H(P ) is equal to the greatest component of H(P ).

  r) (p)) always holds by definition of C B ≤ (O,r) . Proposition 7 is used to limit the search for balls directly covering B ≤ (O, r) in the covering function C B ≤ (O,r) . When a covering cone is detected, only its vertex o or the vectors of the cone (if o = O) have to be tested. P R E P R I N T Proposition 8. If there is an integer j ∈ [1 . .

2 while L is not empty do 3 pop 5 S

 235 (p, q) from 4 take r in S \ {p, q} ← S \ {p, q, r} 6 for s in S do 7 if s is above (p, q, r) then 8 created edge (r, p) and (q, r) in or remove it from 12 if already contained end

13

 13

Figure 3 :

 3 Figure 3: The rational ball B R of the 3D chamfer mask 7,8,11,14 and a unimodular triangulation of B R ∩ G. (b) Unimodular triangulation in 3D cones (represented by triangles) of B R ∩ G (thick lines in (a)). Each Z 3 point (x, y, z) is projected onto (y/x, z/x) (Farey space).

Figure 4 :

 4 Figure 4: Lattice of cone vector sets for the unimodular triangulation of B R ∩ G (chamfer mask 7,8,11,14).

Figure 5 :

 5 Figure 5: The partition of a 3D unimodular cone in 2 3 subcones (from darkest to lightest): 1 point, C 3 1 = 3 1D cones, C 3 2 = 3 2D cones and C 3 3 = 1 3D cone.

1 for r ← 1 3 o 5 end 6 end 7 Algorithm 2 :

 1135672 to R do // inner radius 2 foreach C O,U ∈ L do // cones with increasing dimension ← O + -→ u ∈U -→ u ; 4 visitCone(T , C o,U , r); Computation of the test neighborhood

P 1 / 3 o 5 end 6 end 7 Procedure 1 T 3

 1356713 R E P R I N T Input: T , C o,U , inner radius r, y = H(B ≤ (o, r)) visitPoint (o) / Visit subcones if C o,U is not a covering cone if ∀j, ∀ -→ u ∈ U : y j = max k y k and A Mj • -→ u = max k A Mk • -→ u then 2 foreach subcone C o,U of C o,U do ← o + -→ u ∈U -→ u ; 4 visitCone(T , C o ,U , r); visitCone(T , C o,U , r) Input: T , inner radius r, ball center p to test if ∀ -→ v ∈ T : B ≤ (O + -→ v , C B ≤ (O,r) (O + -→ v )) ⊂ B ≤ (p, C B ≤ (O,r) (p)) thenProcedure visitPoint(T , r, p)

  ]. For instance, with d 14,20,31,44 , Lut (2,1) (291) = 321 and Lut (2,1) (321) = 352, but the smallest open ball of center O that covers B < ((4, 2), 291) is B <

Table 1 :

 1 Run times (in seconds) of reference (ref.) and proposed (H) method implementations for four 3D masks and various volume sizes L × L × L.

		3, 4, 5	3, 4, 5, 7	7, 8, 11, 14	4, 6, 7, 9, 10
	L 10	ref. 0.0006 0.0003 0.0002 0.0004 0.0013 0.0010 0.0004 0.0008 H ref. H ref. H ref. H
	20	0.0049 0.0002 0.0032 0.0004 0.0135 0.0011 0.0047 0.0010
	50	0.1079 0.0003 0.1096 0.0005 0.2886 0.0027 0.1592 0.0015
	100 1.6316 0.0005 1.6831 0.0008 4.4704 0.0054 2.4626 0.0023
	200 30.391 0.0010 31.126 0.0012 81.627 0.0109 53.126 0.0038
	500	3523	0.0022	3537	0.0024	9333	0.0272	6465	0.0088
	1000		0.0042		0.0073		0.0549		0.0220

Acknowledgment

The authors would like to thank the referees for their valuable suggestions and comments, which greatly improved the paper.

P R E P R I N T obtain a constant time covering test by the direct computation of covering radii. Moreover the search space is greatly reduced using covering cones.

While applications always using the same mask can use precomputed test neighborhood T and LUT, other applications that potentially use several masks, adaptive masks, variable input image size can benefit from these algorithms. A faster computation of T is also highly interesting to explore chamfer mask properties. Beyond improved run times, the H-polytope representation helped to prove new properties of chamfer masks. And a new formula of distance which doesn't need to find in which cone lies a point was given.

Examples and source codes for the 3D and 2D cases are available by simple email request to the authors or in the code section of the IAPR-TC18 web site (http://www.cb.uu.se/ ~tc18/).