Aurélie Richard

Gaëlle Skapin

Marc Rodríguez

Eric Andres

Laurent Fuchs

Jean-Serge Dimitri Ouattara

Properties and Applications of the Simplified Generalized Perpendicular Bisector

Keywords:

. The SGPB has some interesting properties that we explore. We show in particular that the SGPB can be used for the recognition and exhaustive parameter estimation of noisy discrete circles. A second application we are considering is the error estimation for a class of rotation reconstruction algorithms.

Introduction

In this paper we are discussing properties and applications of the Simplified Generalized Perpendicular Bisector (SGPB) that has been introduced in [START_REF] Rodríguez | Generalized perpendicular bisector and circumcenter[END_REF][START_REF] Andres | Generalized perpendicular bisector and exhaustive discrete circle recognition[END_REF]. The SGPB is an extension of the classical notion of perpendicular bisector. The Generalized Perpendicular Bisector (GPB) for two 2D regions A and B is defined as the union of the perpendicular bisectors of all the couple of points (p, q) where p and q are respectively points of the regions A and B. The regions we are going to focus on are pixels. The boundary of a GPB between two pixels is composed of line segments and parabola segments. In order to simplify the computational aspects, the SGPB has been introduced. It is only composed of straight line segments. Contrary to a 2D perpendicular bisector, the GPB and SGPB are not lines but surfaces. Bisectors between points and curves or between two curves have also been discussed in detail in the literature [START_REF] Farouki | Computing point/curve and curve/curve bisectors[END_REF][START_REF] Peternell | Geometric properties of bisector surfaces[END_REF] but to our best knowledge, no definition for the bisector between two coplanar surfaces such as pixels has been proposed before [START_REF] Rodríguez | Generalized perpendicular bisector and circumcenter[END_REF].

The SGPB has some interesting properties and can be useful in several different applications mainly when noisy data is considered. The first application we present is the adaptive pixel SGPB and its use for noisy circle recognition. This is similar to the idea presented in [START_REF] Rodríguez | Adaptative pixel resizing for multiscale recognition and reconstruction[END_REF] for discrete straight lines. The second application concerns the parameter estimation of rotations. For this, the Generalized Perpendicular Symmetry transform is introduced. We use this to examine the rotation estimation algorithm recently proposed by Fontjine et al. [START_REF] Fontijne | Reconstructing rotations and rigid body motions from points correspondences as a sequence of reflections[END_REF] and illustrate its behaviour for noisy data.

The paper is organized as follow: section 2 recalls some properties of the SGPB. The section 3 is devoted to the SGBP and the adaptive pixel size SGBP. In particular some of its properties are highlighted and an application to noisy circle recognition is provided. The next section deals with the symmetry relatively to a SGPB and its application to rotation parameter estimation.

2 Definition and properties of the SGBP [START_REF] Rodríguez | Generalized perpendicular bisector and circumcenter[END_REF][START_REF] Andres | Generalized perpendicular bisector and exhaustive discrete circle recognition[END_REF] As mentioned in the introduction, the Generalized Perpendicular Bisector (GPB) between two 2D regions A and B (see Figure .1(a)) is defined as the union of the perpendicular bisectors of all the couples of points (p, q) where p (resp. q) belongs to the regions A (resp. B). In our case we are considering bounded connected regions since we focus on regions that are pixels (that may have different sizes).

By definition, points of the perpendicular bisector are equidistant to both points p and q. Thus a particular point r on the perpendicular bisector is the center of a hypersphere that passes through p and q. Hence, we can set an alternative definition of the GPB. Let S 1 and S 2 be two bounded connected regions and X an Euclidean point of

R n . Let d imin (X) = min Y ∈Si (d(X, Y )), d imax (X) = max Y ∈Si (d(X, Y )) where d is the usual Euclidean distance. Every Euclidean point X ∈ R n such that: [d 1min (X), d 1max (X)] [d 2min (X), d 2max (X)] = ∅ (1)
belongs to the GPB of S 1 and S 2 (see Figure 1(a)). The boundary of the GPB between two pixels is composed of line segments and parabola segments. The parabola segments can be easily removed by extending the line segments. This defines the Simplified GPB (SGPB) (See Figure 1(b,c)). Another way of considering this is simply to state that the minimum distance to a pixel is approximated by the distance to the closest vertex of the pixel rather than to the closest vertex or edge. In other words, the SGPB of two pixels P 1 and P 2 of size respectively λ 1 and λ 2 is given by the following equations:

SGP B(P 1 , P 2 ) = (x, y) ∈ R 2 , (x -C 2x ) 2 + y -C 2y 2 ≤ (x -F 1x ) 2 + y -F 1y 2 ∧ (x -C 1x ) 2 + y -C 1y 2 ≤ (x -F 2x ) 2 + y -F 2y 2 } where C ix , F ix ∈ { x i + λ i 2 , x i -λ i 2 }, C iy , F iy ∈ { y i + λ i 2 , y i -λ i 2 }.
The SGPB between two pixels P 1 and P 2 is bounded by line segments and half-lines as shown in the Figure 2(b).

Proposition 1. The boundary of 2D-Simplified Generalized Perpendicular Bisector between two pixels (x 1 , y 1 ) of size λ 1 and (x 2 , y 2 ) of size λ 2 is composed by at most 10 line segments and half-lines.

The proof of this proposition is straightforward. Let us introduce some concepts and notations useful for the purpose of this section. Let {D i } i∈ [START_REF] Rodríguez | Generalized perpendicular bisector and circumcenter[END_REF][START_REF] Cheng | Joint rotation between two attitudes in the spherical rotation coordinate system[END_REF] be the border segments and half-lines of SGP B(P 1 , P 2 ). Let {P i } i∈ [START_REF] Rodríguez | Generalized perpendicular bisector and circumcenter[END_REF][START_REF] Audin | Geometry[END_REF] be the characteristic points of SGP B(P 1 , P 2 ): i.e. P i belongs to the border of the SGP B(P 1 , P 2 ) and we organize in such a way that

P i = D i ∩ D i+1 if i ∈ [1; 4] and P i = D i+1 ∩D i+2 if i ∈ [5; 8] (see Figure 2(b)). D 1 , D 10 , D 5 , D 6
are the support lines for the half-lines on the border of the SGP B(P 1 , P 2 ). In fact, as the following lemma shows half-lines are linked two by two:

A first immediate property is given by following lemma:

Lemma 1. Let A(a 1 , a 2 ) of size λ A and B(b 1 , b 2 ) of size λ B be two pixels such as a i ± λ A 2 = b i ± λ B 2 .
The border segments D 1 and D 10 are supported by the same line and this is also the case for segments D 5 and D 6 . The proof of this lemma is straightforward using formula [START_REF] Rodríguez | Generalized perpendicular bisector and circumcenter[END_REF]. It is also easy to see that in case a i ± λ A 2 = b i ± λ B 2 then D 1 and D 10 or D 5 and D 6 are not equal but parallel.

Dual of a SGPB

It is easy to determine if a point belongs to a SGPB but not so for a straight line although the SGPB is defined as a union of straight lines. One way of handling this problem is to consider the dual of a SGPB. The dual (see [START_REF] Dexet | Architecture d'un modeleur discret base topologique d'objets discret et methodes de reconstruction en dimensions 2 et 3[END_REF] for more details) of the perpendicular bisector between Euclidean points A(x a , y a ) and B(x b , y b ) is defined as the point ( x b -xa ya-y b ,

x 2 a -x 2 b +y 2 a -y 2 b 2(ya-y b )
) for y Ay B = 0. In the case where y Ay B = 0 we have a point at the infinity.

Proposition 2. The dual of a SGPB is a convex polygon of at most 8 vertices and 8 edges. At most two vertices may be at the infinite and the corresponding edges of the dual polygon are in this case vertical.

Proof. The boundary of the SGPB is formed of straight lines segments and halflines. The dual is therefore a polygon. It is easy to see that this polygon is convex. Let us consider two points A and B on the boundary of the dual of a SGPB. These two points correspond to two straight lines L A and L B that belong to the SGPB. The intersection point I between L A and L B corresponds to the straight line AB in the dual. The point I is inside the SGPB. All the straight lines passing through I form two bundles delimited by L A and L B with one of them containing only lines that are inside the SGPB since L A and L B are inside the SGPB. The other bundle contains lines that are outside the SGPB. A point outside the segment [AB] and not inside the polygon corresponds to a straight line that does not belong to the SGPB and thus to the bundle that contains lines that do not belong to the SGPB. The points that belong to the segment [AB] correspond therefore to straight lines of the bundle that contain only lines inside the SGPB. All the points of [AB] are therefore inside the polygon and this proves that the polygon is convex. The number of vertices and edges of the polygon is an immediate result of Proposition 1.

Determining if a straight line belongs to the SGPB can therefore easily be determined in O(1).

Adaptive Pixel SGPB and noisy circle recognition

In [START_REF] Rodríguez | Adaptative pixel resizing for multiscale recognition and reconstruction[END_REF] some of the authors considered a straight line recognition method that takes into account noise by locally increasing the size of the various pixels. Using a resizing function [START_REF] Rodríguez | Adaptative pixel resizing for multiscale recognition and reconstruction[END_REF] and a local noise estimator (for example the one proposed by Lachaud and al. [?]), pixel sizes are adaptively increased according to the local perturbation on the straight line. The size increased pixels are called Adaptive Pixels [START_REF] Rodríguez | Adaptative pixel resizing for multiscale recognition and reconstruction[END_REF]. The idea here is to do exactly the same with the SGPB for noisy circle recognition purposes.

