Jacques-Olivier Lachaud
email: jacques-olivier.lachaud|xavier.provencal@univ-smb.fr

Xavier Provençal

Tristan Roussillon
email: tristan.roussillon@liris.cnrs.fr

Computation of the normal vector to a digital plane by sampling signicant points

Keywords: Digital geometry, digital plane, recognition, normal vector estimation, lattice reduction

Digital planes are sets of integer points located between two parallel planes. We present a new algorithm that computes the normal vector of a digital plane given only a predicate is a point x in the digital plane or not. In opposition with the algorithm presented in [7], the algorithm is fully local and does not explore the plane. Furthermore its worst-case complexity bound is O(ω), where ω is the arithmetic thickness of the digital plane. Its only restriction is that the algorithm must start just below a Bezout point of the plane in order to return the exact normal vector. In practice, our algorithm performs much better than the theoretical bound, with an average behavior close to O(log ω). We show further how this algorithm can be used to analyze the geometry of arbitrary digital surfaces, by computing normals and identifying convex, concave or saddle parts of the surface.

Introduction

The study of the linear geometry of digital sets has raised a considerable amount of work in the digital geometry community. In 2D, digital straightness has been extremely fruitful. Indeed, properties of digital straight lines have impacted the practical analysis of 2D shapes, especially through the notion of maximal segments [START_REF] Doerksen-Reiter | Convex and concave parts of digital curves[END_REF][START_REF] De Vieilleville | Maximal digital straight segments and convergence of discrete geometric estimators[END_REF], which are unextensible pieces of digital straight lines along digital contours.

In 3D, the main problem is that there is no more an implication between being a maximal plane and being a tangent plane as it is in 2D. This was highlighted in [START_REF] Charrier | Maximal planes and multiscale tangential cover of 3d digital objects[END_REF], where maximal planes were then dened as planar extension of maximal disks. To sum up, the problem is not so much to recognize a piece of plane, but more to group together the pertinent points onto the digital shape. This work has been partly funded by DigitalSnow ANR-11-BS02-009 research grant A rst algorithm was proposed by the authors in [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF]. Given only a predicate is point x in plane P ? and a starting point in P, this algorithm extracts the exact characteristics of P. The idea is to deform an initial unit tetrahedron based at the starting point with only unimodular transformations. Each transformation is decided by looking only at a few points around the tetrahedron. This step is mostly local but may induce sometimes a dichotomic exploration.

At the end of this iterative process, one face of the tetrahedron is parallel to P, and thus determines its normal. The remarkable idea is that the algorithm decides itself on-the-y which points have to be considered for computing the local plane geometry, in opposition to usual plane recognition algorithms [START_REF] Kim | On the recognition of digital planes in threedimensional space[END_REF][START_REF] Veelaert | Digital planarity of rectangular surface segments[END_REF][START_REF] Gérard | An elementary digital plane recognition algorithm[END_REF][START_REF] Charrier | An ecient and quasi linear worst-case time algorithm for digital plane recognition[END_REF].

This approach was thus also interesting for analysing 3D shape boundaries. This paper proposes another algorithm that extracts the characteristics of plane P given only this predicate and a starting conguration. This new algorithm is also an iterative process that deforms an initial tetrahedron and stops when one face is parallel to P. However, this new algorithm both complements and enhances the former approach. It complements it since it extracts triangular facets onto the plane whose Bezout point is above the facet, while the former algorithm extracts the ones whose Bezout point is not above the facet. See g. 2 for an example of execution of both algorithms with the same input. It is also a geometrical algorithm using convex hull and Delaunay circumsphere property, while the former was mostly arithmetic. It enhances it for several reasons. First, its theoretical complexity is slightly better (it drops a log factor) and it is faster also in practice. Second, we control the position of the evolving tetrahedron, which always stays around the starting point. Third, each step is purely local and tests the predicate on only six points. Fourth, it can detect planarity defects onto digital surfaces that are not digital planes. Last, a variant of the proposed algorithm almost always produce directly a reduced basis of the lattice of upper leaning points, without lattice reduction. Of course, this algorithm presents one disadvantage with respect to the former one: the starting conguration must lie at a reentrant corner of the plane, more precisely below the Bezout point of the plane. If it starts at another corner, then the algorithm will stop sooner and outputs only an approximation of the normal of P.

The paper is organized as follows. First, we give basic denitions and present our new algorithm. Second we show its correctness and exhibit worst-case upperbound for its time complexity. Third we study how often this algorithm extracts a reduced basis of the lattice of upper leaning points and present a variant with the same properties that returns almost always a reduced basis. Afterwards we exploit this algorithm to determine the linear geometry of digital surfaces, and we show how to deal with starting congurations that are not under the Bezout point. Finally the pros and cons of this algorithm are discussed and several research directions are outlined.

Algorithm 1: Extracts a triangle aligned with the digital plane P by successive convex hull computations.

Input: a predicate x ∈ P ?, an initial triangle T (0)

1 T ← T (0) ; 2 while Σ T ∩ P = ∅ do 3
Compute the convex hull of T ∪ (Σ T ∩ P) ; 4

Find T , dened as the upper triangular facet intersected by [pq] ; 5 T ← T ; 6 return T;

2 Notations and algorithm

In this section, we introduce a few notations before presenting our new algorithm.

A digital plane is dened as the set

P = {x ∈ Z 3 |µ ≤ x • N < µ + s • N},
where N ∈ Z 3 is the normal vector whose components (a, b, c) are relatively prime, µ ∈ Z is the intercept, s is the shift vector, equal to (±1, ±1, ±1) in the standard case.

By translation, we assume w.l.o.g. that µ = 0. Moreover, by symmetry, we assume w.l.o.g. that the normal of the plane lies in the rst octant, i.e. its components are positive. We also exclude cases where a component is null since then it falls back to a 2D algorithm. Thus, a, b, c > 0 and s = (1, 1, 1). Finally, we denote by ω := s • N = a + b + c the thickness of the standard digital plane P.

The above denition of digital plane suggest to see the space as partitionned into layers of coplanar points, orthogonal to N. The value x • N, called height, is a way of sorting these layers in direction N. Points of height 0 (resp. ω-1), which are extreme in P, are called lower leaning points (resp. upper leaning points ). Points of height ω, the closest ones above P, are called Bezout points. Algorithm 1 nds N given a predicate is x ∈ P ?. The algorithm starts at any reentrant corner as follows: the corner point p is in P; the shifted point q := p + s is not in P since s • N is the thickness of P; the initial triangle T (0) := (v

(0) k ) k ∈{0,1,2} is such that ∀k, v (0) k := p + e k + e k+1 (g. 1.c). It is easy to check that T (0) ⊂ P for any p such that 0 ≤ p • N < min{a, b, c},
which corresponds exactly to reentrant corners onto a standard digital plane.

The algorithm then updates this initial triangle and iteratively aligns it with the plane by calling the above predicate for well-chosen points, until the solution equals N.

At step i ∈ N, the solution is described by a triangle denoted by T (i) . Algo- rithm 1 is designed so that ∀i, T (i) is included in P and is intersected by segment [pq]. Segment [pq] thus forces the exploration to be local.

Let k be an integer taken modulo 3, i.e. k ∈ Z/3Z. The three counterclockwise oriented vertices of T (i) are denoted by v (i) k (g. 1.a). For sake of clarity, we will is dropped for sake of clarity). The starting triangle is illustrated in (c).

p q v 0 v 1 v 2 d 1 d 2 d 0 m 0 m 1 m 2 (a) triangle p q v 0 v 1 v 2 d 1 d 2 d 0 (b) neighborhood • p e 0 e 1 e 2 • q • v 0 0 • v 0 1 • v 0 2 (c) step 0
write ∀k instead of ∀k ∈ Z/3Z. We also introduce the following vectors (g. 1.a):

∀k, m (i) k := q -v (i) k (1) ∀k, d (i) k := v (i) k+1 -v (i) k = m (i) k -m (i) k+1 . (2)
The normal of T (i) , denoted by N(T (i) ), is merely dened by the cross product between two consecutive edge vectors of T (i) , i.e. N(

T (i) ) = d (i) 0 × d (i) 1 .
In order to improve the guess at step i, the algorithm checks if the points of a small neighborhood around q, parallel and above T (i) , belong to P or not. This neighborhood is dened as follows (g. 1.b):

Σ T (i) := q ± d (i) k k ∈{0,1,2} . (3)
The algorithm then computes Σ T (i) ∩ P by making six calls to the predicate. In algorithm 1, the new triangle T (i+1) is simply dened as the other triangular facet of the convex hull of T (i) ∪(Σ T (i) ∩P) which intersects [pq]. See g. 2 for an example. The algorithm stops when Σ T (n) ∩ P is empty. We show in theorem 2 that, when p is a lower leaning point (its height is 0), the vertices of the last triangle T (n) are upper leaning points of P (their height is ω -1). Corollary 1 implies that the normal to T (n) is the normal to P and corollary 2 implies that triangle edges form a basis of the lattice of upper leaning points to P.

Validity and complexity

In this section, we rst prove the two invariants of algorithm 1 and fully characterize its main operation (lines 3-5 of algorithm 1). Then we prove that if p is a lower leaning point, algorithm 1 retrieves the true normal N in less than ω -3 steps.

Let us assume that the algorithm always terminate in a nite number of steps whatever the starting point and let n be the last step (the proof is postponed, see theorem 1).

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (a) i = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (b) i = 1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (c) i = 2 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (d) i = 3 (e) i = 0 (f ) i = 2 (g) i = 5 (h) i = 7

Algorithm invariants

The following two invariants are easy to prove by induction since by construction, each new triangle is chosen as a triangular facet of the convex hull of some subset of P, which is intersected by ]pq].

Invariant

1 ∀i ∈ {0, . . . , n}, ∀k, v (i) k ∈ P.
Invariant 2 ∀i ∈ {0, . . . , n}, the interior of T (i) is intersected by ]pq].

The only diculty in invariant 2 is to show that the triangle boundary is never

intersected by ]pq]. Below, for lack of space, we only show by contradiction that ]pq] does not intersect any edge of T (i+1) , if invariant 2 is assumed to be true for i ∈ {0, . . . , n -1}.

Proof. If ]pq] and some edge [xy] of T (i+1) intersect, points p, q, x and y are coplanar. But, since the boundary of T (i) is not intersected by ]pq] and since points of Σ T (i) are located on lines that are parallel to the sides of T (i) , x and y must be opposite points in Σ T (i) , i.e. x := q + d (i) l

for some l ∈ {0, 1, 2} and y = qd (i) l .

However, since q / ∈ P, x and y cannot be both in P by linearity and thus cannot be the ends of an edge of T (i+1) , which is included in P by invariant 1.

Operation characterization

