Eric Andres
email: eric.andres@univ-poitiers.fr

Mousumi Dutt
email: duttmousumi@gmail.com

Arindam Biswas
email: barindam@gmail.com

Gaelle Largeteau-Skapin

Rita Zrour

Digital Two-dimensional Bijective Reflection and Associated Rotation

Keywords: Digital Reflection, Bijective Digital Rotation, Digital Rotation, Bijective digital rigid motion Original image with Continuous Reflection Line Discretized Image Reflection Bijective Digital Reflection

In this paper, a new bijective reflection algorithm in two dimensions is proposed along with an associated rotation. The reflection line is defined by an arbitrary Euclidean point and a straight line passing through this point. The reflection line is digitized and the 2D space is paved by digital perpendicular (to the reflection line) straight lines. For each perpendicular line, digital points are reflected by central symmetry with respect to the reflection line. Two consecutive digital reflections are combined to define a digital bijective rotation about arbitrary center (i.e. bijective digital rigid motion).

Introduction

Reflection transformation is a one of the most basic linear transforms [START_REF] Goodman | Alice through looking glass after looking glass: The mathematics of mirrors and kaleidoscopes[END_REF] with many potential applications. In particular, they are key to defining n-dimensional rotations [START_REF] Richard | Decomposition of ndrotations: Classification, properties and algorithm[END_REF] and can be used in applications where reflections may play a key role such as lightwave simulations [START_REF] Heckbert | Beam tracing polygonal objects[END_REF][START_REF] Zrour | Parallelization of a discrete radiosity method[END_REF][START_REF] Mora | Pedagogic discrete visualization of electromagnetic waves[END_REF]. There is however surprisingly few works that deal with the problem such a transform in the digital world. The goal of this paper is to propose an algorithm that performs a digital bijective reflection on a two-dimensional image and to use this to propose a digital bijective rotation algorithm that works for an arbitrary rotation center (i.e. bijective digital rigid motion. See [START_REF] Pluta | Bijective digitized rigid motions on subsets of the plane[END_REF] for a first algorithm and a theoretical discussion on the subject of digital rigid motions). While the reflection transform in the digital world has not been studied much, digital rotations and more precisely bijective digital rotations has been a point of interest of the digital community for some years now [START_REF] Andres | The quasi-shear rotation[END_REF]. The characterization of a discretized rotation is stated in [START_REF] Nouvel | Characterization of bijective discretized rotations[END_REF] which is the composition of an Euclidean rotation with rounding operations. Using Gaussian integers the bijective discrete rotation is characterized in [START_REF] Roussillon | Characterization of bijective discretized rotations by gaussian integers[END_REF]. In [START_REF] Nouvel | Incremental and transitive discrete rotations[END_REF], an incremental approach to discretized rotations is explained. In [START_REF] Pluta | Bijective rigid motions of the 2d cartesian grid[END_REF], bijective rigid motion in 2D cartesian Corresponding Author grid is stated. Other mentionable work on characterization of bijective rotations in 3D space is [START_REF] Pluta | Bijectivity certification of 3d digitized rotations[END_REF]. The bijective rotation is also studied on hexagonal grid in discrete space [START_REF] Pluta | Characterization of bijective digitized rotations on the hexagonal grid[END_REF].

A continuous tranform applied to a digital image is not, in general, bijective, nor injective or surjective for that matter, even for an apparently very simple transform such as the reflection transform (see Figure 1). Our idea is to digitize the Euclidean reflection straight line as naive digital lines (i.e. 8-connected digital lines) and to partition the two-dimensional space into naive digital lines that are perpendicular to the reflection line. A central symmetry is performed on the points of each Perpendicular Digital Straight Line (PDSL) according to the Digital Reflection Straight Line (DRSL). This is always possible because we can compute the exact pixel to pixel correspondance on each side of the DRSL for such naive digital lines. The main problem comes from the fact that the digital reflection line and a given digital perpendicular line may or may not intersect and so two different cases have to be considered.

The organization of the paper is as follows: in Sec. 2, we present the preliminaries. Section 3 presents the mathematical foundation of bijective digital reflection. The bijective digital rotation about an arbitrary center can be determined by applying digital reflection twice (see Sec. [START_REF] Zrour | Parallelization of a discrete radiosity method[END_REF]). An error criteria based on the distance between the continuous and the digital rotated points are presented in Sec. 4.2. The concluding remarks are presented in Sec. 5.

Preliminaries

Let {i, j} denote the canonical basis of the 2-dimensional Euclidean vector space. Let Z 2 be the subset of R 2 that consists of all the integer coordinate points.

A digital (resp. Euclidean) point is an element of Z 2 (resp. R 2 ). Two digital points p (p x , p y ) and q (q x , q y ) are said to be 8-connected if |p x -q x | 1 and

|p y -q y | 1.
For x ∈ R, x is the biggest integer smaller or equal to x and x is the smallest integer greater or equal to x.

A naive digital straight line is defined as all the digital points verifying

-max(|a|,|b|) 2 ax -by < max(|a|,|b|) 2
, where a b represents the slope of the Digital Straight Line (DSL) and ω = max(|a|, |b|) the arithmetical thickness [START_REF] Reveillès | Calcul en nombres entiers et algorithmique[END_REF]. A naive DSL is 8-connected such that if you remove any point of the line then it is not 8-connected anymore [START_REF] Reveillès | Calcul en nombres entiers et algorithmique[END_REF]. There are no simple points.

A reflection transformation R θ,(xo,yo) : R 2 → R 2 reflects (or flips) a continuous point, like in a mirror, on a continuous straight line called the reflection line.

The reflection line is defined as the line passing through a point of coordinates In the center, the continuous reflection transform composed with a rounding function. On the right, our proposed bijective digital reflection

(x o , y o ) ∈ R
3 Digital Reflection

Principle

Let us consider the continuous reflection transform, R θ,(xo,yo) . It is not difficult to note that if we simply compose the continuous reflection with a digitization transform (such as for instance D : (x, y) → ( x + 0.5 , y + 0.5 )) then we have a transform that is, in general, neither injective, surjective and of course not bijective (see the middle image in Figure 1).

To avoid this problem and create a digital bijective reflection transform, a completely digital framework is proposed here based on the following digital primitives:

-A naive DSL called digital reflection straight line (DRSL);

and a partition of the digital space with naive DSLs that are perpendicular to the digital reflection straight line. These digital lines are called Perpendicular Digital Straight Lines (PDSL).

The idea of the reflection method is the following: let us consider a point p that belongs to a given digital perpendicular straight line P k (k is a parameter identifying the PDSL in the partition). The point p may belong to the digital reflection line (in case there is an intersection) or lie on either side of that line.

Let the digital reflection of p be determined in the following way (see Figure 2):

let us first suppose that p lies in the intersection between the PDSL P k and the DRSL. Such an intersection does not always exist between the two digital lines but when it does (in this case point p) then the image of p by the digital reflection is p; -Let us now suppose that p does not belong to the intersection between the PDSL and the DRSL (whether such an intersection exists or not). The point p is then located on one side of the DRSL: the digital space is split by the DRSL into three regions, the DRSL and two regions on each side of it. All the digital points of P k on each side of the DRSL can formally be ordered according to the distance to the DRSL with as first point the closest one to the DRSL. As we will see, there actually is no real ordering needed. Let us suppose that p is the n th point of the ordered list on its side and a digital point q is the n th point of the list on the other side of the DRSL. The digital reflection of p will then be the point q and vice-versa.

Figure 2 shows the two possible cases: on the left, the perpendicular digital line has an intersection point (the point marked as point 0) with the digital reflection line, while on the right there is no intersection point. The points marked by the number i (from 1 to 5) are images of each other. The image of the point marked 0 is itself. The idea is to create a digital transformation that is bijective and easily reversible. Let us now, in the next section, present the various definitions and mathematical details to make this work out.

Mathematical Details

The mathematical details of the works is stated here. Let us consider the continuous reflection L line defined by the continuous point (x o , y o ) ∈ R 2 and the direction vector v = (a, b) = (sin θ, cos θ). The analytical equation of the reflection line is given by: L

= (x, y) ∈ R 2 : a(x -x o ) -b(y -y o ) = 0 .
The digital reflection straight line (DRSL) the digitization (as naive DSL) of L defined as all the digital points verifying:

DRSL : - max(|a|, |b|) 2 a(x -x o ) -b(y -y o ) < max(|a|, |b|) 2
The perpendicular digital straight lines (PDSL) P k are naive DSLs defined as all the digital points (x, y) verifying :

P k : (2k -1) max (|a|, |b|) 2 b(x -x o ) + a(y -y o ) < (2k + 1) max (|a|, |b|) 2
It is easy to see that all the PDSLs are perpendicular to the DRSL and that the set of PDSLs, as a set of naive digital lines, partitions the two dimensional digital space.

Let

(2k-1)b 2 b(x -x o ) + a(y -y o ) < (2k+1)b 2 .
The method supposes that we are able to formally order the points of a given perpendicular line that are located on either side of the reflection line. This is always possible because the PDSLs are naive lines and for an angle verifying -π/4 θ π/4, there is one and only one point per integer ordinate y (For angles π/4 θ 3π/4, we will have, symmetrically, one and only one point per integer abscissa x in a given PDSL).

More precisely,

(2k-1)b 2 b(x -x o ) + a(y -y o ) < (2k+1)b 2 means that 2k-1 2 + x o -a b (y -y o ) x < 2k+1 2 + x o -a b (y -y o ). Since 2k+1 2 -2k-1 2
= 1, there is one and only one value x. For a given ordinate y, the abscissa x in the PDSL P k is given by the following function:

X (y) = (2k -1) 2 + x o -(a/b)(y -y o )
Let us note that for a given digital point p(x, y), it is easy to determine the PDSL P k it belongs to:

(2k-1)b 2 b(x -x o ) + a(y -y o ) < (2k+1)b 2 leads to 2k 2 b(x-xo)+a(y-yo) b + 1 < 2k + 2 and thus: k = (x -x o ) + a b (y -y o ) + 1 2
The next question that arises is the determination and the localization of the potential intersection point between the digital reflection line and a given digital perpendicular line. As already mentioned, there may be 0 or 1 digital intersection points between the DRSL and a given PDSL since both are naive digital lines. We propose here a simple criterion to determine the existence of given by (X (2y 1 -y p ), 2y 1 -y p ). Let us note here that the central symmetry

Algorithm 1: Reflection Transform R θ,(xo,yo) Input : (x, y) ∈ Z 2 , (xo, yo) ∈ R 2 , -π/4 θ < π/4 Output: (x , y ) ∈ Z 2 1 k = (x -xo) + a b (y -yo) + 1 2 2 Function X (y) : Z → Z : X (y) = (2k-1) 2 + xo -(a/b)(y -yo) 3 (x1, y1) = (X ( abk + yo ) , abk + yo ) 4 (x2, y2) = (X ( abk + yo ) , abk + yo ) 5 If -b/2 a(x1 -xo) -b(y1 -yo) < b/2 Then 6 (x , y ) = (X (2y1 -y), 2y1 -y) 7 Elsif -b/2 a(x2 -xo) -b(y2 -y -o) < b/2 Then 8 (x , y ) = (X (2y2 -y), 2y2 -y) 9 Else (x , y ) = (X (y1 + y2 -y), y1 + y2 -y) 10 return (x , y )
around the ordinate y 1 does not mean here that there is a central symmetry around the point (x 1 , y 1 ) because space partition would not be guaranteed and thus bijectivity would be lost. The abcissa X (2y 1 -y p ) is computed so that the reflected point still belongs to P k .