First, let us introduce the Simplified Generalized Circumcenter (SGC) of a set of n finite and connected regions S = (S i ) i∈ [1,n] . It is defined as the intersection of the SGPB of every two regions of the set (see Figure 2(a)) :

SGC(S) = i,j∈[1,n],i<j (SGP B(S i , S j )).

Theorem 1. Each point of the Simplified Generalized Circumcenter corresponds to the center of at least one circle that intersects all the Adaptive Pixels.

Proof. The proof is similar to the one of Theorem 1 presented in [START_REF] Rodríguez | Generalized perpendicular bisector and circumcenter[END_REF]: Each point of the SGC obtained by the intersection of adaptive pixel SGBP s corresponds to an intersection of radii intervals (dimension 1) for every two adaptive pixels. A direct application of Helly's Theorem tell us that there exists at least one common radius to all these intervals and thus at least one circle of this radius centered on the SGC point that intersects all the adaptive pixels. Proposition 3. Let us consider a set of adaptive pixels P i of various sizes. The dual of all the straight lines crossing the duals of all the SGPB of every pair of pixels P i and P j is the dual of the Simplified Generalized Circumcenter. This proposition is an immediate consequence of the definition of the SGC and the definition of our notion of dual.

Application to noisy circle recognition

On the Figure 4, we can see a Bresenham circle of radius 5 with misplaced and missing pixels. The SGC can be seen in the middle of the different adaptive pixels. One circle example (the center is marked by a black dot in the SGC) is shown. Each pixel size is increased according to a local noise estimator. The algorithm is the same as for regular circles presented in [START_REF] Rodríguez | Generalized perpendicular bisector and circumcenter[END_REF]. The SGP B of each couple of pixels (with the new sizes) is computed and their intersection provides a set of possible circle centers, the Simplified Generalized Circumcenter (SGC). A further computation provides for each point of the SGC, the interval of possible circle radii that correspond to circles intersecting all the size increased pixels.

Generalized Reflection Symmetry and biased rotation parameter estimation

As a second application, we show how the SGPB helps to visualize the errors that occur when noisy data are used in the rotation reconstruction algorithm recently proposed by Fontijne et al. [START_REF] Fontijne | Reconstructing rotations and rigid body motions from points correspondences as a sequence of reflections[END_REF]. For this we are going to introduce a generalization of the reflection symmetry.

Bundles and Strips

We define the line bundle F i as the set of lines ∆ that pass through P i with slopes in the interval [slope(D i ), slope(D i+1 )] (or [slope(D i+1 ), slope(D i+2 )] according to the values of i) as shown in the Figure 2(b). We also define a line strip strip(P i , P j ) by the set of parallel lines bordered by ∆ i ∈ M passing through P i for i ∈ [1; 4] and ∆ j ∈ M passing through P j for j ∈ [5; 8] as shown in the figure 2

(b).

A line in the SGPB can either be characterized as a point inside the dual of the SGPB or by using bundles F i and strips. Considering a line D in the SGP B two cases are possible. In the first case, the line D passes through P i then D belongs to F i . In the second case the line D does not pass through P i then there exists P j and P k such that a translate D ′ of D passes through P j and another translate D ′′ of D passes through P k . In this latter case, D belongs to strip(P j , P k ) (see figure 2(b)). This leads to a second characterization of a line inside a SGPB: Property 1 (Lines in the SGPB) Every line in the SGPB belongs either to a bundle F i or is included in a strip strip(P k , P l ) Furthermore, we have the additional property:

Property 2 The union of the 8 bundles covers the GPB.

Proof. Let X a point of the GPB. Since the border of the GPB is composed of two convex polylines ∂C 1 and ∂C 2 , there exists two lines ∆ 1 and ∆ 2 which passes through X and respectively p 1 ∈ ∂C 1 and p 2 ∈ ∂C 2 . Obviously, p 1 and p 2 ∈ P i .

Generalized Reflection Symmetry for a pixel. When considering the perpendicular bisector B of two points p and q, it is obvious that p is the image of q and q the image of p through a reflection symmetry of axis B. Let us define the Generalized Reflection Symmetry GRS(P 1 , P 2 ) of a point as the union of the reflection symmetries of axis the straight lines in the SGPB of two pixels P 1 and P 2 . Not surprisingly, the image of a pixel P 3 by a GRS(P 1 , P 2 ) is not a pixel. Let us characterize the obtained region denoted C GRS(P1,P2) (P 3 ) (See figure 6). By construction, a point p of C GRS(P1,P2) (P 3 ) is the symmetric of a point p 3 belonging to the pixel P 3 with respect to a line ∆ that belongs to SGP B(P 1 , P 2 ). More formally, this can be written as:

C GRS(P1,P2) (P 3 ) = {p ∈ R 2 /∃p 3 ∈ P 3 , ∃∆ ∈ SGBP (P 1 , P 2 ), p = sym(p 3 , ∆)}
where sym(p 3 , ∆) denotes the reflection symmetry of the point p 3 with respect to the line ∆. Hence, the following characterization of a point belonging to C GRS(P1,P2) (P 3 ) is obtained: Definition 1 (GRS of a pixel). Let P 1 and P 2 be two pixels centered in (x 1 , y 1 ) and (x 2 , y 2 ) and the associated SGP B(P 1 , P 2 ). Let C GRS(P1,P2) (P 3 ) be the image of the pixel P 3 by GRS(P 1 , P 2 ). A point (x, y) belongs to C GRS(P1,P2) (P 3 ) if the following assertions are true:

(i) ∃(x ′ , y ′ ) such that (x 3 -0.5 ≤ x ′ ≤ x 3 + 0.5) and (y 3 -0.5 ≤ y ′ ≤ y 3 + 0.5) (ii) (x, y) can be written as

( -2ac-a 2 x ′ +b 2 x ′ -2aby ′ a 2 +b 2 , -2bc-2abx ′ +a 2 y ′ -b 2 y ′ a 2 +b 2 ) such that (x 2 , -c-ax2 b ), (x 1 , -c-ax1 b
), ( -by1-c a , y 1 ), ( -by2-c a , y 2 ) satisfy the distance law [START_REF] Rodríguez | Generalized perpendicular bisector and circumcenter[END_REF].

Let us detail a way to construct the GRS of a pixel. The image of a pixel P 3 by a complete 1 line bundle passing through a point p is a ring centered in p defined by the circles centered in p with radii d min (p, P 3 ) and d max (p, P 3 ). In our case, the bundle is not complete 2 , the result is then a piece of ring (see Figure 5(a)). Let P i and P j be the images of the pixel P by the two lines which define a strip. The image of the pixel P by this line strip is the strip bounded by the images of the pixel relatively to the line passing through the centers of P i and P j as shown in the Figure 5(b).

Let P 3 be the pixel of which the image is to be computed. The lines D j are the lines bounding the SGP B(P 1 , P 2 ) and P i the characteristic points of the SGP B(P 1 , P 2 ). Let I j be the images of P 3 by D j . Since D 1 = D 10 and D 5 = D 6 (lemma 1), there are only eight different images. Then we consider the images of P 3 by the whole line bundle F i that passes through P i = D l ∩ D m . These images are the ring pieces starting with I l , ending with I m of internal radius d min (P i , P 3 ) and external radius d max (P i , P 3 ). We therefore obtain eight ring 1 Complete means all the line passing through the point. 2 That means the slopes of the lines are included between slopemin and slopemax.

pieces that form a closed object (since D 1 = D 10 and D 5 = D 6 ) as shown in the Figure 6(b).

The last building step consist in integrating the images of P 3 by all the possible strips that are included in the SGP B(P 1 , P 2 ). This leads to a second characterization of the shape C GRS(P1,P2) (P 3 ) :

Characterization 1 (GRS of a pixel) A point belongs to the image of a pixel P 3 by the GRS(P 1 , P 2 ) between the pixel P 1 and P 2 if and only if it belongs to the set points defined by the eight pieces of rings centered in P i for 1 ≤ i ≤ 4 of radius d min (P i 1≤i≤4 , P 3 ) and centered in P i for 5 ≤ i ≤ 8 of radius d max (P i 5≤i≤8 , P 3 ), where the P i are the characteristic points of the SGBP.

