Isabelle Sivignon
email: isabelle.sivignon@gipsa-lab.grenoble-inp.fr

Fast recognition of a Digital Straight Line Subsegment: two algorithms of logarithmic time complexity

Keywords: Digital geometry, Digital straight segment recognition, subsegment, local convex hull, Farey fan

Given a Digital Straight Line (DSL) of known characteristics (a, b, µ), we address the problem of computing the characteristics of any of its subsegments. We propose two new algorithms that use the fact that a digital straight segment (DSS) can be defined by its set of separating lines. The representation of this set in the Z 2 space leads to a first algorithm of logarithmic time complexity. This algorithm precises and extends existing results for DSS recognition algorithms. The other algorithm uses the dual representation of the set of separating lines. It consists of a smart walk in the so called Farey Fan, which can be seen as the representation of all the possible sets of separating lines for DSSs. Indeed, we take profit of the fact that the Farey Fan of order n represents in a certain way all the digital segments of length n. The computation of the characteristics of a DSL subsegment is then equivalent to the localization of a point in the Farey Fan. Using fine arithmetical properties of the fan, we design a fast algorithm of theoretical complexity O(log(n)) where n is the length of the subsegment. Experiments show that our algorithms are also efficient in practice, with a comparison to the ones previously proposed by Lachaud and Said [1]: in particular, the second one is faster in the case of "small" segments.

Introduction

Digital Straight Lines (DSL) and Digital Straight Segments (DSS) have been known for many years to be interesting tools for digital curve and shape analysis.

The applications range from simple coding to complex multiresolution analysis and geometric estimators (see for instance [START_REF] Kerautret | Meaningful scales detection along digital contours for unsupervised local noise estimation, Pattern Analysis and Machine Intelligence[END_REF] for a recent example). All these applications require to solve the so called DSS recognition problem. Many algorithms, using arithmetics [START_REF] Debled-Rennesson | A linear algorithm for segmentation of digital curves[END_REF], combinatorics [START_REF] Troesch | Interprétation géométrique de l'algorithme d'euclide et reconnaissance de segments[END_REF] or dual-space [START_REF] Dorst | Discrete straight line segments: Parameters, primitives and properties[END_REF] have been proposed to solve this problem, reaching a computational complexity of O(n) for a DSS of length n (see also [START_REF] Rosenfeld | Digital straightness: a review[END_REF] for an overview).

When no further information is known, all these algorithms are actually optimal. They at the same time decide if the set of grid points is a DSS and compute its characteristics (minimal in some sense). However, we sometimes know beforehand that the set of grid points is a DSS: the algorithm does not need to decide anymore and we can hope for a sublinear-in-time recognition algorithm. For instance, this extra information may come from the knowledge of the characteristics of a DSL containing the set of grid points. This occurs for example in [START_REF] Said | Multiscale discrete geometry[END_REF] where the multiresolution geometry of a digital object is considered. Another example concerns the digitization of a straight segment on a grid of given size: we know that the set of grid points is a DSS, but its characteristics may be much smaller than the ones of the input straight segment (and not greater than the grid size).

In [START_REF] Said | Multiscale discrete geometry[END_REF], the authors introduce the following problem: given a DSL of known characteristics and a subsegment of this DSL, compute the minimal characteristics of the DSS. The authors present two algorithms (SmartDSS and ReversedS-martDSS) to solve this problem in [START_REF] Said | Multiscale discrete geometry[END_REF][START_REF] Said | Computing the characteristics of a subsegment of a digital straight line in logarithmic time[END_REF][START_REF] Lachaud | Two efficient algorithms for computing the characteristics of a subsegment of a digital straight line[END_REF]: both use the decomposition into continuous fractions of the DSL slope and reach a logarithmic complexity.

However, a deeper search in the state-of-the-art shows that this problem is not so new. Indeed, in [START_REF] Mcilroy | A Note on Discrete Representation of Lines[END_REF], the author presents a quick sketch of a method that solves it using the Farey Fan. The announced complexity of the method is O(log 2 n) for a segment of length n. Much later, the authors of [START_REF] Charrier | Approximating a real number by a rational number with a limited denominator: A geometric approach[END_REF], try to compute the "reduction" of a straight line, which is a simplification of DSL characteristics over a bounded domain. As we will see, this reduction does not compute the minimal characteristics, but the idea is similar.

Our contribution in this paper is to demonstrate that it is possible to solve the DSL subsegment problem in logarithmic time complexity revisiting and deepening the study of the two state-of-the-art algorithms [START_REF] Mcilroy | A Note on Discrete Representation of Lines[END_REF] and [START_REF] Charrier | Approximating a real number by a rational number with a limited denominator: A geometric approach[END_REF].

The first algorithm is based on the local convex hull algorithm developed in [START_REF] Charrier | Approximating a real number by a rational number with a limited denominator: A geometric approach[END_REF] together with the framework for DSS recognition described in [START_REF] Roussillon | Algorithmes d'extraction de modèles géométriques discrets pour la représentation robuste des formes[END_REF], but we provide the theoretical results which enable to efficiently compute the minimal characteristics of a subsegment from this hull. The second algorithm, detailed in Section brings the method introduced in [START_REF] Mcilroy | A Note on Discrete Representation of Lines[END_REF] up to date: we investigate it further to provide a thoroughly defined algorithm. Moreover, we show how its complexity can be lowered to O(log(n)) with an astute use of arithmetical properties of the Farey Fan. The latter algorithm was first presented in [START_REF] Sivignon | Walking in the farey fan to compute the characteristics of a discrete straight line subsegment[END_REF] and a more detailed description is provided here. Section 2 is dedicated to the presentation of the notions used in this work. Section 3 describes the algorithm based on the local convex hull computation.

The algorithm using the dual representation, and more specifically the Farey Fan is detailed in Section 4. At the end of this section, two extensions for slightly different frameworks are presented : in particular, we show that the second algorithm can be directly and reliably applied when the input data is not a DSL but a straight line with non integer parameters. The last section is classically devoted to experimental validations.

Preliminary definitions

Digital line, segment and minimal characteristics

A Digital Straight Line (DSL for short) of integer characteristics (a, b, µ) ∈ Z 3 is the infinite set of digital points (x, y) ∈ Z 2 such that 0 ≤ axby + µ < max(|a|, |b|), with a and b relatively prime [START_REF] Debled-Rennesson | A linear algorithm for segmentation of digital curves[END_REF]. These DSL are 8-connected and often called naive. The slope of the DSL is the fraction a b and µ b is the shift at the origin. In the following, without loss of generality, we assume that 0 ≤ a ≤ b, such that, on a DSL, there is exactly one pixel for each value of x. In this context, it is easy to see that the set of pixels of a given DSL is defined by a Note that the notions of leaning points and lines are similarly defined for DSSs. DSS recognition algorithms aim at computing the minimal characteristics of a DSS, taking profit of the following fact: (a, b, µ) are the minimal characteristics of a DSS if and only if the DSS contains at least three leaning points [START_REF] Debled-Rennesson | A linear algorithm for segmentation of digital curves[END_REF]. In this case, the minimal characteristics are the characteristics of the DSS upper leaning line.

Minimal characteristics, separating lines and dual space

If we consider the digitization process related to this DSL definition, the points of the DSL L of parameters (a, b, µ) are simply the grid points (x, y) lying below or on the straight line l : axby + µ = 0 (Object Boundary Quantization), and such that the points (x, y + 1) lie above l. We say that L is the digitization of the straight line l. Note that L is also the digitization of all the straight lines of equation axby + ρ = 0 with µ ≤ ρ < µ + 1, where ρ ∈ R.

