Delaunay property and proximity results of the L-algorithm for digital plane probing

When processing the geometry of digital surfaces (boundaries of voxel sets), linear local structures such as pieces of digital planes play an important role. To capture such geometrical features, plane-probing algorithms have demonstrated their strength: starting from an initial triangle, the digital structure is locally probed to expand the triangle approximating the plane parameters more and more precisely (converging to the exact parameters for inﬁnite digital planes). Among the diﬀerent plane-probing algorithms, the L-algorithm is a plane-probing algorithm variant which takes into account a generally larger neighborhood of points for its update process. We show in this paper that this algorithm has the advantage to guarantee the so-called Delaunay property of the set of probing points, which has interesting consequences: it provides a minimal basis of the plane and guarantees an as-local-as-possible computation.


Introduction
In digital geometry, a digital surface is a quadrangular mesh that corresponds to the boundary of a union of regularly spaced unit cubes.The digital nature of such surfaces is a great advantage for computations in many material sciences or medical imaging applications (e.g., [1,2,3]).However, the local geometry is very poor and difficult to analyze since the cubes provide at most six different normal 1 directions.In order to estimate the geometry of the digital surfaces, an approach is to analyze locally the digital surface within a given neighborhood.The size of the neighborhood should be chosen carefully, because we risk blurring sharp features if the neighborhood is too large.The challenge is to seek the right tradeoff between finding an appropriate geometrical estimation and preservation of sharp features.
Purely digital methods have thus emerged and try to perform digital surface analysis.Geometrical properties of continuous planes are translated into digital planes [4].For example, [5] introduces the concepts of leaning points and leaning plane for digital plane recognition.Other works propose digital plane recognition algorithms with low complexity from a finite subset of three-dimensional integer points, e.g., [6,7].To estimate differential quantities, there exist methods that require user-defined parameters such as Voronoi-based methods [8] and integral invariant methods [9].In this context, plane-probing algorithms could lead to normal vector estimators without the need of external parameters.
Plane-probing algorithms are methods which adapt the neighborhood progressively.The first plane-probing algorithm was proposed in [10].It probes some points in the digital surface and the output represent locally an approximation of the digital surface.Other plane-probing algorithms were proposed later [11,12].They consider a tetrahedron whose apex is outside the digital surface and a triangle formed of three points that belong to the digital surface.
The apex of the tetrahedron always stays in the same position so that one can focus on the movement of the triangle.See Fig. 1 for an example on a digital plane.
Despite the fact that plane-probing algorithms return the exact normal on an infinite digital plane, they may encounter difficulties with non-planar surfaces [12].Probing for points not too far away from the initial point can alleviate those difficulties.That is why, we wish here to estimate a minimal neighborhood that provides a good normal estimation.
There are mainly three types of tetrahedron-based plane-probing algorithms: H-algorithm and R-algorithm, first introduced in [12], and L-algorithm [13].The Figure 1: The evolution (from left to right) of a tetrahedra-based plane-probing algorithm on a digital plane of normal (1,2,5).
primary difference between those algorithms is the considered candidate set of points at each iteration.Sometimes, one can observe that the R-algorithm probes points more locally than the H-algorithm (see Fig. 2 and Fig. 3).Fur-40 thermore, one can observe that the L-algorithm usually takes fewer steps than the two other algorithms (see Fig. 4).Nevertheless, in practice, the R-algorithm and the L-algorithm always return the same triangle at the end and that final triangle has only acute or right angles.It is not trivial at all to give a proof of that fact because the triangles may have an obtuse angle throughout the 45 iterations.Here, L-algorithm has the same output as the R-algorithm.Every triangle of the evolution is superimposed.The initial triangle is blue.The last one is red.
In order to find an invariant, we focus in this paper on pairs of consecutive triangles.For all three above-mentioned algorithms, two consecutive triangles share two common vertices and with the two other distinct vertices one can define a unique circumscribing ball.These balls have very strong and interesting 50 properties in the case of the L-Algorithm.For example, the radius of the balls is strictly increasing (see Fig. 4).In addition, the interior of the balls do not include any point of the digital plane.For normal vectors between ( and (80, 80, 80), we have counted the number of points which lie in both the plane and the ball circumscribing two consecutive triangles.There are in total 55 75235972 points detected in balls for the H-algorithm; 424 points are detected for the R-algorithm; while no points for the L-algorithm.We name this invariant property the Delaunay property after the empty-circle condition used to define the Delaunay triangulation of a set of points in computational geometry.In particular, all the triangular facets whose circumcircle contains none of the 60 points of a digital segment have been characterized in [14].This work can be thought of a partial 3D extension of that result.
Indeed, our objective is to demonstrate that the Delaunay property holds for the L-algorithm and to provide a theoretical upper bound for the position of the final triangle.The Delaunay property is also useful for an optimized 65 implementation of the L-algorithm [13].
The outline of the paper is as follows: we start by introducing the general framework of plane-probing algorithms and the L-algorithm in sec.2. The main theorem that mentions the Delaunay property is announced in sec.3, followed by some proximity results for the L-algorithm.We complete the paper with two sections: in sec.4, we prove the important Lemma 2 leaving the technical details to sec.5, which is arranged into three categories: projection-based-results ( , and proximity results (5.3).

Digital plane probing and the L-algorithm
A digital plane is an infinite digital set defined by a normal N ∈ Z 3 \ {0}, a shift value µ ∈ Z and a thickness ω ∈ Z as follows [15]: In this paper, we set ω := N 1 and we assume w.l.o.g. that µ = 0 and that the components of N are positive, i.e., N ∈ N 3 \ {0}.Given a digital plane of unknown normal vector, a plane-probing algorithm computes the normal vector N of P by sparsely probing it with the predicate "is x in P?".We describe below a specific plane-probing algorithm, called Lalgorithm (see Algorithm 1).An extensive description of the L-algorithm and its properties is provided in [13].For the sake of clarity, we briefly describe its critical steps below.
x, x ∈ H (i) + , we say that x is closer to T (i) than x, denoted by x ≤ T x, if and only if (B(T (i) , x ) ∩ H + ) (see Fig. 5).
x ∈ B(T (i) , x ) x ∈ ∂B(T (i) , x) x ∈ ∂B(T (i) , x) The L-algorithm updates a vertex of T (i) with a point of the set N Note that the triple (k, α, β) may not be unique when several points are in a cospherical position and in this case any triple could be picked for the next 95 iteration.
Once a triple (k, α, β) has been selected, the update rule is: As shown in Algorithm 1, lines 5 to 7, equations (3) and ( 4) are used to update the current triangle.
Termination.The algorithm terminates at a step n, when the neighborhood has an empty intersection with the plane, i.e., when N This result can be found in [12,Theorem 1].Even if [12] only considers slightly different candidate sets, all mentioned results are valid for the larger candidate set we consider in this paper.In addition, if o is one of the lowest point in P, i.e., o • N = 0, T (n) lines up with P, as recalled in the following theorem: and any two edges form a basis of the lattice of upper leaning points, i.e., the
This shows that, for all steps i ∈ {0, . . ., n}, m is a basis of Z 3 , which is especially useful in sec.4.
Theorem 3 ([12], Lemma 5).For all i ∈ {0, . . ., n}, ∀k ∈ Z/3Z, ), then we have v k • N by Theorem 3. In other words, the algorithm always replaces a vertex with a higher candidate point in direction N.That property is a key point in the proof of Theorem 1.It also implies that the set N (i) S ∩ P is always finite.Indeed, the scalar product (v )) • N tends to infinity when α or β (or both) tend to infinity.That is to say, when α or β is large enough, the point k+2 ) does not belong to P. As a consequence, Algorithm 1 can be implemented naively by visiting every point in the set N (i) S ∩ P on lines 4 and 5.A more efficient algorithm is proposed in [13], where a smaller subset is shown to be enough to find a closest point.
In addition, we can derive the following small lemma which is needed in sec.4. Lemma 1.For all i ∈ {0, . . ., n}, Proof.By definition p (0) = o and o is assumed to belong to P. As a consequence, p (0) • N ≥ 0.

