Xiuheng Wang
email: xiuheng.wang@oca.eu

Ricardo Augusto Borsoi
email: ricardo.borsoi@univ-lorraine.fr

Cédric Richard
email: cedric.richard@unice.fr

Online change point detection on Riemannian manifolds with karcher mean estimates

Keywords: Online change point detection, nonparametric, Riemannian manifolds, Karcher mean

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Change-point detection (CPD) aims to detect abrupt changes in the data distribution, and is recognized as one of the most significant tasks in time series data analysis. Despite the huge literature on offline CPD, online CPD still suffers from major challenges while it plays a fundamental role in a wide range of applications such as audio [START_REF] Bietti | An online EM algorithm in hidden (semi-) Markov models for audio segmentation and clustering[END_REF] and video [START_REF] St-Charles | Subsense: A universal change detection method with local adaptive sensitivity[END_REF] segmentation, medical condition monitoring [START_REF] Malladi | Online Bayesian change point detection algorithms for segmentation of epileptic activity[END_REF], or human behavior analysis [START_REF] Aminikhanghahi | Real-time change point detection with application to smart home time series data[END_REF] to cite a few. Considering whether prior knowledge about the data distributions is available or not, online CPD approaches can be divided into parametric and nonparametric. Examples of parametric methods include the cumulative sum (CUSUM) [START_REF] Inclan | Use of cumulative sums of squares for retrospective detection of changes of variance[END_REF] and the generalized likelihood ratio test (GLRT) [START_REF] Gustafsson | The marginalized likelihood ratio test for detecting abrupt changes[END_REF]. Such methods assume that the distribution of the data belongs to a known parametric family. However, knowledge of the data distribution is not always available, making the use of non-parametric methods necessary. Classic non-parametric strategies include monitoring changes in the mean or the variance (e.g., the exponentially weighted The work of C. Richard was funded in part by the ANR under grant ANR-19-CE48-0002, and by the 3IA Côte d'Azur Senior Chair program. moving average -EWMA) [START_REF] Costa | A single ewma chart for monitoring process mean and process variance[END_REF] or in a generalized statistic [START_REF] Gretton | A kernel method for the two-sample-problem[END_REF] of the data stream. In [START_REF] Keriven | NEWMA: a new method for scalable model-free online change-point detection[END_REF], a new algorithm called NEWMA was proposed. It consists of comparing two EWMA of the data statistics computed using two distinct forgetting factors to detect change points without requiring the storage of old samples. A non-parametric online algorithm was designed in [START_REF] Ferrari | Online change-point detection with kernels[END_REF] based on an adaptive kernel-based density ratio estimation. Recently, deep learning has been also considered in non-parametric online CPD [START_REF] Atashgahi | Memory-free online change-point detection: A novel neural network approach[END_REF][START_REF] Wang | Change point detection with neural online density-ratio estimator[END_REF].

Modern signal processing tasks increasingly tackle data that does not reside in Euclidean spaces, such as graphs or categorical data. In particular, Riemannian manifolds [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF][START_REF]An introduction to optimization on smooth manifolds[END_REF] have drawn significant attention due to their widespread applications, including diffusion tensor imaging [START_REF] Pennec | A riemannian framework for tensor computing[END_REF] and pedestrian detection [START_REF] Tuzel | Pedestrian detection via classification on riemannian manifolds[END_REF]. However, developing methods that can process manifold-valued data is still challenging since one has to account for the nonlinear geometry of the space. Moreover, manifolds lack a vector space structure, which makes it hard to generalize algorithms originally developed for Euclidean spaces. Aside from some online CPD algorithms that have been extended to particular non-Euclidean domains such as graphs [START_REF] Ferrari | Non-parametric community change-points detection in streaming graph signals[END_REF][START_REF] Borsoi | Online graph-based change point detection in multiband image sequences[END_REF] or categorical data [START_REF] Ienco | Change detection in categorical evolving data streams[END_REF], few works have investigated manifold-valued data. For instance, an online CPD algorithm in [START_REF] Bouchard | Riemannian geometry for compound gaussian distributions: Application to recursive change detection[END_REF] was specifically designed for the compound Gaussian distribution. This technique is however parametric and not broadly applicable. An example of non-parametric techniques can be found in [START_REF] Dubey | Fréchet change-point detection[END_REF], but it is only able to detect a single change point in an offline manner.

In this work, we introduce a unified framework for online CPD on Riemannian manifolds based on Karcher mean estimation. Specifically, a non-parametric strategy is considered by monitoring the Karcher mean of manifold-valued data, which is estimated efficiently in an online way using a Riemannian stochastic gradient descent (SGD) algorithm. To detect abrupt change points, two Karcher mean estimates with different step sizes -one which takes longer to converge and focus on a long-term trend, and one which converges faster to assimilate change points quickly -are compared to form a test statistic. We then illustrate the proposed framework to detect changes on the Riemannian manifold of symmetric positive definite (SPD) matrices. Experimental results on sequences of SPD matrices generated both synthetically and as feature descriptors of video sequences demonstrate its effectiveness.

BACKGROUND

This section introduces some basics of Riemannian geometry with geometrical tools for optimization. Extensive presentations can be found in [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF][START_REF]An introduction to optimization on smooth manifolds[END_REF]. Then this section focuses on a usual example of Riemannian manifolds, that of p × p SPD matrices, denoted as S ++ p , with an invariant metric.

Riemannian geometry and optimization

A Riemannian manifold (M, g) is defined by a constrained set M equipped with a Riemannian metric g x (•, •), that is, a smoothly varying inner product ⟨•, •⟩ x : T x M × T x M → R, defined for each x ∈ M, where T x M is called the tangent space of M at x. The length of a parameterized curve, say, c : [a, b] → M, is given by:

L(c) = b a g( ċ(α), ċ(α))dα = b a ∥ ċ(α)∥dα, (1)
where ċ = ∂c/∂α is the velocity of c. This allows us to define the geodesic γ : [0, 1] → M, which is the unique curve of minimal length linking x and y, with x = γ(0) and y = γ(1).

The Riemannian distance d M (•, •) : M × M → R is defined as follows: d M (x, y) = 1 0 ∥ γ(α)∥dα. (2)
Note that d M satisfies all conditions to be a metric.

The exponential map w = exp x (v) defines the point w of M located on the unique geodesic

γ v (t) such that γ v (0) = x, γ ′ v (0) = v and γ v (1) = w.
The inverse of the exponential map is defined as v = exp -1 x (w). Since the exponential map may be hard to compute, one often resorts to a second-order approximation, called retraction mapping R x :

T x M → M at x ∈ M, and which satisfies d M (R x (tv), exp x (tv)) = O(t 3 ).
Consider f : M → R a smooth function. The Riemannian gradient of f at x ∈ M is the defined as the unique tangent vector ∇f (x) ∈ T x M satisfying:

d dt t=0 f (exp x (tv)) = ⟨∇f (x), v⟩ x (3)
for all v ∈ T x M.

The Riemannian manifold of SPD matrices

The geodesic distance of S ++ p between two SPD matrices Σ 1 and Σ 2 ∈ S ++ p can be computed in closed form [START_REF] Pennec | A riemannian framework for tensor computing[END_REF] as:

d S ++ p (Σ 1 , Σ 2 ) = log(Σ -1 2 2 Σ 1 Σ -1 2 2 ) F , (4)
where ∥•∥ F denotes the Frobenius norm. The Riemannian gradient ∇f at Σ ∈ S ++ p is given by:

∇f (Σ) = Σ sym(G)Σ, (5)
with G ∈ R p×p the Euclidean gradient of function f at Σ and sym(G) = 1 2 (G T + G). In practice, the Euclidean gradient can be easily computed using automatic differentiation tools.

Let ξ ∈ T Σ S ++ p . A retraction R Σ,S ++ p : T Σ S ++ p → S ++ p is defined as R Σ,S ++ p (ξ) = Σ + ξ + 1 2 ξΣ -1 ξ. ( 6
)
This retraction is a second-order approximation of the exponential mapping, that is,

exp Σ,S ++ p (tξ) = R Σ,S ++ p (tξ) + O(t 3 ). (7)

METHODOLOGY

Let us consider a time series of independent random variables {x t } t∈N lying on a Riemannian manifold M. We assume that there exists a time index t r ∈ N with an abrupt change in the probability distribution of x t , that is,

t < t r : x t ∼ P 1 (x) , t ≥ t r : x t ∼ P 2 (x), (8)
where P 1 (x) and P 2 (x) denote two different probability measures on M [START_REF] Pennec | Probabilities and statistics on riemannian manifolds: A geometric approach[END_REF] which represent the distribution of data x t before and after the change point t r . Note that, to simplify the presentation, (8) considers only a single change point. However, the algorithm presented hereafter can handle multiple change points. CPD algorithms aim to estimate a change point tr according to two complementary objectives: 1) minimizing the detection delay, i.e., tr -t r for tr being the first detection after t r ; and 2) minimizing the probability of false alarms, i.e., of flagging some t ̸ = t r as a change point. In this work, we consider a problem setting where manifold-valued data x t are observed sequentially over time and change points must be detected online. This means that we need to decide if each time instant t ∈ IN is a change point based only on past data {x t ′ } t ′ ≤t . Moreover, unlike [START_REF] Bouchard | Riemannian geometry for compound gaussian distributions: Application to recursive change detection[END_REF], we focus on nonparametric strategies, which do not make additional assumptions about the statistical distribution of the data.

Non-parametric statistics and the Karcher mean

As discussed above, we focus on non-parametric strategies, where there is no prior knowledge about probability measures of the data stream. In Euclidean spaces, this can be done by monitoring changes in the mean or the variance [START_REF] Costa | A single ewma chart for monitoring process mean and process variance[END_REF], or in a generalized statistic [START_REF] Gretton | A kernel method for the two-sample-problem[END_REF] of the data stream. In order to generalize such strategies to Riemannian manifolds, we consider the Karcher mean [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF] of the data stream x t ∈ M, which provides a generalization of the center of mass from Euclidean domains to a manifold M.

The Karcher mean is a generalization of the Fréchet mean, which is defined as the set of values that (globally) minimize the expected variance:

f (m) = E x∼P (x) d 2 M (m, x) = d 2 M (m, x)dP (x)
of the Riemannian distance d M , that is,

