Yohan Thibault

Akihiro Sugimoto

Yukiko Kenmochi

3D discrete rotations using hinge angles

Keywords: 3D discrete rotation, half-grid, multi-grid, hinge angle, integer computation

come   L'archive ouverte pluridisciplinaire

Introduction

Rotations in the 3D space are required for computer imagery in image processing [START_REF] Toffoli | Three-dimensional rotations by three shears[END_REF], computer vision [START_REF] Hartley | Global optimization through rotation space search[END_REF][START_REF] Schmidt | Using quaternions for parametrizing 3-D rotations in unconstrained nonlinear optimization[END_REF] and computer graphics [START_REF] Kaufman | 3D volume rotation using shear transformations[END_REF]. A rotation in the 3D Euclidean space can be in general represented in two typical ways. One is to represent a rotation as the combination of three rotations around the three axes of the coordinate system in concern [START_REF] Voss | Discrete Images, Objects and Functions in n[END_REF]. The other is to represent a rotation by a rotation axis together with an angle around the axis [START_REF] Hartley | Global optimization through rotation space search[END_REF][START_REF] Schmidt | Using quaternions for parametrizing 3-D rotations in unconstrained nonlinear optimization[END_REF]. Even if the representations of a rotation are different, 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 is the case of 2D rotations, computing a 3D rotation once in the discrete space brings displacement from that in the continuous space; computing 3D rotations many times causes difficulty 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 popular for 3D rotation estimation from given images, which is one of important problems in computer vision [START_REF] Hartley | Global optimization through rotation space search[END_REF][START_REF] Schmidt | Using quaternions for parametrizing 3-D rotations in unconstrained nonlinear optimization[END_REF]. This paper presents a study of the rotation in the 3D discrete space. Since we admit only integer computation, 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 correspond to the discontinuity caused by discretization 1 of the rotation in the continuous plane [START_REF] Amir | Twodimensional pattern matching with rotations[END_REF][START_REF] Nouvel | Rotations discr tes et automates cellulaires[END_REF][START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF]. Intuitively, hinge angles determine a transit 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 paper, we first define 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 a "multi-grid", that is given by the intersection between a plane normal to the rotation axis and the half-grid, i.e., the boundary between adjacent voxels. The rational multigrids, which are a subset of multi-grids, allow computations of hinge angles using only rational numbers. Using rational multi-grids, we show that, as in 2D, it is possible to use only integers 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 which transforms the first digital image into the second one. Note that in a discrete space, there exists a set of rotations that give the same rotated image from a given digital image. In this paper, we fix 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 the set of admissible rotation angles and its upper and lower bounds are represented by hinge angles. This method is the extension of the method proposed in [START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF] for the 2D cases into 3D cases.

Differently 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] Toffoli | Three-dimensional rotations by three shears[END_REF][START_REF] Kaufman | 3D volume rotation using shear transformations[END_REF]. In particular, to our best knowledge, this is the first work on 3D discrete rotations using integer computations without digitization errors. This paper is an extended version of our previous paper [START_REF] Thibault | Hinge angles for 3D discrete rotations[END_REF]. Sections 3 and 4 about multi-grids and hinge angles are re-organized as more properties of multi-grids are added in this paper (see section 3.4 for the properties). Besides, we newly add an incremental rotation algorithm using 3D hinge angles in Section 5.4 as well as a study for finding a 3D discrete rotation from a pair of 3D digital images in Section 6.

Hinge angles

Hinge angles for 2D rotations are defined 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.

Definition 1. The 2D half-grid is the set of lines in the plane, each of which is represented by one of x = i + 1 2 and y = i + 1 2 where i ∈ . In other words, the 2D half-grid represents the border between two adjacent pixels. From the definition of the 2D half-grid, we define the hinge angles in the plane.

Definition 2. An angle α is called a hinge angle if at least one point in 2 exists such that its image by the Euclidean rotation with α around the origin is on the 2D half-grid.

