Phuc Ngo

Yukiko Kenmochi

Nicolas Passat

Hugues Talbot

Combinatorial structure of rigid transformations in 2D digital images

Keywords: digital rigid transformations, combinatorial structure, discrete algorithm

Introduction

Rigid transformations, (i.e., transformations based on translations and rotations) are frequently involved in the design of computer vision and image processing techniques (e.g., object tracking [START_REF] Yilmaz | Object tracking: A survey[END_REF], image registration [START_REF] Zitová | Image registration methods: A survey[END_REF]), and considered in applications related to 2D or 3D images (e.g., remote sensing, medical imaging). Despite the digital nature of the processed images, such transformations are generally performed by considering the Euclidean space (R n ) associated to the Eulerian space (Z n ) of the data. Such "partially continuous" transformations then need to be interfaced with a subsequent digitization process to finally obtain a result in Z n .

The purpose of this article is to study rigid transformations of digital images as fully discrete processes. Discrete processes for the classes of rotations and affine transformations have been studied in the literature. One can cite the quasi-shear rotation [START_REF] Andres | The quasi-shear rotation[END_REF] and quasi-affine [START_REF] Coeurjolly | Quasi-affine transformation in 3-D: Theory and algorithms[END_REF] for instance. Their approach consists of decomposing transformations into three shears and then obtain the discrete transformations. From the decomposition, the transformations preserve their bijectivity. One major drawback of this approach is that the result obtained after composing the three shears is not always identical to the discretized result of the initial transform. Here, we study an approach which preserves such discretization rather than the bijectivity. Moreover, we consider several issues related to (i) the number of possible rigid transformations given a finite subspace of Z n , (ii) the way to generate all of them, and (iii) the relations between all such transformations. Some combinatorial and algorithmic answers are provided to these questions, in the case of Z 2 . Recently, C. Hundt and M. Liśkiewicz have proposed in [START_REF] Hundt | A combinatorial geometrical approach to two-dimensional robust pattern matching with scaling and rotation[END_REF][START_REF] Hundt | On the complexity of affine image matching[END_REF][START_REF] Hundt | Affine image matching is uniform TC 0 -complete[END_REF][START_REF] Hundt | Combinatorial bounds and algorithmic aspects of image matching under projective transformations[END_REF] a discretization technique under projective transformations and some of its subclasses for the problem of 2D pattern matching. Our approach is inspired by their discretization technique and adapted for the classes of rigid transformations.

The contributions of this article are twofold. The first -more theoreticalconsists of the proposal of a combinatorial structure (namely a graph) modelling the whole space of digital rigid transformations on arbitrary subset of Z 2 of size N × N , and the links between these transformations. These links correspond to the discontinuities induced by the digitization of the continuous transformations in R 2 associated to those defined in Z 2 . This combinatorial structure presents a space complexity O(N 9 ). On the one hand, this first result provides a contribution to the field of combinatorial analysis of image transformations (where previous works, summarised in Table 1, have already been proposed for rotations [START_REF] Amir | Faster two-dimensional pattern matching with rotations[END_REF][START_REF] Thibault | Rotations in 2D and 3D discrete spaces[END_REF], scalings [START_REF] Amir | Efficient pattern matching with scaling[END_REF][START_REF] Amir | Real two dimensional scaled matching[END_REF], combined scalings and rotations [START_REF] Hundt | A combinatorial geometrical approach to two-dimensional robust pattern matching with scaling and rotation[END_REF], affine transformations [START_REF] Hundt | On the complexity of affine image matching[END_REF][START_REF] Hundt | Affine image matching is uniform TC 0 -complete[END_REF], projective and linear transformations [START_REF] Hundt | Combinatorial bounds and algorithmic aspects of image matching under projective transformations[END_REF]). On the other hand, a new concept related to the topology of digital transformations is introduced in this structure, namely a neighborhood relationship between transformations induced by their relations.

The second contribution -of a more methodological nature-is an efficient algorithm enabling the construction of this combinatorial structure, with a computational cost linear with respect to its space complexity. This algorithm, which handles real (i.e., non-rational) values related to the continuous transformations associated to the discrete ones, is however defined in a fully discrete form, allowing for exact computation, and avoiding in particular any approximations related to floating point arithmetic.

Classes of transformations

Complexity

Rotations [START_REF] Amir | Faster two-dimensional pattern matching with rotations[END_REF][START_REF] Thibault | Rotations in 2D and 3D discrete spaces[END_REF] O(N 3 ) Scalings [START_REF] Amir | Efficient pattern matching with scaling[END_REF][START_REF] Amir | Real two dimensional scaled matching[END_REF] O(N 3 ) Rotations and scalings [START_REF] Hundt | A combinatorial geometrical approach to two-dimensional robust pattern matching with scaling and rotation[END_REF] O(N 6 )

Linear transformations [START_REF] Hundt | Combinatorial bounds and algorithmic aspects of image matching under projective transformations[END_REF] O(N 12 )

Affine transformations [START_REF] Hundt | On the complexity of affine image matching[END_REF][START_REF] Hundt | Affine image matching is uniform TC 0 -complete[END_REF] O(N 18 )

Projective transformations [START_REF] Hundt | Combinatorial bounds and algorithmic aspects of image matching under projective transformations[END_REF] O(N 24 )

Table 1: Space complexity of different classes of transformations on a subspace of Z 2 of size N × N .

The article is organised as follows. Sections 2 and 3 restates background notions and useful concepts. Section 4 introduces two notions: those of tipping surfaces and tipping curves, which constitute the basis of the proposed combinatorial structure. Section 5 actually defines this structure, while Sections 6 and 7 propose an algorithm for its construction. Complexity and experiments are discussed in Section 8. Section 9 presents some applications of the proposed graph and Section 10 concludes the article.

Geometric transformations and digitization

2D images and digitization

Let V be a set called value space containing at least two elements, including one, noted ⊥, corresponding to the "background". A 2D image is a function I : R 2 → V such that

I -1 (V \ {⊥}) is finite. (As examples, if
|V | = 2, we say that I is a binary image; if |V | is equipped with a total order, we say that I is a grey-level image). Without loss of generality, we can suppose that