m * ∈ arg min m f (m) . (9)
Note that the existence and uniqueness of the Fréchet mean are not guaranteed. The Karcher mean relaxes this definition by considering the local optima of f (m), instead of only the global one. This allows to establish existence and uniqueness conditions [START_REF] Kendall | Probability, convexity, and harmonic maps with small image I: uniqueness and fine existence[END_REF], and also makes it possible to compute m by locally solving (9) using Riemannian optimization methods [START_REF] Pennec | Probabilities and statistics on riemannian manifolds: A geometric approach[END_REF]. The Karcher mean is unique in many manifolds, such as those that are connected and have nonpositive curvature [START_REF] Afsari | Riemannian ℓ p center of mass: existence, uniqueness, and convexity[END_REF], which includes S ++ p . The CPD strategy on manifolds proposed in the following monitors abrupt changes in the Karcher mean of the data stream. An important requirement is that change points must be detected in an online way, that is, only based on past data. Consequently, we start the presentation with an online Riemannian descent algorithm to estimate the Karcher mean of streaming data. This will be an integral part of the CPD strategy presented afterward.

Online estimation of the Karcher mean

In optimization problem [START_REF] Keriven | NEWMA: a new method for scalable model-free online change-point detection[END_REF], the cost function cannot be computed explicitly because P (x) is unknown. However, observations {x t } are available to compute function d 2 M (m, x t ) for any parameter m and data point x t . That function can be viewed as a stochastic approximation of the loss f (m) updated with new input data x t . Consequently, we consider using the Riemannian SGD algorithm [START_REF] Bonnabel | Stochastic gradient descent on riemannian manifolds[END_REF] to address problem [START_REF] Keriven | NEWMA: a new method for scalable model-free online change-point detection[END_REF]. On manifold M, the update of m with a step size α is given by:

m t+1 = exp mt -αH(m t , x t ) , (10)
where exp m is the exponential map at m, and H(m, x) denotes the Riemannian gradient of the loss such that

E x∼P (x) H(m, x) = H(m, x)dP (x) = ∇f (m).
For computational simplicity, we replace the exponential map in (10) by a retraction R mt . This yields the alternative update:

m t+1 = R mt -αH(m t , x t ) . (11)
Algorithm 1: Online CPD on manifolds Input: {xt}, step sizes λ, Λ, threshold ξ. 1 Initialization: m λ,0 = mΛ,0 = x0 ; 2 for t = 1, 2, 3, . . . do 3 Update the "fast" and "slow" Karcher mean estimates m λ,t and mΛ,t using ( 12) and ( 13) ;

An adaptive CPD with the Karcher mean

We aim to detect change points by monitoring abrupt changes in m over time, that is, a point t ′ is labeled as a change point if m changed abruptly at t ′ . This requires knowledge of two quantities of interest, m bef and m aft , which correspond to the Karcher mean before and after a candidate change point t ′ . First, we propose to compute estimates of these values, say m bef and m aft . Then, as a test statistic, we propose to compare these two quantities using the Riemannian distance, that is, d M ( m bef , m aft ): the larger the Riemannian distance between the Karcher mean estimates before and after instant t ′ , the more likely we are to flag t ′ as a change point.

The question is how to calculate m bef and m aft efficiently and in an online way. Previous work proposed to partition a data stream of length N into two segments, {1, . . . , t ′ -1} and {t ′ , . . . , N } for every t ′ , and testing for differences between their Karcher mean and variance [START_REF] Dubey | Fréchet change-point detection[END_REF]. However, this strategy cannot process data streams on-the-fly or detect multiple change points. In [START_REF] Keriven | NEWMA: a new method for scalable model-free online change-point detection[END_REF], within the realm of Euclidean geometry, these estimates were computed using two exponentially weighted moving averages with different forgetting factors. Nevertheless, this principle cannot be transposed directly to Riemannian manifold setting. Instead, we propose to use two estimates provided by Riemannian stochastic gradient descent algorithms as presented in Section 3.2, with two different fixed step sizes λ < Λ. The Karcher means are updated according to [START_REF] Atashgahi | Memory-free online change-point detection: A novel neural network approach[END_REF] as follows:

m λ,t+1 = R m λ,t -λH(m λ,t , x t ) , (12)
m Λ,t+1 = R mΛ,t -ΛH(m Λ,t , x t ) , (13)
with initialization m λ,0 = m Λ,0 = x 0 . Convergence of the updates ( 12) and ( 13) is directly affected by λ and Λ. Constraint λ < Λ thus means that m Λ,t is more likely to adapt to new data and approximates m aft , while m λ,t has longer memory and is more suitable to estimate baseline trend m bef . Using these two estimates, we can define an adaptive CPD statistic by comparing the difference between m λ,t and m Λ,t using the Riemannian distance on M as: CPD is then performed by comparing g t to a threshold ξ. The full CPD procedure is summarized in Algorithm 1.