These lines separate the points X of the DSL from the points X +(0, 1), denoted by X (as in [START_REF] Roussillon | Algorithmes d'extraction de modèles géométriques discrets pour la représentation robuste des formes[END_REF]), and they are called separating lines. A similar set of lines can be defined if a DSS is considered. Let us denote by X the points of the DSS and by X the points of the DSS translated by the vector (0, 1). The separating lines are the lines which are above the upper part of the convex hull (upper convex hull for short) of the points X and strictly below the lower part of the convex hull (lower convex hull for short) of the points X (see Figure 1 for an illustration). Note that the strict constraint on the lower convex hull makes this definition slightly different from the classical definition in computational geometry. However, geometrically speaking, the set of separating lines is also bounded by the critical support lines of the two convex hulls. We actually have the property that an 8-connected set of grid points is a DSS if and only if its set of separating lines if not empty. This property is used in [START_REF] Roussillon | Algorithmes d'extraction de modèles géométriques discrets pour la représentation robuste des formes[END_REF][START_REF] Roussillon | Robust decomposition of a digital curve into convex and concave parts[END_REF] to design a fast linear in time DSS recognition algorithm. Among the separating lines, the line with integer characteristics (a, b, µ) with minimal b and minimal µ defines the minimal characteristics of the DSS. In the algorithm of [START_REF] Roussillon | Algorithmes d'extraction de modèles géométriques discrets pour la représentation robuste des formes[END_REF], the points of the DSS are added one by one and the set of separating lines is updated accordingly, but the minimal characteristics are not extracted. For a DSS of minimal characteristics (a, b, µ), the structure of the set of separating lines is perfectly known. Indeed, it is defined by the leaning points of the DSS. If U f (resp. U l ) is the leftmost (resp. rightmost) upper leaning point, and L f (resp. L l ) the leftmost (resp. rightmost) leaning point, then the set of separating lines is bounded by the lines (U f , L l +(0, 1)) and (L f +(0, 1), U l ) which are the critical support lines, and the lines (U f , U l ) and (L f + (0, 1), L l + (0, 1))

which are edges of the lower and upper convex hulls respectively. Figure 2(a)

illustrates this structure.

The set of separating lines can also be defined in a dual space, also called parameter space. In this space a straight line l : αxy + β = 0 is represented by the 2D point (α, β).

Given a DSS S, a line l : αxy + β = 0 is a separating line if and only if for all (x, y) ∈ S, 0 ≤ αxy + β < 1. This definition is strictly equivalent to the one given previously. The preimage of S is the representation of the separating lines in the dual space and is defined as

P(S) = {(α, β), 0 ≤ α ≤ 1, 0 ≤ β ≤ 1; | ∀(x, y) ∈ S, 0 ≤ αx -y + β < 1}.
As in [START_REF] Mcilroy | A Note on Discrete Representation of Lines[END_REF], let us define a ray in this space. Definition 2 (Ray). Let x and y be two integers. The ray defined by x and y is defined and denoted as follows:

L f L l U f U l l ↑ l ↓ l U l L (a) A = l ↓ * B = l U * C = l ↑ * D = l L * β α (b)
R(x, y) = {(α, β)|β = -xα + y}
x is called the slope of the ray. Note that x is not the geometrical slope of the ray but its absolute value. In the following, the order on the slopes is to be understood as the order on the absolute values of the geometrical slopes.

Given a DSS, any point (x, y) of the DSS induces two linear constraints on the preimage: αxy + β -1 < 0 and αx + βy ≥ 0. In other words, the preimage is bounded by two parallel rays: it is below the ray R(x, y + 1) and above the ray R(x, y).

This definition enables to prove that the preimage of a DSS is a convex polygon with a well-defined structure that is directly related to the leaning points and lines defined by its minimal characteristics [START_REF] Mcilroy | A Note on Discrete Representation of Lines[END_REF][START_REF] Dorst | Discrete representation of straight lines[END_REF]. Figure 2 from [START_REF] Coeurjolly | On digital plane preimage structure[END_REF] illustrates this point. In Figure 2(b), all the rays induced by the DSS pixels are depicted as dotted lines. The preimage (in gray) is the intersection of all the constraints bounded by the rays. It is for instance well known that the edges of a DSS preimage are segments of rays induced by the first and last lower and upper DSS leaning points. As a consequence, the preimage of a DSS has either three or four edges. Proposition 1 recalls in detail this specific structure.

Proposition 1 ([16]). Let P(S) be the preimage of S. Let ABCD be the polygon defined by this preimage, where A is the upper left most vertex, and the vertices are named counterclockwise. Following the notations of Figure 2 we have :

• The vertex B = l U * maps to the upper leaning line (U f , U l ) ;

• The vertex D = l L * maps to the lower leaning line (L f , L l ) translated by the vector (0, 1) in the digital space ;

• The vertex A = l ↓ * maps to the straight line (U l , L f + (0, 1)) ;

• The vertex C = l ↑ * maps to the straight line (U f , L l + (0, 1)). B is called the characteristic point of P(S).

In the rest of the paper, and especially in Section 4, the edges [AB] and

[BC] are called lower edges. This dual representation of the set of separating lines has also been used to design DSS recognition algorithms with a linear time complexity [START_REF] Dorst | Discrete straight line segments: Parameters, primitives and properties[END_REF][START_REF] O'rourke | An on-line algorithm for fitting straight lines between data ranges[END_REF].

Statement of the problem

Consider now the following problem:

Problem 1. Given a DSL L of characteristics (a L , b L , µ L ) and two points P (x P , y P ) and Q(x Q , y Q ) of this DSL (with x P < x Q ), compute the minimal characteristics (a, b, µ) of the DSS S = {(x, y)

∈ L | x P ≤ x ≤ x Q }.
Some easy cases can be withdrawn rapidly (see [START_REF] Said | Computing the characteristics of a subsegment of a digital straight line in logarithmic time[END_REF][START_REF] Lachaud | Two efficient algorithms for computing the characteristics of a subsegment of a digital straight line[END_REF]): if x Qx P ≥ 2b L , then the DSS contains at least three leaning points of the DSL and the minimal characteristics of the DSS are simply equal to (a L , b L , µ L ).

For the general case, classical DSS recognition algorithm can obviously be used. But the best complexity of such algorithms is linear in the number of pixels of the DSS. The aim here is to take profit of the extra information given by the DSL that contains the DSS to design a sublinear algorithm. In [START_REF] Said | Computing the characteristics of a subsegment of a digital straight line in logarithmic time[END_REF][START_REF] Lachaud | Two efficient algorithms for computing the characteristics of a subsegment of a digital straight line[END_REF] the authors describe two algorithms of logarithmic time complexity to solve this problem using continued fractions and a top-down or bottom-up path in the Stern-Brocot tree (see for instance [START_REF] Hayes | On the Teeth of Wheels[END_REF]). In the following sections, we present

Fast computation of local convex hulls

Rewriting of the problem

Problem 1 can be rewritten as follows:

Problem 2. Given a separating line l of characteristics (a l , b l , µ l ) and two bounding abscissas x P and x Q such that x P < x Q , find the line of minimal characteristics that is separating for the same set of points as l for the grid points between x P and x Q .

In the following, we denote by X the grid points below l on the interval [x P , x Q ] such that the points X + (0, 1), denoted by X are above l. The points X form a DSS by definition.

We present here two properties that show that the minimal characteristics of the DSS can be retrieved from some particular edges of the lower convex hull of X and the upper convex hull of X.

Property 1. Consider the remainder function r(α, β)(x, y) = αxy + β. Let l be a line of slope α and intercept β that is separating for the DSS S of minimal characteristics (a, b, µ). If the points of S are denoted by X, then the function r(α, β) restricted over the points X reaches its smallest positive value for an upper leaning point of S. Similarly, the function r(α, β) reaches its greatest value for a lower leaning point of S. Equivalently, if r(α, β) is restricted to the points X, it reaches its greatest negative value for a point L + (0, 1) where L is a lower leaning point of S, and smallest value for a point U + (0, 1) where U is an upper leaning point of S.

Proof. If α = a b , the result is straightforward. Otherwise, if α > a b , then l can be written as a linear combination of l ↑ and a line l 0 of slope α 0 = a b and intercept β 0 , lying in between l U and l L . This is easy to see in the dual space representation and illustrated in Figure 3.

If l = (1 -t)l 0 + tl ↑ where t ∈ [0, 1], the remainder function r(α, β) is equal to (1 -t)r(α 0 , β 0 ) + tr(α ↑ , β ↑ ).
The smallest positive value of r(α ↑ , β ↑ ) is equal to 0 and reached for the point U f . Similarly, the smallest positive value of r(α 0 = a b , β 0 ) is reached for all the upper leaning points of S, thus for U f . All in all, the smallest positive value of r(α, β) is reached for the point U f . In the same way, the greatest value of r(α 0 = a b , β 0 ) is reached for all the lower leaning points of S, and the greatest value of r(α ↑ , β ↑ ) is reached for the point L , which concludes the proof when α is greater then a b . The case α < a b is similar, replacing the line l ↑ by the line l ↓ .

From this property, we deduce that the grid point closest to l and below or on l is an upper leaning point for the DSS minimal characteristics. Similarly, the grid point closest to l and above l is the translation by (0, 1) of a lower leaning point for the DSS minimal characteristics. These two points are denoted by U and L respectively.

Consider now the upper convex hull of X and the lower convex hull of X.

Then from 2.2, Proposition 1 and Property 1, the DSS minimal characteristics are given by an edge passing through U or L. This leaves us with four edges to check, and the following property is used to conclude. Proof. From Proposition 1, each edge e i lies on a line that links either (i) two leaning points, both upper or both lower translated by (0, 1), of the DSS we are looking for, or (ii) one upper leaning point to the translation by (0, 1) of a lower leaning point. The minimal characteristics (a, b) are given by an edge of type (i).

L f L l U f U l l ↑ l 0 l (a) l ↓ * l U * l ↑ * l L * β α l 0 * l * (b)
Moreover, at most two out of the four edges may lie on the same line, and since there is at most two edges of type (ii), at least one edge is of type (i). If two edges are of type (i), they have the same slope. Since we have no information on the points of the edges apart from U and L, the slopes are used to discriminate between edges of type (i) and edges of type (ii). Indeed, by definition, an edge of type (ii) has a shorter period, i.e. a smallest b i than an edge of type (i), which concludes the first part of the proof. Once the characteristics (a, b) are known, we just use the fact that U is an upper leaning point of the DSS we are looking for to compute µ and thus end the proof.

Algorithm

Given a line l and two bounds x P and x Q , Properties 1 and 2 directly lead to an algorithm to solve Problem 2 in three steps:

1. compute the upper hull of the grid points X -keep the point closest to l and its two neighbour edges;

2. compute the lower hull of the grid points X -keep the point closest to l and its two neighbour edges;

3. among the four edges of slope ai bi , the one with maximal b i gives the solution.

The algorithm of [START_REF] Charrier | Approximating a real number by a rational number with a limited denominator: A geometric approach[END_REF] offers a fast solution to solve the first two steps. Indeed, the authors propose an algorithm of complexity O(log(x Q -x P )) to compute the upper and the lower convex hull of the grid points below and above a straight line and between a minimal and maximal abscissa.

In this article, the authors actually mention an algorithm to compute the "reduction" of a straight line over a domain defined by a minimal and a maximal abscissa. The aim is to reduce the coefficients of a straight line digitized over a finite domain. This algorithm is based on the computation of the critical support lines of the convex hulls computed by the algorithm. However, the "reduction" they compute is not equal to the line of minimal characteristics.

The first reason is because, as we saw, the critical support lines contain a point of X. The second reason is given hereafter and illustrated in Figure 4.

If there is no grid point lying exactly on l in the interval [x P , x Q ] this algorithm computes exactly the hulls we look for (Figure 4 (a)). However, if there exists a point R on l of abscissa x P ≤ x R ≤ x Q , then the lower convex hull computed by this algorithm is erroneous for our purpose: indeed, the computed lower convex hull of X contains the point R which is a point of X and not a point of X (Figure 4 (b)). To solve this problem, we use the fact that, for a line l : a L xb L y + µ L = 0 with integer characteristics, there is no grid point lying strictly between l and the line l : a L xb L y + µ L = 1. The trick is thus to compute the lower convex hull of the points X defined by the line l .

v 1 v 2 = v 4 v 3 v 1 v 2 = v 4 v 3 (a) R (b) v 1 v 2 = v 4 = u v 3 v 1 v 2 = v 4 v 3 (c)
Figure 4: (a) and (b) Illustration of the algorithm of [START_REF] Charrier | Approximating a real number by a rational number with a limited denominator: A geometric approach[END_REF]. In (b), the lower convex hull of X contains a point of X. (c) Modification so that the lower convex hull of X does not contain a point of X.

As in [START_REF] Charrier | Approximating a real number by a rational number with a limited denominator: A geometric approach[END_REF], the upper convex hull of the points X defined by a line l is split in two parts: the left part is composed of the edges of slope greater than the slope of l, while the right part is composed of the edges of slope lower than the slope of l. Then, the computation of the upper convex hull of the points X is completed in two steps, one for each part of the convex hull. Let localCHLeftRight(a L , b L , µ L , n) be Algorithm 1 of [START_REF] Charrier | Approximating a real number by a rational number with a limited denominator: A geometric approach[END_REF] for a line l : a L x-b L y+ µ L = 0 and abscissas in between 0 and n. This algorithm correctly returns the left part of the upper convex hull of the points X. From this convex hull, we only keep the last two points, the last one being the closest to l.

Algorithm 1: upperConvexHullOfX(a L , b L , µ L , P, Q) v 1 , v 2 = localCHLeftRight(a L , b L , r P , x Q -x P ) 1 v 1 = v 1 + P , v 2 = v 2 + P v 3 , v 4 = localCHLeftRight(a L , b L , b L -r Q , x Q -x P ) 2 v 3 = Q + (0, 1) -v 3 , v 4 = Q + (0, 1) -v 4 The closest point u is equal to v 2 or v 4 3 Return V = (v 1 , v 2 , v 3 , v 4
) removing multiple copies, and closest point u

Algorithm 2: localCHDSLSubsegment(a L , b L , µ L , P, Q) V inf , u = upperConvexHullOfX(a L , b L , µ L , P, Q) 1 V sup = upperConvexHullOfX(a L , b L , b L -r Q -1, (0, 0), Q -P ) 2 forall v in V sup do v = Q + (0, 1) -v Let E be the set of edges e, where e = (v i , v i+1 ) with v i ∈ V inf or v i ∈ V sup (b,

Fast walk in the Farey fan

This section is dedicated to the presentation of another algorithm to solve Problem 1 using the dual space representation of the set of separating lines.

Rewriting of the problem

Let us consider all the possible rays R(x, y) as defined in Section 2 with 0 < y ≤ x ≤ n. This is equivalent to considering all the linear constraints induced by all the pixels (x, y) such that 0 < y ≤ x ≤ n. Definition 3 (Farey Fan). The Farey Fan of order n, denoted by F n is defined in the (α, β) space as the arrangement of all the rays R(x, y) such that 0 ≤ y ≤ x ≤ n, and such that 0 ≤ α ≤ 1 and 0 ≤ β ≤ 1.

A facet of F n is a cell of dimension 2 of this arrangement. In the following, a point of F n stands for any point v of the (α, β) space (0 ≤ α ≤ 1 and 0 ≤ β ≤ 1) belonging to a ray, and such that the abscissa of v is a fraction of denominator smaller than or equal to n.

If P and Q are respectively the first and last point of the DSS, from the definition and the previous remarks, key Property 3 follows.

Property 3 ([9]

). For any n, there is a bijection between the facets of F n and the DSSs of length n (composed of n + 1 pixels) such that P = (0, 0) and Q = (n, y q ) with 0 < y q ≤ n. Definition 4. Let S be a DSS of length n. F acet(S) is the facet equal to P(S) in the Farey fan of order n.

Moreover, from Proposition 1 a one-to-one correspondence can be defined between a facet and the characteristic point of the facet. Definition 5. Let f be a facet of the Farey fan of order n. CP oint(f ) is the point v of f such that if v = ( p q , r q ), then (p, q, r) are the minimal characteristics

of the DSS F acet -1 (CP oint -1 (v)).
The Farey Fan of order 6 is depicted in Figure 5(a). The characteristic point of a few facets is depicted. Note that three types of facets can be identified :

• quadrilateral facets (denoted by Q, in orange in Figure 5(a)) ;

• upper triangular facets (denoted by T ↑ , in green in Figure 5(a)) ; • lower triangular facets (denoted by T ↓ , in blue in Figure 5(a)).

Let us now go back to Problem 1. After a translation of the characteristics of L such that P is set to the origin (µ L ← µ L + a L x Pb L y P ), this problem is equivalent to the following one :

Problem 3. Given a point Λ( a b , µ b ), find the point v of the Farey fan of order

n = x Q -x P such that Λ ∈ CP oint -1 (v).
In other words, the problem is to find the characteristic point of the facet of

F n containing Λ.
All in all solving Problem 3 is equivalent to performing a point location in an arrangement of lines. However, the number of facets in the Farey fan of order n (which is equal to the number of DSS of length n) is in O(n 3 ) [START_REF] Klette | Digital geometry geometric methods for digital picture analysis[END_REF][START_REF] Koplowitz | The number of digital straight lines on an n × n grid[END_REF][START_REF] Berenstein | On the number of digital straight line segments[END_REF],

and point location algorithms in such a structure are expensive in term of both time and space complexity [START_REF] Wein | CGAL -2D Arrangements[END_REF]. This brute force approach is then less efficient than classical DSS recognition algorithms [START_REF] Debled-Rennesson | A linear algorithm for segmentation of digital curves[END_REF][START_REF] Troesch | Interprétation géométrique de l'algorithme d'euclide et reconnaissance de segments[END_REF][START_REF] Kovalevsky | New definition and fast recognition of digital straight segments and arcs[END_REF][START_REF] Dorst | Discrete straight line segments: Parameters, primitives and properties[END_REF].

In the following sections, we revisit the approach proposed by [START_REF] Mcilroy | A Note on Discrete Representation of Lines[END_REF] and present an algorithm to solve Problem 3 in time complexity O(log n), without explicitly computing the Farey fan. In the next section, we recall several structural and arithmetical properties of the Farey fan, and derive some very useful corollaries.

These properties are the core of the algorithm detailed in section 4.3.

Properties of the Farey Fan

The Farey series of order n is the set of irreducible fractions in [0, 1] of denominator lower than or equal to n [START_REF] Hardy | An Introduction to the Theory of Numbers[END_REF]. The construction of the Farey series of order n from the Farey series of order n -1 is simply done as follows.

Consider two consecutive fractions p q and p q of the Farey series of order n. All the properties below are illustrated in Figure 5(b) in the Farey fan of order 6.

The first four properties are from [START_REF] Mcilroy | A Note on Discrete Representation of Lines[END_REF] and the reader is invited to consult this reference for the proofs, that are fairly simple.

Property 4 ([9]

). The abscissas of intersections of a ray R(x, y) of F n with other rays are consecutive terms of a Farey series of order max(x, nx).

In Figure 5(b), the abscissas of the intersections between the ray R(2, 1), depicted in red, and the other rays of F 6 are consecutive terms of the Farey series of order 4 = max(2, 6 -2).

Property 5 ([9]

). Let f i and f i+1 be two consecutive fractions of the Farey series of order n. In the interval f i < α < f i+1 , there is no intersection of rays.

Thus, in this interval the Farey fan is a simple ladder of rungs.

In Figure 5(b), two ladders are depicted in blue for f i = 1 3 and f i = 2 3 .

Property 6 ([9]

). Let v( p q , r q ), 0 ≤ p ≤ q ≤ n, be a point of F n . Let R(x 0 , y 0 ) be the ray of minimum slope passing through v. The other rays passing through v have a slope equal to x 0 + kq with k ∈ Z and x 0 + kq ≤ n.

In Figure 5(b), three rays go through the point ( 1 2 , 1 2 ) (in orange). The slopes of these rays are equal to x 0 = 1, 3 and 5. From this property, we can derive the following corollary.

Corollary 1. Let v( p q , r q ), 0 ≤ p ≤ q ≤ n, be a point of F n . Let R(x, y) be a ray passing through p. R is the ray of smallest slope passing through v if and only if xq < 0. It is the ray of greatest slope passing through v if and only if

x + q > n.

The following property is similar to Corollary 1 in [START_REF] Mcilroy | A Note on Discrete Representation of Lines[END_REF], but brings in more information.

Property 7. Let p q be a fraction of the Farey series of order n. The intersection between the line α = p q and F n is exactly the set of points ( p q , r q ) where r takes all the integer values between 0 and q.

Proof. We study the intersection between R(x, y) defined by the equation β = -αx + y and α = p q . We get β = -px+qy q . For 0 ≤ y ≤ x ≤ q ≤ n, the quantity -px + qy takes all the integral values in the interval [|0, q|], which ends the proof.

In Figure 5(b), the intersection between α = 4 5 (depicted in green) and F n is the set of points ( 4 5 , r 5 ) with r ∈ Z, 0 ≤ r ≤ 5. Using Properties 5 and 7, we can prove the following result to compute the ray of smallest slope in a given point.

Corollary 2. Let v( p q , r q ), 0 ≤ p ≤ q ≤ n, be a point of F n . Let p q be the fraction following p q in the Farey series of order n. The ray of smallest slope passing through v is defined by the point v and the point of coordinates v ( p q , r q )

where r is such that r q ≤ r q and r +1 q > r q .

Proof. From Property 5, F n is a ladder in the interval [ p q , p q ], which means there is no intersection of rays in this interval. From Property 7, all the rays passing through v cut the line of equation α = p q in a point v ( p q , r q ), r ∈ Z, 0 ≤ r ≤ q . Among all these rays, the ray of smallest slope is the one that passes through the point v max ( p q , rmax q ) where r max is the maximal value of r such that r q ≤ r q . It remains to prove that the two points v and v max define a ray of F n . Let x and y respectively be the slope and the intercept of the line defined by v and v max . We have to prove that: i) x and y are integers, ii) x is lower than n. For any r , and from a direct computation, we get x = rq -r q p q-pq .