The two following lemmas characterize the operation that transforms a triangle into the next one (lines 3-5 of algorithm 1). Lemma 1. ∀i ∈ {0, . . . , n -1}, T (i) and T (i+1) share exactly one or two vertices.

Proof. Let γ be the number of common vertices between T (i) and T (i+1) . Since T (i) and T (i+1) have three vertices each, 0 ≤ γ ≤ 3, but we prove below that (i) γ = 3 and (ii) γ = 0.

Inequality (i) is trivial. Indeed, since the volume of the convex hull of T (i) ∪ (Σ T (i) ∩ P) is not empty for i ∈ {0, . . . , n -1} (see algorithm 1, l. 2), we necessarily have T (i) = T (i+1) , which implies that γ = 3.

In order to prove inequality (ii), let P be the plane passing by the points of Σ T (i) (points of Σ T (i) are coplanar by ( 2) and ( 3)). Note that q ∈ P by [START_REF] De Vieilleville | Maximal digital straight segments and convergence of discrete geometric estimators[END_REF].

Let us assume that T (i+1) is included in P. By invariant 1, T (i+1) ⊂ P. However, q / ∈ P by denition of q. As a consequence, the upper leaning plane of P, i.e. {x ∈ R 3 |x • N = ω -1}, strictly separates T (i+1) from q in the plane P. In addition, p is located strictly below T (i) by invariant 2 and is thus also below P, which is parallel and above T (i) by denition. As a consequence, T (i+1) is clearly not intersected by ]pq], which contradicts invariant 2. We conclude that T (i+1) is not included in P, which means that γ = 0.

The following lemma fully characterizes the main operation of algorithm 1:

Lemma 2. ∀i ∈ {0, . . . , n -1}, ∀k, either v (i+1) k = v (i) k , or v (i+1) k = v (i) k + m (i) l , l = k.
We know by lemma 1 that T (i) and T (i+1) share one or two vertices. The proof is thus in two parts. In each part however, we use innite cones of apex q to check whether invariant 2 is true or not. Let us dene ∀i ∈ {0, . . . , n}, T ∞ (i) as the set of innite rays emanating from q and intersecting the triangular facet T (i) . Invariant 2 is equivalent to the following invariant 2': ∀i ∈ {0, . . . , n}, the interior of T ∞ (i) contains p.

v 0 v 1 v 2 1.1 1.3 q 1.1 1.2 1.2 1.3 (a) (T (i) ∩ T (i+1) ) = 2 v 0 v 1 v 2 2.1 2.2 2.3 2.2 q (b) (T (i) ∩ T (i+1) ) = 1
Fig. 3. Illustration of the proof of lemma 2. Case where T (i) and T (i+1) share v (i) 1 and

v (i) 2 in (a), but only v (i) 2 in (b) (exponent (i) is omitted in the gures).
Proof. We rst assume that T (i) and T (i+1) share two vertices. Let us assume w.l.o.g. that the index of the vertex of T (i) that is not in T (i+1) is 0.

Since v (i+1) 0 ∈ Σ T (i) ∩ P, there are six cases for v (i+1) 0

, which are grouped below two by two (see g. 3.a). Invariant 1 is true for all cases. We show however by contradition that the rst four cases (items 1.1 and 1.2) are not possible because otherwise, invariant 2' is not true:

1.1 Let us assume that v

(i+1) 0 = v (i) l + m (i) 0 , l ∈ {1, 2}. In these cases, T ∞ (i+1) is adjacent to T ∞ (i) along facet (q, v (i) 1 , v (i)
2 ). We conclude that the interior of T ∞ (i+1) and T ∞ (i) are disjoint and cannot both contains p, which raises a contradiction.

1.2 If we assume now that v

(i+1) 0 = v (i) 1 + m (i) 2 or v (i+1) 0 = v (i) 2 + m (i) 1 . In these cases, v (i+1) 0
lies on the plane passing by q, v (i)

1 and v (i) 2 . We conclude that the interior of T ∞ (i+1) is empty and cannot contains p, a contradiction.

1.3 Hence, v (i+1) 0 = v (i) 0 + m (i) l , l ∈ {1, 2}. In this case, T ∞ (i+1) , which contains
T ∞ (i) , contains p.

We now assume that T (i) and T (i+1) share one vertex. Let us assume w.l.o.g. that its index is 2.

By denition, v

(i+1) 0 , v (i+1) 1
∈ Σ T (i) ∩ P. There are 6 2 = 15 cases to consider for v

(i+1) 0 and v (i+1) 1
, which are grouped below into four dierent items. See g. 3.b for the last three items. Invariant 1 is true for all cases. We show however that invariant 2 is true for only the three cases described at item 2.3:

2.0 (3 cases) Let us assume that v (i+1) 0 and v (i+1) 1
are opposite with respect to q, i.e. one is equal to q + d

(i+1) k , k ∈ {0, 1, 2}, whereas the other is equal to q -d (i+1) k . Since q / ∈ P, q + d (i+1) k and q -d (i+1) k
cannot be both in P by linearity, which raises a contradiction.

2.1 (7 cases) Let us assume that v

(i+1) 0 = v (i) 2 + m (i) l , l ∈ {1, 2} or v (i+1) 1 = v (i) 2 + m (i) l , l ∈ {1, 2}. In all cases, T ∞ (i+1) is adjacent to T ∞ (i) by edge (q, v (i)
2 ), by facet (q, v

(i) 0 , v (i)
2 ) or by facet (q, v

(i) 1 , v (i)
2 ). As a consequence, the interior of T ∞ (i+1) cannot contains p.

2.2 (2 cases) Let us assume that v

(i+1) 0 = v (i) l + m (i) l+1 and v (i+1) 1 = v (i) l + m (i) l+2 , for l ∈ {0, 1}. In both cases, T ∞ (i+1) is adjacent to T ∞ (i) by edge (q, v (i)
l ).

As a consequence, the previous argument also applies.

2.3 (3 cases) T ∞ (i+1) contains T ∞ (i) and thus contains p in the following cases:

v (i+1) 0 = v (i) 0 + m (i) l0 , l 0 ∈ {1, 2} and v (i+1) 1 = v (i) 1 + m (i) l1 , l 1 ∈ {0, 2}, with (l 0 , l 1 ) = (1, 0).
Lemma 3. Let M (i) be the 3 × 3 matrix that consists of the three row vectors (m

(i) k ) k ∈{0,1,2} . Then, ∀i ∈ {0, . . . , n}, det (M (i) ) = 1.
Proof. We can easily check that det (M (0) ) = 1.

We now prove that if det (M (i) ) = 1 for i ∈ {0, . . . , n -1}, then det (M (i+1) ) = 1. By lemma 1 and lemma 2, we have v

(i+1) k = v (i) k + m (i) l ⇔ m (i+1) k = m (i) k -m (i)
l , l = k for at most two rows over three, while the remaining one or two rows correspond to identity. Such matrix operations do not change the determinant, which concludes.

Termination

In the following proofs, we compare the position of the points along direction N. For the sake of simplicity, we use the bar notation • above any vector x to denote its height relative to N. Otherwise said, x := x • N. Even if N is not known, q ≥ ω by denition and for all x ∈ P, 0 ≤ x < ω.

Theorem 1. The number of steps in algorithm 1 is bounded from above by ω-3.

Proof. First, it is easy to see that ∀k, m

(0) k = e k+2 . Thus, k m (0) k = ω.
Let us recall that ∀i ∈ {0, . . . , n}, ∀k, m

(i) k = q -v (i)
k . We have q ≥ ω by denition of q and ∀i ∈ {0, . . . , n}, ∀k, 0 ≤ v (i) k < ω by invariant 1 and by denition of P. As a consequence, ∀i ∈ {0, . . . , n}, ∀k, m

(i) k > 0 and k m (i) k ≥ 3.
Moverover, by lemma 2 and since any m

(i) k is strictly positive ∀i ∈ {0, . . . , n}, we clearly have ∀i ∈ {0, . . . , n -1}, k m (i) k > k m (i+1) k . The sequence ( k m (i) k ) i=0,.
..,n is thus a strictly decreasing sequence of integers between ω and 3.

Remark that this bound is reached when running the algorithm on a plane with normal N(1, 1, r). We now give the main result when the starting point p is a lower leaning point, i.e. p = 0. Note that in this case q = ω. Since we focus below on the last step n, we omit the exponent (n) in the proofs to improve their readability.

Theorem 2. If p is a lower leaning point (i.e. p = 0 and thus q = ω), the vertices of the last triangle are upper leaning points, i.e. ∀k, v

(n) k = ω -1.
Proof. If there exists k ∈ {0, 1, 2} such that d k = 0, then either (i) d k < 0 or (ii) d k > 0. Since q = ω and |d k | < ω, either (i) q + d k ∈ P or (ii) qd k ∈ P. This implies that Σ T ∩ P = ∅, which is a contradiction because Σ T ∩ P = ∅ at the last step (see algorithm 1, l. 2). As a consequence, ∀k, d k = 0 and ∀k, m k = α, a strictly positive integer.

Let us denote by 1 the vector (1, 1, 1) T . We can write the last system as MN = α1. Since M is invertible (because det (M) = 1 by lemma 3), N = M -1 α1 and as a consequence α = 1 (because the components of N are relatively prime and M -1 is unimodular).

We conclude that ∀k, m k = 1 and, straightforwardly, v k = ω -1.

The following two corollaries are easily derived from lemma 3 and theorem 2.

Corollary 1. If p is a lower leaning point, the normal of the last triangle is equal to N, i.e. N(T (n) ) = N.

Corollary 2. If p is a lower leaning point, (

) is a basis of the lattice of upper leaning points {x ∈ P|x = ω -1}.

Lattice reduction and Delaunay condition

Corollary 2 implies that the edge vectors of the last triangle form a basis of the lattice of upper leaning points to P.

Let us recall that a basis (x, y) is reduced if and only if x , y ≤ x-y ≤ x + y . Given (x, y) a basis of a two dimensional lattice, there is an iterative algorithm to compute a reduced basis from (x, y). This algorithm consists in replacing the longest vector among (x, y) by the shortest one among x + y and xy, if it is smaller.

We have found many examples for which the basis returned by algorithm 1 is not reduced, even if we take the two shortest vectors in {d (n) k } k ∈{0,1,2} (see g. 4.a for an example). To cope with this problem, we can either apply the above reduction algorithm to the returned basis, or keep triangles as small as possible so that the last triangle leads to a reduced basis, without any extra reduction steps. To achieve this aim, we use a variant of algorithm 1; see algorithm 2.

Note that in algorithm 2, we compute the convex hull of T ∪ S , but since S ⊂ Σ T ∩ P and since Σ T ∩ P = ∅ ⇒ S = ∅, all the results of sec. 3 remain valid.