g t = d M (m λ,t , m Λ,t ) . (14

EXPERIMENTS

We shall now present some experiments on manifold S ++ p discussed in Section 2. We consider the problem of detecting change points with Algorithm 1 both in a sequence of synthetic SPD matrices, and in a sequence of region covariance descriptors derived from a real video. With {Σ t } t∈N lying on S ++ p and the metric defined in (4), the Karcher means were estimated by minimizing the objective function

f (Σ) = E S∼P (S) log(S -1 2 ΣS -1 2 ) 2 F
using the Riemannian SGD algorithms in [START_REF] Wang | Change point detection with neural online density-ratio estimator[END_REF] and [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF] and the stochastic approximation S ≃ Σ t at each t. They were used to compute the online CPD statistic in [START_REF]An introduction to optimization on smooth manifolds[END_REF].

We compared our method with two baseline algorithms, namely, NEWMA [START_REF] Keriven | NEWMA: a new method for scalable model-free online change-point detection[END_REF] and the Fréchet CPD (F-CPD) [START_REF] Dubey | Fréchet change-point detection[END_REF]. On the one hand, since NEWMA was originally designed for Euclidean spaces, we applied it to the vectorization of the left triangular and diagonal parts of each SPD matrix. On the other hand, F-CPD was designed for manifold-valued data but can only detect a single change point in an offline manner. We tackled these issues by using F-CPD over consecutive sliding windows of length 100. We set λ = 0.01 and Λ = 0.02 for synthetic data, and λ = 0.05 and Λ = 0.06 for real data.

Experiment with synthetic data: The synthetic matrices Σ t ∈ S ++ p with p = 6 were sampled from a Wishart distribution with the scaling matrix V and d degrees of freedom. We generated 800 samples and set a change point at t r = 500 where we reset V . Fig. 1 shows the Receiver Operating Characteristic (ROC) curves of all methods for 5000 Monte Carlo runs. It can be seen that our method achieved a significant improvement in the detection rate with a low rate of false alarms when compared to both NEWMA, which does not take the manifold geometry into account, and to F-CPD, which was designed to operate offline.

Experiment with real data: To further evaluate our approach, we made use of a real video of an outdoor scene from [START_REF] Goyette | Changedetection.net: A new change detection benchmark dataset[END_REF]. This scene contains intermittent object motions including cars and pedestrians as well as subtle, noninformative changes such as tree leaves moving. The video contains color images with 658×491 pixels. Ground truths of segmentation results for larger moving targets are available. With ground truths, the first T = 1150 sequential frames were selected and cropped to 658 × 260 pixels. We considered detecting change points in disjoint, compact regions of this scene by segmenting each frame into superpixels via the SLIC technique [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF]. A snapshot of this scene with its superpixel decomposition is shown in Fig. 2.

We considered the region covariance descriptor features in [START_REF] Tuzel | Pedestrian detection via classification on riemannian manifolds[END_REF] but with the pixel location information removed. Let {s t } t∈N denote the video sequence. Each frame s t was first processed by computing a feature vector for each pixel: where (x, y) represents pixel locations and I x , I y , I xx , I yy denote intensity derivatives. For each superpixel p containing n pixels, the covariance descriptor was an SPD matrix of the feature vectors z t (x, y) computed as follows:

Σ t = 1 n -1 (x,y)∈p z t (x, y) -zt z t (x, y) -zt T , with Σ t ∈ S ++ p and p = 6, and zt the sample mean of z t (x, y) for all (x, y) ∈ p.

We used Algorithm 1 to detect change points in the data stream {Σ t }. The resulting test statistic for one given superpixel and the corresponding change points are illustrated in Figure 3, and compared to the ground truth. We can observe that all peaks of the test statistic are located near change points pointed out by the ground truth, which indicates that the proposed method can reach a low false alarm rate.

CONCLUSION

In this paper, we presented a general approach for the online detection of change points in Riemannian manifolds based on Karcher mean estimation. An adaptive test statistic was computed by comparing two Karcher means estimated with Riemannian SGD algorithms, one converging faster to assimilate new data, and another one converging more slowly to focus on a long-term trend. Experimental results on the Riemannian manifold of SPD matrices illustrated the superiority of our strategy compared with two baseline algorithms that either operate offline or do not take manifold geometry into account.

4 5 if gt > ξ then 6 Flag

 56 Compute the test statistic gt = dM(m λ,t , mΛ,t) ; t as a change point;

Fig. 1 .

 1 Fig. 1. ROC curves for all compared algorithms.

Fig. 2 .Fig. 3 .

 23 Fig.2. Snapshot of the video sequence with its super-pixels decomposed via the SLIC algorithm[START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF].

  z t (x, y) = |I x |, |I y |, I 2 x + I 2 y , |I xx |, |I yy |, arctan I x I y T


Tristan Roussillon

Combinatorial Generation of Planar Sets *

HAL is

Introduction

Digital geometry mainly deals with sets of discrete elements considered to be digitized versions of Euclidean objects. A digital surface may be seen as a mesh of unit square faces whose vertices have integer coordinates. A challenge is to decompose digital surfaces into patches, such as pieces of digital planes.

A digital plane has been analytically defined as a set of points with integer coordinates lying between two parallel planes [START_REF] Reveillès | Géométrie Discrète, calculs en nombres entiers et algorithmique[END_REF]. The topological properties of such object depend on the thickness of the strip. If one considers the graph whose nodes are points of integer coordinates and whose edges link points which are distant to each other less than a given thickness, one can efficiently generate an arbitrary piece of digital plane using a breadth-first search and the analytical definition as a set-membership predicate. Furthermore, given a finite point set, one can decide whether this set belongs to a digital plane or not in linear time using linear programming. See for instance [START_REF] Françon | Recognizing arithmetic straight lines and planes[END_REF] and, for a review on digital planarity, [START_REF] Brimkov | Digital planarity -a review[END_REF]. However, a linear programming solver does not help so much for the analysis of digital surfaces, because one does not know which point set should be tested to obtain patches that approximate the tangent planes of the underlying surface.

In order to cope with this problem, plane-probing algorithms have been developed [START_REF] Lachaud | An optimized framework for plane-probing algorithms[END_REF][START_REF] Lachaud | Two plane-probing algorithms for the computation of the normal vector to a digital plane[END_REF]. Their main feature is to decide on-the-fly how to probe a given point set and locally align a triangle with it, without requiring any user input. However, if probing for points in a sparse way is perfect for digital planes, it is not enough for non-convex parts where the triangles may jump over holes or stab the digital surface. Therefore, that approach also requires to associate pieces of digital planes to the triangles and check whether they fit the digital surface or not.

In this paper, we make a first step towards the incremental generation of planar patches, which can be used during the execution of a plane-probing algorithm (see Fig. 1). To do that, we take advantage of the combinatorial properties of digital planes.

The particular case of digital lines has been studied for a long time in different contexts and has led to many applications. One key result is that digital lines are hierarchical point sets whose structure is exactly described by the continued fraction expansion of their slope and strongly Figure 1: Local approximation of a digital sphere of radius 63 by planar patches (in green) from two starting points (in blue). The implementation combines the plane-probing algorithm H [START_REF] Lachaud | Two plane-probing algorithms for the computation of the normal vector to a digital plane[END_REF] with the generation method of Section 3. relies on the Euclidean algorithm. See for instance [START_REF] Troesch | Interprétation géométrique de l'algorithme d'Euclide et reconnaissance de segments[END_REF] and, for a survey on digital straightness, [START_REF] Klette | Digital straightness -a review[END_REF].

There has been much effort done in order to find similar results in three dimensions despite the lack of a canonical algorithm and the diversity of existing generalizations of the Euclidean algorithm. Some combinatorial results, involving symmetries, piece exchanges and flips, have been stated thanks to an appropriate representation of digital planes [START_REF] Labbé | A d-dimensional Extension of Christoffel Words[END_REF]. However, most of other related works depend on a multi-dimensional generalization of the Euclidean algorithm such as Brun, e.g., [START_REF] Berthé | Brun expansions of stepped surfaces[END_REF], Jacobi-Perron, e.g., [START_REF] Berthé | A study of Jacobi-Perron boundary words for the generation of discrete planes[END_REF], Fully substractive, e.g. [START_REF] Berthé | Fully Subtractive algorithm, tribonacci numeration and connectedness of discrete planes[END_REF], or a mix of several of them [START_REF] Jamet | Generation of digital planes using generalized continued-fractions algorithms[END_REF]. Those algorithms have been used to generate digital planes from a normal vector. There are two different but closely related construction schemes in the literature. The first one is based on union and translation of point sets [START_REF] Berthé | Critical connectedness of thin arithmetical discrete planes[END_REF][START_REF] Domenjoud | Geometric Palindromic Closure[END_REF][START_REF] Jamet | Generation of digital planes using generalized continued-fractions algorithms[END_REF]. It has been used mostly to construct the thinnest digital plane that is connected [START_REF] Berthé | Fully Subtractive algorithm, tribonacci numeration and connectedness of discrete planes[END_REF][START_REF] Berthé | Critical connectedness of thin arithmetical discrete planes[END_REF][START_REF] Domenjoud | Facet connectedness of arithmetic discrete hyperplanes with non-zero shift[END_REF][START_REF] Domenjoud | Facet connectedness of discrete hyperplanes with zero intercept: The general case[END_REF]. The second one is based on a description of standard digital planes as unions of square faces and uses geometrical extensions of substitutions, i.e. rules that replace square faces by unions of square faces. Since the pioneer work of Arnoux and Ito [START_REF] Arnoux | Pisot substitutions and Rauzy fractals[END_REF], that formalism has been used for instance in [START_REF] Berthé | Brun expansions of stepped surfaces[END_REF][START_REF] Berthé | Critical connectedness of thin arithmetical discrete planes[END_REF][START_REF] Berthé | A study of Jacobi-Perron boundary words for the generation of discrete planes[END_REF][START_REF] Fernique | Multidimensional Sturmian Sequences and Generalized Substitutions[END_REF][START_REF] Fernique | Generation and recognition of digital planes using multi-dimensional continued fractions[END_REF]. Both construction schemes incrementally generate sets so that the current set will be included in the next one, but suffer from topological and geometrical limitations. This work is based on the second scheme that generates few elements in comparison with the first scheme, e.g., [START_REF] Jamet | Generation of digital planes using generalized continued-fractions algorithms[END_REF].

This paper is an extension of [START_REF] Meyron | Approximation of Digital Surfaces by a Hierarchical Set of Planar Patches[END_REF] and is organized as follows. We first present a multidimensional generalization of the Euclidean algorithm in Section 2. It includes a far larger class of existing algorithms than the generalization proposed in [START_REF] Meyron | Approximation of Digital Surfaces by a Hierarchical Set of Planar Patches[END_REF]. Furthermore, we do not only show that plane-probing algorithms look like a geometrical version of a three-dimensional Euclidean algorithm, but also show that all the classical three-dimensional Euclidean algorithms like Brun or Jacobi-Perron can be translated into plane-probing algorithms. In Section 3, we explain how one can generate unions of faces from the output of the above-mentioned algorithms. We also prove several combinatorial properties of the generated sets. Several of them are new, others are stronger results or have a more self-contained proof than the ones included in [START_REF] Meyron | Approximation of Digital Surfaces by a Hierarchical Set of Planar Patches[END_REF].

Generalization of the Euclidean Algorithm

We first propose in Subsection 2.1 a multi-dimensional generalization of the Euclidean algorithm: Algorithm 1. Its validity and termination are then established in Subsection 2.2. Algorithm 1 is linked with several other known algorithms in the last two subsections.

Notations, Definitions and Algorithm

Let {e 1 , . . . , e d } be the canonical basis of R d , with d ≥ 2. We denote 0 the origin and 1 = d i=1 e i the vector with all coordinates equal to 1 whatever the dimension d. Let I d be the d × d identity matrix. Matrices and vectors are written in bold, respectively with upper-case and lower-case letters. Vectors are thought as column vectors, while their transpose, in front of which is put the letter t, are thought as row vectors. For instance, the i-th coordinate of x is equal to t e i x. However, for sake of shortness, t e i x will be also denoted x i .

Let N d be the set of vectors of length d with nonnegative integer coordinates. Let ≤ be the coordinatewise partial order defined such that, for all pairs x, y ∈ N d , x ≤ y if and only if x i ≤ y i for all i ∈ {1, . . . , d}. Endowed with the order relation ≤, N d is a lattice.

The Euclidean algorithm is often used to compute the greatest common divisor of two numbers, but in what follows we are only interested in the sequence of operations performed when the algorithm is executed. Since multiplying the two numbers by a constant does not change the sequence of operations, we can assume that the greatest common divisor is equal to 1. More generally, for d numbers, we can similarly focus on vectors of N d whose coordinates are relatively prime, i.e., on the set P d := {x ∈ N d | gcd (x 1 , . . . , x d ) = 1}. Note that P d contains e 1 , . . . , e d and 1, but not 0.

We now define the set of inputs accepted in Algorithm 1:

for each S ⊂ P d , Ω(S) :

= x∈S {y ∈ P d | x ≤ y}.
We will only consider Ω({1}) and Ω({e 1 , . . . , e d }). In the first set, the coordinates of the vectors are coprime and strictly positive, whereas they are coprime and nonnegative in the second one. Algorithm 1 takes as input a set T ⊂ P d whose elements are called terminal and a vector v ∈ Ω(T ). While v is not terminal, it is iteratively transformed into another vector by matrix multiplication. The matrices must belong to U d , defined as the set of d × d matrices having -1 at the intersection between the k-th row and the l-th column with distinct k, l ∈ {1, . . . , d} and the same entries as I d elsewhere. Note that the elements of U d are matrices with determinant 1 and are thus unimodular. Note that multiplication by a unimodular matrix U preserves vectors of coprime coordinates, i.e., v ∈ P d ⇒ Uv ∈ P d . This fact is used in the next section.

Algorithm 1 Multidimensional Euclidean algorithm.

Require: a terminal set T and a vector v ∈ Ω(T ) Ensure: a list of matrices (U n ) n∈{0,...,N } and a terminal t ∈ T such that t

= U N • • • U 0 v. 1: lst ← {I d } 2: while v / ∈ T do 3: Select U ∈ U d such that Uv ̸ = v and Uv ∈ Ω(T ) 4: lst ← lst ∪ {U} 5:
v ← Uv 6: end while 7: return lst, v

Validity, Termination and Complexity

In this subsection, we answer two questions regarding Algorithm 1: is there always a matrix U ∈ U d satisfying the constraints on line 3? Does the algorithm always terminate?

Before answering, we introduce the following technical result:

Lemma 1. Let T be either equal to {1} or {e 1 , . . . , e d }. Let us consider a vector v ∈ Ω(T ) and a matrix U ∈ U d containing -1 at the intersection of the k-th row and the l-th column, with distinct k, l ∈ {1, . . . , d}.

The following relations involving the coordinates of v,

v l ≥ 1 and, ( 1
) v k -v l ≥ 1, if T = {1}, (2) v k -v l ≥ 0, if T = {e 1 , . . . , e d },
imply both Uv ̸ = v and Uv ∈ Ω(T ).

Proof. The k-th coordinate of Uv is the result of substracting the l-th coordinate from the k-th coordinate of v, i.e., t e k (Uv) = v kv l . The other coordinates of Uv are the same as v.

On one hand, we have Uv ̸ = v if v l ≥ 1. On the other hand, we have Uv ∈ Ω(T ) if and only if Uv ∈ P d and there exists t ∈ T such that t ≤ Uv. The first condition is true because v ∈ P d and U is unimodular. The second condition is true if t e k (Uv) ≥ 1 (resp. t e k (Uv) ≥ 0) if T = {1} (resp. T = {e 1 , . . . , e d }).

Since t e k (Uv) = v kv l , putting together the above results concludes.

We are now able to answer the first question about the existence of a matrix U ∈ U d for two specific values of T .

Lemma 2. Let T be either equal to {1} or {e 1 , . . . , e d }. If v ∈ Ω(T ) \ T , then there exists U ∈ U d such that Uv ̸ = v and Uv ∈ Ω(T ).

Proof. Since there are different corner cases for the two possible values of T , we separately deal with them.

Case T = {1}: If v ∈ Ω(T ) \ T , there are two consequences. First, 1 ≤ v, which means that the coordinates of v are greater than or equal to 1. Second, v ∈ P d and v ̸ = 1, which means that the coordinates of v are not all equal. As a result, there exist two disinct k, l ∈ {1, . . . , d} such that v l ≥ 1 and v kv l ≥ 1.

Case T = {e 1 , . . . , e d }: Again, if v ∈ Ω(T ) \ T , there are two consequences, but slightly different. First, there is t ∈ {e 1 , . . . , e d } such that t ≤ v, which means that all coordinates of v are greater than or equal to 0. Second, v ∈ P d and v ̸ = t, which means that v contains at least two nonzero coordinates. As a result, there exist two disinct k, l ∈ {1, . . . , d} such that v l ≥ 1 and v kv l ≥ 0.

In both cases, using Lemma 1, we conclude that there exists a matrix U ∈ U d , whose entry (k, l) equals -1, such that Uv ̸ = v and Uv ∈ Ω(T ).

We are now able to study the termination and complexity of Algorithm 1 for two specific values of T . Theorem 3. Let T be either equal to {1} or {e 1 , . . . , e d }. Let ω be equal to the 1-norm of the input vector v. Algorithm 1 terminates after less than ω iterations.

Proof. The coordinates of v are always nonnegative and the matrix multiplication consists in subtracting a nonzero coordinate from another. The integral number ∥v∥ 1 thus strictly decreases throughout the iterations. In addition, v is always in Ω(T ) by Lemma 2. The algorithm stops when v reaches a bottom element in Ω(T ), i.e., an element of T . Consequently, the algorithm terminates when ∥v∥ 1 = d if T = {1} or when ∥v∥ 1 = 1 if T = {e 1 , . . . , e d }, after less than ω iterations.

In the rest of the paper, we will mainly focus on the following finite sequences:

Definition 1. A sequence of matrices (U n ) 0≤n≤N is valid if and only if U 0 = I d and every U n , n ∈ {1, . . . , N }, is in U d .
In addition, a valid sequence reduces a vector v if and only if there exists

T ⊂ P d such that U n • • • U 0 v ∈ Ω(T ) \ T for all n ∈ {0, . . . , N -1} and U N • • • U 0 v ∈ T . In that case, we set v n := U n • • • U 0 v and a n := U -1 0 • • • U -1 n v N .
By design and by Theorem 3, Algorithm 1 returns a sequence of matrices that reduces the input vector v. Furthermore, the output of the algorithm provides a way of representing the input thanks to the following equalities:

a N = U -1 0 • • • U -1 N v N = U -1 0 • • • U -1 N (U N • • • U 0 v) = v. (1)
Note that a n is generally not equal to v n and, for n ̸ = N , is not equal to v.

Examples in low dimensions

In this subsection, we show how Algorithm 1 relates with the Euclidean algorithm in 2D and some of its well-known generalizations in 3D.

In 2D.

The set U 2 only contains two matrices:

L := 1 -1 0 1
and R := 1 0 -1 1 .

When fed up with the terminal set T = {1} and the input vector v, Algorithm 1 does the following:

if v 1 > v 2 , then v ← Lv and repeat, if v 1 < v 2 , then v ← Rv and repeat,
otherwise, v must be equal to 1, stop and return, which is exactly the Euclidean algorithm in its additive form.

For instance, with T = {1} and t v = (5, 3) as input, Algorithm 1 returns the list (I 2 , L, R, L). Hence the following equivalence: [START_REF] Berthé | Fully Subtractive algorithm, tribonacci numeration and connectedness of discrete planes[END_REF] ). Note that using the terminal set T = {e 1 , e 2 } instead of T = {1} results in replacing 1 with L -1 e 2 .