Figure 1 illustrates all the hinge angles in the first quadrant for the grid point (2, 1) ⊤ . Note that hinge angles in the other quadrants are obtained by symmetry with respect to the x-axis and/or y-axis from those in Figure 1.

To extend the definition of hinge angles into the 3D case, we first define the half-grid in the 3D space. Similarly to the 2D half-grid, the 3D half-grid defines the limit between two adjacent voxels in the 3D discrete space. Definition 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 ∈ . Introducing the definition of the 3D half-grid allows us to define 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. Definition 4. An angle α is called a hinge angle if at least one point in 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 discretised point of the rotation result of p with angle α + ǫ becomes different from that with angle αǫ for any ǫ > 0.

Differently 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 a variety of transitions of a grid point across voxels exist depending 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 2, hinge angles are strongly related to the halfgrid both in 2D and 3D. However, in the 3D case, rotations of a point are always in its rotation plane that is 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 definition 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.

Definition of multi-grids

When a rotation plane in 3D is given, the intersection between the plane and the half-grid in the 3D space is obtained as illustrated in Figure 2(a). Figure 2(b) shows that, the intersection consists of three different sets of parallel lines, except for cases where the normal of the rotation plane is parallel to one of the axes defining the coordinate system of the 3D space. As such exceptional cases provide only two different sets of parallel lines, which are identical with those of the 2D half-grid, we here do not take into account those cases. In other words, in such cases, 3D discrete rotations become identical with 2D discrete rotations. We call the three different sets of parallel lines a multi-grid, which is used for characterizing hinge angles for 3D discrete rotations instead of the 2D half-grid for 2D discrete rotations.

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) ⊤ ∈ Ê 2 |a i x + b i y + c i + k = 0, a 2 i + b 2 i = 0, k ∈ , a i , b i , c i ∈ Ê} (1)
where i = 1, 2, 3. The integer parameter k denotes the index number of each parallel line. Figure 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 ai , y) ⊤ is on the k-th next line, providing that a i = 0. Now we can give a formal definition of a multi-grid. Definition 5. Let L i for i = 1, 2, 3 be each set of parallel lines induced from a given rotation plane and the 3D half-grid. Then the multi-grid M is the union of L i : M = ∪ 3 i=1 L i . Hereafter, we denote by L A,p i for i = 1, 2, 3 the three sets of parallel lines defined by a rotation plane with a normal vector A = (a x , a y , a z ) ⊤ going through point p = (p x , p y , p z ) ⊤ . Using the same idea, we denote by M A,p the multi-grid defined as the union of M A,p = ∪ 3 i=1 L A,p i . The multi-grid in the rotation plane forms various closed convex polygons surrounded by lines, which we call convexels. Depending on the rotation plane, we have a variety of shapes of convexels in this paper. The convexel in a multigrid is the counterpart of the squared pixel defined 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. Note that another study on convexels is reported in 4: The five different 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]. [START_REF] Reveillès | The geometry of the intersection of voxel spaces[END_REF], to clarify the intersection of overlapping voxels for medical image fusion; the required properties there are different for those in this paper. Remark 6. 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) in Figure 4.

We remark that when the normal of the rotation plane is parallel with one of the axes defining the coordinate system of the 3D space, the convexel coincides with the pixel.

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 2 and L A,p 3 . 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 ∈ , we obtain the directional vector v 1 of lines in 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 BB1 = 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 , (2)
where ψ

= 1 √ a 2 x +a 2 y +a 2 z
. We remark that A is transformed to (0, 0) ⊤ by P BB1 :