-Else if (x 2 , y 2 ) belongs to the the DRSL then there is an intersection point, (x 2 , y 2 ). A central symmetry of the PDSL points can be performed around the ordinate value y 2 . The digital reflection of a point p (x p , y p ) ∈ P k is given by (X (2y 2 -y p ), 2y 2 -y p ).

-Otherwise there is no intersection point and (x 1 , y 1 ) and (x 2 , y 2 ) are the first points of the PDSL on each side of the Reflection line. The central symmetry can be performed around the ordinate value y1+y2

2 . The digital reflection of a point p (x p , y p ) ∈ P k is given by (X (y 1 + y 2 -y p ), y 1 + y 2 -y p ).

The important point here is that, since there is only one point per ordinate y, a central symmetry on the ordinate leads directly to the reflection point.

Algorithm 1 presents the digital reflection transform method.

Bijectivity of the Digital Reflection Transform

Algorithm 1 provides a digital reflection method for a digital point. It is not very difficult to see that this defines a bijective digital reflection transform. Let us briefly summarizes the arguments for this:

-The computation of the value k in line 1 of algorithm 1 yields the same value for any point of a given PDSL (and only for those).

-Line 6, 8 and 9 ensure that the image of a digital point of a given PDSL P k At first, we wanted to have an idea on the effect of the rotation center on the rotation error measure with the angle α = 0 and θ = π/8. This means that the first reflection corresponds to a reflection with the continuous reflection line y = y o , for -0.5 y o < 0.5 and the digital reflection y = 0. In Figure 4 we represented the influence of moving the center on the interval [0, 0.5] 2 . As can be seen, the error measures are almost not affected by the shift on the ordinate but greatly by a shift on the abscissa of the rotation center. This can be explained by the fact that the shift on the perpendicular lines has a direct influence on the lateral error that adds up to the global error. Let us note that the surfaces are similar looking for other couple of angles (α, θ) even though the amplitude of the error measure may be different.

Next, we wanted to have an idea on the effect of the angle α on the rotation error measures. So, for each rotation angle 0 θ π/4, the average error measures for 500 rotations with randomly chosen center and respectively angle α = 0 and a randomly chosen angle 0 α π/4. The result can be seen in figure 5. It can be seen that the angle α has only a minor influence on the average error but the maximal error is typically significantly increased for an angle α > 0. It is interesting to note that this is not always so clear: the digital rotation of angle π/6 and center (0, 0) has an average / maximum error of (0.6367, 1.3972) for α = 0 and (0.5726, 1.4186) for α = π/6. What can be noticed as well is that the error ratio between α = 0 and random chosen α is relatively stable over all the rotation angles. The error measure in general are higher than those obtained for Fig. 5. For each rotation angle 0 θ π/4, the average error measures of 500 Rotations with random center and respectively α = 0 and random α (with 0 α π/4) the shear based method [START_REF] Andres | The quasi-shear rotation[END_REF] that is also defined for all angles and centers. The reflection based rotation is however way easier to implement.

Conclusion

A novel bijective digital reflection transform in two dimensions has been proposed. The reflection line is defined by an arbitrary Euclidean point and a straight line passing through this point. A rotation can be defined by the composition of two reflections, and so an associated new bijective digital rotation transform has been proposed. This new bijective digital rotation is defined for all angles and for all rotation centers, defining a rigid motion transform (see figure 6 that illustrates the rotation of an image at various angles and centers). In average the distance between a continuously and a digitally rotated point is about 0.8 (for a pixel of side 1). This reflection based rotation seems well suited for extensions to higher dimensions which is not easily done with previous methods. This is the main perspective for the future: exploring digital reflection based bijective rotations in three and higher dimensions.

Fig. 1 .

 1 Fig. 1. Reflection Transform applied to an image. On the left, the original image. In the center, the continuous reflection transform composed with a rounding function. On the right, our proposed bijective digital reflection

Fig. 2 .

 2 Fig. 2. The image of each point is the mirror image in the perpendicular line according to the digital reflection straight line. The image of a point numbered i is the opposite point numbered i. The point numbered 0 is its own image.

2 a

 2 us now suppose for what follows, w.l.o.g, that we have -π/4 θ π/4, then 0 |a| b and thus max (|a|, |b|) = b. The DRSL is then defined by b (x -x o ) -b(y -y o ) < b 2 and the PDSLs are then defined by P k :

Fig. 3 .

 3 Fig.3. The green point represents the intersection point between the continuous reflection line and the continuous perpendicular line. The orange disks mark the two candidate points for the digital intersection. On the left (x 1 , y 1 ) is the digital intersection point between the DRSL and the PDSL. On the right there is no digital intersection point.

-Fig. 4 .

 4 Fig. 4. Rotations of center (x, y), angle π/8 and α = 0

6 Fig. 6 .

 66 Fig.6. Rotation of the Lena image by various angles and randomly chosen center.

y = 2y 2 -y or y = y 1 +y 2 -y. Let us note first that (x, y) = (X (y), y). For all three cases, we have respectively y = 2y 1 -y = 2y -1 -(2y 1 -y) = y, y = 2y 2 -y = 2y 2 -(2y 2 -y) = y and y = y 1 +y 2 -y = y 1 +y 2 -(y 1 +y 2 -y) = y which proves the point.

Reflection Based Rotation

As mentioned in the preliminaries, a continuous rotation transform Rot θ,(xo,yo) of center (x o , y o ) and angle θ can be defined as the composition of two reflections based on two reflection lines passing through (x o , y o ) with an angle θ/2 between the two lines. In the same way, a digital rotation:

Let us note right away that, since the digital reflection transform is bijective, the digital rotation based on the reflection transforms will be bijective as well.

Furthermore, the inverse transform is easily defined. This last point may seem obvious but it is not because a digital transform is bijective that the inverse transform is easily computed.

Compared to previous digital rotations methods [START_REF] Andres | The quasi-shear rotation[END_REF][START_REF] Jacob | On discrete rotations[END_REF][START_REF] Pluta | Characterization of bijective digitized rotations on the hexagonal grid[END_REF][START_REF] Nouvel | Incremental and transitive discrete rotations[END_REF], there is an extra parameter that comes into play: the angle α.

Rotation Evaluation Criteria

In order to evaluate the "quality" of such a digital rotation, let us present some simple error measures [START_REF] Andres | The quasi-shear rotation[END_REF]. Each grid point has one and only one image through a bijective digital rotation but that does not mean that the digital rotation is a good approximation of the continuous one. The idea is to measure what error is made by using the digital rotation instead of the continuous one in terms of distance to the optimal position (to the continuously rotated point). Let us finish by noting that for the best known digital bijective rotations [START_REF] Jacob | On discrete rotations[END_REF][START_REF] Roussillon | Characterization of bijective discretized rotations by gaussian integers[END_REF][START_REF] Pluta | Bijective rigid motions of the 2d cartesian grid[END_REF] and for the angles where the digitized rotations are bijective (except for the trivial kπ/2 angles), we obtain the best maximum MD values of √ 2 2 ≈ 0.7 and average distance AD of ≈ 0.3. Note also that this is not obtained for all angles.


Eric Andres
email: eric.andres@univ-poitiers.fr

Shear based Bijective Digital Rotation in Hexagonal Grids

Keywords: Hexagonal grid, Digital Rotation, Bijective Digital Rotation, Shear Transform

In this paper, a new bijective digital rotation algorithm for the hexagonal grid is proposed. The method is based on an original decomposition of rotations into shear transforms. It works for any angle with an hexagonal centroid as rotation center and is easily invertible. The algorithm achieves an average distance between the digital rotated point and the continuous rotated point of about 0.42 (for 1.0 the distance between two neighboring hexagons).

Introduction

Applying a continuous geometric transform T : R 2 -→ R 2 to a 2D regular grid (defined by two vectors and an origin) is not trivial, although it is a fundamental problem in imaging. A transformed grid point through a continuous transform is usually not a grid point. One of the ways to deal with this, is to add a digitization transform D, classically the nearest grid point. The digitized transform D•T , has however no reasons, in general, to be injective nor surjective, even when the original continuous transform T is bijective. This leads to information loss. A digitized rotation may cause up to 17% information loss. For this reason, in imaging softwares, interpolation and sometimes oversampling are used to perform digital rotations [START_REF] Danielsson | High-accuracy rotation of images[END_REF], [START_REF] Sterling | Two-Dimensional Image Rotation in: Mathematical Visualization: Algorithms, Applications and Numerics[END_REF]). This is however not always possible when the image has low color resolution or when grid points carry information that are not meant to be interpolated. That is where bijective digital rotations are useful. In our case, we are interested in bijective digital rotations on an hexagonal grid.

Hexagonal grids have enjoyed renewed interest these last couple of years [START_REF] Middleton | Hexagonal Image Processing: A Practical Approach[END_REF] with, in particular applications in 3D printing filling with honeycomb structures [START_REF] Compton | 3d-printing of lightweight cellular composites[END_REF], [START_REF] Gibson | Additive manu-facturing technologies: 3D printing, rapid prototyping, and direct digital manufacturing[END_REF]), networks [START_REF] Kayaturan | Routing in hexagonal computer networks: How to present paths by bloom filters without false positives[END_REF], [START_REF] Stojmenovic | Honeycomb network[END_REF]), image processing (Grigoryan and Agaian (2016), [START_REF] Karkishchenko | Hexagonal images processing over finite eisenstein fields[END_REF], [START_REF] Mostafa | An edge detection method for hexagonal images[END_REF]), computer vision [START_REF] He | Hexagonal structure for intelligent vision[END_REF], [START_REF] Veni | Vision-based hexagonal image processing using hex-gabor[END_REF]), etc. Our motivation for lossless rotations in hexagonal grids came from recent work by K. Pluta andal. (Pluta et al. (2017, 2018)) and the perspectives it opened for 3D bijective rotations of honeycomb structures for 3D printing.

Transformations in hexagonal grids are still a largely open problem with few references (see [START_REF] Her | Geometric transformations on the hexagonal grid[END_REF] for a discussion on transforms on hexagonal grids). Recently K. Pluta et al. showed that, as for the square grid [START_REF] Jacob | On discrete rotations[END_REF], [START_REF] Nouvel | Characterization of bijective discretized rotations[END_REF]), there are angles for which the digitized rotation is bijective [START_REF] Pluta | Honeycomb geometry: Rigid motions on the hexagonal grid[END_REF][START_REF] Pluta | Characterization of bijective digitized rotations on the hexagonal grid[END_REF]) in the hexagonal grid. The problem is that these angles represent only a subset of all angles [START_REF] Pluta | Honeycomb geometry: Rigid motions on the hexagonal grid[END_REF][START_REF] Pluta | Characterization of bijective digitized rotations on the hexagonal grid[END_REF]). This is the limitation we propose to overcome in this paper.