In g. 4, the results of the two algorithms for the digital plane of normal [START_REF] Charrier | Maximal planes and multiscale tangential cover of 3d digital objects[END_REF][START_REF] De Vieilleville | Maximal digital straight segments and convergence of discrete geometric estimators[END_REF]9) are compared. The triangles returned by algorithm 2 are more compact and closer to segment [pq]. The last triangle returned by algorithm 2 leads to a reduced basis, but this is not true for algorithm 1.

Algorithm 2: Variant of algorithm 1 based on a Delaunay condition.

Input: a predicate x ∈ P ?, an initial triangle T (0)

1 T ← T (0) ; 2 while Σ T ∩ P = ∅ do 3
Compute the set S of points s ∈ (Σ T ∩ P) such that the circumsphere of T ∪ {s } does not include any point s ∈ (Σ T ∩ P) in its interior ; 4

Compute the convex hull of T ∪ S ; 5

Find T , dened as the upper triangular facet intersected by [pq] ; 6 T ← T ; 7 return T; (1, 1, 1) to (200,200,200). There are 6578833 vectors with relatively prime components in this range. It turns out that less than 0.01% of the basis returned by algorithm 2 are non-reduced against about 73% for algorithm 1.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (a) alg. 1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (b) alg. 2

Applications to digital surfaces

In this section, we consider a set of voxels, Z, where voxels are seen as unit cubes whose vertices are in Z 3 . The digital boundary, BdZ, is dened as the topological boundary of Z. Since a digital boundary looks locally like a digital plane, it is natural to run our algorithm at each reentrant corner of the digital boundary with predicate Is x in BdZ in order to estimate the local normal vector to the volume Z (see g. 5).

Although we do not go into details here due to space limitations, we can mention the following points:

Our algorithm works well on digitally convex shapes Z. The set of points

x ∈ BdZ located below some facet F of Conv(BdZ) whose normal vector is in the rst octant and such that points x + s are strictly above F, is a piece of a digital plane, P. If the starting reentrant corner is a lower leaning point of P, the last triangle computed by our algorithm is included in F (see g. 5). We call such triangles patterns of P.

When starting from a corner that is not a lower leaning point of a facet of Conv(BdZ), the algorithm returns a triangle called reduced pattern of P which is approximately aligned with the plane.

It is possible to detect and remove reduced patterns with a simple algorithm.

Put in queue every point q of all computed triangles, together with their three vectors m k . Loop on the queue while it is not empty by popping them in sequence. For each popped point q and vectors m k , check if q translated by any of m k is a shifted point q of another triangle. In this case, mark q as reduced pattern and put it back in queue but now with the vectors m k .

Otherwise, do nothing. When the algorithm stops, since reduced patterns have shifted points above the shifted point of the lower leaning point, all reduced patterns are marked.

If the shape is not convex, the algorithm can be adapted to detect planarity defects. The idea is that Σ T ∩ BdZ should contain at most three elements if it is locally a piece of plane. Moreover these elements must always be consecutive neighbors around q in Σ T . We thus stop the algorithm whenever at least one of the two previous conditions fails, meaning that the surface is locally non convex. We also check that each new triangle is separating onto the digital boundary: for any point x ∈ BdZ below a triangle of shift s, the point x + s must be above the triangle.

Conclusion and perspectives

In this paper, we proposed a new algorithm that computes the parameters of a digital plane. In opposition to usual plane recognition algorithms, this algorithm greedily decides on-the-y which points have to be considered like in [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF]. Compared to [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF], this algorithm is however simpler because it consists in iterating one geometrical operation, which is fully characterized in lemma 2. In addition, the returned solution, which is described by a triangle parallel to the sought plane, is close to the starting point. On one hand, the starting point must projects onto the triangle and on the other hand, the two shortest edge vectors of the triangle almost always form a reduced basis (in the second version of our algorithm).

For the future, we would like to nd a theoretically sound way of always getting a reduced basis, without any extra reduction operation. Moreover, we would like to nd a variant of our algorithm in order to retrieve complement triangles whose Bezout point is not above the triangle. For sake of completeness, we are also interested in degenerate cases, where at least one component is null.

After having achieved these goals, we would have a complete working tool for the analysis of digitally convex surfaces. In order to correctly process concave or saddle parts however, we must precisely associate to a given triangle a piece of digital plane containing the triangle vertices and having the same normal vector than the triangle one. Although there are many such pieces of digital plane for a given triangle, we hope that this work will provide a canonical hierarchy of pieces of digital plane.

Fig. 1 .

 1 Fig.1. Main notations are illustrated in (a) and (b) (iteration number as exponent (i)

Fig. 2 .

 2 Fig.2. Illustration of the running of algorithm 1 and algorithm from[START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF] on a digital plane of normal vector[START_REF] Charrier | An ecient and quasi linear worst-case time algorithm for digital plane recognition[END_REF][START_REF] Charrier | Maximal planes and multiscale tangential cover of 3d digital objects[END_REF][START_REF] Gérard | An elementary digital plane recognition algorithm[END_REF]. Images (a) to (d) shows the four iterations of algorithm 1 starting from the origin. For i ∈ {0, 1, 2, 3}, triangles T (i) are in blue, whereas neighborhood points are depicted with red disks (resp. circles) if they belong (resp. do not belong) to P. In (d), the normal of the last triangle is[START_REF] Charrier | An ecient and quasi linear worst-case time algorithm for digital plane recognition[END_REF][START_REF] Charrier | Maximal planes and multiscale tangential cover of 3d digital objects[END_REF][START_REF] Gérard | An elementary digital plane recognition algorithm[END_REF]. Images (e) to (h) shows iterations 0 (initial), 2, 5 and 7 (nal) of the algorithm from[START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF]. The initial tetrahedron (a) is placed at the origin and the nal one (h) has an upper triangle with normal vector[START_REF] Charrier | An ecient and quasi linear worst-case time algorithm for digital plane recognition[END_REF][START_REF] Charrier | Maximal planes and multiscale tangential cover of 3d digital objects[END_REF][START_REF] Gérard | An elementary digital plane recognition algorithm[END_REF].

Fig. 4 .

 4 Fig.4. Triangles returned by algorithm 1 in (a) and algorithm 2 in (b) for a digital plane with normal (2, 3, 9) (the corner point p is set to the origin). In (b), triangles returned by algorithm 2 are smaller and the last one denes a reduced basis.

Fig. 5 .

 5 Fig.5. Our algorithm has been runned at each reentrant corner of a digital plane (left) and an ellipse (right). The last triangle of each run is printed.


Jacques-Olivier Lachaud

Xavier Provençal

Tristan Roussillon

.

Two plane-probing algorithms for the computation of the normal vector to a digital plane *

Introduction

The study of the linear geometry of digital sets has raised a considerable amount of work in the digital geometry community. In 2D, digital straightness has been extremely fruitful. Indeed, digital straight lines present many properties, whether geometric, arithmetic, or combinatoric (e.g. see survey [START_REF] Klette | Digital straightness -a review[END_REF]). Furthermore, these properties have impacted the practical analysis of 2D shapes, especially through the notion of maximal segments [START_REF] Doerksen-Reiter | Convex and concave parts of digital curves[END_REF][START_REF] De Vieilleville | Maximal digital straight segments and convergence of discrete geometric estimators[END_REF], which are unextensible pieces of digital straight lines along digital contours. To sum up, they were shown to be characteristic of convex and concave parts [START_REF] Feschet | Canonical representations of discrete curves[END_REF][START_REF] Doerksen-Reiter | Convex and concave parts of digital curves[END_REF][START_REF] Roussillon | Faithful polygonal representation of the convex and concave parts of a digital curve[END_REF], to induce several multigrid convergence results [START_REF] De Vieilleville | Maximal digital straight segments and convergence of discrete geometric estimators[END_REF][START_REF] Lachaud | Fast, accurate and convergent tangent estimation on digital contours[END_REF] and even to be able to identify noisy parts of the contour [START_REF] Kerautret | Meaningful scales detection along digital contours for unsupervised local noise estimation[END_REF].

A lot was thus expected from the study of 3D digital planes (e.g. see the survey [START_REF] Brimkov | Digital planarity -a review[END_REF]). It is true that several plane recognition algorithms were proposed to determine if a given set of points could be a piece of digital plane. Among them, the most fruitful ones adopt a geometric approach [START_REF] Kim | On the recognition of digital planes in threedimensional space[END_REF][START_REF] Veelaert | Digital planarity of rectangular surface segments[END_REF][START_REF] Gérard | An elementary digital plane recognition algorithm[END_REF][START_REF] Charrier | An efficient and quasi linear worst-case time algorithm for digital plane recognition[END_REF]. A first difficulty is the more complex combinatorial structure of digital planes [START_REF] Labbé | A d-dimensional extension of christoffel words[END_REF], which appears when studying its connectedness [START_REF] Jamet | Minimal arithmetic thickness connecting discrete planes[END_REF] or the link with multidimensional continued fractions [START_REF] Fernique | Generation and recognition of digital planes using multidimensional continued fractions[END_REF][START_REF] Berthé | Brun expansions of stepped surfaces[END_REF]. Their multiple definitions become obvious when looking at the geometry of digital planes: there are multiple ways of approximating them depending on the chosen point of view.

These ambiguities make practical 3D digital shape analysis difficult. The segmentation of the shape boundary into linear parts is thus generally greedy [START_REF] Klette | Digital planar segment based polyhedrization for surface area estimation[END_REF][START_REF] Provot | 3D noisy discrete objects: Segmentation and application to smoothing[END_REF][START_REF] Chica | Pressing: Smooth isosurfaces with flats from binary grids[END_REF], for instance by extracting first the biggest plane, and then repeat the process [START_REF] Zrour | Optimal consensus set for digital line and plane fitting[END_REF]. In 3D, the main problem is that there is no more an implication between "being a maximal plane" and "being a tangent plane" as it is in 2D. This was highlighted in [START_REF] Charrier | Maximal planes and multiscale tangential cover of 3d digital objects[END_REF], where maximal planes were then defined as planar extension of maximal disks. Contrary to the 2D case, the surface topology around the point of interest does not give us sufficient constraints to identify unambiguously the set of points that should define the tangent plane. To sum up, the problem is not so much to recognize a piece of plane, but more to group together the pertinent points onto the digital shape.

We need thus methods that identify locally significant points to test and extract plane parameters at the same time. From now on, we have as input a starting point in some plane P and a predicate that answers to the question "is point x in plane P ?". The objective is to find the exact plane parameters solely from this information, by testing points as locally as possible. We call plane-probing algorithms that class of plane recognition algorithms.

We proposed a first approach to solve this problem in [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF]. Its principle is to deform an initial unit tetrahedron based at the starting point with only unimodular transformations. Each transformation is decided by looking mostly at a few points around the tetrahedron. These points are chosen so that the transformed tetrahedron is in P, with the same volume, and is pushed towards one side of the plane. At the end of this iterative process, one face of the tetrahedron has an extremal position in the plane and is thus parallel to P.

In [START_REF] Lachaud | Computation of the normal vector to a digital plane by sampling significant points[END_REF], we designed a new algorithm for this problem. It shares some features of the previous one, because it is also an iterative process that deforms an initial tetrahedron and stops when one face is parallel to P. It differs from it on several points, as illustrated on fig. 1. One vertex of the evolving tetrahedron is the fixed point lying above the starting point and the opposite triangular facet. The position of the evolving tetrahedron is thus better controlled than in [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF]. Moreover, this new algorithm is mostly a geometrical algorithm using Delaunay circumsphere property, while the former was mostly arithmetic. Its theoretical complexity is slightly better in theory, since it drops a log factor, and in practice also. This paper extends in several ways the DGCI paper [START_REF] Lachaud | Computation of the normal vector to a digital plane by sampling significant points[END_REF]:

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (a) i = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
• • • • • • • (b) i = 1
• • • • • • • (c) i = 2 • • • • • • • • • • • • • • • • • • • • • • • • • • •
• • • • • • • (d) i = 3 (e) i = 0 (f) i = 2 (g) i = 5 (h) i = 7
• The H-algorithm is essentially the same algorithm as Algorithm 2 of [START_REF] Lachaud | Computation of the normal vector to a digital plane by sampling significant points[END_REF].

We give a new variant, called R-algorithm, for finding the next tetrahedron, which guarantees that, at the end of the process, the output facet forms a reduced basis of the plane P. For now, we observe this property but we are not able to fully prove it.

• The running time of this new variant is similar to the previous algorithms, and even better for big ω.

• We detail how this algorithm can be used for digital shape analysis. We show how to detect convex/concave/inflexion zones onto a digital surface.

• For some well-identified starting points, this algorithm stops and outputs only an approximation of the normal to P. We show how to detect such bad starting points and how to connect them to their corresponding facet.

• We give a comprehensive experimental evaluation of its time complexity and compare it extensively with the plane-probing algorithm of [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF].

The paper is organized as follows. First, we give basic definitions and present the Hand the R-algorithms. Then we show their correctness and provide worstcase time complexity. We then conduct an experimental evaluation, which shows that the new variant always output a reduced basis. The last part describes how this algorithm can be used to analyze the linear geometry of digital surfaces.

Two plane-probing algorithms

We introduce a few notations before presenting in a unified manner our new plane-probing algorithms: the Hand the R-algorithms.

Digital plane and initialization

We wish to extract the parameters of an arbitrary standard digital plane P, defined as the set

P = { x ∈ Z 3 | µ ≤ x • N < µ + s • N},
where N ∈ Z 3 is the normal vector whose components (a, b, c) are relatively prime, µ ∈ Z is the intercept, s is the shift vector. In the standard case, the shift vector is equal to (±1, ±1, ±1), where the sign of the components are chosen so that the thickness ω := s • N is equal to |a| + |b| + |c|. Moreover, we extract a basis of P, that is a pair of vectors that forms a basis of the 2D lattice

{ x ∈ Z 3 | x • N = 0}.
A basis of a two dimensional lattice is reduced if and only if it consists of the two shortest non-zero lattice vectors.

By translation, we assume w.l.o.g. that µ = 0. Moreover, by symmetry, we assume w.l.o.g. that the components of the normal vector are all positive. We also exclude cases where a component is null since then it falls back to a 2D algorithm. Thus, a, b, c > 0, which implies that s = (1, 1, 1).

We may see the space as partitioned into layers of coplanar points, orthogonal to N. The height x • N sorts these layers along direction N. Points of height 0 and ω -1 are extreme within P, and are called lower leaning points and upper leaning points respectively.

We propose an algorithm that, given a predicate "is x ∈ P ?", computes the normal vector of a piece of digital plane surrounding a starting point p ∈ P.

The algorithm places an intial sequence of 3 points

T (0) := ( v (0) k ) k ∈{0,1,2} such that ∀k, v (0)
k := p + e k + e k+1 , (for sake of clarity, we write ∀k instead of ∀k ∈ Z/3Z) and requires that T (0) ⊂ P.

• p e0 e1 e2 • q • v 0 0 • v 0 1 • v 0 2
It is easy to check that T (0) ⊂ P for any p such that 0 ≤ p • N < min{a, b, c}, which corresponds to points lying inside reentrant corners (see fig. on the right). The algorithm then iteratively updates this initial point set by calling the above predicate for well-chosen points. We explain in the next subsection how the shifted point q := p + s, which is not in P because s • N is the thickness of P, is used to select those points.

Iteration and termination

At any step i ∈ N, the tangent plane surrounding the starting point is described by a sequence of 3 points denoted by T (i) .

For any finite point set S, let us denote by conv(S) its convex hull and aff(S) its affine hull. Thus, conv(T (i) ) is a triangle, whose three counterclockwise oriented vertices are denoted by ( v 2), whereas aff(T (i) ) is the plane passing by these vertices.

(i) k ) k ∈{0,1,2} (see fig.
In order to update T (i) , the algorithm checks whether the points of some neighborhood around q and above aff(T (i) ) belong to P or not. Before defining two such neighborhoods, let us introduce the following notation:

∀k, m (i) k := q -v (i) k . (1)
The following two neighborhoods are some subsets of a cone of apex qk m

k and base { m

(i)
k } k ∈{0,1,2} (see fig. 2). First, let us define the H-neighborhood at step i as follows:

N (i) H := q -m (i) σ(0) + m (i) σ(1) } σ, a permutation over {0,1,2} , (2)
In other words, N