Since p q and p q are consecutive fractions of a Farey series, we have p qpq = 1 and x is an integer. The same relation is used to show that y is an integer, which proves i). Let us prove ii). From the definition of r max we have r q -rmax q < 1 q , which is equivalent to rqr max q < q. Since q is lower than or equal to n, this ends the proof.

Algorithmically, two solutions are possible to compute the ray of smallest slope through a point using Corollary 2: direct computation or dichotomy. With a direct computation, we get r max = rq q : r max is the result of the integer division. A dichotomy costs O log(q ) and is not interesting since integer numbers only are involved.

Algorithm: a walk in the Farey Fan

Following Problem 3, we look for the characteristic point of the facet containing a given point Λ( a b , r b ). From Proposition 1, Section 4.1 and Property 7

we have the following characterization of the characteristic point.

Property 8. A point v( pv qv , rv qv ) is the characteristic point of a facet if and only if:

1. either v is the intersection of the two lower edges:

(a) the ray supporting the right lower edge is the one of smallest slope in v;

(b) the ray supporting the left lower edge is the one of greatest slope in v;

2. or v is on the unique lower edge and more than one ray passes through the point ( pv qv , rv+1 qv )

As in [START_REF] Mcilroy | A Note on Discrete Representation of Lines[END_REF], the algorithm consists of three steps that are detailed in the following sections :

1. Find the ladder to which Λ belongs;

2. Locate the highest ray that lies on or below Λ: this ray supports a lower edge of the facet (Section 4.3.2, Algorithm 5);

3. Walk along the ray(s) to determine the characteristic point (Section 4.3.3, Algorithm 6).

Find the ladder

Given a point Λ( a b , µ b ), finding the ladder to which Λ belongs in F n is equivalent to finding the two fractions with a denominator smaller than n closest to a b (greater and lower). We look for two fractions f = p q and g = p q such that q ≤ n, q ≤ n, f ≤ a b ≤ g, and there is no fraction of denominator smaller or equal to n neither between f and a b nor between a b and g.

The solution of this problem uses continued fractions representation of the number to be approximated (see [START_REF] Hardy | An Introduction to the Theory of Numbers[END_REF][START_REF] Davenport | The Higher Arithmetic: An Introduction to the Theory of Numbers[END_REF] for instance for an introduction on continued fractions). The solution of our problem is brought by the following Theorem (stated in [START_REF] Charrier | Approximating a real number by a rational number with a limited denominator: A geometric approach[END_REF][START_REF] Harvey | Computing two-dimensional integer hulls[END_REF], with the proof in [START_REF]Algebra: An Elementary Text-book for the Higher Classes of Secondary Schools and for Colleges[END_REF]).

Theorem 1 (as stated in [START_REF] Harvey | Computing two-dimensional integer hulls[END_REF]). Suppose we are required to find the fraction, whose denominator does not exceed n, which most closely approximates, but is no greater than, the quantity a b . If we construct a sequence of fractions containing all the odd principal convergents of a b with their corresponding intermediate convergents (if such convergents exist), then the fraction we desire is the element of this sequence with the largest denominator no greater than n.

