Thibault Yohan

Keywords:

Petite pensé. En anglais, il existe le verbe "to rotate", en français son équivalent serait le verbe rotationner, sa dénition serait eectuer ou faire eectuer une rotation. Cependant le verbe "rotationner" n'existe pas dans la langue française. Et pourtant, il tourne.

The discrete geometry is to classical geometry what the language is to thought, i.e. an imperfect means to represent the reality. It took centuries for the language to evolve in a way almost capable to faithfully describe our though.

Today the discrete geometry tries to do the same thing with the continuous geometry.

The continuous geometry is a mathematical model that cannot be correctly or exactly reproduced in the real world and in computer science. A simple example is the famous number π. The theoretical mathematic model supposes an exact value of this number, however, the representation of a circle on the ground with a rope or on a sheet of paper by a compass can only give an approximated value of π, whatever the diameter of the circle, the size of the rope or the precision of the compass. In computer science, for any approximation of π used during computations, results will always be an approximation.

Today, one of the biggest challenge in computer science is to nd new methods so that computers can represent reality as faithfully as possible. Regarding geometry, we strongly believe that these methods belong to the discrete geometry.

List of Figures

. 2.8

The quasi-shear rotation and its two intermediate results. . . . . . . . . . 2.9 A pixel and (a) its four 4-connected pixels, (b) its eight 8-connected pixels in a square grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 A pixel and its six 6-connected pixels in a hexagonal grid . . . . . . . . . 2.11 A voxel and (a) its six 6-connected voxels, (b) its eighteen 18-connected voxels and (c) its twenty-six 26-connected voxels in a cubic grid. . . . . . . 2.12 Conversion from a square grid into a triangle grid . . . . . . . . . . . . . . 2.13 Conversion from a square grid into a hexagonal grid . . . . . . . . . . . . 2. [START_REF] Louis | Über die Zahl π[END_REF] The percentage of points remaining after a rotation in a triangular grid. . 2.15 Comparison of the average conservation of neighborhood during a rotation angle between the square grid and (a) 4-connectivity on a triangle grid, (b) 8-connectivity on a triangle grid,(c) 4-connectivity on a hexagonal grid and (d) 8-connectivity on a hexagonal grid. . . . . . . . . . . . . . . . . . 2. [START_REF] Hansen | Real-time scene stabilization and mosaic construction[END_REF] The percentage of points remaining after rotation in a hexagonal grid. . . 2.17 Result of four rotations of the image (a) by an angle of 15 • (b) using a square grid, (c) using a hexagonal grid and (d) using a triangle grid. . . . 2.18 Result of four rotations of the image (a) by an angle of 30 • (b) using a square grid, (c) using a hexagonal grid and (d) using a triangle grid. . . . 2. [START_REF] Paeth | A fast algorithm for general raster rotation[END_REF] Result of four rotations of the image (a) by an angle of 45 • (b) using a square grid, (c) using a hexagonal grid and (d) using a triangle grid. . . . xiv xv 3.1 A hinge angle α(P, Q, K) (left) and four symmetrical hinge angles (right).

3.2

The corners of H (⃗ q), namely, four corners of a pixel around ⃗ q. . . . . . . [START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF] The ve dierent shapes of convexels, which are constructed as the intersections between a rotation plane and voxels [START_REF] Andres | New methods in oblique slice generation[END_REF]. . . . . . . . . . . . . . . . 4.5 Two multi-grids generated by a rotation axis but with two dierent points ⃗ p 1 , ⃗ p 2 one of which is obtained by a translation of ±(⃗ p 1 -⃗ p 2 ) from the other. 4.6 A quadruple of grid points {⃗ q 1 = ⃗ q, ⃗ q 2 = (q x + A x , q y + A y , q z ) ⊤ , ⃗ q 3 = (q x + A x , q y , q z -A z ) ⊤ , ⃗ q 4 = (q x + 2A x , q y + A y , q z -A z ) ⊤ } forms two types of triangles, α and β tiled in a multigrid. . . . . . . . . . . . . . . . 6.1 Two representations of a Pythagorean triple (i 1 , i 2 , λ): Pythagorean angle θ and Pythagorean vector ⃗ v in the plane . . . . . . . . . . . . . . . . . . . 6.2 Construction of all Pythagorean triples of the form (i 1 , i 2 , i 1 + 1) . . . . . 6.3 A convex cone in the plane dened by the pair of vector ⃗ v 1 , ⃗ v 2 . . . . . . . 6.4 Geometrical representations of a Pythagorean quadruple (i 1 , i 2 , i 3 , λ) by the Pythagorean vector (i 1 , i 2 , i 3 ) ⊤ and the pair of angles (θ, ψ). . . . . . . 6.5 A 3D convex cone generated by three vectors ⃗ v 1 , ⃗ v 2 , ⃗ v 3 in the space . . . . 6.6 (a) The 3D convex cone C 3 dened by ⃗ v 1 , ⃗ v 2 , ⃗ v 3 and their projected vectors Rotation is a well known function learned at junior high school. For most of the people, rotation has been a well dened function from R 2 into R 2 for centuries.

⃗ v ′ 1 , ⃗ v ′ 2 , ⃗ v ′ 3 on
Recently, with the spread of computer sciences, more and more people are becoming interested in the rotation as a function from Z 2 into Z 2 , commonly denoted by discrete rotation. Indeed, in computer science, the data acquired from digital cameras or other methods are often digital data, i.e. containing only integers, and it is known that using only integers during computations reduce the computation time and avoid approximations errors.

Nonetheless, performing rotations in Z 2 is not obvious. Many problems arise mainly due to the fact that π, which seems to be inseparable of the rotation function, is a transcendental number.

For the last twenty years, some works have been done to understand the particularity of rotations in the discrete space [START_REF] Hartley | Global optimization through rotation space search[END_REF], to dene and design discrete rotation in 2D [START_REF] Andres | The quasi-shear rotation[END_REF][START_REF] Hénon | Good rotations[END_REF] or discrete rotation in 3D [810]. However, the discrete rotation problem is far to be solved.

The comprehension on discrete rotations is far to be the same as comprehension on Euclidean rotation. The apparently obvious properties of the Euclidean rotation such as bijectivity, transitivity are hard to keep in discrete spaces. Today, at our best knowledge, there is no discrete rotation that have all the properties of the Euclidean rotation.

During my thesis, we try to improve the understanding of the discrete rotation. The results obtained are more theoretical than practical. Our main results showed that it is possible to perform rotations in 2D and in 3D using only integers computations, then obtaining exact results under some assumptions. We also showed that every rotations performed with oating number can be performed using the methods proposed in this manuscript. However, even if these results are exact, the associated algorithms have a hight time complexity.

The discrete rotation is a function coming from a new mathematical eld denoted by digital geometry. In the next sections of this introduction, we make a brief history of the geometry and its evolutions that lead us to the digital geometry and then to the discrete rotations. The following sections of this chapter are not required for the understanding of this thesis and can be skipped, except the last one that gives a quick overview of the content of this manuscript.

A quick history of geometry

Geometry, literally "Earth-Measuring" is a part of mathematics concerned with questions of sizes, shapes, relative positions of gures, and spaces properties. Geometry is one of the oldest sciences. Initially a body of practical knowledge concerning lengths, areas, and volumes, in the 3rd century before Christ geometry was put into an axiomatic form by Euclid [START_REF] Vitrac | Euclide : Les Eléments. Vol. 1 : Introduction générale et Livres I à IV[END_REF]. Euclidean geometry set a standard for many centuries to follow.

The introduction of coordinates by Descartes, during the 17th century, and the concurrent development of algebra, marked a new stage for geometry, since geometric gures, such as plane curves, could now be represented analytically, i.e., with functions and equations. This played a key role in the emergence of calculus in the 17th century. The subject of geometry was further enriched by the study of intrinsic structure of geometric objects that originated with Euler and Gauss and led to the creation of topology and dierential geometry.

Since the discovery of non-Euclidean geometry in the 19th-century, the concept of space has undergone a radical transformation. Contemporary geometry considers manifolds, spaces that are considerably more abstract than the familiar Euclidean space. These spaces may be endowed with additional structures, allowing speaking about lengths.

For fty years, geometry has been required in a lot of applications in computer science such as computer vision, computer graphics, medical image analysis, etc. However, in computer science, the classical geometry is not suitable. First, in computer science, it is only possible to represent exactly integers and rational numbers in memory, whereas

Euclidean geometry mainly uses real numbers. Secondly, data used in many application are represented using only integer numbers. Therefore, using the oating numbers, which are the numbers that approximate real numbers in computer science, inherent to the classical geometry is not desirable. Consequently, a new eld of geometry, the digital geometry, have been introduced a few years ago in order to design a geometry using integers and rational numbers only.

Digital geometry

In recent years, with the spread of the usage of images generated or processed by computers, a new mathematical theory called digital geometry has emerged as a subeld of the discrete geometry. This theory aims at studying objects called discrete objects consisting of countable sets of discrete points, whereas Euclidean continuous objects generally consist of non-denumerable sets of points. Besides their names, discrete objects have little in common with Euclidean geometric objects. Indeed, most elementary results of Euclidean geometry do not hold in discrete spaces: fundamental notions like continuity does not hold (what does "continuous" mean in a discrete space ?), even the basic denition of objects become complex (many denitions of a discrete circle that are not equivalent exists). Digital geometry tries to address these problems by devising new results specic to discrete spaces and objects, and transposing the familiar notions of Euclidean geometry to discrete geometry. an innity since the number of pixels in a discrete circle is necessarily nite. Another example can be two circles that have no intersection in the Euclidean space, but their discretized versions have one or more common pixels.

From digital geometry to discrete rotations

The rotation is a geometrical transformation that have probably been introduced by the Greeks one or two centuries B.C. for studying astrology. The rotation is inseparable from the Archimedes constant, better-known since the 17th or 18th century as π. Archimedes was the rst to rigorously study π by realizing that he can obtain a lower and a upper bound by inscribing circles in regular polygons and calculating the outer and inner polygons respective perimeter [START_REF] Archimedes | Measurement of a circle[END_REF].

Between Archimedes, about ve centuries B.C. and the end of the 19 th century, many mathematicians tried to give a good approximation of π or its exact value. The accuracy of the estimation of π quickly increases with the development of mathematics. With the introduction by Liouville of transcendental number in 1844 [START_REF] Liouville | A propos de l'existence des nombres transcendants[END_REF] and the demonstration by Lindemann in 1882 that π is one of these numbers [START_REF] Louis | Über die Zahl π[END_REF], we know that π cannot be exactly represented or computed.

The impact of the belonging of π to the transcendental numbers on digital geometry is very important. Indeed, we know that each computation from Z n to Z n using π will give an approximated result. More generally, π is not the only problem of the rotations in discrete space. Particular angles exist that allow rotations to be performed without using π denoted by Pythagorean angles [START_REF] Anglin | Using pythagorean triangles to approximate angles[END_REF]. Whereas, performing rotations using these angles will imply dierent problems that will be explained in Chapter 2. Some methods have been developed in digital geometry in order to solve these problems and to perform rotations in discrete space that keep some of the properties of the Euclidean rotations.

These rotations will be introduced in Chapter 2. However, the method that have been proposed cannot solve all the problems of rotations in discrete space. Each one of them have some properties but none have all the properties of Euclidean rotation. Therefore, during my thesis, we tried to improve the understanding of rotations in discrete space.

In the next section, we give a quick overview of the content of this manuscript.

Structure of this manuscript

This manuscript is centered on six chapters. Chapter 1 is nished, so it is useless to introduce it.

Chapter 2 presents a survey on 2D discrete rotations. The rotations in the 2D discrete space have been studied during the last 30 years. However many problems are remaining in the application of rotations in the discrete space. Thus, in Chapter 2, we also introduce and explain these problems. We also present a short study on the rotation in dierent regular grid and show that using other regular grid can give better result than orthogonal grid.

In the Euclidean space, for a given point and center of rotations, there is an innity of dierent rotations of the point around the given center. This is not true in a discrete space as the number of dierent rotations is nite and depends on the distance between the point and the center. Chapter 3 presents a discrete rotation using this property of rotations in the discrete space and then gives a method that estimates the angle of rotation from a given pair of digital images that represent the same object before and after a rotation.

Chapter 4 presents a short survey on rotations in the 3D continuous space since, at our best knowledge, there is not discrete 3D rotation. Rotations in the 3D space can be done in many ways such as composing three 2D rotations or by rotating around an arbitrary axis. In this chapter, we will present the extension of the rotation introduced in Chapter 3 to the 3D discrete space using a rotation around an arbitrary axis.

Chapter 5 extends the second part of Chapter 3. It presents a discrete method to nd from a pair of 3D digital images the angle of rotation and the axis of rotation that transforms the rst digital image into the second digital image where the correspondence of points between these two digital images are supposed to be known.

Chapter 6 presents a study on Pythagorean n-tuples that are the extension to dimension n of the well-known Pythagorean triples. They can be represented in any dimension as a vector. In this chapter, we will show that they are dense in every dimension, i.e. any vector in any dimension can be approximated by an innity of Pythagorean n-tuples.

Since all computations on Pythagorean n-tuples can be done using only integers, they are helpful for designing discrete rotation in any dimension.

We will conclude by some future works in Chapter 7. The rst one will be to design a discrete rotation valid in any dimension. As for angles, there is a nite number of rotation axis in the discrete space. However, there is no clear relation established between these two points. Another future work will be to clarify this relation.

In appendix, we add a glossary that recalls the main concepts and denitions used and introduced in this thesis but in a less formal way.

Chapter 2

Rotations in 2D discrete spaces

Introduction

In this chapter, we present a survey on the discrete rotations and the problems arising while we apply rotations in a discrete space. For many applications, the rotation in the discrete space is a required transformation for image computation in computer science such as image matching, construction of mosaic images [START_REF] Hansen | Real-time scene stabilization and mosaic construction[END_REF] but also in physics such as celestial mechanics [START_REF] Hénon | Good rotations[END_REF]. Typically, the input image is a discretized version of the original analog signal and, therefore, stored as a byte image. A rotation is applied to process the image; this process may result in intermediate oating point numbers. However, the nal image displayed or analyzed is similar to the input image stored as a byte image.

As an example, a 2D digital image is represented as a set of N × N discrete points where each discrete point ⃗ p = (i, j), ∈ {0, . . . , n} 2 contains the color value. If we apply the classic rotation matrix

( cos α -sin α sin α cos α )
on this set of points, we obtain a set of points that belongs to R 2 instead of Z 2 . Therefore, the result of the rotation must be discretized in order to obtain an output comparable to the input. During the discretization, some problems occur such as points without pre-image or points with two pre-images.

Usually, in application of discrete rotations, we consider that the discrete plane is tiled by squares, usually denoted by pixel. However, there exist other regular forms that tile the plane such as triangles and hexagons. Performing rotation in the plane using these two forms can give some interesting results.

In this section, we will present the classic denitions and properties expected for the discrete rotation algorithms. The expected properties are the same than the Euclidean rotation, however, we will also explain why it is dicult to keep these properties in the discrete space. Then we will present a survey on the main existing rotation algorithms.

Finally, we propose a study on rotation in triangle and hexagonal grid.

Denitions and Properties

It is generally admitted that a discrete rotation is a rotation algorithm that transforms a set of discrete points into another set of discrete points. More formally:

Denition 2.1. An algorithm is denoted by a discrete rotation algorithm if it performs a rotation of a set of points from

Z 2 into Z 2 .
This denition allows discrete rotation algorithms to use oating numbers 1 during computations. However, using oating numbers during computations involves errors of approximation in the results. A simple example is the sum: sin 2 α+cos 2 α. In mathematics this sum is always equal to 1 although with the use of oating numbers, we can obtain the result 1 ± ϵ, ϵ > 0. Therefore, we believed that Denition 2.1 is not reliable enough.

Then we propose a more restrictive denition for the discrete rotation: Denition 2.2. An algorithm is denoted by a discrete rotation algorithm if it performs a rotation of a set of points from Z 2 into Z 2 and uses only integers during computations.

A rotation algorithm using only rational numbers during computations is also a discrete rotation algorithm since, as explained in introduction, rational numbers and integers are equivalent. Hereafter, the notation discrete rotation will always refers to Denition 2.2.

Denition 2.2 does not ensure any property on discrete rotation algorithms. In the next section, we will present some discrete rotation algorithms that do not have the same properties than the Euclidean rotation. Hereafter, we will introduce the Euclidean rotation adapted to discrete space: discretized Euclidean rotation, abbreviate DER and explain the most common problems arising while applying DER and associate to each of these problems a property that a discrete rotation algorithm should have. Note that DER is not a discrete rotation regarding to Denition 2.2 since it uses sine and cosine functions.

The rst problem denoted by double points is illustrated in distance between two 4-connected pixels is 1 when the diagonal of a pixel is √ 2. Then after a rotation and before the discretization, the centers of two 4-connected pixels can belong to the same pixel. However, using the same idea, we can conclude that triple points cannot exist.

The second problem denoted by blank points is illustrated in Figure 2.1(b). There is sets of four points that form a 2 × 2 square, which is 4-connected, that after a rotation by a given angle become 8-connected. We call the point within the 8-connected square a blank point. This problem has the same origin than the problem of double points. After the application of DER on an image and if we do not crop the rotated image, the number of double points is the same as the number of blank points, however, this is not true on an arbitrary set of points.

Few examples of these two problems are given in Figures 2. These two problems show that a DER is not necessarily a bijective function since there are points resulting of a discrete rotation that have two pre-images and the inverse rotation does not give the original image. To capture these two problems, we introduce the following denition: Another problem similar to the problem of composition of discrete rotations is the inverse rotation as illustrated in Figure 2.1(d). Let ⃗ p 2 be the result of the rotation by DER of ⃗ p 1 by an angle α. We cannot ensure that the result of the rotation by DER of ⃗ p 2 by an angle -α is ⃗ p 1 .

These two problems show that rotations perform in discrete space are not necessarily transitive. Note that the inverse rotation is just a particular case of the transitivity problem were θ = -α. To capture this problem, we introduce the following denition: We say that a rotation is DER-like if the object obtained using the red way is the same than the object obtained using the blue way.

At our best knowledge there is no discrete rotation that is bijective, transitive and DERlike for any angle. We strongly believe that is it not possible to design such a rotation.

Therefore, we propose the following conjecture:

Conjecture. 1. If the function R : Z → Z, is a discrete rotation, this function is not bijective and transitive and DER-like for any angle.

In the next section, we present some discrete rotation algorithms. We will use the three denitions presented above to characterize these discrete rotations.

Figure 2.5: If a discrete rotation is DER-like, then the red way that discretize then rotate will gives the same result than the blue way that rotate then discretize.

Survey on discrete rotations in the plane

In this section, we present some of the existing discrete rotation according to Denition 2.2. The rst rotation is the rotation by discrete circle, then we present the rotation by Pythagorean line already addressed in previous section. We then present the shear rotation that consist in the decomposition of the rotation into three sheer transforms.

Finally, we present the rotation by hinge angles which is the rotation that we mainly studied in this thesis.

Rotation by discrete circle

We rst present the rotation by discrete circles [START_REF] Andres | Cercles discrets et rotations discrètes[END_REF]. A discrete circle of radius r is the

set of n pixels ⃗ p 1 = (p 1x , p 1y ), ⃗ p 2 = (p 2x , p 2y ) . . . , ⃗ p n = (p nx , p ny ) such that (r -1 2 ) 2 ≤ p 2 ix + p 2 iy < (r + 1 2 ) 2 , i = 1 . . . n 3 .
Using this denition of discrete circles, we know that we will obtain a tiling of the plane by discrete circles. In other words, we can ensure that each pixel belongs to a discrete circle and to only one. The main idea of the rotation by discrete circles is to move the pixel in each discrete circle of a number depending on the rotation angle. This number is obtained using the following equation :⌊ α×n 360 ⌋ where α is the given rotation angle and n the number of pixels belonging to the discrete circle. An example of the rotation by discrete circles is given in Figure 2.6. In his thesis, Andres shows that the rotation by discrete circles is transitive and bijective. We can easily be convinced of these two properties as we just displaced pixel in the circle during 3 There are dierent denitions of a discrete such as Bresenham circle [START_REF] Bresenham | A linear algorithm for incremental digital display of circular arcs[END_REF], but this rotation requires a denition of discrete circles that tile the plane. the rotation of a discrete circle. Therefore, there is no blank or double points. In other words, this discrete rotation does not lose any pixel.

The application of DER and the rotation by discrete circles on the same set of points will give dierent results. This remark seems obvious since the rotation by discrete circles is bijective. However, Figure 2.6, which gives an example of a rotation by discrete sphere on an angle of 45 • , shows some shears in the resulting digital image. These shears do not arise during a rotation by DER. Thus, we can conclude that according to Denition 2.5, the rotation by discrete circles is not DER-like.

Rotation by Pythagorean lines

Another discrete rotation is the rotation by Pythagorean lines. A line of slope

i 1 i 2 that veries i 2 1 + i 2 2 = λ 2 is a Pythagorean line if i 1 , i 2 , λ ∈ Z.
In [START_REF] Andres | Cercles discrets et rotations discrètes[END_REF], the author introduce a rotation algorithm using the particular set of Pythagorean lines where λ = i 2 + 1. In [START_REF] Jacob | On discrete rotations[END_REF], authors give the proof that the rotations induced by these set of Pythagorean lines is bijective, transitive and DER-like. However, the rotation by Pythagorean lines has restricted angles. The biggest angle is given by the Pythagorean triple [START_REF] Nagy | An algorithm to nd the number of the digitizations of discs with a xed radius[END_REF][START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF][START_REF] Voss | Discrete Images, Objects and Functions in Z n[END_REF] Note that the results obtained by DER illustrated in Figure 2.3(c) are the same than the results obtained by a rotation by the Pythagorean line associated to [START_REF] Nagy | An algorithm to nd the number of the digitizations of discs with a xed radius[END_REF][START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF][START_REF] Voss | Discrete Images, Objects and Functions in Z n[END_REF]. It is possible to compose rotations by Pythagorean lines in order to obtain any angle, but in case of composition, the nal result will not be consistent since the connectivity between pixels will be lost during the composition 4 . In the nal result, while using composition, some shear eects will appear, which do not arise during a rotation by DER. The twelve biggest angles that can possibly be obtained by Pythagorean lines

Shear rotation and quasi shear rotation

It is also possible to design a discrete rotation using three translations. Indeed, the rotation matrix

( cos α -sin α sin α cos α
) can be rewritten as the product of the three matrices :

( 1 -tan α 2 0 1 ) , ( 1
0 sin α 1 ) , ( 1
-tan α 2 0 1
) that are translation matrix. This method called shear rotation is well studied in the continuous plane [START_REF] Paeth | A fast algorithm for general raster rotation[END_REF]. In [START_REF] Andres | The quasi-shear rotation[END_REF], the author gives a method, called quasi-shear rotation, to discretize these three matrices and then obtain a discrete rotation, working in the discrete space, using three horizontal or vertical shear.

Obviously, the three translations are bijective, then we can conclude that the quasi-shearrotation is bijective. However, the translations are not transitive, indeed the composition of rotations of angles θ and -θ does not necessarily gives the original image and the composition of the two rotations of angles θ 1 and θ 2 does not necessarily gives the same result as the rotation of angle θ 1 + θ 2 . Moreover, this rotation is not DER-like since some shears will appear during the rotation. We can conclude that QSR is bijective, transitive and non DER-like. An example of the quasi-shear rotation where the three translations are decomposed is given in Figure 2.8.

Rotation by hinge angles

More recently, a new discrete rotation has been designed by Nouvel in [START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF] and improved by us Chapter 3. This rotation is based on hinge angles and will be more developed in the next sections. Roughly speaking, an angle is a hinge angle for a discrete point if the result of this point's rotation by this angle belongs to the half-grid. In other words, if the angle of rotation is slightly greater than a hinge angle, the result of the rotation of the associated point will be in a pixel p, if the angle of rotation is slightly lower than a hinge angle the result of the rotation of the same point will be in a pixel that is 4-connected with p. As far as we know, the discrete rotation by hinge angles is the only discrete rotation DER-like for any angle, however, rotation by hinge angles is neither bijective nor transitive.

In our purpose, we assume that discrete data are the result of the acquisition by a digital camera of continuous objects. This supposition enforces the property that the discrete rotation between two dierent sets of discrete data gives the same result as DER. For this reason, during my thesis, we have chosen to improve the rotation by hinge angles.