Main results
For convenience let T (−1) denotes the degenerated triangle whose three vertices are all at o.At every step i ∈ N, let B (i) be the ball uniquely determined by the four distinct points of T (i−1) ∪ T (i) .We have experimentally observed that the following property is verified by the L-Algorithm, but neither by the H-Algorithm, nor by the R-Algorithm: Property 1 (Delaunay property for plane-probing algorithms).For all i ∈ {0, . . ., n}, the ball B (i) does not contain any point of P in its interior.
If we randomly pick a point from the set N (i) S ∩ P at each iteration, our procedure would still terminate and return a triangle whose normal is equal to the normal of the plane.However, that triangle might have a very bad aspect ratio and its vertices might lie far away the starting point.The Delaunay property is a strong geometric result ensuring that the last triangle have only acute or right angles (Corollary 1, in sec.3.1) and that its vertices stay close enough to the starting point (sec.3.2).Lemma 2. For all i ∈ {0, . . ., n − 1}, if the interior of B (i) contains no point of P, then the interior of B (i+1) contains no point of P ∩ H (i) + .
For a fixed step i, we partition the points of H (i) + into different categories according to their position and we treat each case with distinct lemmas before concluding.Since we now focus on a step i ∈ {0, . . ., n−1}, for sake of simplicity, we drop the exponent (i) in the notations of this section.