Thus the fractions we are looking for can be found by searching in the sequence of odd convergents for the greater one, and even convergents for the lower one. A nice geometrical interpretation of the continued fraction of a number was given by Klein in 1895, and can be found in [START_REF] Davenport | The Higher Arithmetic: An Introduction to the Theory of Numbers[END_REF] or [START_REF] Har-Peled | An output sensitive algorithm for discrete convex hulls, Comput[END_REF]: if a fraction p q is represented by the grid point (q, p), then the odd (resp. even) convergents of a number α are the vertices of the lower (resp. upper) convex hull of the grid points lying above (resp. below) the line y = αx. In particular, this leads to a very simple geometrical algorithm to compute the convergents of a rational number. This algorithm called Geometric-GCD is presented in [START_REF] Har-Peled | An output sensitive algorithm for discrete convex hulls, Comput[END_REF] and has a complexity O(log(min(a, b)) for a fraction a b .

In order to compute the closest convergent with a bounded denominator, we use Geometric-GCD algorithm and simply add an upper bounding constraint of the form x ≤ n as in [START_REF] Charrier | Approximating a real number by a rational number with a limited denominator: A geometric approach[END_REF] to get the hybrid Algorithm 3. As in [START_REF] Charrier | Approximating a real number by a rational number with a limited denominator: A geometric approach[END_REF], the

Intersection(P, v,
3: BoundedGeometricGCD(l,n) l right is the vertical line x = n L = (1, 0) U = (0, 1) i = 0 while continue do if i is even then α1 = Intersection(U, L, l), α2 = Intersection(U, L, l right ) α = min(α1, α2) U = U + αL if (α = α2 or U is on l) then continue = false else α1 = Intersection(L, U, l), α2 = Intersection(L, U, l right ) α = min(α1, α2) L = L + αU if (α = α2
or L is on l) then continue = false end return L as the lower fraction, U as the greater one.

Locate a lower edge

At this point, we work in a ladder defined by two fractions f = p q and g = p q of F n and f < g. This step consists in localising Λ in the ladder by computing the highest ray under Λ in F n . In [START_REF] Mcilroy | A Note on Discrete Representation of Lines[END_REF], this step is performed as a binary search among the rays of the ladder. However, each stage of the binary search requires to solve a Diophantine equation with the extended Euclidean algorithm, reaching a total complexity of O(log 2 n).

Our algorithm, presented in Algorithm 5 and illustrated in Figure 7, also performs a dichotomy (line 3), but only on the rays of smallest slope passing through the points of abscissa p q (in red in Figure 7). The basic operation used in this part is thus the computation of the position of a point Λ with respect to a given ray: the function PositionWrtRay(point,ray) returns on, below or above.

Thanks to Property 7, the set of points of abscissa p q can be defined as on line 1, and the rays of smallest slope are computed in time O(1) in the ladder using Corollary 2 (line 2). On line 4, the ray of greatest slope is computed from the ray

Algorithm 4: FindLadder( a b ,n) l : y = a b x if b > n then 1 p q , p q = BoundedGeometricGCD(l, n) else 2 p q = a b (inf L, supL) = BoundedGeometricGCD(l, b -1) supR = (b, a) + ((b, a) -inf L)
if (supL closer to l than supR) or (the abscissa of supR is greater than n) then p q ← supL else p q ← supR return p q and p q of smallest slope thanks to Property 6: its slope is simply equal to n -(nx 0 ) (mod q) where x 0 is the slope of the ray of smallest slope. Property 6 is used in line 6. Two solutions are possible: (i) either a dichotomy is performed on the rays passing through v j+1 , once again using the function PositionWrtRay, (

or a direct computation is done. In the case of (ii), let x be the slope of the line passing through v j+1 and Λ. Let R j+1 (x j+1 , y j+1 ) be the ray of smallest slope passing through v j+1 . Let [x] be the value x j+1 + kq nearest to and lower

than x, k ∈ Z: [x] is equal to x + x j+1 -( x (mod q)) if x = x j+1
, and equal to x j+1 otherwise. If the complexity of solution (i) is straightforwardly logarithmic in the number of rays, which is smaller than q, the complexity of solution (ii) is more complicated to evaluate since it depends on the way floor and modulo functions are implemented. However, we show in Section 4.4 that dichotomy is of help when floating-point input data is considered.

In Figure 7, on the left, the point Λ is located under the ray of greatest slope passing through v j+1 (in green, line 5 in Algorithm 5), R j is returned. On the right, the point Λ is in between the rays passing through v j+1 .

Find the characteristic point

Let us denote by M and N the two points defined as the intersection between the ray R(x, y) returned by Algorithm 5 and the vertical lines defining the ladder, i.e. α = p q and α = p q as defined in Section 4.

3.1. The segment [M N ] is part of a lower edge of the facet of F n containing Λ in F n .
The first step of the algorithm detailed in Algorithm 6 is to compute the extremities of the lower edge containing [M N ]. To do so, the key point is to use Property 4 to characterize the points of intersection between a ray and other rays. Given a ray R(x, y) of the Farey Fan F n and a point v( p q , r q ) on this ray, v is the crossing point of several rays if and only if q ≤ max(x, nx). Thus, the abscissa of the left (resp. right) extremity of the lower edge is given by the term of the Farey series of order max(x, nx) preceding (resp. following p q (resp. p q ) (line 1 of Algorithm 6). This step is simply completed with a call to Proof. We assume a computing model where standard arithmetic operations are done in constant time. Finding the ladder is done using Algorithm 4 that has a complexity of O(log(n)).

Λ v j v j+1 1 4 2 7 Λ v j+1 v j 2 7 1 4
The localization of a lower edge is done with Algorithm 5: the computation of the R i (line 2) does not have to be done as a precomputation since the dichotomy (line 3) can be performed on the indices i. Thus they are computed on the fly and only when necessary, and the complexity of these two lines is in O(log(q)) with q ≤ n. The operations done in lines 4 to 5 are done in constant time. The complexity of line 6 was discussed in Section 4.3.2 and is O(log(q)) in the worst case. All in all, the complexity of Algorithm 5 is in O(log(q)).

Algorithm 6 performs the last step of the algorithm. On line 1, computing the points O and O costs O(log q) with q ≤ n (see Section 4.3.3). The computation of the mediant fraction on line 5 also has a logarithmic worst time complexity if the fraction is not irreducible. However, we can show that p + p and q + q are relatively prime: since p q and p q are succcessive terms of a Farey series, the denominator of the mediant must be strictly greater than max(q, q) ; if p+p q+q was reducible, there would exist an integer k ≥ 2 such that p+p q+q = kp kq , which contradictorily implies q ≤ max(q, q). Thus, the mediant computation is done in constant time.

All the other operations of this algorithm take O(1), which ends the proof. Algorithm 6: findCPoint( p q , p q , R) Let R(x, y) be the ray output by Algorithm 5 ( p q , p q ) = BoundedGeometricGCD(y = p q x, max(x, nx))

M N p q p q R(x, y) O O (a) M N p q p q R(x, y) O O (b) M N O p q p q R(x, y) p q O (c)
1

Let O (resp. O) be the intersection point between α = p q (resp. p q ) and R 2 if xq < 0 then Adapting the algorithm FareyFanDSLSubsegment for the computation of the minimal characteristics of 4-connected DSL subsegments is actually very easy using the following property.

From a theoretical point of view, the algorithm FareyFanDSLSubsegment proposed in Section 4 works the same. However, as often for geometrical algorithms, things get more complicated when the implementation in concerned.

Working on this issue is a research domain by itself, and a huge literature can be found about robustness in geometrical problems. See for instance [START_REF] Kettner | Classroom examples of robustness problems in geometric computations[END_REF][START_REF] Mehlhorn | Robust Geometric Computation[END_REF]. The uncertainty over the floating-point data is handled in a very classical way using an ε parameter. The way this parameter is used is illustrated in Figure 10. In (a), a centered band of height ε is defined around the line l of the Intersection(P, v, l) algorithm. If P + α v lies in the gray area, the point is said to be on the line. If it is above the gray area, the point lies above l, and below otherwise. Similarly, in (b) a vertical interval of height ε is defined around Λ in the function PositionWrtRay: if the ray R crosses the interval, the point is said to be on the ray, if R is below the interval, R is below l, and above otherwise. We could equivalently have represented the uncertainty around the ray R in (b) as we do in (a), but since the uncertainty is carried by the point Λ we find this representation more accurate. The value of the parameter ε depends on the precision of the input data. If the coordinates of Λ are known with a precision of 10 -r , then one can expect that the result of the algorithm on the floating-point point Λ(α L , β L ) is the same as the result obtained for the integer characteristics (a, b, µ) = (α L .10 r , 10 r , β L .10 r ). This suggests that a good value for ε could be 10 -r . Indeed, consider a line of equation l : ax -10 r y + µ = 0, where a and 10 r are relatively prime. Then there is no integer point in the centered band of vertical height 2 * 10 r . Thus, for instance in the case of the Intersection function, any integer point that is in the band of vertical height 10 r is closer from l than from any other line with the same slope and integer characteristics. A similar reasoning can be done for the function PositionWrtRay. An experimental validation is conducted in Section 5.2.

Last, we would like to address rapidly the consequences in the case of failures.

Concerning the function PositionWrtRay, an erroneous answer leads to a bad localization of the lower edge: the result is either the edge below or the edge above the ground truth edge. This means that the final answer will be the characteristic point of the cell just below or just above the ground truth cell.

Nevertheless, from the definition of the preimage of a DSS, it is easy to see that the difference between two DSSs the preimage of which share an edge is of exactly one pixel. The abscissa of this pixel is given by the slope of the common edge. This means that an erroneous answer of the function PositionWrtRay leads to a difference of at most one pixel between the DSS computed and the ground-truth DSS.

Concerning the function Intersection, the analysis is more complicated. Indeed, if the two fractions returned by the function are not consecutive fractions in a Farey series, the algorithm fails to find a solution. If the two fractions are not the correct ones, but are consecutive terms of a Farey series the algorithm will output a result, but the error committed is difficult to estimate. A deeper study could be done if need be.

Experimentation

Implementation and settings

The two algorithms presented in this paper are implemented and available in the generic C++ open-source library DGtal [START_REF]DGtal: Digital Geometry Tools and Algorithms Library[END_REF]. The DGtal library includes several Digital Straight Segment recognition algorithms. Moreover, the authors of [START_REF] Said | Multiscale discrete geometry[END_REF][START_REF] Said | Computing the characteristics of a subsegment of a digital straight line in logarithmic time[END_REF][START_REF] Lachaud | Two efficient algorithms for computing the characteristics of a subsegment of a digital straight line[END_REF] made there algorithms available in this library. Comparing our respective results was then an easy and robust task.

To conduct the experiments detailed below, we also reuse the protocol described in [START_REF] Lachaud | Two efficient algorithms for computing the characteristics of a subsegment of a digital straight line[END_REF] and available as a test file in DGtal. The overall protocol is governed by two parameters : N governs the value of b while n is the length of the subsegment. We recall here this protocol that includes a few minor changes.

1. Input characteristics (a, b, µ) are randomly chosen as follows:

• b is randomly chosen in the interval [N -N 2 , N + N 2 ];
• a is randomly chosen in the interval [0, b] -we also ensure that gcd(a, b) = 1;

• µ is randomly chosen in the interval [0, 2N ];

2. Points P and Q defining the subsegment are chosen as follows:

• x P is randomly chosen in the interval [0, n]; The mean computation time of these draws for each couple (N, n) is reported.

• x Q is equal to x P + n
As stated in Section 2.3, easy cases are withdrawn: when n gets bigger and close to N , the number of easy cases increases and the mean time would decrease if they were kept, bringing no useful information on the efficiency of the core of the algorithms.

Experimental correctness

The first experiment is to validate the correctness of the implementation of our two algorithms. To do so, many tests have been conducted. First, the results of the algorithms FareyFanDSLSubsegment and localCHDSLSubsegment have been directly compared to the results given by the linear-in-time algorithm

ArithmeticalDSS, as implemented in the DGtal library.

Next, the results of the algorithm FareyFanDSLSubsegment in the case of 4-connected DSLs (see Section 4.4) have been compared to the results returned by the algorithm ReversedSmartDSS.

Finally, the correctness of the floating-point implementation of the algorithm FareyFanDSLSubsegment was also evaluated. To do so, we reused the protocol defined in the previous section, but the possible values of b were powers of 10.

For each b = 10 k , random integer values of a and µ were chosen. Then the results of the algorithm on the integer data (a, b, µ) and on the floating-point (decimal) data ( a b , µ b ) were compared. In the experiments, the ε parameter was set to 10 k as explained in Section 4.4, and no errors were reported for the several millions of tests carried out.

General speed contest

The three algorithms ReversedSmartDSS, localCHDSlSubsegment and FareyFanDSLSubsegment all have a theoretical logarithmic time complexity, the logarithm being applied to different values (length of the segment, or difference of depth of the input and output continued fractions). Algorithm SmartDSS has a time complexity which depends on the sum of the quotients of the continued fraction of the output slope and on the number of pattern repetitions. In [START_REF] Lachaud | Two efficient algorithms for computing the characteristics of a subsegment of a digital straight line[END_REF], the authors showed that the algorithm SmartDSS was always slower than ReversedSmartDSS. Consequently, for the sake of clarity, we compare our algorithms with ReversedSmartDSS only and try to exhibit in which cases one algorithm may be faster than the others. First if we have a look at a particular value of N (one color), we note that FareyFanDSLSubsegment is faster than the other two for small values of n: when n becomes bigger ReversedSmartDSS gets better. For very big values of N and very small n, localCHDSlSubsegment is faster than ReversedSmartDSS but still slower than FareyFanDSLSubsegment.

Next, this graph also brings information about the behaviour of each algorithm for increasing values of N . Algorithms FareyFanDSLSubsegment and localCHDSlSubsegment seem to be insensitive to the value of N for small values of n: for a given n, the computation time is similar for all N . However, in both cases, a slight decrease of the mean computation time occurs when n gets bigger than N .

Concerning the algorithm ReversedSmartDSS, the graph reflects the fact that the complexity depends on both the value of N and the value of n: for a given n, the lower the N , the faster the computation.

To conclude on this experimental study, let us replace this work in the context of image analysis, where the DSS length is bounded by the image size.

Considering that best compact consumer cameras provide 10 megapixels images, the maximal length of DSSs in such images is bounded by a few thousands of pixels. If gigapixel images are considered, the length of the DSSs can reach values of several tens of thousands of pixels, not more for now. In both cases, our algorithm FareyFanDSLSubsegment is very competitive for any value of N .

Conclusion

We have proposed two algorithms to compute the characteristics of a DSS which is a subsegment of a DSL of known characteristics. These algorithms use two dual representations of the set of separating lines for a given set of grid points. localCHDSLSubsegment uses local computation of upper and lower convex hull to find the separating line of minimal characteristics. We provide the theoretical proof that only a few edges of the hulls are necessary to find the result. FareyFanDSLSubsegment uses the Farey fan and its numerous arithmetical properties. With this algorithm, the structure of the set of separating lines does not need to be computed since it is known through the Farey Fan, and the problem comes down to a point localization in an arrangement. We also showed that it could be extended to floating-point input data. Another perspective is to use this algorithm in fast digitization algorithms.

Suppose we want to digitize a straight segment given by its two floating-point endpoints on a grid of size n. A fast solution could be to compute the minimal characteristics of the DSS before drawing it using the arithmetical DSS definition.

  unique triplet (a, b, µ). The remainder of a DSL of characteristics (a, b, µ) for a given digital point (x, y) is the value ax-by+µ. The upper (resp. lower) leaning line of a DSL is the straight line axby + µ = 0 (resp. axby + µ = b -1). Upper (resp. lower) leaning points are the digital points of the DSL lying on the upper (resp. lower) leaning lines. A Digital Straight Segment (DSS) is a finite 8-connected part of a DSL. It can be defined by the characteristics of a DSL containing it and two endpoints P and Q. However, a DSS belongs to an infinite number of DSLs. In this context, the minimal characteristics of a DSS are the characteristics of the DSL containing it with minimal b [13, 1]. Since a DSL is defined by a unique triplet (a, b, µ), the values of a and µ are uniquely defined for a given b, and the minimal characteristics of a DSS are also uniquely defined. Definition 1 (minimal characteristics). Let L = {(a, b, µ) ∈ Z 3 , 0 ≤ a ≤ b, gcd(a, b) = 1}. For a given DSS S, we can define L S = {L ∈ L, the pixels of S all belong to the DSL of characteristics L}. Let f : L → Z, f (a, b, µ) = b. Then the minimal characteristics of S is the triplet (a S , b S , µ S ) = arg min L∈L S f (L).

Figure 1 (

 1 a) illustrates the separating lines of a DSL.

Figure 1 :

 1 Figure 1: The separating lines of a DSL(a) and of a DSS (b) are the straight lines lying in the gray area.

Figure 2 :

 2 Figure 2: (a) DSS of minimal characteristics (1, 3, 1) with its leaning points U f , U l , L f , L l . (b) Each vertex of the preimage P(S) maps to a straight line in the digital space. The vertex B( 1 3 , 1 3 ) maps to the upper leaning line, the characteristics of which are the minimal characteristics of the DSS.

  Note that point B or D may be on the line (AC). (a, b, µ) are the minimal characteristics of S if and only if B = l U * = ( a b , µ b ) (a and b relatively prime).

  two new algorithms of logarithmic time complexity to solve this problem. They use the two representations of the set of separating lines presented in Section 2.2 and illustrated in Figure 2. The algorithm of Section 3 uses properties of the upper and lower local convex hulls (Figure 2(a)) to compute the minimal characteristics, while the algorithm presented in Section 4 takes profit of a strong arithmetical structure (in the dual space, see Figure 2(b)) called Farey fan to solve the problem. The experimental section will show that our algorithms, especially the second one, have a very nice behaviour.

Figure 3 :

 3 Figure 3: A separating line l of slope greater than a b can be written as a linear combination of l ↑ and a line l 0 .

Property 2 .

 2 Consider the two edges e 0 and e 1 of the upper convex hull of X passing through U , and the two edges e 2 and e 3 of the lower convex hull of X passing through L + (0, 1). We denote by ai bi (with gcd(a i , b i ) = 1) the slopes of these edges. (a, b, µ) are the minimal characteristics of the DSS S if and only if (a, b) = (a k , b k ) where b k = max(b i ), and µ is such that the remainder of the DSS (a, b, µ) is equal to 0 on U .

  a) = (0, 0) forall e(e x , e y ) ∈ E do 3 if e x > b then (b, a) = e end µ = -au x + bu y 4 Return (a, b, µ) Algorithm 1 computes the upper convex hull of the points X lying right below a line l between two points P and Q. r P (resp. r Q ) stands for the remainder of point P (resp. Q) for the characteristics (a L , b L , µ L ). It consists of two calls to localCHLeftRight on lines 1 and 2. The call on line 2 corresponds to the computation of the hull from right to left. It returns the ordered list of vertices adjacent to the two edges defined in Property 2 and also remembers the closest point denoted by u on line 3. Algorithm 2 is the general algorithm to solve Problem 1, rewritten as Problem 2. Lines 1 and 2 compute the upper hull of the points of X and the lower hull of the points of X respectively. The result of these calls is illustrated in Figure 4(a) and (c): points v i are the one returned by the first call to upperConvexHullOfX on line 1, whereas the points v i are the one returned by the second call on line 2, using the trick presented in the previous paragraph. The loop on line 3 corresponds to the application of Property 2: the edge with maximal x-coordinate defines the sought characteristics.The point u is the point closest to l below l: from Property 1 u is an upper leaning point of the DSS we are looking for and µ is computed from this point on line 4 Concerning the time complexity, algorithm localCHLeftRight(a L , b L , µ L , n) works in O(log(n)) from [10], and is called four times. All other operations are done in constant time, so that the complexity of Algorithm localCHDSLSubsegment is O(log(n)) for a subsegment of length n.

Figure 5 :

 5 Figure 5: (a) The Farey Fan of order 6. (b) Illustration of properties 4 to 6.

  l) function computes the intersection point between the straight line defined by the point P and the vector v, and the straight line l. This point is of the form P + α v, and the function Intersection returns α . Algorithm 4 implements to computation of the ladder around the fraction a b in the Farey Fan of order n. If b is greater than n, then it consists in a simple call to the BoundedGeometricGCD algorithm. Otherwise, a direct call to the BoundedGeometricGCD algorithm would return the fraction a b for both the lower and the upper fraction, which is not the result sought. However, as shown in Figure 6, since we suppose b < 2n, a simple call to BoundedGeometricGCD also does the trick.

Algorithm

  a) + ((b, a)inf L)

Figure 6 :

 6 Figure 6: Computation of the ladder when b ≤ n: the smallest fraction greater than a b of denominator lower than or equal to n is given by one of the surrounded points.

Figure 7 :Lemma 1 .

 71 Figure 7: Illustration of Algorithm 5: the dichotomy is performed on the red rays only.

Figure 8 :

 8 Figure 8: Three cases for the lower edge [OO]: (a) all the rays passing through O (in blue) have a slope greater than x and O is the characteristic point; (b) all the rays passing through O (in blue) have a slope smaller than x and O is the characteristic point; (c) neither O nor O is the solution, and the characteristic point is found with a mediant computation.

3 return O else if x + q > n then 4 return O else Let p q = p+p q+q 5 return

 45 the intersection point between α = p q and R This algorithm solves Problem 3 in O(log(n)) where n is the order of the Farey fan. From the equivalence of Problems 1 and 3, this algorithm also solves Problem 1 in logarithmic time where n is the length of the DSS.

4. 4 .

 4 Extensions 4.4.1. 4-connected DSL subsegment In the framework presented above, the DSL and DSS considered are 8connected sets of pixels. In [8, 1] the authors consider the same problem but with 4-connected digital straight lines and segments. Their definition is similar to the 8-connected lines: a 4-connected DSL of integer characteristics (a, b, µ)is the infinite set of digital points (x, y) ∈ Z 2 such that 0 ≤ axby + µ < a + b assuming that, as before, 0 ≤ a ≤ b.

  Without going to deep in these considerations, we propose a solution to get a robust algorithm for floating-point input data. The robustness is evaluated in Section 5.2.The only functions that may cause some problems are the ones involving directly the point Λ(α L , β L ). A careful analysis of the algorithm shows that this concerns only two particular functions: the function Intersection(P, v, l) called in Algorithm 4.3.1 and the function PositionWrtRay(point,ray) called on lines 4 and 5 of Algorithm 5. In the first case, the computation of the position of an integer point with respect to a line with floating-point characteristics in involved, while in the second case, the computation of the position of a point with floating-point coordinates with respect to a line with integer characteristics is done. Let us also come back to the line 6 of Algorithm 5. We saw in Section 4.3.2 that two algorithmic choices are possible: either (i) a dichotomy using the function PositionWrtRay or (ii) a direct computation involving the computation of the slope of a line going through Λ. As we see in the few next lines, it is possible to make the function PositionWrtRay robust, while the computation of a slope involving floating-point coordinates plus its rounding seems more difficult to control, on an uncertainty point of view. Thus choice (i) seems to be the better one for floating-point input data.

Figure 10 :

 10 Figure 10: Use of the ε parameter for the Intersection (a) and the PositionWrtRay (b) functions.

  where n is the length of the subsegment; Concerning the number of draws, we randomly draw 4000 couple of values for b and a, for which 5 values of µ are chosen. For each of the 20000 triplets (a, b, µ), we draw 10 random values for x P and x Q . We repeat this process for values of N equal to 10 k with k ∈ [1, 9]. For each N , the values of n are in the interval [10, 2N ] with an increment n ← 4 3 n. Thus, for each couple of values (N, n) we proceed to 4000 * 5 * 10 = 200000 draws.

Figure 11 :

 11 Figure 11: Runtime comparison of our algorithms and one algorithm of[START_REF] Said | Multiscale discrete geometry[END_REF][START_REF] Said | Computing the characteristics of a subsegment of a digital straight line in logarithmic time[END_REF][START_REF] Lachaud | Two efficient algorithms for computing the characteristics of a subsegment of a digital straight line[END_REF]

Figure 11

 11 Figure 11 presents the results obtained for three values of N for the three algorithms, for n varying as defined in the previous section. Each value of N is represented by a color (red for 10 3 , green for 10 6 and blue for 10 9 ) and each algorithm is identified by a point type (square for FareyFanDSLSubsegment, disk for ReversedSmartDSS and triangle for localCHDSlSubsegment).

  Both algorithms have a logarithmic time complexity. Moreover, they are efficient in practice, and easy to implement. The results have been thoroughly compared to existing algorithms, both in terms of correctness of the result (to validate the implementation) and in terms of computation time.There now exists four algorithms of logarithmic time complexity to solve the DSL subsegment problem. However, no lower bound on the complexity has been proven so far. Is is possible to compute the DSL subsegment minimal characteristics in sub-logarithmic time ? Is it possible in constant time ? This is still an open question.

(Isabelle Sivignon) 1 This work has been partially supported by the French Agence Nationale de la Recherche (Grant agreement ANR-11-BS02-009)

Algorithm 5: localizeLowerEdge( p q , p q , Λ) Let v i = ( p q , i q ), i ∈ Z, 0 ≤ i ≤ q 1 Let R i (x i , y i ) be the ray of smallest slope passing through v i 2

Perform a dichotomy on the R i to compute j ∈ [|0, q -1|] such that Λ is 3 above R j and below R j+1 Let R j+1 be the ray of greatest slope passing through v j+1

Among all the rays passing through v j+1 , find the ray R which is 6 right under Λ and return R the function BoundedGeometricGCD for the line of slope p q and with an upper bounding constraint set to max(x, nx) (Algorithm 6, line 1). We get the two fractions p q and p q preceding and following p q in the Farey series of order max(x, nx). Note that since p q and p q are consecutive fractions of the Farey series of order n, the fraction p q is also greater than p q . From these two fractions we compute the two points O of R with abscissa equal to p q and O of R with abscissa equal to p q (line 2).

At this point, [OO] is a lower edge of the facet containing Λ. Then, the three cases illustrated in Figure 8 can occur: either O or O is the characteristic point (case (a) and (b)), or not (case (c)). We use Property 8 to distinguish between these cases:

the condition line 3 refers to Corollary 1;

• if R is the ray of greatest slope in O, then O is the characteristic point:

the condition line 4 refers to Corollary 1;

• otherwise, the facet is lower triangular, and the abscissa of the characteristic point is given by the mediant of the abscissas of the lower edge extremities, i.e. O and O (direct consequence of Property 4): on line 5, the mediant is computed, and the point of R with this abscissa is the characteristic point.

General algorithm and Complexity

The general algorithm FareyFanDSLSubsegment gathering all the functions presented before is summed up in Algorithm 7. It solves Problem 3, equivalent

When the DSL characteristics are floating-point numbers

Let us now consider the case when the characteristics of the DSL given as input data are not rational numbers anymore, but real numbers. We now have Λ(α L , β L ) with α L and β L in R.


Isabelle Sivignon
email: isabelle.sivignon@gipsa-lab.grenoble-inp.fr

Algorithms for Fast Digital Straight Segments Union

Keywords: Digital geometry, Union, Digital straight segment

Given two Digital Straight Segments (DSS for short) of known minimal characteristics, we investigate the union of these DSSs: is it still a DSS ? If yes, what are its minimal characteristics ? We show that the problem is actually easy and can be solved in, at worst, logarithmic time using a state-of-the-art algorithm. We moreover propose a new algorithm of logarithmic worst-case complexity based on arithmetical properties. But when the two DSSs are connected, the time complexity of this algorithm is lowered to O(1) and experiments show that it outperforms the state-of-the art one in any case.

Introduction

Digital Straight Lines (DSL) and Digital Straight Segments (DSS) have been used for many years in many pattern recognition applications involving digital curves. Whether it be for polygonal approximation or to design efficient and precise geometric estimators, a basic task is the so-called DSS recognition problem: given a set of pixels, decide whether this set is a DSS and compute its characteristics. Many linear-in-time algorithms have been proposed to solve this problem through the years. Furthermore, Constructive Solid Geometry-like operations have been considered for these objects: the intersection of two DSL has been studied in [START_REF] Said | Multiscale analysis of digital segments by intersection of 2d digital lines[END_REF][START_REF] Sivignon | Digital intersections : minimal carrier, connectiviy and periodicity properties[END_REF][START_REF] Jacob | Applications quasi-affines[END_REF][START_REF] Debled-Rennesson | A new approach to digital planes[END_REF], algorithms for the fast computation of subsegments were described in [START_REF] Lachaud | Two efficient algorithms for computing the characteristics of a subsegment of a digital straight line[END_REF][START_REF] Sivignon | Walking in the farey fan to compute the characteristics of a discrete straight line subsegment[END_REF]. Surprisingly enough, the union of DSSs has not yet been studied. The problem was raised in [START_REF] Damiand | A generic and parallel algorithm for 2d digital curve polygonal approximation[END_REF] in the context of parallel recognition of DSSs along digital contours. The recognition step was followed by a merging step where the problem of DSSs union appeared. In this work, we show how to solve this problem, both using state-of-the-art algorithm, and proposing a new and faster algorithm.

General considerations

Preliminary definitions

A Digital Straight Line (DSL for short) of integer characteristics (a, b, µ) is the infinite set of digital points (x, y) ∈ Z 2 such that 0 ≤ axby + µ < max(|a|, |b|) [START_REF] Debled-Rennesson | A linear algorithm for segmentation of digital curves[END_REF]. These DSL are 8-connected and often called naive. The fraction a b is the slope of the DSL, and µ b is the shift at the origin. In the following, without loss of generality, we assume that 0 ≤ a ≤ b. The remainder of a DSL of characteristics (a, b, µ) for a given digital point (x, y) is the value axby + µ. The upper (resp. lower) leaning line of a DSL is the straight line axby + µ = 0 (resp. axby + µ = b -1). Upper (resp. lower) leaning points are the digital points of the DSL lying on the upper (resp. lower) leaning lines. A Digital Straight Segment (DSS) is a finite 8-connected part of a DSL.

If we consider the digitisation process related to this DSL definition, the points of the DSL L of parameters (a, b, µ) are simply the grid points (x, y) lying below or on the straight line l : ax-by +µ = 0 (Object Boundary Quantization), and such that the points (x, y +1) lie above l. Otherwise said, line l separates the points X of the DSL from the points X + (0, 1) [START_REF] Roussillon | Algorithmes d'extraction de modèles géométriques discrets pour la représentation robuste des formes[END_REF], and is called separating line. More generally, for an arbitrary set of digital points X, the separating lines are the lines that separate the points X from the points X + (0, 1). In other words, the separating lines separate the upper convex hull of X from the lower convex hull of X + (0, 1). Computing the set of separating lines of two polygons is a very classical problem of computational geometry. It is well known that specific lines called critical support lines can be defined: there are the separating lines passing through a point of each polygon boundary. Critical support points are the points of the polygons belonging to critical support lines [START_REF] O'rourke | Computational Geometry in C[END_REF].

All the separating lines of a DSL have the same slope, but this is not true for arbitrary sets of digital points. The minimal characteristics of a set of digital points X are the characteristics of the separating line of minimal b and minimal µ. The set of separating lines of a DSS is well known, and the critical support points are exactly defined by the DSS leaning points: they define the minimal characteristics of the DSS.

The set of separating lines of a set of points X can also conveniently be defined in a dual space, also called parameter space. In this space a straight line l : αxy + β = 0 is represented by the 2D point (α, β). Given a set of digital points X, a line l : αxy + β = 0 is a separating line if and only if for all (x, y) ∈ X, 0 ≤ αxy + β < 1. This definition is strictly equivalent to the one given previously. The preimage of X is the representation of its separating lines in the dual space and is defined as

P(X) = {(α, β), 0 ≤ α ≤ 1, 0 ≤ β ≤ 1| ∀(x, y) ∈ X, 0 ≤ αx -y + β < 1}.
The set of separating lines of a set of pixels is an open set in the digital space, but it is a convex polygon in the dual space. In this work, this dual space will not be used explicitly in the algorithms, but we will see that this representation is convenient in some proofs. Moreover, the arrangement of all the constraints for any pixel (x, y) with y ≤ x ≤ n is called Farey Fan [START_REF] Mcilroy | A note on discrete representation of lines[END_REF] of order n: each cell of this arrangement is the preimage of a DSS of length n. Figure 1 is an illustration of the separating lines of a DSS, both in the digital space and in the dual space: they separate the points X of the DSS in black, from the points X + (0, 1) in white. Note that the edges of the preimage of a DSS are exactly supported by the dual representation of its leaning points, or equivalently its critical support points. ) maps to the upper leaning line, the characteristics of which are the minimal characteristics of the DSS.