The angular arc between the extrema of the region (as defined below) C GRS(P1,P2) (P 3 ) gives us an upper bound of "noise" that is added to a pixel when its image by a GRS is computed. Let consider the lines ∆ 1 and ∆ 2 that are the supports of the half-lines D 1 and D 5 (see Figure 6 Then, the angular length between the extrema of C GRS(P1,P2) (P 3 ) is given by (β + 2θ)( R1+R2

2

), β and θ respectively denote the angles between the lines Oa and Od and the lines ∆ 1 and ∆ 2 .

Rotation reconstruction using the SGPB

In this section, the rotation reconstruction algorithm proposed by Fontijne et al. [START_REF] Fontijne | Reconstructing rotations and rigid body motions from points correspondences as a sequence of reflections[END_REF] is adapted to noisy data using the SGP B. This algorithm is based on the Cartan-Dieudonné theorem [START_REF] Audin | Geometry[END_REF] that decomposes nD rotations into reflections. Methods based on this approach [START_REF] Aragon-Gonzalez | Reflections, rotations, and pythagorean numbers[END_REF][START_REF] Fontijne | Reconstructing rotations and rigid body motions from points correspondences as a sequence of reflections[END_REF] need exact point correspondences and suffer from noisy data.

However, rotation reconstruction methods from point correspondences are used in many application domains such as Computer Vision or body movement analysis. Different approaches can be found in the literature [10, 11, ?]. Their drawbacks are the impossibility to extend them in dimension higher than two or the subtle geometrical model and computations they need. A simple algorithm based on the decomposition of nD rotations into planar rotations has been developed by some of the authors and it has been shown to be rather robust against noisy data [START_REF] Richard | An algorithm to decompose n-dimensional rotations into planar rotations[END_REF][START_REF] Richard | Decomposition of nDrotations: classification, properties and algorithm[END_REF]. But, contrary to the Fontijne et al. [START_REF] Fontijne | Reconstructing rotations and rigid body motions from points correspondences as a sequence of reflections[END_REF] algorithm, this approach is not well suited when the data are acquired incrementally.

The purpose of this section is to give some hints on how the SGPB could help to develop an adapted version of the Fontijne et al. [START_REF] Fontijne | Reconstructing rotations and rigid body motions from points correspondences as a sequence of reflections[END_REF] algorithm to noisy data. For this first attempt, only dimension two is explored but as SGPB is defined for any dimension, extension to nD is expected.

by n-dimensional pixels and to introduce the SGPB into the reflection determination process. The expected result is a class of rotations that fit the input data.

Principle of the algorithm (see Figure 7(b)). First the points p i and their corresponding points p ′ i are replaced by n-dimensional pixels P i and P ′ i . Then, the SGPB M 1 = SGP B(P 1 , P ′ 1 ) (corresponding to R 1 in the exact case) is determined. Using M 1 , the region C GRS(P1,P ′ 1 ) (P 2 ) is computed. This region corresponds to R 1 (p 2 ) in the exact case.

In order to restrict C GRS(P1,P ′ 1 ) (P 2 ), we use the property that a circle centered on p ′ 1 passing through p ′ 2 and R 1 (p 2 ) exists (see figure 7(a)). In the pixel case this means that we have to consider the set S of all circles centered on points of P ′ 1 passing through P ′ 2 . So region of interest for the next step is C GRS(P1,P ′ 1 ) (P 2 ) ∩ S. Now, the SGBP (C GRS(P1,P ′ 1 ) (P 2 ) ∩ S, P ′ 2 ) is computed (see Figure 7(b)). The obtained region must also be restricted with the set L of lines that pass through P ′ 1 while in the exact case the bisector R 2 passes through p ′ 1 (see Figure 7(a)). This finally determines the region M 2 = SGBP (C GRS(P1,P ′ 1 ) (P 2 ) ∩ S, P ′ 2 ) ∩ L corresponding R 2 in the exact case. In higher dimension, this process is continued with next point correspondences.

At the end of the algorithm we obtain a class of rotations that are determined by lines chosen in the sequence of the M i , 1 ≤ i ≤ n.

At present, we only have experienced this process in dimension 2 and we face efficiency problems for computing the SGBP between a pixel and a complex region such as the C GRS(P1,P ′ 1 ) (P 2 ). Currently, we are working on this efficiency problem and also on the extension of this algorithm to higher dimension.

Conclusion

In this paper, we studied the Simplified Generalized Perpendicular Bisector and its properties: whereas the classical euclidean 2D perpendicular bisector (and all the commonly used medial axis notions), the SGPB of two regions is a 2D surface that collects all the perpendicular bisector of any couple of points of both regions. We have particularly defined the SGPB between two pixels of different sizes. The SGPB is useful to obtain an exhaustive parameter estimation of noisy circles. We have also characterized the lines in the SGPB by showing that the dual of a SGPB is a convex polygon and we have defined a new operation: the Generalized Reflection Symmetry which is a symmetry relatively to a SGPB. This operation allows us to adapt the algorithm of Fontijne et al. [START_REF] Fontijne | Reconstructing rotations and rigid body motions from points correspondences as a sequence of reflections[END_REF] to reconstruct rotations while taking noisy data into account. In future works, the GPB will be further studied. There are various notions of discrete bisector that have been proposed. They are, among other applications, used to analyze and filter medial axis [START_REF] Couprie | Discrete bisector function and euclidean skeleton in 2d and 3d[END_REF][START_REF] Talbot | Euclidean skeletons and conditional bisectors[END_REF] where the medial axis of a Jordan curve is in any point equidistant to its borders. These notions of discrete bisectors are discrete curves. How these notions are linked to the Generalized Bisector is one the question we would like to investigate. Especially the links with Voronoï diagrams, medial axis and skeletons seems promising. An extension in higher dimension will also be investigated.

Fig. 1 .

 1 Fig. 1. The exact perpendicular bisector (a) and the approximation (b) where the parabolic pieces have been dropped by extending the straight lines (c). The generalized bisector is slightly reduced.

Fig. 2 .

 2 Fig. 2. (a) Simplified Generalized Perpendicular Bisector and Simplified Generalized Circumcenter of three adaptive pixels ; (b) The 6 border segments and the 4 half-lines of a SGPB {Di} i∈[1,10] and its characteristic points {Pi} i∈[1,8] (D1 = D10 and D5 = D6). The line bundle F8 and the line beam beam(P4, P5).

Fig. 3 .

 3 Fig. 3. Dual of Three SGPB corresponding to three pixels.

Fig. 4 .

 4 Fig. 4. Circle parameters exhaustive estimation with noisy data.

P 1 dm 1 = (x 1 , y 1 )Fig. 5 .

 11115 Fig. 5. (a) Image of the pixel P2 by a partial bundle from m1; (b) image of the pixel P by the strip(∆i, ∆j)

Fig. 6 .

 6 Fig. 6. (a) Generalized Reflection Symmetry of a pixel (12,16) relatively to a SGPB between the pixel (0,0) and (3,4); (b) Image of a pixel by a GRS. The obtained shape is composed of ring pieces; (c) Symmetric of a pixel P3 with respect to the lines ∆1 and ∆2.

  (c)). They intersect at the point O = ( x1+x2 2 , y1+y2 2 ). On Figure 6(c), the symmetric squares ABDC and EF GH of the square abdc by the lines ∆ 1 and ∆ 2 are the extrema of the region C GRS(P1,P2) (P 3 ). As the lines ∆ 1 and ∆ 2 belong to the bundle of lines that pass through the point O, the squares ABDC and EF GH are contained in the ring centered in O formed by the circles C 1 and C 2 of radii R1 = d min (O, P 3 ) and R2 = d max (O, P 3 ). Since B and H belong to C 1 and C, F belong to C 2 , the points A, D, G and E belong to the circle centered in O of radius R3 = R1+R2 2 .

In this section, Ri denotes indistinguishably the reflection and the line that determines the reflection.

Rotation reconstruction algorithm. To reconstruct the rotation R from n points p i and their correspondences p ′ i = R(p i ) the algorithm proposed in [START_REF] Fontijne | Reconstructing rotations and rigid body motions from points correspondences as a sequence of reflections[END_REF] works by finding successive reflections using the Cartan-Dieudonné theorem [START_REF] Audin | Geometry[END_REF]. Hence, for each i, we have

, where R i are the reflections 3 .

The R i are determined incrementally as perpendicular bisectors. The first reflection R 1 is determined as the perpendicular bisector of p 1 and p ′ 1 , the second reflection R 2 is the perpendicular bisector of R 1 (p 2 ) and p ′ 2 . By construction, the reflections have the property that their composition does not move previously aligned points. For example, as 7 illustrates the algorithm in dimension two. It must be noted that p ′ 1 = R 1 (p 1 ) is on the line R 2 and it is the center of circle passing through R 1 (p 2 ) and p ′ 2 . Now, replacing the perpendicular bisectors with the SGPB could help to adapt the algorithm of Fontijne et al. [START_REF] Fontijne | Reconstructing rotations and rigid body motions from points correspondences as a sequence of reflections[END_REF] to noisy data.

Adaptation of the rotation reconstruction algorithm to noisy data. The algorithm described in the previous section supposes that exact point correspondences are provided. However, in the case of noisy data, the reflections became hard to determine. To deal with noisy data, the idea is to replace points


A Richard

L Fuchs

E Andres

G Largeteau-Skapin

Decomposition of nD-rotations: classification, properties and algorithm

Keywords: nD-dimensional rotations, planar rotations, Cartan-Dieudonné theorem, Schur decomposition, classification of nD-rotation

In this paper, the decomposition of nD-rotations is studied. Using this decomposition, nD-rotations are classified and properties are underlined. For each class, an algorithm previously presented by the authors to decompose nD-rotation into planar rotations is proposed and a generalized algorithm is presented. Moreover, an alternate algorithm based on the Schur decomposition is investigated. A comparison between these different algorithms is finally provided.

Introduction