Outline of the proof and notations
Remind that p is equal to q − k m k .We conveniently describe any integer point y ∈ Z 3 as a linear combination of m 0 , m 1 and m 2 , which form a basis of To check that any y ∈ H + is in one of the previous cases, it is enough to consider the partition of Z 3 into eight octants depending on the signs of the coefficients and with a convention for null coefficients (see Fig. 6).The negative octant, in red, does not intersect H + and is therefore discarded.The positive octant is itself divided into two regions, the interior, in yellow, corresponds to (Case 1), whereas the boundary faces, in green, corresponds to (Case 2).Among the last six octants, three of them, in blue, correspond to (Case 3), whereas the other three, in purple, correspond to (Case 4).The proofs of the following lemmas require a lot of technical details that are postponed in sec.5 for the sake of readability.They also require the introduction of new notations (see Fig. 7): For sake of clarity, we use the bar notation whenever a scalar product with N is required, i.e., y instead of y • N for any vector y ∈ Z 3 .
That is why we will only check if p + k c k m k < ω, whenever we want to determine whether such a point is in P or not.
Finally, let Σ be the set of all permutations over {0, 1, 2}.Permutations will be useful to describe in a uniform way the various sign combinations of the coefficients.

(Case 1)
The following lemma indicates that the points y corresponding to (Case 1) do not need to be considered because they are not in P.

(Case 3)
This subsection contains Lemma 4 and Lemma 5 that focus on (Case 3).
More precisely, they indicate that the points y corresponding to (Case 3) do not need to be considered because if they are in P, then there is at least one specific point x ∈ N S ∩ P (Lemma 4) such that x ≤ T y (Lemma 5).Proof.We assume w.l.o.g. that σ is the identity, i.e., σ(0) = 0, σ(1) = 1 and Since y ∈ P, we have Since q = ω, the last inequality is equivalent to k (c k − 1)m k < 0.
With h set to min (m 1 , m 2 ) and noticing that c 0 < 0 ⇔ −(c 0 − 1) > 1, we equivalently have In addition, we have which means that h < m 0 .
As a consequence, Those results are used in a second step to complete the proof: (11), ( 12), ( 13), ( 14) are used in cases (i) and (ii), while ( 12) and ( 15) are used in case (iii).We can similarly get ( 12) and ( 15).To explain why, we focus on the case where p + m 1 + 2m 2 is assumed to be in P because the other case is symmetric.
Note first that p + m 0 ∈ P (using the same arguments as in the previous paragraph for p + 2m 2 ).As a consequence, both p + m 1 + 2m 2 and p + m 0 are not in B by hypothesis.We can then apply Lemma 14 with the two vectors m 2 , (−u) and the point v 0 as origin.Since v 0 and v 0 + m 2 − u = v 1 are indeed on the boundary of B, we get m 2 • (−u) < 0 and thus (12).13) is actually a simple consequence of (12).
It remains (14), whose proof is separated into two distinct cases.
are on the boundary of B, we deduce that the point necessarily in the interior of B.Moreover, since no point of B belongs to P, we deduce that q − u is not in P.