L f L l U f U l l ↑ l ↓ l U l L (a) A = l ↓ * B = l U * C = l ↑ * D = l L * β α (b)

Setting the problem and useful properties

Consider now the following problem :

Problem 1. Given two DSSs S 1 = [P 1 Q 1 ] and S 2 = [P 2 Q 2 ]
of known minimal characteristics, decide if there exists a DSL containing both S 1 and S 2 . If yes, compute the minimal characteristics of S 1 ∪ S 2 .

If S 1 and S 2 do not belong to the same octant, then it is easy to conclude that there is no DSL containing both S 1 and S 2 . Thus, in the following S 1 and S 2 belong to the first octant, i.e. we have

S 1 (a 1 , b 1 , µ 1 ) and S 2 (a 2 , b 2 , µ 2 ) with 0 ≤ a 1 ≤ b 1 and 0 ≤ a 2 ≤ b 2 . We denote by r i (x, y) = a i x -b i y + µ i , i ∈ {1, 2} the remainder function of each DSS.
In what follows, we moreover suppose that the leaning points of S 1 and S 2 are known as input data. This is not a strong requirement since the most efficient recognition algorithms actually compute this data on the fly.

By convention, we also suppose that the abscissa of Q 2 is greater than the abscissa of Q 1 . Note that we make no assumption on the connectivity of S 1 and S 2 : the abscissa of P 2 can be lower than, equal to or greater than the abscissa of Q 1 . If the abscissa of P 2 is lower than the abscissa of P 1 , then the problem is trivial since S 1 is either a subsegment of S 2 or the union is impossible. Consequently, we also assume that the abscissa of P 2 is greater than the abscissa of P 1 .