1 = LRL( 5 3 ) ⇔ L -1 R -1 L -1 1 = ( 5

In 3D.

In dimension 2, it is clear to decide which coordinate has to be subtracted to the other, whereas in dimension 3, it is not the case, hence the diversity of existing generalizations. See [START_REF] Labbé | 3-dimensional continued fraction algorithms cheat sheets[END_REF] for a list as well as a detailed analysis of several of them. Several existing algorithms stick to a convention. For instance, assuming w.l.o.g. 0 ≤ v 1 ≤ v 2 ≤ v 3 , we have:

Brun: v → 1 0 0 0 1 0 0 -1 1 v, Selmer : v → 1 0 0 0 1 0 -1 0 1 v.
In both strategies, the subtraction is encoded in a matrix U ∈ U 3 . Furthermore, for any vector of coprime and strictly positive coordinates, Selmer can be reproduced using Algorithm 1 and T = {1}. In order to be able to deal with null coordinates, it is enough to choose T = {e 1 , e 2 , e 3 }. Brun can also be reproduced using Algorithm 1, but only with T = {e 1 , e 2 , e 3 } because a zero coordinate may appear whenever the two largest coordinates are equal.

For instance, with T = {e 1 , e 2 , e 3 } and t v = (1, 2, 2) as input, using Brun in line 3 for the matrix selection and taking the index order in case of ties, Algorithm 1 returns the terminal e 1 and the following list:

U 0 = I 3 , U 1 = 1 0 0 0 1 -1 0 0 1 , U 2 = 1 0 0 0 1 0 -1 0 1 , U 3 = U 2 .
Hence the following equivalence:

e 1 = U 3 U 2 U 1 U 0 1 2 2 ⇔ U -1 0 U -1 1 U -1 2 U -1 3 e 1 = 1 2 2
.

Other less elementary strategies may be used as well. For instance,

Poincaré: let 0 ≤ v 1 ≤ v 2 ≤ v 3 , v → 1 0 0 -1 1 0 0 -1 1 v, Jacobi-Perron: let 0 < v 1 , v 2 ≤ v 3 and q 2 = [ v 2 v 1 ], q 3 = [ v 3 v 1 ], v → 1 0 0 -q 2 1 0 -q 3 0 1 v.
One can check that we can choose T = {e 1 , e 2 , e 3 } as terminal set for these two strategies by possibly applying the two-dimensional Euclidean algorithm once a coordinate has reached the value 0 in the Jacobi-Perron case. Even if the above matrices are not in U 3 , one can check that they can be replaced by a product of elementary matrices, which are all in U 3 . Algorithm 1 is thus also a generalization of such strategies.

Relation with Digital Geometry

The goal of this subsection is to relate Algorithm 1 to algorithms recently developed in the context of digital geometry [START_REF] Lachaud | Two plane-probing algorithms for the computation of the normal vector to a digital plane[END_REF][START_REF] Lu | A New Lattice-based Plane-probing Algorithm[END_REF]. This relation is based on the transpose of the partial product U n . . . U 0 , i.e., t U 0 • • • t U n . Even if the above-mentioned algorithms are three-dimensional, the following propositions are valid in a d-dimensional space. They gather several properties involving the matrices t U 0 • • • t U n as well as v n and a n introduced in Definition 1: Proposition 4. For all n ∈ {0, . . . , N } and i ∈ {1, . . . , d}:

(i) t v 0 t U 0 • • • t U n e i = t e i v n , (ii) t a n t U 0 • • • t U n e i = t e i v N ,
Proof. In both cases, it is enough to use the definitions of v n and a n :

(i) t v 0 t U 0 • • • t U n e i = t e i U n • • • U 0 v 0 = t e i v n . (ii) t a n t U 0 • • • t U n e i = t e i U n • • • U 0 (U -1 0 • • • U -1 n v N ) = t e i v N .
Definition 2. For all n ∈ {0, . . . , N }, let Λ n := {x ∈ Z d | t a n x = 0} and let B n be defined as follows: Proof. (Sketch) For both definitions of B n , it is enough to notice that the dot product between a n and each vector of the tuple is equal to 0 by Proposition 4, item (ii). The fact that the tuple is a basis of the lattice Λ n comes from the fact that the matrix

B n := t U 0 • • • t U n (e i -e i-1 ) i∈{2,...,d} if v N = 1, B n := t U 0 • • • t U n e i ) i∈{2
t U 0 • • • t U n has determinant 1.
Fig. 2 illustrates the action of the matrix t U 0 • • • t U N . Note that for each n ∈ {0, . . . , N },

the matrix t U 0 • • • t U n transforms the basis {e i | i ∈ {1, 2, 3}} into { t U 0 • • • t U n e i | i ∈ {1, 2, 3}}.
In other words, it deforms an orthonormal basis to a basis that is more and more aligned with the plane of normal v 0 because the sum of the quantities

{ t v 0 t U 0 • • • t U n e i | i ∈ {1, 2,
3}} is smaller and smaller by Proposition 4, item (i). At the last step, the i-th column of t U 0 • • • t U N is perfectly aligned with the plane of normal v 0 when t e i • v N = 0, but not when t e i • v N = 1. Despite this, one can always deduce d -1 vectors perfectly aligned with the plane of normal v 0 by Proposition 5.

The next proposition requires the following definition.

U 0 • • • t U N v 0 Figure 2: Geometrical interpretation of the matrix t U 0 • • • t U N .
We have implicitly represented the scalar projection of the basis vectors in the direction of v 0 with the help of dotted lines. The scalar products

{ t v 0 t U 0 • • • t U N e i | i ∈ {1
, 2, 3}} are equal to the coordinates of v N and are thus all equal to 0 or 1. Definition 3. The digital hyperplane of normal v ∈ P d is the infinite set:

P v := {x ∈ Z d | 0 ≤ t vx < ∥v∥ 1 }.
(2) Proposition 6. For all n ∈ {0, . . . , N } and i ∈ {1, . . . , d}:

(i) if 1 -t U 0 . . . t U n e i ∈ P v 0 , then t e i v n ≥ 1. (ii) if 1 -t U 0 . . . t U n e i / ∈ P v 0 and 1 + t U 0 . . . t U n e i / ∈ P v 0 , then t e i v n = 0.
Proof. We first have

t v 0 (1 ± t U 0 . . . t U n e i ) = t v 0 1 ± t v 0 ( t U 0 . . . t U n e i ) = ∥v 0 ∥ 1 ± t e i v n ,
where the last equality comes from t v 0 1 = ∥v 0 ∥ 1 and Proposition 4, item (i). By Definition 3,

1 -t U 0 . . . t U n e i ∈ P v 0 implies t v 0 (1 -t U 0 . . . t U n e i ) < ∥v 0 ∥ 1 ⇔ ∥v 0 ∥ 1 -t e i v n < ∥v 0 ∥ 1 ⇔ t e i v n > 0,
which proves (i).

We now focus on (ii). By Definition 3, 1 ± t U 0 . . . t U n e i / ∈ P v 0 implies either

t v 0 (1 ± t U 0 . . . t U n e i ) ≥ ∥v 0 ∥ 1 or t v 0 (1 ± t U 0 . . . t U n e i ) < 0.
The second inequality is equivalent to

± t e i v n > ∥v 0 ∥ 1 , which is false because ± t e i v n ≤ ∥v n ∥ 1 ≤ ∥v 0 ∥ 1 .
As a result, only the first inequality, which is equivalent to ± t e i v n ≤ 0, is true. We thus clearly get t e i v n = 0 if we have both

1 -t U 0 . . . t U n e i / ∈ P v 0 and 1 + t U 0 . . . t U n e i / ∈ P v 0 .
Let us now consider the predicate InPlane := "is x in P?" for any hyperplane P ∈ {P v | v ∈ P d } of unknown normal vector. Proposition 6 shows that we are able to check if a given coordinate of the unknown normal vector is strictly positive or equal to 0 by probing the hyperplane with the predicate InPlane. By Lemma 1, such tests on coordinates are a way of selecting a possible matrix on line 3 of Algorithm 1. As a consequence, a geometrical version of Algorithm 1 can be obtained by replacing the vector v by the predicate InPlane and changing the way a new matrix is selected accordingly. The three-dimensional plane-probing algorithms H, R [START_REF] Lachaud | Two plane-probing algorithms for the computation of the normal vector to a digital plane[END_REF] and L [START_REF] Lu | A New Lattice-based Plane-probing Algorithm[END_REF] rely on exactly that idea. The difference between the three algorithms lies in their set of possible operations. Those operations can be represented up to symmetries by a matrix:

v → 1 -a -b 0 1 0 0 0 1 v.
In H, a is set to 1 and b to 0, which means that the matrix belongs to U 3 . In R, a is also set to 1, but b can be any nonnegative integer. In L, both a and b can be any nonnegative integers, provided that they are not both null. The matrices used in R and L are not necessarily in U 3 . However, if they are not, they can be replaced by a product of elementary matrices, which are all in U 3 . We can therefore assume that H, R and L produce a valid sequence of matrices (U) 0≤n≤N . In addition, that sequence reduces the normal vector of the underlying hyperplane because:

∀n ∈ {0, . . . , N -1}, U n • • • U 0 v ∈ Ω({1}) \ {1} by item 1 of [20, Lemma 1], U N • • • U 0 v = 1 by [20, Theorem 2].
As in (1), we have a N = v [20, Corollary 4]. In other words, those algorithms, which compute a N only from the predicate InPlane, are indeed able to retrieve the normal vector v of P.

To end the section, note that one can even further generalize our approach if one replaces InPlane by another set-membership predicate involving a given digital surface as shown in Fig. 1 and [20, Section 5].

Discussion

In this subsection, we compare the plane-probing algorithms H, R and L to classical threedimensional Euclidean algorithms, in particular, Brun, Selmer, Poincaré and Jacobi-Perron presented above.

What makes those two kinds of algorithms very different from each other lies in the way they select a matrix in line 3 of Algorithm 1. On one hand, plane-probing algorithms use a geometrical criterion involving the column vectors of t U 0 • • • t U n . On the other hand, classical three-dimensional Euclidean algorithms, even when run with the predicate InPlane, use an arithmetical criterion involving the coordinates of

U n • • • U 0 v.
Non-uniqueness. In general, several matrices can be chosen. The arithmetical criterion aims at reducing the set to one, but if v n has equal coordinates, there is still more than one possible matrix to get v n+1 . An often-used convention to break ties is to compare the indices of the equal coordinates. The geometrical criterion is based on an in-sphere test [START_REF] Lachaud | Two plane-probing algorithms for the computation of the normal vector to a digital plane[END_REF] and usually reduces the set to one, but in case of co-spherical points, again, there is more than one possible matrix. In such cases, one can resort to a lexicographic order.

Fig. 3 illustrates the first steps when applying Selmer (top row) and H (bottom row) on

t v = (2, 3, 5). If the point 1 -t U 0 • • • t U n (e k -e l )
, with distinct k, l ∈ {1, 2, 3}, belongs to Psuch a point is depicted with a red disk in Fig. 3 -, then there exists a unique matrix U ∈ U 3 containing -1 at the intersection of the k-th row and the l-th column such that

1 -t U 0 • • • t U n (e k -e l ) = 1 -t U 0 • • • t U n t Ue k ∈ P.
As a result, by Proposition 6 and Lemma 1, U is a possible candidate for U n+1 . With that geometrical point of view, in Fig. 3, choosing a red disk is like choosing a matrix. Selmer first selects x (see (a) and (b) in Fig. 3). Then, an equality between two coordinates of v 1 is detected, because there are two points symmetric about 1, namely y and y ′ , that do not belong to P. From the two candidate points

x ′ = 1 -t U 0 • • • t U n e 2 and x ′′ = 1 -t U 0 • • • t U n e 3 ,
the second one is chosen because it has the greatest index: 3 (see (b) and (c) in Fig. 3).

Regarding H, the red and blue disks are in a co-spherical position at the very beginning, but z is chosen according to the lexicographic order on points (see (d) and (e) in Fig. 3). Then, z ′′ is chosen because the sphere passing by the triangle and z ′′ does not contain z ′ (see (e) and (f) in Fig. 3).

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • e 1 e 2 e 3 • • • • • • • • • • x (a) t v0 = (2, 3, 5) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • x x ′ x ′′ y ′ y (b) t v1 = (2, 3, 3) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • x ′′ (c) t v2 = (2, 3, 1) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • e 1 e 2 e 3 • • • • • • • • • • z (d) t v0 = (2, 3, 5) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • z z ′ z ′′ (e) t v1 = (2, 3, 2) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • z ′′ (f) t v2 = (2, 1, 2)
-t U 0 • • • t U n e i ) i∈{1,2,3}
is depicted as a blue triangle, while the set of differences (1

-t U 0 • • • t U n (e k -e l )) k,l∈{1,2,3}
,k̸ =l is depicted as a red hexagon, which is a usual representation for plane-probing algorithms. The vertices are depicted with disks (resp. circles) if they belong (resp. do not belong) to the digital plane.