In order to propose a bijective rotation in the hexagonal grid that works for all angles, we propose to revisit an idea the author of the present paper has already used for bijective rigid motions in the classical square grid: decomposing a rotation into a sequence of shear transforms [START_REF] Andres | The quasi-shear rotation[END_REF] developed out of an original idea by [START_REF] Reveillès | Calcul en Nombres Entiers et Algorithmique[END_REF]). A.W. Paeth used such a decomposition in order to propose a fast anti-aliasing method for image rotations [START_REF] Paeth | A fast algorithm for general raster rotation[END_REF]). The idea proposed by J-P. Reveilles is to use the same decomposition to push rows of pixels by an integer number of pixels, ensuring bijectivity and trivial reversibility. In this paper we propose to consider three directions (the three symmetry axes of an hexagon) for the shear transforms. With three such shear transforms in the hexagonal coordinate system, one for each hexagonal symmetry direction, a novel rotation decomposition is proposed. The idea behind the shear transforms is to push whole rows of hexagons successively in those three directions by integer numbers of hexagons, ensuring bijectivity. As for the shear based rotation proposed by A.W. Paeth, we have angles where the formula diverges. There is however a very simple way to overcome this problem. This leads to a digital rotation algorithm for all angles in the hexagonal grid, with the limitation that the center is a grid point (the centroid of an hexagon). A future work could consist in lifting this limitation and proposing a rigid motion transform (rotation with an arbitrary center).

The organization of the paper is as follows: in section two, we present the preliminaries. In particular, we present shear transforms and how shear transforms have been used in the classical square grid to define digital rotations. An error criteria based on the distance between the continuous and the digital rotated points are presented. In section three, we introduce our method of digital bijective rotation for the hexagonal grids. We conclude and present perspectives in the last section.

Preliminaries

Hexagonal grid

We are considering an hexagonal cell (centered on the grid point) with "pointy top" hexagons (two sides of the hexagons are parallel to the ordinate axis of the classical Euclidean coordinate system). It should be easy to transpose this work for other hexagon orientations or hexagon sizes. The hexagons are regular with a side length of 1 √ 3 . This means that between the centroids of two neighboring hexagons, there is a distance of 1. There are various ways of creating a coordinate system for hexagonal grids. We chose a simple 2D coordinate system where the first coordinate, the hexagonal hx-axis, is given by the Cartesian

{-1, -1} {-1, 0} {-1, 1} {0, -1} {0, 0} {0, 1} {1, -1} {1, 0} {1, 1} 1 hx -axis hy -axis (-x, +y) Direction A Direction B Direction C 1  3
Figure 1: coordinate system, hexagon size and A, B and C directions vector (1, 0) and the second coordinate, the hexagonal hyaxis, is defined by the Cartesian vector 1 2 , √ 3 2 [START_REF] Snyder | Coordinate system for hexagonal pixels[END_REF]). The coordinate transforms from hexagonal grid to the Cartesian grid and vice versa are given by the following transforms [START_REF] Snyder | Coordinate system for hexagonal pixels[END_REF]):

Cart2Hex : (x, y) → (x h , y h ) = x - y √ 3 , 2y √ 3
Hex2Cart : (x h , y h ) → (x, y) =       x h + y h 2 , y h √ 3 2      
We define three directions: direction A defined by vector (1, 0), direction B defined by vector (0, 1) and direction C defined by vector (-1, 1). These vectors are given in the hexagonal coordinate system (see Figure 1). A A-row is a set of hexagonal grid points with hexagonal coordinates (x, y) ∈ Z 2 such that they have all the same y value. All the hexagons of an A-row are aligned in the A-direction and are generated by hexagonal coordinate vector (1, 0) and one of its points. A B-row is a set of hexagonal grid points with hexagonal coordinates (x, y) ∈ Z 2 such that they have all the same x value. All the hexagons of a Brow are aligned in the B-direction and are generated by hexagonal coordinates vector (0, 1) and one of its points. A C-row is a set of hexagonal grid points with hexagonal coordinates (x, y) ∈ Z 2 such that they have all the same x + y value. All the hexagons of a C-row are aligned in the C-direction and are generated by hexagonal coordinates vector (-1, 1) and one of its points.

Shear Transforms

The idea behind the method proposed in this paper is to decompose a rotation transform into a sequence of shear transforms in the hexagonal coordinate system. A shear transform is a linear mapping that translates a point in a given direction by an vector proportional (by a shear factor) to the signed distance to a line parallel to that direction. A typical shear transform uses the axis-lines as shear lines. Shear transforms preserve areas. This is why it is quite natural to decompose isometries into sequences of shear transforms (as atomic transforms). As an example, the decomposition of a rotation into three shear transforms with the Cartesian axes as shear lines (an 'ULU' decomposition with '1's on the diagonals) used by A.W. Paeth [START_REF] Paeth | A fast algorithm for general raster rotation[END_REF]):

cos θ -sin θ sin θ cos θ = 1 -tan θ 2 0 1 1 0 sin θ 1 1 -tan θ 2 0 1
This rotation decomposition leads to a fast and simple antialised rotation that is still in use in some image libraries such as ImageMagick.

Bijective Digital Rotations in the square grid

In this paper we are interested in bijective digital rotations on hexagonal grids. Let us first see how this problem has been tackled in the classical square grid. Let us consider a digitized rotation DR(α): DR(θ) :

x y → x cos θ + y sin θ + 0.5 -x sin θ + y cos θ + 0.5 with u the biggest integer smaller or equal to u (i.e. floor function). A digitized rotation is not, in general, surjective or injective, except for some angles where the transform is actually bijective. For the interested reader, please refer to [START_REF] Jacob | On discrete rotations[END_REF], [START_REF] Nouvel | Characterization of bijective discretized rotations[END_REF], [START_REF] Roussillon | Characterization of bijective discretized rotations by Gaussian integers[END_REF]. The problem with the angles for which the digitized rotation is bijective is that they do not cover all angles. J.P. Reveilles took the same decomposition than A.W. Paeth to propose a bijective digital rotation that works for all angles [START_REF] Reveillès | Calcul en Nombres Entiers et Algorithmique[END_REF]). E. Andres improved on this idea and gave formulas for an improved bijective digital rotation and a bijective rigid motion in the classical square grid [START_REF] Andres | The quasi-shear rotation[END_REF]). The idea that leads to a bijective digital rotation is the following: the shear coefficients multiplied by x and y respectively, -x tan θ 2 and y sin θ, are approximated by their closest integer -x tan θ 2 + 0.5 and y sin θ + 0.5 . Moving a grid point in the x-axis direction and y-axis direction by an integer displacement is obviously a reversible operation. One could think that this leads to a very coarse approximation of the continuous rotation, but this is not the case. The way to measure this "approximation" has lead E. Andres to propose an error measure for digital bijective rotations [START_REF] Andres | Cercles Discrets et Rotations Discretes[END_REF][START_REF] Andres | The quasi-shear rotation[END_REF])) which we'll present in the next subsection.

Error Measure

Each grid point has one and only one image through a bijective digital rotation but that does not mean that the digital rotation is a good approximation of the continuous one. To measure how "wrong" we are by choosing the digital rotation over the continuous one, we are considering two distance criterii [START_REF] Andres | Cercles Discrets et Rotations Discretes[END_REF][START_REF] Andres | The quasi-shear rotation[END_REF]). Let us first denote G = ax g + by g |(a, b) ∈ Z 2 ⊂ R 2 a grid defined by the point (0, 0) and two vectors x g and y g . Let us denote R θ (p) the continuous rotation of center (0, 0) and angle θ of a grid point p ∈ G and R θ (p) its digital rotation of center (0, 0) (R θ (p) is a grid point in this case).