The first part of Problem 1 consists in deciding if there exists a separating line for the set S 1 ∪ S 2 : we will say that the union is possible in this case. If so, then among all the separating lines, the final goal is to find the one with minimal characteristics.

Property 1. The preimage of S 1 ∪ S 2 is equal to the intersection between the preimages of S 1 and S 2 .

Proof. The proof is straightforward since the lines that are separating for S 1 ∪S 2 are the ones that are separating for S 1 and S 2 .

Corollary 1. The critical support points of the set of separating lines of S 1 ∪ S 2 are either upper leaning points or lower leaning points translated by (0, 1) of S 1 and S 2 . Thus, to compute the set of separating lines of S 1 ∪ S 2 , it is enough to update the set of separating lines of S 1 with the leaning points of S 2 (or conversely).

Proof. The critical support points are, in the dual space, lines supporting the edges of the preimage. From Property 1, the lines supportting the edges of P(S 1 ∪ S 2 ) are lines supporting the edges of P(S 1 ) or P(S 2 ). However, since S 1 and S 2 are DSSs, the edges of their preimages are supported by the dual representation of either upper leaning points or lower leaning points translated by (0, 1).

3 Fast Union of DSSs: an arithmetical algorithm

Fast computation of the set of separating lines

A first straightforward solution to compute the set of separating lines of S 1 ∪ S 2 is to use the state-of-the-art algorithm of O'Rourke [START_REF] O'rourke | An on-line algorithm for fitting straight lines between data ranges[END_REF], re-interpreted in the digital space by Roussillon [START_REF] Roussillon | Algorithmes d'extraction de modèles géométriques discrets pour la représentation robuste des formes[END_REF]. Whether it be in the dual space or in the digital space, these algorithms update the critical support points iteratively for each point added. Since at most four points have to be considered in our case, the algorithm is already quite efficient compared to the classical arithmetical recognition algorithm for instance. However, we propose an algorithm that is both faster and simpler to implement, in the spirit of the arithmetical recognition algorithm.