Impact of previous choices. In the arithmetical approach, the choice of U n+1 only depends on the coordinates of v n , regardless of how v n has been obtained. If, reducing two different input vectors v and v ′ , one has v n = v ′ n ′ , the next matrices will be exactly the same for the two sequences. This is not the case in the geometrical approach, because the choice of U n+1 depends on the whole sequence. In that approach, even if v n equals v ′ n ′ , the next matrices will be in general different for the two reductions.

Geometrical properties. The column vectors of t U 0 • • • t U n exhibit interesting properties in the geometrical approach. One such property is that the basis of the rank-two lattice Λ n returned by R (resp. H) is experimentally always (resp. almost always) reduced [START_REF] Lachaud | Two plane-probing algorithms for the computation of the normal vector to a digital plane[END_REF]. It has even been proven that the basis returned by L is always reduced [START_REF] Roussillon | Delaunay property and proximity results of the L-algorithm[END_REF]. Another property concerns the points selected during the execution of the algorithm, i.e., the set:

{1 -t U 0 . . . t U n e i | n ∈ {0, . . . , N }, i ∈ {1, . . . , d}}. (3)
Experimentally, those points are in average much closer to 1 in the geometrical approach than in the arithmetical approach. This observation is illustrated in Fig. 4 and supported by the experiments described below.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (a) Selmer • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (b) H
- t U 0 • • • t U n e i ) i∈{1,2,3}
is depicted as a blue triangle.

Experimental results.

As in [20, Table 1], we ran all algorithms on vectors ranging from (1, 1, 1) to (200,200,200). There are 6578833 vectors with relatively prime components in that range. For all of them, we checked whether the returned basis is reduced or not. In the latter case, we also computed the number of reductions needed to reduce it a posteriori. Results are reported in Table 1. The returned bases are generally not reduced with the arithmetical approach, while they generally are with the geometrical approach. Furthermore, we computed for all vectors the maximal distance between 1 and every point of the set described in (3). In Fig. 5, we plotted that maximal distance against the magnitude of the input vector. In order to smooth the lines, we discretized the magnitude of the input vector into bins of size 10 and took the average for every bin. other because those three algorithms often make the same choices for the considered vectors. More importantly, those lines lie quite below the other ones, which means that the geometrical approach is much more local than the arithmetical one, as suggested by Fig. 4. The maximal distance between 1 and every point selected during the execution of the algorithm is plotted against the magnitude of the input vector (see the text for the experimental settings). The line labeled by P is truncated, because of its extremely large slope.

Pattern Generation with Extensions of Substitutions

As discussed in the previous section, Algorithm 1 provides a terminal and a sequence of matrices that reduces the input vector v. A variant consists in passing in the predicate InPlane instead of v as with plane-probing algorithms. In this section, we show how to generate a piece of digital hyperplane of normal v from a sequence of matrices returned by Algorithm 1 or its variant.

Our method is based on a description of digital hyperplanes as unions of faces, recalled in Subsection 3.1. The mathematical framework we use is presented in Subsection 3.2. It provides us with a clear definition from which one can derive a generation method in Subsection 3.3. The properties of the generated sets are then discussed in the last two subsections.

Digital Hyperplanes as Unions of Faces

For x ∈ Z d and i ∈ {1, . . . , d}, we define the (d -1)-dimensional face in x perpendicular to the direction i as the following subset of R 3 :

(x, i * ) := x + j∈{1,...,d}\i λ j e j | λ j ∈ [0, 1] . (4)
The point x will be called the origin of the face and i the type of the face. A three-dimensional illustration is given in Fig. 6 (a).

S v := {(x, i * ) ∈ Z d × {1, . . . , d} | 0 ≤ t vx < t ve i }. (5)
The link between digital hyperplanes (Definition 3) and stepped hyperplanes (Definition 4) is made explicit in the following proposition.

Note that the vertices of a face (x, i

* ) ∈ S v form the set {x + t | t ∈ T i }.
We first show that all the vertices of a face (x, i * ) ∈ S v are in P v . By definition (x, i * ) ∈ S v is equivalent to 0 ≤ t vx < t ve i . By adding t vt, t ∈ T i , we get t vt ≤ t v(x + t) < t v(e i + t). In addition, we have by definition 0 ≤ t vt and t ve i + t vt < t ve i + t v j∈{1,...,d}\i e j = t v1 = ∥v∥ 1 . We have therefore 0 ≤ t v(x + t) < ∥v∥ 1 , i.e., x + t ∈ P v for all t ∈ T i by [START_REF] Arnoux | Pisot substitutions and Rauzy fractals[END_REF]. Now, we show that for all points y ∈ P v , there is a face of type i in S v such that y is one of its vertices. More precisely, we choose i such that t ve i is maximal, i.e., i = arg max j∈{1,...,d} t ve j . In addition, we number the K partial sums of T i according to the following order:

t vt 1 ≤ • • • ≤ t vt k ≤ • • • ≤ t vt K .
Note that t K = j∈{1,...,d}\i e j and, because of our choice for i, t vt k+1 -t vt k ≤ t ve i for all k ∈ {1, . . . , K -1}.

If t vt K ≤ t vy < ∥v∥ 1 , then 0 ≤ t v(yt K ) < ∥v∥ 1 -t vt K = t ve i , i.e. (yt K , i * ) ∈ S v by [START_REF] Berthé | Critical connectedness of thin arithmetical discrete planes[END_REF].

For all k ∈ {1, . . . , K -1}, if t vt k ≤ t vy < t vt k+1 , then 0 ≤ t v(y-t k ) < t vt k+1 -t vt k ≤ t ve i , i.e. (yt k , i * ) ∈ S v by (5).

Words, Substitutions and their Extensions

In this subsection, we present the framework we use for our generation method. It is based on geometrical extensions of substitutions: rules that replace faces by unions of faces in the same way that substitutions replace letters by words.

We consider a d-letter alphabet A d := {1, . . . , d}. A word is an element of the free monoid A ⋆ d generated by A d . The empty word is denoted by ϵ and the concatenation operation is denoted by • or is left implicit. The abelianization mapping f :

A ⋆ d → N d is such that t f (w) = (|w| 1 , |w| 2 , |w| 3 )
, where |w| i denotes the number of occurrences of the letter i in w.

Example 1. The word w = 1 • 2 • 3 • 1 • 2 = 12312 is an element of {1, 2, 3} ⋆ . t f (w) = (2, 2, 1).
A substitution σ over A d is a non-erasing endomorphism of A ⋆ d , completely defined by its image on the letters of A d by the relation σ(w

• w ′ ) = σ(w) • σ(w ′ ). The incidence matrix M σ of σ is the d × d matrix whose i-th column is f (σ(i)) for every i ∈ A d . For any word w ∈ A ⋆ d , applying σ then f is the same as applying f then M σ , i.e., M σ f (w) = f (σ(w)). Example 2. For σ: 1 → 1, 2 → 21, 3 → 32, σ(12312) = 1 • 21 • 32 • 1 • 21 = 12132121. M σ = 1 1 0 0 1 1 0 0 1 , M σ f (w) = f (σ(w)) = 4 3 1 .
For two substitutions σ ′ , σ ′′ , we denote by • their composition: (σ ′′ • σ ′ )(w) := σ ′′ (σ ′ (w)). Note that the above commutation relation implies that

M σ ′′ •σ ′ = M σ ′′ M σ ′ .
A substitution is said to be unimodular if det (M σ ) = ±1. In that case, M σ is invertible and there is a formula involving M -1 σ for the (d -1)-extension of σ [1, Proposition 2.11]. From now on, we assume that all substitutions are unimodular. The (d -1)-extension of σ acting on (d -1)-dimensional faces is then defined as follows:

E d-1 (σ)(x, i * ) := s|σ(j)=p•i•s j∈A d (M -1 σ (x + f (s)), j * ). ( 6
)
The subscript means that the union is done over all pairs (j, s) such that j ∈ A d and σ(j) = p•i•s, where i is a letter, whereas p and s are words.

Example 3. Let us consider the substitution σ: 1 → 1, 2 → 12, 3 → 1213 over A 3 and let us consider the image by E 2 (σ) of (0, 1 * ). Since the type of that face is 1, the union in (6) must be done over all pairs (j, s) such that j ∈ A 3 and σ(j) = p • 1 • s. There are four such pairs:

(1, ϵ), (2, 2), (3, 213), (3, 3)
. Furthermore, since the origin of the face is 0, every pair (j, s) gives the face (M -1 σ (f (s)), j * ). Let us consider the third pair as an example: (3, 213) corresponds to the face (M -1 σ (f (213)), 3 * ), where f (213) = e 1 + e 2 + e 3 . The image of (0, 1 * ) thus consists of the following four faces:

(0, 1 * ), (M -1 σ e 2 , 2 * ), (M -1 σ (e 1 + e 2 + e 3 ), 3 * ), (M -1 σ (e 3 ), 3 * ).
Geometrically, the action of E d-1 (σ) can be considered as a digitization of the action of M -1 σ (see Section 2 and Fig. 3 in [START_REF] Fernique | Generation and recognition of digital planes using multi-dimensional continued fractions[END_REF]). Example 3 is illustrated in Fig. 7 with that point of view. We extend [START_REF] Berthé | A study of Jacobi-Perron boundary words for the generation of discrete planes[END_REF] to unions of faces:

E d-1 (σ)(F) := (x,i * )∈F E d-1 (σ)(x, i * ). (7)
It is convenient to use the following notation for translations of faces: if (x, i * ) is a face and y is a vector, then (x, i * ) + y := (x + y, i * ), which extends in a natural way to union of faces. One can then check that σ and E 2 (σ), where σ is defined as in Example 3. The faces of type 1, 2, 3 are respectively colored with red, green and blue. because M -1 σ x does not depend on the union in [START_REF] Berthé | A study of Jacobi-Perron boundary words for the generation of discrete planes[END_REF]. The images by E d-1 (σ) of the faces (0, i * ) for 1 ≤ i ≤ d thus suffice to define E d-1 (σ).

E d-1 (σ)(x, i * ) = M -1 σ x + E d-1 (σ)(0, i * ), (8)
For two substitutions σ ′ , σ ′′ , we also denote by • the composition of their extensions:

E d-1 (σ ′′ )• E d-1 (σ ′ ) (F) := E d-1 (σ ′′ ) E d-1 (σ ′ )(F) .
We have the following key property [15, Proposition 1.2.4, item [START_REF] Arnoux | Algebraic numbers, free group automorphisms and substitutions on the plane[END_REF]]:

E d-1 (σ ′′ ) • E d-1 (σ ′ ) = E d-1 (σ ′ • σ ′′ ). (9)
From ( 6), ( 7), ( 8) and ( 9), one can derive the following proposition, which shows that one can incrementally compute the image by E d-1 (σ ′ • σ ′′ ) of a face (0, i * ) as a union-translation scheme:

Proposition 8. We have for all i ∈ A d , E d-1 (σ ′ • σ ′′ )(0, i * ) = s|σ ′ (j)=pis j∈A d (M σ ′ •σ ′′ ) -1 f (s) + E d-1 (σ ′′ )(0, j * ) .
Proof.

E d-1 (σ ′ • σ ′′ )(0, i * ) = E d-1 (σ ′′ ) E d-1 (σ ′ )(0, i * ) = E d-1 (σ ′′ ) σ ′ (j)=pis (M -1 σ ′ f (s), j * ) = σ ′ (j)=pis E d-1 (σ ′′ ) M -1 σ ′ f (s), j * ) = σ ′ (j)=pis (M σ ′′ ) -1 (M σ ′ ) -1 f (s) + E d-1 (σ ′′ )(0, j * ) = σ ′ (j)=pis (M σ ′ •σ ′′ ) -1 f (s) + E d-1 (σ ′′ )(0, j * ) .
Example 4. Let us consider the following substitutions:

σ ′ : 1 → 1 2 → 2 3 → 31 , σ ′′ : 1 → 1 2 → 21 3 → 321 , σ ′ • σ ′′ : 1 → 1 2 → 21 3 → 3121 .
As shown in Fig. 8, σ ′ describes how the images by E 2 (σ ′ • σ ′′ ) relate to the images by E 2 (σ ′′ ).

First, let us consider E 2 (σ ′ • σ ′′ )(0, 1 * ), which is red in (b). It has been obtained as the union of two previous sets: E 2 (σ ′′ )(0, 1 * ), which is red in (a), and E 2 (σ ′′ )(0, 3 * ), which is blue in (a), because letter 1 belongs to both σ ′ (1) and σ ′ (3) and is also in the last position, which means with no suffixe and thus no translation. Let us now consider

E 2 (σ ′ • σ ′′ )(0, 3 * ), which is blue in (b).
It is a copy of E 2 (σ ′′ )(0, 3 * ), which is blue in (a), translated by M σ ′ •σ ′′ e 1 = e 1 , because letter 3 appears only once in σ ′ (3), followed by the single-letter word 1.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (a) E 2 (σ ′′ ) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (b) E 2 (σ ′ • σ ′′ )
Figure 8: The images of (0, 1 * ), (0, 2 * ) and (0, 3 * ) by E 2 (σ ′′ ) (a) and E 2 (σ ′ •σ ′′ ) (b) are displayed respectively in red, green and blue (see Example 4).

Sequences of Substitutions and Patterns

A sequence of matrices (U n ) 0≤n≤N that reduces a vector v is introduced in Definition 1. In this subsection, we define its associated sequence of substitutions.

Definition 5 (Associated substitutions). Let (U n ) 0≤n≤N be a sequence of matrices that reduces a vector v ∈ P d . Its associated sequence of substitutions (σ n ) 0≤n≤N verifies, for all n ∈ {0, . . . , N }, M σn = ( t U n ) -1 and for all i ∈ A d , σ n (i) starts with letter i.

Example 5. Let us consider the following matrices:

U 1 = 1 0 0 0 1 -1 0 0 1 , ( t U 1 ) -1 = 1 0 0 0 1 0 0 1 1 .
Remark that ( t U 1 ) -1 is the incidence matrix of two disinct substitutions:

σ ′ : 1 → 1, 2 → 23, 3 → 3, σ ′′ : 1 → 1, 2 → 32, 3 → 3.
Actually, the inverse of the transpose of all matrices in U d is the incidence matrix of two distinct substitutions. According to Definition 5, σ 1 := σ ′ because σ ′ (i) starts with letter i for all i ∈ A d , which is not the case for σ ′′ .

To save space, we set

σ n•••0 := σ n • • • • • σ 0 for 1 ≤ n ≤ N . In accordance with Definition 1, we have thus: ∀n ∈ {0, . . . , N }, v n = ( t M σ n•••0 ) -1 v, a n = t M σ n•••0 v N . (10)
Table 2 shows, for a given sequence, how

σ n , M n , σ n•••0 , M σ n•••0 , v n
and a n are related.

We will now focus on the image by E d-1 (σ n•••0 ) of faces with origin at 0. The result is a finite set of faces that we call pattern. A full example is given by Table 2 and Fig. 9.

Definition 6 (Pattern). Let (U n ) 0≤n≤N be a sequence of matrices that reduces a vector v ∈ P d and let (σ n ) 0≤n≤N be its associated sequence of substitutions. Let W 0 be the faces of type i and origin 0 such that t v N e i = 1, i.e., ∪ i∈A d , t v N e i =1 (0, i * ). For all n ∈ {1, . . . , N }, let W n be the image of