Therefore, in the next chapter, we will dene hinge angles and give their properties.

Rotation in other grids

All the rotations presented in Section 2.3 are based on the square grid. However, other regular grids exist, in this section, we present a short study on rotations in the two other kinds of regular grids that are triangles and hexagonal grids. In discrete geometry, the square grid is the most commonly used, however, this grid is not the only regular grid that tiles the plane. There are two other regular geometrical forms that tile the plane:

the triangle and the hexagonal.

In the eld of discrete geometry, the square grid is almost the only grid used to perform transformations. During my thesis we exclusively used the square grid. However, using other kind of grids can present interesting properties such as conservation of the neighborhood, that helps to keep a visual consistency of the digital image, or loosing (a) (b)

Figure 2.9: A pixel and (a) its four 4-connected pixels, (b) its eight 8-connected pixels in a square grid.

less points during rotations, that allows keeping more information on the transformed digital image. Thus in this subsection, we present some experiments on rotations done using these two other kinds of grids. Here we only consider the two other regular grids:

the triangles and hexagons, that tile the plane with an area of 1 and with regular edges.

First, we introduce some notions required to understand the studies done in this chapter, then, we present the methodology used to perform our test and nally, we present our results on the rotation on triangle grid then on hexagonal grid.

Connectivity on the discrete grids

In this section, we present the notion of connectivity in discrete grids.

The notion of neighborhood in the continuous and in the discrete space is totally dierent.

In the continuous space a neighborhood is the set of all the points whose distance from a given point is shorter than a given value. This set contains an innity of points if the given value is not null. However, this is not the case in the discrete space as the neighborhood of a point only contains a nite number of points if the given value is nite.

Moreover, for any discrete point, we can always nd at least one point so that there is no other points between them. In discrete geometry, instead of talking of neighborhood between two points, we use the connection between pixels/voxels. Usually, we say that two pixels/voxels are m-connected. The variable m depends on the dimension of the discrete space and the grid that we are considering. Considering the 2D discrete space tile by squares, the variable m can take the values of 4 or 8 as illustrated in Figures 2.9 A pixel and its six 6-connected pixels in a hexagonal grid an edge and they are 26-connected if they share a face, an edge or a vertex. Note that in 2D, there exist three shape of regular polygons that tile the plane. However, in 3D, the only regular polygon that tile the space is the cube.

The notions of connectivity explained in this section is required for the study presented hereafter.

Methodology

In this section, we analyze two aspects of the discrete rotations: the number of points lost during rotations and the conservation of the neighborhood after rotations.

The rst step to study the number of points lost during a rotation on a non square grid is to establish a bijection between the square grid and the other regular grid. The aim is to transform a digital image into a set of points belonging to the regular grid or to transform a set of points of the regular grid into a digital image. Figures 2.12 and 2.13 illustrate the bijections between the square, triangle and hexagonal grid. By using the bijection between the two grids, we convert a digital image of size 1000 × 1000 containing white points only into a set of 10 6 oating points centered on a polygon of the regular grid.

Then we apply the standard matrix rotation

( cos α -sin α sin α cos α
) on these set of points and discretize them by using a function, which depends on the regular grid used, that returns a point centered on a regular grid's polygon. We then apply the bijection function that transforms the regular grid into a squared grid, count the number of remaining white points and compare that number to the original white points number.

To study the conservation of the neighborhood after a rotation, we rst generate a set of points containing all the points of a square grid. Then, we generate a table that contains for each points four or eight corresponding neighbors. The next step converts the set of points into an-other regular grid, applies the rotation transformation on it and converts the result into a square grid. Finally, from the obtained results, we generate a second the 8-connectivity after a rotation in every kind of regular grids because the comparison of remaining connectivity between a 4/8-connected grid and a 6 connected grid is not relevant.

Experimental result

Results on triangle grid We can conclude that rotation on triangle grid has no interest regarding the number of points lost or in the connectivity's preservation.

Results on Hexagonal grid

Figure 2.16 presents the results obtained on rotations using the hexagonal grid. We compare them with the result obtained in Figure 2.3. First, Figure 2.16 seems to be the contrary of Figure 2.14, this impression is probably due to the fact that the triangle grid is the dual of the hexagonal grid.

We can see that the rotations on hexagonal grid rarely lose more than 9, 5% of points.

This result is really better than the results obtained on square grid since on the square grid we lost more that 15% of points between 33 • and 57 • . However, for angles lower than We can see that the results obtained on neighborhood using hexagonal grid are globally worse than the results on square grid. This is probably due to the fact that the hexagonal grid counts only six neighbors instead of four or height. However, for the angle between 33 • and 58 • for the 4-connectivity and between 37 • and 54 • for the 8-connectivity, rotations on hexagonal grid give slightly better results.

Finally, regarding the number of points lost during rotations, we see that rotations on hexagonal grid globally give better results, than the two other regular grids. We also see, that for 4 and 8-connectivity, results around the angle of 45 • are better. Thus, we can conclude that rotations on hexagonal grid globally give better results than rotations on square grid.

Visual result

In this paragraph, we present the result of rotations for three kind of regular grids. They illustrate the results given above that are theoretical results that cannot give a good idea to the problems of double points and connectivity encountered. These pictures are the results of the composition of four rotations of a given angle for each regular grid.

Figure 2.17 shows that for an angle of 15 • , the rotation on square grid keep a little bit more points than the rotation on hexagonal grid. The curve and line bone between

French lands and the sea seem to be better preserved by the rotation on square grid.

Figures 2.18 and 2.19 shows that rotations on hexagonal grid keep more points than rotations on other grids. However, in both gures, we see that the lines and the curves are better preserved by rotation on square grid.

We can conclude that rotations on hexagonal grids are useful to keep more points during rotations, however rotations on square grids are better for visual consistency.

Conclusion

In this Chapter, we have introduced the problematic of rotations in the discrete plane.

We have shown that existing discrete rotations cannot have all the properties of the Chapter 3

2D rotations in a discrete space using hinge angles

Introduction

Rotations in the discrete plane are required in many applications for image computation such as image matching, construction of mosaic images [START_REF] Hansen | Real-time scene stabilization and mosaic construction[END_REF]. At the moment, the most popular method to estimate the rotation angle is to approximate the rotation matrix by minimizing errors [START_REF] Hansen | Real-time scene stabilization and mosaic construction[END_REF].

In the continuous plane, the Euclidean rotation is well-dened and possesses the property of bijectivity. This implies that for two angles γ 1 , γ 2 and a set of points A, if the Euclidean rotation with angle γ 1 applied to A gives the same result as the Euclidean rotation with angle γ 2 applied to A, then we have

γ 1 = γ 2 .
In the discrete plane, however, two points resulting of two dierent Euclidean rotations (with two dierent angles) of the same grid point can be discretized in the same pixel, i.e. discretized as the same grid point. For this reason, two dierent angles can give the same rotation result for a set A of grid points 1 . In other words, we can dene a set of admissible rotation angles S so that any angle in S gives the same rotation result for A. Note that S depends on A. The two most interesting angles in S are the lower and the upper bounds. Indeed, with only these two angles we can deduce the other angles in S. To identify the exact bounds, we should not involve any computation error. Thus, we work with discrete geometry tools, which guarantee computation avoidance with real numbers. Moreover, because we assume that our data are discretized from continuous images of an object, we enforce the property that the discrete rotation2 between two dierent sets of grid points gives the same result as DER.

Some works on discrete rotations already exist. The most widely used discrete rotation at the beginning is the CORDIC algorithm [START_REF] Jack | The cordic trigonometric computing technique[END_REF]. The CORDIC algorithm uses a sequence of fractions for an approximation of π. It multiplies/adds fractions in this sequence to approximate values of sine and cosine. Thus multiple approximations lead to little dierences between results of DER and those of CORDIC. Note that CORDIC algorithm is not a discrete rotation according to Denition 2.2. As presented in Chapter 2, Andres described in [START_REF] Andres | Cercles discrets et rotations discrètes[END_REF][START_REF] Andres | The quasi-shear rotation[END_REF] some discrete rotations such as the rotation by discrete circles, the rotation by Pythagorean lines or the quasi-shear rotation. Computation executed during these rotations are exact. However, we have shown in Chapter 2 that these rotation algorithms do not give the same results as DER.

On the other hand, Nouvel and Rémila proposed in [START_REF] Nouvel | Incremental and transitive discrete rotations[END_REF] another discrete rotation based on hinge angles, which gives the same results as DER. It is known that a sequence of hinge angles is a set of particular angles determined by a digital image in the sense that any angle between two consecutive hinge angles gives the identical rotated digital image. This means that hinge angles correspond to the discontinuity of DER. Nouvel and Rémila showed that each hinge angle is represented by an integer triplet, so that any discrete rotation of a digital image is realized only with integer calculation. Because their algorithm gives the same rotation results as DER, we see that hinge angles represented by integer triplets give sucient information for executing any digital image rotation.

In this chapter, we introduce the hinge angles and give their properties. Then we present a method to perform a discrete rotation using hinge angles. Finally, we propose a discrete method to nd the lower and upper bounds of admissible rotation angles. Our method uses hinge angles because this algorithm is DER-like and because they allow exact computations. The input data of our method is two sets of grid points where point correspondences across the two sets are known. The output is two hinge angles that give the lower and the upper bounds of the admissible rotation angles for the two sets.

Hinge angles

We consider grid points in Z 2 as the centers of pixels and rotate them in such a way that the rotation center has integer coordinates. Hinge angles are particular angles that make some points in Z 2 rotated to points at the frontier between adjacent pixels. In this section, we give the denition of hinge angles and their properties related to Pythagorean angles.

Denition of hinge angles

Let ⃗ p = (x, y) be a point in R Because H can be seen as the discontinuity of the rounding function, hinge angles can be regarded as the discontinuity of the discretized Euclidean rotation. In other words, hinge angles determine a transit of a grid point from a pixel to its adjacent pixel during the rotation.

The following theorem is important because it shows that we can represent every hinge angle with three integers.

Theorem 3.2 ([21]

). An angle α is a hinge angle for a grid point (P, Q) ∈ Z 2 if and only if there exists K ∈ Z such that

2Q cos α + 2P sin α = 2K + 1. (3.1)
Geometrically, a hinge angle α is formed by two rays going through (P, Q) and a halfgrid point (K + 1 2 , λ) where the two rays share the origin as their endpoints as shown in Figure 3.1 (left). This theorem indicates that all calculations related to hinge angles can be done only with integers. Hereafter, α indicates a hinge angle.

We denote by α(P, Q, K) the hinge angle generated by an integer triplet (P, Q, K).