The idea is the following: if we know that the slopes of the separating lines of S 1 ∪ S 2 are greater/lower than the slopes (given by the minimal characteristics) of S 1 and S 2 respectively, then we can conclude that some leaning points of S 1 or S 2 cannot be critical support points for S 1 ∪ S 2 .

Property 2. Let S 1 be a DSS of minimal characteristics (a 1 , b 1 , µ 1 ). Let L 1f , L 1l , U 1f , U 1l be its first and last, lower and upper leaning points. If all the separating lines of S 1 ∪ S 2 have a slope greater (res. lower) than a1 b1 , then U 1l (resp. U 1f ) and L 1f (resp. L 1l ) are not critical support points for S 1 ∪ S 2 .

Proof. If all the separating lines of S 1 ∪ S 2 have a slope lower than a1 b1 , then, in the dual space and from Property 1, P(S 1 ∪ S 2 ) is a subpart of the triangle defined by the vertices ABD (see Figure 1(b)). In particular, the edges [BC] and [DC] of P(S 1 ) supported by U * 1f and L * 1l respectively cannot be edges of P(S 1 ∪ S 2 ). Therefore, the leaning points U 1f and L 1l are not critical support points for S 1 ∪ S 2 . The proof is the same if we suppose that the separating lines all have a slope greater than a1 b1 . Note that if S 1 has three leaning points only, let's say for instance only one lower leaning point L 1 , then setting L 1l and L 1f to L 1 (L 1 is "duplicated"), the property is also valid. A similar result holds when the leaning points of S 2 are considered. However, guessing the slope of the union can be tricky, and taking into account only the DSS slopes is not enough. For example, it is easy to exhibit cases where the slope of S 2 is greater than the slope of S 1 , and the slope of S 1 ∪ S 2 is nevertheless lower than both the slope of S 1 and the slope of S 2 (see Figure 2(a)). We establish hereafter some properties linking the remainder of the leaning points of S 2 and the slope of the separating lines for S 1 ∪ S 2 if they exist.

S 1 (2, 5, 0) S 2 (1, 2, -3) S(3, 8, 1) (a) α β a 1 b 1 µ 1 b 1 Preimage of S 1 Constraints related to P ∈ S 2 with r 1 (P ) ≥ b 1 A B C D (b) Fig. 2. (a)
The slope of S1 is equal to 2 5 and lower than the slope of S2, which is equal to 1 2 . However, the slope of S1 ∪ S2 is equal to 3 8 and smaller than both. (b) Illustration of Property 4.

Let's start with a very simple consideration. Property 3. If for all the leaning points P of S 2 (resp. S 1 ), we have 0 ≤ r 1 (P ) < b 1 (resp. 0 ≤ r 2 (P ) < b 2 ), then there exists a DSL containing S 1 ∪ S 2 and its minimal characteristics are the one of S 1 (resp. S 2 ).

Proof. In the dual space, the points B and D of the preimage of S 1 satisfy all the constraints related to the points P defined as above, which ends the proof.

With the following property, we investigate the other cases.

Property 4. Let P be a leaning point of S 2 .

if r 1 (P) ≥ b 1 , the slope of all the separating lines for S 1 ∪S 2 , if any, is lower than the slope of S 1 -if r 1 (P) < 0, the slope of all the separating lines for S 1 ∪S 2 , if any, is greater than the slope of S 1

Let P be a leaning point of S 1 .

if r 2 (P) ≥ b 2 , the slope of all the separating lines for S 1 ∪ S 2 , if any, is greater than the slope of S 2 -if r 2 (P) < 0, the slope of all the separating lines for S 1 ∪ S 2 , if any, is lower than the slope of S 2

Proof. We prove the first item, for a leaning point of S 2 with a remainder greater than or equal to b 1 . Proving the other cases is similar. Consider a point P ∈ S 2 such that r 1 (P ) ≥ b 1 . Let us consider the stripe defined by the constraints related to this point in the dual space. It is very simple to see that the point ( a1 b1 , µ1 b1 ) is above this stripe (see Figure 2(b)). Since P does not belong to S 1 , and with the assumptions made in Section 2.2 on the relative position of S 1 and S 2 , its abscissa is greater than the abscissas of all the leaning points of S 1 . This means that the intersection, if not empty, between the stripe and P(S 1 ) lies in the subspace α < a1 b1 .

Table 1 summarises the computation of the four possible critical support points combining Properties 2 and 4. Figure 6 in Appendix illustrates the first line of this table.

remainder value < 0 ≥ b P ∈ S2 r1(P ) (U f , L f ) = (U 1f , L 1l ) (U f , L f ) = (U 1l , L 1f ) P ∈ S1 r2(P ) (U l , L l ) = (U 2l , L 2f ) (U l , L l ) = (U 2f , L 2l ) Table 1.

Possible critical support points according to remainder values

At this point, we have identified four points denoted by U f , U l , L f and L l , that may be critical support points for S 1 ∪ S 2 . However, they may not be all critical support points. Since the preimage of S 1 ∪ S 2 is a convex polygon, it has at least three edges and thus, at least three out of the four possible points are indeed critical support points. Property 5 gives a way to decide whether the four points are critical support points or not. Property 5. U f and U l (resp. L f and L l ) are both critical support points if and only if L f and L l (resp. U f and U l ) belong to the DSL of directional vector U l -U f (resp. L l -L f ) and upper (resp. lower) leaning points U f and U l (resp. L f and L l ).

Proof. U f and U l are both critical support points is equivalent to say that the straight line (U f U l ) is separating for S 1 ∪ S 2 . This is also equivalent to the fact that L f and L l belong to the DSL as defined in the property statement.

If the four points are not all critical support points, the three critical support points are identified using Property 6. Property 6. Let U f , U l , L f and L l be the four possible critical support points for S 1 ∪ S 2 . If they are not all critical support points, then:

if U f and U l are both critical support points then:

• if the slope of (U f U l ) is lower than the slope of (L f L l ), L f is the third critical support point ; • otherwise, L l is the third critical support point ; -if L f and L l are both critical support points then:

• if the slope of (L f L l ) is greater than the slope of (U f U l ), U f is the third critical support point ; • otherwise, U l is the third critical support point ;

Proof. We write the proof for the case where U f and U l are both critical support points. The other case is similar.

Consider the dual representation of the points U f and U l , denoted by U * f and U * l . By hypothesis, these two lines support edges of P(S 1 ∪ S 2 ). Consider now the dual representation of the points L f + (0, 1) and L l + (0, 1), denoted by L * f + and L * l+ . The third edge of P(S 1 ∪ S 2 ) is a segment of either L * f + or L * l+ . We suppose now that the slope of (U f U l ) is lower than the slope of (L f L l ) and illustrate the rest of the proof with Figure 3. Then, the abscissa of point D = L * f + ∩ L * l+ is greater than the abscissa of point B = U * f ∩ U * l (it lies in the light-gray half-space on Figure 3). It is now easy to see that if D is above the line U * f , then both L * f + and L * l+ support edges of P(S 1 ∪ S 2 ), which is not possible by hypothesis. Then, D is below the line U * f , which implies that the third edge of the preimage is a segment of L * f + , and equivalently L f is the third critical support point.

Pulling out the minimal characteristics

In the previous section, we showed how to efficiently compute the three or four critical support points of S 1 ∪S 2 . These points also define the preimage of S 1 ∪S 2 . Until now, the results were valid whether the two DSSs were connected or not. In order to compute the minimal characteristics, we have to consider several cases.

Input DSSs are connected

We consider here the case where the first point of S 2 is either a point of S 1 or 8-connected to the last point of S 1 . In this case, if there exists a DSL containing S1 ∪ S 2 , then S 1 ∪ S 2 is a DSS of length n, the difference of abscissa between the first point of S 1 and the last point of S 2 . As a consequence, P(S 1 ∪ S 2 ) is a cell of the Farey Fan of order n, with very wellknown properties. In particular, the critical support points computed in Section 3.1 are exactly the leaning points of the DSS. S 1 last point and S 2 first point have the same ordinate We show that this case is actually as easy as the previous one. Property 7. If the last point of S 1 and the first point of S 2 have the same ordinate, then P(S 1 ∪ S 2 ) is a unique cell of the Farey Fan of order n, the difference of abscissa between the last point of S 2 and the first point of S1.