W 0 by E d-1 (σ n•••0 ), i.e., W n := E d-1 (σ n•••0 )(W 0 ). n 0 1 2 3 σ n 1 → 1 2 → 2 3 → 3 1 → 1 2 → 2 3 → 32 1 → 1 2 → 21 3 → 3 1 → 1 2 → 2 3 → 31 M n I 3 1 0 0 0 1 1 0 0 1 1 1 0 0 1 0 0 0 1 1 0 1 0 1 0 0 0 1 σ n•••0 1 → 1 2 → 2 3 → 3 1 → 1 2 → 2 3 → 32 1 → 1 2 → 21 3 → 321 1 → 1 2 → 21 3 → 3121 M σ n•••0 I 3 1 0 0 0 1 1 0 0 1 1 1 1 0 1 1 0 0 1 1 1 2 0 1 1 0 0 1 t a n 1 ( 1 1 2 ) ( 1 2 3 ) ( 1 2 4 ) t v n ( 1 2 4 ) ( 1 2 2 ) ( 1 1 2 ) 1
Table 2: A sequence of substitutions and other related quantities. 2. It is associated with a sequence of matrices obtained from the reduction of (1, 2, 4) using algorithm H [START_REF] Lachaud | Two plane-probing algorithms for the computation of the normal vector to a digital plane[END_REF]. In the top row, the images of (0, 1 * ), (0, 2 * ) and (0, 3 * ) by E 2 (σ n•••0 ) are displayed respectively in red, green and blue. In the bottom row, the faces of type 1, 2, 3 are respectively colored with red, green and blue.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (a) W0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (b) W1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (c) W2 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (d) W3 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (e) W0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (f) W1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (g) W2 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • (h) W3
From a sequence of substitutions (σ n ) 0≤n≤N , there are at least two ways of generating patterns. In the first method, we compute the substitution σ n•••0 and then apply the formula given by [START_REF] Berthé | A study of Jacobi-Perron boundary words for the generation of discrete planes[END_REF] to get W n . The second method is incremental: we compute W n from W n-1 as in a union-translation scheme thanks to Proposition 8. One can check that the time-complexity of both methods linearly depends on the number of faces in W n . By Proposition 10 or Theorem 11, item (iii), proven in the next subsection, this number is equal to ∥a n ∥ 1 . As a result, the overall time-complexity of both methods is O(∥a n ∥ 1 ).

From a vector v, one has to first compute a sequence of matrices that reduces v and then compute the associated sequence of substitutions. The overall time-complexity thus depends on the reduction algorithm. For most of them, e.g., Brun, Selmer, H [START_REF] Lachaud | Two plane-probing algorithms for the computation of the normal vector to a digital plane[END_REF], every matrix choice takes a constant time, while there are at most O(∥v∥ 1 ) choices in total (Theorem 3). In that case, the reduction stage takes O(∥v∥ 1 ) time. Then, the generation of W N also takes O(∥v∥ 1 ) time, because a N = v by (1). However, note that there are reduction algorithms, e.g., R [START_REF] Lachaud | Two plane-probing algorithms for the computation of the normal vector to a digital plane[END_REF], L [START_REF] Lu | A New Lattice-based Plane-probing Algorithm[END_REF], which does not compute new matrices in constant time, but in logarithmic time, thus increasing the whole time-complexity by a logarithmic factor.

Properties of Patterns

In this subsection, we show several properties for the patterns defined in Definition 6. Even if a pattern is a set of faces, we will say that a point x is in a pattern W, denoted by x ∈ W by abuse of notation, whenever there exists a face w ∈ W such that x is a vertex of w.

The following proposition shows that the patterns always contain some very specific faces and points.

Proposition 9. Let us assume that v N ∈ {1} ∪ {e 1 , . . . , e d } and let A v N be equal to {i ∈ A d | t v N e i = 1}. For all n ∈ {0, . . . , N } and for all i ∈ A v N , (e i -M -1 σ n•••0 e i , i * ) ∈ W n and 1 -M -1 σ n•••0 e i ∈ W n .
Proof. Since every word σ 0 (i), . . . , σ n (i) starts with i, σ n•••0 (i) also starts with i. Let s be the suffix such that

σ n•••0 (i) = i • s. The face (M -1 σ n•••0 f (s), i * ) is by definition included in E d (σ n•••0 )(0, i * ) and thus in W n for all i ∈ A v N . In that case, f (s) = f (σ n•••0 (i)) -e i , where f (σ n•••0 (i)) is the i-th column of M σ n•••0 and is thus equal to M σ n•••0 e i .
We have therefore:

W n ∋ (M -1 σ n•••0 f (s), i * ) = (M -1 σ n•••0 f (σ n•••0 (i)) -e i , i * ) = (M -1 σ n•••0 M σ n•••0 e i -e i , i * ) = (e i -M -1 σ n•••0 e i , i * ).
In addition, the point 4) and thus is in W n .

e i + k̸ =i e k -M -1 σ n•••0 e i = 1-M -1 σ n•••0 e i is a vertex of (e i -M -1 σ n•••0 e i , i * ) by (
In plane-probing algorithms, the point used to select a new matrix at a step n is exactly Proposition 6 and (3)) and thus lies in W n by Proposition 9. The convention used in Definition 5 in order to associate only one substitution to a given matrix turns out to guarantee that the generated patterns contain those points.