Setting λ = √ P 2 + Q 2 -(K + 1
2 ) 2 , we easily derive the following equations from (3.1) and Figure 3.1 (left), Note that we have a case where a half-grid point is (λ, K + 1 2 ) instead of (K + 1 2 , λ). In such a case, the above equations become To distinguish the case of (K + 1 2 , λ) from that of (λ, K + 1

cos α = P λ + Q(K + 1 2 ) P 2 + Q 2 , sin α = P (K + 1 2 ) -Qλ P 2 + Q 2 . ( 3
cos α = Qλ + P (K + 1 2 ) P 2 + Q 2 , sin α = P λ -Q(K + 1 2 ) P 2 + Q 2 . ( 3
2 ), we change the sign of K; we use α(P, Q, K) for the case of (K + 1 2 , λ), and α(P, Q, -K) for the case of (λ, K + 1 2 ). Because the symmetries allow us to restrict α to the range [0, π 2 ], as mentioned above, we may assume that K is positive.

Properties related to Pythagorean angle

Because hinge angles are strongly related to Pythagorean angles, certain properties of Pythagorean angles are required to prove some properties of hinge angles. Thus, we rst give the denition of Pythagorean angles and their properties.

|c| = max{| a |, | b |, | c |}. If gcd(a, b, c) = 1, then c is odd. If gcd(a, b, c) = i, then gcd( a i , b i , c i ) = 1 and the triplet of integers ( a i , b i , c i )
generates the same Pythagorean angle as (a, b, c).

Theorem 3.6. Let E h be the set of hinge angles and E p be the set of Pythagorean angles.

Then we have E

h ∩ E p = ∅.
Proof. 1. Assume that there is an angle

α such that α ∈ E h and α ∈ E p . Since α ∈ E p ,
we can nd an integer triplet (a, b, c) generating α where gcd(a, b, c) = 1. By substitution of (3.4) in (3.1), we obtain This theorem shows that it is not possible to rotate a point (i, j) ∈ Z 2 to a point (x, y)

2 Qa + P b c = 2K + 1, (3.5
such as x = i + 1 2 , y = j + 1 2
, if the angle of the rotation is a hinge angle.

The next theorem shows an interesting relation between hinge angles and Pythagorean angles.

Theorem 3.7 ([4]). Let θ be a Pythagorean angle and α be a hinge angle. The angle α ′ = α + θ is a hinge angle.

Rotation by hinge angles

As explained in Chapter 2, as our work is focused on the discrete space, the rotation of a grid point by two dierent angles can give the same result. Namely, two dierent angles give the same result after the rotation of a grid point followed by discretization.

Generally, there is a range of angles in which the same result is obtained. We thus dene admissible rotation angles, abbreviated hereafter by ARA, to represent this range of angles.

In this section, we propose a method for computing the lower bound α inf of ARA for a given digital image from a given angle. Note that with minor modications, this method can also nd the upper bound α sup of ARA. Thus applying any rotation to the given Algorithm 1 Function for the lower bound rotation angle for a point.

Require: A point ⃗ p(P, Q), a Pythagorean angle θ

Ensure: A hinge angle α(P, Q, K)

K max ← ⌊ √ P 2 + Q 2 -1⌋; K min ← 0; K ← ⌊ Kmax+K min 2 ⌋; while K max -K min ̸ = 1 do if α(P, Q, K) > θ then K max = K; else K min = K; end if K = ⌊ Kmax+K min 2 ⌋;
end while return α(P, Q, K);

digital image with an angle between α inf and α sup gives the same result. We note that both α inf and α sup are hinge angles.

Our input is a Euclidean angle. However, we can replace it by a Pythagorean angle as there is a method with linear time complexity O(m) to approximate a given Euclidean angle with a Pythagorean angle with a precision of 1 10 m [15], where m is a xed integer that represents the quality of approximation. Below, we assume that a Pythagorean angle is given as in [START_REF] Nouvel | Incremental and transitive discrete rotations[END_REF].

Nouvel and Rémila presented a method for computing all possible hinge angles for a grid point or a pixel in a digital image [START_REF] Nouvel | Incremental and transitive discrete rotations[END_REF]. Their method can be used for nding the hinge angle that is the lower bound of the admissible rotation angles. Its time complexity is O(n log(n)) where n is the number of all hinge angles for a given grid point. Note that n depends on the coordinates of the grid point.

Computing the lower bound rotation angle for grid point

For each grid point [START_REF] Nouvel | Incremental and transitive discrete rotations[END_REF]. We can compare in magnitude any pair of hinge angles. This means that we have a totally ordered set {α(P, Q, K 1 ), α(P, Q, K 2 ), ..., α(P, Q, K n )} of hinge angles in the ascending order where K i ∈ Z. Given a Pythagorean angle θ, in order to nd the lower bound rotation angle α(P, Q, K i ) such that α(P, Q, K i ) < θ < α(P, Q, K i+1 ), we use a binary search. The binary search allows us to nd α(P, Q, K i ) in O(log(n)), providing that we can compare a hinge angle with a Pythagorean angle in a constant time. The algorithm is described in Function 1. Note that thanks to Theorem 3.9 a Pythagorean angle for the input can be replaced by a hinge angle.

⃗ p = (P, Q) ∈ Z 2 , there are less than n = ⌊ √ P 2 + Q 2 + 1 2 ⌋ dierent hinge angles in each quadrant
The following theorem shows that the comparison between a hinge angle and a Pythagorean angle is executed in a constant time.

Theorem 3.8. Let α be a hinge angle and θ be a Pythagorean angle. We can check whether α > θ in a constant time with integer calculation.

Proof. 2. Let α(P, Q, K) be a hinge angle in [0, π 2 ] and θ(a, b, c) be a Pythagorean angle in [0, π 2 ]. From (3.2) and (3.4), we obtain

cos α -cos θ = Q(K + 1 2 ) + P λ P 2 + Q 2 - a c .
If θ is greater than α, cos α -cos θ > 0. Thus

cQ(2K + 1) -2a(P 2 + Q 2 ) > -2cP λ. (3.6)
Since we know that c, P, λ are positive, the right-hand side of (3.6) is always negative.

Thus, if the left-hand side of (3.6) is not negative, then θ > α. Otherwise, we take the square of (3.6), so that we only have to check whether the following inequality holds:

[ cQ(2K + 1) -2a(P 2 + Q 2 ) ] 2 < 4c 2 P 2 λ 2 . ( 3.7)
Note that because λ =

√ P 2 + Q 2 -(K + 1
2 ) 2 , we see that 4λ 2 in the right-hand side of (3.7) contains only integer values. Therefore, we can verify (3.7) with integer calculation.

If it is true, θ > α; otherwise, α > θ. Note that because of Theorem 3.6, it is impossible to obtain θ = α.

We claim that this comparison of a hinge angle and a Pythagorean angle is executed in a constant time because even in the worst case, we only have to check two equations (3.6) and (3.7).

We mention the importance of the rotation with angle π 2 and its multiplications. In fact, if the angle of a rotation is equal to π 2 , π, 3π 2 , we just have to ip x and/or y-coordinates by changing their signs. It gives the justication that we can restrict the input angle θ to 0 < θ < π 2 .

Computing the lower bound rotation angle for a set of grid points

In this subsection, we present an algorithm for computing the lower bound rotation angle from a given Pythagorean angle θ for a digital image consisting of m grid points A. The output is a triplet of integers that represents the lower bound rotation angle for A. We note that the lower bound rotation is a hinge angle.

The algorithm computes all hinge angles for all points in A, and sorts them to keep the largest one. More precisely, we rst compute the lower bound rotation angle for the rst point of A, and store it as a reference. Then, we compute the lower bound rotation angle for the second point in A and compare it with the reference to keep the larger one. After repeating this procedure for all points in A, our algorithm returns the lower bound rotation angle α such that α < θ. The time complexity of this algorithm is O(m log(n)) because we call m times the binary search (Function 1) whose time com- plexity is O(log(n)). Function 2 illustrates our algorithm. As shown in Theorem 3.9, the comparison between two hinge angles is realized in a constant time, so that our algorithm does not change the global complexity.

Theorem 3.9. Let α 1 , α 2 be two hinge angles. We can check whether α 1 > α 2 in a constant time and with integer calculation.

Proof. 3. Let α 1 (p, q, k) and α 2 (r, s, l) be two hinge angles in [0,

π 2 ]. From (3.3) we obtain cos α 1 -cos α 2 = p(k + 1 2 ) + qλ 1 p 2 + q 2 - r(l + 1 2 ) + sλ 2 r 2 + s 2 . (3.8) If α 2 is greater than α 1 , cos α 1 -cos α 2 > 0. Thus (r 2 + s 2 )p(2k + 1) -(p 2 + q 2 )r(2l + 1) > 2(p 2 + q 2 )sλ 2 -2(r 2 + s 2 )qλ 1 . ( 3.9)
If the left-hand side of (3.9) is negative and the right-hand side of (3.9) is positive, then

α 1 > α 2 .
If the left-hand side of (3.9) is positive and the right-hand side of (3.9) is negative, then α 2 > α 1 . We can easily check the signs of the left-hand side and the right-hand side of (3.9) with integer computation. Note that p, q, k, r, s, l are all positive, and that (2(p 2 + q 2 )sλ 2 ) 2 and (2(r 2 + s 2 )qλ 1 ) 2 contain only integer values.

If the signs of the left-hand side and the right-hand side of (3.9) are the same, we rst compute the square of each side and then compare the values to identify which is the greater. For simplicity, we assume that the signs of both sides of (3.9) are positive, and

let A = (r 2 + s 2 )p(2k + 1), B = (p 2 + q 2 )r(2l + 1), C = (r 2 + s 2 )q and D = (p 2 + q 2 )s.

Now (3.9) is rewritten by

A -B > 2Dλ 2 -2Cλ 1 . (3.10)
Then we take the square of equation (3.10) to obtain If the sign of the left-hand side of (3.11) is positive, we can deduce that α 2 > α 1 .

(A -B) 2 -4(C 2 λ 2 1 + D 2 λ 2 2 ) > -8CDλ 1 λ 2 . ( 3
Otherwise, taking the square of each side gives us

[ (A -B) 2 -4(C 2 λ 2 1 + D 2 λ 2 2 ) ] 2 < 64C 2 D 2 λ 2 1 λ 2 2 .
(3.12)

We note that we can easily verify whether (3.12) is satised with integer computation alone. If (3.12) is true, α 1 < α 2 ; otherwise α 2 < α 1 . The same logic can be applied to the case where the signs of both sides of (3.9) are negative.

This comparison of a pair of hinge angles is executed in a constant time because in the worst case, we only have to check three equations (3.9), (3.11) and (3.12).

Digital image rotation by a lower bound rotation angle

This section uses the results obtained in Section 3.3 to present an algorithm for rotating a set of points with a given lower bound rotation angle.

It is already proved in [START_REF] Nouvel | Incremental and transitive discrete rotations[END_REF] that we can obtain the same result as the DER with respect to the original rotation angle. Note that our input is a lower bound rotation angle and the input of the algorithm presented in Function 1 is a Pythagorean angle. In spite of this dierence, we can apply the same algorithm thanks to Theorem 3.9, since we are end for return A ′ ; looking for K, which gives the arriving pixel (K, ⌊λ + 1 2 ⌋), for each pixel (P, Q). The algorithm is presented in Function 3. It supposes that the center of rotation is the origin.

For each point, our algorithm calls the binary search (Function 1) to nd the corresponding hinge angle, which designates its new position. If we consider n as the biggest coordinate of all points in A, we can assume that there are less than 4n 2 points in A.

Thus we can conclude that the complexity of our algorithm is O(n 2 log(n)). The rst advantage of our method is that it does not require any oating number calculation. The second advantage is that the exact rotation of the set of points is obtained with only an integer triplet. We need neither matrices nor angles to realize the rotation.

Obtaining admissible rotation angles from two digital images

Let us assume that a set of grid points in the rst image and its corresponding set in the second image are given: A = {⃗ p 1 , ⃗ p 2 , ..., ⃗ p l } and B = {⃗ q 1 , ⃗ q 2 , ..., ⃗ q l } are given where ⃗ p i corresponds to ⃗ q i . Given A and B, we obtain a hinge angle pair {α inf , α sup }. This pair of hinge angles is the lower and the upper bounds of the ARA. Therefore each γ such that α inf ≤ γ < α sup is consistent with the point correspondences between A and B. Hereafter, we assume that A is the original point set and B is the rotated point set by angle γ. In this section, we show how to obtain the ARA from A and B.

Setting Rotation Centers

For any rotation, we need to set a rotation center. Without a loss of generality, we may choose any grid point in a digital image for the rotation center. Indeed, if we apply two rotations on the same image with the same angle but two distinct centers of rotation, the dierence between the two result is captured by a translation. Assuming that rotation centers for A and B are ⃗ p 1 and ⃗ q 1 3 respectively, we dene two translation functions T A and T B such that

T A (⃗ p i ) = ⃗ p i -⃗ p 1 , T B (⃗ q i ) = ⃗ q i -⃗ q 1 , for all ⃗ p i ∈ A, ⃗ q i ∈ B.
We can regard the origin as the rotation centers after these translations. Hereafter, we assume that these translations have already been applied in order to obtain A = {⃗ p 1 , ⃗ p 2 , ..., ⃗ p l } and B = {⃗ q 1 , ⃗ q 2 , ..., ⃗ q l }. The corners of H (⃗ q), namely, four corners of a pixel around ⃗ q.

Computing lower and upper bounds of rotation angles from two corresponding point pairs

In this subsection, we consider a special case of two corresponding point pairs.

We let

A = {⃗ p 1 , ⃗ p 2 } and B = {⃗ q 1 , ⃗ q 2 } where ⃗ p i = (P i , Q i ) and ⃗ q i = (R i , S i ). We then dene a circle C (⃗ p 2 ) going through ⃗ p 2 whose center is ⃗ p 1 . Thus the radius of C (⃗ p 2 ) is r = d(⃗ p 1 , ⃗ p 2 ) where d(⃗ p 1 , ⃗ p 2 ) is the Euclidean distance between ⃗ p 1 and ⃗ p 2 .
Let us dene the half-grid H (⃗ q 2 ) around ⃗ q 2 :

H (⃗ q 2 ) = {(x, y) ∈ H : S 2 - 1 2 ≤ y ≤ S 2 + 1 2 if x = R 2 ± 1 2 , R 2 - 1 2 ≤ x ≤ R 2 + 1 2 if y = S 2 ± 1 2 }.
We set ⃗ p 1 and ⃗ q 1 to be the rotation centers. Then, we need to detect intersections between C (⃗ p 2 ) and H (⃗ q 2 ) in order to nd a hinge angle pair. We thus investigate which corners of H (⃗ q 2 ) are inside of C (⃗ p 2 ).

Setting four corners of

H (⃗ q 2 ) to be C 1 (⃗ q 2 ) = (R 2 -1 2 , S 2 -1 2 ), C 2 (⃗ q 2 ) = (R 2 -1 2 , S 2 + 1 2 ), C 3 (⃗ q 2 ) = (R 2 + 1 2 , S 2 + 1 2 ), C 4 (⃗ q 2 ) = (R 2 + 1 2 , S 2 -1 2
) as shown in Figure 3.2, we dene a binary function F :

F (C i (⃗ q 2 )) =    1 if C i (⃗ q 2 ) is inside of C (⃗ p 2 ), 0 otherwise. In order to obtain F (C i ( ⃗ q 2 )) with integer calculation, we compare each of ∥(2(R 2 ± 1 2 ), 2(S 2 ± 1 2
))∥ 2 with (2r) 2 . Note that we may assume that C (⃗ p 2 ) and H (⃗ q 2 ) always intersect with each other. This is because no intersection between C (⃗ p 2 ) and H (⃗ q 2 ) indicates that ⃗ p 2 and ⃗ q 2 are not corresponding.

The following lemmas are needed to prove Theorem 3.12 below.

Lemma 3.10. For a circle

C (⃗ p 2 ) centered on ⃗ p 1 , any C i (⃗ q 2 ) = (R 2 ± 1 2 , S 2 ± 1 2 ) cannot be on C (⃗ p 2 ) for i ∈ {1, 2, 3, 4} where R 2 , S 2 ∈ Z. Proof. 4. Let r be the radius of C (⃗ p 2 ). Because C (⃗ p 2 ) goes through ⃗ p 2 , r 2 ∈ Z. Let us assume that C i (⃗ q 2 ) is on C (⃗ p 2 ).
Thus the Euclidean distance between the origin and

C i (⃗ q 2 ) is equal to r. This indicates that r 2 = (R 2 ± 1 2 ) 2 + (S 2 ± 1 2 ) 2 .
However, this contradicts the above fact that r 2 ∈ Z.

Lemma 3.11. Let D be a line that belongs to

H . If C (⃗ p 2 ) is a circle centered on ⃗ p 1 .
Then, the number of distinct intersections of C (⃗ p 2 ) and D is two or zero.

Proof. 5. Let ⃗ p 2 = (P 2 , Q 2 ) ∈ Z 2 and the equation representing D be x = i + 1 2 where i ∈ Z. Letting (x 1 , y 1 ) be the coordinates of the intersecting point of D and C (⃗ p 2 ). Then

we have    x 2 1 + y 2 1 = P 2 2 + Q 2 2 , x 1 = i + 1 2 .
From these equations, we obtain

y 2 1 = P 2 2 + Q 2 2 -(i + 1 2 ) 2 . Because P 2 2 + Q 2 2 -(i + 1 2 ) 2
does not belong to Z, y 1 cannot be equal to 0. Therefore there are two distinct solutions for y

1 if P 2 2 + Q 2 2 -(i + 1
2 ) 2 > 0; no solution otherwise. Similar discussion can be done for the case of y 1 = j + 1 2 where j ∈ Z. Theorem 3.12. If two points ⃗ p 2 and ⃗ q 2 are corresponding, the circle C (⃗ p 2 ) and the half-grid H (⃗ q 2 ) always have two distinct intersections.

Proof. [START_REF] Andres | The quasi-shear rotation[END_REF]. In general, if a circle intersects with a square and the center of the circle is not inside the square, we have 1, 2 or 4 intersections. Having just one intersection means that the circle goes through one corner of the square or the circle is tangential to a half-grid. Lemma 3.10 shows that the circle cannot go through one corner. Lemma 3.11 shows that the circle cannot be tangential to any half-grid.

Let C be a circle centered on the origin and going through the point

⃗ p = (p x , p y ) ⊤ ∈ Z 2 .
We assume that there is a pixel with four intersections with C. If such a pixel exists, one of its coordinates must be equal to zero and the intersections are symmetric regarding to the x-axis or the y-axis (according to the null coordinate). We assume that ⃗

p ′ = (0, p ′ y ) ⊤ ∈ Z 2
is the center of the pixel. Then the four points of intersection 1 2 ]. Because these four points belong to C, we have

⃗ p 1 , ⃗ p 2 , ⃗ p 3 , ⃗ p 4 have the coordinates ⃗ p 1 = (δ 1 , p ′ y +ϵ 1 ) ⊤ , ⃗ p 2 = (δ 2 , p ′ y + 1 2 ) ⊤ , ⃗ p 3 = (-δ 2 , p ′ y + 1 2 ) ⊤ , ⃗ p 1 = (-δ 1 , p ′ y +ϵ 1 ) ⊤ where δ 1 , δ 2 , ϵ 1 ∈ [0,
p 2 x + p 2 y = δ 2 2 + (p ′ y + 1 2 ) 2 . (3.13)
By denition, p 2 x , p 2 y , p ′2 y are integers, thus δ 2 + 1 4 must be an integer that is impossible since δ 2 ∈ [0, 1 2 ]. We can deduce that ⃗ p 2 and ⃗ p 3 cannot exist.

Therefore we have only 2 intersections.

(a) ΣF (C i (⃗ q 2 )) = 0 (b) ΣF (C i (⃗ q 2 )) = 1 (c) ΣF (C i (⃗ q 2 )) = 2 (d) ΣF (C i (⃗ q 2 )) = 2 Figure 3.3: Illustration of cases ΣF (C i (⃗ q 2 )) = 0,1,2 or 3.
From Theorem 3.12, we always have two distinct intersections between C (⃗ p 2 ) and H (⃗ q 2 ), and we see that there are four cases corresponding to dierent possibilities to have 0,1,2 or 3 corners inside of C (⃗ p 2 ), as illustrated in Figure 3.3.

Remark 3.13. The case represented in (a) in Figure 3.3 never happen unless two conditions are satised.

1. ⃗ q 2 is on the x-axis or on the y-axis.

2. The radius r of C (⃗ p 2 ) is suciently close to a half-integer.

Supposing that ⃗ q 2 is neither on the x-axis nor the y-axis, we see that neither R 2 nor S 2 is zero. We assume that they are positive. In the rst quadrant, the y-coordinate (respectively the x-coordinate) of points in C (⃗ p 2 ) is strictly decreasing with respect to

x (respectively y). Thus it cannot intersect twice with a line parallel to the x-axis (respectively the y-axis).

Supposing that ⃗ q 2 is on the y-axis, we see that the distance between the origin and

C 1 (⃗ q 2 ) (respectively C 3 (⃗ q 2 )
) is greater (respectively smaller) than r. Thus we have r <

√ (R 2 -1 2 ) 2 + ( 1 2 ) 2 (respectively r > √ (R 2 + 1 2 ) 2 + ( 1 2 ) 2 ). Letting ϵ = r -(R 2 -1 2 ) (respectively ϵ = r -(R 2 + 1
2 )) where ϵ < 1 2 , we obtain ϵ(8r -4ϵ) < 1 (in both cases). We experimentally observe that as far as r ≤ 10 5 , ϵ(8r -4ϵ) < 1 never holds. On the other hand, if r is suciently large enough to have 8r -4ϵ ≈ 8r then we have ϵ < 1 8r , which means that r is suciently close to a half-integer. We can thus conclude that [START_REF] Andres | Cercles discrets et rotations discrètes[END_REF] in Remark 3.13 is satised only if r is suciently large. Note that a similar discussion can be done when supposing that ⃗ q 2 is on the x-axis. Corners of H (⃗ q 2 ) inside of C (⃗ p 2 ) and ARA.

Value of I

⃗ q 2 α inf α sup I ⃗ q 2 = 1 or I ⃗ q 2 = 14 α(P 1 , Q 1 , R 2 -1) α(P 1 , Q 1 , 1 -S 2 ) I ⃗ q 2 = 2 or I ⃗ q 2 = 13 α(P 1 , Q 1 , R 2 -1) α(P 1 , Q 1 , -S 2 ) I ⃗ q 2 = 3 or I ⃗ q 2 = 12 α(P 1 , Q 1 , 1 -S 2 ) α(P 1 , Q 1 , -S 2 ) I ⃗ q 2 = 4 or I ⃗ q 2 = 11 α(P 1 , Q 1 , R 2 ) α(P 1 , Q 1 , -S 2 ) I ⃗ q 2 = 6 or I ⃗ q 2 = 9 α(P 1 , Q 1 , R 2 -1) α(P 1 , Q 1 , R 2 ) I ⃗ q 2 = 7 or I ⃗ q 2 = 8 α(P 1 , Q 1 , R 2 ) α(P 1 , Q 1 , 1 -S 2 ) I ⃗ q 2 = 0 and R 2 = 0 α(P 1 , Q 1 , 1 -S 2 ) α(P 1 , Q 1 , 1 -S 2 ) I ⃗ q 2 = 0 and S 2 = 0 α(P 1 , Q 1 , R 2 -1) α(P 1 , Q 1 , R 2 -1)
The main function of our algorithm for nding the lower and the upper bounds of admissible rotation angles, consists of three steps. The rst step sets the rotation center at ⃗ p 1 and ⃗ q 1 , as described in Section 3.4.1. The second step computes the corners that are in C (⃗ q 2 ) and then compute the index

I ⃗ q 2 = ∑ i 2 i × F (C i ( ⃗ q 2 )
). Therefore we can easily identify which corners are in C (⃗ p 2 ) from I ⃗ q 2 . The third step calls a function that returns hinge angles corresponding to I ⃗ q 2 . There are fourteen possible values for I ⃗ q 2 from 0 till 15 except for 5 and 10. Note that geometrically I ⃗ q 2 can be neither 5 nor 10. The value 15 of I ⃗ q 2 implies an error such that all corners are inside of C (⃗ q 2 ). Since I ⃗ q 2 whose value is 0 corresponds to the case ΣF (C i (⃗ q 2 )) = 0, we should verify whether H (⃗ q 2 ) really intersects with C (⃗ p 2 ). Note that for the other values for I ⃗ q 2 , we can make a pair (d, e) such that d + e = 15. The two indices of each pair design the same pair of lower and upper bounds of ARA. Table 1 gives the corresponding lower and upper bound rotation angles for each value of I ⃗ q 2 . Each step of this algorithm has the constant time complexity. Thus the global complexity of this algorithm is also O(1).

Incremental computing lower and upper bounds of rotation angles

In general, the corresponding point sets contain more than two points. Therefore, in this section, we extend our algorithm in Section 3.4.2 to two sets of corresponding point pairs, A and B, each of them having l points where l > 2.

To simplify the notation, we denote by ARA(⃗ p i , ⃗ q i ) = (α i inf , α i sup ) the pair of angles that gives the lower and the upper bounds of admissible rotation angles for the pair of points (⃗ p i , ⃗ q i ). Note that α i inf , α i sup are hinge angles. ARA(A n , B n ) denotes the two most restrictive angles for all points i such as i ≤ n. We recursively dene it by A new algorithm handles all points incrementally. This algorithm is divided into two parts. The rst part is to initialize the algorithm by computing ARA(⃗ p 2 , ⃗ q 2 ). The second part computes ARA(A i , B i ) for i = l. Note that ARA(⃗ p 1 , ⃗ q 1 ) cannot be computed because ⃗ p 1 and ⃗ q 1 are the centers of the rotation.

ARA(A n , B n ) = ARA(A n-1 , B n-1 ) ∩ ARA(⃗ p n , ⃗ q n ).
The time complexity of this algorithm is O(l). As explained in Section 3. Given input data of the algorithm as shown in Figure 3.4 (A), we rst obtain the result of the translation described in Section 3.4.1 as illustrated in (B). We then compare, for each pair of points (⃗ p i , ⃗ q i ) with i ≥ 2, the distance of ⃗ p i from the origin with that of each corner from H (⃗ q i ) to deduce the corresponding hinge angle as explained in Section 3.4.2.

Finally, we obtain (D) that shows the intersection of all ARA(⃗ p i , ⃗ q i ) obtained in (C).

Experimental Evaluation using Synthetic Data

We evaluated our algorithm using synthetic data and veried the quality of obtained admissible rotation angles. To test our algorithm, we need two sets of points. We

Discussion on practical application to digital images

In the previous section, we start with grid points in the discrete plane. In other words, all coordinates of the points are integers. But in reality, points in the Euclidean space are not represented by integers but real numbers. Thus we start here with the set of oating points to see how our algorithm works. Then we test our algorithm with real data acquired by a digital camera.

Synthetic data

To test our algorithm using synthetic data more similar to real data, we applied the Euclidean rotation with angle θ = 50 • directly to the set F that is used in Section 3.5.

Then we discretized the rotated image to obtain the second set I ′′ . To the two sets I and I ′′ , we applied our algorithm. We note that the set I comes from Section 3.5. γ. Let ⃗ p 2 and ⃗ q 2 denote the discretization of ⃗ p 1 and ⃗ q 1 . If we apply our algorithm to the pair of points {⃗ p 2 , ⃗ q 2 }, we obtain ARA(⃗ p 2 , ⃗ q 2 ) = (α inf , α sup ). In this case, γ does not belong to the interval [α inf , α sup ] (Figure 3.8 (right)). We see that this is caused by discretization of the oating points. We can give a bound for this error. In fact, this bound directly depends on the distance between the rotation center and the pair of points. If we denote by d this distance, the maximum error is equal to sin -1 √ 2 2d . Note that the pair of points {⃗ p 3 , ⃗ q 3 } and the angle ρ bring the same problem.

To avoid this problem, we take the example of modifying our algorithm as follows. We keep all the pairs of bound rotation angles determined by all pairs of points. After sorting all lower and upper bound rotation angles into two lists L inf , L sup , we remove all hinge angles α inf from L inf such that for ∃α sup ∈ L sup , α inf > α sup . We also remove all hinge angles α sup from L sup such that for ∃α inf ∈ L inf , α sup < α inf .

With this procedure, we can guarantee that the lower bound hinge angle is smaller than the rotation angle applied to the points and that the upper bound hinge angle is larger than the rotation angle. The complexity O(n log n) is required in sorting the remaining hinge angle pairs, which does not increase the computational cost of the algorithm as a Both methods return a valid ARA. It is preferable to use the second method when there are a few pairs of points because it does not remove any hinge angles. The rst method is preferable when data contains many pairs of points in correspondence because the obtained ARA is more restrictive than the ARA obtained by the second method.

Real data

We applied our algorithm to see its practical usefulness. We used a turntable that is rotated with respect to the vertical axis with respect to a digital image plane. The precision of rotations in control is 10 -3 degrees.

We put a toy block on the turntable and then took its image using a standard digital camera where the camera was xed so that its optical axis is parallel with the rotation axis (Figure 3.9(a)). Next we rotated the turntable with the angle of 44.99 • and then took another image of the toy block by the xed camera (Figure 3.9(b)).

We manually selected ve points in the rst image and their corresponding points in the second image. Then we applied our algorithm to the ve pairs of corresponding points.

Our algorithm might return the empty result (Figure 3.9(c)). This is because no intersection between C (⃗ p i ) and H (⃗ q i ) was found for all i = 1, 2, . . . , 5.

In the case of real data, we cannot always guarantee to detect correct corresponding pairs of points even manually. This indicates that a problem dierent from the discretization problem (see Section 3.6.1) arises 4 . Namely the maximum dierence of distances from the 4 We assume here that we are to nd ARA from given correspondences. rotation center between two points in correspondence can become greater than √ 2. This is illustrated in Figure 3.10 where points ⃗ p and ⃗ q are in correspondence and

ϵ d = |d ⃗ p -d ⃗ q | is greater than √ 2
where d ⃗ p , d ⃗ q are respectively the distances from ⃗ p, ⃗ q to the origin. We see that no intersection exists between the circle C (⃗ p) going through ⃗ p and H (⃗ q).

To avoid this problem on corresponding points, we can change the resolution of the images. Namely, we degrade the image resolution until we nd an intersection between C (⃗ p) and H (⃗ q) 5 . Because of the change in image resolution, we cannot accept the obtained hinge angles as they are. Instead, we add (subtract) a Pythagorean angle θ to (from) the obtained upper (lower) hinge angle. Let assume that d ⃗ p < d ⃗ q and that H (⃗ q) is dened for the image whose resolution is degraded with 2 -n from the original one.

Then we have to choose θ satisfying θ > n We can see that the real rotation angle (44.99 • ) is included in the ARA obtained by our algorithm with this modication. Accordingly, we can conclude that our modication is eective. 5 There are other criteria for nding a reasonable image resolution. For example, Brimkov presents criteria of faithfully digitization for continuous objects in [START_REF] Valentin | Scaling of plane gures that assures faithful digitization[END_REF].

Conclusion

In this chapter, we have presented a new method to use rotation in discrete space.

Because rotations are mathematically dened for the continuous space, it is necessary to develop new operational tools for such discrete rotations.

Our problem was from two digital images to nd all possible rotation angles to rotate the rst image into the second. Namely, to nd the admissible rotation angles between two digital images. Since we need to have only exact computation without any approximation, we decide to adopt hinge angles, which allow us rotations using only integer computation. By using their properties, we have proposed a method for incrementally computing the lower/upper bound rotation angles from a Pythagorean angle for a set of discrete points. Based on this method, we have also proposed a linear algorithm for nding from two digital images the lower and the upper bound rotation angles representing the admissible rotation angles between them.

In addition, we have shown with experiments that our algorithm is ecient with synthetic data and gives consistent results. Then we also have done some experiments to real data.

Results given by our algorithm were not consistent because of the loss of correct point correspondences in our input. We then explained why these problems appear and gave necessary modications to obtain consistent results. These modications did not increase the time complexity while keeping the integer computation.

The range size of ARA can be considered as the unreliability of the point correspondence.

In other words, larger the range of ARA between two points, less is the reliability of their point correspondence. Therefore, it may be interesting if we can propose a new method for evaluating "good" point correspondence from a given pair of digital images by using our results.

In this paper we always assume that rotations centers are on integer points. But with real data it is not always the case. One of our future work will be to adapt the presented method to such rotations. In this purpose, studies on digital discs whose centers are not only on integer points [START_REF] Nagy | An algorithm to nd the number of the digitizations of discs with a xed radius[END_REF][START_REF] Huxley | The number of n -point digital discs[END_REF][START_REF] Pham | Digital circles with non-lattice point centers[END_REF] would help. Chapter 4 3D Rotations in discrete space using hinge angles

Introduction

Rotations in the 3D space are required in many applications for computer imagery, such as image processing [START_REF] Tooli | Three-dimensional rotations by three shears[END_REF], computer vision [START_REF] Schmidt | Using quaternions for parametrizing 3-d rotations in unconstrained nonlinear optimization[END_REF][START_REF] Hartley | Global optimization through rotation space search[END_REF] and computer graphics [START_REF] Chen | 3d volume rotation using shear transformations[END_REF]. A rotation in the 3D Euclidean space can be in general represented in two dierent typical ways. One way is to represent a rotation as a combination of three particular rotations around the three axes of the coordinate system [START_REF] Voss | Discrete Images, Objects and Functions in Z n[END_REF]. The other way is to represent a rotation by a rotation axis together with an angle around the axis [START_REF] Schmidt | Using quaternions for parametrizing 3-d rotations in unconstrained nonlinear optimization[END_REF][START_REF] Hartley | Global optimization through rotation space search[END_REF]. Even if the representations of a rotation are dierent, computed rotation results are the same as far as the space is continuous. However, this is not the case for the discrete space. Namely, depending on the rotation representation, the computed rotation result can change in the discrete space [START_REF] Andres | Cercles discrets et rotations discrètes[END_REF]. As this is the case of 2D rotations, computing a 3D rotation once in the discrete space brings displacement from in the continuous space; computing 3D rotations many times causes diculty in analyzing inherited displacements during the computation. Accordingly, representing a 3D rotation by a rotation axis together with an angle around the axis is more preferable in the 3D discrete space. Besides, it is known that such axis-angle representation is useful for 3D rotation estimation from given image sets, which is one of the important problems encountered in computer vision [START_REF] Schmidt | Using quaternions for parametrizing 3-d rotations in unconstrained nonlinear optimization[END_REF][START_REF] Hartley | Global optimization through rotation space search[END_REF].

This chapter presents a study of the rotation in the 3D discrete space. Since we admit only integer computations, we assume that our rotation center is a grid point such as the origin, and that a rotation axis has integer coordinates. In the 2D case, hinge angles are known to be corresponding to the discontinuity caused by discretization of the rotation in the continuous plane as explained in Chapter 3 and in [START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF][START_REF] Amir | Two-dimensional pattern matching with rotations[END_REF]. Intuitively, hinge angles determine a transfer of a grid point from a pixel to its adjacent pixel during the rotation. In other words, two rotations with nearby angles transform the same grid point to two adjacent pixels because discrete rotations around a given center are locally continuous with regard to the angle. Hinge angles are their discontinuity angles.

Computing hinge angles using integers alone allows us to compute 2D discrete rotations without any approximation errors, which designs the 2D discrete rotation. Extending these to the 3D case, we design a 3D discrete rotation. In the 3D case, however, depending on the rotational axis, we have a variety of transitions of a grid point across voxels. How to capture these transitions systematically is a big issue.

In this chapter, we rst dene hinge angles for 3D rotations so that they determine a transit of a grid point from a voxel to its adjacent voxel. To compute the hinge angles for 3D rotations, we introduce a notion, called "multi-grids", that is given by the intersection between a plane normal to the rotation axis and the half-grids that mark the boundary between adjacent voxels. The rational multi-grids, which are a subset of multi-grids, allow computations of hinge angles using only rational numbers. Using rational multigrids, we show that, as in 2D, it is possible to only use integer during computation and to have an integer representation of hinge angles. Then, we give a method to sort all the possible hinge angles in concern to design a 3D discrete rotation. We also propose a method to obtain from a pair of 3D digital images in correspondence a set of 3D rotations, each of them transforming the rst digital image into the second one. Note that in a discrete space, there are sets of rotations that give the same rotated image of a given digital image. In this paper, we x a rotation axis and look for all the possible rotation angles from a given pair of 3D digital images. The set of all possible rotation angles is called admissible rotation angles and its upper and lower bounds are represented by hinge angles. This method is the extension of the method proposed in Chapter 3 for the 2D cases into 3D cases.

Dierently from 2D discrete rotations [START_REF] Andres | Cercles discrets et rotations discrètes[END_REF][START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF][START_REF] Andres | The quasi-shear rotation[END_REF], few attempts on 3D discrete rotations have been reported [START_REF] Chen | 3d volume rotation using shear transformations[END_REF][START_REF] Tooli | Three-dimensional rotations by three shears[END_REF]. In particular, to our best knowledge, this is the rst work on 3D discrete rotations using integer computations without digitization errors.

Hinge angles

Hinge angles for 2D rotations are dened to represent the discontinuities of rotations in the discrete plane [START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF]. Hinge angles determine a transit of a grid point from a pixel to its adjacent pixel during the rotation. To characterize those hinge angles, the 2D half-grid plays an important role. In other words, the 2D half-grid represents the border between two adjacent pixels. From the denition of the 2D half-grid, we dene the hinge angles in the plane.

Denition 4.2. An angle α is called a hinge angle if at least one point in Z 2 exists such that its image by the Euclidean rotation with α around the origin is on the 2D half-grid. Note that hinge angles in other quadrants are obtained by symmetry with respect to the

x-axis and/or y-axis from those in Figure 4.1.

To extend the denition of hinge angles into the 3D case, we rst dene the half-grid in the 3D space. Similarly to the 2D half-grid, the 3D half-grid denes the limit between two adjacent voxels in the 3D discrete space.

Denition 4.3. The 3D half-grid is the set of planes in the 3D space, each of which is represented by one of x = i + 1 2 , y = i + 1 2 and z = i + 1 2 where i ∈ Z.

Introducing the denition of the 3D half-grid allows the denition of hinge angles in 3D as a natural extension of hinge angles in 2D. As mentioned in the introduction, we only consider here 3D rotations whose rotation axes have directional vectors with integer coordinates and go through the origin. Hereafter, we call such an axis an integer-axis.

Denition 4.4. An angle α is called a hinge angle if at least one point in Z 3 exists such that its image by the Euclidean rotation with α around an integer-axis is on the half-grid.

Similarly to the 2D case, for a grid point ⃗ p in 3D, an angle α is a hinge angle if and only if the discretized point of the rotation result of ⃗ p with angle α + ϵ becomes dierent from that with angle α -ϵ for any ϵ > 0.

Dierently from the case of 2D rotations, we need not only a rotation angle but also a rotation axis in order to specify a 3D rotation. This requires investigation of the intersection between voxels and a plane determined by a given rotation axis because of a variety of transitions of a grid point across voxels existing on the plane. To capture this variety, we introduce the multi-grid in the next section, which allows us to study hinge angles in the 3D rotation plane.

Multi-grids

In this section, we introduce a notion that is required to extend hinge angles from 2D to 3D and thus required to perform discrete rotations in the 3D discrete space. As explained in Section 4.2, hinge angles are strongly related to the half-grid both in 2D and 3D. However, in the 3D cases, rotations of a point are always in its rotation plane normal to the rotation axis and goes through the point. In fact, the 3D half-grid is not well adapted to describe rotations for a grid point. We thus consider the intersection between a rotation plane and the 3D half-grid, which is a planar grid consisting of three sets of parallel lines, called a multi-grid instead of the half-grid.

In this section, we give a formal denition of multi-grids. Then we show how to obtain the line equations of a multi-grid from a grid point and a rotation axis. Then, we restrict multi-grids to the rational multi-grids that form a set of multi-grids useful to perform discrete rotations in the 3D discrete space. Finally we present some useful properties of rational multi-grids.

Denition of multi-grids

When a rotation plane in 3D is given, the intersection between the plane and the halfgrid in the 3D space is obtained as illustrated in In a multi-grid, the interval between parallel lines having the same directional vector is regular. Normalizing the interval allows us to represent each set of parallel lines having the same normal vector (a i , b i ) ⊤ as

L i = {(x, y) ⊤ ∈ R 2 |a i x + b i y + c i + k = 0, a 2 i + b 2 i ̸ = 0, k ∈ Z, a i , b i , c i ∈ R} (4.1)
where i = 1, 2, 3. The integer parameter k denotes the index number of each parallel line. Figure 4.3 gives a geometrical explanation of the parameters of L i . For example, if a point (x, y) ⊤ is on one of the parallel lines of L i , (x -k a i , y) ⊤ is on the k-th next line, providing that a i ̸ = 0. Now we can give a formal denition of a multi-grid. The ve dierent shapes of convexels, which are constructed as the intersections between a rotation plane and voxels [START_REF] Andres | New methods in oblique slice generation[END_REF].

The multi-grid in the rotation plane forms various closed convex polygons surrounded by lines, that we call convexels. Depending on the rotation plane, we have a variety of shapes of convexels in this paper. The convexel in a multi-grid is the counterpart of the squared pixel dened in the 2D half-grid. The shapes of convexels are investigated in [START_REF] Andres | New methods in oblique slice generation[END_REF] under the context of the intersection of a voxel and a plane, and the number of vertices of a convexel can be 3, 4, 5 or 6 as illustrated in Figure 4.4. We remark that when the normal of the rotation plane is parallel with one of the axes dening the coordinate system of the 3D space, the notion of convexels coincides with the notion of pixels.

Multi-grid line equations

To simplify the derivation of line equations for each L A,⃗ p i , we introduce a new base B i where lines in L A,⃗ p i are parallel with the x-axis. In the following, we derive the equations for L A,⃗ p 1 . Note that the same discussion can be applied to L A,⃗ p

L A,⃗ p 1 as ⃗ v 1 = A ∧ ⃗ e 1
where ⃗ e 1 = (1, 0, 0) ⊤ . We set ⃗ v 2 = ⃗ v 1 ∧A ∥A∥ , which is orthogonal to ⃗ v 1 . Note that both ⃗ v 1 and ⃗ v 2 are orthogonal with respect to A.

We introduce a new base B 1 in such a way that ⃗ v 1 and ⃗ v 2 respectively become ⃗ u 1 = (1, 0) ⊤ and ⃗ u 2 = (0, 1) ⊤ in P. The transformation from B to B 1 is realized by

P BB 1 =   0 az a 2 y +a 2 z -ay a 2 y +a 2 z -ψ ψaxay a 2 y +a 2 z ψaxaz a 2 y +a 2 z   , ( 4.2)
where ψ

= 1 √ a 2 x +a 2 y +a 2 z
. We remark that A is transformed to (0, 0) ⊤ by P BB 1 ; the rotation center in P thus becomes the origin of B 1 .

We remark that if A is collinear with one of the axes dening the coordinate system of the 3D space, P BB i degenerates: the rank of M BBi becomes 1. In such cases, 3D rotations become identical with 2D rotations. However, in Section 4.3.1, we explained that we do not consider the cases where the vector A is collinear with an axis dening the coordinate system in concern, so that lines in L A,⃗ p i are not orthogonal to those in

L A,⃗ p j
where i ̸ = j. Thus, we do not take these particular cases in consideration.

Applying

P BB 1 to the plane x = k for k ∈ Z induces a line in M A,⃗ p 1 whose equation is ψ(a 2 y + a 2 z )y + k -A.⃗ pψ 2 a x = 0. (4.3)
Changing the roles between L A,⃗ p 1 and L A,⃗ p 2 (resp, L A,⃗ p 3 .) and between ⃗ e 1 and ⃗ e 2 = (0, 1, 0) ⊤ (resp, ⃗ e 3 = (0, 0, 1) ⊤ .), we obtain the transformation matrices P BB i for i = 2, 3 such that

P BB 2 =   -az a 2 y +a 2 z 0 ax a 2 y +a 2 z ψaxay a 2 y +a 2 z -ψ ψayaz a 2 y +a 2 z   , ( 4.4
)

P BB 3 =   ay a 2 y +a 2 z -ax a 2 y +a 2 z 0 ψaxaz a 2 y +a 2 z ψayaz a 2 y +a 2 z -ψ   , ( 4.5)
and the line equations for L A,⃗ p 2 and L A,⃗ p 3 such that:

ψ(a 2 x + a 2 z )y + k -A.⃗ pψ 2 a y = 0, (4.6
)

ψ(a 2 x + a 2 y )y + k -A.⃗ pψ 2 a z = 0, (4.7)
where k ∈ Z.

We note that (4.3), (4.6) and (4.7) correspond to the equation of (4.1) for i = 1, 2, 3

respectively. All a i of (4.1) for i = 1, 2, 3 are null since each B i is set such that the lines of L A,⃗ p i are parallel to the x-axis. We also remark that every b i of (4.1) depends only on A, but not on ⃗ p. Indeed, c i is the only parameter depending on A and ⃗ p. This implies that all the multi-grids generated from the same normal vector A with dierent point

Rational multi-grids

In Section 4.3.2, we obtained (4.3), (4.6) and (4.7) for the lines in the multi-grid M A,⃗ p .

In general, parameters of these equations belong to R. In order to use only integers during computation, we need these parameters to belong to Q.

If all elements of P BB i belong to Q, then all the parameters of (4.3),(

become rational as well. In order to obtain rational parameters in P BB i we set A to be a Pythagorean vector; a vector ⃗ v = (i 1 , i 2 . . . , i n ) ⊤ , i 1 , i 2 , . . . , i n ∈ Z, is a Pythagorean vector if ∥⃗ v∥ = λ where λ ∈ Z. Note that a rotation axis whose directional vector is such a Pythagorean vector is called a Pythagorean axis. This assumption ensures that ψ becomes a rational value. We call a multi-grid dened from a Pythagorean axis a rational multi-grid. Note also, that every axis can be approximated by a Pythagorean axis. A study on Pythagorean axis is presented in Chapter 6.

The rational multi-grids presented in this section are a subset of the multi-grid. In addition of the possibility to use only integers during the computations, rational multigrids also oer some useful properties that are valid only for rational multi-grids.

Properties of rational multi-grids

In the case of multi-grids that are not rational, there is an innity of dierent convexels.

However, for rational multi-grids, the number of convexels is nite and depend only on the coordinates of A that is the normal vector of rotation plane. Here, we consider any grid point ⃗ p for M A,⃗ p . A prime Pythagorean axis is a Pythagorean axis A = (a x , a y , a z ) ⊤ such that gcd(a x , a y , a z ) = 1. We dene the arithmetical rest of a voxel ⃗ v as follow:

R(⃗ v) = a x v x + a y v y + a z v z + b. (4.8)
Theorem 4.7. Let M A,⃗ p be a rational multi-grid associated to the prime Pythagorean vector A = (a x , a y , a z ) ⊤ . The number of dierent convexels in M A,⃗ p for any ⃗ p is either

|a x | + |a y | + |a z | or |a x | + |a y | + |a z | + 1.
Proof. 7. Let P be a plane of equation a x x + a y y + a z z + b = 0. We know [START_REF] Andres | New methods in oblique slice generation[END_REF] that the voxel around the grid point ⃗ v = (v x , v y , v z ) ⊤ is intersected by P if and only if:

- |a x | + |a y | + |a z | 2 ≤ R(⃗ v) ≤ |a x | + |a y | + |a z | 2 . ( 4
- |a x | + |a y | + |a z | 2 < R(⃗ v) < |a x | + |a y | + |a z | 2 (4.10)
and the number of dierent values admissible for R(⃗ v) is a x + a y + a z .

Note that, if A is a Pythagorean vector but not prime, the number of dierent convexels will be either Such a repetition of convexels in M A,⃗ p is shown in the next theorem, which is easily derived from [START_REF] Kenmochi | On combinatorial properties of discrete planar surfaces[END_REF] given in the context of discrete planar surfaces.

Theorem 4.8. Let M A,⃗ p be a rational multi-grid associated to the prime Pythagorean vector A = (a x , a y , a z ) ⊤ and a grid point ⃗ p. There are three integers such that A x = L ax , A y = L ay , A z = L az where L = lcm(a x , a y , a z ) such that for any grid point ⃗ q = (q x , q y , q z ) ⊤ , we have a quadruple of grid points {⃗ q 1 = ⃗ q, ⃗ q 2 = (q x +A x , q y -A y , q z ) ⊤ , ⃗ q 3 = (q x -A x , q y , q z + A z ) ⊤ , ⃗ q 4 = (q x , q y + A y , q

z -A z ) ⊤ } satisfying R(⃗ q 1 ) = R(⃗ q 2 ) = R(⃗ q 3 ) = R(⃗ q 4 ).
This indicates that the triple of integers A x , A y , A z describes the frequency of the repetition of convexels in a rational multi-grid M A,⃗ p . Let P be the support plane of M A,⃗ p . For each grid point ⃗ v = (v x , v y , v z ) intersected by P, the grid point ⃗ v ′ obtained by translation of v by a linear combination of two vectors (A x , A y , 0) ⊤ and (A x , 0, A z ) ⊤ has the same formed convexel as that of ⃗ v in P.

Another consequence of Theorem 4.8, illustrated in Figure 4.6, is the existence of triangles that tile the multi-grid. Indeed, let ⃗ q 1 = (q x , q y , q z ) be a grid point and ⃗ q 2 = (q x +A x , q y + A y , q z ), ⃗ q 3 = (q x + A x , q y , q z -A z ) and ⃗ q 4 = (q x + 2A x , q y + A y , q z -A z ) be three points resulting of a translation of ⃗ q. The two triangles of vertex {⃗ q 1 , ⃗ q 2 , ⃗ q 3 } and {⃗ q 2 , ⃗ q 3 , ⃗ q 4 } Figure 4.6: A quadruple of grid points {⃗ q 1 = ⃗ q, ⃗ q 2 = (q x + A x , q y + A y , q z ) ⊤ , ⃗ q 3 = (q x + A x , q y , q z -A z ) ⊤ , ⃗ q 4 = (q x + 2A x , q y + A y , q z -A z ) ⊤ } forms two types of triangles, α and β tiled in a multigrid.

contain the same set of convexels and are mirror images. Note that triangles formed by four triangles, three α and one β or one α and three β, also tile the plane and contain all convexels. If A x , A y , A z are chosen to be prime, they describe the smallest triangle that tiles M A,⃗ p .

Hinge angles characterized by a multi-grid

The goal of this section is to show the relation between multi-grids and hinge angles, and then to show how to obtain the unique representation of an integer quintuplet for a 3D hinge angle, namely an injective map from hinge angles to quintuplets. In 2D, any hinge angle can be uniquely represented by a triple of integers as shown in Section 3.2.1. In order to ensure that the representation is unique, some properties on the multi-grid are required. For 3D hinge angles, similar properties on multi-grid are required. Therefore, we rst show some properties on the 3D half-grid, then we represent 3D hinge angles using ve integers and then explain how to decode them to obtain their hinge angles.

In Section 4.2, we dene the 3D hinge angles in the framework of the 3D half-grid.

As manipulation of 3D hinge angles is more convenient in the framework of multi-grids and rational multi-grids, we propose an alternative denition of 3D hinge angles in the following Proposition 4.9. Let ⃗ p be a grid point and ⃗ p ′ be the result of the rotation of ⃗ p by an angle α around an integer-axis. If α is a hinge angle, then ⃗ p ′ is on the multi-grid.

Hinge angles and rational multi-grids

In 2D, there is a property on hinge angles ensuring that the locus of rotation of a grid point cannot contain the intersection of two lines belonging to the half-grid [START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF]. In this section, we show that a similar property for 3D hinge angles in the rotation plane holds if the rotation axis is an integer axis. In multi-grids, the locus of rotation of a grid point in a rotation plane cannot contain the intersection of two or three lines of the multi-grid as illustrated in Proof. 8. The equation of the rotation plane P of M A,⃗ p is a x x + a y y + a z z -A.⃗ p = 0.

Let ⃗ p ′ = (p ′ x , p ′ y , p ′ z )
⊤ be a point which belongs to the locus of the rotation of ⃗ p in P.

Let us assume that ⃗ p ′ is also a vertex of a convexel of M A,⃗ p , so that it belongs to two planes of the 3D half-grid. Thus we can set, without loss of generality, that p ′

x = k x + 1 2 and p ′ y = k y + 1 2
where k x , k y ∈ Z, and then

⃗ p ′ = (k x + 1 2 , k y + 1 2 , p ′ z ) ⊤ .
The locus of the rotation of ⃗ p is the intersection between P and the sphere S:

x 2 + y 2 + z 2 -(p 2 x + p 2 y + p 2 z ) = 0.
Thus, by the assumption that ⃗ p ′ belongs to P and S, we have

a x ( k x + 1 2 ) + a y ( k y + 1 2 ) + a z p ′ z -A.⃗ p = 0, (4.11) ( k x + 1 2 ) 2 + ( k y + 1 2 ) 2 + p ′2 z -(p 2 x + p 2 y + p 2 z ) = 0. (4.12)
From (4.11) we see that p ′ z must be a rational number, so that there is a pair of integers λ 1 , λ 2 such that p ′ z = λ 1 λ 2 and gcd(λ 1 , λ 2 ) = 1. From (4.12) we see that p ′2 z + 1 2 must be an integer. Thus we have

p ′2 z = λ 2 1 λ 2 2 = k + 1 2 , (4.13)
where k, λ 1 , λ 2 ∈ Z. Since gcd(λ 1 , λ 2 ) = 1 and 2λ 2 1 = (2k + 1)λ 2 2 , λ 2 must be even. Setting λ 2 = 2λ ′ 2 where λ ′ 2 ∈ Z, we then obtain 2λ 2 1 = 4(2k + 1)λ ′2 2 and deduce that λ 1 must be also even, which contradicts the assumption that gcd(λ 1 , λ 2 ) = 1. Therefore we can conclude that there is not such a point ⃗ p ′ .

This lemma shows that two adjacent convexels, whose transition of ⃗ p during its rotation is done between them, always share a convexel edge. In other words, the rotation locus of ⃗ p passes through a sequence of voxels such that any successive voxels are connected by their common face. Thus, Lemma 4.10 shows that cases (a) and (b) illustrated in Figure 4.7 cannot happen in a rational multi-grid. Note that the above proof only needs the assumption that A has integer coordinates.

Note that in 2D, there is also another property on hinge angles ensuring that the locus of rotation of a grid point cannot intersect twice on the same line of the half-grid without intersecting another line of the half-grid between the two intersections. In 3D, this property does not hold with rational multi-grids. More details will be discussed in Section 5.2.2.

Quintuplet integer representation of hinge angles

A hinge angle in a given rotation plane is represented by a quintuplet of integers (p x , p y , p z , i, k)

given by α(p x , p y , p z , i, k). The rst three integers p x , p y , p z represent, in the 3D basis B, the coordinates of ⃗ p. The fourth integer i indicates the index number for the set

L A,⃗ p i , i = 1, 2, 3
, where the hinge angle α is dened. The last integer k represents the index number of the line in L A,⃗ p i . Therefore, a quintuplet of integers keeps the coordinates of ⃗ p and the information required for obtaining the coordinates of the arriving point ⃗ q after the rotation of ⃗ p by α.

From these ve integers, we can obtain, in the basis B i whose projection matrix P BB i is given by (4.2), (4.4) and (4.5), the coordinates (P x , P y ) ⊤ and (Q x , Q y ) ⊤ corresponding to ⃗ p and ⃗ q, as represented in Since ⃗ q belongs to the locus of the rotation of ⃗ p, we have Q 2

x + Q 2 y = P 2 x + P 2 y . We then notice that Q x is not a rational number. However, since P x , P y and Q y are rational, all computations involving Q x can be done using only integer computations. In general, Q x can take two values which dene two points. To discriminate two dierent hinge angles corresponding to these two points by the integer quintuplets, we add the positive sign to the fourth integer for the greater Q x and the negative sign for another value of Q x . If Q x = 0, then the two points merge into one. This particular case where the locus of the rotation of ⃗ p is tangent with the k-th line of M A,⃗ p i dene no hinge angle since there is no transition of convexel. Hereafter, the representation of hinge angles by integer quintuplets will be denoted by α(p x , p y , p z , ±i, k)

Note that according to Lemma 4.10 we know that any hinge angle cannot rotate a grid point to the intersection of more than one lines of a multi-grid. Therefore, we have the following theorem:

Theorem 4.11. Let α p (p x , p y , p z , i p , k p ) and α q (q x , q y , q z , i q , k q ) be two hinge angles with their integer quintuplet representations. α p = α q if and only if p x = q x , p y = q y , p z = q z , i p = i q and k p = k q .

Theorem 4.11 is rephrased as: two dierent integer quintuplets cannot represent the same hinge angle.

We remark that in [START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF], hinge angles for 2D discrete rotations are represented with a triple of integers that represents the 2D coordinates of the point and the index number of the line that is intersected by the locus of the point's rotation. There is no integer that decides whether the intersected line is parallel to the x-axis or the y-axis. To dierentiate both cases, the last integer of the triple is positive if the line is parallel to

x-axis or otherwise, it is negative.

3D discrete rotations around a Pythagorean axis

In this section, we develop a 3D discrete rotation based on hinge angles. This rotation is the extension to the 3D space of the 2D discrete rotation given in Section 3.3. In order to obtain the complexity of the algorithm described in this section, we need to enumerate all the hinge angles existing for an image to rotate. In the case of 2D hinge angles, the upper bound of the number of dierent hinge angles for an image is O(m 3 ) where m is the largest side of the image [START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF][START_REF] Amir | Two-dimensional pattern matching with rotations[END_REF]; we will use a similar method for obtaining the upper bound of hinge angles for a 3D image. Our algorithm, to be ecient, needs to compare a pair of hinge angles in constant time. Besides, in order to keep integer computations, we need that the comparison is done using only integers. Note that the hinge angle comparison in constant time is not trivial indeed due to our integer computation constraint. After

showing how to compare a pair of hinge angles in constant time with integer computation, we will give the upper bound of the number of hinge angles for a 3D digital image and then present a 3D discrete rotation algorithm for such a given image.

In Section 4.5 we assume that A is given to be a prime Pythagorean vector.

Comparing hinge angles with integer computations

From the integer representation of hinge angles we can obtain the sine and cosine of a hinge angle α(p x , p y , p z , i, k) characterized in the base B i by the points ⃗ p that have the coordinates (p x , p y , p z ) ⊤ in B and (P x , P y ) ⊤ in B i and ⃗ q that have the coordinates We remark that if the multi-grid is a rational multi-grid, P x , P y , Q y and Q 2 x are rational. Thus it is possible to compare the sine and the cosine of two hinge angles using only integer computations as follows.

(Q x , Q y ) ⊤ in B i . P x ,
Proposition 4.12. Let α 1 and α 2 be two hinge angles dened for A. Then it is possible to decide if α 1 > α 2 using only integer computations.

Proof. 9. Let α 1 = α(p 1x , p 1y , p 1z , i 1 , k 1 ) and α 2 = α(p 2x , p 2y , p 2z , i 2 , k 2 ). Comparing α 1 and α 2 is equivalent to comparing their sines and cosines which are given in (4.14) and (4.15). First we compare the signs of both sine and cosine between α 1 and α 2 . If the sines and cosines of both angles have dierent signs, then we can conclude whether α 1 > α 2 without other computation. Otherwise, without loss of generality, we can assume that both α 1 and α 2 belong to [0, π 2 ], so that cos α i ≥ 0 and sin α i ≥ 0 for both i = 1, 2. As the method for comparing two sines is similar to the one for comparing two cosines, we will only show the later one.

A 1 + B 1 > A 2 + B 2 . (4.17)
Squaring both sides of (4.17) since they are not negative, and moving rational values to the left-hand side and the rest to the right-hand side, we obtain

A 2 1 + B 2 1 -A 2 2 -B 2 2 > 2A 2 B 2 -2A 1 B 1 . (4.18)
If the left-hand side and the right-hand side of (4.18) do not have the same sign, then we can conclude whether α 1 > α 2 or α 2 > α 1 . We can check the sign of both sides of (4.18)

with integer computations since the left-hand side contains only rational numbers and the sign of the right side is the same as

A 2 2 B 2 2 -A 2 1 B 2
1 which also contains only rational numbers. If signs of both sides are the same, assuming that they are positives, we square both sides of (4.18) to obtain

(A 2 1 + B 2 1 -A 2 2 -B 2 2 ) 2 -4A 2 1 B 2 1 -4A 2 2 B 2 2 > -8A 1 B 1 A 2 B 2 . (4.19)
If the sign of the left-hand side of (4. [START_REF] Paeth | A fast algorithm for general raster rotation[END_REF]) is positive, we can deduce that α 2 > α 1 .

Otherwise, taking the square of each side gives us

[ (A 2 1 + B 2 1 -A 2 2 -B 2 2 ) 2 -4A 2 1 B 2 1 -4A 2 2 B 2 2 ] 2 < 64A 2 1 B 2 1 A 2 2 B 2 2 . (4.20)
We note that we can easily verify whether (4.20) is satised with integer computation alone. If (4.20) is true, α 1 < α 2 ; otherwise α 2 < α 1 . The same logic can be applied to the case where the signs of the both sides of (4.18) are negative.

Thanks to Proposition 4.10, a hinge angle cannot have two quintuplet integer representations. Thus we can conclude that the comparison of a pair of hinge angles α 1 , α 2 is always possible with integer computation if they have dierent quintuplets.

Note that, if the comparison of two hinge angles is done using oating numbers, only one comparison is required. However, if we want to keep integer computations, then in the worst case we have to check (4.18), (4. [START_REF] Paeth | A fast algorithm for general raster rotation[END_REF]) and (4.20). From Proposition 4.12, we now have a guarantee of a constant number for each comparison. Then we conclude that the comparison of a pair of hinge angles is done in constant time.

The following proposition is required for our algorithm of 3D discrete rotation. In the 2D case, is it known that the comparison between a hinge angle and a Pythagorean angle can also be done using only integers during its computation and in constant time. An angle is a Pythagorean angle if its sine and cosine belong to Q [START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF]. In the 3D case, a similar proposition will be still valid.

Note that the proof of Proposition 4.13 is similar to the proofs of Proposition 4.12 and on comparison for a hinge angle and a Pythagorean angle in the plane presented in Section 3.3.1.

Proposition 4.13. Let θ be a Pythagorean angle and α be a hinge angle dened for A.

Then it is possible to decide if α > θ in constant time with integer computations.

Proof. 10. Let α = α(p x , p y , p z , i, k) be an hinge angle and θ be a Pythagorean angle associated to the Pythagorean triple (i 1 , i 2 , λ) such that cos θ = i 1 λ and sin θ = i 2 λ where

i 1 , i 2 , λ ∈ Z.
If α is greater than θ, cos θ -cos α > 0 thus we have from (4.14),

i 1 (P 2 x + P 2 y ) > λ(P x Q x + P y Q y ). (4.21)
Using similar methods to those used for the proof of Proposition 4.12, we obtain the equation:

(i 1 (P 2 x + P 2 y ) -λP y Q y ) 2 > λ 2 P 2 x Q 2 x . (4.22)
We note that we can easily verify whether (4.22) is satised with integer computation alone. If (4.22) is true, α < θ; otherwise θ < α.

Upper bound of the number of 3D hinge angles

In the 2D case, the upper bound of the number of hinge angles for a given image of size

m × m is known to be O(m 3 ) [4]
. This is obtained by computing the bound for the furthest point from the origin in the image and multiplying this bound by the number of points in the image. To compute the number of hinge angles in the 3D case, we will use a similar method.

In the 3D case, we assume that an image of size m × m × m is given. The number of hinge angles for a given point ⃗ p depends on the distance between ⃗ p and the axis of

3D discrete rotations induced by hinge angles

In this section, we explain how to design a discrete rotation of a 3D digital image using hinge angles for a given Pythagorean axis of rotation. This method is the 3D extension of the 2D discrete rotation described in Section 3.3.3. As input of such discrete rotation, we have a digital image I of size m × m × m, a vector A, supposed to be a prime Pythagorean vector, and an angle θ supposed to be Pythagorean. The assumption that the rotation axis is a Pythagorean axis and the angle is a Pythagorean angle does not restrict the eld of possible rotations. Indeed, it is proved in [START_REF] Schmutz | Rational points on the unit sphere[END_REF] that the Pythagorean vectors are dense on the 3D unit sphere and it is proved in [START_REF] Anglin | Using pythagorean triangles to approximate angles[END_REF] that any angle can be approximated with a small dierence ϵ > 0 by a Pythagorean angle. The output of our algorithm is a rotated digital image I ′ .

The rotation algorithm is described in Algorithm 4. For each point ⃗ p = (p x , p y , p z ) ⊤ in the image I, the algorithm computes the corresponding multi-grid M A,⃗ p and searches for each line k-th in M A,⃗ p i , i = 1, 2, 3, a pair of hinge angles α(p x , p y , p z , ±i, k). Then we stock and sort all hinge angles corresponding to ⃗ p using Proposition 4.12. The algorithm searches in the sorted list L the greatest hinge angle α which is smaller than θ using Proposition 4.13. This operation can be done using only integer computations thanks to our assumption that our input angle θ is a Pythagorean angle. Finally the new point after the rotation of ⃗ p by α is generated in I ′ .

The time complexity of this algorithm is O(m 4 log m). The computation and the sorting of all hinge angles for each point is done in O(m log m) operations because the comparison between two hinge angles is done in constant time according to Proposition 4.12.

Searching the largest hinge angle α smaller than θ is done in O(log m) operations because the comparison between a hinge angle and a Pythagorean angle can be performed in constant time according to Proposition 4.13. Therefore, the time complexity of this algorithm is O(m 4 log m) because we repeat m 3 times the previous operations.

3D incremental discrete rotation

In this section, we present an incremental algorithm that performs a 3D discrete rotation.

In [START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF], the author describes a 2D incremental rotation that allows obtaining all possible congurations of 2D discrete rotations for a given image. Similarly, 3D hinge angles on rotation planes also allow us to design a 3D incremental discrete rotation. This algorithm may help us to understand the congurations of 3D discrete rotations even thought we do not see yet its practical uses.

For the incremental rotation algorithm, we consider the input data to be a digital image I of size m × m × m and a given rotation axis.

The incremental algorithm consists in four main steps. The rst step generates a table of m × m × m pointers that refers to every voxel of I. Then, for each point in I it computes and sorts all hinge angles associated to this point, and stores them into a list. The third step merges the lists of all points into the nal list and then sorts all the hinge angles in the list. The nal step browses in ascending order the list containing all the hinge angles and displaces the voxel corresponding to the current hinge angle to its neighborhood.

Firstly, the algorithm initializes a table T of m × m × m pointers, such that each pointer (i, j, k), where i, j, k ∈ [0, m], refers to the grid point (i, j, k). T is used to track the voxel location during the incremental rotation. Then for each point ⃗ p = (p x , p y , p z ) ⊤ in I, the algorithm computes and stores in the list T ⃗ p all hinge angles α(p x , p y , p z , i, k) for i = 1, 2, 3 and k = -l, -l + 1, . . . , l -1, l where l = ⌊ √ p 2 x + p 2 y + p 2 z ⌋. Hinge angles are stored using their integer representations. Then each list T ⃗ p is made in ascending order of hinge angles. The second step merges all the lists T ⃗ p into a unique list. In order to avoid the sorting of the nal list T this operation should be done using the merge sort algorithm [START_REF] Kronrod | Optimal ordering algorithm without operational eld[END_REF]. The last step of the algorithm performs the incremental rotation. Reading every hinge angle α(p x , p y , p z , i, k) in T , we displace the grid point referred by T (p x , p y , p z ) to one of its 6-connected voxel according to i and k. The pointer in T (p x , p y , p z ) is updated to the new position so that the actual position of ⃗ p will be known using T (p x , p y , p z ).

We repeat this operation while hinge angles are remaining in T . When all the hinge angles in T have been processed, every point has described its full rotation locus around a given axis of rotation. Then all the possible congurations of rotations of I around the rotation axis have been reached.

The time complexity for the rst step is O(m 4 log m) since there is m 3 points to consider and for each point, there are at most O(m) hinge angles. Thus the sorting of each list is done in O(m log m). The second step requires the same time complexity than the rst step. Note that if the merging does not use the method presented in [START_REF] Kronrod | Optimal ordering algorithm without operational eld[END_REF], the time complexity remains the same, but it will practically increase the computation time. The last step browses the list T that contains at most O(m 4 ) hinge angles. Each displacement is done in constant time and thus it requires O(m 4 ) operations. We conclude that the time complexity for the incremental rotation algorithm is O(m 4 log m).

Conclusion

In this chapter, we extended the notion of hinge angles, introduced Chapter 3 to the 3D. Extension of hinge angles from the 2D to 3D space involves many problems because most of the 2D hinge angles properties are not valid for 3D hinge angles. In order to regard hinge angles in the 3D space similarly to the 2D ones, we introduced the notion of a multi-grid that is the intersection of the 3D half-grid and a rotation plane. By redening the hinge angles on the rotation plane, which are the extension of hinge angles for the rotation in 2D, we showed a subgroup of the multi-grids where all parameters are rational, called rational multi-grids. This subgroup allows us to compare two hinge angles on rotation planes in constant time by using integer computations. It also allows us to design a 3D discrete rotation and a 3D incremental rotation.

The multi-grids introduced in this paper can be extended in any dimension to perform discrete rotations. Roughly speaking, a rotation in nD requires ⌊ n 2 ⌋ angles and rotation planes. However, rotations in nD are not well understood, for example, in Euclidean space we do not know yet if a given nD rotation can be uniquely decomposed into ⌊ n 2 ⌋ planar rotations [START_REF] Richard | An algorithm to decompose n-dimensional rotations into planar rotations[END_REF]. Therefore, the extension of the work presented in this paper to nD discrete rotations requires a study of decomposition of nD discrete rotation into 2D discrete rotation.

Chapter 5

Estimation of rotation between a pair of 3D digital images

Introduction

This chapter is the extension to the 3D space of our algorithm to nd the admissible rotation angles for a pair of 2D digital images presented in Chapter 3. We present a discrete method to nd from a pair of 3D digital images admissible rotations that transform the rst digital image into the second digital image. We assume that the points correspondence between these two 3D digital images are given.

The estimation of a rotation in 3D discrete space is done by estimating a rotation angle and a rotation axis. It is also possible to estimate a 3D discrete rotation by estimating three angles and the rotation center. However, as explained in Chapter 2 and Chapter 4, it is preferable to avoid rotations composition.

Hereafter, we assume that the given data are two or n pairs of points in correspondence.

The points p 1 = (p 1x , p 1y , p 1z ) and p 2 = (p 2x , p 2y , p 2z ) are said to be in correspondence if p 1 in the rst image corresponds to p 2 in the second image.

In this chapter, we show how to deduce information on the rotation axis and the rotation angle that transform the rst image in the second image from n pairs of points in correspondence. First, we give a method that approximates the rotation axis from n pairs of grid points in correspondence. Secondly, we give a method to obtain the admissible rotation angles knowing a pair of discrete points in correspondence and the rotation axis.

Finding a 3D discrete rotation from a pair of digital images

In this section, we present a method to nd a discrete rotation between a pair of 3D digital images that represents the same object from two dierent points of view. Suppose that we have a set of n grid points A = (⃗ a 1 , ⃗ a 2 , . . . , ⃗ a n ) in the rst image and its corresponding

set B = ( ⃗ b 1 , ⃗ b 2 , . . . , ⃗ b n ) in the second image. Each pair of points (⃗ a i , ⃗ b i ) corresponds to
the same point of the object. We say that such a pair of grid point sets, A and B, are in correspondence.

In 2D, nding a discrete rotation between such a pair of A and B mainly consists in identifying a set of angles that give the same rotated image as explained in Chapter 3.

In Euclidean space, such a rotation is unique. This is not the case in discrete space where two slightly dierent angles α 1 , α 2 may transform A into B. Moreover, for such pair α 1 , α 2 , each angle α 3 such that α 1 ≤ α 3 ≤ α 2 also rotate A into B. Then, we dene the admissible rotation angle, abbreviated by ARA, that is the set of all angles that give the same rotated digital image. It is bounded by the pair of hinge angles α inf and α sup .

In 3D, it is necessary to nd not only a set of admissible rotation angles but also an admissible rotation axis.

We note that we identify a set of admissible rotation angles but only one admissible rotation axis. Indeed, if we consider two rotation axis, the two multi-grids associated to these axes and a point will provide two distinct sets of hinge angles. Therefore, we set a unique rotation axis.

Approximation of the rotation axis by a 3D Pythagorean vector

The determination of an admissible rotation axis from a pair of digital images in correspondence requires two steps. The rst step is to obtain a rotation axis that is consistent with all pairs of points in correspondence. The second step is to approximate the axis obtained through the rst step by a Pythagorean axis in order to obtain rational multigrids.

First step

Given a pair of Euclidean points in correspondence, it is known that any rotation axis such that the rst point is transformed into the second one belongs to the bisection plane of these two points. If we consider two pairs of Euclidean points in correspondence, their rotation axis is the intersection between their two bisection planes. If we add a third pair of Euclidean points in correspondence, the intersection between the three bisection plane still give the same rotation axis. However, if we consider three pairs of grid points in correspondence instead of Euclidean points, the intersection between the three bisection plane may not be unique. In other words, we may have three rotation axes because grid points approximate Euclidean points. Therefore, in order to obtain a unique rotation axis, we consider the rotation axis that is the average of the three obtained rotation axis. Now, we consider voxels instead of Euclidean or grid points. For a given pair of voxels in correspondence, the bisection plane is not unique. Indeed, each pair of Euclidean points that belongs to the pair of voxels in correspondence dene a bisection plane. Therefore, there is an innity of admissible bisection plane for a pair of voxels. Computing the bounds of all the bisection planes for a pair of voxels in correspondence is a big issue. If n pairs of voxels in correspondence are given, and to dene an admissible rotation axis for the n pairs of voxels, we have to compute the intersections between all the obtained bounds of admissible bisection planes. In 2D, a study on the intersection of the bound of bisection for n pairs of pixels have been proposed [START_REF] Rodriguez | Generalized perpendicular bisector and circumcenter[END_REF] and can be extended to the 3D discrete space. However, this work is not nished yet. Therefore, in order to obtain a rotation axis for n pairs of voxels, we consider the n pairs of grid points, which are the centers of voxels, and compute for each pair the bisection plane. Then we compute ⌊ n+1 2 ⌋ axes from the n obtained bisection planes. The rotation axis is the average of the ⌊ n+1 2 ⌋ obtained axes. We strongly believe that this method gives a rotation axis that belongs to the intersections between all the bounds of admissible bisection planes. However, we do not have any formal proof yet.

Second step

In this section, we show a method to approximate the rotation axis obtained in the previous section by a Pythagorean axis that is required to compute a rational multigrid. Based on [START_REF] Anglin | Using pythagorean triangles to approximate angles[END_REF], we can derive a naive method to approximate any 3D vector by a 3D Pythagorean vector; its main idea is to construct two 2D Pythagorean vectors to approximate a 3D Pythagorean vector. For this naive method, some simple denitions are required. We call a Pythagorean triple a set of three integers (i 1 , i 2 , λ) such that

i 2 1 + i 2 2 = λ 2 .
Each Pythagorean triple is associated with the Pythagorean vector (i 1 , i 2 ) ⊤ and a Pythagorean angle θ such that cos θ = i 1 λ . Pythagorean quadruples are the sets of four integers (i

1 , i 2 , i 3 , λ) such that i 2 1 + i 2 2 + i 2 3 = λ 2 .
Each Pythagorean quadruple is associated with the Pythagorean vector (i 1 , i 2 , i 3 ) ⊤ .

The method presented by Anglin in [START_REF] Anglin | Using pythagorean triangles to approximate angles[END_REF] allows the approximation of any angle α with a Pythagorean angle θ such that |α -θ| < ϵ for any ϵ > 0. It uses the theorem, given in [START_REF] Niven | An Introduction to the Theory of Numbers[END_REF], that for each Pythagorean triple (i 1 , i 2 , λ), there is a pair of integers (u, v) such

that v < u, i 1 = 2uv, i 2 = u 2 -v 2 and λ = u 2 + v 2 . Setting X = tan(α -ϵ) + sec(α -ϵ)
and Y = tan(α + ϵ) + sec(α + ϵ), Anglin proved that for every pair of integers (u, v) such that X < u v < Y , the angle θ associated with the Pythagorean triple generated by (u, v) approximates α in such that α -ϵ < θ < α + ϵ.

Let ⃗ r = (r 1 , r 2 , r 3 ) ⊤ be a 3D real vector to be approximated by a Pythagorean vector such that the angle between the two vectors is smaller than ϵ. Such a pair of vectors

is said ϵ-close. From ⃗ r, we dene two 2D vectors ⃗ r 2 , ⃗ r 3 such that ⃗ r 2 = (r 1 , r 2 ) ⊤ , ⃗ r 3 = ( √ r 2 1 + r 2 2 , r 3 ) ⊤ . Each vector ⃗ r i = (x i , y i ) ⊤ is associated with an angle α i satisfying cos α i = x i |⃗ r i | .
Using the algorithm presented in [15] for these two angles α i with a precision ϵ, we obtain two Pythagorean angles associated with the two Pythagorean triples (i 2 , j 2 , λ 2 ), (i 3 , j 3 , λ 3 ). Now, we consider the Pythagorean vectors associated with the two Pythagorean triples. We remark that if (i 2 , j 2 , λ 2 ) is a Pythagorean triple, then (ki 2 , kj 2 , kλ 2 ) is also a Pythagorean triple for any k ∈ Z * . Note that, these two Pythagorean triples are associated with the same Pythagorean angle. This remark allows us to generate from the two Pythagorean triples, (i 2 , j 2 , λ 2 ), (i 3 , j 3 , λ 3 ), a Pythagorean quadruple (ki 2 , kj 2 , lj 3 , lλ 3 ) with two integers k, l such that k(i 2 + j 2 ) = li 3 . The 3D vector ⃗ r ′ associated with this Pythagorean quadruple (ki 2 , kj 2 , lj 3 , lλ 3 ) is our approximation result of ⃗ r.

In [START_REF] Anglin | Using pythagorean triangles to approximate angles[END_REF], the author introduced ϵ that corresponds to the maximum angle between the original vector and the Pythagorean vector that approximates it. With the above method, to construct a 3D Pythagorean vector that approximates a given vector by ϵ, we apply the Anglin's method twice. An example of the method that approximates a 3D vector is given in Figure 5.1. First we decompose the 3D vector ⃗ r to approximate into the two 2D vectors ⃗ r 2 and ⃗ r 3 . ⃗ r 2 and ⃗ r 3 can be approximated by Anglin's method with a precision of ϵ by two 2D Pythagorean vectors ⃗ r ′ 2 and ⃗ r ′ 3 that belong to the two 2D convex cones illustrated in Figure 5.1 (a) and (b), respectively. We generate ⃗ r ′ from ⃗ r ′ 2 and ⃗ r ′ 3 , which is the approximation of ⃗ r, so that it belongs to the 3D convex cone represented in Fig- 2ϵ. Therefore, we can deduce that to reach a precision of ϵ ′ while approximating a 3D vector, we need to give to Anglin's method a precision of

ϵ ≤ ϵ ′ √ 2 .
The algorithm given in this section computes very quickly because, according to Anglin, the computation of a Pythagorean triple is done in constant time. Accordingly, we can give the approximation of any 3D vector in O(1) operations. However, we cannot give any bounds on the size of integers that belongs to the nal Pythagorean quadruple.

(

Approximation of the rotation angle using rational multi-grids

After obtaining an admissible rotation axis, the next step is to obtain the admissible rotation angles. For the 2D case, we have developed, in Section 3.4, a method to obtain from n pairs of grid points the admissible rotation angles where n ≥ 1 in time complexity O(n). This method is designed to work in the 2D half-grid. With minor modications, this method can also be applied on rational multi-grids and keeps the same complexity. a lower bounds such that any angle between them gives the same discrete rotation from A to B. We denote by ARA(⃗ a i , ⃗ b i ) = (α inf i , α sup i ) the pair of hinge angles that gives the lower and the upper bounds of ARA for the pair of points (⃗ a i , ⃗ b i ).

For a given pair of points (⃗ a i , ⃗ b i ) in correspondence, we rst compute the rational multigrid M A,⃗ a i . Then we search for the convexel c containing ⃗ q and compute the two hinge angles between c and the circle centered at the origin of the rational multi-grid and going through ⃗ a i . These two intersections dene the ARA for this pair of points. We note that each operation are done in constant time.

Generally, the given input contains n pairs of points instead of one pair of points. We then incrementally compute the upper and lower bounds of the ARA for these n pairs of points. We rst compute the ARA corresponding to the two rst pair of points (⃗ a 1 , ⃗ b 1 ) and (⃗ a 2 , ⃗ b 2 ). We compare the two pairs of hinge angles (α inf 1 , α sup 1 ) and (α inf 2 , α sup 2 ) obtained and keep the two most restrictive ones so that we have the new pair of hinge angles (max(α inf 1 , α inf 2 ), min(α sup 1 α sup 2 )). We incrementally repeat this operation for the n -2 remaining pairs of points and obtain ARA(A, B) = ∩ i ARA(⃗ a i , ⃗ b i ) for the n pairs of points. Since the comparison of a pair of hinge angles is done in constant time, the time complexity of the incremental algorithm is O(n).

In the above algorithm, we assumed that there are always two intersections between a convexel and the locus of the rotation of ⃗ p. Practically, it may be possible that the locus of rotation of ⃗ p intersects twice a line of M A,⃗ p without intersecting any other line of M A,⃗ p as illustrated in Figure 5.2. In other words, four intersections between a convexel and the rotation locus may exist. We perform experiments in order to evaluate the frequency of such a case. We randomly generate ten thousands Pythagorean axes and for each axis we randomly generate ten thousands grid points. For each grid points ⃗ p, we compute all intersections between M A,⃗ p and the locus of rotation of ⃗ p and search if four of these intersections belong to the same convexel. About 3% of generated points have such particular convexels and the probability that a convexel intersected by the locus of rotation of ⃗ p have four intersections is 1 62500 . Since this particular case does not often occur, we chose not to take it into consideration in our algorithm. Note that in 2D such a case does not happen and the proof can be found in Chapter 3.

Conclusion

In this chapter, we have introduced a method to obtain the rotation angle and rotation axis from a pair of digital images. This work is the extension of the work presented in Section 3.4 of Chapter 3.

In Section 5.2, we introduce the admissible rotation axis that is the set of all possible rotation axis for a pair of digital image. The shape of such a set is not yet studied, however in [START_REF] Rodriguez | Generalized perpendicular bisector and circumcenter[END_REF] the authors studied in discrete plane the shape of admissible center of rotation for n pairs of points. Thus, one of the future work will be to extend it in 3D and then in nD.

Chapter 6

Study on Pythagorean n-tuples

Introduction

In Chapter 4, we introduced a method that performs discrete rotations in 3D using the multi-grid. We explained that to be discrete, this method needs that the rotation axis to be Pythagorean. But in order to approximate all possible 3D rotations it is necessary that Pythagorean vectors, that are the direction vector of Pythagorean axis, are dense in 3D.

The Pythagorean triples are the triple of integers (i 1 , i 2 , λ) such that i 2 1 + i 2 2 = λ 2 . They are well known since Pythagoras, and many studies have been done such as studies on their repartition in the plane [START_REF] Kühleitner | An omega theorem on dierences of two squares[END_REF], how to generate them [START_REF] Teigen | On generating pythagorean triples[END_REF], how to approximate them [START_REF] Anglin | Using pythagorean triangles to approximate angles[END_REF] and the density of their representation [START_REF] Pythagoras | How to found a sect based on three numbers, 520BC[END_REF]. The natural extension of the Pythagorean triples is Pythagorean quadruples that are the quadruple of integers

(i 1 , i 2 , i 3 , λ) such that i 2 1 + i 2 2 + i 2 3 = λ 2
. Some of the studies on Pythagorean triples have been done for Pythagorean quadruples such as generation [START_REF] Cass | Matrix generation of pythagorean n-tuples[END_REF][START_REF] Oliverio | Self-generating pythagorean quadruples and n-tuples[END_REF]. The property of density have been proved arithmetically in [START_REF] Schmutz | Rational points on the unit sphere[END_REF], here we propose a constructive proof that Pythagorean vectors are dense in 3D.

In this chapter, we rst introduce the Pythagorean triples and a sketch of the proof of their density. Then we extend this result to Pythagorean quadruples and to Pythagorean n-tuples. Finally we present a naive algorithm to approximate any sets of (n-1) integers by a Pythagorean n-tuple.

Pythagorean triples

In this section, we rst dene the Pythagorean triples and recall some known properties.

Then we give the proof of the density of Pythagorean triples, which may originally be proved by Pythagoras [START_REF] Pythagoras | How to found a sect based on three numbers, 520BC[END_REF][START_REF] Maor | The Pythagorean theorem: a 4,000-year history[END_REF]. This result will be used in the subsequent sections for the proof of density of Pythagorean quadruples and n-tuples. Denition 6.1. A Pythagorean triple is an integer triple (i 1 , i 2 , λ) such that i 2

1 +i 2 2 = λ 2 .
We denote by P 2 the set of Pythagorean triples such that

P 2 = {(i 1 , i 2 , λ) ∈ Z 3 | i 2 1 + i 2 2 = λ 2 }.
There are two geometrical representations for each element in P 2 : the angle and the vector.

The rst representation of a Pythagorean triple (i 1 , i 2 , λ) can be done by the angle between the line i 1 x -i 2 y = 0 and the abscissa axis (Figure 6.1). Denition 6.2. A Pythagorean angle associated with a Pythagorean triple (i 1 , i 2 , λ) of P 2 is an angle θ such that sin θ = i 2 λ and cos θ = i 1 λ .

Remark 6.3. If a triple of integers (i 1 , i 2 , λ) belongs to P 2 , the triples (ki 1 , ki 2 , kλ) for all k ∈ Z also belong to P 2 ; they are associated with the same Pythagorean angle.

We denote by A P 2 the set of Pythagorean angles. Roughly speaking, Pythagorean angles are the representation of Pythagorean triples in R.

Before showing that A P 2 is dense in R, we give the classic denition of density in R. Denition 6.4. [START_REF] Glenn | Mathematic dictionary, chapter Density[END_REF] A set S is dense in R, if for any element e 1 in R and for any ϵ ∈ R, there is an element e 2 ∈ S such that e 1 < e 2 < e 1 + ϵ.

Regarding Denition 6.4, we have the following theorem.

Theorem 6.5. The set of Pythagorean angles A P 2 is dense in R. One of simple proofs of Theorem 6.5 uses the two following lemmas. Lemma 6.6. [START_REF] Eckert | The group of primitive pythagorean triangles[END_REF] The sum of any two Pythagorean angles is a Pythagorean angle.

In other words, A P 2 is closed by addition. It is easy to admit this lemma, because the sum of two angles consists of additions and multiplications of the sines and cosines of the two given angles. Thus, the sine and cosine of the angle resulting from the addition still belong to rational numbers. Lemma 6.7. For any angle α ∈ [0, 2π], there is an angle θ ∈ [0, 2π] of A P 2 such that θ < α.

The proof of Lemma 6.7 is similar to a proof of the statement, known since Pythagoras, that the number of Pythagorean triples is innity. The idea of this proof is illustrated in Figure 6.2; each time we nd a square of an integer on the third line, we can form a Pythagorean triple with the two corresponding integers on the rst line. Since the third line contains all the odd number squares, we deduce that there is an innity of Pythagorean triples of the form (i 1 , i 2 , i 1 + 1). It is obvious that i 1 > i 2 and lim

i 1 →∞ i 2 i 1 = 0.
Therefore, when i 1 → ∞, the associated angle θ → 0. This induces Lemma 6.7.

The main idea to prove Theorem 6.5 is to show that for any angle α ∈ [0, 2π] and for any ϵ > 0 there is a Pythagorean angle θ such that α < θ < α + ϵ. To obtain such θ, we rst nd a Pythagorean angle θ ′ such that θ ′ < α and a Pythagorean angle γ such that γ < ϵ, thanks to Lemma 6.7. From Lemma 6.6, we know that the angles resulting from the sums θ ′ + γ, θ ′ + 2γ,. . . ,θ ′ + nγ are Pythagorean angles as well. Therefore, we conclude that there is an integer n such that θ = θ ′ + nγ and α < θ < α + ϵ.

The second representation of a Pythagorean triple is a vector as illustrated in Figure 6.1. Denition 6.8. A 2D Pythagorean vector is a vector ⃗ v = (i 1 , i 2 ) ⊤ associated with the element (i 1 , i 2 , λ) of P 2 .

Note that the Euclidean norm of a Pythagorean vector is always an integer. We denote by V P 2 the set of 2D Pythagorean vectors.

To each Pythagorean triple, we have its corresponding Pythagorean angle and Pythagorean vector respectively. Thus, intuitively, we know that V P 2 is dense because A P 2 is dense. To dene the density of a set of such vectors, we need the notion of convex cone, which is dened in [START_REF] Boyd | Convex Optimization[END_REF]. A 2D convex cone in the plane is represented by a pair of linearly independent 2D vectors ⃗ v 1 , ⃗ v 2 bounding of the cone as illustrated in Figure 6.3. Then a vector ⃗ v belongs to the 2D convex cone if there is a pair of scalars

s 1 , s 2 ∈ R * + 1 such that ⃗ v = s 1 ⃗ v 1 + s 2 ⃗ v 2 .
Similarly an nD convex cone in the nD space is represented by a set of n linearly independent nD vectors ⃗ v 1 , ⃗ v 2 , . . . , ⃗ v n bounding the cone as illustrated for the 3D case in Figure 6.5.

The following denition introduces the notion of density for a set of vectors in any dimensions. This denition is needed to prove that Pythagorean triples and Pythagorean n-tuples are dense. Denition 6.9. Let U be a set of nD vectors in R n . We say that U is dense in R n if and only if for the nD convex cone C dened by any set of n linearly independent vectors

⃗ v 1 , ⃗ v 2 , . . . , ⃗ v n ∈ ⃗ R n ,
∈ R * + such that ⃗ u = s 1 ⃗ v 1 + s 2 ⃗ v 2 + • • • + s n ⃗ v n .
With Denition 6.9, it is possible to show that V P 2 is dense in R 2 ; however, this proof is omitted since it is similar to the proof of Theorem 6.5.

Theorem 6.10. The set of 2D Pythagorean vectors V P 2 is dense in R 2 .

Pythagorean quadruples

In Section 6.1, we explained that Pythagorean quadruples can be considered as 3D rational unit vectors. 3D unit vectors appear in many applications in computer vision and computer graphics. Here, we highlight Pythagorean quadruple. For the practical use, those vectors associated with Pythagorean quadruples must be dense on the unit sphere. In this section, we give a constructive proof that they are dense in R 3 using Theorem 6.10. Denition 6.11. A Pythagorean quadruple is an integer quadruple (i

1 , i 2 , i 3 , λ) such that i 2 1 + i 2 2 + i 2 3 = λ 2 .
We dene the set containing all Pythagorean quadruples such that

P 3 = {(i 1 , i 2 , i 3 , λ) ∈ Z 4 | i 2 1 + i 2 2 + i 2 3 = λ 2 }.
Similarly to the case of Pythagorean triples, we present two geometrical representations for each element of P 3 .

Pythagorean quadruples can be represented by angles, similarly to Pythagorean triples.

Each Pythagorean quadruple can be associated with a pair of angles (θ, ψ) where cos θ =

i 1 √ i 2 1 +i 2 2 , sin θ = i 2 √ i 2 1 +i 2 2 and cos ψ = √ i 2 1 +i 2 2 λ
, sin ψ = i 3 λ . Note that, with a Pythagorean triple, we can ensure that the values of the sine and cosine generated by the associated Pythagorean angle belong to Q, and this is not the case with any Pythagorean quadruple because in most cases

√ i 2 1 + i 2 2
does not belong to Z.

We can also represent a Pythagorean quadruple (i 1 , i 2 , i 3 , λ) by a vector ⃗ v = (i 1 , i 2 , i 3 ) ⊤ in the 3D space as seen in Figure 6.4. Denition 6.12. A 3D Pythagorean vector is a vector

⃗ v = (i 1 , i 2 , i 3 ) ⊤ associated with the element (i 1 , i 2 , i 3 , λ) of P 3 .
As well as 2D Pythagorean vectors, each 3D Pythagorean vector has an integer Euclidean norm. We denote by V P 3 the set of 3D Pythagorean vectors. Generally speaking, such geometrical representation for Pythagorean n-tuples is done in (n -1)D space.

Similarly to Pythagorean triples with their 2D geometric representation, the question of density for the representation of Pythagorean quadruples arises. In 2D the proof of the density was made by using Pythagorean angles. However, in 3D, as explained above, the representation by angles is not convenient because it uses irrational numbers. Therefore,

⃗ v ′ 1 , ⃗ v ′ 2 , ⃗ v ′ 3 on the OXY plane, and (b) the 2D convex cone C 2 dened by {⃗ v ′ 1 , ⃗ v ′ 2 }.
we will prove the density of 3D Pythagorean vectors. The denition of density used here is the one given in Denition 6.9.

Theorem 6.13. The set of 3D Pythagorean vectors V P 3 is dense in R 3 .

Proof. 11. To prove Theorem 6.13, it is sucient to show that for every set of three vectors, we can construct a Pythagorean vector belonging to the 3D convex cone dened by the three vectors.

Let ⃗ v 1 , ⃗ v 2 , ⃗ v 3 be three linearly independent 3D vectors. We denote by C 3 the 3D convex cone dened by position vectors ⃗ v 1 , ⃗ v 2 , ⃗ v 3 . We also denote by

⃗ v ′ 1 , ⃗ v ′ 2 , ⃗ v ′ 3 the projections of ⃗ v 1 , ⃗ v 2 , ⃗ v 3 in
the OXY plane (Figure 6.6(a)).

Because ⃗ v 1 , ⃗ v 2 , ⃗ v 3 are linearly independents, at least one of these three pairs

{⃗ v ′ 1 , ⃗ v ′ 2 }, {⃗ v ′ 1 , ⃗ v ′ 3 }, {⃗ v ′ 2 , ⃗ v ′ 3 } is linearly independent in the plane OXY . Let us assume that the pair of vectors {⃗ v ′ 1 , ⃗ v ′ 2 } is linearly independent. Thus {⃗ v ′ 1 , ⃗ v ′ 2 }
denes in the OXY plane a 2D convex cone, which we denote by C 2 as shown in Figure 6.6(b). Thanks to Theorem 6.10, there is a 2D Pythagorean vector ⃗ u ′ that belongs to C 2 (Figure 6.6(b)). ⃗ u ′ is associated with the Pythagorean triple (i 1 , i 2 , λ ⃗ u ). We set the vector ⃗ u = (i 1 , i 2 , 0) ⊤ .

Let P be a plane dened by the two vectors ⃗ u and (0, 0, 1) ⊤ and the point (0, 0, 0) ⊤ (Figure 6 cone C ′ 2 that belongs to P (Figure 6.7(b)). Thanks to Theorem 6.10 again, there is a 2D Pythagorean vector ⃗ w ′ that belongs to C ′ 2 (Figure 6.7(b)); ⃗ w ′ is associated with the Pythagorean triple (j 1 , j 2 , λ j ).

As said in Remark 6.

Approximation of a given 3D vector by a Pythagorean 3D vector

In this section we explain how to approximate a given vector by a Pythagorean vector.

Because we proved in Section 6.3 that 3D Pythagorean vectors are dense, we know that we nd a Pythagorean vector as close as possible to the given vector considering the distance d between the two vectors

⃗ u = (u x , u y , u z ) ⊤ , ⃗ v = (v x , v y , v z ) ⊤ as follow d(⃗ v, ⃗ u) = √ ( ux ∥⃗ u∥ -vx ∥⃗ v∥ ) 2 + ( uy ∥⃗ u∥ - vy ∥⃗ v∥ ) 2 + ( uz ∥⃗ u∥ -vz ∥⃗ v∥ ) 2 . Two vectors are said to be ⃗ v, ⃗ u are ϵ-close if d(⃗ v, ⃗ u) < ϵ.
In this section, we rst explain how to generate all Pythagorean quadruple, then we give two methods to approximate a given vector by an ϵ-close Pythagorean vector. Finally, we show some experiment on approximation using both methods.

Generation of Pythagorean quadruples

In this section we explain how to obtain all Pythagorean quadruples using twelve matrices. This method is the extension to the Pythagorean quadruples of the method used for Pythagorean triples and presented in many papers [START_REF] Teigen | On generating pythagorean triples[END_REF][START_REF] Cass | Matrix generation of pythagorean n-tuples[END_REF][START_REF] Eckert | The group of primitive pythagorean triangles[END_REF]. Here we will rst describe the method used to generate Pythagorean triples, then we will extend it to Pythagorean quadruple and show what problems appear during the generation of the tree.

Note that in this section we only consider Pythagorean triples (respectively Pythagorean quadruples) of the form

< i 1 , i 2 , λ > (respectively < i 1 , i 2 , i 3 , λ >) where 0 < i 1 ≤ i 2 ≤ λ (respectively 0 ≤ i 1 ≤ i 2 ≤ i 3 ≤ λ).
This assumption does not restrict the Pythagorean triples (respectively quadruples) generated, since all other Pythagorean triples (quadruples) can be obtained by switching two values or multiplying some values by -1.

It is generally admitted that the smallest Pythagorean triple is < 3, 4, 5 >. We say that the Pythagorean triple < i 11 , i 12 , λ 1 > is smaller than the Pythagorean triple < i 21 , i 22 , λ 2 > if λ 1 < λ 2 or if λ 1 = λ 2 and i 12 < i 22 . With this order, we can also consider the triple < 0, 1, 1 > as the smallest Pythagorean triple, but this Pythagorean triple is degenerated since it contains a 0.

The method to generate the tree that contains all Pythagorean triples uses the three following matrices: An important property of the construction is that there is a proof that each Pythagorean triple in the tree is a prime Pythagorean triple and each Pythagorean triple appears only in the tree.

M 1 =     1 2 2 2 1 2 2 2 3     , M 2 =     -1 2 2 -2 1 2 -2 2 3     , M 3 =     1 -2 2 2 -1 2 2 -2 3     . ( 6
Note that there is also a fourth matrix:

M 4 =     -1 -2 2 -2 -1 2 -2 -2 3     (6.2)
that gives the father of a Pythagorean triple in the tree.

We expect that there is a Pythagorean tree for Pythagorean quadruples that has the same properties than the tree of a Pythagorean triple, however, we will see that this is not the case. Now we consider Pythagorean quadruples. The naive adaptation of the above method to obtain the Pythagorean quadruple tree is to add a dimension to matrices 6.1 and to apply it on all quadruples (i 1 , i 2 , i 3 , λ), (i 2 , i 1 , i 3 , λ) and (i 3 , i 2 , i 1 , λ) to obtain all sons. As the smallest Pythagorean quadruple, there are two candidates < 0, 0, 1, 1 > and < 1, 2, 2, 3 >. The Pythagorean quadruple < 0, 0, 1, 1 > can be considered as degenerated, thus the naive solution is to consider the smallest Pythagorean quadruple as < 1, 2, 2, 3 >. However, in the Pythagorean triples case, applying the three matrices 6.1 to < 0, 1, 1 > always gives < 3, 4, 5 >, in the Pythagorean quadruple case, application of Matrix 6.3 on < 0, 0, 1, 1 > will give two dierent Pythagorean quadruples: < 1, 2, 2, 3 > and < 0, 3, 4, 5 >. Even if < 0, 3, 4, 5 > is degenerated, some of its sons are not, thus we have to consider it to create the Pythagorean quadruple tree.

The extension of Matrix 6.1 to the 4D gives us the following matrix:

M 1 =        1 0 0 0 0 1 2 2 0 2 1 2 0 2 2 3        , M 2 =        1 0 0 0 0 -1 2 2 0 -2 1 2 0 -2 2 3        , M 3 =        1 0 0 0 0 1 -2 2 0 2 -1 2 0 2 -2 3        . ( 6.3)
In the Pythagorean triples case, we just have to apply each of these matrices once to obtain all possible sons. However, to obtain all sons of a Pythagorean triple, it is necessary to apply Matrix 6.3 and then to switch i 1 and i 2 and apply again Matrix 6.3. We repeat this operation switching i 1 and i 3 . The switch operation is needed because Matrix 6.3

does not modify the rst value of the Pythagorean quadruple. Therefore the switching operation allows us to obtain nine sons.

When we generate the Pythagorean quadruples tree using the naive method, some problems arise. As explained before, we should consider < 0, 0, 1, 1 > as the smallest Pythagorean quadruple even if its a degenerated one. However, considering this quadruple as the root of the tree will remove the guarantee that each Pythagorean quadruple only appears in the tree. Applying the Matrix 6. Pythagorean quadruples. There is an innity of Pythagorean quadruples that cannot be reached using this method. The solution of this problem is to use the fourth matrix From all these problems, we can derive a new method of generation of the Pythagorean tree that removes all repetitions. We set the root as < 0, 0, 1, 1 > and use the four matrices M 1 , M 2 , M 3 , M 4 . The rst level of the tree contains only the two Pythagorean quadruples: < 0, 3, 4, 5 > and < 1, 2, 2, 3 >. To obtain the next level, we proceed as follow. For each Pythagorean quadruple of the last level, we compute the twelve sons that correspond to the three applications of the four matrices. For each son that results on the application of M 4 we check if it is the same as its grandfather, if so we remove it. For the other sons, we remove all the repetitions.

M 4 =        1 0 0 0 0 -1 -2 2 0 -2 -1 2 0 -2 -2 3        , ( 6

Approximation methods

In this section, we present two methods to approximate a given vector ⃗ v = (v x , v y , v z ) ⊤ by an ϵ-close Pythagorean vector. In few words, the rst method constructs the Pythagorean tree level by level and checks in each level if there is a Pythagorean vector that is ϵclose to the given vector. The second one selects for each level, the three Pythagorean

Experiments

In this section, we present some experiment on the two methods presented in Section 6.4.2. In order to evaluate our method, we generate for both methods n vectors to approximate. For the rst method, we dene a depth d and compute the corresponding Pythagorean tree. Then, for each vector to approximate, we compute the distance between this vector and all Pythagorean vectors in the tree and keep the lowest distance.

Then we compute the average distance for all generated vectors. For the second method, we dene a maximum depth d, and compute incrementally the d convex cones. We keep the closest Pythagorean vector computed. For the rst method, we set d = 4, 5, 6, 7, 8 and generated 1000 vectors to approximate.

These experiments are illustrated in Figure 6.12(a-e). For the second method, we set d = 10, 20, 50, 100, 200, 300 and use the same 1000 vectors to approximate than the rst method. These experiments are illustrated in Figure 6.13(af). For all the experiments, computations are fast, approximately 5ms when d = 10 to 0.15s when d = 300 for each Pythagorean vector. The computation time is linear regarding to the depth of the tree. The average distance between the original vector and its approximation evolves from 0.046 when d = 10 to 0.042 when d = 300. These average distances are not really improved between d = 10 and d = 300. The average distance obtained by the second method is always worse than the average distance obtained by the rst method for all depth of the tree. However, even if the average distance is bad, many vectors are well approximated. The analysis of the result shows that many vectors that are not correctly approximated by the rst method when d = 4, 5 obtain a better approximation using the second method for d = 10.

The second method cannot be used without the rst one since many Pythagorean vectors are not well approximated. However, it is interesting to use the second method in parallel of the rst one when d is high since the rst method has a high computation time and the second method can give a good result. 15. An nD Pythagorean vector is a vector ⃗ v = (i 1 , i 2 , . . . , i n ) ⊤ associated with the element (i 1 , i 2 , . . . , i n , λ) of P n .

The Euclidean norm of an nD Pythagorean vector is an integer. We denote by V Pn the set of nD Pythagorean vectors and we are going to prove the density of this set. The following theorem is the generalization of Theorems 6.10 and 6.13.

Theorem 6.16. The set of nD Pythagorean vectors V Pn is dense in R n .

Since we already proved the case of n = 2, 3 in Theorems 6.10 and 6.13, we give the proof for n ≥ 4 using the induction. The proof for n ≥ 4 is similar to the proof for n = 3. The main idea is to construct an nD Pythagorean vector in a given nD convex cone.

Proof. 12. First we assume that the set V P n-1 is dense in R n-1 .

Let ⃗ v 1 , ⃗ v 2 , . . . , ⃗ v n be n linearly independent nD vectors. We denote by C n the nD convex cone dened by position vectors ⃗ v 1 , ⃗ v 2 , . . . , ⃗ v n .

Let ⃗ v ′ 1 , ⃗ v ′ 2 , . . . , ⃗ v ′ n be the projections of ⃗ v 1 , ⃗ v 2 , . . . , ⃗ v n on the hyper-plane H orthogonal to the n-th axis, respectively. Because ⃗ v 1 , ⃗ v 2 , . . . , ⃗ v n are linearly independent in nD, there is j ∈ {1, . . . , n} such that ⃗ v ′ 1 , ⃗ v ′ 2 , . . . , ⃗ v ′ j-1 , ⃗ v ′ j+1 , . . . , ⃗ v ′ n are linearly independent in H. Thus, in H, the set of vectors ⃗ v ′ 1 , ⃗ v ′ 2 , . . . , ⃗ v ′ j-1 , ⃗ v ′ j+1 , . . . , ⃗ v ′ n denes an (n -1)D convex cone, which we denote by C n-1 . Because the set V P n-1 is dense in (n -1)D, there is an (n -1)D Pythagorean vector ⃗ u ′ that belongs to C n-1 such that ⃗ u ′ is associated with the Pythagorean n-tuple (i 1 , i 2 , . . . , i n-1 , λ ⃗ u ). We denote by ⃗ u the nD vector (i 1 , i 2 , . . . , i n-1 , 0) ⊤ .

Let P be a 2D plane dened by the two directional nD vectors: (0, . . . , 0, 1) ⊤ and ⃗ u and the point (0, . . . , 0) ⊤ . By the construction, P intersects C n . This intersection is a 2D convex cone, which we denote by C 2 that belongs to P. Thanks to Theorem 6.10, there is a 2D Pythagorean vector associated with a Pythagorean triple (j 1 , j 2 , λ ⃗ w ) that belongs to C 2 .

As said in Remark 6.3, if we multiply all values of a Pythagorean triple by k ∈ Z, the obtained triple is still a Pythagorean triple. Note that Remark 6.3 is still valid for Pythagorean k-tuples even for k ≥ 3. Therefore, it is possible to design from a Pythagorean triple and a Pythagorean n-tuple a Pythagorean (n + 1)-tuple as follows.

Let a, b ∈ Z * such that aλ ⃗ u = bj 1 . Then we obtain the Pythagorean (n + 1)-tuple (ai 1 , ai 2 , . . . , ai n-1 , bj 2 , bλ ⃗ w ). By construction the vector ⃗ u = (ai 1 , ai 2 , . . . , ai n-1 , bj 2 ) ⊤ , associated with the Pythagorean (n + 1)-tuple (ai 1 , ai 2 , . . . , ai n-1 , bj 2 , bλ ⃗ w ), belongs to the nD convex cone C n .

Approximation of an nD vector by an nD Pythagorean vector

In this section, we explain how to approximate a given nD vector by an nD Pythagorean vector. Our method is based on the algorithm proposed by Anglin in [START_REF] Anglin | Using pythagorean triangles to approximate angles[END_REF] that allows us to approximate any angle α with a precision ϵ by a Pythagorean triple and its associated angle θ. It uses the theorem, given in [START_REF] Niven | An Introduction to the Theory of Numbers[END_REF], that for each Pythagorean triple (i 1 , i 2 , λ), there is a pair of integers (u, v) such that v < u, i 1 = 2uv, i 2 = u 2 -v 2 and λ = u 2 + v 2 .

Setting X = tan(α -ϵ) + sec(α -ϵ) and Y = tan(α + ϵ) + sec(α + ϵ), Anglin proved that for every pair of integers (u, v) such that X < u v < Y , the angle θ associated with the Pythagorean triple generated by (u, v) approximates α such that α -ϵ < θ < α + ϵ.

Based on [START_REF] Anglin | Using pythagorean triangles to approximate angles[END_REF], we can derive a naive method to approximate any nD vector by an nD Pythagorean vector; its main idea, which is similar to the previous proof of density, is to construct m 2D Pythagorean vectors to approximate an nD Pythagorean vector.

Let ⃗ v = (r 1 , r 2 , . . . , r n ) ⊤ be an nD real vector to be approximated by a Pythagorean vector that is ϵ close 2 . From ⃗ v, we dene (n -1)

2D vectors ⃗ v 2 , ⃗ v 3 , . . . , ⃗ v n such that ⃗ v 2 = (r 1 , r 2 ) ⊤ , ⃗ v 3 = ( √ r 2 1 + r 2 2 , r 3 ) ⊤ , . . . , ⃗ v n = ( √ r 2 1 + r 2 2 + • • • + r 2 n-1 , r n ) ⊤ .
Each vector ⃗ v i = (x i , y i ) ⊤ is associated with an angle α i satisfying cos α i = x i |⃗ v i | . Using the algorithm presented in [START_REF] Anglin | Using pythagorean triangles to approximate angles[END_REF] for each of the (n -1) angles α i with a precision ϵ, we obtain (n -1) Pythagorean angles associated with the (n -1) Pythagorean triples (i 2 , j 2 , λ 2 ), (i 3 , j 3 , λ 3 ), . . . , (i n , j n , λ n ). Now, we consider the Pythagorean vectors associated with the n Pythagorean triples. Remark 6.3 allows us to generate from the two rst Pythagorean triples (i 2 , j 2 , λ 2 ), (i 3 , j 3 , λ 3 ), a Pythagorean quadruple (k 1 i 2 , k 1 j 2 , l 1 j 3 , l 1 λ 3 ) with two integers k 1 , l 1 such that k 1 (i 2 + j 2 ) = l 1 i 3 . Similarly, from this Pythagorean quadruple and the third Pythagorean triple (i 4 , j 4 , λ 4 ), we can generate the Pythagorean quintuple (k 1 k 2 i 2 , k 1 k 2 j 2 , l 1 k 2 j 3 , l 2 j 4 , l 2 λ 4 ) with two integers k 2 , l 2 such that k 2 (k 1 i 2 + k 1 j 2 + l 1 j 3 ) = l 2 i 4 . We repeat this combination operation between a Pythagorean (i + 2)-tuple and a Pythagorean triple by adding a pair of integers k i and l i until obtaining a Pythagorean (n + 1)-tuple. The nD vector associated with this Pythagorean (n + 1)-tuple is our approximation result of ⃗ v.

In [START_REF] Anglin | Using pythagorean triangles to approximate angles[END_REF], the author introduced ϵ that corresponds to the maximum angle between the original vector and the Pythagorean vector that approximates it. With the above method, to construct an nD Pythagorean vector that approximates a given vector by ϵ, we apply (n -1) times Anglin's method. An example of the method that approximates a 3D vector is given in Figure 5.1. First we decompose the 3D vector ⃗ v to approximate into the two 2D vectors ⃗ v 2 and ⃗ v 3 . ⃗ v 2 and ⃗ v 3 can be approximated by Anglin's method with 2 Two vectors are ϵ close if the angle between them is smaller than ϵ. a precision of ϵ by two 2D Pythagorean vectors ⃗ v ′ 2 and ⃗ v ′ 3 that belong to the two 2D convex cone illustrated in Figure 5.1 (a) and (b) respectively. Thanks to Remark 6.3, we generate ⃗ v ′ from ⃗ v ′ 2 and ⃗ v ′ 3 , which is the approximation of ⃗ v, such that it belongs to the 3D convex cone represented in Figure 5.1 (c) forming a square pyramid. The minimum circular cone including the square pyramid has the solid angle 2 √ 2ϵ. Therefore, we can deduce that to reach a precision of ϵ ′ while approximating a 3D vector, we need to give to Anglin's method a precision of ϵ ≤ ϵ ′ √ 2 . In general, to reach a precision of ϵ ′ while approximating an nD vector, it is necessary to give to Anglin's method a precision of ϵ ≤ ϵ ′ √ n for each 2D vector to approximate.

The algorithm given in this section is fast because, according to Anglin, the computation of a Pythagorean triple is done in constant time. Accordingly, we can give the approximation of any nD vector in O(n) operations. However, we cannot give any bounds on the size of integers that belong to the nal Pythagorean (n + 1)-tuple.

Conclusion

In this chapter, we introduce the Pythagorean n-tuple that are required for hinge angles in 3D. We gave a constructive proof of the density on 3D Pythagorean vector and then on nD Pythagorean vectors. This proof is required to show that for any given vector, we can construct a Pythagorean vector that approximates the given vector.

The algorithm presented in Section 6.6 does not give an upper bound for the size of integers composing the Pythagorean n-tuple. A future work is to give an upper bound.

For reaching this goal, we have two directions. One is to use a method that computes all primitive Pythagorean n-tuples using a tree structure [START_REF] Cass | Matrix generation of pythagorean n-tuples[END_REF]. The other is to study the repartition of nD Pythagorean vectors in the n dimensional space.

Chapter 7 Conclusion

As explained in the introduction, the main goal of this thesis was to improve the understanding of the rotations in the 2D and 3D discrete space and to develop new tools to compute discrete rotations using only integers. To conclude this manuscript, we resume the main result obtained during my thesis and introduce some future works and developments.

Discrete rotation in 2D

During my thesis, we studied the particularity of the 2D rotations in discrete space. In Chapter 2, we present a survey on 2D discrete rotations and the problems that arise during computations. We show that each existing discrete rotation solves some of the problems presented, but not all. We also present a study on rotations in the two other regular grids: triangle and hexagonal grid. We show that the rotations in hexagonal grid give better results than rotation in orthogonal grid regarding points lost and connectivity. However, even if the results obtained are better, the properties of bijectivity and transitivity are lost and the rotation is not discrete regarding our denition. Then one of the future work will be to develop a new discrete rotation based on hexagonal grid that improves the results on points lost and connectivity presented in Chapter 2. A possible method to obtain such results is to mix the results obtained with a rotation on hexagonal and orthogonal grids. Another possible method is to adapt the rotation by circle to the hexagonal grid.

Another property, which is not discussed in this manuscript, is the conservation of topology. To sum up, topology is the mathematical study of the properties that are preserved through deformations, twisting, and stretching of objects. A coee cup and a donuts are topologically equivalent since they both have one hole but they are not equivalent to a glass that has no hole. Obviously, Euclidean rotation keeps the topology of an object and discrete rotation does not keep it since it is required to keep connectivity and to not create holes while losing points. Another future work will be to design a discrete rotation that ensures that the topology of object is kept during discrete rotations.

In Chapter 3, we studied hinge angles and their applications to discrete rotations. Hinge angles represent the discontinuities of the rotations on discrete space. They allow the design of a discrete rotation having the same results than the discretized Euclidean rotation that is the Euclidean rotation applied on a discrete space. The 2D discrete rotation using hinge angles is neither bijective or transitive. A future work is to try to give these two properties to this rotation. A method based on multi-scale can give such results. A promising idea is to divide the pixel into many little pixels with dierent weight. The discretization will be done regarding the number of each little pixel in the normal pixel and their respective weight.

Discrete rotation in 3D

As we explained in the rst section of Chapter 4, discrete rotations in 3D are not as well studied as discrete rotations in 2D. There was at our best knowledge, no rotation that is discrete under our denition. Then, in Chapter 4, we have designed a 3D discrete rotation based on hinge angles. For this rotation, we have extended the hinge angles in 3D, dened the multi-grids that are the intersection between the rotation plane and the 3D half-grid and present a study on Pythagorean quadruple. This discrete rotation gave the same results than the discretization of a 3D continuous rotation around a given axis.

From this discrete rotation, we have extended the research of a rotation done between a pair of digital images in correspondence presented in Chapter 3 to the 3D discrete space. The extension done in Chapter 5 searches for an admissible axis of rotation and x it in order to obtain admissible rotation angle regarding to this axis. However, another approach is to obtain all admissible rotation axis and to obtain rotation angles regarding to the rotation axis. To obtain such results, we have to compute the bound of all admissible rotation axis for the given n pairs of points in correspondence. In 2D, such bounds have been computed for an admissible center of rotation [START_REF] Rodriguez | Generalized perpendicular bisector and circumcenter[END_REF]. The extension of this work to 3D is required to obtain the bound of the admissible rotation axis. As explain in above section, a future work will be to design a discrete rotation that keep the topology. Such a rotation is also required for 3D space. However, in 3D a new artifact called tunnel, appear in 3D objects. Thus the simple extension of a 2D discrete rotation that keep topology to 3D may not keep topology in 3D.

In 2D we have presented a short study on rotation in dierent regular grids. In 3D, according to the Russian mathematician Fyodorov, the only regular polygon that tile the plane is the cube. However, some other semi-regulars polygons tile the plan such as truncated octahedron or rhombic dodecahedron. The shape of these two polygons is closer to the shape of a sphere than the shape of the cube. Therefore, study the rotation in 3D space tiled by one of these polygons may give interesting result.

Another future work will be to design a 3D rotation that extends the rotation by discrete circles presented in Chapter 2. The rotation by discrete spheres will be a bijective 3D discrete rotation. A study is required to determine which denition of discrete spheres will give the best results. This denition requires at least to tile the 3D discrete space.

Discrete rotation in nD

Commonly, the rotation is required in 2D or in 3D. However, for some particular elds, rotations in higher dimension can be required such as color image analysis where 4D rotations are needed or global visibility where 5D or 6D rotation are needed [START_REF] Batard | Cliord fourier transform for color image processing[END_REF] 1 . In 2D, the rotations are computed in a unique plane. In 3D, to each rotation is associated a unique rotation axis and a unique plane, denoted by rotation plane in this thesis, to compute the rotation. It is important to notice that in 2D and 3D, two given rotations with dierent parameters (angle of rotation and/or axis of rotation) will give dierent result. When the dimension is higher than three, this is not necessarily true. Indeed, when the dimension n is higher than three, the number of angles required to perform a rotation is ⌊ n 2 ⌋. Moreover, to perform a rotation in dimension n, we need to compute ⌊ n 2 ⌋ rotation planes that are associated to the ⌊ n 2 ⌋ angles. Two dierent sets of rotation planes and angles can dene the same rotation. Another remark, which is strongly related to the previous one, is that in n dimensions, a rotation with a unique angle will leave (n-2) stationary dimensions. Then in a planar rotation we have a 0D stationary dimension, in other word a point, in 3D we have a stationary axis, in 4D we have a stationary plane and so on.

In nD some continuous method to perform rotation are known such as method using Cliord algebra. At our best knowledge, there exists no method to perform discrete rotation in nD. The hinge angles and multi-grids presented in this thesis can be extended to perform rotation in any dimension. However, some problems appear. Firstly, as explain above, a rotation in nD must be decomposed in ⌊ n 2 ⌋ planar rotation. In 3D, we know that this decomposition in unique, in nD the question is open. Another problem is the uniqueness of representation of hinge angle by n-tuple of integers. In 2D and 3D,

1. 1

 1 Intersection between two circles and between their discretization. The discrete circles have four pixels for the intersection while the continuous circles have only one point of intersection. That shows one of the main dierences between Euclidean geometry and digital geometry. . . . . . . . 2.1 The four main problems of discrete rotation algorithms: (a) double point, (b) blank point, (c) errors due to composition and (d) inverse rotation. . . 2.2 The DER rotation of a white image by the four angles (a)α = 15 • , (b)α = 30 • , (c)α = sin -1 3 5 ≈ 36.87 • and (d)α = 37 • . . . . . . . . . . . . . . . . . 2.3 Percent of point lost during the rotation of a 200 × 200 picture for all angles between [0, π 2 ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 During the composition of DER, the number of points lost due to double point is cumulative. For an angle α = 40 • we lost approximately 15% of points for (a), 28% for (b), 39% for (c) and 50% for (d). . . . . . . . . . . 2.5 If a discrete rotation is DER-like, then the red way that discretize then rotate will gives the same result than the blue way that rotate then discretize. 2.6 The rotation by circles by an angle of 45 • . . . . . . . . . . . . . . . . . . . 2.7 The twelve biggest angles that can possibly be obtained by Pythagorean lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3 . 3

 33 Illustration of cases ΣF (C i (⃗ q 2 )) = 0,1,2 or 3. . . . . . . . . . . . . . . . . 3.4 Running of the incremental algorithm. . . . . . . . . . . . . . . . . . . . . 3.5 Randomly generated points and their corresponding rotated points. . . . . 3.6 Result of our algorithm applied on the sets of points in Figure 3.5 . . . . . 3.7 Result of our algorithm applied on a set of oating points. . . . . . . . . . 3.8 Examples of errors for computing lower and upper bounds of rotation angles introduced by rotation of oating points. . . . . . . . . . . . . . . . 3.9 Result of our algorithm applied on real data. . . . . . . . . . . . . . . . . 3.10 Problem of corresponding points with real data. . . . . . . . . . . . . . . . 4.1 All hinge angles in the rst quadrant for the grid point ⃗ p = (2, 1) ⊤ , such that α 1 ≈ -13.68 • , α 2 ≈ 15.59 • , α 3 ≈ 21.31 • , α 4 ≈ 50.52 • . The hinge angle α 5 is obtained by symmetry such that α 5 = π 2 -α 1 ≈ 76.32 • . . . . 4.2 The 3D half-grid cut by a plane (a), and its multi-grid (b). . . . . . . . . . 4.3 Parallel lines of a set L A,⃗ p i and geometric interpretation of their parameters 4.

4 . 7

 47 the two impossible cases of intersection between a circle centered on a grid point and a rational multi-grid. . . . . . . . . . . . . . . . . . . . . . . . . 4.8 A hinge angle α for a point ⃗ p in a rational multi-grid. . . . . . . . . . . . . 5.1 Example of the approximation of a 3D vector using our method. (a) and (b) represent the projections of ⃗ r into ⃗ r 2 and ⃗ r 3 in 2D planes respectively. Each 2D convex cone in (a) and (b) is the admissible approximation for ⃗ r 2 and ⃗ r 3 regarding to ϵ. (c) represents the square pyramid constructed from the two 2D convex cones that contains the sets of admissible approximation for ⃗ r. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 In a rational multi-grid M A,⃗ p a convexels can have four intersections with the locus of rotation of ⃗ p. . . . . . . . . . . . . . . . . . . . . . . . . . . .

6 . 8 1 Introduction 1 . 1

 68111 the OXY plane, and (b) the 2D convex cone C 2 dened by {⃗ v ′ 1 , ⃗ v ′ 2 }. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 (a) The plane P dened by OZ and ⃗ u, and (b) the vector ⃗ w that belongs to C ′ 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi The three rst level of the Pythagorean tree. . . . . . . . . . . . . . . . . 6.9 The Pythagorean quadruple < 23, 24, 24, 41 > can only be reached using Matrix 6.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.10 In most cases, M 4 gives the father of a Pythagorean quadruple. . . . . . . 6.11 There exists an innity of Pythagorean quadruples that have themselves as grandfather. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.12 (a) Distance between 1000 vectors and their approximated Pythagorean vector for a Pythagorean tree of depth 4 (a), of depth 5 (b), of depth 6 (c), of depth 7 (d) and of depth 8 (e) . . . . . . . . . . . . . . . . . . . . . 6.13 Distance between 1000 vectors and their approximated Pythagorean vector for a partial Pythagorean tree of depth 10 (a), of depth 20 (b), of depth 50 (c), of depth 100 (d), of depth 200 (e) and of depth 300 (f) . . . 6.14 Distance between 1000 vectors and their approximated Pythagorean vector for a random Pythagorean tree of depth 10 (a), of depth 20 (b), of depth 50 (c), of depth 100 (d), of depth 200 (e) and of depth 300 (f) . . . Chapter Introduction to the rotation problem in Z 2

Figure 1 . 1 :

 11 Figure 1.1: Intersection between two circles and between their discretization. The discrete circles have four pixels for the intersection while the continuous circles have only one point of intersection. That shows one of the main dierences between Euclidean geometry and digital geometry.

Figure 1 .

 1 Figure 1.1 gives a simple example of the dierence between the Euclidean geometry and the digital geometry. In Euclidean geometry, two circles have zero, one, two or an innity of intersections (if they are merged). However, in digital geometry, two discrete circles can have any number of common pixels (four common pixels on Figure 1.1) but not

Figure 2 .

 2 1(a). During the DER, is it possible that two points, which are 4-connected, have their images in the same pixel. Then during the discretization of the rotation, these two points are discretized in the same pixel. To give an idea of this problem, we must notice that the Euclidean1 Floating numbers are the representation in computer science of the real numbers. They have a limited precision that implies an approximation of the real numbers.

Figure 2 . 1 :

 21 Figure 2.1: The four main problems of discrete rotation algorithms: (a) double point, (b) blank point, (c) errors due to composition and (d) inverse rotation.

2

 2

Figure 2 . 2 :Figure 2 . 3 :

 2223 Figure 2.2: The DER rotation of a white image by the four angles (a)α = 15 • , (b)α = 30 • , (c)α = sin -1 3 5 ≈ 36.87 • and (d)α = 37 •

  (a),(b),(c) and (d). A short study of these patterns have been done in [4]. Figures 2.2(c) shows that some particular angles do not create any blank points. These particular angles are used to create a discrete rotation presented Section 2.3.2.

Denition 2 . 4 .

 24 A discrete rotation is transitive if for all angles α and γ and for all points ⃗ p the rotations of ⃗ p by α then by γ give the same result than the rotation of ⃗ p by α + γ. Note that, the problems of double points and blank points are cumulative with the problem of compositions. An example is given in Figure 2.4. These gures represent the four rotations of a white image of size 450 × 450 by an angle of α = 40 • . Figure 2.4(a)is the result of the rst rotation of the white image by α.

Figure 2 .Figure 2 . 4 :Figure 2 . 5

 22425 Figure 2.4(c) and Figure 2.4(d). After the last rotation, it remains approximately 50% of the white points of the original image.

Figure 2 . 6 :

 26 Figure 2.6: The rotation by circles by an angle of 45 • .

Figure 2 . 7 :

 27 Figure 2.7: The twelve biggest angles that can possibly be obtained by Pythagorean lines

Figure 2 . 8 :

 28 Figure 2.8: The quasi-shear rotation and its two intermediate results.

  (a) and (b). Two pixels are 4-connected if they share an edge and they are 8-connected if they share an edge or a vertex. If we consider the hexagonal grid, two pixels are 6-connected or not connected as illustrated on Figure 2.10. Regarding the 3D case, we only consider the orthogonal grid. In this case, the variable m can take the values of 6, 18 or 26 as illustrated on Figures 2.11

  (a) (b) and (c). Two voxels are 6-connected if they share a face, they are 18-connected if they share a face or

Figure 2 .

 2 Figure 2.10: A pixel and its six 6-connected pixels in a hexagonal grid

Figure 2 .

 2 Figure 2.11: A voxel and (a) its six 6-connected voxels, (b) its eighteen 18-connected voxels and (c) its twenty-six 26-connected voxels in a cubic grid.

Figure 2 . 12 :

 212 Figure 2.12: Conversion from a square grid into a triangle grid

Figure 2 . 13 :

 213 Figure 2.13: Conversion from a square grid into a hexagonal grid

Figure 2 . 14 :

 214 Figure 2.14: The percentage of points remaining after a rotation in a triangular grid.

Figure 2 .

 2 Figure2.14 presents the results obtained on a rotation using the triangle grid. We compare them with the results obtained in Figure2.3. First, the rotation using triangle grid always lose at least 13% of points after a rotation. This result is globally worse than the results obtained on a square grid. However, for angles between 32 • and 58 • , the results are almost the same than the results obtained on a square grid.

Figures 2. 15

 15 Figures 2.15(a) and (b) show the average conservation of neighborhood depending on

  Figures 2.15(a) and (b) show the average conservation of neighborhood depending on the rotation angle on a triangle grid and compare them with the result obtained on

Figure 2 . 17 :

 217 Figure 2.17: Result of four rotations of the image (a) by an angle of 15 • (b) using a square grid, (c) using a hexagonal grid and (d) using a triangle grid.

Figure 2 . 18 :

 218 Figure 2.18: Result of four rotations of the image (a) by an angle of 30 • (b) using a square grid, (c) using a hexagonal grid and (d) using a triangle grid.

Figure 2 . 19 :

 219 Figure 2.19: Result of four rotations of the image (a) by an angle of 45 • (b) using a square grid, (c) using a hexagonal grid and (d) using a triangle grid.

. 2 )Figure 3 . 1 :

 231 Figure 3.1: A hinge angle α(P, Q, K) (left) and four symmetrical hinge angles (right).

. 3 )

 3 The symmetries on hinge angles are important because it allows us to restrict rotations in the rst quadrant of the circle such that α ∈ [0, π 2 ]. Corollary 3.3. Each triplet (P, Q, K) corresponds to four symmetrical hinge angles α + πk 2 where k = 0, 1, 2, 3.

Figure 3 . 1 (

 31 Figure 3.1 (right) gives an example of Corollary 3.3.To distinguish the case of (K +1 2 , λ) from that of (λ, K + 12 ), we change the sign of K; we use α(P, Q, K) for the case of (K + 1 2 , λ), and α(P, Q, -K) for the case of (λ, K + 1 2 ). Because the symmetries allow us to restrict α to the range [0, π 2 ], as mentioned above, we may assume that K is positive.

Denition 3 . 4 .

 34 An angle θ is called Pythagorean if and only if both its cosine and sine belong to the set of rational numbers Q. We can deduce from Denition 3.4 that a Pythagorean angle θ is represented by an integer triplet (a, b, c) such that cos θ = a c , sin θ = b c . (3.4) In the following, θ indicates a Pythagorean angle. The lemma below for Pythagorean angles is well known. Lemma 3.5. Let (a, b, c) be an integer triplet generating a Pythagorean angle where

) from which we derive 2

 2 Qa+P b c ∈ Z. Because we know that c is odd according to Lemma 3.5, we obtain Qa+P b c ∈ Z. This means that 2(Qa+P b) c is even, while 2K + 1 is always odd, which leads to a contradiction.

Algorithm 3

 3 Discrete rotation algorithm by a lower bound rotation angle. Require: A set of points A, a lower bound rotation angle α (hinge angle) Ensure: A set of points A ′ for all ⃗ p ∈ A do α 1 = Function 1 (⃗ p, α); Move ⃗ p to (K, ⌊λ + 1 2 ⌋) or (⌊λ + 1 2 ⌋, K), depending on the sign of K and store the new point in A ′ ;

3 ⃗

 3 p1 and ⃗ q1 can have dierent coordinates since they are a pair of point in correspondence.

Figure 3 . 2 :

 32 Figure 3.2:The corners of H (⃗ q), namely, four corners of a pixel around ⃗ q.

Figure 3 . 4 :

 34 Figure 3.4: Running of the incremental algorithm.

  4.2, the function giving a pair of lower and upper bound rotation angles from a pair of points is realized in a constant time O(1). Moreover, as explained in Section 3, we can compare two hinge angles in a constant time O(1). Therefore, the computation of this algorithm for l points takes the time complexity of l × (O(1) + O(1)) = O(l) as a whole.

Figure 3 .

 3 Figure 3.4 gives an example of the incremental algorithm for two sets of three points.

Figure 3 . 5 :

 35 Figure 3.5: Randomly generated points and their corresponding rotated points.

Figure 3 . 6 :

 36 Figure 3.6: Result of our algorithm applied on the sets of points in Figure 3.5 .

  randomly generated the set F of 100 oating points in a 200 × 200 square. The rst set I is obtained by discretizing F . The second set I ′ is obtained by applying the DER with angle θ = 50 • to the rst set I (Figure 3.5). We assumed that the point correspondences across the two sets I and I ′ are known.

Figure 3 .

 3 Figure 3.6 shows the hinge angles obtained by our algorithm. The green and the red

Figure 3 . 2 d

 32 Figure 3.6 also indicates that most pairs of points give the admissible rotation angle range smaller than 1 • (for example, pair #39). Pair #34 and pair #66, however, give a range greater than 4 • . In fact, the maximum range for a pair of points is directly related to the distance between the center of rotation and the pair of points. We denote by d the distance between the rotation center and the pair of points. Then, the maximum dierence between the lower and the upper bound rotation angles is sin -1 √ 2 d . For the pair #34, d = √ 173; thus the maximum range is approximately 6.1 • . For the pair #39, d = √ 16505 and thus the maximum range is approximately 0.63 • . This is consistent with our experimental result.

Figure 3 .Figure 3 . 7 :

 337 Figure 3.7 shows the hinge angles obtained. Because we started with the set of oating points, we observe some errors in bounding admissible rotation angles. In fact, we see that some pairs of points, the pairs #64 and #73 for instance, do not contain θ.

Figure 3 .

 3 Figure 3.8 (left) explains how errors arise in rotating oating points. Let assume that the point ⃗ q 1 is obtained after the Euclidean rotation of the oating point ⃗ p 1 with angle

Figure 3 . 8 :

 38 Figure 3.8: Examples of errors for computing lower and upper bounds of rotation angles introduced by rotation of oating points.

Figure 3 . 9 :

 39 Figure 3.9: Result of our algorithm applied on real data.

√ 2 2dp.

 2 As mentioned in Section 3.6.1, thanks to Theorem 3.7, the addition of a Pythagorean angle does not change the nature of the upper (lower) bound rotation angle.

Figure 3 .

 3 Figure 3.9 (d) shows the results obtained by the modication above to the same data.

Figure 3 . 10 :

 310 Figure 3.10: Problem of corresponding points with real data.

Figure 4 . 1 : 1 2

 411 Figure 4.1: All hinge angles in the rst quadrant for the grid point ⃗ p = (2, 1) ⊤ , such that α 1 ≈ -13.68 • , α 2 ≈ 15.59 • , α 3 ≈ 21.31 • , α 4 ≈ 50.52 • . The hinge angle α 5 is obtained by symmetry such that α 5 = π 2 -α 1 ≈ 76.32 •

Figure 4 .

 4 Figure 4.1 illustrates all the hinge angles in the rst quadrant for the grid point (2, 1) ⊤ .

Figure 4 .

 4 2(a).

Figure 4 .

 4 2(b) shows that the intersection consists of three dierent sets of parallel lines, except for cases where the normal of the rotation plane is parallel to one of the axes, dening the coordinate system of the 3D space. As such exceptional cases provide only two dierent sets of parallel lines, which are identical with those of the 2D half-grid, we do not take into account those cases here. In other words, in such cases, 3D discrete rotations become identical with 2D discrete rotations. We call the three dierent sets of parallel lines a multi-grid,

Figure 4 . 2 :

 42 Figure 4.2: The 3D half-grid cut by a plane (a), and its multi-grid (b).

Figure 4 . 3 :

 43 Figure 4.3: Parallel lines of a set L A,⃗ p i and geometric interpretation of their parameters

Denition 4 . 5 .Figure 4 . 4 :

 4544 Figure 4.4:The ve dierent shapes of convexels, which are constructed as the intersections between a rotation plane and voxels[START_REF] Andres | New methods in oblique slice generation[END_REF].

Remark 4 . 6 .

 46 On any M A,⃗ p , the convexel containing ⃗ p always has a symmetric shape, and ⃗ p is at the center of the convexel. Moreover, the convexel that contains ⃗ p is necessarily of the type (b),(d) or (e) according to Figure 4.4.

2 and L A,⃗ p 3 .

 23 Let A = (a x , a y , a z ) ⊤ and ⃗ p = (p x , p y , p z ) ⊤ in the standard orthonormal base B and the plane P with normal vector A that goes through ⃗ p. Assuming lines in L A,⃗ p 1 come from the intersection between P and the planes of the 3D half-grid that are parallel to yz-plane, denoted by x = k where k ∈ Z, we obtain the directional vector ⃗ v 1 of lines in

Figure 4 . 5 :

 45 Figure 4.5: Two multi-grids generated by a rotation axis but with two dierent points ⃗ p 1 , ⃗ p 2 one of which is obtained by a translation of ±(⃗ p 1 -⃗ p 2 ) from the other.

  |ax|+|ay|+|az| gcd(ax,ay,az) or |ax|+|ay|+|az| gcd(ax,ay,az) + 1.

Figure 4 . 7 :

 47 Figure 4.7: the two impossible cases of intersection between a circle centered on a grid point and a rational multi-grid.

Figure 4 .

 4 7(a) or (b) presented as Lemma 4.10 in the following. Lemma 4.10. Let M A,⃗ p be a multi-grid where A ∈ Z 3 and ⃗ p ∈ Z 3 . Then, the locus of the rotation of ⃗ p on the rotation plan P does not go through any vertex of the convexels on M A,⃗ p .

Figure 4 . 8 .

 48 The coordinates P x and P y , of ⃗ p in B i , are obtained by applying P BB i to the coordinates of ⃗ p in B. Note that P x and P y are rationals. The coordinate Q y is obtained from one of (4.3), (4.6) and (4.7) depending on the values of i and k. All the values in (4.3), (4.6) and (4.7) are rationals, and thus Q y is also rational.

Figure 4 . 8 :

 48 Figure 4.8: A hinge angle α for a point ⃗ p in a rational multi-grid.

  rotation A. Therefore, we dene the distance function d(⃗ p) that is the Euclidean distance between A and ⃗ p. Then, the rotation of ⃗ p around A intersects at most 3⌊d(⃗ p)⌋ planes of

ure 5 .

 5 1 (c) forming a square pyramid. The minimum circular cone including the square pyramid has the solid angle 2 √

Figure 5 . 1 :

 51 Figure 5.1: Example of the approximation of a 3D vector using our method. (a) and (b) represent the projections of ⃗ r into ⃗ r 2 and ⃗ r 3 in 2D planes respectively. Each 2D convex cone in (a) and (b) is the admissible approximation for ⃗ r 2 and ⃗ r 3 regarding to ϵ. (c) represents the square pyramid constructed from the two 2D convex cones that contains the sets of admissible approximation for ⃗ r.

Firstly, we remind

  the denition of admissible rotation angles abbreviated by ARA and introduced in Section 3.3. The ARA for the two sets of n points A = (⃗ a 1 , ⃗ a 2 , . . . , ⃗ a n ) and B = ( ⃗ b 1 , ⃗ b 2 , . . . , ⃗ b n ) in correspondence is the set of angles dened by an upper and

Figure 5 . 2 :

 52 Figure 5.2: In a rational multi-grid M A,⃗ p a convexels can have four intersections with the locus of rotation of ⃗ p.

Figure 6 . 1 :

 61 Figure 6.1: Two representations of a Pythagorean triple (i 1 , i 2 , λ): Pythagorean angle θ and Pythagorean vector ⃗ v in the plane

Figure 6 . 2 :

 62 Figure 6.2: Construction of all Pythagorean triples of the form (i 1 , i 2 , i 1 + 1)

Figure 6 . 3 :

 63 Figure 6.3: A convex cone in the plane dened by the pair of vector ⃗ v 1 , ⃗ v 2

1 4 *

 14 + is equivalent to ]0, +∞[.

Figure 6 . 4 :

 64 Figure 6.4: Geometrical representations of a Pythagorean quadruple (i 1 , i 2 , i 3 , λ) by the Pythagorean vector (i 1 , i 2 , i 3 ) ⊤ and the pair of angles (θ, ψ).

Figure 6 . 5 :Figure 6 . 6 :

 6566 Figure 6.5: A 3D convex cone generated by three vectors ⃗ v 1 , ⃗ v 2 , ⃗ v 3 in the space

Figure 6 . 7 :

 67 Figure 6.7: (a) The plane P dened by OZ and ⃗ u, and (b) the vector ⃗ w that belongs to C ′ 2 .

. 1 )

 1 Applying these three matrices to the rst Pythagorean triple < 3, 4, 5 > will give the second level of the tree, then applying these three matrices to the three children of < 3, 4, 5 > will give nine children that are the third level of the tree. Iterating this operation will give us all possible prime Pythagorean triples, the three rst levels of the Pythagorean tree are given in example in Figure6.8. Note that in[START_REF] Nouvel | Rotations discrètes et automates cellulaires[END_REF] a dierent Pythagorean tree is presented, in his tree the author sorts the Pythagorean triple according to their angles.

Figure 6 . 8 :

 68 Figure 6.8: The three rst level of the Pythagorean tree.

. 4 ) 4 on

 44 that is the extension of the matrix that gives the father for a Pythagorean triple. An example of a Pythagorean quadruple that requires M 4 to be reached is illustrated on Figure6.9. The Pythagorean quadruple < 23, 24, 24, 41 > is reached using M Pythagorean quadruples < 10, 10, 23, 27 > or < 11, 12, 24, 29 >. Using M 4 on < 23, 24, 24, 41 > will give < 10, 10, 23, 27 > or < 11, 12, 24, 29 >, then we deduce that the fathers of < 23, 24, 24, 41 > are < 10, 10, 23, 27 > or < 11, 12, 24, 29 > and that we cannot reach it using matrix M 1 , M 2 , M 3 . Even if for most Pythagorean quadruples, M 4 gives the father of the Pythagorean quadruple as illustrated on Figure 6.10, there is an innity of Pythagorean quadruple that are their own grandfather. In other words, there is an innity of Pythagorean quadruples such that the application of M 4 gives another Pythagorean quadruple and

Figure 6 . 9 :

 69 Figure 6.9: The Pythagorean quadruple < 23, 24, 24, 41 > can only be reached using Matrix 6.4.

Figure 6 . 10 :

 610 Figure 6.10: In most cases, M 4 gives the father of a Pythagorean quadruple.

Figure 6 . 11 :

 611 Figure 6.11: There exists an innity of Pythagorean quadruples that have themselves as grandfather.

Figure 6 .

 6 Figure 6.12: (a) Distance between 1000 vectors and their approximated Pythagorean vector for a Pythagorean tree of depth 4 (a), of depth 5 (b), of depth 6 (c), of depth 7 (d) and of depth 8 (e)

Figure 6 . 13 :Figure 6 . 14 :

 613614 Figure 6.13: Distance between 1000 vectors and their approximated Pythagorean vector for a partial Pythagorean tree of depth 10 (a), of depth 20 (b), of depth 50 (c), of depth 100 (d), of depth 200 (e) and of depth 300 (f)

  Denition 2.3. A function F : E 2 → E 2 is bijective if for all elements e ∈ E 2 , there is exactly one element e ′ ∈ E 2 such that F -1 (e ′ ) = e.

	Another problem is arising during composition of discrete rotations. Let ⃗ p 1 be a point
	and α, θ be two angles. The point ⃗ p 2 is the result of the application of DER on ⃗ p 1 by
	α and ⃗ p 3 is the result of the application of DER on ⃗ p 2 by θ. We cannot ensure that
	the result of the application of DER on ⃗ p 1 by α + θ is ⃗ p 3 as illustrated on Figure 2.1(c).

  The set of points each of which has a semi-integer coordinate is called the half-grid and is denoted by H . Thus, H represents the set of points on the frontiers of all pixels whose centroids are points in Z 2 .

2

. We say that ⃗ p has a semi-integer coordinate if x + 1 2 ∈ Z or y + 1 2 ∈ Z. Denition 3.1. An angle α is called a hinge angle if at least one point in Z 2 exists such that its image by the Euclidean rotation with α belongs to H .

  .11) Algorithm 2 Function for nding the lower bound rotation angle for a digital image. Require: A set of points A, a Pythagorean angle θ

	Ensure: A hinge angle α
	α = Function 1 (rst point ⃗ p 1 of A, θ); for all ⃗ p ∈ A\{⃗ p 1 } do
	α temps = Function 1 (⃗ p, θ); if α < α temps then
	α = α temps ; end if end for return α;

Table 3 . 1 :

 31

  .9) Since all values in (4.7) are integers , we can deduce that for each ⃗ v, R(⃗ v) is an integer. Thus, if |a x | + |a y | + |a z | is even, we can conclude that there are |a x | + |a y | + |a z | + 1 dierent values for the arithmetical rest R(⃗ v).

	If

|a x | + |a y | + |a z | is odd then, (4.9) becomes

  If α 1 is greater than α 2 , cos α 2 -cos α 1 > 0. Thus we have, from (4.14),

	(P 2 1x + P 2 1y )(P 2x Q 2x + P 2y Q 2y ) > (P 2 2x + P 2 2y )(P 1x Q 1x + P 1y Q 1y ).	(4.16)
	For simplicity, let A 1 = (P 2 1x + P 2 1y )P 2x Q 2x , B 1 = (P 2 1x + P 2 1y )P 2y Q 2y , A 2 = (P 2 2x +
	P 2 2y )P 1x Q 1x and B 2 = (P 2 2x + P 2 2y )P 1y Q 1y . Note that A 2 1 , B 1 , A 2 2 , B 2 ∈ Q. Now (4.16) is
	rewritten as

  Algorithm 4 Rotation of a 3D image around an integer-axis whose direction is A by a Pythagorean angle θ.

	Require: An image I, a vector A, a Pythagorean angle θ Ensure: A rotated image I ′ 1: for all points ⃗ p in I do
	2:	Set T to be an empty list;
	3:	Compute the three generic equations (4.3), (4.6) and (4.7) for each M A,⃗ p i , i =
		1, 2, 3;
	4:	for all lines in M A,⃗ p i , i = 1, 2, 3 do
	5:	Compute all hinge angles corresponding to ⃗ p and the current line,
	6: 7:	and add α(p x , p y , p z , ±i, k) to the list T ; end for
	8:	Sort all the hinge angles corresponding to ⃗ p in T ;
	9:	Search in T the greatest hinge angle α which is smaller than θ;
	10:	Copy the voxel color from ⃗ p in I to the rotated point ⃗ q with α in I ′
	11: end for 12: return I

′ ; the half-grid and denes at most 6⌊d(⃗ p)⌋ dierent hinge angles. Because d(⃗ p) ≤ √ 3m, the upper bound of the number of hinge angles for any point in the image is 6 √ 3m. Accordingly, we can conclude that the upper bound of the number of hinge angles for a given image of size m × m × m is 6 √ 3m 4 ; thus O(m 4 ).

  there is a vector ⃗ u ∈ U that belongs to C.In other words, U is dense if for any set of n vectors ⃗ v 1 , ⃗ v 2 , . . . , ⃗ v n ∈ V , there exists a vector ⃗ u ∈ U and s 1 , s 2 , . . . , s n

  3, if we multiply all integers of a Pythagorean triple by k ∈ Z, the obtained triple is still a Pythagorean triple. This enables us to design from two Pythagorean triples a Pythagorean quadruple as follows. Let a, b ∈ Z * such that aλ ⃗ u = bj 1 . Then we obtain the Pythagorean quadruple (ai 1 , ai 2 , bj 2 , bλ ⃗ w ). By construction, it is obvious that the vector ⃗ w = (ai 1 , ai 2 , bj 2 ) ⊤ , associated with the Pythagorean quadruple (ai 1 , ai 2 , bj 2 , bλ ⃗ w ), belongs to the 3D convex cone C 3 .

  3 on this quadruple will only give two dierent sons, that will lead to repetitions. Pythagorean quadruple of the form < i 1 , i 2 , i 2 , λ > like < 1, 2, 2, 3 > and < 4, 4, 7, 9 > are also problematic since they have many repetitions in their sons. The solution may be to remove all repeti-

tions after generating the sons of a Pythagorean quadruple. However, this does not solve the problem since some Pythagorean quadruples can be generated by two distinct Pythagorean quadruples. For example, the two Pythagorean quadruples < 2, 3, 6, 7 > and < 1, 4, 8, 9 > generate the Pythagorean quadruple < 3, 4, 12, 13 > when we apply M 3 . One other problem of the naive method is that it does not generate all possible

  The computation time depends on the depth d of the Pythagorean three. When d = 8, it takes almost one hour to nd the 1000 closest Pythagorean vectors or around three seconds by vector to approximate, but when d = 4, it takes around half a second to nd the 1000 closest Pythagorean vectors or 0.5ms for each Pythagorean vector. The evolution of computation time between two levels of the tree have a factor 9. As we can see, the average distance between the original vector and its approximation decreases between depth 4 and depth 8 from 0.034 to 0.0077. However, the maximum distance is still high in both cases, around 0.15 when d = 4 and 0.073 when d = 8. Analyzing the results shows that all the vectors with a bad approximation have two small coordinates and one a lot greater. The ve vectors with the worst approximation for d = 8 are < 1628, 1706, 30389 >, < 179, 214, 3539 >, <

	358, 437, 7611 >, < 547, 668, 11872 >, < 1057, 1106, 21398 > for error between 0.07 and
	0.075.

Accordingly, DER is not bijective.

A discrete rotation is a rotation designed for the discrete space. It transforms a set of grid points into another set of grid points.

Acknowledgements

quadruples that are the closest one to the given vector and compute only their sons.

Roughly speaking, the rst method is a force brute algorithm. The second one cuts the Pythagorean tree. We know that the rst one will give us a result since Pythagorean vectors are dense and it tries all possible Pythagorean quadruple until it nds an ϵ-close Pythagorean vector. For the second one, we still do not have any proof of its convergence.

The rst approximation method starts from the Pythagorean quadruple < 0, 0, 1, 1 > and its two sons < 0, 3, 4, 5 > and < 1, 2, 2, 3 >. For this three Pythagorean quadruples, we compute the distance between ⃗ v and these three Pythagorean vectors. If none of them is ϵ-close, we compute the sons of < 0, 3, 4, 5 > and < 1, 2, 2, 3 > as explained in Section 6.4.1 and for each sons, we compute the distance between them and ⃗ v. If none of them is ϵ-close, we iterate the operation until we nd a matching vector. Since Pythagorean vectors are dense, we know that this algorithm will end in a nite time, however, we cannot give any upper bound of the number of level to compute, then we cannot give any upper bound of the computation time. The problem of this method is the number of Pythagorean quadruple to compute and to compare with ⃗ v. As explained in Section 6.4.1, most of the Pythagorean quadruples have nine sons, thus the number of Pythagorean quadruples for a Pythagorean tree with n levels can be approximated by n ∑ i=1 9 i considering that < 0, 3, 4, 5 > and < 1, 2, 2, 3 > are the level 0 of the tree. Thus reaching the 10th level of the tree requires about 3 billions comparisons.

The second approximation method cuts the Pythagorean tree. The rst step is the same as the rst method. After computing the rst level of the Pythagorean tree, we compute the distance between ⃗ v and all computed sons and check if one of the son is ϵ-close to ⃗ v.

If not, we sort the sons regarding their distance to ⃗ v and search a convex cone formed by three of the sons. Because there are many convex cones, we consider the one that minimize the sum of the distances between the three sons and ⃗ v. To decide if the convex cone formed by three vectors contains ⃗ v, we compute the three planes formed by two of these three vectors and the origin and check if ⃗ v and the vector that does not belong to the plane belongs to the same half-space. The next step is to compute the sons of the three vectors that form the convex cone and compute the distance between ⃗ v and them. If none of them is ϵ-close to ⊑, we compute a new convex cone and iterate the last operation until we found and ϵ-close Pythagorean vector. For each iteration we compute a maximum of 36 Pythagorean quadruples (12 for each vector of the convex cone), then the number of Pythagorean quadruples computed may be seriously reduced in comparison with the rst method. However, unlike the rst method, we cannot ensure that this method will nish since we do not have any proof on the repartition of the Pythagorean quadruple depending on the father.

The next section presents some experiments of the two methods presented here.

The second method leads to an interesting observation on the repartition on Pythagorean quadruples. Even when we know the father, we cannot predict the direction of the sons.

Then the study of the repartition of the Pythagorean quadruples on the unit sphere will be dicult. Therefore, from this observation, we can design a third method. The third method is almost the same than the second one; the main dierence is that instead of keeping the three Pythagorean vectors that form a convex cone around the vector to approximate, we select three vectors randomly among the sons computed from the three previous vectors.

For the third method, we set d = 10, 20, 50, 100, 200, 300 and use the same 1000 vectors to approximate than the rst method. These experiments are illustrated in Figure 6.14(af) . The computation times to approximate each vector are almost the same than the computation times obtained with the second method. However, the results are a lot better when d gets high. The obtained result when d = 300 is really close to the result obtained by the rst method when d = 6 since the maximum distance between the vector to approximate and its approximated Pythagorean vector is lower than 0.1 and the average distance is 0.015 for the third method and 0.014 for the rst method.

The third method gives a good result; however, this is a theoretical method since the integers which form the Pythagorean quadruples when d = 300 are huge (≈ 10 40 ) and need a special library to be coded in all languages.

From our experiments, we can conclude that the only practical method is the rst one.

Pythagorean n-tuples

In this section, we rst present Pythagorean n-tuples and their representation in n -1 dimensions. Then, we extend Theorem 6.13 to any dimensions. The following denition is the extension of Denitions 6.1 and 6.11. Denition 6.14. A Pythagorean n-tuple is an integer n-tuple

In order to simplify the reading, hereafter, n denotes the dimension of their representation, instead of the numbers of integers of the n-tuples. In other words, we will use notations nD Pythagorean vectors and Pythagorean (n + 1)-tuples.

Similarly to the cases of n = 2, 3 in Sections 6.2 and 6.3, we dene

which is the set of Pythagorean (n + 1)-tuples.

We represent each element of P n by an nD Pythagorean vector dened as follows.

Appendix 96 we have the proof that these representations are unique, but by carefully reading the proof, we see that for dimension higher than four, the uniqueness disappear. The work presented in [START_REF] Rodriguez | Generalized perpendicular bisector and circumcenter[END_REF] that must be extended to nd admissible rotation axis in 3D must also be extended in nD to nd the set of ⌊ n 2 ⌋ admissible plane of rotation or at least one set if the decomposition is not unique. We can conclude that the description of discrete rotation methods in nD will requires lot of work.

Chapter 8 Appendix

This last chapter is a recall of all denitions used in this manuscript. Here the denitions are not formal, just intuitive. It may help the reader to have a better understanding of the notions or nd a denition that he forgot. Convex hull: Imagine a land containing a nite number of trees. You attach a rope to the most eastern tree, surround all the trees with the rope and attach again the rope to the same tree. The rope is the convex hull of the trees.

Convexel: Take a cube, slice it one time. The new surface obtained is a convexel. This Discrete geometry: It is a subeld of mathematics. Simply put, Discrete geometry is the part of mathematics considering only nite set. Instead of considering an innity of numbers, we just select some of these numbers such as all the integers smaller than one billion and try to nd the geometrical properties of this set.

Discrete grid: Can be seen as a computer screen. On a continuous plane, there is always an innity of points between two points. On a discrete grid, there is a nite number of points between two points. If you consider two neighbor pixels on a computer screen, there is no pixel between them. This dierence implies one of the main issue in computer science.

Discrete point: A point in any dimension with integer coordinates. On a computer screen, it can be considered as a pixel.

Floating numbers: In computer science, they are an approximation of the real numbers. In math, there are dierent sets of numbers, the two most used sets are natural numbers and real numbers. In computer science, it is possible to represent natural numbers but, since most of the real number have an innity of digit after the point, they cannot be exactly represented in the memory of a computer. Thus, we use the oating numbers instead. The dierence between a real number and its oating approximation necessarily implies errors during computations.

Grandfather (tree): In computer science, a tree look like a family-tree. The notion of grandfather is exactly the same.

Half-grid: In any dimension, the half-grids are the union of all hyper-plans of equations

2 where x represents one of the dimensions and i ∈ Z. The half-grid, is the border between the hypercube of dimension n which tiles the space. On a paper, the half-grid can be compared to the line on a graph paper.

Multi-grid: In any dimension, the multi-grids are the intersection between the half-grid and a given plane.

Pixel: Pixels can be seen in the discrete space as points in the continuous space. On a chess board, the pixel will be a black or white square, since a piece belongs to only one case, no matter where it is. On a computer screen, pixels are the smallest points on the screen.

Points in correspondence: Two points from two dierent digital images (in any dimension) are said to be in correspondence if they represent the same things in both digital images. For two digital images of the Eiel tower, the two pixels which represent the top of the tower in both images are said to be in correspondence.

Prime:

In other words, a set of element is prime, if there is no element that divide each element of the set.

Pythagorean angles: An angle is Pythagorean if its sine and its cosine are rational. Each Pythagorean angle is associated to a Pythagorean triple. They are useful to compute rotations using only integers.

Pythagorean n-tuples: Read rst Pythagorean triples. Instead of a right triangle in the plane, we imagine a right triangle in the (n-1)D space. We then add n-3 integers to the triple to obtain a n-tuples of integers i

Pythagorean quadruples: Read rst Pythagorean triples. Instead of a right triangle in the plane, we imagine a right triangle in the 3D space. We then add one integer to the triple to obtain a quadruple of integers i 1 , i 2 , i 3 , λ such that i

Pythagorean triples: Firstly, we have to remember our old mathematic classes with Pythagoras of Samos [START_REF] Pythagoras | How to found a sect based on three numbers, 520BC[END_REF]. Hypotenuse and so on ... So if the three sides of a right triangle (Pythagorean triangle) have integer norms, then the three norms form a Pythagorean triple. In other words, each triple of integers i 1 , i 2 , λ such that i 2 1 + i 2 2 = λ 2 is a Pythagorean triple.

Pythagorean vectors: Consider a Pythagorean triple, the associated Pythagorean vector is the vector with the coordinates (i 1 , i 2 ) ⊤ and with the norm λ.

Rational computation: Equivalent to integer computation since a rational number can be seen as a pair of integers.

Tree: It is a structure mostly used in computer science in order to represents data. Like a real tree, it has a root and many branches denoted by nodes. A node belongs to only one branch but one branch can have many nodes.

Topology: It's the study of the spatial properties of objects under continuous deformations. Roughly speaking, if we consider an object done with modeling clay, we can scratch Appendix 100 it, extend it without modifying its topology. However if we cut it, create a hole, etc the topology is modied. A donut and a coee cup are topologically equivalent since they both have only one hole. According to [START_REF] Chaussard | Topologie ???!!![END_REF], topology is a useless branch of mathematics and computer science.

Voxel: Read rst Pixel. Remember the Lego toys. Consider an object made only by using smallest colored cubes. Each cube in the object can be seen as a voxel. They are the smallest and indivisible parts of the object (except if you introduce a three year old child, but we do not consider bugs in this manuscript).