U * f U * l L * l+ L * f + B D (a) U * f U * l L * l+ L * f + B D (b)
Proof. In a DSS, the edges of the preimage are defined by the leaning points only. Actually, the preimage of a DSS is equal to the preimage of its leaning points, all the other points make no contribution. If the last point of S 1 and the first point of S 2 have the same ordinate, then all the missing points between these two points also have this same ordinate. Consequently, they cannot be neither lower nor upper leaning points for any DSL containing S 1 ∪ S 2 . This proves that P(S 1 ∪ S 2 ) is the same as the preimage of the set of pixels composed of S 1 , S 2 and all the missing points between the two. Then, P(S 1 ∪ S 2 ) is the preimage of a DSS of length n, the difference of abscissa between the last point of S 2 and the first point of S1, which is similar to the previous case.

Disconnected DSSs This case is trickier since P(S 1 ∪ S 2 ) may not be a unique cell but can be a union of adjacent cells of a Farey Fan (see Figure 4 for an example). The characteristics given by the critical points may not be the minimal ones. However, from the critical support points we can easily infer the range of slopes of the separating lines. If we denote s low and s up the minimum and maximum slopes of the separating lines, the slope of the line of minimal characteristics is given by the fraction of smallest denominator between s low and s up . It is finally easy to decide which one of either U f or U l is an upper leaning point of the line of minimal characteristics (see Algorithm 1 for more details).

General algorithm

All the properties presented above are put together to design the fast union algorithm described in Algorithm 1. The algorithm returns the minimal characteristics of S 1 ∪ S 2 if the union is possible. The result is given as a directional (a) vector (b, a) and an upper leaning point. The algorithm can be decomposed into three main parts. Between line 1 and 2, the four possible critical support points are computed. The function initCriticalSupportPoint is the implementation of Table 1 and is detailed in Algorithm 2 presented in Appendix. The string parameter given in input is there to discriminate between the two lines of Table 1. At the same time, easy cases where S 1 ∪ S 2 has the same minimal characteristics as S 1 or S 2 are treated: in such cases, the variable inDSL is set to true by initCriticalSupportPoint and we can conclude directly. Then, between lines 3 and 4, the exact critical support points are computed. The function isSolution? implements Property 5 and is detailed in Algorithm 3 in Appendix: here, the string parameter tells if the upper leaning points are tested w.r.t the lower leaning points or conversely. The function thirdPoint implements Property 6: if the variable solU is true, then the first item of the property is concerned, otherwise solL is true, and the second item is concerned. The last part, between lines 5 and 7 returns the minimal characteristics of S 1 ∪ S 2 . On line 6, if S 1 and S 2 are connected or if the last point of S 1 and the first point of S 2 have the same ordinate, the result is straightforward from the critical support points. Otherwise, as explained in Section 3.2 the line of minimal characteristics is computed among all the separating lines. Function minimalSlope is then called. It can be implemented in several way, using for instance the decomposition into continued fractions, or, like in [START_REF] Sivignon | Digital intersections : minimal carrier, connectiviy and periodicity properties[END_REF]-Algorithm 3 (see also [START_REF] Harel | Fast algorithms for finding nearest common ancestor[END_REF][START_REF] Graham | Concrete Mathematics[END_REF]) using the Stern-Brocot tree.

S 1 S 2 (b) S 1 S 2 (c) S 1 S 2 (d) S 1 S 2 (e)

Complexity analysis

Lemma 1. The complexity of Algorithm 1 is O(1) when S 1 and S 2 are connected or when the last point of S 1 and the first point of S 2 have the same ordinate. Its

Algorithm 1 was implemented in C++ using the open-source library DGtal [START_REF]DGtal: Digital Geometry Tools and Algorithms Library[END_REF]. We refer to it as the FastArithmetical algorithm in the following. We compare our algorithm with two other ones. The first one is the well-known arithmetical recognition algorithm [START_REF] Debled-Rennesson | A linear algorithm for segmentation of digital curves[END_REF] a DSL of characteristics (a, b, µ) is picked up at random ; -the abscissas x 1 and x 2 of the first and last points of S 1 are randomly selected ; -the abscissa x 3 of the first point of S 2 is either equal to x 2 +1 in the connected case, or randomly selected and greater than x 2 in the disconnected one ; -the abscissa of the last point is set at a fixed distance from x 3 ; Two parameters govern this setup : maxb is the maximal value of b ; valX is the length of S 2 . b is randomly picked in the interval [1, maxb], a is drawn in the interval [1, b] and such that a and b are relatively prime, and µ in the interval [0, 2maxb]. The value of x 1 is drawn in the interval [0, maxb]. The length of S 1 (i.e. x 2x 1 ) is randomly selected in the interval [valX, 2valX], so that S 1 is always longer than S 2 . In our test, maxb is set to 1000 and valX varies from 10 to 2maxb. For each value of valX, 2000 pairs of values (a, b) are drawn. For each of them, 5 different values of µ are picked up, and then 10 different values of x 1 are tested, for a total of 10 5 tests.

When the two DSSs are connected, the first test we perform consists in verifying that the three algorithms actually compute the same minimal characteristics. Then, the performances in terms of computation time are compared. Figure 4 shows the results (logarithmic scale for both axis): the x-axis represents the value of valX, the y-axis is the mean CPU computation time for a pair of DSSs, and for each algorithm.

First we can observe that the experimental behaviour of FastArithmetical algorithm confirms the constant-time complexity. Unsurprisingly, the Arithmetical algorithm has a linear-time complexity. Concerning the StabbingLine algorithm, its performances are slightly worse than the FastArithmetical algorithm, and a slight increase of the mean computation time is observed for larger DSS lengths: this is due to the fact that a post-treatment has to be done on the result returned by this algorithm in order to compute the minimal characteristics. This post-treatment involves a gcd computation, which explains the plot. However, the main information is that the FastArithmetical algorithm gets faster than the classical Arithmetical one when the length of the smallest DSS is greater than 20. In comparison, the StabbingLine algorithm becomes faster for lengths greater than 200 only. This means that what could appear as a small gain on a constant term in comparison to the StabbingLine algorithm makes the FastArithmetical relevant in practice compared to the Arithmetical algorithm. Last, the FastArithmetical remains faster than the StabbingLine algorithm even when the two DSSs are not connected. The slight decrease of the mean computation time for long DSSs is related to the fact that the longer the DSSs, the more easy cases appear.

Conclusion

In this work, we have shown that the union of two DSSs can be very efficiently computed since it is enough to "update" the minimal characteristics of the first segment with the leaning points of the second one. To do so, we have demonstrated that a state-of-the-art algorithm -the stabbing line algorithm -can be used to compute the union in logarithmic time. Moreover, we have exhibited a number of simple arithmetical properties to design an even faster algorithm. This algorithm runs in O(log(n)) worst-time complexity, and O(1) for easy cases and the experiments have shown that the implementation concretises this complexity.

Now, an interesting question remains: what if the union is not possible ? Can we measure the "distance" between the two DSSs ? A solution would be to consider "thicker" DSSs and to compute the thickness necessary for the union to be possible. This problem seems actually to be very close the the blurred DSS recognition algorithms [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF][START_REF] Roussillon | Computation of binary objects sides number using discrete geometry, application to automatic pebbles shape analysis[END_REF], and this trail seems worthy to be followed.

Fig. 1 .

 1 Fig. 1. (a) DSS of minimal characteristics (1, 3, 1) with its leaning points U f , U l , L f , L l . (b) Each vertex of the preimage maps to a straight line in the digital space. The vertex B( 1 3 , 1 3 ) maps to the upper leaning line, the characteristics of which are the minimal characteristics of the DSS.

Fig. 3 .

 3 Fig. 3. Illustration of the proof of Property 6.

Fig. 4 .

 4 Fig. 4. When S1 and S2 are not connected, like the ones depicted in (b), P(S1 ∪ S2) may be the union of several cells: in (a), P(S1) is depicted in light gray, and the two constraints related to the leaning points of S2 are depicted in red and blue. The intersection is bordered by a dotted black line: it is composed of three cells, each one being the preimage of a DSS containing S1 ∪ S2, depicted in (c), (d), and (e).

Fig. 5 .

 5 Fig. 5. Experimental results

  , implemented in DGtal (called Arithmetical algorithm in what follows). As stated at the beginning of Section 3.1, the algorithm of O'Rourke[START_REF] O'rourke | An on-line algorithm for fitting straight lines between data ranges[END_REF] can be used to compute the set of separating lines. It was implemented in DGtal by T. Roussillon as the StabbingLine algorithm. The Arithmetical algorithm works only when the two DSSs are connected and is used as follows : the minimal characteristics are initialised with the ones of S 1 , and updated as the points of S 2 are added one by one. Concerning the StabbingLine algorithm, the preimage is initialised with the one of S 1 and updated as the leaning points of S 2 are added (Corollary 1). The result is the set of critical support points of S 1 ∪ S 2 .The experimental setup is the following:

This work was partially founded by the French Agence Nationale de la Recherche (Grant agreement ANR-11-BS02-009)

Algorithm 1: FastArithmeticalDSSUnion(DSS S 1 , DSS S 2 ) U f , U l , L f , L l critical support points of S1 ∪ S2; boolean inDSL ← false; boolean solU, solL; connected ← true if S1 and S2 are connected or the last point of S1 and the first point of S2 have the same ordinate, false otherwise;

// Four possible critical support points are known 3 (solL,aL, bL) ← isSolution? (L f ,L l ,U f ,U l ,"lower") (solU,aU , bU ) ← isSolution? (U f ,U l ,L f ,L l ,"upper") if solU = false and solL = false then return DSS(0,0,Point(0,0)) // Union of S1 and S2 is not possible else if solU = false or solL = false then // Three points only are critical

The four points are critical // At this point, the exact critical support points are known 5

otherwise, where n is the difference of abscissa between the last point of S 2 and the first point of S 1 .

Proof. If we assume a computing model where the standard arithmetic operations are done in constant time, then all the operations from line 1 to line 5 are also done in constant time. Whichever the algorithm chosen, the function minimalSlope on line 6 always requires, in a more or less direct way, the computation of the continued fractions of two fractions p q with p ≤ q ≤ n, and n is the difference of abscissa between the last point of S 2 and the first point of S 1 . This is done in O(log(n)) time (see [START_REF] Graham | Concrete Mathematics[END_REF] for instance). 6. Illustration of the first line of Table 1: the leaning points of S1 marked with a cross cannot be critical support points if the point U of S2 is added.

Appendix

// compute the characteristics defined by the points P f and P l a ← P l .y -P f .y, b ← P l .x -P f .x if type = "lower" then µ ← b -1 -aP f .x + bP f .y else µ ← -aP f .x + bP f .y // check the position of Q f and Q l w.r.t these characteristics Let r(Q) = aQ.x -bQ.y + µ if 0 ≤ r(Q f ) < b and 0 ≤ r(Q l ) < b then return (true,a, b, µ) else return (false,a, b, µ)