(Case 4)
This subsection contains Lemma 6 and Lemma 7 that focus on (Case 4).
More precisely, they indicate that the points y corresponding to (Case 4) do not need to be considered because, as in the previous section, if they are in P, then there is at least one specific point x ∈ N S ∩P (Lemma 6) such that x ≤ T y (Lemma 7).
Proof.We assume w.l.o.g. that σ is the identity.
Since y ∈ P, we have With h set to max (m 0 , m 1 ) and noting that (c 2 − 1) ≥ 2 (since k c k ≥ 3 and c 0 , c 1 ≤ 0), we equivalently have In addition, we have which means that m 2 < h.
We conclude that if h = m 0 (resp.h = m 1 ), p + m 1 + 2m 2 (resp.p + m 0 + 2m 2 ) is strictly smaller than p + k m k = q = ω and thus, the point p + m 1 + 2m 2 (resp.p + m 0 + 2m 2 ) is in P. A fortiori and whatever h is, If y ∈ P and if the interior of B contains no point of P, then there exists a point x ∈ N S ∩ P such that x ≤ T y.
Proof.We assume w.l.o.g. that σ is the identity.
Since y ∈ P, p + 2m 2 ∈ P by Lemma 6.That point, which is also at We assume below that p + m 1 + 2m 2 ∈ P, because the case where only p + 2m 1 + m 2 ∈ P can be proven symmetrically.
One can check that

Conclusion of the proof
Now we have all the material required to prove Lemma 2: Proof.For all i ∈ {0, . . ., n − 1}, the interior of B (i) is assumed to contain no point of P.
Let x be the point chosen by the algorithm at step i, i.e., x = T (i+1) \T (i) .
We want to show that x ≤ T y, for all y ∈ P ∩ H (i) + .Let y be denoted as + .Moreover, since y ∈ P, the coefficients cannot be all strictly positive by Lemma 3. if one coefficient is zero and the others are strictly positive (Case 2), then x ≤ T y by the design of the algorithm, if one coefficient is strictly negative and the others are strictly positive (Case 3), then there exists a point x ∈ N S ∩ P such that x ≤ T y by Lemma 5.
Then, x ≤ T x by the design of the algorithm and x ≤ T y by transitivity.if one coefficient is strictly positive and the others are strictly negative or null (Case 4), then, similarly, there exists a point x ∈ N S ∩ P such that x ≤ T y by Lemma 7.Then, x ≤ T x by the design of the algorithm and x ≤ T y by transitivity.
Since there is no other possibility, the proof is complete.

Technical details
The proof of Lemma 2 refers to several technical details that we elaborate in this section.The results are organized into three categories.First, we present in sec.5.1 several useful angle relations in the tetrahedron formed by the current triangle and the fixed point q.These results come from the fact that q always projects inside the current triangle.Then, we present several general and purely geometrical circumsphere-based properties in sec.5.2, because the relation ≤ T and the selection of a closest point according to ≤ T involves circumspheres.
Finally, in sec.5.3, we derive in an algebraic way several other results about the comparison of specific points according to ≤ T .These results are used in Lemma 5 and Lemma 7, which are the main ingredients in the proof of Lemma 2.

Projection-based results
Remind that k is taken modulo 3. To keep notations short, we simply write ∀k instead of ∀k ∈ Z/3Z in this section.Let us introduce the following extra notations (see Fig. 7): ∀i ∈ {0, . . ., n}, Note that the following equality also holds for the estimated normal vector, which is normal to the current triangle: Lemma 8.For all i ∈ {0, . . ., n}, ∀k, Base Case.The triangle T (0) and q forms a trirectangular tetrahedron.We have ∀k, Induction case.We now assume that for any i ∈ {0, . . ., n−1}, ∀k, By the update rule, equation 4, we straightforwardly have: Since we have ∀k, ), the induction hypothesis implies the result.