Rotation characterization is very useful in many application domains such as in computer vision [START_REF] Shapiro | Computer Vision[END_REF], in medicine [START_REF] Pl | Determination of axial rotation angles of limb segments -a new method[END_REF][START_REF] Pl | Joint rotation between two attitudes in the spherical rotation coordinate system[END_REF], in Procrustes analysis [START_REF] Gower | Procrustes problems[END_REF] or in global visibility computation [START_REF] Charneau | Exact, robust and efficient full visibility computation in the Plücker space[END_REF] for example. In many cases, estimating 3D rotations is sufficient but higher dimensions are needed in some cases (color image analysis, global visibility computation using Plücker space for example). rajouter une phrase sur le modele conforme From a set of points and their images, our aim is to characterize the rotation parameters (determination of the angles and planes) in nD. Since, our framework is experimental, the method must be somewhat robust to noise. This work is supposed to be a basis for the description of nD discrete rotations by extension of hinge angles for 3D discrete rotations [START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF]. In this paper, we suppose without loss of generality that the rotations are centred on the origin of the frame and that the angles are in [0; π].

Unfortunately, as explained in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], in the available literature there are not many methods that determine the characteristics of a nD-rotation. Some methods can be used only in 3D [START_REF] Pl | Determination of axial rotation angles of limb segments -a new method[END_REF][START_REF] Pl | Joint rotation between two attitudes in the spherical rotation coordinate system[END_REF][START_REF] Ga | Computing Helmert transformations[END_REF]. Other methods [START_REF] Perwass | Numerical evaluation of versors with Clifford Algebra[END_REF][START_REF] Dorst | Determining a versor in nD geometric algebra from known transformation of n vectors[END_REF] estimate the rotor which describes the rotation. In [START_REF] Perwass | Numerical evaluation of versors with Clifford Algebra[END_REF], the rotation planes cannot be computed easily.

Moreover, in [START_REF] Dorst | Determining a versor in nD geometric algebra from known transformation of n vectors[END_REF], the versor can be determined only with exact data. The constructive proof of the Cartan-Dieudonné theorem proposed in [START_REF] Aragon-Gonzalez | Reflections, rotations, and pythagorean numbers[END_REF] could provide an algorithm to characterize nD rotations by their decomposition into planar rotations. Unfortunately it has been shown that this solution is not practicable [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF].

In [START_REF] Fontijne | Reconstructing rotations and rigid body motions from points correspondences as a sequence of reflections[END_REF], Fontijne et al. proposed a method to reconstruct versors using vector correspondences. This method is based on the construction of a reflection sequence which aligns the points with their correspondences one by one. As explained by its authors, this method is not well suited for noisy data [START_REF] Richard | Properties and applications of the simplified generalized perpendicular bisector[END_REF]. In this case, retrieving the reflexion axis could be very difficult. Even if the found axis is in accordance with the data, due to the multiple necessary decisions, the obtained versor could be significantly different from the correct one.

Another method can be deduced from the Schur decomposition [START_REF] Golub | Matrix computations[END_REF]15] which factorizes a matrix A into two matrices (an orthogonal matrix U and a blockdiagonal matrice T such that A = U T U t ). As far as we know, the Schur decomposition has not been used to decompose nD-rotations into planar rotations. Since few methods exist, the authors have proposed a new algorithm to determine nD rotations from n points and their images [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF]. From the Cartan-Dieudonné theorem [START_REF] Audin | Geometry[END_REF], we know that a nD rotation can be decomposed into an even number of planar rotations. In [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], some of the authors have proposed an algorithm that finds these planar rotations and the angle for each of them. This fully characterizes the rotation.

The paper is organized as follow. In Section 2, our method presented in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] for determining rotation parameters is summarized and its limits are underlined. In Section 3, nD-rotations are classified and an improvement of our algorithm is proposed. Another method based on the Schur decomposition is investigated in the next section. Finally, a comparison between our algorithm, its improvements and the Schur decomposition is provided.

2. The nD rotation decomposition algorithm presented in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] In [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], the authors have proposed an algorithm to decompose the nD-rotations. It works in the Geometric Algebra framework [START_REF] Perwass | Numerical evaluation of versors with Clifford Algebra[END_REF][START_REF] Dorst | Geometric algebra for computer science: an object oriented approach to geometry[END_REF][START_REF] Hildebrand | Geometric algebra and its application to computer graphics[END_REF][START_REF] Sommer | Geometric computing with Clifford Algebra[END_REF][START_REF] Charneau | Étude et application des algèbres géométriques pour le calcul de la visibilité globale dans un espace projectif de dimension n ≥ 2[END_REF]. However, there are particular classes of rotations that are not handled by the proposed algorithm: the isoclinic rotations and what we called the pseudo-isoclinic rotations. We will therefore extend this algorithm to integrate those particular cases. We briefly present the useful notions for the understanding; for more details the reader may refer to [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF].

In the following, the n dimensional space E (usually R n ) is considered with a basis (not necessarily an orthogonal one) denoted by (x i ) i=1...n . The space E must be equipped with a non-degenerated quadratic form and the associated bilinear form is denoted by B. In standard cases, B corresponds to the usual scalar product. Thus, the space E is a metric space upon which isometries can be defined.

Some results about nD-rotations

Let R be a rotation in SO(E), the rotation group of E. The fixed-point subspace of R is denoted by F . Its orthogonal subspace in E is denoted by P or F ⊥ . In other words,

F = {x ∈ E|R(x) = x} and F ⊥ = {x ∈ E | ∀y ∈ F, B(x, y) = 0}.
The respective images of (x i ) i=1...n by R are denoted by (y i ) i=1...n . Since B is non-degenerate, the space E can be decomposed into E = P ⊕ F .

In the following of this paper, the following notations are used:

* the orthogonal subspace to F is denoted by P * the rotation planes are denoted by P i * P i is used in the R-isogonal subspace decomposition provided in Section 3.3

Moreover, E can also be decomposed into another direct sum given by the following well-known proposition [START_REF] Audin | Geometry[END_REF]: Proposition 1. Let f be a linear isometry of E. The vector space E is an orthogonal direct sum :

E = V ⊕ W ⊕ P 1 ⊕ • • • ⊕ P r
where the subspaces V , W and P i are stable under f , and

f | V = Id V , f | W = -Id W and every P i for i = 1, . . . , r is a plane such that f | Pi is a rotation.
In term of matrices, Proposition 1 shows us that there exists an orthonormal basis of E such that the matrix of f can be written as

F =        I p -I q R θ1 . . . R θr       
where I p denotes the identity matrix of order p and R θi (i = 1, . . . , r) are 2 × 2 planar rotation matrices. As we deal with rotations, W has an even dimension q. Thus, we can consider the submatrix -I q as the concatenation of q/2 rotations of angle π. So, we can set V = F , W = 0 and

P = P 1 ⊕ • • • ⊕ P r .

First version of our algorithm [7]

The method we proposed in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] consists of three steps. We begin by the determination of P , the orthogonal subspace to the rotation fixed-points set. From it, we can construct the rotation planes and finally we can retrieve the rotation angles.

Determination of P .

The invariant subspace F of a transformation R and its orthogonal P are in direct sum. Moreover, the null space (denoted by Ker) and the range (denoted by Im) of an application are in direct sum too. From the two equalities E = P ⊕ F and E = Ker(R -Id) ⊕ Im(R -Id) we deduce that P is generated by (y i -x i ) i=1...n . The complete proof is provided in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF].

x 1 x 2 y 2 P 1 θ x 2 -y 2 x 1 -y 1 x 1 y 1 y 1 θ
Figure 1: Graphical illustration of rotation plane and angle. x 1 , y 1 (resp. x 2 , y 2 ) represent a vector and its image by a rotation. The rotation plane P 1 is generated by x 1 -y 1 and x 2 -y 2 . The angle between x 1 and y 1 is the same than those between x 1 and y 1 in P 1 .

Determination of the rotation planes.

The goal here is to build a basis of P in order to be able to determine the rotation planes. Using Proposition 1, the subspace P can be decomposed as

P = P 1 ⊕ • • • ⊕ P r , i = 1, . . . , r.
Then, if we compute an orthogonal basis of P , two elements of basis will determine

P i . Actually, let (b 1 , b 2 , • • • , b 2r-1 , b 2r
) be an orthogonal basis of P . In other words,

P = b 1 ∧ b 2 ∧ • • • ∧ b 2r-1 ∧ b 2r .
The plane P i is generated by b 2i-1 ∧ b 2i for i = 1...r. The most obvious choice to compute it would be to use the Gram-Schmidt orthogonalization process over the vectors y i -x i . However, in case of noisy data, it is difficult to obtain a matrix with correct rank1 . This is very important in case of real data. It is one of the main reasons why we worked on this method rather than simply using the Schur method as it can be seen in Section 4. We decided to use an alternative which is to use a SVD decomposition over the matrix formed by the vectors y i -x i . The approximation property of the SVD [START_REF] Golub | Matrix computations[END_REF] ensures that we obtain a matrix with rank 2r that is the dimension of P even in case of slightly noisy data. Hence, following the usual presentation of the SVD decomposition, the first 2r columns form the desired basis.

Determination of rotation angles.

The last step of the algorithm finds the rotation angle of each planar rotation. For each rotation plane P i , the rotation angle between x and y is the same than the angle between their projections x and y in P i (cf Figure 1).

x ⊥ x x B x B
Then, for each P i we have to determine these two vectors x and y . To do so, we suppose that the vectors x and y are chosen such that their projections onto P i are not null.

A collinear vector to x is given by (x P i ) P i where (. .) denotes the orthogonal contraction2 (see [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF]). A graphical illustration of the orthogonal contraction of a vector into a plane is provided in Figure 2.

Then, the rotation angle θ in P i is the angle between x and y and is given by: cos(θ) =