1 -M -1 σ n•••0 e i (see
The following proposition shows that all faces in W n have a specific location in space with respect to direction a n . It is key for several other properties stated in Theorem 11.

Proposition 10. Let us assume that v N ∈ {1} ∪ {e 1 , . . . , e d }. For all n ∈ {0, . . . , N } and for all j ∈ A d , for every face (x, j * ) ∈ W n , the quantity t a n x appears once and only once in the integer interval {0, . . . , t a n e j -1}.

Proof. By Definition 6 and according to [START_REF] Berthé | A study of Jacobi-Perron boundary words for the generation of discrete planes[END_REF], all faces of type j in W n are located along a path encoded by σ n•••0 (j).

If there is only one letter i in σ n•••0 (j), there is only one face in E d-1 (σ n•••0 )(0, i * ) whose origin is 0 and the proposition is obviously true.

Otherwise, let us denote by a and b two consecutive letters in σ n•••0 (j). In other words

σ n•••0 (j) = p • a • b • s. Letter a (resp. b) corresponds to a face in E d-1 (σ n•••0 )(0, a * ) (resp. E d-1 (σ n•••0 )(0, b * ))
. By abuse of terminology, we call the faces a and b and say that they are consecutive. The difference between the origins is equal to M

-1 σ n•••0 f (bs)-M -1 σ n•••0 f (s) = M -1 σ n•••0 e b . However, the quantity t a n M -1 σ n•••0 e b is equal to t v N M σ n•••0 M -1 σ n•••0 e b = t v N e b .
Remind that if t v N e b = 1, the face of (M -1

σ n•••0 e b , j * ) is kept in W n . Otherwise, t v N e b =
0 and the face is discarded. As a consequence, the dot product between a n and the difference between the origins of a and b is 1 if b is in W n and 0 otherwise. The set { t a n x | (x, j * ) ∈ W n } thus form an integer interval whose lower bound is 0. The upper bound plus one is given by the total number of faces of type j in W n , which is equal to the number of letters

i ∈ σ n•••0 (j) such that t v N e i = 1: t v N f (σ n•••0 (j)) = t v N M σ n•••0 e j = t a n e j .
We are now able to gather several properties of patterns into our main theorem: Two extra properties hold if v N = 1:

(iv) ∀n ∈ {1, . . . , N }, W n-1 ⊂ W n , (v) ∀n ∈ {0, . . . , N }, W n contains points from which one can compute the vectors of a basis of Λ n .

Proof. (i) By Proposition 10, if W n contains at least one face of type j ∈ A d , then, for every face (x, j * ) ∈ W n , we have 0 ≤ t a n x < t a n e j , which means that (x, j * ) is in the stepped plane of normal a n according to [START_REF] Berthé | Critical connectedness of thin arithmetical discrete planes[END_REF].

In addition, if n < N , note that the coordinates of a n are always nonnegative and the multiplication by t M σn consists in adding a coordinate to another, so that a n ≤ a n+1 . By induction, we get a n ≤ a N and in particular t a n e j ≤ t a N e j . As a consequence, all faces of W n are also in the stepped plane of normal a N .

(ii) By Proposition 10, if W n contains at least one face of type j ∈ A d , then, for every face (x, j * ) ∈ W n , the quantity t a n x appears once and only once in the integer interval {0, . . . , t a n e j -1}. Since the copies of W n are translated along the lattice Λ n , i.e., by a translation vector whose dot product with a n is zero, the union of all copies is the set {(x, j * ) | 0 ≤ t a n x < t a n e j }. Taking into account all face types, we get the stepped plane of normal a n according to [START_REF] Berthé | Critical connectedness of thin arithmetical discrete planes[END_REF].

It remains to prove that two distinct copies of W n cannot share a common face (x, i * ). Let us denote by W ′ and W ′′ the two copies and by t ̸ = 0 the translation vector such that W ′′ = W ′ +t. If (x, i * ) ∈ W ′′ , then (xt, i * ) ∈ W ′ . By definition of Λ n , we have t a n t = 0, which implies that t a n x = t a n (xt). However, by Proposition 10, there are no two faces in W n whose origin has the same dot product with t a n . Since that fact is obviously translation invariant, the same is true for W ′ , which means that (x, i * ) and (xt, i * ) cannot be both in W ′ . Consequently, (x, i * ) cannot belong to both W ′ and W ′′ .

(iii) Direct consequence of Proposition 10.

(iv) if v N = 1, we have for all n ∈ {1, . . . , N }:

W n = ∪ i∈A d E d-1 (σ n•••0 )(0, i * ) (Definition 6) = ∪ i∈A d ∪ σn(j)=pis (M σ n•0 ) -1 f (s) + E d-1 (σ n-1•••0 )(0, j * ) (Proposition 8) ⊃ ∪ j∈A d E d-1 (σ n-1•••0 )(0, j * ) = W n-1 (Definition 6),
where the inclusion comes from the trivial fact that for all j ∈ A d , σ n (j) may always be written p • i • ϵ, where i is any letter in A d , i.e., the word σ n (j) ends with a letter i ∈ A d . (v) By Proposition 5, the vectors of a basis of Λ n are M -1 σ n•••0 (e ie i-1 ) i∈{2,...,d} .

In addition, by Proposition 9, W n contains the points 1 -M -1 σ n•••0 e i i∈{1,...,d} .

Consequently, for all i ∈ {2, . . . , d}, we have:

1 -M -1 σ n•••0 e i-1 -(1 -M -1 σ n•••0 e i ) = M -1 σ n•••0 (e i -e i-1
), which concludes the proof.

Theorem 11 shows that our method provides a sequence of patterns, all included in a given stepped plane (i). The last pattern periodically generates the underlying stepped plane without hole or overlap (ii) and is of minimal size (iii) because if one sum the normal vector of all its faces, we get exactly a N , i.e., the normal of the stepped plane, and one cannot expect to find a smaller pattern with the same normal.

In addition, if v N = 1, the patterns form a hierarchical set by inclusion (iv) and they contain points from which one can compute the period vectors (v).

Geometrical and Topological Remarks

In this subsection, we discuss two additional and desirable properties: shape compactness and connectivity.

Our first remark is that the choice of the algorithm has a great impact on the shape of the pattern (see Fig. 10).

Experiments were conducted to measure how far the patterns are from the origin. We generated patterns using the reduction algorithms H, Brun (B), Selmer (S) and Jacobi-Perron (JP) on 811214 vectors with relatively prime components ranging from (1, 1, 1) to (100, 100, 100). We discarded R and L, which behave almost exactly like H in that range, as well as Poincaré, which probes for points very far away from the origin as shown in Fig. 5. For all input vectors, we computed the maximal distance between every point of the pattern and the origin. In Fig. 11, we plotted that maximal distance against the magnitude of the input vector. We smoothed the lines as in Fig. 5. The patterns generated from the geometrical approach (H) are on average much closer to the origin than those generated from the arithmetical approach (B, S and JP).

In our opinion, this is because the basis of the lattice Λ n returned by the algorithm R and L (resp. H) is always (resp. almost always) reduced [START_REF] Lachaud | Two plane-probing algorithms for the computation of the normal vector to a digital plane[END_REF][START_REF] Roussillon | Delaunay property and proximity results of the L-algorithm[END_REF]. There are no such results with classical three-dimensional Euclidean algorithms, because they are not geometry-aware. A perspective is to bound the distance of the pattern boundary to the origin. For the algorithms For all vectors with relatively prime components ranging from (1, 1, 1) to (100, 100, 100), the maximal distance between the origin and every point of the pattern is plotted against the magnitude of the vector. method has a greater practical interest in the context of digital surface analysis.

Conclusion

We have introduced a d-dimensional generalization of the Euclidean algorithm, which covers a lot of algorithms: not only Euclid's algorithm itself and the classical three-dimensional Euclidean algorithms like Brun, Selmer, Poincaré or Jacobi-Perron, but also plane-probing algorithms recently developed in the context of digital geometry. All the classical three-dimensional Euclidean algorithms can even be translated into a plane-probing version due to Proposition 6 and the discussion that follows.

All those algorithms share in common that they return a sequence of elementary matrices that transforms a given normal vector into a terminal element. That output can be used together with a geometrical extension of substitutions, to generate sets of (d -1)-dimensional faces.

The generated sets have mainly three combinatorial properties listed in Theorem 11: they can form a hierarchical set of patterns, all included in a given digital hyperplane. The pattern of the highest level is of minimal size and periodically generates the underlying digital hyperplane without hole or overlap, hence the name pattern. In addition, those patterns experimentally have an additional geometrical property when used together with a plane-probing algorithm: they are rather compact and close to the origin as illustrated in Fig. 10. It is however a perspective to quantify compactness and proximity and to provide upper bounds on such measures. To the best of our knowledge, no other method, whether it was based on digitization or based on union-translation schemes, has all the above properties.

Plane-probing algorithms compute a normal vector by iteratively updating a triangle that locally fits a digital plane or, with slight changes, a digital surface. If the digital surface is not convex, the triangles may jump over holes or stab the digital surface, because the set of probes, which is too sparse, does not faithfully represent the digital surface. In order to correctly approximate the tangent plane in such cases, our idea is to incrementally compute a pattern instead of a triangle: new probes provide the next matrix in the sequence, we then use the proposed method to generate the next pattern and we finally check whether it fits the digital surface or not, as illustrated in Fig. 1. We plan to thoroughly study how accurately do such patterns approximate tangent planes. We hope it will improve the local analysis of digital surfaces using plane-probing.

Figure 3 :

 3 Figure 3: The digital plane of normal t v = (2, 3, 5) is depicted in gray. The vector v has been reduced with Selmer in (a)-(c) and H in (d)-(f). Only the steps n = 0, 1, 2 are shown. For every step, the three-point tuple (1t U 0 • • • t U n e i ) i∈{1,2,3}is depicted as a blue triangle, while the set of differences (1t U 0 • • • t U n (e ke l )) k,l∈{1,2,3},k̸ =l is depicted as a red hexagon, which is a usual representation for plane-probing algorithms. The vertices are depicted with disks (resp. circles) if they belong (resp. do not belong) to the digital plane.

Figure 4 :

 4 Figure 4: The digital plane of normal t v = (1, 13, 17) is depicted in gray. The vector v has been reduced with Selmer in (a) and H in (b). For every step n, the three-point tuple (1 -

  the unit vector and the selected points

Figure 5 :

 5 Figure5: The maximal distance between 1 and every point selected during the execution of the algorithm is plotted against the magnitude of the input vector (see the text for the experimental settings). The line labeled by P is truncated, because of its extremely large slope.

Figure 6 :

 6 Figure 6: In (a), illustration of two-dimensional faces embedded in a three-dimensional space.The faces (0, 1 * ), (0, 2 * ), (0, 3 * ) are respectively in red, green and blue. In (b), the stepped plane of normal t (2, 6, 15) is displayed.

Proposition 7 .

 7 Every face in S v have its vertices in P v . Conversely, every point in P v is a vertex of a face in S v .Proof. Let us introduce the following set of partial sums: ∀i ∈ {1, . . . , d}, T i := j∈J e j | J ⊆ {1, . . . , d} \ i .

Figure 7 :

 7 Figure 7: Image of (0, 1 * ) by M -1σ and E 2 (σ), where σ is defined as in Example 3. The faces of type 1, 2, 3 are respectively colored with red, green and blue.

Figure 9 :

 9 Figure 9: The patterns are the images by E 2 (σ n•••0 ) of the faces with origin at 0. The sequence of substitutions is given in Table2. It is associated with a sequence of matrices obtained from the reduction of (1, 2, 4) using algorithm H[START_REF] Lachaud | Two plane-probing algorithms for the computation of the normal vector to a digital plane[END_REF]. In the top row, the images of (0, 1 * ), (0, 2 * ) and (0, 3 * ) by E 2 (σ n•••0 ) are displayed respectively in red, green and blue. In the bottom row, the faces of type 1, 2, 3 are respectively colored with red, green and blue.

Theorem 11 .

 11 Let us assume that v N ∈ {1} ∪ {e 1 , . . . , e d }. The three following properties hold on the patterns defined in Definition 6: (i) ∀n ∈ {0, . . . , N }, W n ⊂ S an and W n ⊂ S a N . (ii) ∀n ∈ {0, . . . , N }, the copies of W n translated along the lattice Λ n (Definition 2 and Proposition 5) is a partition of S an .(iii) ∀n ∈ {0, . . . , N }, ∀j ∈ A d , W n has t a n e j faces of type j.

  the unit vector and the pattern points

Figure 11 :

 11 Figure 11:For all vectors with relatively prime components ranging from (1, 1, 1) to (100, 100, 100), the maximal distance between the origin and every point of the pattern is plotted against the magnitude of the vector.

  ,...,d} if v N ∈ {e 1 , . . . , e d } and assuming w.l.o.g. that v N = e 1 . The (d -1) tuple B n is a basis of the lattice Λ n .

	Proposition 5.

Table 1 :

 1 The lines for H, R and L lie upon each All algorithms were run on all vectors ranging from (1, 1, 1) to(200, 200, 200) with relatively prime coordinates. Besides H, R and L, other algorithms are denoted by their initials. The average number of reductions is computed only among non-reduced bases.

	Algorithm	H R L	B	S	P	JP
	non-reduced (%) 0.01 0 0 93.37 84.16	96.18 97.29
	avg. reductions	1.00	5.27	5.43 481.68	7.73

* This work has been funded by PARADIS ANR-18-CE23-0007-01 research grant.

H, R and L, such bound may involve geometrical arguments based on the empty-circumsphere criterion.

Our second remark is about the connectivity of the patterns. In general, the patterns may be neither vertex-nor edge-connected. For instance, with Selmer, an example of non-connected pattern is given in Fig. 10 (c). On the contrary, with Brun and Jacobi-Perron, the patterns are always edge-connected as proven in [START_REF] Berthé | A study of Jacobi-Perron boundary words for the generation of discrete planes[END_REF][START_REF] Jolivet | Combinatorics of Pisot Substitutions[END_REF]. Note that the tie-break rule mentioned in Subsection 2.5 is important to guarantee connectivity.

For all vectors with relatively prime components ranging from (1, 1, 1) to (100, 100, 100), we checked whether the pattern is edge-connected or not. Results are reported in Table 3.

Algorithm

H S B JP Edge-connected (%) 73.11 43.04 100 100 Table 3: Patterns were generated from the reduction algorithms H, Selmer (S), Brun (B) and Jacobi-Perron (JP) on all vectors ranging from (1, 1, 1) to (100, 100, 100) with relatively prime coordinates.

Unfortunately, there is an opposition between proximity and connectivity. Indeed, Brun and Jacobi-Perron provide patterns that are edge-connected, but also much less close to the origin than the patterns obtained from the geometrical approach, which may be not connected however. It seems that one cannot have all geometrical, topological and combinatorial properties at the same time, but that statement has yet to be rigorously proven.

The lack of guarantee of connectivity is clearly a drawback. It seems indeed natural to require that patterns must be connected. However, even unconnected patterns may be good candidate to locally approximate tangent planes in the context of digital surface analysis. We tend to think that shape compactness is a more determinant property than connectivity for an accurate and convergent approximation of tangent plane but this has to be further studied.

To end the section, let us mention the method developed in [START_REF] Jamet | Generation of digital planes using generalized continued-fractions algorithms[END_REF], which generates, from an input vector v, a set of points included in a digital plane of normal v. It is incremental and based on several three-dimensional continued-fraction algorithms. It generates point sets that are provably connected but that may have holes. In addition, those sets are not thick enough, e.g., they do not contain any point x such that t vx = ∥v∥ 1 -1, and are composed of a great amount of points that cover the digital plane with overlaps. That is why we believe that our