The Maximum Distance error criteria (MD) consists in computing max p∈G (d(R θ (p), R θ (p)). The Average Distance error criteria (AD) consists in computing avg p∈G (d(R θ (p), R θ (p)), where avg p∈G is the average distance over the grid.

Bijective Digital Rotation in the Hexagonal Grid

3.1. Rotation decomposition into shears in the hexagonal grid K. Pluta showed that there exists a subset of angles for which the digitized rotation of an hexagonal grid is bijective [START_REF] Pluta | Honeycomb geometry: Rigid motions on the hexagonal grid[END_REF][START_REF] Pluta | Characterization of bijective digitized rotations on the hexagonal grid[END_REF]). Our idea is to use shear transforms to define a digital rotation is the hexagonal grid that works for all angles by pushing grid points into specific directions. Since we are working in an hexagonal grid, we have three privileged directions to push grid points: direction A defined by vector (1, 0), direction B defined by vector (0, 1) and direction C defined by vector (-1, 1) (vectors expressed in the hexagonal coordinate system. See Figure 1). The shear transforms for the three directions, expressed in the hexagonal grid coordinate system, correspond to the following matrices:

matA hex = 1 a 0 1 ; matB hex = 1 0 b 1 matC hex = 1 -c -c c 1 + c
matA hex and matB hex are classical shear transforms with the axis as shear lines. For the third matrix, it is the same idea: We have a direction C defined by the vector (-1, 1) with an orthogonal direction (1, 1) (in the hexagonal coordinate system) . This leads to matrix matC hex .

Since the rotation matrix is simplest expressed in the classical Cartesian coordinate system, we switch back to the Cartesian coordinate system, which leads to:

matA cart = 1 2a √ 3 0 1 matB cart =        1 + b 2 -b 2 √ 3 b √ 3 2 1 -b 2        matC cart =        1 -c 2 -c 2 √ 3 c √ 3 2 1 + c 2       
Each of the matrices matX cart corresponds to matX hex expressed in the classical Cartesian coordinate system. The idea of our method is to push grid points in those three directions in order to approximate a rotation. Therefore, we solved the following equation:

matA cart .matB cart .matC cart = cos θ -sin θ sin θ cos θ
which has a unique solution for a, b, c:

          a b c           =               -1 + √ 3-2 sin θ √ 3 cos θ-sin θ 1 -cos θ + sin θ √ 3 1 - √ 3-2 sin θ √ 3 cos θ-sin θ              
It is easy to see that a = -c. Let us note that from here on, we are only going to work in the hexagonal coordinate system with matrices matA hex , matB hex and matC hex .

Algorithm 1: RotCBA(x, y, θ): Point Rot. CBA of center (0, 0) and angle θ

Input : (x, y) ∈ Z 2 , 0 ≤ θ < 2π/3 Output: (x 3 , y 3 ) ∈ Z 2 If θ = π 3 Then x 3 y 3 ← -y x + y 1 Else a b ←         -1 + √ 3-2 sin θ √ 3 cos θ-sin θ 1 -cos θ + sin θ √ 3         2 x 1 y 1 ← x * (1 + a) + y * a + 0.5 -x * a + y * (1 -a) + 0.5 3 x 2 y 2 ← x 1 y 1 + x 1 * b + 0.5 4 x 3 y 3 ← x 2 + y 2 * a + 0.

Dealing with the divergence and the inverse

The proposed solution for (a, b, c) is not universal: the denominator for a (and c),

√
3 cos θsin θ, is equal to 0 for θ equal to π/3 or 4π/3 (all angles are considered to be between 0 and 2π). As one can see in Figure 2, the zeros for π/3 and 4π/3 are not of the same nature: lim t→ π 3 a = 0.5 and lim t→ 4π 3 a = ±∞. Angle π/3 is a point singularity that can be easily dealt with since it is trivial to rotate an hexagonal image by such an angle bijectively: R π 3 (x, y) = (-y, x + y) is still an hexagonal grid point if (x, y) is one. This leads to the digital Rotation RotCBA (Algorithm 1) that rotates an hexagonal grid point for angles from 0 to 2π/3. Algorithm 2: RotCBANeg(x, y, θ): Inverse Rot. CBA of center (0, 0) and angle θ

Input : (x, y) ∈ Z 2 , 0 ≤ θ < 2π/3 Output: (x 3 , y 3 ) ∈ Z 2 If θ = π 3 Then x 3 y 3 ← x + y -x -Rotation -π/3 1 Else a b ←         1 - √ 3-2 sin θ √ 3 cos θ-sin θ -1 + cos θ -sin θ √ 3         2 x 1 y 1 ← x + y * a + 0.5 y 3 x 2 y 2 ← x 1 y 1 + x 1 * b + 0.5 4 x 3 y 3 ← x 2 * (1 + a) + y 2 * a + 0.5 -x 2 * a + y 2 * (1 -a) + 0.5 5 return (x 3 , y 3 ) 6
Let us see how we can handle the problems around 4π/3. In Figure 2, we can see that all the values around angle 4π/3 are affected by the divergence. Actually, starting at 2π/3 the average (AD) and maximal (MD) distance error values start to increase. The idea here is to consider only rotations for angles between 0 and 2π/3 and map all the other angles to this interval. This is not that difficult since, in the hexagonal grid, rotations by angles kπ/3 (with a grid point center) are naturally bijective so that for an angle θ, we can always decompose it as θ = kπ/3 + θ where 0 ≤ θ ≤ 2π/3.

There is however another problem that needs to be dealt with if we want an easy to use digital rotation algorithm: invertibility. Indeed, with Algorithm 1, RotCBA(x, y, θ) -1 is not equal to RotCBA(x, y, -θ). Having a bijective transform does not automatically mean that we have an easy way of performing the inverse transform. The approach we adopted to be sure that our rotation for angle θ is easily invertible is to perform all hexagon push operations in the reverse order for an angle -θ (or more precisely 2π -θ). For this, we introduce Algorithm 2 which is literally the inverse mapping of the mapping of algorithm 1, with shear factors -(a, b) and reversed shear transforms application order. With all this, we can now propose a rotation for an hexagonal point. The rotation of a point (x, y) by an angle 0 ≤ θ < 2π is performed in the following way (see Algorithm 3):

Algorithm 3: RotPoint(x, y, θ): Rotation of center (0, 0) and angle θ of point (x,y) in the hexagonal Grid.

Input : (x, y) ∈ Z 2 , θ ∈ R 2 Output: (x , y ) ∈ Z 2 θ = θ Mod 2π 1 If 0 ≤ θ ≤ 2π 3 Then (x , y ) ← RotCBA(x, y, θ) 2 Elsif 2π 3 < θ ≤ π Then 3 (x , y ) ← RotCBA -x -y, x, θ -2π 3 Elsif π < θ < 4π 3 Then 4 (x , y ) ← RotCBANeg x, y, 4π 3 -θ ; (x , y ) ← (y , -x -y ) Else (x , y ) ← RotCBANeg(x, y, 2π -θ) 5 Return(x , y ) 6 1.
For an angle 0 ≤ θ2 ≤ 2π/3: apply algorithm 1 to the point (x, y) with angle θ; 2. For an angle 2π/3 < θ ≤ π: apply algorithm 1 to the point (-xy, x) with angle θ -2π/3 (the point (-xy, x) is the rotation of (x, y) by 2π/3); 3. For an angle π < θ < 4π/3: apply algorithm 2 to the point (x, y) with angle 4π/3 -θ and then perform a rotation of angle -2π/3 on the resulting point (x , y ) which maps the point (x , y ) to (y , -xy ) (This ensures that operations are exactly performed in the reverse order of point 2 (angle 2π/3 < θ ≤ π)). 4. For an angle 4π/3 ≤ θ < 2π: apply algorithm 2 to the point (x, y) with an angle 2π -θ. This ensures that operations are exactly performed in the reverse order of point 1 (angle 0 ≤ θ ≤ 2π/3).

This leads to point rotation RotPoint (Algorithm 3) and to the final bijective rotation RotHexa (Algorithm 4). We'll explain why Rotation Algorithm 4 is bijective in the next subsection.

Let us quickly comment Figure 3 which illustrates the distance error criterii for the bijective digital rotation algorithm. In order to have a point of comparison, let us note that the best MD and AD error values are obtained when considering digitized rotations in the hexagonal grid (that are not necessarily bijective). Except for the angles kπ/3, we obtain mostly MD error values of unit hexagon to its center). As a reminder, the distance between two neighboring grid points is 1. In figure 3, we see that for our algorithm, the average distance (AD criteria) between a point rotated continuously and a point rotated discretely is very stable and close to 0.42 (compared to 0.35). The maximal distance (MD criteria) does not exceed 1.1 (compared to 0.57). This means basically that an hexagon lands usually in the immediate grid neighborhood of the spot where he should be. Let us note that we obtain slightly better error values than for the similar type of algorithm in the square grid [START_REF] Andres | The quasi-shear rotation[END_REF]).

Bijectivity

All has been done inorder to have a bijective transform RotHexa. Let us show that formally by examining the mappings of Algorithm 1. Let us define the mapping A(a) : (x, y) → (x + ay + 0.5 , y), a ∈ R 2 (line 5 in Algorithm 1). For a A-row, we have A(a)(S ) = S (y does not change). All hexagons on a given A-row are pushed in direction A by the same integer number of grid points.

In the same way, for the mapping B(b) : (x, y) → (x, y + bx + 0.5 , y), b ∈ R 2 (line 4 in Algorithm 1), if S is a B-row, we have B(b)(S ) = S (x does not change). All hexagons on a given B-row are pushed in direction B by the same integer number of grid points.

Let us now consider the mapping C(a) : (x, y) → ( x(1 + a) + ay + 0.5 , -ax + y(1a) + 0.5 ), a ∈ R 2 (line 3 in Algorithm 1). Let us show that for a C-row S , we have C(a)(S ) = S (i.e. that x + y does not change). Let us consider a point (x, y) ∈ S and (x , y ) = C(a)(x, y). We have x + y = x(1 + a) + ay + 0.5 + -ax + y(1a) + 0.5 = x + y + a(x + y) + 0.5 + -a(x + y) + 0.5 . Now x + 0.5 is the rounding function and it is easy to see that -u + 0.5 =u + 0.5 which proves that x + y = x + y and thus that C(a)(x, y) ∈ S . All hexagons on a given C-row are pushed in direction C by the same integer number of grid points. All this is of course verified as well for Algorithm 2. All the other mappings (rotations by kπ/3 angles) are trivially bijective.

This proves that the transform presented in Algorithm 4 is bijective (as sequence of bijective mappings). By design, we have made sure that the transform is also easily invertible such that, for ImageIn an image in the hexagonal grid, we have RotHexa(ImageIn, θ) -1 = RotHexa(ImageIn, -θ). Figure 4 shows how each matrix acts on the hexagons. In order to illustrate the action of each shear transform, we set a = b = c = 1. The rows containing the rotation center (0, 0) (marked by an ellipse) do not move. On the right of Figure 4, we show an actual rotation of the proposed image by an angle of π/8. Lastly, Figure 5 shows the result of the bijective rotation on a 'Lena' image of size 32x32 and 128x128 for different angles.

Conclusions

In this paper, we have proposed a bijective digital rotation method for hexagonal grids that works for all angles. For this rotation, the average distance, on an image, between a point rotated by a continuous rotation and the point rotated by the digital rotation is about 0.42 (with a distance of 1.0 between two neighboring hexagons), while the maximal distance is bounded by 1.1. This algorithm extends the direct mapping proposed by K. Pluta in the sense that it works for all angles. However, contrary to his work, we have only an algorithm that works for hexagonal grid points as center. Extending our method to arbitrary centers would be a very interesting extension. There are many other questions that are raised by this work. Using shears for directions that are not the usual axis directions can of course be applied to the regular square grid. How does that type of rotation measure up to the previously published bijective rotations [START_REF] Andres | The quasi-shear rotation[END_REF]? the triangular grid is a dual grid of the hexagonal grid however so may be there is something that can be done there. What about bijective rotations in arbitrary grids? One aspect that is of special for us is the extension to higher dimensions especially for honeycomb type grids with possible applications in 3D printing.

Figure 2 :

 2 Figure 2: Average (AD) and Maximal (MD) Distance error criteria for Rotation CBA for angles 0 to 2π

Algorithm 4 :Figure 3 :

 43 Figure 3: Average (AD) and Maximal (MD) Distance error criteria for the hexagonal grid Rotation Algorithm for angles 0 to 2π

ABBBCAC

  of angle 𝜋/8, center (0,0) A Row pushed by -1 Row pushed by +1 A Row pushed by +2 B Row pushed by -1 Row pushed by -2 Row pushed by -3 Row pushed by +3 B Row pushed by +2 B Row pushed by +1 C Row pushed by +2 C Row pushed by +1 C Row pushed by -2 Row pushed by -1 Row pushed by -2 Row pushed by -3 C Row pushed by +3

Figure 4 :

 4 Figure 4: shears in directions A,B,C with a=b=c=1. The rows containing (0, 0) (marked by an ellipse) do not move. On the right, rotation by π/8 and center (0, 0).

Figure 5 :

 5 Figure 5: Examples with a Lena image size 32x32 and 128x128 in a hexagonal grid.

≈ 0.57 and AD error values of ≈ 0.35 (which can be analytically calculated as the distance from a randomly point inside a


Eric Andres
email: eric.andres@univ-poitiers.fr

Gaëlle Largeteau-Skapin
email: gaelle.largeteau.skapin@univ-poitiers.fr

Rita Zrour
email: rita.zrour@univ-poitiers.fr

Shear based Bijective Digital Rotation in Triangular Grids

Keywords: Triangular grid, Rotation, Bijective Rotation, Shear Transforms

In this paper we are proposing a way to perform bijective digital rotations on a triangular cell grid. The method is based on a decomposition of a rotation into shear transforms and on a way to transform the original triangle centroids in order to have a regular grid. The rotation method works for any angle and achieves an average distance between the digital rotated point and the continuous rotated point of about 0.5 (with 1.0 the side of a triangle and a distance between two neighboring triangle centroids of 0.57).

Introduction

The aim of this paper is to propose an algorithm for performing a bijective digital rotation on triangular grids. On one hand, rotation transforms are one of the most basic and fundamental operations in fields like computer graphics, computer vision, pattern recognition, etc. On the other hand, triangular grids [START_REF] Freeman | Algorithm for generating a digital straight line on a triangular grid[END_REF], Nagy (2003a[START_REF] Nagy | Generalized triangular grids in digital geometry[END_REF]) have enjoyed a continuous interests over the years in various fields like geometric modelling [START_REF] Botsch | Geometric modeling based on triangle meshes[END_REF]), computer vision Cooper et al. (2005), [START_REF] Ogawa | Labeled point pattern matching by delaunay triangulation and maximal cliques[END_REF], image processing [START_REF] Sun | Triangle mesh-based edge detection and its application to surface segmentation and adaptive surface smoothing[END_REF], digital tomography [START_REF] Nagy | Memetic algorithms for reconstruction of binary images on triangular grids with 3 and 6 projections[END_REF], networks [START_REF] Koster | AGraphs and Algorithms in Communication Networks: Studies in Broadband, Optical, Wireless and Ad Hoc Networks[END_REF], etc. To the best of the author's knowledge, there are very few papers yet on triangular cell grid transforms in general and none on digital rotations for such grids.