H := q ± ( m (i) k -m (i) k+1 ) k ∈{0,1,2} (i)
. This set consists of six points arranged in an hexagon and that is why it is called H-neighborhood (see red disks in fig. 2). We define below the R-neighborhood, using the notion of ray, with:

R (i) σ := q -m (i) σ(0) + m (i) σ(1) + λ m (i) σ(2) ) λ≥0 , (3)
N (i) R := R (i) σ } σ, a permutation over {0,1,2} . (4)
Looking at the first definition of H-neighborhood, it is clear that the R-neighborhood contains it and extends it along rays that go out of the hexagon, hence the name R-neighborhood (see green squares in fig. 2). We thus have

N (i) H ⊂ N (i) R . Let N (i) be any neighborhood in {N (i) H , N (i) R }.
Our algorithm selects a point x of N (i) ∩ P such that the circumsphere of T (i) ∪ { x } does not include any point of N (i) ∩ P in its interior. The upper triangular facet of the tetrahedron conv(T (i) ∪ { x }) intersected by the straight line passing by p and q is the new triangle conv(T (i+1) ).

Let us introduce the edge vectors of conv(T (i) ):

∀k, d (i)
k := v (i) k+1 -v (i) k = m (i) k -m (i) k+1 . ( 5
) q v 0 m 0 v 1 m 1 v 2 m 2
H , called the hexagon at step i, is depicted with red disks, whereas the set

N (i)
R is depicted with green squares, located along rays coming out of the hexagon vertices. We can see that

N (i) H ⊂ N (i)
R . (Iteration number is dropped for sake of clarity).

The normal of aff(T (i) ), denoted by N(T (i) ), is merely defined by the cross product between two consecutive edge vectors of conv(T (i) ), i.e. N(

T (i) ) = d (i) 0 × d (i) 1 .
The algorithm stops at a step n, when N (n) has an empty intersection with P. The output of the algorithm is T (n) . We will prove in sec. 3.3 that if p is a lower leaning point, then the points of T (n) are upper leaning points of P and that the normal N(T (n) ) of the triangle conv(T (n) ) is aligned with the normal N of P.

Unified presentation of plane-probing algorithms

Algorithm 1 summarizes our two variants for recognizing on-the-fly a digital plane. The predicate "is x in P ?" is used to compute the intersection between the neighborhood and P (see lines 2 and 3). Clearly, Algorithm 1 remains around its starting point by construction, since every triangle has a non-empty intersection with the straight line passing by p and q. It also stays as local as possible with the empty circumsphere property.

Algorithm 1: Unified plane-probing algorithm: for any neighborhood definition (H or R), it extracts a 3-point sequence of upper leaning points Input: a shift vector s, a point q and an initial 3-point sequence T (0) i ← 0 ;

1 while N (i) ∩ P = ∅ do 2
Compute a point x ∈ (N (i) ∩ P) such that the circumsphere of 3 T (i) ∪ { x } does not include any point x ∈ (N (i) ∩ P) in its interior ; Find T (i+1) , defined as the vertex sequence of the upper facet of 4 conv(T (i) ∪ { x }) intersected at a single point by the straight line of direction s passing by q ; i ← i + 1 ;

5 return T (i) ; 6
The algorithm variant using the H-neighborhood is called the H-algorithm, whereas the one using the R-neighborhood is called the R-algorithm. We show in the next section that both variants always extract the exact parameters of plane P in O(ω) iterations in worst cases. As we will see in the experimental section, the H-algorithm often extracts the reduced basis of P, while the Ralgorithm always extracts the reduced basis of P in our experiments. In most cases, it falls back to the H-algorithm, but sometimes it looks for points further away along rays that go out of the hexagon.

To end, note that the H-algorithm is slightly different from [16, Algorithm 2] because only one point of the hexagon is selected at each iteration, instead of possibly several ones in case of cospherical points in [START_REF] Lachaud | Computation of the normal vector to a digital plane by sampling significant points[END_REF]Algorithm 2]. This choice leads to a simpler characterization (by lemma 1-item 5, two consecutive triangles must share exactly two vertices, instead of one or two in [16, lemma 1]), while the number of iterations is not changed (O(ω) in both case). The result is even better, since the H-algorithm more often leads to a reduced basis of upper leaning points than [START_REF] Lachaud | Computation of the normal vector to a digital plane by sampling significant points[END_REF]Algorithm 2] (for all 6578833 vectors with relatively prime components ranging from (1, 1, 1) to (200,200,200), the H-algorithm returns 480 non reduced basis against 924 with [16, Algorithm 2], see sec. 4).

Implementation details

In algorithm 1, line 3, at a step i, we have to compute a point x ∈ (N (i) ∩ P) such that the sphere passing by T (i) ∪ { x } does not contain any other point x ∈ (N (i) ∩ P) in its interior. We say below that x is a closest point to T (i) , since the sphere passing by T (i) ∪ { x } has minimal radius (over the set of spheres passing by

T (i) ∪ { x}).
Searching for a closest point to T (i) is trivial and in O( 1) for the H-algorithm, because the H-neighborhood is finite and its intersection with P has at most 6 points. The algorithm is then similar to finding the minimum element of a sequence: we take an arbitrary point y ∈ (N (i) ∩P) as a current closest point and for each remaining point x of this set, if the sphere passing by T (i) ∪ { y} strictly contains x, then x becomes the new current closest point (see algorithm 2). Searching for a closest point to T (i) is more tricky for the R-algorithm, because the R-neighborhood is unbounded. Two questions may be raised: how to compute its intersection with P and how to determine a closest point without visiting all the points of the intersection?