x .y ||x ||||y || and sin(θ) = x ∧ y ||x ||||y ||P i .

The computation is applied to all rotation planes in order to find all the rotation angles. Unfortunately, as shown in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], this algorithm is not adapted in some particular cases.

Limits of this algorithm

In 4D, a rotation can be decomposed into one or two planar rotations. These 4D-rotations are said to be respectively simple and double rotations. They leave fixed respectively a plane and one point (the rotation center). If the angles θ 1 and θ 2 of a double rotation have the same magnitude (|θ 1 | = |θ 2 | ∈ [0, π]) then the 4D rotation is said to be isoclinic. This is also called a Clifford displacement [START_REF] Mebius | Applications of quaternions by dynamical simulation, computer graphics and biomechanics[END_REF].

More generally, we have the following definition:

Definition 1 (Isoclinic rotations). In a n-dimensional space, with n even, a rotation is said isoclinic if all its n 2 angles are equal (up to the sign). These particular rotations are only defined in an even dimension space. Indeed, in an odd dimension space, the unique isoclinic rotation is the identity transformation. In [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], the term pseudo-isoclinic was introduced as follows:

Definition 2 (Pseudo-isoclinic rotations). A nD-rotation with at least two equal angles is called pseudo-isoclinic. This definition works both in an even and in an odd dimension space.

In [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], we showed that our algorithm is not adapted to these particular cases. Actually, if the data represents these two kinds of rotations, the rotation parameters are not well determined by the algorithm such as it is presented in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF]. In [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], numerical experiments have been conducted with dimension strictly smaller than 7. In this case, four groups of pseudo-isoclinic rotations can be distinguished:

(a) In 4D, the rotation has 2 equal angles.

(b) In 6D, the rotation has 3 equal angles.

⇒ the rotation is said isoclinic.

(c) In 5D, the rotation has 2 equal angles.

(d) In 6D, the rotation has exactly 2 equal angles.

An adaptation of the algorithm to these cases is proposed in Section 3. There exist previous works on isoclinic rotations [START_REF] Mebius | Applications of quaternions by dynamical simulation, computer graphics and biomechanics[END_REF]: the rotation planes are given by orthogonal planes generated by a vector and its image. Our algorithm uses this to handle cases (a,b). To the best of the knowledge of the authors case (c) has not yet been addressed. Nevertheless, some rotations in an odd dimension space (see Section 3.1.1) can be handled like isoclinic rotations in an even dimension space. An enhancement is also provided in case (d). It is based on the projections onto the range. These improvements are detailed case-by-case in the following section.

This list is only exhaustive for rotations in a space of dimension strictly smaller than 7. Actually, in higher dimension, other cases have to be considered. For example, a 9D-rotation can be decomposed into four rotations having four angles θ 1 , θ 2 , θ 3 , θ 4 with θ 1 = θ 2 = θ 3 = θ 4 . This pseudo-isoclinic rotation is composed of a couple of isoclinic rotations. This underlines that the term pseudo-isoclinic, as it is defined in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], gathered various cases. A classification of these pseudo-isoclinic rotations can thus be proposed. In Section 3.3 a general method is proposed.

Extension to the pseudo-isoclinic rotations

In this section, a classification of the nD-rotations as well as some new definitions are proposed. According to the best of the knowledge of the authors, it is the first time that such a classification is proposed. For each class of rotations, a decomposition algorithm is presented. Finally, a generalized decomposition algorithm which can be applied to all pseudo-isoclinic rotation is supplied.

Classification and properties of nD-rotations

In the following, R denotes a nD-rotation which is not the identity transform.

Classification of nD-rotations

The nD-rotations can be classified with respect to their decomposition. Let us recall or even introduce some new terms useful for the proposed classification.

The isoclinic notion, explained in the previous section, is only defined in an even dimension space. In an odd dimension space, a similar notion can be defined.

Definition 3 (Pure pseudo-isoclinic rotations). In an odd dimension space, the nD-rotations, with all its n-1 2 angles equal are called pure pseudo-isoclinic rotation.

The (pseudo) isocliny property is defined for angle equality. Since the data is noisy, it is necessary to decide when two angles are equal. This problem is a recurrent one when working with real numbers. In order to have a partition of the finite set of the angles we introduce the following relation.

Let θ 1 and θ 2 be two angles. Let β be the approximation parameter and let =β be the following relation:

θ 1 =β θ 2 ⇔ |θ 1 | -|θ 2 | ≤ β
For our experiments we arbitrarily fix the value of β to π 12 . Moreover, since β is given and to avoid cumbersome indexes, we simply set = in order to denote =β in the following.

Let us introduce the new term pseudo-isogonal to characterize the relation =. This relation is reflexive and symmetric but not transitive so = is not an equivalence relation. The idea is to use this relation to form a partition of the angle set. The following definition is required: Definition 4 (Pseudo-isogonal component). The pseudo isogonal component (PIC) of an angle θ denoted by [θ] is defined as the set of the angles α ≤ θ in relation with θ. In other words :

[θ] = {α ∈ [0; θ] | α =θ}.
The pseudo-isogonal component is not an empty set. Actually, from the reflexivity property, an angle belongs always to its pseudo-isogonal component.

Let S be the ordered set of the angles (in decreasing order).

Definition 5 (Pseudo-isogonal partition).

A pseudo-isogonal partition (PIP) is a set of pseudo-isogonal components such that it forms a partition of the angle set; that means:

∀θ i ∈ S, ∃![θ j ] ∈ P IP | θ i ∈ [θ j ].
Algorithm 1 Compute the pseudo-isogonal partition of a rotation Input: S (ordered set of angles of the rotation R)

Output: C (set of the pseudo-isogonal component of R) C = ∅ while S = ∅ do Θ : first element of S Compute [Θ] C ← C ∪ [Θ] S ← S\[Θ] end while
Remark 1. The following property holds: P IP ⊂ P IC. There is no equality between both sets because of the non transitivity of the relation =.

The method to compute a pseudo-isogonal partition of a rotation R is summarized by the algorithm 1.

Definition 6 (Length of the pseudo-isogonal component). The length of the pseudo-isogonal component denoted by l([θ]

) is defined as its cardinal number.

This length is included between 1 and n 2 . For a given rotation, the component number is as well included in the interval [1; n 2 ].

Proposition 2. The pseudo-isogonal component of an isoclinic or a pure pseudoisoclinic rotation is unique and its length is maximal.

Definition 7 (Isogonal space). An even dimensional subspace of E (of dimension at least 4) is said R-isogonal if the restriction of the rotation R to this subspace is an isoclinic rotation.

A nD-rotation R (for n ≥ 8) can have many R-isogonal subspaces. For example, a 9D-rotation with two pseudo-isogonal components of length 2 has two R-isogonal subspaces. An isoclinic rotation R (resp. pure pseudo-isoclinic rotation) has only one R-isogonal subspace. Its dimension is n (respectively n-1

2 ). The R-isogonal subspace is the maximal subspace in P (the orthogonal subspace to the invariant subspace).

Example 1. Let us consider a 11D-rotation which can be decomposed into 5 rotations with angles (θ 1 , θ 2 , θ 3 , θ 4 , θ 5 ) = (130 o , 130 o , 72 o , 68 o , 60 o ). With these new definitions, the following relations and equalities are true:

θ 1 Rθ 2 , θ 3 Rθ 4 , θ 3 Rθ 5 , θ 4 Rθ 5 , [θ 1 ] = {θ 1 , θ 2 } = [θ 2 ], [θ 3 ] = {θ 3 , θ 4 , θ 5 } = [θ 4 ] = [θ 5 ], l([θ 1 ]) = l([θ 2 ]) = 2, l([θ 3 ]) = l([θ 4 ]) = l([θ 5 ]) = 3.
Remark 2. The partition based on the pseudo-isogonal components can be calculated from the input data of our algorithm. Actually for a given matrix A, its SVD decomposes A as U SV t where S is diagonal and its elements 3 , called singular values, are the square root of the eigenvalues of AA t . Here, the SVD is applied to the matrix (A -Id). An elementary calculation shows that the eigenvalues of (A -I)(A -I) t with A be a rotation matrix are given by (2 -2 cos θ i ). Each eigenvalue is of order two. Thereby, the number of angles is the half of the number of non-null singular value number (see Section 3.2.1).

The nD-rotations can be classified according to the number of angles N a and to the number of pseudo-isogonal components N c . This classification is given by the following theorem.

Theorem 1 (Classification of nD rotations). Let R be a nD-rotation. If all the angles of R are different then R is a general rotation else R is a pseudo-isoclinic rotation. In this case, if all the n 2 angles of R are equal and n is even (resp. odd) then R is an isoclinic rotation (resp. pure pseudoisoclinic rotation) in other cases R is an ordinary pseudo-isoclinic rotation.

Properties of nD-rotations

In [START_REF] Mebius | Applications of quaternions by dynamical simulation, computer graphics and biomechanics[END_REF], Mebius defines two particular 4D-isoclinic rotations. The 4D-space is equipped with an arbitrary standard orientation. A 4D-isoclinic rotation has two angles of same magnitude (the angles have the same or an opposite sign). If the two angles have the same orientation then the rotation is said to be left-isoclinic else it is said to be right-isoclinic. In other words:

Rot 4D = Rot 2D (θ 1 )Rot 2D (θ 2 ) with θ 1 = θ 2 or θ 1 = -θ 2 .
Using these definition, we have the following theorem: Theorem 2 (Isoclinic decomposition theorem [START_REF] Mebius | Applications of quaternions by dynamical simulation, computer graphics and biomechanics[END_REF]). Any 4D rotation can be decomposed into two 4D-rotations: one left-and one right-isoclinic rotation. This decomposition can be processed in two ways, differing by a central reversion.

In other words, review 1: since a 4D rotation stands in two orthogonal planes P 1 and P 2 we have:

Rot 4D = [Rot 2D (θ 1 )Rot 2D (θ 2 )] P1 [Rot 2D (θ 3 )Rot 2D (θ 4 )] P2 = [Rot 2D (θ 1 )Rot 2D (θ 3 )] Left-isoclinic [Rot 2D (θ 2 )Rot 2D (θ 4 )] Right-isoclinic with θ 1 = θ 3 and θ 2 = -θ 4 .
The proof of this theorem is provided in [START_REF] Mebius | Applications of quaternions by dynamical simulation, computer graphics and biomechanics[END_REF]. Moreover, Mebius decomposes the 3D rotation into the same way:

Theorem 3.
[21] Any 3D-rotation of angle 2θ can be decomposed in a rightand a left-isoclinic rotation of angle θ.

In other words, review 1: a 3D-rotation can be seen as the composition of rotation in two planes P 1 and

P 2 Rot 3D (2θ) = [Rot 2D (θ)Rot 2D (θ)] P1 [Rot 2D (θ)Rot 2D (-θ)] P2 = R P1 (2θ).
The proof of this second theorem is also provided in [START_REF] Mebius | Applications of quaternions by dynamical simulation, computer graphics and biomechanics[END_REF].

The notion of left-and right-isoclinic rotations make no sense in higher dimensions. Actually, two angles can be said to be equal or opposite but relations between three angles can not be characterized in this way. In order to by-pass the problem, another notion needs to be introduced: Definition 8 (Standard rotations). A nD-rotation with all its angles of same magnitude is said to be a standard rotation (n not necessarily even).

To decompose a rotation into standard rotations, one needs to order the rotation in a convenient way. In 3D, rotations do not commute but there is no isoclinic rotation. Whereas, in 4D, under some special condition on orthogonality, rotations can be permuted and thus placed in the convenient order. Let us recall the notion of completely orthogonal planes: Definition 9 (Completely orthogonal planes). Two planes P 1 and P 2 are said to be completely orthogonal if each line in P 1 is orthogonal to any line in P 2 .

The isoclinic decomposition theorem of Mebius can be extended to any nDrotations in completely orthogonal planes: Theorem 4 (Decomposition of the nD-rotations). Any nD-rotation in planes completely orthogonal can be decomposed into standard rotations. The decomposition is not unique.

Let us introduce the notion of simple, double and k-rotations used in the proof of Theorem 4.

Definition 10 (Simple and double rotations).

A nD-rotation with only one (resp. two) angle is called simple (resp. double) rotation.

Definition 11 (k-rotation).

A nD-rotation which can be decomposed into k planar rotations with k angles (3

≤ k ≤ n 2 ) is called by k-rotations.
More precisely, Theorem 4 can be reformulated as : a simple rotation can be decomposed into two standard rotations, a k-rotation can be decomposed into k standard rotations. Let us prove Theorem 4.

Proof. In the proof, the notation R P (α) denotes the rotation of angle α in the plane P . The proof is in two parts: for k = 1 (simple rotation) and for k ≥ 2 (k-rotation).

* Simple rotations:

Let R be the nD-rotation of an angle θ in the rotation plane P . The rotation R is the composition of the rotation R 1 of angle θ in the plane P and the rotation R 2 of angle 0 in P with P and P completely orthogonal. The rotation R 1 can be written as:

R P (θ) = R P θ + 0 2 R P θ -0 2 .
In the same way, the rotation R 2 can be written as:

R P (0) = R P θ + 0 2 R P 0 -θ 2 .
With these decompositions, the rotation R can be written:

R = R P (θ)R P (0) = R P θ + 0 2 R P θ -0 2 R P θ + 0 2 R P 0 -θ 2 .
Since the rotation planes are completely orthogonal, the factors commute:

R = R P θ + 0 2 R P θ + 0 2 R P θ -0 2 R P 0 -θ 2 .
The rotation R is the composition of two standard double rotations of angle of magnitude |θ|/2 and | -θ|/2.

* the k-rotations:

Let R be a k-rotation with k angles θ 1 , θ 2 , ..., θ k . The rotation planes are supposed to be completely orthogonal. In each plane, the rotation can be expressed as:

R P1 (θ 1 ) = R P1 θ1+θ2 2 ...R P1 θ1+θ k 2 R P1 θ 1 -(θ1+θ2) 2 -... -(θ1+θ k ) 2 R P2 (θ 1 ) = R P2 θ1+θ2 2 R P2 -(θ1+θ3) 2 ...R P2 -(θ1+θ k ) 2 R P2 -θ 1 + (θ1+θ2) 2 + ... + (θ1+θ k ) 2 and for j = 2, ..., k -1 R Pj+1 (θ j+1 ) = R Pj+1 -(θ1+θ2) 2 ...R Pj+1 -(θ1+θj-1) 2 R Pj+1 θ1+θj 2 R Pj+1 -(θ1+θj+1) 2 ...R Pj+1 -(θ1+θ k ) 2 R Pj+1 -θ 1 + (θ1+θ2) 2 + ... + (θ1+θ k ) 2
If we denote by Θ 1 , Θ 2 , ..., Θ k the k angles of the standard double rotations then the angles are related by the following relations:

Θ i = θ1+θi 2 i = 2, ..., k. Θ k = θ 1 - k-1 i=1 Θ i
The first (k -1) th standard rotations are of angle of magnitude . The common factor is θ 1 but we can choose another one so the decomposition is not unique. Table 1 sums up all the different nD-rotations according to these criteria. The algorithms which permit their decomposition are also provided. They are developed in the following sections.

Rotations

Extension of the algorithm from [7]

In this section, for each pseudo-isoclinic rotations, a method to extend the algorithm provided in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] is proposed. A generalized method available for all pseudo-isoclinic is then provided.

First, we extend our algorithm to the isoclinic rotations, the pure pseudoisoclinic rotations and the ordinary pseudo-isoclinic rotations only in 6D with two equal angles. These cases are particular cases of the generalized method which deals with all the pseudo-isoclinic rotations provided in Section 3.3.

Isoclinic rotations

The first improvement of our algorithm presented in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] is the extension of the isoclinic rotations.

This particular case can be identified from the SVD. If the rotation is isoclinic then the eigenvalue is unique and is of order n. Thus, isoclinic rotations can be easily detected (cf Remark 2) and treated separately. In case of noisy data, the eigenvalues are compared using an error margin (|v

1 -v 2 | ≤ ⇒ v 1 = v 2 ).
In this case, all the half-lines from the origin point are rotated by the same angle and all planes generated by a vector and its image are invariant. In other words, a 4D-isoclinic rotation is completely characterized by any point of the hypersphere centered in 0 and its image which belongs to this hypersphere too [START_REF] Mebius | Applications of quaternions by dynamical simulation, computer graphics and biomechanics[END_REF]. In 4D, it is sufficient to find two couples (x i , y i ) and (x j , y j ) with y i = R(x i ) and y j = R(x j ) such that P i = x i ∧ y i and P j = x j ∧ y j are orthogonal. The planes P i and P j can be used as rotation planes. In a general way, two planes generated by two vectors and their images are not mutually orthogonal. Consequently, it is not judicious to randomly choose vectors among our input data. In order to correctly build those planes using only our input data (n points and their images) we apply the following method (based on the Gram-Schmidt orthogonalization). The first plane P 1 is arbitrarily defined by x 1 ∧y 1 . Then a couple (x i , y i ) which does not belongs to P 1 can be chosen. The projection of (x i ,y i ) on P 1 is given by (x i , y i ). The vectors x i⊥ = x i -x i and y i⊥ = y i -y i can be computed. They define a plane P 2 generated by a vector and its image 4 . Moreover the plane P 2 is orthogonal to P 1 (by construction). So we can choose P 1 and P 2 as rotation planes.

In nD, n 2 planes are required. In the same way, P 1 and P 2 are built. The third plane is generated by x k -proj P1 (x k ) -proj P2 (x k ) and y k -proj P1 (y k )proj P2 (y k ) with (x k , y k ) which belongs to neither P 1 nor P 2 . In a general way, the (i+1) th plane is generated by x i -i-1 j=1 proj Pj (x i ) and y i -i-1 j=1 proj Pj (y i ) with (x i , y i ) not belongs to P 1 , ..., P i . This method can be used in any even dimension space to characterize the isoclinic rotations.

Algorithm 2 Find P i (rotation plane) for isoclinic rotation in all even dimension space Input: (x i ) i=1,...,n , (y i ) i=1...n Output:

P i P 1 = x 1 ∧ y 1 for i = 2 to dim/2 do Choose (x i , y i ) / ∈ P 1 , ..., P i-1 x i = i-1 j=1 proj Pj (x i ) y i = i-1 j=1 proj Pj (y i ) P i = (x i -x i ) ∧ (y i -y i ) end for
The algorithm of [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] can be extended to the isoclinic rotations. It is sufficient to detect these rotations from the SVD and then to apply Algorithm 2.

Pure pseudo-isoclinic rotations