Let us first note that there are two different types of triangular grids in the literature, based on a triangle plane tessellation: the triangular vertex grid defined by the triangle vertices and the triangular cell grid formed by the centroids of the triangles. It is the latter, the triangular cell grid, that we are considering since we are interested in triangular grid image transforms. In the classical triangular cell grid, each triangle is identified by its centroid. The problem with this is that there are two types of cells (odd and even triangles) and, more importantly, the centroids do not form a regular 2D grid (defined by a point and two vectors). One of the first steps of our method is therefore to propose a small modification of the triangular cell grid in order to create a new grid that is regular. Each triangle is now identified by a point that is slightly different from the classical centroid.

Different methods have been used to rotate an image. There exists angles for which digitized rotations (a continuous rotation followed by a rounding operator) are bijective in the square grid [START_REF] Jacob | On discrete rotations[END_REF], [START_REF] Nouvel | Characterization of bijective discretized rotations[END_REF], [START_REF] Roussillon | Characterization of bijective discretized rotations by Gaussian integers[END_REF]. Recently, [START_REF] Pluta | Characterization of bijective digitized rotations on the hexagonal grid[END_REF] showed that there also exists angles for which such a digitized rotation is bijective in the hexagonal grid. The problem is that the digitized rotations are not bijective for all angles. Paeth-Tanaka proposed an approach [START_REF] Paeth | A fast algorithm for general raster rotation[END_REF], [START_REF] Tanaka | A rotation method for raster image using skew transformation[END_REF], that consists in decomposing a rotation matrix into three shear matrices. They used this decomposition to perform antialiased image rotations (bijectivity was not considered). In this paper we are following the idea of Paeth-Tanaka by considering a rotation decomposition into shear matrices that are each combined with a rounding function. The idea is to push grid cells (triangles) rows (aligned cells) by an integer number of cell spots in a given direction. Pushing a cell row by an integer number of cells means that the whole process is reversible and that the resulting transform will be bijective. The idea of pushing rows of pixels has been originally introduced by J-P [START_REF] Reveillès | Calcul en Nombres Entiers et Algorithmique[END_REF], improved by [START_REF] Andres | The quasi-shear rotation[END_REF] for square grids and recently adapted to the hexagonal grid by [START_REF] Andres | Shear based bijective digital rotation in hexagonal grids[END_REF].

The rotation method works for any angle and achieves an average distance between the digital rotated point and the continuous rotated point of about 0.5 (with 1.0 the side of a triangle and a distance between two neighboring triangle centroids of 0.57).

The organization of the paper is as follows, Section 2 presents the basic notions and notations. Section 3 describes the rotation in the triangular grid. Section 4 provides examples as results and qualitative comparison with continuous euclidean rotation. The last section gives some conclusions and perspectives.

Triangular grid and Shear transforms

Description of the Triangular grid

We are considering a triangular cell grid with equilateral triangles of side one. One side of the cells is parallel to the Euclidean abscissa axis of the regular coordinate system. As shown in Figure 1, we consider a two coordinates system, where the tx-axis is the horizontal axis defined by the vector (1, 0) (in the triangular grid coordinate system) and the ty-axis the vertical axis defined by the vector (0, 1) (in the triangular grid coordinate system). The origin cell has the coordinates (0, 0) (in the triangular grid coordinate system). B. [START_REF] Nagy | Cellular topology on the triangular grid[END_REF][START_REF] Nagy | Cellular topology and topological coordinate systems on the hexagonal and on the triangular grids[END_REF] proposed a coordinate system with three coordinates instead of two. It is also the coordinate system we used in a previous paper on triangular grids [START_REF] Dutt | Characterization and generation of straight line segments on triangular cell grid[END_REF]. Let us note that our two coordinate system and Nagys' three coordinate system are linked by simple formulas (a%b is the Modulo of a by b):

• (x, y, z) → (x -z, -y); • (x, y) → x+y+(x+y)%2 2 , -y, -x+y+(x+y)%2 2 .
There are two types of triangles. The up-pointing triangles (yellow triangles in our figures), are called Even, the down-pointing triangles (blue ones), are called Odd. Using our triangular coordinate system, we can see that an Even (resp. Odd) triangle of coordinates (i, j) has an even (resp. odd) sum i + j.

Traditionally, for cell transforms, the centroid is considered as representative of the cell. If we consider the triangle centroids, we can see that Even and Odd triangles centroids do not form a regular grid defined by two vectors and a center (see Figure 1). The centroids of the triangles form two intertwined grids. One solution could be to rotate each grid separately but this poses the problem of the rotation center which cannot be a grid point for both grids. To avoid this problem, we propose to change the point that will represent the cell in the transform. The new triangular grid points are defined as follows:

• For even triangles, gridpoint = centroid + √ 3 12 , 0 ;

• for odd triangles, gridpoint = centroid -√ 3 12 , 0 .

In Figure 1, the red circles show the triangle centroids while the red disks show the grid points we are considering. As one can see, the centroids do not form a regular grid while the new grid points form a (rectangular) regular grid defined by the center (0, 0) and the Euclidean vectors 1 2 , 0 and 0, √ 3

2 . The coordinate transform from our triangular grid points to cartesian grid and vice versa are given by the following transforms :

T ri2Cart : (x t , y t ) → (x, y) = x t 2 , √ 3y t 2
(1)

Cart2T ri : (x, y) → (x t , y t ) = 2x, 2y √ 3 3 (2)

Shear Transforms for the Triangular Grid

The idea is to propose a rotation decomposition into a sequence of shear transforms in the (new) triangular coordinate system. A shear transform is a linear mapping that translates a point in a given direction by a vector proportional (by a shear factor) to the signed distance to a line parallel to that direction. Shear transforms preserve areas and because of this, it is quite natural to decompose isometries into sequences of shear transforms. As an example, the decomposition of a rotation into three shear transforms with the Cartesian axes as shear lines (an 'ULU' decomposition with '1's on the diagonals) used by [START_REF] Paeth | A fast algorithm for general raster rotation[END_REF] and [START_REF] Tanaka | A rotation method for raster image using skew transformation[END_REF]:

cos θ -sin θ sin θ cos θ = 1 -tan θ 2 0 1 1 0 sin θ 1 1 -tan θ 2 0 1
Our goal is to propose a bijective rotation and for that, following the original idea of [START_REF] Reveillès | Calcul en Nombres Entiers et Algorithmique[END_REF], we are going to push rows of triangular cells. The choice we made in this paper is to consider shear transforms corresponding to the directions of angle 0, π/3 and 2π/3. We have called those directions, for the sake of simplicity, directions A, B and C (See Figure 1). Direction A corresponds to a translation along the tx-axis: vector (1, 0). Direction B corresponds to a translation of vector (1, 1) in the tri-angular coordinate system. Direction C corresponds to a translation of vector (-1, 1) in the triangular coordinate system. The triangles can be pushed in those three directions without tearing the grid apart. Those transforms are linked to the following matrices in the triangular coordinate system: Note that a row corresponds to a set of triangles with aligned grid points. As such, a given B-row or a given C-row is composed of the same type of triangles (odd or even). Only an A-row is composed of odd and even triangles.

matA tri = 1 a 0 1 matB tri = 1 + b -b b 1 -b matC tri = 1 -c -c c 1 + c

Bijective Digital Rotation

Let us now find values for a, b and c such that the composition of three such shear transforms corresponds to a rotation of center (0, 0) and angle θ. For the sake of simplicity, we will write simply a, b, c instead of a(θ), b(θ) and c(θ). Since it is simpler to express the rotation matrix in the classical Cartesian coordinate system, we solved the following equation:

T ri2Cart•matA•matB•matC•Cart2T ri = cos θ sin θ -sin θ cos θ
which has a unique solution for a, b and c (function of θ):

a = -2 + -6 + 4 √ 3 sin θ -3 cos θ + √ 3 sin θ (3) b = 1 6 3 -3 cos θ + √ 3 sin θ (4) c = - a 4
(5)

Dealing with divergence

Let us first note that we consider here, that all angles are between -π and π. The solutions (a, b, c) do not represent a complete solution since the denominator of a (and equivalently c) has a value zero for angles θ equal to -2π/3 and π/3. The zeros for both angles are not of the same nature.

As can be seen in Figure 3, the zero in π/3 is a 0/0 singularity where lim θ→ π 3 (a) = -1 and lim θ→ π 3 (c) = 0.25. To avoid this singularity problem, one can simply set the values of coefficients a and c by their limit value for angle θ = π/3 (see line 1 in Algorithm 1 and Algorithm 2).

The problem for angle -2π/3 has to be dealt in an entirely different way. As one can see in Figure 4, the error measures (see section 4) start to increase rapidly when approaching the angle -2π/3. The idea here is to avoid altogether the problem of the interval from -π to 0 and simply consider that a digital rotation of such an angle is the inverse of the digital rotation of angle -θ (that falls then in the interval 0 to π). That way we avoid the problems around -2π/3 and ensure that the digital rotation of angle θ is the inverse of the digital rotation of angle -θ (note that it is an important side effect as a bijective transform is not always easily invertible). Note that the decomposition proposed by [START_REF] Paeth | A fast algorithm for general raster rotation[END_REF] and [START_REF] Tanaka | A rotation method for raster image using skew transformation[END_REF] has a similar divergence problem around angle π.

Algorithms

Let us now focus on the algorithm. We have found a shear based decomposition for a rotation that preserves A, B and C-rows respectively. Of course, a triangular cell is only transformed into another triangular cell if the coefficients a, b and c are integers. For this, we simply considered the Floor(x + 0.5) function (Floor(x) is the biggest integer smaller or equal to x). Careful here because the classical rounding function might lead to errors since it is not translation invariant. Some attention has A simple way around this problem is to consider that for all angles close to π/3, we consider a b value equal to 0.5 and a c value equal to 0.25. The term close here can be defined based on the size of the image that will be rotated. That is also the reason why this algorithm is set up for angles Algorithm 1: RotCBA(x, y, θ): Point Rotation of center (0, 0) and angle θ