We omit the iteration exponent (i) below to simplify the presentation. Algorithm 3 processes each ray of the R-neighborhood separately, because over one ray, we can efficiently find a closest one to T in P and add it to the candidate point set. Then, we compute a closest point over the candidate point set by algorithm 2 in O(1) since this set contains at most 6 points (one per ray).

Along a ray, we have two different tasks to perform: find the points that belong to P and over such points, find a closest one to T.

First, the starting point of a given ray is guaranteed to be in P (see algorithm 3, line 3). Second, knowing that q / ∈ P and T ⊂ P (see lemma 1-item 1) and by equations ( 1) and ( 3), there always exist some points of a ray that are not in P. Third, due to the convexity of P, all the points not in P must follow the points in P on a given ray (see lemma 7). Locating the furthest point in P Algorithm 3: Search for a closest point to a 3-point sequence T over the R-neighborhood.

Input: a 3-point sequence T and a shifted point q CandidateP ointSet ← ∅ ;

1

Compute the set of rays {R σ } σ, a permutation over {0,1,2} from T and q ;

2 foreach ray R σ do 3 if the starting point of R σ is in P then 4 x ← ClosestPointInRay(T, R σ ) ; 5 CandidateP ointSet ← CandidateP ointSet ∪ { x}; 6 return ClosestPointInSet(T, CandidateP ointSet) ; 7
on a ray can be done in two stages. In the first stage, one advance in the ray direction by doubling the step at each iteration until a point not in P is found. The last two points determine a range whose lower bound is in P, but upper bound is not in P. In the second stage, a binary search is performed on this range. The overall time complexity is logarithmic in the number of points in P on a given ray. Finding a closest point to T may be done in a similar way. Indeed, the function that maps a point x on the ray to the radius of the sphere passing by T ∪ { x} is convex and has a global minimum, because the ray does not intersect aff(T) (see lemma 1-item 3 and corollary 2).

In fig. 3, we provide a two-dimensional illustration of this function. For each point x on the ray, depicted with green boxes, the spheres passing by T ∪ { x} are depicted with blue circles. We can see that in the ray direction, the function is first decreasing, reaches its minimum at the middle point before increasing.

Thus, in algorithm 4, we use two predicates in the exponential march and the binary search to find a closest point to T in P: the predicate "is x in P" because we are only interested in points in P, and the in-sphere test: "does the sphere passing by T ∪ { x} contain y?" in order to find a closest point. This in-sphere test is made with a determinant computation.

Correctness of plane exploration algorithms

We start by giving some properties that are valid for both algorithms. We denote by n the last iteration. The proof that the algorithm always terminates is postponed to section 3.2, theorem 1.

Algorithm invariants and characterization of update operations

Lemma 1 The following properties are true:

1. ∀i ∈ {0, . . . , n}, ∀k, v

(i) k ∈ P.
Algorithm 4: ClosestPointInRay(T,R σ , q) Input: a 3-point sequence T, a ray R σ , a shifted point q x ← qm σ(0) + m σ(1) ; // starting point of the ray R σ (must be in P) y ← m σ(2) ; // direction vector of the ray R σ // exponential march κ ← 0, λ ← 1 ; // lower and upper bound respectively while x + κ y ∈ P and the sphere passing by T ∪ { x + κ y} contains

x + λ y do κ ← λ, λ ← 2λ ; κ ← κ 2 ;
// binary search while (λκ) > 4 do 3. ∀i ∈ {0, . . . , n}, q is strictly above aff(T (i) ) in direction s.

// invariant: x + κ y ∈ P α = (3κ+λ)
4. ∀i ∈ {0, . . . , n -1}, for a point x ∈ (N (i) ∩ P), the convex hull of T (i) ∪ { x } has non-zero volume.

5. ∀i ∈ {0, . . . , n -1}, T (i) and T (i+1) have two vertices in common.

Proof. We prove the first item by induction. Others follow. The first property is obviously true for i = 0. Let us assume now that it is true for i < n. The set (N (i) ∩ P) contains at least one point because i < n (when the set is empty, the algorithm stops and i = n). For a point x ∈ (N (i) ∩ P), let us consider the tetrahedron defined as the convex hull of T (i) ∪ { x }. Since T (i) ⊂ P by the induction hypothesis, any three-point subsets of T (i) ∪ { x } (and therefore T (i+1) ) belong to P, which proves item 1.

We can now prove items 2 and 3. Indeed, let us assume that { m

(i) k } k ∈{0,1,2} are
coplanar. Then, by construction (algorithm 1, line 4), q belongs to conv(T (i) ). However, since q / ∈ P and T (i) ⊂ P, this contradicts the convexity of P, which proves item 2 by contradiction.

The same argument of convexity may be used to prove that q is strictly above aff(T (i) ) in direction s for all i ∈ {0, . . . , n} (item 3).

To end, we prove items 4 and 5. Since q is strictly above aff(T (i) ) and since the neighborhoods are defined in ( 2) and ( 4) by vectors { m (i) k } k ∈{0,1,2} going from points of T (i) to q (by ( 1)), all the points of N (i) are strictly above aff(T (i) ) in direction s. As a consequence, the convex hull of T (i) ∪ { x } has non-zero volume (item 4) and T (i) = T (i+1) . Since the convex hull of T (i) ∪ { x } is a non degenerate tetrahedron, there are three triangular facets different from conv(T (i) ). All of them (including conv(T (i+1) )) share exactly two vertices with T (i) and have point x as third vertex (item 5).

The following lemma fully characterizes the main operation of algorithm 1 (lines 3-5):

Lemma 2 ∀i ∈ {0, . . . , n -1}, ∃ k s.t.          v (i+1) k = v (i) k + α m (i) k +1 + β m (i) k +2 , with α, β ∈ N, α = 1 or β = 1, α + β ≥ 1, v (i+1) k +1 = v (i) k +1 , v (i+1) k +2 = v (i) k +2
Proof. Let us assume w.l.o.g. that k = 0. By algorithm 1, ( 2) and ( 4), for a permutation σ over {0, 1, 2}, ∀i ∈ {0, . . . , n -1}, v 2) , and by (1),m 2) . We will prove below by contradiction that σ(0) = 0. Let us assume that σ(0) = 1 (the case where σ(0) = 2 is similar).