A nD-rotation R (with n odd) is said to be a pure isoclinic rotation if its n 2 angles are equal. A pure pseudo-isoclinic rotation has an invariant subspace F of dimension 1 and the restriction of R to P (of even dimension n -1) is an isoclinic rotation. Therefore, the subspace P is R-isogonal.

Let the matrix A be the rotation matrix. Let the matrix U be the matrix which comes from the SVD of the matrix (A-Id). The subspace P is generated by the (n -1) first column vectors of the matrix U and the subspace F is generated by its last vector.

One way to determine the rotation planes is to project all the n vectors and their images onto the subspace P using orthogonal contractions. Then, Algorithm 2 is applied. Naturally, its input data are these projections. This method can be extended to all the pure pseudo-isoclinic rotations in any dimension.

The rotation parameter determination for pure pseudo-isoclinic rotation in any dimension space is summarized in Algorithm 3.

Algorithm 3 Find P i (rotation plane) for pure pseudo-isoclinic rotation in all odd dimension space Input:

(x i ) i=1...n , (y i ) i=1...n Output: P i z i = y i -x i (S, U, V ) = SV D[z 1 ...z n ] P = n-1 i=1 U (:, i) for i = 1 to n do x i = proj P (x i ) y i = proj P (y i ) end for Choose (x k , y k ) = (0, 0) P 1 = x k ∧ y k for i = 2 to (dim -1)/2 do Choose (x i , y i ) / ∈ P 1 ...P i-1 x i = i-1 j=1 proj Pj (x i ) y i = i-1 j=1 proj Pj (y i ) P i = (x i -x i ) ∧ (y i -y i ) end for

Ordinary pseudo-isoclinic rotations in 6D with exactly 2 equal angles

To better understand the following, we start to detail the 6D-rotation with three angles θ 1 , θ 2 and θ 3 where exactly two are equal. In the following, the angles θ 1 and θ 2 are supposed equal. 5 The rotation R has got two pseudoisogonal components. This is an example which can be generalized to any dimension and any pseudo-isogonal component number (see Section 3.3).

The space E can be decomposed as E = P ⊕ F . The space F is reduced to 0. The space P can be decomposed as P = P 1 ⊕ P 2 ⊕ P 3 . We suppose that the rotation angle θ i belongs to P i , for i = 1, 2, 3.

Let A be the rotation matrix. The SVD of the matrix (A -Id) gives two distinct singular values: s 1 of order 4 and s 2 is of order 2. Let (j, j + 1) with j odd be the indexes of the matrix S such that S(j, j) = S(j + 1, j + 1) = s 2 . The j th and (j + 1) th columns of the matrix U (denoted by U (:, j) and U (:, j + 1) in the following) span P 3 and the four other columns span P 1 ⊕ P 2 .

In other words, P 3 = U (:, j) ∧ U (:, j + 1).

P 1 ⊕ P 2 = 6 k=1 U (:, k), k = j, j + 1. (1)
The restriction of the rotation R to P 1 ⊕ P 2 is an isoclinic rotation, that means the space P 1 ⊕ P 2 is R-isogonal. In this space all planes generated by a vector and its image are invariant. Thus, the spaces P 1 and P 2 are generated by a vector and its images.

If the six vectors and their images are projected onto the space P 1 ⊕P 2 , these projections can be used as input data (instead of the x i and y i ) of Algorithm 2.

The algorithm provided in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] can be extended using these three previous methods as shown by Algorithm 4. Nevertheless, a generalized method can be now proposed which includes the previous three algorithms. This method works with all the pseudo-isoclinic rotations.

General method

The aim of this section is to provide an algorithm which can be applied to the three pseudo-isoclinic rotations previously stated in any dimension space. It is a generalization of the three algorithms presented in Sections 3.2.1, 3.2.2, 3.2.3. In particular, it is based on the identification of the R-isogonal subspaces of the nD-rotation R. Actually, any pseudo-isoclinic rotations have got at least one R-isogonal subspace.

Let R be a nD-rotation in E and let A be its matrix in the canonical basis of E. The space E can be decomposed into E = P ⊕ F where F is the invariant subspace and P its orthogonal. The space F is eventually reduced to 0.

Let us introduce some notations: * Let θ 1 , ..., θ k for 1 ≤ k ≤ n 2 be the k rotation angles of R. * Let j be the number of distinct pseudo-isogonal components of the partition.

* Let C αi i for 1 ≤ i ≤ j be the i-th component of the partition with length α i .

By convention, since singular values are in decreasing order, the components are in angle value decreasing order.

Since, the pseudo-isogonal components form a partition of the angle set, we have

j i=1 α i = k.
Consequently, the space P can be decomposed into j subspaces as

P = P 1 ⊕ ... ⊕ P j where each P i is of dimension 2α i with 1 ≤ α i ≤ n 2 .
This decomposition is called R-isogonal subspace decomposition of the rotation. 6Obviously, if α i ≥ 2, the subspace P i is R-isogonal. These subspaces can be calculated from the SVD of the matrix A -Id. Actually, the matrix A can be decomposed into U SV t . The matrix S is diagonal and its diagonal elements s i for i = 1, ..., n are in decreasing order. The components are in decreasing order too. Thus, the first (2α 1 ) th singular values are equal, the next (2α 2 ) th are equal too, and so on. Consequently the first (2α 1 ) th columns of the matrix U span P 1 , the next (2α 2 ) th span P 2 and so on. In a general way, we have the following equality:

P k = 2α k i=2α k-1 +1 U (:, i). (2)
By convention, α 0 = 0. The formula (2) extends the formula (1) given in the previous part. Using the formula (2), rotation planes are easy to retrieve. Actually, for each P k , two cases are possible:

* Either 2α k = 2 then P k is a plane. It is determined by the formula (2).

* Or 2α k > 2 then P k is not a plane but a direct sum of planes. All the n vectors and their images are projected on R-isogonal subspace P k . The planes are determined by application of Algorithm 2.

This can be summarized by Algorithm 5.

Algorithm 4 Find P i (rotation plane) in the pseudo-isoclinic rotation case.

Input: (x i ) i=1...n , (y i ) i=1...n Output:

P i z i = y i -x i (S, U, V ) = SV D[z 1 . . . z n ]
Calculate the pseudo-isogonal partition composed with the C αi i for 1 ≤ i ≤ j (i-th component which length α i of the partition, in conventional order).

P = P 1 ⊕ ... ⊕ P j P k = 2α k i=2α k-1 +1 U (:, i) with α 0 = 0 for all P k do if 2α k = 2 then P k = 2α k i=2α k-1 +1 U (:, i) else
Calculate the projections into P k of (x i ) i=1...n , (y i ) i=1...n P k =Output of Algorithm 2 with these projections as input data end if end for

An alternative: the Schur decomposition

It is well-known that every rotation matrix can be block-diagonalized in an orthogonal basis [START_REF] Audin | Geometry[END_REF]. One way to block-diagonalize a matrix is to compute Algorithm 5 General algorithm to find P i (rotation plane) and θ i (rotation angle). The lines 3-9 can be replaced by the algorithm 4. Input: (x i ) i=1...n , (y i ) i=1...n Output: P i , θ i // Compute rotation planes 1: Find x j = x j , P i > ⊥ , P i > ⊥ and y j = y j , P i > ⊥ , P i > ⊥ non-null 18:

z i = y i -x i 2: (S, U, V ) = SV D[z 1 . . .
Calculate cos(θ j ) =

x j .y j ||x j ||||y j || and sin(θ j ) =

x j ∧y j ||x j ||||y j ||Pi . 19: end for the Schur decomposition [START_REF] Golub | Matrix computations[END_REF]15]. This decomposition is very used for example for the resolution of Riccati equations [START_REF] Laub | A Schur method for solving algebraic Riccati equations[END_REF], for the definition of the matrix sign fonction [START_REF] Sun | Perturbation analysis of the matrix sign function[END_REF] and in control theory [START_REF] Bojanczyk | The periodic Schur decomposition. algorithms and applications[END_REF]. However, according to the best knowledge of the authors, the Schur decomposition has never been used to characterize rotations that means to find rotation angles and rotation planes. In the following, we explain how to do this.

The principle of the Schur decomposition is to factorize a given matrix A into A = U T U t where U is orthogonal and T is block-diagonal with blocks of size 1 × 1 or of size 2 × 2. In order to become block-diagonal, a rotation matrix A undergoes orthogonal transformations. These transformations are stored into the matrix U (which is a basis change matrix). The blocks 2 × 2 of the matrix T represent 2 × 2 planar rotation matrices and the blocks 1 × 1 represent the identity matrices of order 1. From T , the rotation angles are thus retrieved and from the corresponding columns of U the rotation planes can be computed. This method is summarized by Algorithm 6.

Algorithm 6 Find P i (rotation plane) and θ i (rotation angle) Input: (x i ) i=1...n , (y i ) i=1...n Output: P i , θ i A = [y 1 ...y n ] (U, T ) = Schur-decomposition(A) // U is orthogonal, T is block-diagonal for each 2 × 2-block in T do (j, j + 1) = indexes of columns which compose the block θ i = acos(T jj ) P i = T (:, j) ∧ T (:, j + 1) end for

Comparison of the different algorithms

Numerical experimentations

Our algorithm, its improvement and the Schur decomposition have been implemented and compared in higher dimensions with the library developed in [START_REF] Charneau | Étude et application des algèbres géométriques pour le calcul de la visibilité globale dans un espace projectif de dimension n ≥ 2[END_REF] using the OCaml language [25].