Input : (x, y) ∈ Z 2 , 0 ≤ θ ≤ π Output: (x 3 , y 3 ) ∈ Z 2 If θ = π 3 Then a c ← -1 0.25 1 Else a ← -2 + -6+4 √ 3 sin θ -3 cos θ+ √ 3 sin θ 2 c ← -a/4 3 b = 1 6 (3 -3 cos θ + √ 3 sin θ). 4 x 1 y 1 ← x -Floor(0.5 + c * (x + y)) y + Floor(0.5 + c * (x + y)) 5 x 2 y 2 ← x 1 + Floor(0.5 + b * (x 1 -y 1 )) y 1 + Floor(0.5 + b * (x 1 -y 1 )) 6 x 3 y 3 ← x 2 + Floor(0.5 + a * y 2 ) y 2 7 return (x 3 , y 3 ) 8
between -π and π rather than for angles between 0 and 2π. This leads to the first two algorithms that correspond to digital cell rotations: Algorithm 1 for rotations of angles between 0 and π. The second algorithm, Algorithm 2, is exactly the same algorithm but all the operations of Algorithm 1 are performed in reverse order with negative push coefficients: x + Floor(0.5 + a * y) becomes x -Floor(0.5 + a * y). This ensures that the digital rotation is invertible and that the rotation of angle θ is the inverse of rotation of angle -θ. Algorithm 3 combines both Algorithm 1 and Algorithm 2 to perform a digital rotation of a triangular cell grid. Finally Algorithm 4 shows the final digital algorithm for a triangular grid image.

Bijectivity

All has been done in order to have a bijective transform RotT ri. Let us show that formally by examining the mappings of Algorithm 1. Let us define the mapping A(a) : (x, y) → (x + Floor(0.5 + ay), y), a ∈ R 2 (line 7 in Algorithm 1). For an A-row, we have A(a)(S ) = S (y does not change). All the triangular cells on a given Arow are pushed in direction A by the same integer number of grid points. Since A-rows partition the grid points, it proves that mapping A(a) is bijective over the grid.

Algorithm 2: RotCBANeg(x, y, θ): Inverse Point Rotation of center (0, 0) and angle θ

Input : (x, y) ∈ Z 2 , 0 ≤ θ ≤ π Output: (x 3 , y 3 ) ∈ Z 2 If θ = π 3 Then a c ← -1 0.25 1 Else a ← -2 + -6+4 √ 3 sin θ -3 cos θ+ √ 3 sin θ 2 c ← -a/4 3 b = 1 6 (3 -3 cos θ + √ 3 sin θ). 4 x 1 y 1 ← x -Floor(0.5 + a * y) y 5 x 2 y 2 ← x 1 -Floor(0.5 + b * (x 1 -y 1 )) y 1 -Floor(0.5 + b * (x 1 -y 1 )) 6 x 3 y 3 ← x 2 + Floor(0.5 + c * (x 2 + y 2 )) y 2 -Floor(0.5 + c * (x 2 + y 2 )) 7 return (x 3 , y 3 ) 8
Algorithm 3: RotPoint(x, y, θ): Rotation of center (0, 0) and angle θ of point (x,y) in the triangular Grid. The proof for C-rows is similar. All this is, of course, verified as well for Algorithm 2.

Input : (x, y) ∈ Z 2 , θ ∈ R 2 Output: (x , y ) ∈ Z 2 θ = θ Mod 2π 1 If 0 ≤ θ ≤ π Then (x , y ) ← RotCBA(x,
This proves that the transform presented in Algorithm 4 is bijective (as sequence of bijective mappings). By design, we have made sure that the transform is also easily invertible such that, for an image ImageIn in the triangular grid, we have RotT ri(ImageIn, θ) -1 = RotT ri(ImageIn, -θ). Figure 2 shows how each matrix acts on the triangular grid. In order to illustrate the action of each shear transform, we set a = b = c = 1. The rows containing the rotation center (0, 0) do not move. Note that mapping B and mapping C conserve the type (even or odd) of the triangles which is not the case for mapping A that changes the types of the triangles depending if Floor(0.5 + ay) is odd or even.

Results

In this section, we present some results and error measurements. The algorithms have all been implemented in Mathematica.

Error Measure

In order to evaluate the quality of the rotation, two distance criterii [START_REF] Andres | The quasi-shear rotation[END_REF][START_REF] Andres | Shear based bijective digital rotation in hexagonal grids[END_REF] were considered: the Maximum Distance criteria (MD) and the Average Distance criteria (AD). We measure the average and the maximal distance between the digital points obtained after the digital rotation and the Euclidean points obtained after the corresponding continuous rotation. The goal is to measure the error committed by considering the digital point instead of the continuous point. Let us recall here that the triangular cells have sides of size 1 and two neighboring triangles centroids are at a distance of 0.57.

For the first error measure, we simply considered the average and maximal distance criterii for our triangular grid. We obtained an average distance bounded by 0.52 (mean average distance over all the angles is 0.40). The maximal distance is bounded by 1.47 (mean maximal distance over all the angles of about 0.93). In order to be consistent with a more classical way of representing a triangular grid, we recomputed our error measures with the centroids rather than our corrected grid points. The measure becomes therefore: at what average and maximal distance is the continuously rotated point from the centroid of the digitally rotated corresponding triangle (see Figure 5). Let us call Rot θ (x, y) the continuous rotation transform of center (0, 0) and angle θ for point (x, y).

For a point (x, y) of the triangular grid:

• Continuous rotation: (x, y) → (x , y ) = T ri2Cart(x, y) ± 0, √ 3 12 → (x , y ) = Rot θ (x , y ). • Digital rotation: (x, y) → (x 1 , y 1 ) = RotPoint(x, y, θ) → (x 2 , y 2 ) = T ri2Cart(x 1 , y 1 ) ± 0, √ 3 12
The error measure is therefore defined as : for all (x, y) ∈ Z 2 : average and maximal distance value between the points (x , y ) and (x 2 , y 2 ). In order to compute the final error (Figure 5), we considered a triangular grid of size 1000 2 for angles from -π to π with steps of kπ/600. Figure 5 shows an average error of about 0.5 (maximum average error is 0.56 near. Mean average error over all angles is around 0.42). The maximal distance is bounded by 1.6.

One last question was to verify that there is no rotation center problem since the rotation center for the Euclidean rotation is Euclidean point (0, 0) while the rotation center used for our triangular grid rotation is the one with triangular grid coordinates (0, 0) and therefore with Euclidean coordinates 0, √ 3 12 . There is however no reason for rotation center error since the error between both ways of considering the point representing a cell creates an additional small error localized around a cell. To validate these values and ensure that there is no drift when considering triangles at a long distance, we took 10 6 randomly chosen points at a distance of 10 9 from the center and obtained similar results than the one presented in Figure 5.

Examples

Lastly, Figure 6 shows the result of the bijective rotation on a 'Lena' image of size 64x64 for different angles. The holes that appear on the boundary of the rotated image is not a problem of bijectvity but simply because the rotation does not guarantee that a straight line is transformed into a straight line.

Conclusions

In this paper, we have presented a method to build bijective digital rotations on a triangular grid using three shear transforms. The main advantages of the method is that it works for all angles, is one to one, easy to implement and provides a good approximation of the continuous rotation. In order to evaluate the error committed using such a digital rotation when compared to a continuous rotation, the average distance, on an image, between the point rotated by continuous rotation and the one rotated by digital one has been computed. Experiments showed that this average distance is not too high and oscillates around 0.5.

Figure 1

 1 Figure 1: a) Coordinate system on the triangular grid. The yellow triangles are called Even, the blue ones are called Odd. The centroid of each triangle (empty circle) on the grid doesn't form a regular grid which brought us to translate the centroids by + √ 3/12 and -√ 3/12 for respectively Even and Odd triangles on the Cartesian grid to obtain a regular grid (full red circles). b) Triangle size.

Figure 2 Figure 2 :

 22 Figure2shows the effect of the shear transforms for coefficients a = b = c = 1. In order to explain what those matrices mean, let us define three type of cell rows: let us call an A-row a set of cells whose triangular ty-coordinate is constant: A s = {(x, y) ∈ R : y = s} is an A-row. For a point (x, y), matA tri • (x, y) = (x + ay, y), so matA tri • A y =

Figure 3 :Figure 4 :

 34 Figure 3: The graph of coefficient a = -2 + -6+4 √ 3 sin θ -3 cos θ+ √ 3 sin θ .

4

  y, θ) 2 Else (x , y ) ← RotCBANeg(x, y, 2π -θ) RotT ri(ImageIn, θ): Bijective Rotation for the triangular Grid of center (0, 0) and angle θ Input : ImageIn, θ Output: ImageOut For all (x, y) in ImageIn 1 (x , y ) ← RotPoint(x, y, θ) 2 ImageOut (x , y ) ← ImageIn(x, y) Let us now consider the mapping B(b) : (x, y) → (x + Floor(0.5 + b(xy)), y + Floor(0.5 + b(xy))), b ∈ R 2 (line 6 in Algorithm 1). Let us show that for a B-row S , we have B(b)(S ) = S (i.e. that xy does not change). Let us consider a point (x, y) ∈ S and (x , y ) = B(b)(x, y). We have xy = xy + Floor(0.5 + b(xy)) -Floor(0.5 + b(xy)) = xy. All the triangular cells on a given B-row are pushed in direction B by the same integer number of grid points. Since B-rows partition the grid points, it proves that mapping B(b) is bijective over the grid.

Figure 5 :

 5 Figure 5: Average (AD) and Maximal (MD) Distance error criteria for the triangular grid Rotation Algorithm for angles 0 to 2π

Acknowledgments

The authors gratefully acknowledge the help of Ranita Biswas, Assistant-Professor from IIT Roorkee, India.


Lei Mo
email: lei.mo@inria.fr

Angeliki Kritikakou
email: angeliki.kritikakou@irisa.fr

Olivier Sentieys
email: olivier.sentieys@irisa.fr

Imprecise Computation Task Mapping on Multi-Core Wireless Sensor Networks

Keywords: Imprecise Computations, Multicore Architectures, Task Mapping. Definitions

) considers that a task can be logically decomposed into a mandatory subtask and an optional subtask. The mandatory subtask should be completed before the task's deadline in order to generate the minimum acceptable Quality of Service (QoS). The optional subtask is to be executed after the mandatory subtask, and still before the deadline, if there are resources in the system that are free to execute the optional subtask. The longer the optional subtask is executed, the better is the QoS of the obtained result.

1. Most of the sensor nodes have energy constraints, as they are powered by a battery that has a limited energy budget. 2. Real-time responsiveness is required by many WSN applications, e.g., mobile target tracking. Missing of the task deadline can cause serious results.

In the past, many research works have focused on optimizing the energy efficiency and real-time execution of applications on WSNs, e.g., data aggregation [START_REF] Dhasian | Survey of data aggregation techniques using soft computing in wireless sensor networks[END_REF], topology control [START_REF] Esch | A survey on topology control in wireless sensor networks: taxonomy, comparative study, and open issues[END_REF], sink mobility [START_REF] Silva | Mobility in wireless sensor networks: a survey and proposal[END_REF] and delay-aware routing [START_REF] Lai | Energy efficient link-delay aware routing in wireless sensor networks[END_REF]. These methods improved significantly the energy efficiency and real-time execution of the WSNs. However, further improvements exist:

1. Most of the past works focused on the software optimization methods without including the hardware optimization strategies. 2. Most of the past works focused on reducing the energy consumption of the wireless transmission, but ignored the energy consumption of the processor modules.