the rotation center in P thus becomes the origin of B 1 . We remark that if A is collinear with one of the axes defining the coordinate system of the 3D space, P BBi degenerates: the rank of M BBi becomes 1. In such cases, 3D rotations become identical with 2D rotations. However, in Section 3.1, we explained that we do not consider the cases where the vector A is collinear with an axis defining 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 BB1 to the plane x = k for k ∈ 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. (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 BBi for i = 2, 3 such that

P BB2 = -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)
P BB3 = 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 -ψ , (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, (6)
ψ(a 2 x + a 2 y )y + k -A.pψ 2 a z = 0, (7)
where k ∈ . We note that (3), ( 6) and ( 7) correspond to the equation in (1) for i = 1, 2, 3 respectively. All a i in (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 in (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 different point p are similar; more precisely, they are obtained simply by translations as illustrated in Figure 5.

Rational multi-grids

In Section 3.2, we obtained (3), ( 6) and ( 7) for the lines in the multi-grid M A,p . In general, parameters of these equations belong to Ê. In order to use only integers during computation, we need these parameters to belong to É. If all elements of P BBi belong to É, then all the parameters in ( 3),( 6) and ( 7)

become rational as well. In order to obtain rational parameters in P BBi we set A to be a Pythagorean vector:

a vector v = (i 1 , i 2 . . . , i n ) ⊤ , i 1 , i 2 , . . . , i n ∈ , is a Pythagorean vector if v = λ where λ ∈ .
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 defined from a Pythagorean axis a rational multi-grid in the following.

The rational multi-grids presented in this section are in a subset of the multigrid. In addition to allowing us to use only integers during the computations, rational multi-grids also offer some useful properties which are valid only for rational multi-grids. Note that studies on more general multi-grids generated by planes whose normal vectors are integer-axes can be found in [START_REF] Andres | New methods in oblique slice generation[END_REF][START_REF] Reveillès | The geometry of the intersection of voxel spaces[END_REF][START_REF] Feschet | A generic approach for n-dimensional digital lines[END_REF].

Properties of rational multi-grids

In the case of multi-grids that are not rationals, there is an infinity of different convexels. However, for rational multi-grids, the number of convexels is finite and depending only on the coordinates of A that is the normal vector of rotation planes. 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 define the arithmetical rest of a voxel v = (v x , v y , v z ) ⊤ for a rotation plane P such that a x x + a y y + a z z + A.p = 0 as follows:

R(v ) = a x v x + a y v y + a z v z + A.p. (8)
Theorem 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 different 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. We know from [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 a plane P if and only if:

- |a x | + |a y | + |a z | 2 ≤ R(v ) ≤ |a x | + |a y | + |a z | 2 . ( 9
)
Since all values in (8) are integers, we can deduce that for any v , R(v ) is an integer. Thus, if |a x | + |a y | + |a z | is even, we can conclude that there Figure 6: A quadruple of grid points {q 1 = q , q 2 = (qx + Ax, qy + Ay, qz) ⊤ , q 3 = (qx + Ax, qy, qz -Az) ⊤ , q 4 = (qx + 2Ax, qy + Ay, qz -Az) ⊤ } forms two types of triangles, α and β tiled in a multigrid.

are |a x | + |a y | + |a z | + 1 different values for the arithmetical rest R(v ). If |a x | + |a y | + |a z | is odd then, (9) becomes - |a x | + |a y | + |a z | 2 < R(v ) < |a x | + |a y | + |a z | 2 (10)
and the number of different 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 different convexels will be either |ax|+|ay|+|az| gcd(ax,ay,az) or |ax|+|ay|+|az| gcd(ax,ay,az) + 1. 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 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 exist 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 +2A 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 8, illustrated in Figure 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 vertices {q 1 , q 2 , q 3 } and {q 2 , q 3 , q 4 } contain the same set of convexels and are mirror images. If A x , A y , A z are chosen to be prime, they describe the smallest triangle that tile 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 [START_REF] Nouvel | Rotations discr tes et automates cellulaires[END_REF][START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF]. In order to ensure the uniqueness of the representation, some properties on the multi-grid are required. For 3D hinge angles similar properties on multi-grid are required. We first show them on the 3D half-grid, and then represent 3D hinge angles using five integers and explain how to decode them to obtain their hinge angles.

In Section 2, we defined the 3D hinge angles in the framework of the 3D halfgrid. As manipulation of 3D hinge angles is more convenient in the framework of multi-grids and rational multi-grids, we propose an alternative definition of 3D hinge angles in the following. Proposition 9. Let p be a grid point and p ′ be the result of the rotation of p by an angle α around an integer-axis. α is a hinge angle if and only if p ′ is on the multi-grid.

Hinge angles and rational multi-grids

In 2D, there exists 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 axis of rotation 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 Figure 7(a) or (b), presented as Lemma 10 in the following. Lemma 10. Let M A,p be a multi-grid where A ∈ 3 and p ∈ 3 . Then, the locus of the rotation of p on the rotation plane P does not go through any vertex of the convexels on M A,p .

Proof. 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 ∈ , 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 Figure 7: the two impossible cases of intersection between a circle centered on a grid point and a rational multi-grid.

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, (11)
k x + 1 2 2 + k y + 1 2 2 + p ′2 z -(p 2 x + p 2 y + p 2 z ) = 0. (12)
From [START_REF] Thibault | Hinge angles for 3D discrete rotations[END_REF] we see that p ′ z must be a rational number, so that there exists a pair of integers λ 1 , λ 2 such that p ′ z = λ1 λ2 and gcd(λ 1 , λ 2 ) = 1. From [START_REF] Andres | New methods in oblique slice generation[END_REF] 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 , ( 13
)
where k, λ 1 , λ 2 ∈ . Since gcd(λ 1 , λ 2 ) = 1 and 2λ 2 1 = (2k + 1)λ 2 2 , λ 2 must be even. Setting λ 2 = 2λ ′ 2 where λ ′ 2 ∈ , 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, between which the transition of p during its rotation is done, 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 10 shows that cases (a) and (b) illustrated in Figure 7 cannot happen in a rational multi-grid. Note that the above proof need only the assumption that A has integer coordinates.

In 2D there also exists another property on hinge angles ensuring that the locus of rotation of a grid point cannot intersects twice 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. The more details will be discussed in Section 6.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). The first 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 defined. 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 five integers, we can obtain the coordinates (P x , P y ) ⊤ and (Q x , Q y ) ⊤ corresponding to p and q in the basis B i , as represented in Figure 8. Note that the projection matrix P BBi from B to B i is given by ( 2), ( 4) and ( 5). The coordinates P x and P y , of p in B i , are obtained by applying P BBi to the coordinates of p in B. Note that P x and P y are rational. The coordinate Q y is obtained from one of (3), ( 6) and [START_REF] Amir | Twodimensional pattern matching with rotations[END_REF] depending on the values of i and k. All the values in (3), ( 6) and ( 7) are rationals, and thus Q y is also rational. 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 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 define two points. To discriminate two different 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 defines no hinge angle since there is no transition of convexels. 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 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 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. Then α 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 11 is rephrased as: two different 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 rotation of the point. This indicates that there is no integer that decides whether the intersected line is parallel to the x-axis or the y-axis. To differentiate the cases, the last integer of the triple is set to be positive if the line is parallel to the x-axis; otherwise, it is set to be 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 [START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF]. 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 different hinge angles for an image is O(m 3 ) where m is the largest side of the image [START_REF] Amir | Twodimensional pattern matching with rotations[END_REF][START_REF] Nouvel | Rotations discr tes et automates cellulaires[END_REF]; we will use a similar method for obtaining the upper bound of hinge angles for a 3D image. Our algorithm, to be efficient, needs to compare a pair of hinge angles, which can be made 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 indeed trivial 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 5 we assume that A is a given 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 (Q x , Q y ) ⊤ in B i . P x , P y , Q x and Q y can be obtained from (p x , p y , p z , i, k) as explained in Section 4.2. The following equations are then derived from Figure 8:

cos α = P x Q x + P y Q y P 2 x + P 2 y , (14)
sin α = P x Q y -P y Q x P 2 x + P 2 y . ( 15
)
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 sines or cosines of two hinge angles using only integer computations as follows.

Proposition 12. Let α 1 and α 2 be two hinge angles defined for A. Then it is possible to decide if α 1 > α 2 using only integer computations.

Proof. 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 or cosines which are given in [START_REF] Feschet | A generic approach for n-dimensional digital lines[END_REF] and [START_REF] Kenmochi | On combinatorial properties of discrete planar surfaces[END_REF]. First we verify the signs of both of the sines and cosines of α 1 and α 2 . If the sines and cosines of both angles have different signs, then we can conclude whether α 1 > α 2 without additional 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.

If α 1 is greater than α 2 , cos α 2 -cos α 1 > 0. Thus we have, from [START_REF] Feschet | A generic approach for n-dimensional digital lines[END_REF],

(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 ). ( 16
)
For simplicity, let A 1 = (P

A 1 + B 1 > A 2 + B 2 . ( 17
)
Squaring both sides of (17) since they are not negative, and moving rational values to the left-hand side and the rests 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 . (18)
If the left-hand side and the right-hand side of (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 (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 (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 . ( 19
)
If the sign of the left-hand side of ( 19) 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 . ( 20
)
We note that we can easily verify whether (20) is satisfied with integer computation alone. If (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 (18) are negative. Thanks to Theorem 11, 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 different quintuplets.

Note that, if the comparison of two hinge angles is done using floating num-bers, only one comparison is required. However, if we want to keep integer computations, then in the worst case we have to check [START_REF] Kronrod | Optimal ordering algorithm without operational field[END_REF], [START_REF] Rodriguez | Generalized perpendicular bisector and circumcenter[END_REF] and [START_REF] Niven | An Introduction to the Theory of Numbers[END_REF].

From Proposition 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 be also done using only integers during its computation and in constant time. An angle is a Pythagorean angle if its sine and cosine belong to É [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 13 is similar to that of Proposition 12 and the one for comparison between a hinge angle and a Pythagorean angle in the plane presented in [START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF]. Proposition 13. Let θ be a Pythagorean angle and α be a hinge angle defined for A. Then it is possible to decide if α > θ in constant time with integer computations.

Proof. Let α = α(p x , p y , p z , i, k) be a hinge angle and θ be a Pythagorean angle associated to the Pythagorean triple (i 1 , i 2 , λ) such that cos θ = i1 λ and sin θ = i2 λ where i 1 , i 2 , λ ∈ . If α is greater than θ, namely cos θ -cos α > 0, then we have from ( 14)

i 1 (P 2 x + P 2 y ) > λ(P x Q x + P y Q y ). (21)
Using similar methods to those used for the proof of Proposition 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 . ( 22
)
We note that we can easily verify whether [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] is satisfied with integer computation alone. If [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] 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 ) [START_REF] Nouvel | Rotations discr tes et automates cellulaires[END_REF]. This is obtained by computing the bound for the furthest point from the origin in the image which is the rotation center 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 rotation A. Therefore, we define 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 the half-grid and defines at most 6⌊d(p)⌋ different hinge angles. Set T to be an empty list;

3:

Compute the three generic equations ( 3), ( 6) and ( 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, 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 ′ ;

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 [START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF]. 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 field 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 difference ǫ > 0 by a Pythagorean angle. The output of our algorithm is a rotated digital image I ′ .

The rotation algorithm is described in Algorithm 1. For each point p = (p x , p y , p z ) ⊤ in the image I, the algorithm computes the corresponding multigrid M A,p and search 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 12. The algorithm searches in the sorted list L the greatest hinge angle α which is smaller than θ using Proposition 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 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 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], a 2D incremental rotation is presented that allows obtaining all possible configurations 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 configurations of 3D discrete rotations even though 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 first 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 second step merges the lists of all points into the final list and then sorts all the hinge angles in the list. The final 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 final list T this operation should be done using the merge sort algorithm [START_REF] Kronrod | Optimal ordering algorithm without operational field[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 ) 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 configurations of rotations of I around the rotation axis have been reached.

The time complexity for the first step is O(m 4 log m) since there is m 3 points 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 as the first step. Note that if the merging is realized without using the method presented in [START_REF] Kronrod | Optimal ordering algorithm without operational field[END_REF], the time complexity remains the same, but practically it will increase the computation time. The last step browses the list T that contain 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).

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

In this section, we present a method to find a discrete rotation between a pair of 3D digital images that represents the same object from two different points of view. Suppose that we have a set of n grid points A = (a 1 , a 2 , . . . , a n ) in the first 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, finding 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 [START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF]. In Euclidean space such a rotation is unique. However, this is not the case in the discrete space where either of two slightly different angles α 1 , α 2 may transform A into B. Moreover, for such a pair (α 1 , α 2 ), every angle α 3 such that α 1 ≤ α 3 ≤ α 2 also rotates A into B. Thus, we define 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 find not only a set of admissible rotation angles but also an admissible rotation axis. However, the later problem is not treated in this paper, and is set to be our future work.

In this section, we identify a set of admissible rotation angles but only one admissible rotation axis. Indeed, if we consider two rotation axis, the two multigrids associated to these axes and a point will provide two distinct sets of hinge angles. In order to avoid such a situation, we set a unique rotation axis.

Approximation 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 first step is to obtain an axis of rotation that is consistent with all pairs of points in correspondence. The second step is to approximate the axis obtained through the first step by a Pythagorean axis in order to obtain rational multi-grids.

First step

Given a pair of Euclidean points in correspondence, it is known that any rotation axis such that the first 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 gives 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 planes may not meet at a 3D line. This is because grid points approximate Euclidean points.

Therefore, in order to obtain an appropriate rotation axis, we consider the rotation axis that is the average of the three lines, namely the intersections of all pairs of the three bisection planes. 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, each of which belongs to each voxel of the pair, defines a bisection plane. Therefore there is an infinity of admissible bisection planes 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 to define an admissible rotation axis, we have to compute the intersections between all the obtained bounds of admissible bisection planes. In 2D, a study on the intersections between all the bounds of bisections for n pairs of pixels have been proposed [START_REF] Rodriguez | Generalized perpendicular bisector and circumcenter[END_REF], and its extension to the 3D discrete space is still in a future work. Therefore, as an alternative, 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. We compute for each pair the bisection plane, and obtain ⌊ n+1 2 ⌋ axes from the n obtained bisection planes. The rotation axis is the average of the ⌊ n+1 2 ⌋ obtained axes. We expect that this method gives a rotation axis that belongs to the intersections between all the bounds of admissible bisection planes. The verification will be one of our future work; even though no formal proof will be found, empirical evidence will be given.

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 multi-grid. Based on the method presented in [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 definitions 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 θ = i1 λ . 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 to approximate 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 exists 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 define two 2D vectors r 2 , r 3 such that with an angle α i satisfying cos α i = xi |r i| . Using the algorithm presented in [START_REF] Anglin | Using pythagorean triangles to approximate angles[END_REF] 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 ∈ * . 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 .

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 (a) (b) (c)
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 twice Anglin's method. An example of the method that approximates a 3D vector is given in Figure 9. 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 9 (a) and (b), respectively. We generate r ′ from r ′ 2 and r ′ 3 , which is the approximation of r , such that it belongs to the 3D convex cone represented in Figure 9 (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 . 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 3D vector in O(1) operations. However, we cannot give any bounds on the size of integers that belong to the final 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 a method to obtain from n pairs of grid points the admissible rotation angles where n ≥ 1 in time complexity O(n) [START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF]. This method is designed to work in the 2D halfgrid. With minor modifications, this method can be also applied on rational multi-grids and keep the same complexity.

Firstly, we remind the definition of admissible rotation angles abbreviated by ARA introduced in [START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF]. 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 defined by upper and 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 a pair of points (a i , b i ).

For a given pair of points (a i , b i ) in correspondence, we first compute the rational multi-grid M A,a i . Then we search for the convexel c containing b i 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 define the ARA for this pair of points. We note that each operation is done in constant time.

Generally, a 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 first compute the ARA corresponding to the two first 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 such 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 a i . Practically, it may be possible that the locus of rotation of a i intersects twice a line of M A,a i without intersecting any other line of M A,a i as illustrated in Figure 10. In other words, there may exist four intersections between a convexel and the locus of rotation. We perform experiments in order to evaluate the frequency of such a case. We randomly generated ten thousands Pythagorean axes and for each axis we randomly generated ten thousands grid points. For each grid point 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 has four intersections is 1 62500 . Since this particular case does not often occur, we choose to not take it in consideration in our algorithm. Note that in 2D such a case does not happen and the proof can be found in [START_REF] Thibault | Rotations in 2D and 3D discrete spaces[END_REF].

Conclusion and future works

In this paper, we extended the notion of hinge angles, introduced for 2D discrete rotation in [START_REF] Nouvel | Rotations discr tes et automates cellulaires[END_REF][START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF] to the 3D. Extension of hinge angles from the 2D to 3D space involves many problems because most of properties of 2D hinge angles 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 multi-grid that is the intersection of the 3D half-grid and a rotation plane. By redefining the hinge angles on the rotation plane, which are the extension of hinge angles for the rotation in 2D, we showed a subset of the multi-grids where all parameters are rational, called rational multi-grids. This rational multi-grid 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 to extend the search of a 3D rotation from pairs of discrete points introduced for 2D in [START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF].

In Section 6, we introduced the admissible rotation axis that is the set of all possible rotation axes for a pair of digital images. 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 centers of rotations for n pairs of points. Thus one of our future work will be to extend it in 3D and then in nD. Since the extension of [START_REF] Rodriguez | Generalized perpendicular bisector and circumcenter[END_REF] is not done yet, we have proposed an empirical method to find an admissible rotation axis. This method may not offer any formal proof, and we will show that this method is reliable in future.

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; 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 ndimensional rotations into planar rotations[END_REF]. Indeed, a study of such decomposition of an nD discrete rotation is required for the extension of the work presented in this paper to nD discrete rotations.

Figure 1 :

 1 Figure 1: All hinge angles in the first 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 2 :

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

Figure 3 :

 3 Figure 3: Parallel lines of a set L i and geometric interpretation of their parameters

Figure

  Figure4: The five different 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].

Figure 5 :

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

Figure 8 :

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

6 :

 6 and add α(p x , p y , p z , ±i, k) to the list T ;

Figure 9 :

 9 Figure 9: 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 set of admissible approximations for r .

Figure 10 :

 10 Figure 10: In a rational multi-grid M A,a i a convexel can have four intersections with the locus of rotation of p.

  • , α 2 ≈ 15.59 • , α 3 ≈ 21.31 • , α 4 ≈ 50.52 • . The hinge angle α 5 is obtained by symmetry such that α 5 = π 2α 1 ≈ 76.32 •

  Rotation of a 3D image around a Pythagorean axis whose direction is A by a Pythagorean angle θ. Input: An image I, a Pythagorean vector A, a Pythagorean angle θ Output: A rotated image I ′ 1: for all points p in I do

	Algorithm 1 2:
	Because d(p) ≤ hinge angles for any point in the image is 6 √ 3m, the upper bound of the number of √ 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 ).

The discretization is done by applying the rounding function to the result.

Acknowledgment

We would like to thank Laurent Fuchs and Éric Andrès at University of Poitier, France, and Bertrand Nouvel at the National Institute of Informatics, Japan, for the useful discussions.