(i+1) 0 = q -m (i) σ(0) + α m (i) σ(1) + β m (i) σ(
(i+1) 0 = -m (i) σ(0) + α m (i) σ(1) + β m (i) σ(
Let us consider det(-m

(i) 1 , -m (i) 2 , -m (i+1) 0 ). Replacing -m (i+1) 0 by -m (i) 1 + α m (i) 0 + β m (i) 2
in the previous determinant, we obtain the following identity:

det(-m (i) 1 , -m (i) 2 , -m (i+1) 0 ) = α det(-m (i) 1 , -m (i) 2 , m (i) 0 ). If α = 0, then det(-m (i) 1 , -m (i) 2 , -m (i+1) 0 ) = 0, which is in contradiction with lemma 1-item 2.
Otherwise, since α > 0, det(-

m (i) 1 , -m (i) 2 , -m (i+1) 0
) and det(-

m (i) 1 , -m (i) 2 , -m (i)
0 ) must have opposite signs. By [START_REF] Berthé | Brun expansions of stepped surfaces[END_REF], it follows that the plane passing by v

(i) 1 , q, v (i) 2 separates v (i) 0 from v (i+1) 0
(see fig. 4.a). In this case, the straight line passing by p and q can intersect both T (i) and

T (i+1) only if it intersects segment [ v (i) 1 v (i) 2 ].
However, a ray that goes inside the convex hull of

T (i) ∪ T (i+1) through edge [ v (i) 1 v (i)
2 ] must exit the convex tetrahedron through a single point of a facet that is not T (i+1) (see fig. 4.b), which raises a contradiction (T (i+1) is not the upper facet intersected by such a ray in this case).

We conclude that σ(0) = 1 and similarly that σ(0) = 2, which implies that σ(0) = 0.

This characterization leads to extra properties: lemma 3 and corollary 2 are stronger versions of lemma 1, items 2 and 3, respectively.

v i 1 v i 2 v i 0 q v i+1 0 (a) v i 1 v i 2 v i 0 v i+1 0 (b)
(i) 1 , q, v (i) 2 sep- arates v (i) 0 from v (i+1) 0
. In (b), a ray that goes inside the tetrahedron through edge [ v

(i) 1 v (i) 2 ] must exit through a point of a facet that is not T (i+1) .
Let M (i) be the 3 × 3 matrix formed by the three vectors joining the vertices of the current triangle to q. Otherwise said, it consists of the three row vectors ( m

(i) k ) k ∈{0,1,2}
. We prove below that M (i) is unimodular, which is an important property to show that the algorithm returns a basis of upper leaning points at termination (see theorem 2 and corollary 5).

Lemma 3 ∀i ∈ {0, . . . , n}, det (M (i) ) = 1.
Proof. It is easily checked that det (M (0) ) = 1. We now prove that if det (M (i) ) = 1 for ∀i ∈ {0, . . . , n -1}, then det (M (i+1) ) = 1. As proven in lemma 1-item 5, only one vertex changes at each step, for instance v (i) k . If we look at lemma 2, we have v

(i+1) k = v (i) k + α m (i) k+1 + β m (i) k+2 , which is equivalent to m (i+1) k = m (i) k -α m (i) k+1 -β m (i)
k+2 , for some non negative integers α, β such that α or β equals to 1, α + β ≥ 1. The other vertices are not modified so the remaining two rows of M (i+1) are not modified. We get

det(M (i+1) ) = det( m (i) k -α m (i) k+1 -β m (i) k+2 , m (i) k+1 , m (i) k+2 ) = det( m (i) k , m (i) k+1 , m (i) k+2 ) (by linearity) = det(M (i) ) = 1.
(by induction hypothesis)

The height of each vector m k in the direction given by the estimated normal is equal to 1.

Corollary 1 ∀i = 0, . . . , n, M (i) • N(T (i) ) = 1. Proof. M (i) • N(T (i) ) = 1 because ∀k, ( d (i) 0 × d (i) 1 ) • m (i) k = (( m (i) 0 -m (i) 1 ) × ( m (i) 1 -m (i) 2 )) • m (i) k = ( m (i) k+1 × m (i) k+2 ) • m (i) k = det (M (i)
), which is equal to 1 by lemma 3.

The height of q in the direction given by the estimated normal is equal to 1.

Corollary 2 ∀i = 0, . . . , n, q • N(T (i) ) = 1.

Proof. Since ∀i = 0, . . . , n, q = v

(i) 0 + m (i) 0 , we compute on one hand v (i) 0 • N(T (i)
), which is equal to 0 by definition, and on the other hand m (i) 0 • N(T (i) ), which is equal to 1 by corollary 1.

To end, the following lemma leads to a strong geometrical property : the straight line passing by p and q intersects the interior of every triangle.

Lemma 4 ∀i ∈ {0, . . . , n}, ∀k, ( m

(i) k × m (i) k+1 ) • s > 0. Proof. By definition, ∀k, ( m (0) k × m (0) k+1 ) • s = 1.
Let us now assume that ∀k, ( m

(i) k × m (i) k+1 ) • s > 0 and let us prove that ∀k, ( m (i+1) k × m (i+1) k+1 ) • s > 0.
Let k be the index of the vertex of T (i+1) , that is not a vertex of T (i) . By lemma 2, we have m

(i+1) k = m (i) k -α m (i) k +1 -β m (i)
k +2 , for some non negative integers α, β such that α or β equals to 1, α + β ≥ 1, while m

(i+1) l = m (i) l for l ∈ {0, 1, 2} \ k .
We must check two expressions involving m (i+1) k (the remaining one does not change).

The first one is:

( m (i+1) k × m (i+1) k +1 ) • s = ( m (i) k × m (i) k +1 ) • s -α( m (i) k +1 × m (i) k +1 ) • s -β( m (i) k +2 × m (i) k +1 ) • s.
We conclude that ( m

(i+1) k × m (i+1) k +1 ) • s > 0 because • ( m (i) k × m (i) k +1 ) • s > 0 by induction hypothesis • -α( m (i) k +1 × m (i) k +1 ) • s = 0 • -β( m (i) k +2 × m (i) k +1 )• s = β( m (i) k +1 × m (i) k +2 )• s ≥ 0 by induction hypothesis (and β ≥ 0).
The second expression is similar.

By (1), lemma 4 is equivalent to:

∀i ∈ {0, . . . , n}, ∀k, ( q -m (i) k ) × ( q -m (i) k+1 ) • ( q -p) > 0,
which means that p is strictly in the first octant of the frame ( q;m

(i) 0 , -m (i) 1 , -m (i)
2 ). This implies that the straight line passing by p and q intersects the interior of the triangle whose vertices are ( qm (i)

k ) k ∈{0,1,2} , i.e. ( v (i) k ) k ∈{0,1,2} .
This result guarantees that there is no ambiguity in the computation of T (i+1) in line 4 of algorithm 1, since the straight line passing by p and q never crosses an edge of conv(T (i) ∪ T (i+1) ), but only the interior of T (i) and T (i+1) .

Termination

In the following proofs, we compare the position of the points along direction N. For the sake of simplicity, we use the bar notation • above any vector x to denote its height relative to N. Otherwise said, x := x • N. Even if N is not known, q ≥ ω by definition and for all x ∈ P, 0 ≤ x < ω. By (1) and lemma 1-item 1 we straightforwardly get the following lemma:

Lemma 5 ∀i ∈ {0, . . . , n -1}, ∀k, m (i) k > 0.
As a consequence, any operation strictly increases the height of the updated vertex.

Lemma 6 ∀i ∈ {0, . . . , n -1}, let k be the index of the updated vertex such that v

(i+1) k = v (i) k . Then, v (i+1) k > v (i) k . A corollary is m (i+1) k < m (i) k .
Proof. By lemma 2, we have ∀i ∈ {0, . . . , n -1}, v

(i+1) k = v (i) k + α m (i) k +1 + β m (i) k +2 , with two non negative integers α, β such that α = 1 or β = 1, α + β ≥ 1. Since ∀k, m (i)
k > 0 by lemma 5 and α, β ≥ 0, but α and β are not both equal to 0, we clearly have v

(i+1) k > v (i) k and, by (1), m (i) k < m (i+1) k .

The termination theorem follows:

Theorem 1 The number of steps in algorithm 1 is bounded from above by ω-3.

Proof. The result comes from the fact that the sequence ( k m (i) k ) i=0,...,n is a strictly decreasing sequence of integers between ω and 3 because:

• ∀k, m (0) k = e k+2 and k m (0) k = ω.
• by lemma 5, ∀i ∈ {0, . . . , n}, ∀k, m

(i) k > 0 and k m (i) k ≥ 3. • by lemma 6, ∀i ∈ {0, . . . , n -1}, k m (i) k > k m (i+1) k .
Remark that this bound is tight since it is reached when running the algorithm on a plane with normal N(1, 1, r).

Correctness

We show that our two plane-probing algorithms extract the correct normal of the input digital plane. Let us begin with a small technical lemma: Lemma 7 For any permutation σ over {0, 1, 2}, ∀i ∈ {0, . . . , n}, if there is a point x of ray R (i) σ that is not in P, then no point further than x on the ray is in P.

Proof.

For two non negative integers λ and λ such that λ < λ , let

x := q -m (i) σ(0) + m (i) σ(1) + λ m (i) σ(2) and y := q -m (i) σ(0) + m (i) σ(1) + λ m (i) σ(2) be two points of R (i) σ (see equation 3). If x /
∈ P, then either x < 0 or x ≥ ω. However, it cannot be the former because x = v (i)

σ(0) + m (i) σ(1) + λ m (i) σ(2) and • v (i) σ(0) > 0 by lemma 1-item 1, • m (i) σ(1) , m (i) σ(2) > 0 by lemma 5,
• λ is assumed to be non negative.

We can now bound from below the height of y 2) > 0 (by lemma 5) and (λλ) > 0, we have y ≥ ω, which implies that y / ∈ P.

= x + (λ -λ) m (i) σ(2) . Since x ≥ ω, m (i) σ(
Corollary 3 ∀i ∈ {0, . . . , n}, (N (i)

H ∩ P) = ∅ ⇒ (N (i) R ∩ P) = ∅.
Proof. Due to the neighborhood definitions (2) and ( 4), ∀i ∈ {0, . . . , n}, any point y in

N (i) R but not in N (i)
H is located in a ray whose starting point x is in

N (i) H and y = x + λ m (i)
k for some index k ∈ {0, 1, 2} and non negative integer λ.

If (N (i)

H ∩ P) = ∅, then x / ∈ P and the result follows by lemma 7. This corollary implies in particular that at the last step we can focus on the H-neighborhood. Since we focus below on the last step n, we omit the exponent (n) in the proofs to improve their readability.

We now give the correctness result when the starting point p is a lower leaning point, i.e. p = 0.

Theorem 2 If p is a lower leaning point (i.e. p = 0 and thus q = ω), the vertices of the last triangle are upper leaning points, i.e. ∀k, v

(n) k = ω -1.
Proof. The first step of the proof is to show that the vertices of the last triangle are all at the same height, i.e. m 0 = m 1 = m 2 . If not, then there exists k ∈ {0, 1, 2} such that d k = 0. In this case, either (i) d k < 0 or (ii) d k > 0.

Since q = ω and | d k | < ω, either (i) q + d k ∈ P or (ii) qd k ∈ P. This implies that N ∩ P = ∅, which is a contradiction because N ∩ P = ∅ at the last step (see algorithm 1, l. 2). As a consequence, ∀k, d k = 0 and ∀k, m k = γ, a strictly positive integer.

The second step of the proof is to show that γ = 1.

Let us denote by 1 the vector (1, 1, 1) T . We can write the last system as MN = γ1. Since M is invertible (because det (M) = 1 by lemma 3), N = M -1 γ1 and as a consequence γ = 1 (because the components of N are relatively prime and M -1 is unimodular).

We conclude that ∀k, m k = 1 and, straightforwardly,

v k = ω -1.
The following two corollaries are derived from lemma 3 and theorem 2.

Corollary 4 If p is a lower leaning point, the normal of the last triangle is equal to N, i.e. N(T (n) ) = N.

Proof. On one hand, M N(T) = 1 because ∀k, (

d 0 × d 1 ) • m k = (( m 0 - m 1 ) × ( m 1 -m 2 )) • m k = ( m k+1 × m k+2 ) • m k = det (M), which is equal to 1 by lemma 3.
On the other hand, MN = 1 by theorem 2. Since M is invertible, we have

N(T) = N. Corollary 5 If p is a lower leaning point, ( d (n) 0 , d (n) 1 ) is a basis of the lattice of upper leaning points { x ∈ P | x = ω -1}.
Proof. The unit parallelepiped in the lattice

{( q + α m 0 , q + β m 1 , q + γ m 2 )|(α, β, γ) ∈ Z 3 }
does not contain any integer point because it is equivalent to Z 3 (det (M) = 1 by lemma 3). It follows that the facet conv(T) does not contain any integer point. Since the points of T are upper leaning points by theorem 2, ( d 0 , d 1 ) is a basis of the lattice of upper leaning points.

We end the section by providing the worst-case time complexity of both algorithms in a computation model where the evaluation of the predicate "is x in P" only requires a constant time: Theorem 3 If p is a lower leaning point, the H-algorithm (resp. R-algorithm) returns three upper leaning points of P in O(ω) (resp. O(ω log ω)), where ω is the arithmetical thickness of the digital plane.

Proof. The time complexity and correctness of the H-algorithm straightforwardly comes from theorem 1 and theorem 2 respectively, because each iteration runs in O(1) (N (i) H ∩ P contains at most 6 points and algorithm 2 runs in linear time, see sec. 2.4).

However, the time complexity and correctness of the R-algorithm depend also on the time complexity and correctness of algorithm 4, run at each iteration to find the closest point in P to the current triangle. Algorithm 4 uses an exponential march (lines 4-5) followed by a binary search (lines 7-14). The relevance of such an approach comes from lemma 5 and lemma 7 (along a ray, points in P are followed by points lying above P) and corollary 2 (which implies, together with ( 1) and ( 3), that the function that maps a point x in the ray sequence to the radius of the sphere passing by T ∪ { x} is convex and has a global minimum).

Let us focus now on the complexity of algorithm 4. After line 3, κ = 0 and x + κ y ∈ P (precondition, see algorithm 3, line 3). Therefore, x + κ y ≤ ω. Moreover, we know by lemma 5 and lemma 7, that there must be a greater value for κ, such that x + κ y > ω. It is clear that this value is reached after at most log ω + 2 iterations in the exponential march (lines 4-5). After line 6, the size of the range [κ, λ] is at most 3ω 2 . In the binary search (lines 7-14), there are O(log ω) iterations because each iteration shrinks the range to half its size until a size less than 4. The last line takes a constant time since algorithm 2 runs in linear time and the cardinality of the input point set is at most 4. The overall complexity of algorithm 4 is thus O(log ω).

Since there are O(ω) iterations in algorithm 1 (theorem 1) and since algorithm 4 is used at each iteration at most 6 times in algorithm 3, we conclude that the overall complexity of the R-algorithm is O(ω log ω).

It is worth noting that in both cases, the time complexity corresponds to the number of calls to the predicate "is x in P ?". This means that the time taken by a call to the predicate impacts directly the multiplicative constant in O(ω) (resp. O(ω log ω)).

Experimental evaluation

In this section, we conduct an experimental evaluation of both H-and Ralgorithms. Furthermore, we also compare these new algorithms to the planeprobing algorithm presented in [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF] called FindNormal. We evaluate the number of steps, the number of calls to the predicate "is x in P ?" as a function of the norm of the normal vector of P. We also check the ability of the algorithms to produce reduced lattice basis, and, in the case where the basis is not reduced, the number of lattice reduction operations necessary to transform it into a reduced basis.

The graphics on fig. 5 and fig. 6 shows that the three algorithms have quite a similar behavior. For the three algorithms, the number of steps is clearly sublinear on average. That being said, there are still some cases that reach the linear bound of theorem 1. Unlike the FindNormal algorithm, at each step the H-and R-algorithms select a point based on geometrical criteria. These criteria are stronger in the case of the R-algorithm which explains that in general, it terminates with less steps than the others.

Regarding the number of points tested or, equivalently, the number of calls to the predicate "is x in P ?", the H-algorithm shows a better behavior on average. Of course, the systematic exploration of 6 rays using algorithm 2 generates extra calls to the predicates.

Corollary 5 states that the edge vectors of the last triangle form a basis of the lattice of upper leaning points to P. FindNormal from [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF], middle: H-algorithm and bottom: R-algorithm. All graphics are made using the same vectors picked randomly in such way that their 1-norm is located in the interval displayed by the width of each column.

Figure 6: Number of calls to the predicates "is x ∈ P" versus the 1-norm of the normal vector N. Top: the algorithm is FindNormal from [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF], middle: H-algorithm and bottom: R-algorithm. All graphics are made using the same vectors picked randomly in such a way that their 1-norm is located in the interval displayed by the width of each column.

Let us recall that a basis ( x, y) is reduced if and only if x , y ≤ xy ≤ x + y . Given ( x, y) a basis of a two dimensional lattice, there is an iterative algorithm to compute a reduced basis from ( x, y). This algorithm consists in replacing the longest vector among ( x, y) by the shortest one among x + y and xy, if it is smaller. This operation is called a reduction.

In [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF], a solution is proposed in order to generate reduced basis. This method can be summarized as if at one step, a non-reduced basis is formed and that it may be corrected using a reduction, then do it. In H-and R-algorithm, we return the two shortest vectors in { d

(n) k } k ∈{0,1,2}
as a basis. We do not perform any reduction because such a basis is almost always reduced. We ran all three algorithms (FindNormal, H-algorithm, R-algorithm) on all vectors ranging from (1,1,1) to (200,200,200). There are 6578833 vectors with relatively prime components in this range. Tab. 1 shows that less than 0.01% of basis computed by the H-algorithm were non-reduced. Regarding the R-algorithm, not only all bases computed in the range (1,1,1) to (200,200,200) were reduced but we have performed billions of tests on different normal vectors and we have never found a basis that was not reduced.

Digital surface analysis

In this section, we consider a set of voxels, Z, where voxels are seen as axisaligned unit cubes whose vertices belong to Z 3 . The digital boundary, BdZ, is defined as the topological boundary of the union of the voxels of Z. Since a digital boundary looks locally like a digital plane, it is natural to run our plane probe algorithms at each reentrant corner of the digital boundary with predicate "is x in BdZ ?" in order to estimate a local tangential facet to the volume Z (see fig. 7). This facet also defines naturally a local normal vector to the volume Z.

Since the predicate "is x in BdZ ?" defines only locally a plane, our planeprobing algorithms must be slightly adapted to this new context. We discuss first the case where Z is the digitization of a convex set, and we address the problem of initializing the algorithms at bad reentrant corners, that is a corner that does not correspond to a lower leaning point. Then we present how Figure 7: Our H-algorithm has been run at each reentrant corner of a digital plane (left) and an ellipse (right). The last triangle of each run is printed and each triangle is associated one-one to a pattern. On the left, the large triangles that regularly tile the digital plane are Bezout patterns; the small triangles hidden underneath are reduced patterns.

we can use the H-neighborhood to detect planarity defects. Last we explain how to make facets follow closely the local planar geometry along the digital surface. Presented experiments use H-algorithm, but would be identical with R-algorithm.

From now on, we call pattern any tuple ( q, s, m 0 , m 1 , m 2 ), such that q and s define a reentrant corner of BdZ and qm 0 , qm 1 , qm 2 is the output facet of the H-algorithm run at this reentrant corner. In proofs, since vector s will be (1, 1, 1) for all considered patterns, we will omit vector s in the tuples.

Pattern on convex shapes

Let us assume that Z is a digitally convex shape, i.e. the digitization of its convex hull is Z itself. Let F be a facet of the convex hull Conv(BdZ) of BdZ, oriented so that its normal points outside. The shift vector s F of F is the vector (±1, ±1, ±1) whose component signs match the sign of the normal vector to F. We define the set of boundary points of BdZ below F as

BdZ(F) := x ∈ BdZ ∩ Z 3 s.t. [ x, x + s F [ ∩ F = ∅ .
By convexity, for any point y ∈ BdZ(F), the point y + s F does not belong to BdZ.

It is clear that BdZ(F) is a piece of some digital plane P: it suffices to define P as the digital plane with normal vector identical to the normal vector of F and with intercept such that the vertices of F are all upper leaning points.

If a pattern ( q, s, m 0 , m 1 , m 2 ) rooted at point qs ∈ BdZ(F) and with s = s F has its triangle aligned with F, then qs is a lower leaning point of the digital plane carrying BdZ(F). In this case, q is called a Bezout point and ( q, s, m 0 , m 1 , m 2 ) is called a Bezout pattern of F. Unfortunately not every facet of Conv(BdZ) has a Bezout pattern. This is illustrated on fig. 7. We can see that approximately half of the facets of the convex hull of BdZ are extracted.

The other half of the facets are not extracted because these facets have no Bezout point above them ! Their Bezout point is above the other half of the fundamental domain.

An interesting consequence is that we can use the first plane-probing algorithm of [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF] to extract the complementary facets. This is illustrated on fig. 8. Almost all the geometry of the convex shape is captured. Missing parts are related to facets with normal vectors with a null component.

Reduced patterns

If the pattern with corner point ps ∈ BdZ(F) is not aligned with the facet F then it is called a reduced pattern of F. In this case the pattern is only approximately aligned with the plane. You can see some examples of reduced patterns on the digital plane of fig. 7, left. There are indeed several reentrant corners on such planes that are not located beneath a Bezout point. Starting from such corners make H-algorithm stops prematurely and outputs a smaller triangle that is only approximately aligned with the local tangent plane.

We propose to detect a posteriori the reduced patterns with algorithm 5. Its principle is simple. Let us define Plane( q, m 0 , m 1 , m 2 ) as the digital plane of upper leaning points m 0 , m 1 , m 2 and Bezout point q. Since there are Bezout patterns among all the input patterns, their vectors m k climb along the normal to the facet one layer at a time (their height is 1). Reduced patterns correspond to starting shifted point q that are not Bezout point. Their height (along the normal to some facet) is thus greater than the height of Bezout points, i.e. ω. So, in the first iterations, the algorithm finds the reduced patterns whose shifted point q has height ω + 1. Those shifted points are put back in the queue, but with the Bezout vectors that climb one by one. So, in the following iterations, the reduced patterns whose shifted point q has height ω + 2 are reached, and so on. This is illustrated on fig. 9.

// Push back q but with updated vectors. If it is already inside, update vectors only. Q.push( ( q , s , m 0 , m 1 , m 2 ) ); A Bezout pattern of this digital plane is displayed in dark green, its shifted point q is in green. The plane-probing algorithm was also run at the reentrant corner of height 1, its output is the reduced pattern displayed in dark red, with shifted point q in red. It is easily seen that q = q + m 2 and is thus located by algorithm 5.

We also remark that reduced patterns can be reached sooner (e.g. a shifted point of height ω + 3 can be reached from a shifted point of height ω + 1 with a vector m k climbing two by two). This is fine, we just want to mark reduced patterns. On a digital surface, it may also happen that one pattern can be reached from another pattern, but they do not correspond to the same digital plane. This is also tested in the algorithm: in order to be marked, the reduced pattern must be included in the digital plane carried by the Bezout pattern (line 10).

We now prove that algorithm 5 terminates on arbitrary BdZ. To make the explanation simpler, we suppose that the shift vector is equal to 1 := (1, 1, 1) for each pattern (see also line 8 of the algorithm, where equality of shift vectors is tested). The two following lemmas show how to build an order on patterns.

Lemma 8 Given a pattern P = ( q, m 0 , m 1 , m 2 ), denoting by M the matrix composed of rows of m k , then the digital plane Plane(P) has the following parameters (resp. normal N, intercept µ and thickness ω):

N = M -1 1, µ = q • N -1 • N, ω = 1 • N.
To lighten the exposition, we will speak of the normal, intercept and thickness of the pattern P, to refer to these characteristics defined on the digital plane Plane(P). We say that a pattern P = ( q , m 0 , m 1 , m 2 ) is included in the pattern P = ( q, m 0 , m 1 , m 2 ), and we write P ⊂ P, whenever every q -1, qm 0 , qm 1 , qm 2 belongs to Plane( q, m 0 , m 1 , m 2 ). We have: Lemma 9 If P ⊂ P then it holds that N ≤ N, where ≤ means less or equal for all components. A corollary is that ω ≤ ω.

Proof. Since we have ( q -1) ∈ Plane(P), we have ( q -1) • N ≥ µ and it follows from lemma 8 that q • N ≥ q • N. Since we have ( qm k ) ∈ Plane(P), we have ( qm k )

• N ≤ µ + ω -1. Using lemma 8, this is equivalent to m k • N ≥ 1 + q • N -q • N.
The first relation entails that M N ≥ 1 when writing it in matrix form. Since it holds that M N = 1 and M is invertible, the result follows.

Proposition 1 Algorithm 5 terminates in less than nω iterations, where n is the number of patterns in L and ω is the maximal thickness of the patterns in L.

Proof.

First this algorithm manages patterns in the queue. Then it is clear that the number of patterns cannot increase, since patterns inserted back into the queue are associated to existing shifted points of L. Since only patterns with the same shift vectors are compared, we limit our reasoning to patterns with the same shift vector, here the positive orthant s = (1, 1, 1).

The key argument is that whenever some pattern is pushed back into the queue, its thickness strictly grows as well as its intercept. At some point, there won't be any pattern that can include another one, so the algorithm stops.

More precisely, looking at lines 8 and 10, we denote by P = ( q, m 0 , m 1 , m 2 ) the pattern in the list L and by P = ( q , m 0 , m 1 , m 2 ) the pattern popped from the queue. The pattern that is pushed back into the queue is denoted by P . Using the above notations for thickness and intercept of digital planes associated to patterns (Lemma 8), we obtain straightforwardly that:

N = N, ω = ω, µ = µ + 1.
Furthermore, since a success of the test at line 10 means P ⊂ P, we have ω < ω = ω (the "≤" comes from Lemma 9, the strictness "<" comes from q = q + m k ). We may also state that P ⊂ P and P ⊂ P since their intercepts are different while their thickness are equal.

To sum up, the queue contains patterns with strictly increasing thicknesses. However, since the new inserted pattern P has the same thickness as some already existing pattern, it is clear also that the thickness of all patterns is upperbounded by the thickness of the thickest pattern of L, i.e. ω. If we follow one pattern of L, it may thus be reassigned in Q at most ω -1 times before having maximal thickness. At this moment, it cannot be included in any pattern of L, since the tests at lines 8 and 10 may be true only if the thickness of P is strictly lower than the thickness of a pattern of L. So this pattern is definitively popped out of the queue. It follows than the algorithm terminates in at most nω iterations.

In practice, as shown on tab. 2, this algorithm runs extremely fast and much fewer iterations than nω are needed. As expected, the algorithm is slower when analyzing shapes made of wide planar sides with important arithmetic thickness (e.g. fandisk 512 3 ).

Detecting planarity defects

If the shape is not convex, the algorithm can be adapted to detect planarity defects. The idea is that, at a step i, N

H ∩ BdZ must contain at most three elements if it is locally a piece of plane. Moreover these elements must always be consecutive neighbors around q in N (i) H . We thus stop the algorithm whenever at least one of the two previous conditions fails, meaning that the surface is locally non convex. We sum up the possible cases in the following table, where the predicate values of the points of the H-neighborhood are symbolized with white circles for ∈ BdZ and black circles for ∈ BdZ. When our plane probe algorithm stops, it returns in addition to the output triangle the stopping criterion among "convex or planar" and "non convex". These various stopping criteria are illustrated in fig. 10. Convex and planar zones contain green triangles, which are patterns which stops on a "convex and planar" H-neighborhood. Inflexion zones, saddle points, or places with at least one negative principal curvature contain mostly magenta triangles, meaning "non-convex" H-neighborhood, or yellow triangles, which refer to non-separating patterns.

H-neighborhood configurations Stop

Local planarity yes convex or planar no (still probing) yes non-convex

Following closely digital surface geometry

Our algorithm probes for points in a plane in a sparse way. This works well when identifying a true digital plane, but it may badly identify pieces of plane on digital surfaces. For instance the algorithm may jump over holes or cracks in the surface. Therefore, we have to check at each iteration that the current triangle fits closely the digital surface. We proceed as follows:

• the current triangle is denoted by T = ( v 0 , v 1 , v 2 ), the shifted point q and the shift vector s;

• the corresponding digital plane is D := Plane( q, qv 0 , qv 1 , qv 2 );

• we compute the digital triangle under T on D as

A := { x ∈ D, [ x, x + s[ ∩ conv(T) = ∅},
by a breadth-first traversal from q on D;

• we compute also the three digital segments B k , k ∈ {0, 1, 2}, along the border of T as the 6-connected path in D from v k to v k+1 whose points are closest to the straight line segment [ v k , v k+1 ];

• then the current triangle is said separating iff A ∪ B 0 ∪ B 1 ∪ B 2 ⊂ BdZ.

When running our plane-probing algorithm, we check at each step if the current triangle is separating. If not, the algorithm exits with stopping criterion "non-separating". An illustration of such cases is given in fig. 10.

The sets B k are important to connect points of A. Should we not consider them, the current triangle might have a needle shape that is separating while being far away from the surface (A is reduced to the three vertices of T).

Some timings and statistics of our algorithms to analyze digital surface geometry are given in tab. 2.

Conclusion and perspectives

In this paper, we proposed two new algorithms that compute the parameters of a digital plane. In opposition to usual plane recognition algorithms, these Figure 10: Illustration of H-algorithm run on every reentrant corner of digital surfaces "cube and sphere 128 3 " (top row) and "fandisk 128 3 " (bottom row). Extracted triangles are colored according to the stopping criterion: "convex or planar" is green, "non-convex" is magenta, "non-separating" is yellow. All patterns are displayed on left column. On right column, reduced patterns have been removed with algorithm 5. algorithm greedily decide on-the-fly which points have to be considered like in [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF]. We have called plane-probing algorithms that kind of plane recognition algorithms. Compared to [START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF], these algorithms are however simpler because they consist in iterating one geometrical operation. Furthermore the returned solution, which is described by a triangle parallel to the sought plane, always lies above the starting corner. Besides, the two shortest edge vectors of the triangle almost always form a reduced basis for the purely local H-algorithm. For the sometimes less local R-algorithm, it has always returned a reduced basis.

Shape

For the future, we would like to prove that the R-algorithm always returns a reduced basis. Moreover, we would like to find a variant of our algorithm in order to retrieve complement triangles whose Bezout point is not above the triangle. For sake of completeness, we are also interested in degenerate cases, where at least one component is null. We would also like to recognize a piece of digital planes, i.e. to extract the digital plane with minimal characteristics containing this piece, but with theoretical guarantees. The difficulty is that our plane-probing algorithms need few but specific points of the plane to fully recognize it. We feel that we still need a correct definition of what is a valid piece of digital plane in 3D. Although connectedness was a sufficient condition for digital segment in 2D, it is not enough in the 3D case and we are currently working on it.

After having achieved these goals, we would have a complete working tool for the analysis of digitally convex surfaces. General digital surfaces are even more complex to analyze. We have shown that our plane-probing algorithms are able to detect non-convex parts. To go further, it appears necessary to precisely associate to a given output triangle a subset of points of the digital surface. Then we will be able to analyze how these subsets of points overlap each other in order to delineate convex, concave and saddle zones. The segmentation of a digital surface into such zones is certainly an essential step in digital shape analysis, which would greatly facilitate higher-level processing, like digital shape matching, indexing or recognition.

Figure 1 :

 1 Figure 1: Illustration of the running of[START_REF] Lachaud | Computation of the normal vector to a digital plane by sampling significant points[END_REF] Algorithm 2] and the algorithm from[START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF] on a digital plane of normal vector[START_REF] Berthé | Brun expansions of stepped surfaces[END_REF][START_REF] Brimkov | Digital planarity -a review[END_REF][START_REF] Chica | Pressing: Smooth isosurfaces with flats from binary grids[END_REF]. Images (a) to (d) show the four iterations of [16, Algorithm 1] starting from the origin. Each triangle is intersected by the small dashed red segment. In (d), the normal of the last triangle is[START_REF] Berthé | Brun expansions of stepped surfaces[END_REF][START_REF] Brimkov | Digital planarity -a review[END_REF][START_REF] Chica | Pressing: Smooth isosurfaces with flats from binary grids[END_REF]. Images (e) to (h) show iterations 0 (initial), 2, 5 and 7 (final) of the algorithm from[START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF]. The initial tetrahedron (a) is placed at the origin and the final one (h) has an upper triangle with normal vector (1, 2, 5).

Figure 2 :

 2 Figure 2: The triangle conv(T (i) ) is depicted in grey. The set N

Algorithm 2 : 1 foreach x ∈ S do 2 if 3 y ← x ; 4 return y ; 5

 212345 ClosestPointInSet(T,S)Input: a 3-point sequence T, a point set S Let y be an arbitrary point of S ; the sphere passing by T ∪ { y} strictly contains x then

4 , β = (κ+λ) 2 , γ = (κ+3λ) 4 ; 13 κ 14 return 1 v 2 Figure 3 :

 4241314123 Figure 3: Two-dimensional illustration of the function that maps a point x on a ray (in green) to the radius of the sphere passing by T ∪ { x} (blue circles).

Figure 4 :

 4 Figure 4: Illustration of lemma 2. In (a), the plane passing by v

Figure 5 :

 5 Figure 5: Number of iterative steps performed by the plane-probing algorithms as a function of the 1-norm of the normal vector N. Top: the algorithm isFindNormal from[START_REF] Lachaud | An output-sensitive algorithm to compute the normal vector of a digital plane[END_REF], middle: H-algorithm and bottom: R-algorithm. All graphics are made using the same vectors picked randomly in such way that their 1-norm is located in the interval displayed by the width of each column.

Figure 8 :

 8 Figure 8: Left: a digital surface. Middle: triangular faces (in green) produced by the H-algorithm starting from each reentrant corners. Right: triangular faces (in blue) produced by the FindNormal algorithm are showed together with the previous ones.

Figure 9 :

 9 Figure9: Locating reduced patterns with algorithm 5. A Bezout pattern of this digital plane is displayed in dark green, its shifted point q is in green. The plane-probing algorithm was also run at the reentrant corner of height 1, its output is the reduced pattern displayed in dark red, with shifted point q in red. It is easily seen that q = q + m 2 and is thus located by algorithm 5.

  •

  •

Table 1

 1

	Algorithm	FindNormal	H-alg.	R-alg.
	nb. non-reduced	6197115 ( 94.2%) 480 ( 0.007%)	0
	avg. nb. reductions	6	1	0
	max. nb. reductions	117	1	0
	: Algorithms FindNormal, H-algorithm and R-algorithm were used on
	all 6578833 vectors ranging from (1, 1, 1) to (200, 200, 200) with relatively prime
	coordinates. The average number of reduction is computed only among non-
	reduced basis.

Table 2 :

 2 Total timings of H-algorithm (algorithm 1) runned at each point of BdZ for several shapes. Note that this time includes the time for checking the separation of triangles. The time taken by algorithm 5 to remove reduced patterns is listed. The total number of extracted patterns and the number of patterns marked as reduced are listed.

	Z	#BdZ	Time (ms) H-algorithm	Time (ms) reduced	#patterns #reduced
	cube and sphere 128 3	34034	1513	30	7041	3266
	fandisk 128 3	48220	2732	69	12106	6984
	fandisk 256 3 186760	13954	791	45999	31948
	fandisk 512 3 734658	71660	18201	178751	142293
	octaflower 512 3	692738	47299	944	179905	98816

* This work has been partly funded by DigitalSnow ANR-11-BS02-009 research grant and CoMeDiC ANR-15-CE40-0006 research grant.

Algorithm 5: Compute reduced patterns from the list of all patterns.

Input: L: the list of tuples ( q, s, m 0 , m 1 , m 2 ), one for each output triangle, where q is the shifted point, s the shift and qm k are the three output vertices. Var: Q: a queue. Output: M : the set of reduced patterns. begin

, m 0 , m 1 , m 2 ) ∈ L such that s = s and q + m k = q then if qs , qm 0 , qm 1 , qm 2 ∈ Plane( q, m 0 , m 1 , m 2 )

then // Mark q as a reduced pattern M ← M ∪ q ;