However, processing and transmitting large amount of sensed data in real-time applications exceeds the capabilities of traditional WSNs. Since single-core embedded sensor nodes will soon be unable to meet the increasing requirements of data intensive applications, e.g., Wireless Multimedia Sensor Networks (WMSNs), the next generation sensor nodes must possess enhanced computation and communication capabilities, i.e., multi-core embedded sensor nodes [START_REF] Munir | Multi-core embedded wireless sensor networks: architecture and applications[END_REF]. For instance, the fire detection WMSN node MiLive-v2 in ( 13) is equipped with three types of processors:

1. A 8-bit low-power AVR processor (ATmega128rfa1), which can be used to run simple tasks. 2. A 32-bit powerful ARM processor (ARM1176JZF), which can be used to run more complicated tasks. 3. A Digital Signal Processor (DSP) unit, which can be specially used to perform image processing.

The requirements of WSNs are usually high performance computation, low energy consumption and low task execution delay. The efficient mapping of the tasks on multi-core sensor nodes (13) is required to balance these contradictory requirements. It can simultaneously optimize task allocation (on which processor each task is executed) and task scheduling (when each task starts and ends its execution). On the one hand, in order to map a set of applications on multi-core system, the applications need to be partitioned (parallelized), whenever possible, into multiple tasks that can be executed concurrently on different processors. Therefore, task execution delay can be further reduced. On the other hand, multi-core systems have better energy savings than traditional single-core systems in two ways:

1. Multi-core embedded sensor nodes can extract the desired information from the sensed data and process this information. It leads to reduction in the data transmission volume sent to the sink node (base station). By replacing a large percentage of communication with in-network computation, the multi-core system realizes large energy savings that increase the sensor network's overall lifetime.

2. A multi-core embedded sensor node allows the computations to be split across multiple processors, while running each processor at a lower voltage and frequency, as compared to a single-core system, which results in energy savings. Utilizing a single-core embedded sensor node for information processing in data intensive applications demands the sensor node run at a high voltage and frequency in order to meet the application's real-time requirements. It will increase the power dissipation of the processor.

The tradeoff between energy consumption and performance has been extensively studied using the following methods:

1. Dynamic Voltage and Frequency Scaling (DVFS): by lowering the supply voltage, quadratic savings in dynamic energy consumption can be achieved, while the performance is approximately degraded in linear manner. 2. Dynamic Power Management (DPM): by exploring idle intervals of the processors and switching off a processor, when the idle interval is longer than a certain threshold, which is called break-even time.

Existing works that focus on energy-aware task mapping problem aim at minimizing the energy consumption under resource and application constraints.

1. When the voltage level is discrete, the task mapping problem is usually formulated as Integer Programming (IP), e.g., (14; 22; 11). To efficiently solve the IP problem, a hybrid Genetic Algorithm (GA) is proposed in ( 14), a polynomialtime two-step heuristic is designed in [START_REF] Xu | Energy minimizing for parallel real-time tasks based on levelpacking[END_REF], and the IP problem is relaxed to a Linear Programming (LP) in [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF]. Combining DVFS and DPM, a Mixed-Integer Linear Programming (MILP)-based task mapping problem is considered in (3) and the problem is solved by CPLEX solver. 2. When the voltage level is continuous, a convex task mapping problem is formulated in [START_REF] Kong | Energy-efficient scheduling of real-time tasks on cluster-based multicores[END_REF] and the problem can be solved by using polynomial-time methods. In [START_REF] Leung | Simultaneous task allocation, scheduling and voltage assignment for multiple-processors-core systems using mixed integer nonlinear programming[END_REF], Mixed-Integer Non-Linear Programming (MINLP) is used to formulate the task mapping problem. The problem is relaxed to an MILP by linear approximation and solved by Branch and Bound (B&B) method.

As long as the resource and application constraints (e.g., task deadline) allow, the methods with DVFS or DPM achieve significant energy reductions compared to the basic methods, where no DVFS or DPM is applied.

Key Applications

In several application domains, an approximate-but still in time-result is preferable than an exact result obtained too late. For example, in audio and video streaming, frames with a lower quality are better than missing frames. In radar tracking, an estimation of target's location in time is better than an accurate location arriving too late. In control loops, an approximate result produced by a control scheme is more preferable as long as the controlled system, e.g., indoor temperature control system, remains stable. This type of applications can be modeled as Imprecise Computation (IC) tasks. The more cycles of the optional subtasks are executed, the higher is the generated QoS. An interesting question is how should the function of QoS be defined in order to represent realistic application areas:

1. A linear function (25; 23; 26; 21; 16; 17) models the case where the QoS of the overall system increases uniformly during the optional subtask execution. 2. A concave function (2; 24; 4; 20) addresses the case where the increase of QoS exhibits a continuously nondecreasing rate as the optional subtask execution goes on.

Linear and general concave functions are considered as the most realistic and typical QoS representation in the literature, as they adequately capture the behavior of many application areas, such as image and speech processing, control engineering, and automatic target recognition (1; 2; 12).

Foundations

The three constraints energy, deadline, and QoS play an important roles in the QoS-aware task mapping. This is different from the energy-aware task mapping, where no exploration of the QoS improvement through the optional subtasks adjustment take place. The QoS-aware task mapping problem uses the IC task model and has a goal to maximize the QoS under a set of real-time and energy supply constraints, e.g., (2; 4; 24; 25; 20; 23; 7; 26; 21; 16; 17). Specifically, QoS-aware task mapping determines on which processor should the task be executed (task-toprocessor allocation), its frequency (frequency-to-task assignment), the start time (task scheduling), and the end time (optional subtask adjustment) of each task, such that the system QoS is maximized, while meeting the energy supply and the task deadlines. Usually, the task-to-processor allocation and frequency-to-task assignment decisions are binary variables, while the task scheduling and optional subtask adjustment decisions can be either integer (24; 25; 23; 7; 26; 21) or continuous (2; 4; 20; 16; 17) variables. Based on different constraints and variables introduced into the problem formulation, the variables may be coupled with each other nonlinearly, i.e., MINLP. However, it is not desirable since solving MINLP is much more complex than solving MILP. In order to linearize the nonlinear items, the common methods include:

1. Linear approximation (10); 2. Variables replacement (3; 17).

The QoS-aware task mapping problem is similar to Quadratic Assignment Problem, a well-known NP-hard problem. Therefore, finding an optimal solution satisfying all the given constraints (e.g., energy efficiency, deadline, QoS, task dependency, and DVFS) is very difficult and time consuming [START_REF] Mo | Controllable QoS for imprecise computation tasks on DVFS multicores with time and energy constraints[END_REF]. Existing works usually focus on different contexts, as summarized in Table 1. The platforms considered in (4; 24; 20) are single-core platforms, i.e., the task allocation problem is not taken into account. Although the works in (2; 25; 23; 7; 26; 21; 16) target at multi-core platforms, some assumptions are made during the problem formulation, e.g., the task-to-processor allocation is fixed and given in advance for all the tasks in (25; 23), the tasks considered in [START_REF] Isabel | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF] are independent, and in (2; 26; 21; 16) each processor has a predefined frequency. The methods used to solve the aforementioned problems can be classified into two main classes:

1. The first class includes the methods based on heuristics, e.g., (24; 25; 23; 7; 26; 21). 2. The second class includes the methods that always produce an optimal solution, e.g., (2; 4; 20; 16; 17).

Although heuristic methods can find feasible solutions in a short amount of time, they do not provide the bounds on the solution quality. When new assumptions or constraints are taken into account they must be redeveloped, and, thus, they are not always easily extensible. The heuristic methods mentioned above usually adopted a multi-step optimization, i.e., to decouple the variables and to determine their values in sequence. For instance, a two-step heuristic called Adaptive Task Allocation (ATA) is proposed in [START_REF] Zhou | Energy-adaptive scheduling of imprecise computation tasks for QoS optimization in real-time MPSoC systems[END_REF]. The aim of the first step is to find a task-to-processor allocation, such that the energy consumption is minimized. With the given task allocation decision, the energy consumed to execute one task is proportional to the length of its optional subtask. Based on the given task allocation decision, the aim of the second step is to adjust the optional subtasks so as to maximize OoS under the energy supply constraint. On the other hand, in order to find the optimal solution, the common methods include: a lower bound for the original problem along with the temporary values of the binary variables. And then the SP is solved to obtain an upper bound by utilizing the MP solution. The bounds are updated if the stopping criterion, i.e., the gap between the upper and lower bounds is smaller than a predefined threshold, is not met, and a new constraint (Benders cut) is generated using the SP solution and is added to MP in next iteration.

1 .Fig. 1

 11 Fig.1The structure of Benders decomposition.

Table 1

 1 Task Mapping Method

	Energy-aware

Non-optimal√ √ √ √ √ √ √ √ √ √ √


Abderrahman Mokni
email: abderrahman.mokni@mines-ales.fr

Marianne Huchard
email: huchard@lirmm.fr

Christelle Urtado
email: christelle.urtado@mines-ales.fr

Vauttier + Sylvain
email: sylvain.vauttier@mines-ales.fr

Yulin Zhang
email: yulin.zhang@u-picardie.fr

A three-level formal model for software architecture evolution

Keywords: Software architecture, architecture levels, reuse, software evolution, B formal models

This papers gives an overview of our formal approach to address the architecture-centric evolution at the three main steps of component-based software development: specification, implementation and deployment. We illustrate our proposal with an example of software evolution that leads to erosion and we demonstrate how our evolution process can resolve this problem.

Introduction

Software evolution has gained a lot of interest during the last years [START_REF]Evolving Software Systems[END_REF]. Indeed, as software ages, it needs to evolve and be maintained to fit new user requirements. This avoids to build a new software from scratch and hence save time and money. Handling evolution in component-based software systems is non trivial since an ill-mastered change may lead to architecture inconsistencies and incoherence between design and implementation. Many Adls (Architecture Description Languages) were proposed to support architecture modeling and analysis. Examples include C2SADL [START_REF] Medvidovic | ADLs and dynamic architecture changes[END_REF], Wright [START_REF] Allen | A formal basis for architectural connection[END_REF] and Darwin [START_REF] Magee | Dynamic structure in software architectures[END_REF]. Although, some Adls integrate architecture modification languages, handling and controlling architecture evolution in the overall software lifecycle is still an important issue. In this paper, we attempt to provide a solution to the architecture-centric evolution that preserves consistency and coherence between architecture levels. We propose an architecture evolution process based on the formal foundations of our three-level Adl Dedal [START_REF] Zhang | Architecture-centric component-based development needs a three-level ADL[END_REF]. The process is then illustrated by an evolution scenario on a simplified Home Automation Software architecture. The remainder of this paper is organized as follows: Section 2 briefly discusses examples of existing Adls. Section 3 gives an overview of Dedal and its formal foundations. Section 4 describes the evolution process based on Dedal. Section 5 presents an evolution scenario that illustrates the proposed evolution process before Section 6 concludes and discusses future work.

Over the two last decades, a number of Adls were proposed [START_REF] Medvidovic | A classification and comparison framework for software architecture description languages[END_REF]. Most of them provide textual notations to describe architectural entities (i.e. components, interfaces and connections). Initially, Adls were domain-specific. Examples include C2SADL [START_REF] Medvidovic | ADLs and dynamic architecture changes[END_REF] for the design of concurrent systems and Wright [START_REF] Allen | A formal basis for architectural connection[END_REF] and Darwin [START_REF] Magee | Dynamic structure in software architectures[END_REF] for the design and analysis of distributed architectures. Later on, attempts to unify Adls and make them general-purpose were made. For example, ACME [START_REF] Garlan | Acme: An architecture description interchange language[END_REF] was designed for such purpose. It consists in a common interchange description language that offers annotation facilities to support architecture descriptions in other languages. Another relevant example is xADL 2.0 [8]. It was designed to support various types of systems. The strength of xADL 2.0 resides in its extensibility since it is XML-based. It offers then an easy way for architects to adapt its use to any kind of architectures.