415
From now on, we omit once again the exponent (i) for clarity.We go on with this purely geometrical result (see Fig. Proof.We first expand (N×d)•m < 0, which is equivalent to (d×m)•(d×d ) > 0, using the scalar quadruple product rule: We then similarly expand ( Multiplying both sides by d 2 and (d • d ) leads to: We now combine the two preceding lemmas to find angular relations in the tetrahedron formed by the current triangle and q, i.e., involving the vectors (m k ) and (d k ).See Fig. 10.These results are used in Lemma 5 and in sec.5.3.We show below that the three terms are positive, so is the whole sum.

Conclusion
We introduce the Delaunay property for plane-probing algorithms and prove that the L-algorithm verifies the Delaunay property.We invoke several geometry properties related to projections and spheres in order to proceed to the proof by recurrence.Since L-algorithm verifies the Delaunay property, a direct consequence is that the output triangle contains a minimal basis for the underlying rank-two lattice.We also show that such minimal basis provides a raw estimation of the bound for the distance between the vertex of the last triangle to the fixed apex of tetrahedron.
As for future work, our study serves as a basis for the study of the locality for the L-algorithm.The progress of the algorithm depends on every point that it visits during the iteration.We hope that we can find a bound for the convex hull of every point visited by the algorithm in order to optimize the performance of the algorithm on digital surfaces.On the other hand, we also wish to understand theoretically how the R-algorithm always returns the

Figure 4 :
Figure 4: Evolution of the radius of the balls circumscribing two consecutive triangles for a digital plane of normal (2, 5, 156).

P
that is a closest one according to ≤ T .More precisely, if N (i) S P = ∅, there is an index k ∈ Z/3Z and there are numbers (α, β) ∈ N 2 \ (0, 0) such that ∀x ∈ N (i) (n) S P = ∅ (Algorithm 1, line 3).The number of steps, n, is less than or equal to ω − 3, which is a tight bound reached for any normal of components (1, 1, r) with r ∈ N \ {0}.

Z 3 (
by Theorem 2), i.e., y := p + k c k m k , with c k ∈ Z for all k ∈ Z/3Z.By construction, the bounding plane of H + is defined by the vertices {p + m 0 + m 1 , p + m 1 + m 2 , p + m 0 + m 2 }.All lattice points y on this plane are such that k c k = 2. Hence, for any lattice point y, we have y ∈ H + ⇔ k c k ≥ 3.In this section, we always assume that y ∈ H + .We consider several cases: (Case 1) the coefficients c 0 , c 1 , c 2 are all strictly positive (see Lemma 3), (Case 2) one coefficient is zero and the others are strictly positive; these points are exactly the ones probed in the L-algorithm (see also the definition of the candidate points, equation (2)), (Case 3) one coefficient is strictly negative and the others are strictly positive (see Lemma 4 and Lemma 5), (Case 4) one coefficient is strictly positive and the others are strictly negative or null (see Lemma 6 and Lemma 7).

Figure 6 :
Figure 6: The discrete space Z 3 (intersected with the box [−5, 5] 3 for the illustration) is partitioned into five regions: the yellow, green, blue and purple regions respectively correspond to (Case 1), (Case 2), (Case 3) and (Case 4), the red one is discarded because none of its points lie in H + (the three black arrows indicate the direction of the grid axes).

Figure 7 :
Figure 7: Notations for v k (black), m k (red) and d k (blue).Note that Nk (grey) is used only in sec.5.1.

9 ): Lemma 9 .Figure 9 :
Figure 9: Illustration of Lemma 9. Note that m does not belong to the span of d and d .However, it projects along N into the interior of the convex combination of d and d (hatched area).

)
If d • d = 0, we can conclude from (17) for d • m and from (18) for d • m.If not, then d • d > 0 and we can derive lower and upper bounds for d • m, respectively from (17) and (18):