I -1 (V \ {⊥}) ⊆ E =] -1 2 , N -1 2 [ 2 for a given N ∈ N.
The set E is called the support of I and N is called the size of I.

A 2D digital image is defined as a function I : Z 2 → V such that I -1 (V \ {⊥}) is finite. Without loss of generality, we can suppose that

I -1 (V \{⊥}) ⊆ S = [0, N -1] 2 for a given N ∈ N. The set S is called the support of I and N is the size of I.
In the sequel, by abuse of notation, we will sometimes set I : E → V and

I : S → V , instead of I : R 2 → V and I : Z 2 → V .
The 2D digital image I can be seen as a numerical representation of the 2D image I obtained by a digitization process that associates I : S → V to

I : E → V , such that S = E ∩ Z 2 .
It then maps each point of E to exactly one point of S according to a square grid structure. This digitization process is defined via the following function:

D : E -→ S x = (x, y) -→ p = (p, q) = ([x], [y]) (1)
where [ • ] is the rounding function. Note in particular that we have I = I |S , that is I(x) = I(x) for any x ∈ S. In other words, we have

D |S = Id S .
Broadly speaking, the real points x ∈ E are mapped onto the integer point p ∈ S with respect to the Voronoi diagram of S in E. This mapping is deterministic everywhere except on the frontier of the Voronoi cells, that is, for the real points x ∈ E presenting at least one semi-integer coordinate.

We call the subset H of E that consist of all such points the half-grid :

H = E ∩ R × (Z + 1 2 ) ∪ (Z + 1 2 ) × R . (2)
Broadly speaking, the digitization function D decomposes E into unit (open) squares, namely the pixels, whose centres are points of Z 2 . In other words, the half-grid H represents the union of the boundaries of all pixels.

The correct handling of this half-grid, which can be seen as the "interpixel" space in digital imaging is sometimes of great importance, for instance in (digital) topology [START_REF] Kovalevsky | Finite topology as applied to image analysis[END_REF][START_REF] Kong | Digital topology: Introduction and survey[END_REF]. However, in the digital geometry context of this article, the way the points of H are digitized is not crucial, since they represent an infinitesimal part of E. Nevertheless, H remains of first importance in the understanding of digitization, since the discontinuities induced by this process will occur on the points located in this set.

Geometric transformations for digital images

We call geometric transformation any bijective function

T : R 2 → R 2 .
Such a geometric transformation applied on an image

I : R 2 → V provides a new transformed image I • T -1 : R 2 → V .
However it is not possible to apply directly T on a digital image I in order to define I • T -1 , since there are no guarantees that T -1 (p) ∈ Z 2 for all p ∈ Z 2 . The handling of geometric transformations for digital images then requires the definition of a function noted T -1 : Z 2 → Z 2 (which is not necessarily bijective) such that we can define the digital transformed image

I •T -1 : Z 2 → V being correct with respect to T -1 .
This can be conveniently obtained by setting:

T -1 = D • T -1 (3)
as illustrated in Figure 1. Such a transformation T is called a 2D digital image transformation.

Remark 1. In the "nearest neighbour" digitization paradigm considered in

S S E E V V D Id Id D T T -1 I I • T -1 I I • T -1 T T -1

Rigid transformations for digital images

Digital rigid transformations

From this point on, we only consider 2D rigid transformations which are composed of translations and rotations. Formally, a rigid transformation

T : R 2 → R 2 transforms any point p = (p, q) ∈ R 2 into a point p ′ ∈ R 2
according to the following relation:

T (p) = p ′ =   cos θ -sin θ sin θ cos θ     p q   +   a b   (4)
where a, b, θ ∈ R.

If we consider T : Z 2 → Z 2 , from Equations (1), ( 3) and ( 4), the transformed point T (p) = p ′ = (p ′ , q ′ ) associated to p ∈ Z 2 is then defined as:

T (p) =   p ′ q ′   =   [p cos θ -q sin θ + a] [p sin θ + q cos θ + b]   . ( 5
)
This transformation is called a 2D digital rigid transformation. Note that any rigid transformation (resp. digital rigid transformation) can be represented by a triplet of parameters (a, b, θ), noted T abθ (resp. T abθ ). We note T (resp.

T d ) the set of all rigid transformations (resp. digital rigid transformations): In the following discussion, as the complexity is expressed as an order of N , the multiplicative constant factor √ 2 does not matter. Without loss of correctness, we can assume that both S and T (S) (resp. E and T (E)) are of size N such that N is "sufficiently large" with respect to S (resp. E) and the √ 2 factor.

T = {T abθ | (a, b, θ) ∈ R 3 }, (6)
T d = {T abθ | (a, b, θ) ∈ R 3 }. (7

Discontinuities of digital rigid transformations

Let us define the function T (resp. T d ), which associates to each triplet of parameters (a, b, θ) the rigid transformation (resp. digital rigid transformation) modelled by these parameters (see ( 4) and ( 5)):

T : R 3 -→ T (a, b, θ) -→ T (a, b, θ) = T abθ , (8)
T d : R 3 -→ T d (a, b, θ) -→ T d (a, b, θ) = T abθ . (9)
The function T is continuous on R 3 . In other words, for any value δ > 0, there exists ǫ > 0 such that for all (a

1 , b 1 , θ 1 ), (a 2 , b 2 , θ 2 ) ∈ R 3 (a 1 , b 1 , θ 1 ) -(a 2 , b 2 , θ 2 ) < ǫ ⇒ T(a 1 , b 1 , θ 1 ) -T(a 2 , b 2 , θ 2 ) < δ (10)
for given norms • on R 3 and (R 3 ) R 3 = {F : R 3 → R 3 } (for instance the N ∞ and pointwise N ∞ norms, respectively).

Due to the digitization involved in the definition of T in Equation ( 5), the function T d is piecewise constant on R 3 . For the same reason, T d presents discontinuities. More precisely, there exist some points c ∈ R 3 such that for any ball B(c, ǫ) of center c and radius ǫ > 0, there exist

(a 1 , b 1 , θ 1 ), (a 2 , b 2 , θ 2 ) ∈ B(c, ǫ
) such that for some values δ > 0 (for instance δ = 1), we have

(a 1 , b 1 , θ 1 ) -(a 2 , b 2 , θ 2 ) < ǫ ∧ T d (a 1 , b 1 , θ 1 ) -T d (a 2 , b 2 , θ 2 ) ≥ δ. ( 11
)
Such points c ∈ R 3 are critical points in the parameter space of the digital rigid transforms. They characterize, in particular, some transformations which map at least one point onto the half grid (see Figure 3). These considerations lead to the following definition.

Definition 1. Let (a, b, θ) ∈ R 3 , and T abθ be the associated rigid transformation. We say that T abθ is a critical transformation if

∃p ∈ Z 2 , T abθ (p) ∈ H. (12)
We denote by T c the set of all critical transformations.

The set T c of the critical transformations obviously corresponds to the part of R 3 inducing discontinuities in T d .

Discrete rigid transformations

As a result of the discontinuity of digital rigid transformations, it is possible for the digitization D of two rigid transformations of a digital image to Note that the term discrete thus refers to the non-continuous structure of transformations for digital images. In this way, the parameter space of the rigid transformations is partitioned into the disjoint sets of the equivalence classes, such that each of which associates exactly to one DRT.

(p, q) ∈ Z 2 T abθ (λ, l + 1 2 ) ∈ H (a) (p, q) ∈ Z 2 T abθ (k + 1 2 , λ) ∈ H (b)

Tipping surfaces and tipping curves

In this section, we introduce two notions, that allow us to study the decomposition of the parameter space into the DRTs. As stated in the preceding section, the digital rigid transformations present a discontinuous property, which is characterized by the critical transformations. The later transformations are actually related to the structure of DRTs in the parameter space, since they represent the boundaries of the equivalence classes. We propose the two following notions: tipping surfaces and tipping curves which are, in fact, analytical expressions of critical rigid transformations.

Tipping surfaces

According to Definition 1, a critical rigid transformation T abθ moves at least one integer point p = (p, q) ∈ Z 2 into a half-grid point which can be either a vertical half-grid point (k + 1 2 , λ) or a horizontal half-grid point (λ, l + 1

2 ) for k, l ∈ Z and λ ∈ R, or both. The set of values (a, b, θ) ∈ R 3 associated to this critical transformation forms a surface in the parameter space, called a tipping surface. Note that a tipping surface can be represented by an integer triplet associated to such a vertical or horizontal half-grid point. Definition 4. Given an integer triplet (p, q, k) (resp. (p, q, l)), we define a vertical (resp. horizontal) tipping surface as the function Φ pqk (resp. Ψ pql ) such that The sets of all tipping surfaces defined by Formulae ( 13) and ( 14

Φ pqk : R 2 -→ R (b, θ) -→ a = k + 1 2 + q sin θ -p cos θ, (13)
Ψ pql : R 2 -→ R (a, θ) -→ b = l + 1 2 -p sin θ -q cos θ. (14)

Tipping curves

We remark that the tipping surface Φ pqk in Formula ( 13) is independent of b. Consequently, the cross-sections of Φ pqk orthogonal to the b-axis orthogonal are identical for any b ∈ R, and each has the form of a trigonometric curve.

We call these tipping curves, and we denote them by φ pqk . This means that Formula (13) defines a tipping curve as well if it is observed in the plane (a, θ). Similarly, Formula ( 14) is also considered as a tipping curve on the plane (b, θ), denoted by ψ pql .

Definition 5. Given an integer triplet (p, q, k) (resp. (p, q, l)), we define a vertical (resp. horizontal) tipping curve as a function φ pqk (resp. ψ pql ) such that

φ pqk : R -→ R θ -→ a = k + 1 2 + q sin θ -p cos θ, (15)
ψ pql : R -→ R θ -→ b = l + 1 2 -p sin θ -q cos θ. (16)
The sets of all tipping curves φ pqk and ψ pql provide two families of critical rigid transformations projected on the planes (a, θ) and (b, θ) respectively. Definition 6. Two families of tipping curves are defined for θ ∈ R by

F φ (θ) = {φ pqk (θ) | p, q, k ∈ Z}, (17)
F ψ (θ) = {ψ pql (θ) | p, q, l ∈ Z}. (18)
Note that an operation applied to a family of curves means that it is applied to all tipping curves belonging to the family. For instance,

-F φ (θ) = {-φ pqk (θ) | φ pqk (θ) ∈ F φ (θ)}.

Properties of tipping curves

So far we have studied the definitions of tipping surfaces and tipping curves. We now discuss the combinatorial structure of discrete rigid transformations by using the sets of tipping surfaces and curves. Before this, we state some properties about tipping curves. From Definition 5, we remark that F φ and F ψ have similar forms. Indeed, if we translate F ψ by π 2 with respect to θ, then we obtain F φ . For simplicity, we will thus show the properties for F φ that are valid for F ψ as well.

The first property concerns the unique representation of tipping curves. This is a direct result of the fact that the three functions {1, cos, sin} involved in Formula (15) are linearly independent.

Property 1. There exists a unique integer triplet (p, q, k) ∈ Z 3 for each tipping curve.

The next property shows the relationships between two tipping curves.

This will be used in the algorithm of construction of the combinatorial structure of discrete rigid transformations. Geometrically, two tipping curves can relate to each other in four different ways (see Figure 6): they can be identical, intersecting1 , tangent2 or not intersecting. Analytical interpretation of these relationships is expressed as follows.

Property 2. let φ pqk and φ p ′ q ′ k ′ be two tipping curves.

Setting K = k -k ′ , P = p -p ′ , Q = q -q ′ , ∆ 1 = P 2 (P 2 + Q 2 -K 2 ), ∆ 2 = Q 2 (P 2 + Q 2 -K 2 )
and ∆ = ∆ 1 + ∆ 2 , we have the following relations between φ pqk and φ p ′ q ′ k ′ :

• if P = Q = 0, (i) if K = 0 then φ pqk and φ p ′ q ′ k ′ are identical, (ii) if K = 0 then they have no intersection, • otherwise, (iii) if ∆ = 0 and |KP | ≤ P 2 + Q 2 and |KQ| ≤ P 2 + Q 2 then they are tangent, (iv) if ∆ > 0 and |KP ± √ ∆ 1 | ≤ P 2 + Q 2 and |KQ ± √ ∆ 2 | ≤ P 2 + Q 2
then they intersect, (note that ± means "+ or -") (v) otherwise, they have no intersection.

Proof (i) is easily induced by Property 1. Assume that two tipping curves φ pqk and φ p ′ q ′ k ′ are different. The solution set satisfying both equations φ pqk and φ p ′ q ′ k ′ determines the nature of the relationships between two curves; it provides either the tangent point (iii), the intersection points (iv) or the empty set (ii)(v). Supposing the equation system has a solution, then we can find a value θ satisfying the following:

k + 1 2 -p cos θ + q sin θ = k ′ + 1 2 -p ′ cos θ + q ′ sin θ.
Replacing K = kk ′ , P = pp ′ and Q = qq ′ , then we have

K -P cos θ = -Q sin θ. (19)
Squaring both sides of Equation ( 19), we obtain

K 2 -2KP cos θ + P 2 cos 2 θ = Q 2 sin 2 θ.
Because cos 2 θ + sin 2 θ = 1, we have

(P 2 + Q 2 ) sin 2 θ -2KQ sin θ + K 2 -P 2 = 0, (20)
(P 2 + Q 2 ) cos 2 θ + 2KP cos θ + K 2 -Q 2 = 0, (21)
which are quadratic equations if P 2 + Q 2 = 0. As we assume that φ pqk and 19) has no solution (ii). Otherwise, P 2 + Q 2 = 0 so that the discriminants ∆ 1 and ∆ 2 of Equations ( 20) and (21) determine the number and nature of the roots. There are three cases; the curves are:

φ p ′ q ′ k ′ are different, either P , Q or K are not equal to 0. If P = Q = 0 and K = 0, then Equation (
• tangent (iii): each of Equations ( 20) and ( 21) has one solution, if ∆ = 0

(namely, ∆ 1 = ∆ 2 = 0) and |KP | ≤ P 2 + Q 2 and |KQ| ≤ P 2 + Q 2 , since -1 ≤ sin θ ≤ 1 and -1 ≤ cos θ ≤ 1,
• intersecting (iv): at least one of Equations ( 20) and ( 21) has more

than one solution, if ∆ > 0 (namely, either ∆ 1 ≥ 0 ∧ ∆ 2 > 0 or ∆ 1 > 0 ∧ ∆ 2 ≥ 0) and |KP ± √ ∆ 1 | ≤ P 2 + Q 2 and |KQ ± √ ∆ 2 | ≤ P 2 + Q 2 , since -1 ≤ sin θ ≤ 1 and -1 ≤ cos θ ≤ 1,
• otherwise, there are no intersections (v).

Given two curves φ pqk and φ p ′ q ′ k ′ , we can find out their relationship by evaluating the integer triplets (p, q, k) and (p ′ , q ′ , k ′ ) exactly, i.e., by using only integers during computation. Indeed, the values P, Q, K, ∆ 1 , ∆ 2 and ∆ are all integers, while the inequalities As the tipping curves have the form of trigonometric functions with sin θ and cos θ, they inherit the properties of the trigonometric functions as follows.

|KP ± √ ∆ 1 | ≤ P 2 + Q 2 and |KQ ± √ ∆ 2 | ≤ P 2 + Q 2 can
Property 3. The family of tipping curves F φ is symmetric:

F φ (θ) = -F φ (θ), F φ (θ) = F φ (-θ + m π 4 ), ∀m ∈ Z.
Property 4. The family of tipping curves F φ is periodic:

F φ (θ) = F φ (θ + m π 2 ), ∀m ∈ Z, F φ (θ) = F φ (θ) + m, ∀m ∈ Z.
So far, we have considered a digital image support as a set of points of Z 2 . Under this assumption, the subdivision of the parameter space into DRTs is infinite. However a digital image is finite in practice, and we can then consider S instead of Z 2 . Given a digital image I of size N × N , we define the set of tipping surfaces with respect to I as Φ pqk and Ψ pql for 0 ≤ p, q ≤ N -1, 0 ≤ k, l ≤ N and p, q, k, l ∈ Z. Thus we have 2N 2 (N + 1)

tipping surfaces partitioning the parameter space. Tipping curves φ pqk and ψ pql are correlated to the tipping surfaces, since they are the projections of Φ pqk and Ψ pql respectively on the planes (a, θ) and (b, θ) of the parameter space.

Graph representation of discrete rigid transformations

As discussed in the previous section, the parameter space of rigid transformations is divided into DRTs, whose boundaries are the tipping surfaces. This subdivision can be represented by a graph which is defined as follows.

Definition 7. Let V be a set of vertices and E be a set of labelled edges, such that

• a vertex v ∈ V corresponds to a DRT, and

• an edge e = (u, w, f ) ∈ E, where f = Φ pqk or Ψ pql , connects a pair of DRTs {u, w} ∈ V separated by the tipping surface f , which is considered as the label of the edge.

The graph G = (V, E) is called a discrete rigid transformation graph (DRT graph) (see Figure 7).

In order to distinguish between the integer triplets of vertical and horizontal tipping surfaces, we use an integer quadruple (p, q, k, i) for modelling each tipping surface, where i = 0 or 1 indicates respectively the vertical or horizontal set.

Strictly speaking, each vertex in V represents one digital transformed image generated by the corresponding DRT. Each labelled edge in E that links two vertices represents a critical transformation. This transformation links a digital image to a neighbouring one, such that they differ by only one pixel. Note that the DRT graph G = (V, E) does not contain any geometric information, such as the parameter (a, b, θ) of rigid transformations, but only the integer quadruples (p, q, k, i) representing the tipping surfaces.

The construction of a DRT graph can be performed with the help of surface arrangements [START_REF] Sharir | Recent developments in the theory of arrangements of surfaces[END_REF]. However this implies a complexity of Ω(n 4 ) for general surfaces, where n is the number of surfaces. In this article, we propose an algorithm with a better complexity of O(n 3 ) by using the properties of tipping surfaces described in Section 4. We know that while projecting two families of tipping surfaces on the planes (a, θ) and (b, θ), we obtain the corresponding families of tipping curves defined by Formulae ( 13) and [START_REF] Kong | Digital topology: Introduction and survey[END_REF].

The combinatorial structure of DRTs in a 2D parameter space is called a 2D

DRT graph. This graph is built from tipping curves on the planes (a, θ) or (b, θ). The DRT graph can be reconstructed by combining these two 2D DRT graphs. For this, we will first describe an algorithm for building a 2D DRT graph. Then we extend it to 3D to reconstruct the complete DRT graph.

Construction of 2D discrete rigid transformation graph

This section presents a method for constructing a 2D DRT graph of a set of tipping curves generated from a given digital image I of size N × N . • each vertex v ∈ V C corresponds to a face, and

• each edge e = (u, w, φ) ∈ E C corresponds to an arc that connects two faces {u, w} ∈ V C sharing a boundary tipping curve φ.

Note that the tipping curve φ is considered as a label of the edge e ∈ E C .

Figure 8 illustrates a 2D DRT graph.

The problem of the construction of G C is related to curve arrangements [START_REF] Edelsbrunner | Arrangements of curves in the plane-topology, combinatorics, and algorithms[END_REF]. Various methods have been proposed including incremental constructions [START_REF] Edelsbrunner | Arrangements of curves in the plane-topology, combinatorics, and algorithms[END_REF], sweeping lines [START_REF] Snoeyink | Sweeping arrangements of curves[END_REF] and others. A comprehensive discussion on arrangements can be found in [START_REF] Edelsbrunner | Arrangements of curves in the plane-topology, combinatorics, and algorithms[END_REF][START_REF] Halperin | Arrangements[END_REF]. Since our curves are tipping curves, many degenerate cases can present themselves, such as tangent and multiple intersections. In addition, we are only interested in the information about faces and arcs in the arrangement. Therefore, instead of using the basic algorithm of curve arrangement, we propose a variation of the sweeping line method for constructing the graph G C . The main idea of the algorithm is that a (vertical) cut is swept throughout tipping curves, and stops at some points to construct G C incrementally. The details of algorithm and its implementation are explained in the sequel.

Principles of incremental 2D DRT graph reconstruction

Let C be a set of n tipping curves and G C denotes the partition graph of C. We define a (vertical) cut, noted γ, as a monotonic line intersecting each tipping curve in C exactly once [START_REF] Snoeyink | Sweeping arrangements of curves[END_REF]. Note that γ is an unbounded simple curve and is represented by a sequence of tipping curves which γ intersects from top to bottom, as illustrated in Figure 9. We assume that γ starts at θ = 0 and ends at θ = 2π. While moving γ, its sequence of tipping curves changes, but not continuously. Indeed, γ changes only at intersection points, called event points. When γ reaches an event point, the algorithm updates γ and constructs a part of G C . This is called an elementary step of the algorithm. As a set of event points forms a series of elementary steps, we need to maintain a sorted set of event points and make the curve γ go through them in their increasing order.

In fact, the cut γ can be also represented as a directed graph such that each edge corresponds to a tipping curve in γ and each vertex corresponds to a face bounded by two consecutive tipping curves in γ.

Definition 9. Let γ = (φ 1 , φ 2 , ..., φ n ) be the cut. The partial graph δG C = (δV C , δE C ) with respect to γ is defined as a directed graph, such that

• δV C = {v 0 , v 1 , ..., v n } is the set of vertices,

• δE C = ((v 0 , v 1 , φ 1 ), (v 1 , v 2 , φ 2 ), ..., (v n-1 , v n , φ n ))
is the ordered set of edges.

In practice, elements of δE C are also ordered in the same way as γ.

At each elementary step, the partial graph δG C is updated with respect to the change of γ. If such an operation is applied, then the sweep progresses such that δG C is partially modified and integrated in the final graph G C .

Proposition 1. Let G C be the 2D DRT graph. Then we have

G C = m i=0 δG C i ( 22
)
where δG C i is the partial graph at the i-th elementary step and m is the total number of ordered event points.

Note that a partial graph is a directed graph because we need edge directions during the update. However the final graph G C is not directed, so that we do not keep directions while integrating δG C i into G C .

Initial graph construction

The initialization step provides the graph δG C 0 with respect to the initial cut γ 0 . In fact, γ 0 is a sequence of tipping curves in C with the order defined by the following relation ≺ 0 . Definition 10. For any pair of tipping curves φ pqk (θ) and

φ p ′ q ′ k ′ (θ), a rela- tion ≺ 0 is defined as φ pqk (θ) ≺ 0 φ p ′ q ′ k ′ (θ) iff • φ pqk (0) < φ p ′ q ′ k ′ (0), or • φ pqk (0) = φ p ′ q ′ k ′ (0) and the first derivatives verify φ ′ pqk (0) < φ ′ p ′ q ′ k ′ (0), or • φ pqk (0) = φ p ′ q ′ k ′ (0) and φ ′ pqk (0) = φ ′ p ′ q ′ k ′ (0) and φ pqk (π) < φ p ′ q ′ k ′ (π).
From Definition 9, we can generate the initial graph δG C 0 of γ 0 . Note that δG C 0 corresponds to the cut at θ = ǫ, where ǫ is some very small positive value. The 2D DRT graph G C is then initialized by δG C 0 .

Incremental 2D DRT graph construction for simple cases

We first present the construction algorithm in simple cases, i.e., such that any intersection consists only of two crossing tipping curves (Figure 10(a)).

We call the other cases degenerate or non-simple (see Figure 10(b)). We will discuss how to deal with such degeneracies in Section 6.5.

Detecting and ordering event points

The following question arises in the sweeping algorithm: how to detect event points, or how to discover when an elementary step is applied? Since event points are intersections of two tipping curves, the answer is linked to Property 2 (iv). More precisely, given two tipping curves φ pqk and φ p ′ q ′ k ′ modelled by the integer triplets (p, q, k) and (p ′ , q ′ , k ′ ), they intersect if the following relation is satisfied:

∆ 1 + ∆ 2 > 0 and |KP ± √ ∆ 1 | ≤ P 2 + Q 2 and |KQ ± √ ∆ 2 | ≤ P 2 + Q 2 , where P = p -p ′ , Q = q -q ′ , K = k -k ′ , ∆ 1 = P 2 (P 2 + Q 2 -K 2 ) and ∆ 2 = Q 2 (P 2 + Q 2 -K 2 ).
All event points need to be sorted and then stored in a queue, to be handled one by one. We call this queue of event points an event queue. It is defined as Q = (E, ≺ E ) where E is a set of all event points and ≺ E is a binary relation defined on E. The event points are sorted by ≺ E as follows.

Definition 11. For any pair of event points u = (u x , u y )

and v = (v x , v y ) ∈ E, a relation ≺ E is defined as u ≺ E v iff u x < v x or u x = v x and u y < v y .
Sorting event points can be performed with exact computation. We know that the coordinates (θ, a) of event points are typically irrational numbers and that we generally cannot compute an exact value for θ. Nevertheless, in order to sort event points with respect to θ ∈ [0, 2π[, we can use the values cos θ and sin θ which are calculated from Corollary 1. We can also easily obtain the value a from Formula (15) with cos θ and sin θ. Note that a, cos θ and sin θ are all quadratic irrationals3 . Therefore, sorting event points then relies on the capacity to compare quadratic irrationals. It is known in [START_REF] Rosen | Elementary Number Theory and its Applications, 3rd Edition[END_REF] that two quadratic irrationals can be compared by an exact method. In fact, a quadratic irrational can be represented exactly using a periodic continued fraction modelled by a sequence of integers, and this representation is unique.

Moreover the comparison of periodic continued fractions can be performed in constant time [START_REF] Flajolet | Continued fractions, comparison algorithms, and fine structure constants[END_REF] (see Appendix A for more details), so that sorting all Each event point is represented by the intersecting tipping curves at that point, but not its coordinates, since they are more important for constructing G C . In simple cases, any event point in Q is thus stored as two tipping curves, represented by the integer triplets.

Elementary step

An elementary step corresponds to a transposition of two curves in a cut γ around an event point, as illustrated in Figure 11. Following a change in γ, the partial graph δG C is modified.

Formally, given an event point q = {φ u , φ v }, if the cut on the left of q is denoted by γ i = (φ i 1 , . . . , φ u , φ v , . . . , φ in ), then after q we have γ i+1 = (φ i 1 , . . . , φ v , φ u , . . . , φ in ). Let δG C i and δG C i+1 denote the partial graph of γ i and γ i+1 respectively. We can generate δG C i+1 from δG C i according to the following four steps:

1. finding the current vertex w bounded by φ u and φ v ; 2. deleting two edges that are adjacent to w; Procedure 1: Elementary step for simple cases.

Input: A partial graph δG C i and an event point q = {φ u , φ v }. Output: A partial graph δG C i+1 .

1 e u ← ε(φ u ); e v ← ε(φ v );

2 {w} ← ϑ(e u ) ∩ ϑ(e v ); This procedure is called an i-th elementary step, by which the partial graph of a cut is modified. The implementation is given in Procedure 1, which requires the following two functions:

3 {w u } ← ϑ(e u ) \ {w}; {w v } ← ϑ(e v ) \ {w}; 4 ∆V C -← {w}; // w is a removed vertex 5 ∆E C -← {(w u , w, φ u ), (w, w v , φ v )}; 6 ∆V C + ← {w ′ }; // w ′ is a new vertex 7 ∆E C + ← {(w u , w ′ , φ v ), (w ′ , w v , φ u )}; 8 δG C i+1 ← δG C i \ ∆G C -∪ ∆G C + ; /* ∆G C -= (∆V C -, ∆E C -)
• ϑ(e) returns two adjacent vertices of the edge e in δV C i ;

• ε(φ) returns the edge corresponding to the tipping curve φ in δE C i .

Algorithm

We now present an algorithm for incremental 2D DRT graph construction.

The algorithm builds G C by picking event points in Q one by one. Each iteration consists of modifying the partial graph δG C according to the current cut γ (see Section 6.4.2), and then integrating δG C into G C .

Algorithm 1: Incremental construction of a 2D DRT graph in simple cases.

Input: A tipping curve set C = {φ 1 , φ 2 , ..., φ n }.

Output: A 2D DRT graph G C = (V C , E C ).
1 Initialize G C and δG C ; (see Section 6.3)

2 Generate an event queue Q; (see Section 6.4.1)

3 while Q = ∅ do 4 q ← dequeue(Q); 5 δG C ← Procedure 1 (δG C , q); 6 G C ← G C ∪ δG C ;

Incremental 2D DRT graph construction for degenerate cases

As real data with tipping curves have degeneracies, we now discuss how to deal with such cases. The algorithm for constructing a 2D DRT graph with degeneracies is similar to Algorithm 1 except for two modifications. The first of these -in step 2-consists of detecting and sorting degenerate event points to generate an event queue Q (see Section 6.5.1). The second modification -in step 5-is, despite dealing with a pair of tipping curves in simple cases, dealing with a family of tipping curves; thus the elementary step needs to be modified (see Section 6.5.2).

Detecting and sorting event points

According to their natures, the degeneracies can be classified into the following three cases, illustrated in Figure 12. More than two tipping curves can be: (i) intersecting at a single point (multiple intersection);

(ii) tangent at a single point (multiple tangent point);

(iii) tangent and/or intersecting at a single point (multiple mixed point).

In fact, those degeneracies can be detected by using Property 2 (iii)(iv)

and Corollary 1 as mentioned in Section 6.4.1. More precisely, comparing the intersections allows us to detect degenerate event points since they have the same coordinates4 at the point. Such a degenerate event point is now represented by a family of tipping curves that go through the intersection point.

Elementary step for degenerate event points

Based on the classification of degenerate event points, we can handle each case explicitly. For a multiple intersection, we swap the order of all intersecting tipping curves in the cut γ before and after this event point.

A multiple tangent is not considered as an event point because there is no change of γ around this point. The last case, a multiple mixed point, is more complicated. Observing carefully Figure 12, we remark that tipping curves are decomposed into sets sorted by tangent values; each set contains tipping curves with the same tangent. When γ passes this point, only the order of these sets of curves are reversed while the order of curves in each set is preserved. In fact, this provides the general procedure for multiple intersection cases in which each tipping curve is seen as a set of equal tangent curves.

At each event point, the elementary step consists of modifying the partial graph δG C according to the change of γ in a similar way to Procedure 1.

Note that degenerate cases, each event point contains a family of tipping curves instead of a pair as in simple cases. Let q be an event point, then q = {τ 0 , τ 1 , ..., τ m } is a family of tipping curves where each τ j = {φ j 1 , φ j 2 , ...} is a set of tipping curves with the same tangent at q. Let δG C i and δG C i+1 be respectively the partial graph with respect to γ i and γ i+1 , which go through on the left and right of the event point q. The construction of δG C i+1 from δG C i proceeds following these steps:

1. Generating two lists of tipping curves which give the order of tipping curves before and after q and storing them into two Last-In First-Out (LIFO) stacks S 1 and S 2 respectively. Details are given below.

2. Finding the initial and terminal vertices u and v for an event point q, between which δG C i changes. 3. Finding the current vertices and edges between u and v.

Replacing the current vertices and edges by the new vertices and edges between u and v.

While creating new edges in the last step, each edge is given a label of a tipping curve taken from S 2 .

We now explain how to fill two LIFO stacks S 1 and S 2 according to the order of tipping curves before and after the event point q. For this, we first need to sort the sets τ j and the tipping curves φ in each set τ j of q with the order obtained from δE C i . After sorting, we assume that q = (τ 0 , τ 1 , . . . , τ m ) and each τ j = (φ j 1 , φ j 2 , . . .) with respect to δE C i . The stacks S 1 and S 2 are generated by using the reverse sequences of τ j , which is τ j = (. . . , φ j 2 , φ j 1 ), as follows.

for j = 0 → m do for each φ ∈ τ m-j do S 1 ← push(φ); end for for each φ ∈ τ j do S 2 ← push(φ); end for end for Note that in S 1 , tipping curves have the same order as in δE C i , and in S 2 the order of sets of tipping curves are reversed while the order of curves in each set is preserved. Procedure 2 generates δG C i+1 from δG C i at the event point q.

Construction of a discrete rigid transformation graph

In this section, we present an algorithm to construct a DRT graph from a set of tipping surfaces. The basic idea is quite similar to the algorithm for constructing a 2D DRT graph, which is based on sweeping a cut and updating the DRT graph G at event points. At each elementary step, we build an additional part ∆G and then integrate it into G.

Principles of incremental DRT graph construction

As explained before, the graph of the parameter space can be constructed from its projections into the planes (a, θ) and (b, θ). We denote

G a = (V a , E a ) (resp. G b = (V b , E b )) the projection DRT graph into the planes (a, θ) (resp.
(b, θ)). From Formulae ( 13) and ( 14) of tipping curves, we derive the following proposition: Proposition 2. Let G a (resp. G b ) be a 2D DRT graph constructed from a set of tipping curves F φ (resp. F ψ ). G a and G b are isomorphic, denoted by Procedure 2: Elementary step for degenerate cases.

Input: A partial graph δG C i and an event point q. Output: A partial graph δG C i+1 .

1 Fill S 1 and S 2 as explained above.

2 φ u ← pop(S 1 ); φ v ← pop(S 1 );

3 e u ← ε(φ u ); e v ← ε(φ v ); 4 {u} ← ϑ(e u ) \ ϑ(e u ) ∩ ϑ(e v ); // u is the initial vertex 5 S 1 ← push(φ v ); S 1 ← push(φ u ); 6 φ u ← pop(S 2 ); φ v ← pop(S 2 ); 7 e u ← ε(φ u ); e v ← ε(φ v ); 8 {v} ← ϑ(e u ) \ ϑ(e v ) ∩ ϑ(e v ); // v is the terminal vertex 9 S 2 ← push(φ v ); S 2 ← push(φ u ); 10 while S 1 = ∅ do 11 φ u ← pop(S 1 ); φ v ← pop(S 1 ); 12 e u ← ε(φ u ); e v ← ε(φ v ); 13 {w} ← ϑ(e u ) ∩ ϑ(e v ); 14
{w u } ← ϑ(e u ) \ {w}; 15 ∆V C -← ∆V C -∪ {w}; // w is a removed vertex 16 ∆E C -← ∆E C -∪ {(w u , w, φ u )}; 17 if S 1 = ∅ then S 1 ← push(φ v ); 18 j ← 0; 19 φ ← pop(S 2 ); 20 ∆V C + ← {w j }; // w j is a new vertex 21 ∆E C + ← {(u, w j , φ)}; 22 while S 2 = ∅ do 23 j ← j + 1; 24 φ ← pop(S 2 ); 25 ∆V C + ← ∆V C + ∪ {w j }; // w j is a new vertex 26 if S 2 = ∅ then ∆E C + ← ∆E C + ∪ {(w j-1 , w j , φ)}; 27 else ∆E C + ← ∆E C + ∪ {(w j , v, φ)}; 28 δG C i+1 ← δG C i \ ∆G C -∪ ∆G C + ; G a ∼ G b .
Proof Translating Formula ( 14) by π 2 with respect to θ, we obtain a set of tipping curves that corresponds to the set of Formula (13), so that there exists exactly one correspondence between φ pqk and ψ pql :

b(θ + π 2 ) = l + 1 2 -p sin(θ + π 2 ) -q cos(θ + π 2 ) = l + 1 2 -p cos θ + q sin θ = a(θ).
As the sets F φ and F ψ are periodic with period π 2 (Property 4), the 2D DRT graphs in the two planes (a, θ) and (b, θ) are isomorphic.

Since G a ∼ G b , we need to construct only one graph, G a for example, and by the correspondence of tipping curves between φ pqk and ψ pql we can induce the other. The proof of Proposition 2 implies the following lemma.

Initial graph construction

The initial DRT graph G 0 = (V 0 , E 0 ) is generated from δG a 0 = (δV a 0 , δE a 0 ) and δG b 0 = (δV b 0 , δE b 0 ) as follows:

• V 0 = {(v a , v b ) | v a ∈ δV a 0 , v b ∈ δV b 0 }; • E 0 = {((u 1 , v), (u 2 , v), φ u ) | u 1 , u 2 ∈ δV a 0 , v ∈ δV b 0 , (u 1 , u 2 , φ u ) ∈ δE a 0 } ∪ {((u, v 1 ), (u, v 2 ), φ v ) | v 1 , v 2 ∈ δV b 0 , u ∈ δV a 0 , (v 1 , v 2 , φ v ) ∈ δE b 0 }.
Therefore G 0 contains (n + 1) 2 vertices and 2n(n + 1) edges, where n is the number of tipping curves of each of F φ and F ψ (see Figure 14).

Elementary step

The additional part ∆G i at each elementary step is constructed by Pro- 1 Initialize ∆G i = ∅, with ∆V i = ∅ and ∆E i = ∅;

2 foreach u ∈ ∆V a + and v ∈ δV b i do 3 ∆V i ← ∆V i ∪ {(u, v)}; 4 foreach u ∈ δV a i and v ∈ ∆V b + do 5 ∆V i ← ∆V i ∪ {(u, v)}; 6 foreach e u = (u 1 , u 2 , φ u ) ∈ δE a i do 7 foreach v ∈ ∆V b + do 8 ∆E i ← ∆E i ∪ {((u 1 , v), (u 2 , v), φ u )}; 9 foreach e v = (v 1 , v 2 , φ v ) ∈ δE b i do 10 foreach u ∈ ∆V a + do 11 ∆E i ← ∆E i ∪ {((u, v 1 ), (u, v 2 ), φ v )};
the additional parts ∆V a + and ∆V b + of δG a i and δG b i respectively. We modify Procedure 2 and obtain Procedure 2', such that it returns not only a partial graph δG C i+1 but also an adding vertex set ∆V C + , which is already calculated in Procedure 2 as an intermediate parameter.

Algorithm

The final algorithm (Algorithm 2) builds a DRT graph G by taking two event points q a ∈ Q a , q b ∈ Q b , whose coordinates are identical. At each iteration, we generate ∆G simply by using Procedure 3 and then integrates (iii) at most 2n 2n edges.

Proof (i) The number of event points is the number of intersections of two curves in C. Since two tipping curves meet at most in two points (Corollary 2), the number of event points is less than or equal to 2( n 2 ) = n(n -1). (ii) The vertices of G C correspond to the faces of the arrangement of tipping curves. If n = 1, the number of faces is 2 = 1 2 + 1, the curve forms the boundary between two faces. Let us now assume that there are (n -1) 2 + 1 faces for n -1 tipping curves. When adding the n-th curve, this curve will be divided into at most 2(n -1) + 1 arcs by the n -1 other curves, and each of these arcs will split at most one face into two. Therefore, at most 2(n -1) + 1 new faces will be created. Thus, the total number of faces (i.e., that of vertices in the 2D DRT graph) is at most (n -1)

2 + 1 + 2(n -1) + 1 = n 2 + 1.
The result follows by induction. (iii) If n = 1, there is one curve and thus we obtain 1 = 2 • 1 2 -1 edge. Let us now assume that for n -1 curves there are at most 2(n -1) 2 -(n -1) edges. When adding the n-th curve, this curve will intersect at most the n -1 previous curves. Since there are at most two intersections for each one, this creates at most 2(n -1) edges.

Moreover, the n-th curve itself will create at most 2(n -1) + 1 new edges as it has intersected at most 2(n -1) points. Thus, the total number of created edges is at most 2(n -1) 2 -(n -1) + 2(n -1) + 2(n -1) + 1 = 2n 2n. The result follows by induction.

The following property allows us to study the complexity of F φ for a digital image I with finite size. This property induces the fact that the number of possible DRTs is finite and bounded by the size of I. We say two tipping curves are vertically offset if one can be obtained by translating the other vertically. In fact, two vertically offset curves is a special case of two tipping curves with no intersection, which corresponds to (ii) in Property 2. Proof (i) The number of tipping curves φ pqk is simply the possible combinations of integer triplets (p, q, k). Since 0 ≤ p, q ≤ N -1 and 0 ≤ k ≤ N , there are N 2 (N + 1) tipping curves. (ii),(iii) From (ii) of Property 2, two trigonometric curves φ p i q i k i and φ p j q j k j are vertically offset if and only if they have p i = p j , q i = q j and k i = k j . We thus obtain N 2 sets of vertically offset tipping curves, such that each set contains N + 1 curves with different values of k. (iv) Thanks to Property 4, we only need to evaluate Formula (15) at θ = 0 and we have a(0) = k + 1 2 + q. The number of intersections is the number of different sums of q and k, which is 2N .

In practice, the various numbers of event points, vertices and edges are all lower than their upper-bound O(n 2 ), due to the degenerated cases in the arrangement of n tipping curves. From Property 5, we know that n = O(N 3 ) for an image of size N × N . Then, those complexities such as the numbers of event point, vertices and edges can be re-expressed as O(N 6 ).

As mentioned in Section 7, the construction of a DRT graph G is obtained from its projections on the planes (a, θ) and (b, θ) using two cuts. Then we notice that the initial graph has a complexity O(N 3 ) × O(N 3 ). We also know that at each elementary step, there are O(N 3 ) vertices generated. As the number of event points is O(N 6 ), in total there are O(N 6 ) × O(N 3 ) vertices added in G. This justifies the following theorem.

Experiments

We have implemented our algorithm in C++. From the experiments, the numbers of vertices (and edges) of the DRT graphs have been computed the DRT graph models the whole space of digital rigid transformations on any subset of Z 2 of size N × N . It should be mentioned that this graph does not contain any geometric parameters (a, b, θ) for the rigid transformations but only the topological information, which gives the relationship between any neighbouring transformed images. Indeed, each edge of the DRT graph is labelled by an integer quadruple (p, q, k, i) indicating that the pixel (p, q) ∈ Z 2 will cross the half-grid line, either x = k + 1 2 if i = 0 or y = k + 1 2 if i = 1. Using this information, when we move from one vertex -a transformed image-to its neighbour, only one pixel (p, q) changes. Thus, we can provide incrementally all possible images under any rigid transformations of I. Given an image I and its transformed image I ′ , we define a transition path between I and I ′ as the (possibly non unique) locus of rigid transformations transitioning I to I ′ in the parameter space (a, b, θ). By using the DRT framework, such a transition path can be represented as a connected path in the DRT graph G, called a discrete transition path. Indeed, we can find a discrete transition path in G between I and I ′ using the topological information of G providing the neighbouring relations of DRTs. Figure 17 shows a discrete transition path between two vertices of G and the sequence of transformed images corresponding to the vertices through the path. This representation of such a sequence illustrates the fact that only one pixel is changed between two incident transformed images. This allows us to define a metric on the DRT graph for measuring similarity between different transformed images, simply by counting the number of pixel changes, i.e., the length of the path.

Evaluation of discrete transition paths

Given two DRT graph vertices, it is observed in Figure 18 that their discrete transition paths are not unique. Assuming the existence of several different transition paths, we would like to find the best path which preserves some additional criteria such as topology and geometry of our object of interest in images. We know that for any two consecutive images in a discrete transition path only one pixel changes and thus it may be easy to verify our criteria incrementally along the path by checking only this pixel. By using the DRT framework, it may be possible to propose a new evaluation strategy for discrete transition paths.

Conclusion

In this article, we have introduced a combinatorial structure represented as a graph for modelling the parameter space of digital rigid transformations. This graph consists of finite sets of vertices and edges. In this graph, each vertex represents a digital transformed image, and each edge linking two vertices represents a transition changing only one pixel between two transformed images. This structure presents a space complexity of O(N 9 ), where N × N is the size of any considered subspace of Z 2 . An algorithm has also been proposed in order to define this structure in linear time with respect to this space complexity.

Experiments performed on a standard computer emphasize both the correctness of the algorithm, and the estimated time/space polynomial complexities. Due to these complexities, it remains however hardly tractable to compute the proposed combinatorial structure for large images.

However, this size limitation is not a crippling default in the case of several applications. Indeed, image processing techniques based on subimage/sample analysis can take advantage of the proposed approach, e.g., in the context of pattern matching, non-local image processing [START_REF] Buades | A review of image denoising algorithms with a new one[END_REF], or markerbased registration [START_REF] Pennec | Landmark-based registration using features identified through differential geometry[END_REF].

From a methodological point of view, further work will now involve studying ways to use the proposed combinatorial structure in multiscale strategies, most importantly in order to process large images without computing the whole data structure. Furthermore, for image registration, most existing methods [START_REF] Zitová | Image registration methods: A survey[END_REF] provide no guarantee to find a global optimal solution in general.

With our approach, we may define a new graph based metric using neighbouring relations between discrete rigid transformations, which may lead to a global optimal solution. From a theoretical point of view, extensions of the presented results to 3D digital images (following some connected works related to 3D pattern matching [23]) will also be investigated.

Figure 1 :Figure 2 :

 12 Figure 1: Geometric transformation and the associated digital image transformation (see text). First row: image support before transformation. Second row: image support after transformation.

Figure 3 :Definition 2 .Definition 3 .

 323 Figure 3: Examples of critical transformations T abθ , which map at least one integer value onto a "horizontal" (a) or "vertical" (b) half-grid point. The digital support S ⊂ Z 2 is depicted by white dots. The Voronoi cells of E associated to the points of S are depicted by squares.

Figure 4 :

 4 Figure 4: Discrete rigid transformations (DRTs) are observed as regions separated by tipping surfaces.

  ) clearly correspond to the discontinuities of digital rigid transformations. They divide the parameter space (a, b, θ) into the equivalence classes. In each of these classes, any rigid transformation leads to an identical digital transformed image. We restate that each equivalence class is represented by a DRT. As the number of critical transformations is limited by the size of a given digital image I, the number of tipping surfaces and DRTs in the parameter space are both finite. Figure4visualizes the DRTs as regions separated by tipping surfaces Φ pqk and Ψ pql in the parameter space, for p, q ∈ [0, 2] and k, l ∈ [0, 3].

Figure 5 Figure 5 :

 55 Figure 5 illustrates a part of F φ for p, q ∈ [0, 4] and k ∈ [0, 5].

Figure 6 :

 6 Figure 6: Relationships between two tipping curves, which are identical (a), intersecting (b), tangent (c) and not intersecting (d,e).

Figure 7 :

 7 Figure 7: Four tipping surfaces in the parameter space (a, b, θ) (left) and its DRT graph (right).

6. 1 .Definition 8 .

 18 Problem formalizationA finite set of n tipping curves C partitions the plane into three types of regions: a vertex is an intersection point of curves, an arc is a largest connected portion of a curve that is not intersected by any other curve, and a face is a largest connected region that is not intersected by any other curve in C. We define as follows an incident graph represented by the partition C. Given a set of tipping curves C, the 2D DRT graph G C = (V C , E C ) of C consists of a set V C of vertices and a set E C of labelled edges, such that:

Figure 8 :

 8 Figure 8: Four tipping curves on the plane (a, θ) (left) and the associated 2D DRT graph (right).

Figure 9 :

 9 Figure 9: Event points and a vertical cut γ = (φ 2 , φ 1 , φ 3 ) .

Figure 10 :

 10 Figure 10: Illustration for simple (a) and non-simple (b) cases.

  event points requires O(|Q| log |Q|) times of such a comparison. Using integer arithmetic avoids the technical problems due to the use of floating point representations. Most importantly, this allows us to prevent detecting false event points that would arise from the limited precision of floating-point representation supported in computer implementations. This would be likely to occur in degenerate cases such as multiple tangent and/or intersecting points.

1 1Figure 11 :

 111 Figure 11: Illustration of a progress of a cut by which the partial graph δG C i+1 is modified from δG C i in simple case.

/ 3 .

 3 is the subtracting part from δG C i and ∆G C + = (∆V C + , ∆E C + ) is the adding part to δG C i *replacing w by a new vertex w ′ ; 4. creating two new edges that are linked to w ′ .

Figure 12 :

 12 Figure 12: Degenerate cases: multiple intersection (a), multiple tangent point (b) and multiple mixed point (c).

1 1Figure 13 :

 113 Figure 13: Illustrating the progress of a cut by which the partial graph δG C i+1 is modified from δG C i in a degenerate case.

Lemma 7 . 1 .

 71 Let E a ⊂ R 2 (resp. E b ⊂ R 2 ) be the set of event points for the set of tipping curves F φ (resp. F ψ ). We haveE a = E b . If Q a (resp. Q b ) denotes the event queue corresponding to E a (resp. E b ), then |Q a | = |Q b |.Note that we store an event point as a list of tipping curves which generates this event point, but not as its coordinates; thus Q a = Q b even if their event points have the same coordinates.In order to construct a DRT graph G we use two vertical cuts, each of which sweeps the plane either (a, θ) or (b, θ) respectively. Thus an elementary step is performed for each pair of event points on the two planes to generate the additional part ∆G of G by combining the partial graph δG a and δG b of G a and G b respectively. The construction of δG a and δG b was described in Section 6.

Figure 14 :

 14 Figure 14: Illustration of construction of an initial DRT graph G 0 from two 2D DRT partial graphs δG a 0 and δG b 0 on the planes (a, θ) and (b, θ).

cedure 3 ,Procedure 3 :

 33 from δG a i and δG b i in a similar way to the initialization. The difference being that a new vertex is generated from a pair of a new vertex in δV a i and one of the vertices in δV b i or vice-versa. Thus a new edge connects two vertices sharing an edge in δE a i and δE b i . Procedure 3 requires Elementary step for DRT graph construction. Input: Two 2D partial graphs δG a i and δG b i and their additional vertex sets ∆V a + and ∆V b + . Output: The additional part ∆G i .

Property 5 .

 5 Given a digital image I of size N × N , the family F φ has (i) N 2 (N + 1) tipping curves, (ii) N 2 sets of vertical offset tipping curves, (iii) N + 1 vertically offset tipping curves in each set, and (iv) 2N intersections at θ = π 2 d for d ∈ Z.

Theorem 1 .

 1 The DRT graph G associated to an image of size N × N has a space complexity of O(N 9 ).

8. 2 .

 2 Run-time complexity of construction algorithmFirst, in Step 1 of Algorithm 2, the partial graphs δG a 0 and δG b 0 are initialized with a time complexity of O(N 3 log N ), when tipping curves are sorted as explained in Section 6.3. Note that the number of tipping curves is O(N 3 ). Then in Step 2, since each of δG a 0 and δG b 0 has O(N 3 ) vertices and edges, from section 7.2 we can generate G 0 in O(N 6 ). For Step 3, firstly a O(N 6 ) time complexity is needed to detect all the intersections of tipping curves, Secondly, from Proposition 3, we know that the total number of event points is O(N 6 ), thus sorting the event points takes O(N 6 log N ). As explained in Section 6.4.1, to compare any two event points we use their corresponding continued fractions. Practically, this comparison is executed in constant time[START_REF] Flajolet | Continued fractions, comparison algorithms, and fine structure constants[END_REF] so that this does not increase the time complexity of sorting event points. The sweep of a cut, in Steps 4 to 10, has O(N 6 ) iterations, for each of which Procedure 2' and Procedure 3 are executed in O(N 3 ) operations respectively. Finally, the DRT graph is constructed in O(N 9 ) time.

Figure 15 :

 15 Figure 15: The relation between image size and number of elements, (i.e., vertices or edges) in Table 2 for 2D DRT graphs (left) and DRT graphs (right).

Figure 16 shows

 16 all generated digital transformed images from a part of the DRT graph for an image of size 3 × 3 where a ∈ ] -

Figure 16 :

 16 Figure 16: A set of 231 images generated from the part of a DRT graph for a given image of size 3 × 3 where a∈ ] -1 2 , 1 2 [, b ∈ ] -1 2 , 1 2 [ and θ ∈ ]0, π4[. The first image is the identity transformation.

9. 2 .

 2 Discrete representation of the transformation of an image: discrete transition path

Figure 17 :

 17 Figure 17: A discrete transition path of the DRT graph in the parameter space (a, b, θ) (left) and a sequence of transformed images (right). Each vertex of the path represents a DRT which gives a digital transformed image. The first image corresponds to the identity transformation and the images from left to right and from top to bottom correspond to the vertices ordered in the path.

Figure 18 :

 18 Figure 18: Illustration of two different discrete transition paths with the same starting and ending vertices in a DRT graph: The two paths (left) and their sequences of images (right). The first row images correspond to the red path and the bottom ones correspond to the blue path.

  )

	Remark 2. A digital rigid transformation T (resp. a rigid transformation
	T ) maps the support S (resp. E) -of size N -of a digital image I (resp. an
	image I) onto a set T (S) (resp. T (E)) that is included in a finite (resp. bounded) part of Z 2 (resp. R 2 ) corresponding to a disk of diameter N √ 2.

  +Q 2 . Note that there are four possible combinations for sin θ and cos θ, but only two of them are valid for intersections (see Corollary 2).

	Corollary 2. Two distinct tipping curves have at most two intersections in
	[0, 2π[.
	√ P 2 Since Equation (19) has at most two solutions in [0, 2π[, we also have the ∆ 2 KP ±
	following corollary.

be verified by using squaring procedures. The proof of Property 2 implies the following corollary.

Corollary 1. Let us consider the same notations as in Property 2. Given two tipping curves φ pqk and φ p ′ q ′ k ′ , if they are • tangent, the tangent point satisfies sin θ = KQ P 2 +Q 2 and cos θ = KP P 2 +Q 2 ,

• intersecting, the intersection points satisfy sin θ = KQ± √ ∆ 1 P 2 +Q 2 and cos θ =

Two tipping curves intersect if they cross each other at a point.

Two tipping curves are tangent if both share the same tangent line at a point, i.e., the curves touch but do not cross each others.

A quadratic irrational is an irrational number that is a solution of some quadratic equations.

Note that we do not know the exact value of θ but the values sin θ and cos θ.

Algorithm 2: Construction of DRT graph.

Input: A set of tipping surfaces, i.e., two sets of tipping curves.

Output: A DRT graph G = (V, E).

1 Initialize δG a and δG b ; (see Section 6.3)

2 Initialize G = G 0 from δG a 0 and δG b 0 ; (see Section 7.2)

3 Generate Q a and Q b ; (see Section 6.5.1)

Complexity analysis and experiments

Space complexity of DRT graph

The complexity of a 2D DRT graph, i.e., the numbers of its vertices and edges, is obtained by counting the number of event points. for images of sizes varying from 1 × 1 to 9 × 9. The experiments were carried out on a personal computer equipped with a processor 3.0GHz Intel R

Core TM 2 Duo and 4GB of memory. The results, shown in Table 2 and Figure 15, validate the theoretical O(N 9 ) space complexity stated in the previous theorem.

Applications

This section illustrates the practical applicability of our proposed combinatorial structure for DRTs.

Generation of all digital transformed images

Given a digital image I of size N × N , we firstly show how to generate all possible transformed images of I under rigid transformations. We know that

Appendix A. Exact comparison of quadratic irrationals

This appendix describes an algorithmic process enabling to compare two quadratic irrationals without numerical approximation. In [START_REF] Rosen | Elementary Number Theory and its Applications, 3rd Edition[END_REF][START_REF] Khinchin | Continued Fractions[END_REF], it was proved that a quadratic irrational can be rewritten as a periodic continued fraction. More formally, for any quadratic irrational

where p, q, r ∈ Z, q > 0, r = 0 we have a periodic continued fraction:

where a 0 ∈ Z and a i ∈ N for i = 1 . . . n. Such a periodic continued fraction is unambiguously modelled by a finite sequence of integers (a 0 , a 1 , . . . , a n ), denoted by [a 0 ; a 1 , a 2 , . . . a n ]. Based on this formulation, the comparison between two quadratic irrationals can be performed as follows.

Given

let k ∈ N be the smallest index for which

can be easily checked by comparing the values p i , q i and r i ), the order between

Algorithm 3: Comparison of two quadratic irrationals Input: (p 1 , q 1 , r 1 ), (p 2 , q 2 , r 2 ) representing two quadratic irrationals Q 1 and Q 2 .

Output: Value in {<, =, >} denoting the relation between Q 1 and Q 2 .

if (p 1 , q 1 , r 1 ) = (p 2 , q 2 , r 2 ) then return =;

E ← (-1) i-1 (a i 1a i 2 ); i ← i + 1; if E > 0 then return >; else return <; in [START_REF] Podsypanin | Length of the period of a quadratic irrational[END_REF] that the length of repeating block for Q is O( √ q ln q). In the worst case, all terms of the repeating block are compared, i.e., we have O( √ q ln q)

comparisons. However, it is proved in [START_REF] Flajolet | Continued fractions, comparison algorithms, and fine structure constants[END_REF] that the comparison of continued fractions has an average-case complexity of O(1).