Although a lot of effort was dedicated to improve the expressiveness of Adls and promote their use to model software architectures, several important issues are not taken into account. First, existing Adls hardly support all the development steps of component-based software architectures (i.e. specification, implementation and deployment). Most of them cover only one or at most two levels which harden their integration in development processes. Second, they hardly support architecture evolution and handle architecture inconsistencies such as drift or erosion [START_REF] Perry | Foundations for the study of software architecture[END_REF]. C2SADL and Darwin are exceptions. They include language to describe changes. However, they do not support reverse evolution and they do not cover all steps of component-based development either.

In this work, we address architecture evolution in the whole component-based development process. We show how architectural erosion could be avoided thanks to reverse evolution.

Overview of Dedal

The three architecture levels

Dedal is a novel ADL that covers the whole life-cycle of a component-based software. It proposes a three-step approach for specifying, implementing and deploying software architectures in a reuse-based process [START_REF] Zhang | A three-level component model in component-based software development[END_REF].

To illustrate the three architecture levels of Dedal, we propose an example of a Home Automaton Software (Has). Figure 1 presents the Has architecture at three abstraction levels:

The abstract architecture specification (cf. Figure 1-a) is the first level of architecture software descriptions. It represents the architecture as designed by the architect and after analyzing the requirements of the future software. In Dedal, the architecture specification is composed of component roles and their connections. Component roles are abstract and partial component type specifications. They are identified by the architect in order to search for and select corresponding concrete components in the next step. The instantiated architecture assembly (cf. Figure 1-c) describes software at runtime and gathers information about its internal state. The architecture assembly results from the instantiation of an architecture configuration. It lists the instances of the component and connector classes that compose the deployed architecture at runtime and their assembly constraints (such as maximum numbers of allowed instances).

The formal foundations of Dedal

Dedal is formalized using B [START_REF] Abrial | The B-book: Assigning Programs to Meanings[END_REF], a set-theoretic and first order logic formalism. The formalization [START_REF] Mokni | Towards automating the coherence verification of multi-level architecture descriptions[END_REF] covers all the concepts of Dedal and includes a set of rules that defines the relations between the different artifacts into and over each architecture level of Dedal (cf. Figure 2). These rules are classified into two categories: the intra-level rules and the inter-level rules.

Fig. 2. Inter-level relations in Dedal

Intra-level rules in Dedal consist in substitutability and compatibility between components of the same abstraction level (component roles, concrete component types, instances). Defining intra-level relations is necessary to check the architecture consistency. For instance, the components must be correctly connected to each other (i.e. each required interface is connected to a compatible provided one).

Inter-level rules are specific to Dedal and consist in relations between components at different abstraction levels as shown in Figure 2. Defining inter-level rules is mandatory to decide about coherence between two architecture descriptions at different abstraction levels. For instance, the realization rule is used to check that a given configuration is a valid implementation of a given specification and the instantiation rule is used to check if an assembly correctly instantiates a given configuration.

Handling software evolution in Dedal is quite advantageous. Indeed, Dedal covers the whole life-cycle of software systems and hence all descriptions can be kept upto-date for further reuse, reimplementation and deployment in different contexts.

To keep architecture descriptions coherent, change must be propagated from where it is initiated to the other abstraction levels descriptions. As a solution, we propose an evolution process based on Dedal to enable change in software systems in a manner that preserves architecture consistency and coherence.

The evolution process in Dedal lies on two kinds of rules: (1) static rules to check architecture consistency and coherence between the different descriptions and (2) evolution rules to trigger the change at each abstraction level and propagate it to the other levels descriptions. Figure 3 presents the condition diagram of the evolution process. Consistency rules are name uniqueness, completeness, connection correctness and graph consistency. These properties are verified at the same abstraction level by the evolution manager to check weather the architecture is structurally consistent or a change must be triggered to restore consistency. Coherence rules are used to check if software descriptions at the three abstraction levels of Dedal are coherent. If an incoherence is detected, the evolution manager propagates change to the other levels to restore coherence. Coherence rules include the verification of the following properties:

-A configuration Conf is an implementation of a given specification Spec.

-A specification Spec is a documentation of a given configuration Conf.

-An assembly Asm is an instantiation of a given configuration Conf.

-A configuration Conf is instantiated by a given assembly Asm.

Evolution rules

An evolution rule is an operation that makes change in a target software architecture by the deletion, addition or substitution of one of its constituent elements (components and connections). Each rule is composed of three parts: the operation signature, preconditions and actions. Specific evolution rules are defined at each abstraction level to perform change at the corresponding formal description. These rules are triggered by the evolution manager when a change is requested. Firstly, a sequence of rule triggers is generated to reestablish consistency at the formal description of the initial level of change. Afterward, the evolution manager attempts to restore coherence between the other descriptions by executing the adequate evolution rules. The following role addition rule is an example of evolution rules at specification level: 5 Evolution scenario

Motivation

To illustrate the evolution approach, we propose an example of evolving the Has architecture. The objective is to enable the control of the house through a mobile device running under Android OS. The change is initiated at the configuration level and attempts to adapt the current HAS implementation to an android device.

Figure 4-a shows the initial implementation of HAS while Figure 4-b shows the evolved one. Two main changes are noticed in the new configuration: the orchestrator is substituted for a new one compatible with android. Since a new service to control the intensity of lamp is required, the component Lamp is replaced by the component AdjustableLamp with additional provided interface (IIntensity) to adjust the luminosity.

Tool support overview

At this stage of work, the evolution process is assisted using ProB [START_REF] Leuschel | Prob: An automated analysis toolset for the b method[END_REF], an animation tool of B models. We manually instantiate the B formal models corresponding to the Has architecture and execute evolution rules at each abstraction level. We control the evolution process by checking consistency and coherence properties thanks to the evaluation console of ProB. This first step provides a proof feasability of our work. Ongoing work is to automate the generation of Dedal formal models and use the ProB solver to automate the evolution process.

Evolving the HAS configuration

The change is initiated by disconnecting and deleting the old orchestrator. Then, the one compatible with android is added and connected. The old Lamp is replaced by the adjustable one and connected to the android orchestrator. The evolution manager performs the following operations: disconnect(HAS conf ig, (cl4, rintIP ower), (cl1, pintIP ower)) disconnect(HAS conf ig, (cl4, rintIClock), (cl2, pintIClock)) disconnect(HAS conf ig, (cl4, rintIT herm), (cl3, pintIT herm)) disconnect(HAS conf ig, (cl4, rintIAirCon), (cl3, pintIAirCon)) deleteClass(HAS conf ig, cl4) addClass(HAS conf ig, cl4a) replaceClass(HAS conf ig, cl1, cl1a) connect(HAS conf ig, (cl4a, rintIP ower2), (cl1a, pintIP ower2)) connect(HAS conf ig, (cl4a, rintIIntensity), (cl1a, pintIIntensity)) connect(HAS conf ig, (cl4a, rintIClock2), (cl2, pintIClock)) connect(HAS conf ig, (cl4a, rintIT herm2), (cl3, pintIT herm)) connect(HAS conf ig, (cl4a, rintIAirCon2), (cl3, pintIAirCon2))

We note that cl1, cl2, cl3 and cl4 refer respectively to the component classes Lamp, Clock, AirConditioner and Orchestrator. cl1a and cl4a refer respectively to the new component classes AdjustableLamp and AndroidOrchestrator. Their Interface names are prefixed with rint(for a required interface) and pint(for provided interface) followed by the interface type name and eventually a number when there is several instances of the same interface

Propagating change to the HAS specification

The current Has specification is no longer a good documentation of the new version of the HAS configuration and thus, we have a problem of erosion. Indeed, the control of the light intensity is not included in the current specification. Hence, a new documentation version is required to keep both descriptions coherent. Figure 5 shows the initial and evolved version of the HAS specification after the change propagation. The change is propagated to the HAS specification by replacing the role HomeOrchestrator (cr4) with HomeOrchestrator2 (cr4a) and the role Light(cr1) by Lunminosity(cr1a). The following operations are performed by the evolution manager: disconnect(HAS spec, (cr4, rintILight), (cr1, pintILight)) disconnect(HAS spec, (cr4, rintIT ime), (cr2, pintIT ime)) disconnect(HAS spec, (cr4, rintIT herm1), (cr3, pintIT herm)) disconnect(HAS spec, (cr4, rintICon), (cr5, pintICon)) deleteRole(HAS spec, cr4) addRole(HAS spec, cr4a) replaceRole(HAS spec, cr1, cr1a) connect(HAS spec, (cr4a, rintILight2), (cr1a, pintILight2)) connect(HAS spec, (cr4a, rintIIntensity), (cr1a, pintIIntensity)) connect(HAS spec, (cr4a, rintIT ime), (cr2, pintIT ime)) connect(HAS spec, (cr4a, rintIT herm1), (cr3, pintIT herm)) connect(HAS spec, (cr4a, rintICon), (cr5, pintICon))

Propagating change to the HAS assembly

The current HAS assembly violates the instantiation rule according to the new version of the HAS configuration. This violation is detected by the evolution manager and change is triggered at the assembly level to restore coherence. Figure 6 illustrates the changes applied on the assembly architecture. In this paper, we give an overview of our three-level Adl Dedal and its formal model. At this stage, a set of evolution rules is proposed to handle architecture change during the three steps of software lifecycle: specification, implementation and deployment. The rules were tested and validated on sample models using a B model checker. As future work, we aim to manage the history of architecture changes in Dedal descriptions as a way to manage software system versions. Furthermore we are considering to automate evolution by integrating Dedal and evolution rules into an eclipse-based platform.

Fig. 1 .

 1 Fig. 1. Illustrative Example

Fig. 3 .

 3 Fig. 3. Condition diagram of the evolution process

/

  * Operation signature takes as arguments an instance of the architecture specification(spec) and the instance of the new role(newRole)*/ addRole(spec, newRole) = /* preconditions */ PRE spec ∈ arch spec ∧ newRole ∈ compRole ∧ newRole ∈ spec components(spec) ∧ /* spec does not contain a role with the same name*/ ∀ cr.(cr ∈ compRole ∧ cr ∈ spec components(spec) ⇒ comp name(cr) = comp name(newRole)) THEN /* actions */ /* update the set of clients (required interfaces), the set of servers (provided interfaces) and the set of component roles */ spec servers(spec) := spec servers(spec) ∪ servers(newRole) || spec clients(spec) := spec clients(spec) ∪ clients(newRole) || spec components(spec) := spec components(spec) ∪ {newRole} END;

Fig. 4 .

 4 Fig. 4. Evolving the HAS configuration

Fig. 5 .

 5 Fig. 5. Evolving the HAS specification by change propagation

Fig. 6 .

 6 Fig. 6. Evolving the HAS assembly