The following numerical experiments have been conducted. The rotation planes and angles are randomly chosen and the rotation matrices A i are then generated. To simulate the noise, the obtained matrices are biased that means the second or the third decimal of each floating number value is randomly changed. The biased matrices are denoted by B i .

Afterward, the different methods are applied to these matrices. All the methods are experimented with the same matrices. The estimated rotation planes and angles are thus estimated (resp. P i and θ i ). In [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], in order to estimate the error of the method, we successively compute the rotor R corresponding to these parameters and the image of each basis vector x k using Rx k R -1 . The rotor R is given by R = cos(θ i ) -P i sin(θ i ). The matrix such that the k th column is the vector Rx k R -1 is denoted by E i and is compared with the initial matrix A i . In [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], we have (arbitrarily) considered that the plane P i is generated by the (2i -1) th and the 2i th columns of the matrix U which results from the SVD. Since there exists k!

2 k 2 k 2 !
possibility of permutations of the vectors of U (where k is the number of columns of U ), it is not obvious to retrieve the same matrix. Here, we evaluate the transformation rather than to calculate the images of basis vectors. To do so, the rotor R is calculated, then a vector v is randomly generated. Its image v theo is calculated using the matrix of y i and is compared to v prat = RvR -1 . The error is thus equal to the norm of (v theo -v prat ). As in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], the tests have been conducted in 4D, 5D and 6D (not in 3D because there exists only one rotation plane so the change of evaluation do not modify the result) with 500 matrices in each case. Dim.

Prec.

Algo [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] ( (2) algorithm presented in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] with the new evaluation; (3) adaptation to the isoclinic case with new evaluation; (4) the Schur decomposition, (5) adaptation to all the pseudo-isoclinic cases.

Result analysis

The results are given in Tables 2,3, 4 and 5. As in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], experiments have been conducted with a precision 10 -2 and 10 -3 for each dimension (4D, 5D and 6D).

Table 2 gives for each dimension and each precision the mean error for each algorithm. The quantities have no unit since it is a norm. The error is calculated for the algorithm presented in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] with the previous (1) and the new evaluation (2), its extension to the isoclinic rotations only (3)(no result in odd dimension), its extension to all the pseudo-isoclinic rotations (5) and for the algorithm deduced of the Schur decomposition [START_REF] Gower | Procrustes problems[END_REF].

The comparative study between all the algorithms is given by Table 3. Percentages indicate the improvements between the various algorithms. The mean improvement is given in the last row.

The error is smaller using the new evaluation than the previous one. Since parameter rotations are the same in both cases, this proves that our evaluation used in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] is not a good evaluation. The rotation parameters are well estimated but due to their non-unicity, the error is very high. The evaluation method used in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] must therefore be used only if the parameters are unique. [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] using the new evaluation) and after improvement (in isoclinic and pseudo-isoclinic case), using algorithm 5 for all the matrices and only those corresponding to an isoclinic or pseudo-isoclinic rotations.

The algorithm extended to the isoclinic case have an edge (7%). This can be explained by the number of isoclinic rotations. These numbers are given in the table [START_REF] Gower | Procrustes problems[END_REF]. In this case, the mean over all the 500 matrices is not representative when only 2 matrices are concerned as for example in 6D. A better comparison in thus provided by the table 5 where the mean over the 500 matrices and only matrices which correspond to isoclinic rotations and pseudo-isoclinic rotations is calculated. By considering only these matrices, the mean is improved by almost 69%. In a general way, our generalized algorithm is very efficient when data represent isoclinic or pseudo-isoclinic rotations. The improvement between the initial version proposed in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] and the extended version is very significant. The improvement is almost 63%. Moreover, our new algorithm is better than those deduced of the Schur decomposition (+50%).

Numerical example in 3D

rotation is given by the angle between x 1 and y 1 , the projection of x 1 and y 1 in P . For more numerical examples in particular in 4D, the reader can refer to [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF].

Conclusion

In [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF], we have proposed an algorithm that computes the decomposition into planar rotations of an nD rotation. The input data are n points and their images. Using the geometric algebra framework this algorithm is easy to implement.

If the data do not represent a pseudo-isoclinic rotation, the calculated planes and angles are very close to the exact ones even if the input data are biased. However, the algorithm was not adapted to the pseudo-isoclinic rotations.

In this paper, we have proposed an improvement which takes into account this particular case of rotations. The evaluation method is changed. The new method is more adapted since the parameters are not unique.

Our generalized algorithm is very efficient when the data represent isoclinic or pseudo-isoclinic rotations. The improvement between the initial version proposed in [START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF] and the extended version is significant. Moreover, our new algorithm is better than those deduced of the Schur decomposition.

A classification of nD-rotations is also presented. Moreover, some properties for decomposition of these rotations are also provided.

As an application of this work, we plan to use it in order to describe nD discrete rotations by extending the hinge angle notion developed for 3D rotations in [START_REF] Thibault | Computing upper and lower bounds of rotation angles from digital images[END_REF].

Figure 2 :

 2 Figure 2: Contraction of a vector x by a plane B.

θ1+θj 2 for j = 2 , 2 -

 222 ..., k. The last one is of angle of magnitude θ 1 -(θ1+θ2) ... -(θ1+θ k ) 2

Table 1 :

 1 The nD-rotation classification and the associated decomposition algorithms. The letter E (resp. O) denotes the even (resp. odd) dimension, Na, the number of angles and Nc the number of pseudo-isogonal components.

			N a	N c	Parity of n	Algo
	Pseudo-isoclinic	General Isoclinic Pure pseudo-isoclinic Ordinary pseudo-isoclinic	≤ n 2 n 2 n 2 x	N a (of length 1) 1 (of length n 2 ) 1 (of length n 2 ) y < x	E & O E O E & O	Algo of [7] Algo 2 Algo 3 Algo from section 3.2.3	Algo 4	Algo 5

  z n ] 3: match case with 4: case [ (n even) & (One singular value (contained in S) of order n ) ] // Isoclinic rotations 5: |P i = Output of Algorithm 2 6: case [(n odd)& (One non-null singular value (contained in S) of order n -1) ] //Pure pseudo-isoclinic rotations 7: |P i = Output of Algorithm 3 8: case [ at least one singular value of order at least 4 ] //Ordinary pseudo-isoclinic rotations 9: |P i =Output of the Algorithm of section 3.2.3 10: case: other cases //General case 11: {k = rank(U ) 12: for all i such that 1 ≤ i ≤ k 2 do

	13:	P i = U (:, 2i -1) ∧ U (:, 2i)
	14: end for}
	15: end match
		// Compute rotation angles
	16: for all P i do
	17:

Table 2

 2

	: Result of the various algorithms: (1): algorithm presented in [7] without modification;
	(2) algorithm presented in [7] with the new evaluation; (3) adaptation to the isoclinic case
	with new evaluation; (4) the Schur decomposition; (5) adaptation to all the pseudo-isoclinic
	cases.
	Dim. P rec. (1)vs(2) (2)vs(3) (2)vs(4) (2)vs(5) (4)vs(5)
	4D	0.001 0.01	52% 50%	15% 8%	43% 40%	15% 8%	52% 45%
	5D	0.001 0.01	58% 56%	--	47% 42%	24% 7%	59% 46%
	6D	0.001 0.1	59% 59%	0% 5%	32% 31%	32% 20%	54% 45%
	Mean	55.6%	7.0%	39.2%	17.7%	50.2%

Table 3 :

 3

Comparative study of the different algorithms: (1): algorithm presented in

[START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF]

without modification;

Table 4 :

 4 Number of isoclinic and pseudo-isoclinic rotations detected by the comparison of the singular values.

	Dimension Precision Isoclinic rotations Pseudo-isoclinic rotations
		4D	0.001 0.01	31 17		--
		5D	0.001 0.01	--		35 24
		6D	0.001 0.01	4 2		97 57
	Dim. Prec.	Before	All matrices After	Comp	Iso. & Pseudo-iso. cases Before After Comp
	4D	0.001 0.013 0.01 0.12	0.011 0.11	(15%) (8%)	0.054 0.37	0.018 0.10	(67%) (73%)
	5D	0.001 0.017 0.01 0.14	0.013 0.13	(24%) (7%)	0.085 0.40	0.026 0.09	(69%) (78%)
	6D	0.001 0.028 0.01 0.20	0.019 0.16	(32%) (20%)	0.078 0.35	0.032 0.12	(59%) (66%)

Table 5 :

 5 Mean of the error before (algorithm presented in

The (y i -x i ) (i=1...n) family generates P . So, the obtained matrix must be of rank

2r ≤ n that is at most 2r = n -1 if n is odd.

If A is an a-blade and B a b-blade then < A, B > ⊥ is the subspace of dimension (b -a) in B which is orthogonal to A

By convention, singular values are in decreasing order.

From[START_REF] Richard | An algorithm to decompose ndimensional rotations into planar rotations[END_REF]:y i⊥ + y i = R(x i⊥ ) + R(x i ) and y i⊥ = R(x i⊥ ) ∈ F, y i = R(x i ) ∈ P

In other words, θ 1 =θ 2 , θ 3 =θ 2 , θ 3 =θ 1 .

If for all i, α i = 1, then each P i are planes and P = P 1 ⊕ ... ⊕ P j