Nathan Crombez
email: nathan.crombez@u-bourgogne.fr

Guillaume Caron
email: guillaume.caron@u-picardie.fr

Takuya Funatomi
email: funatomi@is.naist.jp

Yasuhiro Mukaigawa
email: mukaigawa@is.naist.jp

Reliable Planar Object Pose Estimation in Light Fields From Best Subaperture Camera Pairs

Keywords:

come

I. INTRODUCTION

A. Motivation

T He pose estimation of an object is a fundamental task for a variety of purposes in the robotics field such as 3-D tracking [START_REF] Marchand | Pose estimation for augmented reality: A hands-on survey[END_REF], visual servoing [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF], and motion estimation [START_REF] Fraundorfer | Visual odometry : Part ii: Matching, robustness, optimization, and applications[END_REF]. When the target is a known planar object such as a checkerboard, a single camera can be used for this purpose, and a lot of geometric algorithms have been proposed. The Zhang's algorithm [START_REF] Zhang | A flexible new technique for camera calibration[END_REF] is a standard algorithm of intrinsic parameters estimation, based on corner feature points, and it also provides the poses of the checkerboard for each input views. A binocular stereo pair might also be considered for increasing the accuracy [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] and relaxing all the knowledge about the object dimensions, except it is planar.

However, Hartley and Zisserman [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] highlight that some feature points locations with respect to the pair of images leads to poor estimates of the plane parameters and pose, usually combined as a homography. More precisely, "image points [...] close to collinear with the epipole [leads to] poorly conditioned estimate" [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]Sec. 13.2.1,p. 330]. The latter issue is jointly raised by the points configuration and the cameras configuration, i.e. their relative pose. The problem of points configuration for the robust pose estimation of a planar object has very recently been studied in the literature in order to consider configurations that ensure a precise estimation [START_REF] Acuna | Robustness of control point configurations for homography and planar pose estimation[END_REF].

In this work, we propose a complementary approach in which we aim at selecting the best pair of cameras to ensure a precise estimation of the planar object pose. The selection of a pair of cameras among many is done in this work thanks to light-field (LF) imaging, contrary to the classical use of a stereo rig.

B. Related works

Recently, a LF camera can be used to obtain richer information of a scene, such as 3-D reconstruction [START_REF] Monteiro | Depth range accuracy for plenoptic cameras[END_REF], [START_REF] Hahne | Baseline and triangulation geometry in a standard plenoptic camera[END_REF], instead of usual 2-D cameras. The information obtained by the LF camera is inherently equivalent to a set of images captured from different viewpoints, so-called sub-aperture images (SAI). SAIs, made from the transformation of the image taken by the LF camera or made from several actual camera arrays are full of interest in robotics [START_REF] Tsai | Image-based visual servoing with light field cameras[END_REF], [START_REF] Dong | Plenoptic cameras in real-time robotics[END_REF]. This is due to the high redundancy and strong constraints brought by the configuration of sub-aperture cameras (SAC) in a LF camera, different from general multi-view settings.

The constraints inherent to SAIs have been considered in general Structure-From-Motion [START_REF] Johannsen | On linear structure from motion for light field cameras[END_REF] with the point feature. Zhang et al. [START_REF] Zhang | Ray space features for plenoptic structure-from-motion[END_REF] recently extended the research to the line feature. The latter work also introduces the modeling of planar constraint on feature point manifolds. However, that theoretical modeling is not considered in any motion or structure estimation.

C. Contributions

In this paper, we propose the detailed modeling of relationships between feature points belonging to a planar area of the scene, observed by one and two LF cameras. Then, linear and non-linear estimation algorithms of planar object pose are proposed, exploiting the minimum number of SAIs. To reach precise planar object pose estimation, we propose a new algorithm to select the best pair of SAIs. This selection avoids poor conditioning in the estimation and saves processing time, with respect to a greedy estimation that would consider every available SAI.

The contributions of the article are summed up below:

• The proposition of a method to select which pair of SAIs is the best for plane estimation.

• The exploitation of the best pairs of SAIs in two different

LFs to estimate the planar object motion. • The efficiency of the approach with a low processing time while being precise. Combining the latter contributions leads to the reliable estimation of planar object poses over time in LF sequences.

Theoretical and methodological contributions are deeply evaluated on real LF images with respect to the robotic armbased acquired ground truth. Comparisons are also done with respect to the extrinsic parameters estimation of a LF camera calibration toolbox [START_REF] Dansereau | Calibration and rectification for lenselet-based plenoptic cameras[END_REF], when considering a checkerboard as planar object, and with respect to the standard binocular stereo pair (e.g.: two horizontally aligned side-by-side cameras), when considering a textured planar object.

II. LIGHT-FIELD CAMERA MODEL

This section recalls how the image data acquired by a LF camera is transformed to the 4-D structure that parameterizes the acquisition. Throughout this paper, we use x ∈ R n d , a vector of Euclidean coordinates in the n d -dimensions (n d -D) space and x ∈ P n d , its homogeneous equivalent. Transformations from Euclidean to homogeneous and vice versa are omitted since they are very common.

In this paper, we consider the model of a LF camera introduced by Dansereau et al. [START_REF] Dansereau | Calibration and rectification for lenselet-based plenoptic cameras[END_REF]. The latter model is also considered in several computer vision applications [START_REF] Johannsen | On linear structure from motion for light field cameras[END_REF], [START_REF] Tsai | Image-based visual servoing with light field cameras[END_REF], [START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF]. This LF camera model describes

φ = [m, w, x, y, 1] , (1)
which fully describes both spatial (o) and angular (p) information about the incident light ray. However, physically, the LF camera has a unique photosensitive matrix. So, the four parameters of (1) are acquired in the same raw real image plane. Figure 1 shows a subpart of such a raw image on which one can see that it contains a set of small lenslet images.

Hence, a transformation is needed to get the 4-D φ from pixels of the raw acquired image. To do so, two main stages are required. First, omitting image processing of the raw image to segment lenslets [START_REF] Dansereau | Calibration and rectification for lenselet-based plenoptic cameras[END_REF], the initial common 2-D raw pixels structure is transformed to a 4-D structure (Fig. 2a). The latter organizes raw image pixel coordinates in order to build subsets containing one pixel, of coordinates i = [i, j] , of each lenslet of central coordinates k = [k, l] , with respect to which i is expressed. Stacking i and k leads to the decoded LF image "pixel" coordinates:

r = [i, j, k, l, 1] . (2)
Then, the 4-D structure is easily re-organized from r coordinates in order to build SAIs, as shown in Figure 2b.

As, in the rest of the article, we focus on considering φ's and not r's, r is transformed to φ thanks to the LF camera intrinsic parameters matrix K such that:

      m w x y 1       =       K 1,1 0 K 1,3 0 K 1,5 0 K 2,2 0 K 2,4 K 2,5 K 3,1 0 K 3,3 0 K 3,5 0 K 4,2 0 K 4,4 K 4,5 0 0 0 0 1             i j k l 1       . (3)
The twelve non-zero elements of the intrinsic matrix K are estimated thanks to the LF camera calibration toolbox 1 implementing [START_REF] Dansereau | Calibration and rectification for lenselet-based plenoptic cameras[END_REF].

Finally, since a LF camera captures several light rays emitted by a unique 3-D point P, there are as much φ's corresponding to P as there are SACs seeing P. Hence, we indicate the SAC which, virtually, acquires the ray as superscript left to φ, i.e. ca φ, for SAC c a of frame F ca , c b φ, for SAC c b of frame F c b , and so on. More precisely, we have: considering the coordinates of SAC optical centers are expressed in the global frame F 0 (e.g. 0 o ca = [ 0 m ca , 0 w ca ], regarding camera c a ). In (4), ca p = [ ca x, ca y] , is nothing but a precision, since coordinates x and y are expressed relatively to 0 o ca (1). Remark 1 (Coplanarity of SACs): By construction from the L2PP parameterization, optical centers of SACs belong to the same plane Π (the focal plane). That is why these optical centers coordinates are expressed in 2-D in F 0 .

ca φ = [ 0 m ca , 0 w ca , ca x, ca y, 1] , (4)

III. MULTIPLE VIEW GEOMETRIC RELATIONSHIPS IN A LIGHT-FIELD UNDER PLANAR CONSTRAINT

A LF camera captures several projections of the same 3-D point in one shot. Thus, we propose to consider an observed planar object and to study how its parameters are involved in SAIs relationships. First, we describe how feature correspondences in LF structures are expressed. Then, we introduce for the first time, to our knowledge, the homography matrix related to light rays.

A. Light ray correspondences

In an image captured with a regular camera, a single 3-D point ca P = [ ca X, ca Y , ca Z] , expressed in the c a camera frame F ca , is projected onto a unique 2D point ca p = [ ca x, ca y] in the sensor frame ( ca k = [ ca k, ca l] in the 4-D structure, see [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF], adding superscript c a to show to which SAC it belongs to). The usual input to many computer vision algorithms, for object pose estimation, is a set of 2D points correspondences ca k and c b k (from a second camera c b ) between two (or more) images.

In a calibrated LF camera, correspondences between two SACs c a and c b images are rays, parameterized as 4-D (Sec. II), below written in the L2PP modeling and as homogeneous coordinates:

ca φ ↔ c b φ . (5)
As in classical multiple view computer vision, many correspondences are usually available, i.e. ca φ1 ↔ c b φ1 , ca φ2 ↔ c b φ2 , and so on. Indices are not used in every equation of the paper to limit heavy mathematical writings.

B. Rays homography

Corresponding φ's are issued from the same 3-D point P. When P is lying on a planar object, its rays captured by the LF camera follows a homography projective transformation. Considering the common multiple view computer vision, we have, for two camera frames F ca and F c b , the projective homography matrix c b H ca (3×3) transforming ca P to c b P, up to a scale factor, is:

c b P ∝ c b H ca (3×3) ca P . (6)
Then, as a 3D point ca P is projected in the sensor frame of camera c a following the perspective projection:

ca p = λ ca ca P ⇐⇒   ca x ca y 1   = λ ca   ca X ca Y ca Z   , (7)
where λ ca = 1 ca Z , and similarly done for c b P, corresponding p's are related by the same c b H ca :

c b p ∝ c b H ca ca p   c b x c b y 1   ∝   H 1 H 2 H 3 H 4 H 5 H 6 H 7 H 8 1     ca x ca y 1   . (8)
The analytical expression of c b H ca is known to involve the rotation matrix c b R ca , belonging to the SO(3) group of the Lie algebra, the translation vector

c b t ca (3×1) = [ c b t X ca , c b t Y ca , c b t Z ca ]
, belonging to R 3 between the frame F ca and the frame F c b . It also involves the parameters of the planar surface on which the 3-D point belongs to [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]:

c b H ca = c b R ca - 1 ca d c b t ca ca n , (9)
where

ca n = [ ca n X , ca n Y , ca n Z ]
is the normal vector of the planar surface expressed in F ca and ca d is the orthogonal distance between the planar surface and the origin of the camera c a . Taking ( 8) and ( 9) together leads to the general analytical expression of each H {1,...,8} .

Whereas homographies above are, as often considered, 2-D, φ's are 4-D. Thus, as ca φ and c b φ are rays captured by SACs, crossing at P, we extend (8) to rays as:

c b φ ∝ c b H ca ca φ       0 m c b 0 w c b c b x c b y 1       ∝       H 1 H 2 H 3 H 4 H
            0 m ca 0 w ca ca x ca y 1      
. [START_REF] Dong | Plenoptic cameras in real-time robotics[END_REF] We can identify in [START_REF] Dong | Plenoptic cameras in real-time robotics[END_REF] the presence of the (8) which gives the following nine elements of c b H ca :

H 13 = H 1 H 14 = H 2 H 15 = H 3 H 18 = H 4 H 19 = H 5 H 20 = H 6 H 23 = H 7 H 24 = H 8 H 25 = 1
The latter elements of That translations between F ca and F c b can be written w.r.t the global frame F 0 :

c b t X ca = c b t X 0 + 0 t X ca = 0 t X ca + 0 t X c b = 0 m ca -0 m c b (11) c b t Y ca = c b t Y 0 + 0 t Y ca = 0 t Y ca + 0 t Y c b = 0 w ca -0 w c b . ( 12
)
Equations ( 11) and ( 12) provide alone the useful missing relationships to set the remaining sixteen elements of c b H ca (5×5) :

c b H ca (5×5) =         1 0 0 0 c b t X ca 0 1 0 0 c b t Y ca 0 0 H 1 H 2 H 3 0 0 H 4 H 5 H 6 0 0 H 7 H 8 1         . ( 13
)
Furthermore, assumptions considered to establish [START_REF] Johannsen | On linear structure from motion for light field cameras[END_REF] and ( 12) lead to the expressions of c b R ca and c b t ca :

c b R ca =   1 0 0 0 1 0 0 0 1   , c b t ca =   c b t X ca c b t Y ca 0   =   0 m ca -0 m c b 0 w ca -0 w c b 0   , (14)
thus simplifying c b H ca (9) as:

    1 - c b t X ca ca n X ca d - c b t X ca ca n Y ca d - c b t X ca ca n Z ca d - c b t Y ca ca n X ca d 1 - c b t Y ca ca n Y ca d - c b t Y ca ca n Z ca d 0 0 1     . (15)
Therefore c b H ca becomes:

         1 0 0 0 c b t X ca 0 1 0 0 c b t Y ca 0 0 1 - c b t X ca ca n X ca d - c b t X ca ca n Y ca d - c b t X ca ca n Z ca d 0 0 - c b t Y ca ca n X ca d 1 - c b t Y ca ca n Y ca d - c b t Y ca ca n Z ca d 0 0 0 0 1          . (16)
Remark 2 (Parallelism of SAC optical axes): From (11), SAC optical axes are considered to be parallel. Although it allows reaching ( 16) simple as it is, the assumption reduces the rest of the modeling and planar object pose estimation (Sec. IV and V) to the class of LF camera devices following a design having the same optical properties as the Lytro Photo LF camera device [START_REF] Knight | Light field data acquisition devices, and methods of using and manufacturing same[END_REF]. In the latter, lenses of the considered micro-lenses array, to make the camera a LF one, are coplanar and have their optical axes parallel. The latter design fact leads to the parallelism of SAC optical axes. Several other LF acquisition designs, including multi-camera ones [START_REF] Wilburn | The light field video camera[END_REF], can also benefit from the method we propose.

IV. PLANE PARAMETERS ESTIMATION

A. Linear estimation

Considering the analytical expression of c b H ca in [START_REF] Wilburn | The light field video camera[END_REF] and after some manipulations and rearrangements of (10), we can collect the plane parameters ca n and ca d on one side:

       c b x 1 -ca x 1 c b y 1 -ca y 1 . . . c b x * -ca x * c b y * -ca y *        =        c b t X ca ca x 1 c b t X ca ca y 1 c b t X ca c b t Y ca ca x 1 c b t Y ca ca y 1 c b t Y ca . . . . . . . . . c b t X ca ca x * c b t X ca ca y * c b t X ca c b t Y ca ca x * c b t Y ca ca y * c b t Y ca            ca n X ca d ca n Y ca d ca n Z ca d     (17)
where 11) and [START_REF] Zhang | Ray space features for plenoptic structure-from-motion[END_REF]. Therefore, the plane parameters can be computed directly from a set of at least 3 light rays correspondences ( * >= 3) following:

c b t X ca = ( 0 m ca -0 m c b ), c b t Y ca = ( 0 w ca -0 w c b ) as seen in (
ca η=        c b t X ca ca x 1 c b t X ca ca y 1 c b t X ca c b t Y ca ca x 1 c b t Y ca ca y 1 c b t Y ca . . . . . . . . . c b t X ca ca x * c b t X ca ca y * c b t X ca c b t Y ca ca x * c b t Y ca ca y * c b t Y ca        +        c b x 1 -ca x 1 c b y 1 -ca y 1 . . . c b x * -ca x * c b y * -ca y *        (18)
where superscript + denotes the pseudo-inverse and

ca η = [ ca n X ca d ca n Y ca d ca n Z ca d ] .
Then the parameters of the planar object w.r.t. the camera c a can be obtained from ca η following:

ca d = 1 ca η ( 19
) ca n = ca d ca η. (20)

B. Non-linear optimization

To estimate ca η more robustly and to allow the method to be extended to relative pose estimation between two LFs (Sec. V), we formulate its estimation as a non-linear optimization of which the initial guess is provided by the linear estimation (Sec. IV-A).

As c b H ca in (15) depends on ca η, we rewrite it as c b H ca ( ca η). Considering Q corresponding points, the optimization problem to solve is:

ca η = min ca η Q q = 0 || c b p * q c b pq ||, (21)
where the operator denotes the difference after resetting the homogenous coordinate at 1 and with:

c b p * q = c b H ca ( ca η) ca pq , (22)
c• pq being the measured point coordinates in the image of SAC c • and c• p * q a computed one. A gradient descent-like algorithm (Gauss-Newton in this paper) solves (21). It requires to express the Jacobian matrix of image points with respect to plane parameters, i.e.:

J F1,F2 = ∂ F2 p * ∂ F1 η , ( 23
)
where the index q is omitted for compactness. F 1 and F 2 denote generic coordinate systems of the homography relationship:

F2 p * ∝ F2 H F1 F1 p (24)
similarly to [START_REF] Hahne | Baseline and triangulation geometry in a standard plenoptic camera[END_REF]. Thus, the detailed generic expression of J F1,F2 is (with F 1 omitted, again for compactness):

J F1,F2 = 1 (R 3 -η t Z )p 2 (R 1 t Z -R 3 t X )p p (R 2 t Z -R 3 t Y )p p , (25)
with R • , the •-th row of the rotation matrix between F 1 and F 2 , and t X , t Y and t Z the translation vector components, between F 1 and F 2 as well.

Considering only two SACs c a and c b of the same LF, J F1,F2 is simplified as J ca, c b :

J ca, c b = - c b t X ca ca p c b t Y ca ca p , (26)
due to the same SACs orientation and the coplanarity of their optical centers [START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF]. As much J ca, c b are computed as there are point pairs between SAIs of c a and c b . These Jacobians are stacked as matrix J, considered in the plane parameters increment ca η computation of an optimization iteration:

ca η = -J + e, ( 27
)
where e is the error vector stacking every c b p * q -c b p q , for each point pair q. Then, ca η updates ca η by addition, allowing to compute a new c b H ca and new c b p * q values, for all q, and this process loops until convergence.

C. The best SACs pair to estimate plane parameters

As introduced in Section I, classical stereo-vision would have, by definition, only two views at once of the same scene, or planar object, to be considered in the estimation of plane parameters. The LF camera, by providing many views of the scene at once, offers a new possibility: choosing which pair of SACs to consider in two views-based computer vision algorithms. Two criteria are considered below to make that choice.

First, as well known, the widest baseline between SACs, however allowing feature points matching, should be considered to solve [START_REF] Bay | Speeded-up robust features (surf)[END_REF] and (21) in order to reduce uncertainties due to slight imprecisions of φ values allowing to get reliable estimates of the plane parameters.

Furthermore, one should also avoid solving [START_REF] Bay | Speeded-up robust features (surf)[END_REF] and (21) if φ's are collinear, or close to be, with the epipoles of both SAIs [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. The key idea is, then, to choose a pair of SACs for which the set of φ pairs, to be considered in [START_REF] Bay | Speeded-up robust features (surf)[END_REF] and (21), is far from being collinear with SACs epipoles.

Combining both baseline and non-epipolar collinearity constraints leads to following best SACs pair selection strategy to get reliable plane parameters estimations:

• if the center of gravity of the SAIs φ pairs is in the topleft corner of the LF camera field of view, then choose top-right and bottom-left SACs to solve ( 18) and (21) (similar way for the three other field of view corners) • if the center of gravity of the SAIs φ pairs is in the top-center of the LF camera field of view, then choose bottom-right and bottom-left SACs to solve ( 18) and (21) (similar way for the bottom-center, left-middle and rightmiddle areas of the LF camera field of view) The latter selection scheme is experimentally validated in Section VI-B1 and estimations compared to the greedy approach of stacking every available SAC.

V. PLANAR OBJECT POSE ESTIMATION FROM TWO LFS

LF 1 and LF 2 are two light-fields. A SAC c a of LF 1 is written c a1 and we write similarly for LF 2 as well as c b .

The goal is to exploit the planar geometry of the observed object to compute the frame change between two of its poses or, equivalently, ca 2 M ca 1 , the rigid transformation matrix between F ca 1 and F ca 2 . The constant frame change ca M 0 from SAC frame F ca to the LF frame F 0 easily leads to the frame change from F LF1 to F LF2 .

We consider c a1 as the reference frame, without loss of generality since ca M 0 is constant, so the plane parameters are estimated as ca 1 η. Thus, LF 1 only contributes to the computation of ca 1 η whereas LF 2 contributes to both ca 1 η and ca 2 M ca 1 . This is written by extending (21) and considering m, a vector representation of ca 2 M ca 1 , to:

ca 1 η m = min ca 1 η, m q 1 C q ( ca 1 η) + 2 C q ( ca 1 η, m), ( 28
)
where

1 C q ( ca 1 η) = || c b 1 H ca 1 ( 1 η) ca 1 p * q c b 1 pq ||, (29)
and [START_REF] Knight | Light field data acquisition devices, and methods of using and manufacturing same[END_REF]. In (30), we introduce SACs c e and c f of LF 2 in order to get enough generality in the mathematical writings to be able to consider different pairs of SACs between LFs. Indeed, as choosing the best SACs pair leads to more reliable results (Sec. IV-C), it is also a key idea in the planar object pose estimation. Thus, ce 2 H ca 1 ( ca 1 η, m) is computed similarly to c b 1 H ca 1 ( ca 1 η) but by, first, changing the frame of ca 1 η to F ca 2 , then to F ce 2 , using the constant frame change ce M ca , and, second, considering the varying frame change ca 2 M ca 1 . The latter's translation and rotation (axis-angle representation) elements are merged together

2 C q ( ca 1 η, m)= || cc 2 H ca 1 ( ca 1 η, m) ca 1 p * q ce 2 pq || + || c d 2 H ca 1 ( ca 1 η, m) ca 1 p * q c f 2 pq || . (30) In (29), c b 1 H ca 1 ( ca 1 η) is computed using
as m = ca 2 t X ca 1 , ca 2 t Y ca 1 , ca 2 t Z ca 1 , ca 2 θ X ca 1 , ca 2 θ Y ca 1 , ca 2 θ Z ca 1 . c f 2 H ca 1 ( ca 1 η
, m) follows similarly. Note that if both LFs have the same best SACs pair at some instant, ce 2 H ca 1 reduces to ca 2 H ca 1 and c f 2 H ca 1 reduces to c b 2 H ca 1 .

To optimize simultaneously both ca 1 η and m, following the same methodology that leads to (27), one must compute increments:

ca 1 η ṁ = -J 2 + e 2 , (31)
used for updating ca 1 η and m in a looping process until convergence. ca 1 η updates ca 1 η by addition and ṁ updates m by composition (see the use of the exponential map of se(3) from the Lie Algebra in [START_REF] Ma | An invitation to 3D vision[END_REF] for more details). In (31), e 2 is the stacking of differences between computed and measured corresponding point coordinates, similarly to (27), but for three pairs of SAC images. Jacobian J 2 first three columns contain the stacking of Jacobians J F1,F2 ( see ( 23)), computed for SAC pairs (c a1 ; c b1 ) (thus computed exactly as in (26), i.e. J ca 1 , c b 1 = J ca, c b , for ca 1 p = ca p), (c a1 ; c e2 ) and (c a1 ; c f2 ), in the general case about pairs of SACs. Both latter are computed using (25), since, in these pairs, SACs belong to different LFs so R • -s and t • -s may have any values. Juxtaposing the stacking of interaction matrices Lc• p * q , relating c• ṗ * q to ṁ [2], to J 2 finishes its building:

J 2 =   J ca 1 , c b 1 0 J ca 1 , ce 2 L ce 2 p * ce V ca J ca 1 , c f 2 L c f 2 p * c f V ca   , ( 32
)
where Lc• p * is, itself, the stacking of every Lc• p * q . ce V ca , the twist transformation matrix, directly computed from ce M ca , to locally transform L ce 2 p * to F ca 2 , so that it contributes to the ca2 M ca 1 frame change computation (same for c f V ca ).

VI. EXPERIMENTS

The proposed linear and non-linear plane parameters estimations (Section IV) have been validated on simulation data before being applied on actual LF acquisitions. The estimation of a planar object pose transformation between two LFs (Section V) has similarly been evaluated. Due to lack of space, only real experiments are presented in the following.

First of all, the LF sequences acquisition setup is detailed. Then, we study the influence of the selected pairs of SACs on the plane parameters estimation accuracy to validate the selection criteria proposed in this paper. Finally, the proposed planar object pose estimation is used for tracking a textured plane.

A. Data acquisition setup

The LF camera considered in our experiments is a Lytro Photo (first generation) camera. The LF camera is static and faces a planar surface moved by the end-effector of a Stäubli TX-60 6-axis robotic arm (Fig. 3).

First, a checkerboard of 19×19 squares of 3.6mm side is set on the planar surface of frame F O . Twelve LF s for s = 1...12 are acquired, varying the pose of the checkerboard (Fig. 4). Each checkerboard pose is measured w.r.t the robot base, F B , by the TX-60 internal software to precisely know the 12 B M Os (for s = 1...12) transformations. These 12 LF are set as the input of the Dansereau's MATLAB Toolbox [START_REF] Dansereau | Calibration and rectification for lenselet-based plenoptic cameras[END_REF] to calibrate the Lytro camera. As output, the calibration gave us, inter alia, the 5×5 intrinsics matrix (3), rectification parameters and the 12 poses of the checkerboard w.r.t. the LF camera: LF M Os for s = 1...12. We used the 12 B M Os and the 12 LF M Os to estimate the fixed pose of the camera w.r.t. the robot base: B M LF 2 . Therefore it is now straightforward to obtain a ground truth of the pose for any planar object put on the robot Fig. 3: Acquisition setup -A Lytro Photo LF camera facing a planar object that is moved by a Stäubli TX60 robotic arm. end-effector w.r.t the camera to evaluate the accuracy of our method.

The calibration toolbox output is also the decoded 4-D LF data structure (Section II). After re-organization, 9 × 9 SAIs of size 380 × 380 pixels are obtained from one LF. In other words, the Lytro camera can be considered as a grid of 9 × 9 SACs. However, due to the undistortion, the outer SAIs may be cropped or blurred. For this reason, during the following experiments, we are only considering SAIs that belong to the 5 × 5 centered SACs of the grid. More precisely, the top-left usable SAI is noted as SAI (3×3) , the top-right usable SAI is SAI (7×3) , the bottom-left is SAI (3×7) , the bottom-right is SAI (7×7) and so on.

B. Best pairs of SACs study 1) Selection strategy evaluation:

The purpose of this study is to experimentally validate the best SACs pair selection strategy introduced in Section IV-C.

For this experiment, we consider the 4 LFs shown in Fig. 4. We note these LFs: LF 1 , LF 2 , LF 3 and LF 4 . As it can be seen on these 4 SAIs, the visual features (the checkerboard corners) are located in a specific part of the field of view in each LF, respectively in the top-left, top-right, bottom-right and bottom-left. These 4 LFs are then interesting candidates to validate our best SACs pair selection strategy that combines baseline and epipolar non-collinearity constraints.

We estimate the plane parameters of the 4 checkerboards following the approach described in Section IV and using different pairs of SACs: {SAC (7×3) , SAC (3×7) } and {SAC (3×3) , SAC (7×7) }. Recall that these two pairs of SACs correspond to the two widest baselines of the usable SACs. TABLE I contains the results of this experimentation. We would like to remind that the maximum baseline between usable SACs is equal to 3mm. The best pairs of SAC for LF 1 , LF 2 , LF 3 and LF 4 are highlighted in TABLE I. To quantitatively evaluate the estimation we compute the mean distance between the ground truth plane and the estimated one. More precisely, this metric corresponds to the average distance between a grid of corresponding points equally distributed on the two planes. We can see that there is a link between the pair of SACs used in the estimation of the plane 2) The best SACs pair versus 25 SACs: As previously mentioned (Section IV), more than two SACs may be used in our plane parameters estimation scheme. Indeed, using the 25 available SACs seems, intuitively, interesting. Any SAC may be used as a reference, thus there are 25 combinations. We have performed the same experiment as before on the same 4 LFs (Fig. 4) but using these 25 combinations of 25 SACs. TABLE II contains for each LF, the best estimation among the 25 and the SAC that has been used as a reference to obtain this result.

SAC (7×3) SAC (3×7) SAC (3×7) SAC (7×3) SAC (3×3) SAC (7×7)
First, it is interesting to note that, apart from the LF 2 , the SAC used as a reference that provides the best result is the same (or is very close, as for LF 1 ) to the SAC selected following our proposed best SACs pair selection strategy (TABLE I). Considering the 25 SACs in LF 2 not only led to a different reference SAC than for our proposed selection scheme but also to a lower estimation precision. Thus, it shows that using every available SAC does not always provide a better estimation result than using a pair, i.e. the best pair under the criteria proposed in this paper.

For the other considered LFs, thanks to the redundancy of information the estimations are more accurate than using the best pairs of SACs (TABLE I). However, the computation using 25 SACs is very time-consuming. An estimation using 25 SACs is performed in approximatively 9.0s when it takes 0.7s using a single pair (MATLAB implementation running on a laptop with an Intel i7-4900MQ CPU with 16GB RAM). So, our SAC pair selection strategy shows being an interesting tradeoff between precision and computation time.

C. Evaluation on textured planar object acquisitions

For the following experiments, we changed the planar object held by the robot (Fig. 3b) and acquired several LF sequences. The Lytro camera is not able to record a video. To overcome this issue, we capture LF sequences by taking a series of LF, slightly moving the planar object between each shots using the robotic arm. Fig. 5 shows some central SAIs of LFs extracted from two acquired sequences.

As for the first experiment, we need to build a set of light ray correspondences between two pairs of SAIs extracted from two consecutive LFs to estimate the displacement of the planar object. However, the visual features detection and matching are more difficult in the case of a textured planar scene rather than a checkerboard. In other words, we must be able to find visual features that each of the four SAIs (two SAIs in the first LF and two others in the second) have in common.

To perform this, we detect and match SURF features [START_REF] Bay | Speeded-up robust features (surf)[END_REF] in the four SAIs. The matching is done by setting an experimentally found hard threshold to avoid outliers while keeping enough matches (5 minimum and well balanced on the object surface) to simultaneously compute the plane parameters and its pose. The data provided by this matching has just the purpose to provide feature points with realistic noise and location uncertainties to look forward for robotics applications. True, more optimal methods of features matching or tracking could be used but the current experiment focuses on multi-LF geometry estimation, not about image processing.

The first pose of the acquired planar object is supposed to be known. Then, our method estimates its displacement, shot after shot, incrementally, and is compared to the ground truth. This experiment evaluates the accuracy of this planar object incremental pose computation depending on the pairs of SACs considered for each LF. More precisely, we compare the estimated poses under the five following SACs configurations:

• D1 -the same pair of SACs for every LF: For more generality of results, we consider two LF sequences (Fig. 5) on which the latter five settings are considered for the planar object incremental pose estimation. The two sequences are composed of nine LFs each. During the acquisition of both sequences the plane has been displaced on approximatively 20cm. In the first sequence, the plane follows pure translations. In the second sequence, translations and rotations are applied.

{SAC (3×3) , SAC (
To take into account the estimation drift, inherent to any incremental estimation, the last planar object estimated pose of each path is compared with the ground truth (TABLE III).

Estimation results of both experiments clearly show that the pose estimation of the planar object is highly more accurate when taking advantage of the LF camera acquisition structure, as we proposed, than considering a constant horizontal pair of cameras as in classical stereo-vision: precision improvement of up to an order of 14 for the first angle of DHV, w.r.t. the first angle of H, in sequence 2, for instance. The comparisons between the estimated paths and the ground truth (Fig. 6) also highlight these observations. To compare with a LF-based approach from the state-ofthe-art, Figure 6 also shows an additional estimate of the planar object path (path "J") using [START_REF] Johannsen | On linear structure from motion for light field cameras[END_REF]'s Matlab scripts 3 . The latter Structure-From-Motion approach, while producing precise results in non-planar scenes, produces poor results when considering a planar one due to its essential matrix-based modeling [5, Sec. 11.9.2, p. 296]. Our approach outperforms the latter, thanks to its consideration of the fact that the observed scene is planar.

VII. CONCLUSION

In this paper, we have proposed a planar object pose estimation from light rays captured by a LF camera. To develop this method, we have first expressed the 5 × 5 homography between corresponding rays of a LF. The LF camera geometry and the 5 × 5 homography have been considered to implement a plane parameters estimation from one LF. A linear and a more robust estimation based on non-linear optimization have been proposed, considering the best pair of SACs. Based on all of this, we have introduced a method to estimate the rigid transformation of a planar object between two LFs. Results show that the combination of these contributions leads to the reliable planar object tracking over time from multiple LF acquisitions.

  how light rays emitted by a 3-D point P = [X, Y, Z] are observed by the camera. The local two-plane parameterization (L2PP) describes a light ray as its intersections with two parallel planes, i.e. the focal plane Π and the image plane. The latter is considered the closest to the 3D point and the focal plane the farthest. The LF camera coordinates frame F 0 is attached to the focal plane. Then, o = [m, w] are the coordinates, expressed in the global frame F 0 , of the intersection between the ray and the focal plane. The intersection of the light ray with the image plane gives two more coordinates, p = [x, y] , but expressed relatively to o. Stacking both coordinates leads to the ray:

Fig. 1 :

 1 Fig. 1: LF camera acquisition: a complete raw LF image (a) and a zoomed subpart of this LF (b).

1

  http://dgd.vision/Tools/LFToolbox (a) (b)

Fig. 2 :

 2 Fig. 2: (a) Decoded data structure (4-D). (b) SAIs (each bloc of figures from 0 to 24 is a SAI) generated from the decoded 4-D data structure. This figure is a rewriting of Fig. 6 in [8].

  c b H ca only relate the third to fifth coordinates of φ's. Their first and second coordinates, 0 o ca and 0 o c b , correspond to optical centers coordinates of SACs in the LF camera frame F 0 . Since, optical centers of SACs belong to the same plane (Remark 1), and assuming optical axes of SACs are parallel (assumption discussed in Remark 2, below), the transformation between two SACs is a translation of length c b t X ca and c b t Y ca along the two axes of the focal plane.

Fig. 4 :

 4 Fig. 4: Central SAI of relevant LF to validate our best SACs pair selection strategy.

  7×7) } (first diagonal baseline) • D2 -the same pair of SACs for every LF: {SAC (7×3) , SAC (3×7) } (second diagonal baseline) • H -the same pair of SACs for every LF: {SAC (3×3) , SAC (7×3) } (horizontal baseline) • V -the same pair of SACs for every LF: {SAC (7×3) , SAC (7×7) } (vertical baseline) • DHV -the best pair of SACs for each LF following the strategy selection proposed in Section IV-C.

Fig. 5 :

 5 Fig. 5: Some central SAIs that compose two LF sequences.

  5 H 6 H 7 H 8 H 9 H 10 H 11 H 12 H 13 H 14 H 15 H 16 H 17 H 18 H 19 H 20 H 21 H 22 H 23 H 24 H 25

TABLE I :

 I Plane estimation errors [mean distance in mm] for various LFs using the four pairs of SACs of widest baselines. Bold SAC is the reference one.parameters and the position of the visual features in the LF. The best pairs of SACs for LF 1 , LF 2 , LF 3 , and LF 4 are respectively {SAC (7×3) , SAC (3×7) }, {SAC (3×3) , SAC (7×7) }, {SAC (7×3) , SAC (3×7) }, and {SAC (3×3) , SAC (7×7) }. These best SACs pairs correspond to the misalignment of the checkerboard center of gravity in the LF camera field of view w.r.t. both SACs, as awaited. Thus, these results validate the best SACs pair selection strategy proposed in Section IV-C.

• The expression of the homography between corresponding rays of a LF. • The exploitation of the LF camera geometry to simplify the homography content between two SAIs of a LF, used for plane parameters estimation.

LF images and robot poses are available in the LFMIS dataset at: mis.u-picardie.fr/∼g-caron/pub/data/LFMIS dataset.zip (a) (b)


David Guyon
email: david.guyon@irisa.fr

Anne-Cécile Orgerie
email: anne-cecile.orgerie@irisa.fr

Christine Morin
email: christine.morin@inria.fr

Energy-Efficient IaaS-PaaS Co-design for Flexible Cloud Deployment of Scientific Applications

Keywords: Cloud computing, Infrastructure-as-a-Service, Platform-as-a-Service, energy

Reducing the massive amount of energy consumed by cloud datacenters becomes of major importance. In the usual approach where resources are consolidated into fewer servers in order to power down the others, it still remains periods of time when servers are not fully utilized. Consequently, it exists unused resources that are not exploited although they could be used to execute applications compatible with the variable availability of these resources. In this work, we propose a cloud system where the Platform-as-a-Service (PaaS) and Infrastructure-as-a-Service (IaaS) layers interact to find execution trade-offs that exploit the unused resources at IaaS level. PaaS users are involved in the energy optimization by proposing to delay their executions and adapt resource sizes in order to fit with the available unused resources. Our evaluation by simulation is based on real data and expresses a realistic large scale cloud scenario. Results show that according to the proportion of energy-aware users, this system is able to reduce the amount of servers by using resources that would have been wasted otherwise. Therefore, our solution allows datacenters to consume less energy than with usual resource managers where all applications start their execution at submission time with their initial resource size.

I. INTRODUCTION

Estimations show that by 2030 the energy consumed by datacenters will represent 13% of the worldwide energy consumption [START_REF] Andrae | On Global Electricity Usage of Communication Technology: Trends to 2030[END_REF]. Due to climate change and the increasing success of Internet services, priorities should be given to the reduction of their energy consumption.

In cloud computing, studies have been conducted to increase datacenters energy efficiency [START_REF] Orgerie | A Survey on Techniques for Improving the Energy Efficiency of Large-Scale Distributed Systems[END_REF]. Their efforts focus on optimizing energy consumption at Infrastructureas-a-Service (IaaS) level, as it is closer to hardware resources than the Platform-as-a-Service (PaaS) layer, and thus closer to the electricity consumption itself. However, several studies already showed that higher energy savings are possible when enhancing interactions between these two cloud layers [START_REF] Carpen-Amarie | Towards energy-aware IaaS-PaaS co-design[END_REF], [START_REF] Dupont | Energy Efficient Data Centres Within Smart Cities: Iaas and Paas Optimizations[END_REF]. IaaS knows about the availability of hardware resources and can deliver energy-related information that could help the PaaS layer to make energyaware decisions. In return the PaaS layer could inform IaaS providers on users' applications flexibility in order to help the consolidation process. Typically, when PaaS users ask for virtual resources, requested resources are made available as soon as possible. However, some users could accept their request to be handled differently if it saves energy by either delaying the deployment or changing the size of the requested resources (and consequently the duration). While this approach is not compatible with applications that continuously execute (web jobs), time-bound scientific applications could exhibit flexibility on starting time and resource size as long as results arrive before a deadline and if total cost does not increase. This scenario is realistic as scientists running HPC applications are looking at clouds as a cost effective alternative to HPC [START_REF] Netto | HPC Cloud for Scientific and Business Applications: Taxonomy, Vision, and Research Challenges[END_REF].

Our objective is to reduce IaaS datacenter energy consumption with a cooperation allowing PaaS to express the flexibility of its applications and IaaS to inform on when and how many resources are predicted to be unused. This way, energy savings are achievable by shifting and resizing some applications on these otherwise unused resources.

Our proposition is a PaaS-IaaS co-design that offers to PaaS users energy-efficient execution trade-offs. This codesign allows users to decide the flexibility level they agree to give to their execution. Our evaluation by simulation, based on real datacenter traces and real application execution logs, shows that when all users deploy applications that are flexible in time and size, it is easier for IaaS providers to allocate resources to applications without the need to power up additional servers. Up to 5.49% of energy is saved compared to a scenario already consolidating the workload, powering down unused servers and where all users prefer their application execution to start at submission time.

The rest of the article is organized as follows. Section II presents the context and list the assumptions. The system architecture is described in Section III and evaluated by simulation in Section IV. Discussion on limitations is given in Section V. Section VI discusses the related work. Finally, we conclude in Section VII.

II. CONTEXT AND ASSUMPTIONS

In this study, we consider a PaaS cloud with users deploying scientific applications. These applications are CPUintensive, execute in a single Virtual Machine (VM) and have a finite execution time. Applications with a continuous execution (i.e. web application) are not considered here as they cannot be delayed or changed in size without impacting users' satisfaction.

The PaaS negotiates with multiple IaaS providers to provide virtual resources to its users. We assume that PaaS users can choose to be tolerant to delaying their application execution. This flexibility in time has a limit defined in the PaaS Service Level Agreement (SLA) as the time window (T W ). Its value defines the maximum number of hours the execution of an application can be delayed. As application executions can be delayed, we assume that IaaS resources can be reserved in advance (within T W ).

VMs have a fixed size in terms of vCPUs, RAM and disk usage (called flavor), that remains unchanged during their lifetime (no dynamic resource scaling). However, here we simplify the problem by defining flavors only by their number of vCPUs, the most important component for CPUintensive applications.

VMs energy consumption can be estimated in advance and computed more precisely after their execution. For this calculation to be possible, we make the assumption that IaaS providers know the power profile of their servers. IaaS providers optimize reservation placements with an allocation algorithm that targets the use of a minimum number of servers. In addition, each provider tries to reduce its energy consumption by shutting down unused servers [START_REF] Orgerie | A Survey on Techniques for Improving the Energy Efficiency of Large-Scale Distributed Systems[END_REF] and putting them back online when required.

Based on utilization traces and because workloads tend to repeat themselves over time (usually on a daily basis), we assume that IaaS providers can predict the workload of their servers [START_REF] Wamba | Cloud Workload Prediction and Generation Models[END_REF]. The load is defined as the amount of reserved vCPUs per server. Additionally, the power state of each server is also known. In other words, it is possible to record when servers are online or not (powered up/down) during the day, thus providing information on when resources are active or not. With this, it is possible to compute the perserver Resource Availability (RA). This value represents the amount of unbooked vCPUs when servers are powered up. As shown in Figure 1, RA fluctuates from 0 to the vCPU capacity of the server according to the arrival and departure of resource reservations. When RA reaches 0, it is either because the server is powered down or because all vCPUs are allocated. When all vCPUs are available, it means the server is not used. In this case, the server is powered down, thus making resources unavailable (RA falls back to 0 at the end of Figure 1).

Additionally, we make the following assumptions: • users know the execution time of their applications for each VM flavor. Typically, scientific applications are CPU-intensive and executed multiple times, so users can conveniently know their execution duration depending on the vCPUs number and the data input size. Moreover, predicting application's execution time for varying number of vCPUs in the VM configuration has been shown to be possible [START_REF] Li | Application Execution Time Prediction for Effective CPU Provisioning in Virtualization Environment[END_REF] ; • each IaaS provider has its own incoming workload in addition to the one induced by the studied PaaS ; • virtual resources are given for a specific duration and are released at the end of this period ; • requested virtual resources are considered to be fully booked, regardless of whether the application really use them or not. CPU overcommitment is not considered.

III. PROPOSITION A. General Idea

In this work, we propose to benefit from the exchange of information between PaaS and IaaS layers, as well as the users' willingness to change resource size and to delay their application execution in order to achieve energy savings.

The general idea to achieve energy savings is presented in Figure 2. Based on the RA information, IaaS providers are able to propose to their users execution contracts (EC) with a delayed execution in order to fit with the availability of resources. An EC corresponds to a spatial and temporal placement of a VM on servers that stays within the maximum T W hours tolerated delay. Orange rectangles in Figure 2 depict three different ECs consuming different amounts of energy. Two servers owned by an IaaS provider host a total of three VMs. Server 1 is able to host all VMs, which explains why server 2 remains powered down. A user desires to execute her 8 vCPUs application for a duration of 30 minutes. When the VM request occurs, the only way to execute the application directly is to power up server

(App[C1]

). This is the most energy consuming contract. Delaying the execution of the application of about 1 hour avoids to use server 2 (App[C2]) and reduces the energy consumption. It is also possible to execute the application at submission time without powering up the second server by changing the size to 4 vCPUs (App[C3]). In this simple use-case, we consider that attributing half of the requested vCPUs doubles the execution time of the application. This third contract has the same energy cost as the second contract but execution results are available sooner.

B. System Architecture Overview

In this section, we present the architecture of our system. Details on each component are provided in the following sections. The complete system architecture is presented in Figure 3. Each step of the overall process given in the following list is linked to the labels attached to the figure:

1) a user sends a request to the PaaS provider to execute her application ; 2) the PaaS provider sends several requests adjusted to the user flexibility to multiple IaaS providers ; 3) each IaaS provider proposes an EC for each request it received ; 4) the PaaS provider filters the ECs based on their energy cost and delay in order to propose three contracts to the user ; 5) the user selects the contract she wants to execute her application with ; 6) the PaaS provider informs the corresponding IaaS provider that its contract has been selected ; 7) this IaaS provider plans the execution of the application in a VM as defined in the EC. PaaS users send the following details: the application, the size S in number of vCPUs (within the available flavors listed in Table I), and a table of execution lengths according to each possible size (D[S]).

In a previous work [START_REF] Guyon | Energy-Efficient User-Oriented Cloud Elasticity for Data-Driven Applications[END_REF], we already showed that proposing to IaaS users to variate their resource size helps to reduce datacenters energy consumption. Here, the PaaS provider acts like an IaaS user and sends 4 different execution requirement (ER) requests to its collection of known IaaS providers:

• asks for a VM with S vCPUs (f lavor base ) for a duration of D[f lavor base .cpu] seconds without delay (as soon as the application submission arrives) ; • asks for the same size and duration but allows a maximum delay of T W hours ; • asks for the flavor one step larger than f lavor base for a duration of D[f lavor base+1 .cpu] seconds within a maximum of T W hours delay ; • same as previous requirement but with a flavor one step smaller than f lavor base and a duration of

D[f lavor base-1 .cpu] seconds.
This mechanism, handled by the trade-off size/duration component of our architecture, is based on the duration table given by the user. Whenever an IaaS provider receives an ER request, its consolidation component searches in its collection of servers for the most energy efficient EC. If no delay is tolerated, it searches for a server with just enough available vCPUs at the current time. This is a bin packing problem. We use a greedy algorithm that gives priority to servers that are already active. The scheduling is optimized in order to have the highest amount of reserved vCPUs on the server. It avoids small VMs to "pollute" servers with large free spaces that could be used by larger VMs. If there are not enough available vCPUs on active servers, the system powers up a new one.

If delay is tolerated, the system uses the resource availability component that allows to know the availability of unused vCPUs in the next T W hours. It provides the delay value of a possible time slot for a VM of the requested size that stays within the tolerated deadline. If no time slot can be found, the process informs the PaaS layer and falls back to a scheduling where no delay is tolerated. More details on the resource availability are given in Section III-C.

Each EC comes with an estimation of the VM energy cost. This estimation is based on the power profile of servers, the VM size and its duration. Details on this calculation are given in Section III-D.

The consolidation component may find several ECs for a single ER request. Instead of sending them all to the PaaS layer, each IaaS provider executes a first filtering round to only send the most relevant. Contracts are projected in a graph plotting the energy prediction for the various delay. Data normalization is used in order to have axes with equivalent weight. We select the contract with the shortest Euclidean distance with the (0, 0) coordinates because this is the one with an equally balanced trade-off between delay and estimated energy cost.

To summarize, IaaS providers include in their EC the starting time t start that is between t submission and t submission + T W and the energy E V M that the VM is estimated to consume.

For each ER request sent by the PaaS, IaaS provider that has been contacted submits its proposition of EC. At this stage, the contract filtering component selects the most energy-efficient EC received for a given ER. It searches for the contract that is the most equally balanced between delay and energy consumption. This filtering feature is the same as the one used for the first filtering round executed at the IaaS layer (trade-off between delay and energy).

Finally, end-users are given a choice between three ECs where they can choose the one they prefer: C1: application execution starts at t submission without changing the requested amount of resources ; C2: application execution is delayed of up to T W hours and the requested resources are not changed ; C3: application execution is delayed of up to T W hours and the given resources are one size larger or smaller than originally requested. Each EC displays its starting and completion time (t start + D). Delayed contracts (C2 and C3) display the percentage of energy they save in comparison with the contract starting directly (C1). A percentage of energy saved is provided instead of real values given by IaaS providers because revealing real values could break the separation of concerns between cloud layers. In the case of C3, it is also required to show the new resource configuration.

Each EC is valid for acceptance during a certain amount of time (few minutes) defined in the IaaS SLA as the user's response time. The PaaS user has to select the contract that suits the most her need within this time. All contracts that are not accepted within the allowed response time are released. Her choice is made in terms of execution time and depending on her motivation in reducing her impact on energy consumption. When an EC is accepted, the PaaS layer informs the corresponding IaaS provider to schedule the VM creation. If the selected contract is still valid, the IaaS provider executes two actions: the VM creation is scheduled at time t start as specified in the EC and resources taken by the VM are removed from the RA schedule (resources pre-allocation). Updating RA is required because the resources that have just been pre-allocated will not be available when the VM will be created. In the case where the selected contract is not valid anymore, the PaaS has to restart its negotiation.

At time t start , the VM creation component creates the VM on the server given by the consolidation component.

Next, we detail each component of our architecture.

C. Resource Availability Management

The resource availability informs on when the vCPUs of a server are expected to be available over time based on its historical workload. The computation of RA is based on the IaaS incoming workload excluding the PaaS workload.

Each IaaS provider generates a RA prediction for each server it owns. A per-server RA differs from a RA at the datacenter's level because it allows to know the largest VM a server can host. As an example, with an aggregated RA at infrastructure level, when 10 vCPUs are predicted to be available, it is impossible to know if it is 5 servers with 2 available vCPUs each or a single server with 10 unused vCPUs.

Each time a VM is instantiated or removed from a server, the new amount of vCPUs that are available on this server is recorded. We define n s as the total number of vCPUs on a server s and the function usage s (t) that returns its number of vCPUs reserved at time t. The RA function can then be defined as:

RA(t) = n s -usage s (t).
The consolidation algorithm goes through the RA of each server in order to find a time slot for a requested VM. A valid time slot corresponds to a period of time in a server's RA when at least S vCPUs are available for a duration of minimum D[S] seconds. At first, the algorithm searches on each server for time slots with an availability duration that is the closest to the application duration. Then, following a greedy algorithm, in this sub-list of possible time slots, only the one with the smallest difference between the amount of available vCPUs and the requested size S is retained. This time slot, along with its energy estimation cost (described in Section III-D) are included in the EC returned to the PaaS.

An IaaS provider has to update the corresponding server's RA when one of its proposed contract is selected, because of resources pre-allocation. The amount of vCPUs at the time slot specified in the contract is subtracted from the vCPUs that were predicted to be available. The RA schedule evolves over time as contracts are accepted.

D. Execution Contract Energy Estimation

In order to let the PaaS layer know which contract is the less energy consuming, each EC includes an estimation of the VM energy cost.

Our VM energy estimation is the sum of two parts. The first one is a portion of the idle energy consumption of the server. This energy is shared with the VMs currently hosted on the server according to their size. The second energy value is the maximum dynamic energy that the VM can consume. The dynamic energy is based on the application execution duration, the requested number of vCPUs and the power consumption of a single vCPU. We consider this CPU-based power model because the processor is known to be the most energy consuming component of servers [START_REF] Fan | Power Provisioning for a Warehouse-sized Computer[END_REF] and this work targets CPU-intensive applications. If a server needs to be powered up for a VM, this energy cost is not included in the VM energy calculation.

This energy estimation requires two power values for a given server s: P s idle and P s CP U . P s idle is the power consumed by a server while not being used (i.e. not hosting any VMs) and P s CP U is the power consumption of 1 vCPU used at its maximum capacity. Power consumption values can be measured with wattmeters attached to the servers. P s idle is reached each time there are no VMs running on the server. The measurement of P s CP U relies on P s max the maximum consumption with all vCPUs. The calculation of the CPU power is defined as follows:

P s CP U = (P s max -P s idle )/n s (1)
Following, the equation to calculate the idle energy consumption of a VM.

E V M idle = P s idle × tstart+D[S] tstart S n s -RA(t) + S dt (2)
The ratio refers to the requested number of vCPUs (S) over the total number of allocated and requested vCPUs. It reaches its closest value to 0 when the smallest VM size runs in a fully used server, and equals to 1 when any VM size is running alone on the server. The integral gives the average ratio over the VM lifetime.

The maximum dynamic energy the VM is expected to consume is calculated with:

E V M dyn = P s CP U × S × D[S]
. Finally, the total amount of joules a VM is estimated to consume corresponds to the sum of idle and dynamic energies:

E V M = E V M idle + E V M dyn
The power estimation component allows to have an estimation of the VM energy consumption. However, this estimation can differ from the energy consumed after execution. As a matter of fact, the ratio used when calculating E V M idle evolves according to the update of the RA(t) function. Therefore, it causes the energy cost to be sometimes overestimated.

IV. EVALUATION

A. Experimental Setup

To validate our approach a simulator has been developed. The simulated PaaS application is the Montage scientific application (detailed in Section IV-A2) that is executed on real VMs of all possible sizes in order to get the different execution time lengths. Our system considers the flavors listed in Table I. It is based on the types of instances that can be found on the Amazon EC2 platform. Finally, workloads used at both IaaS and PaaS levels are traces from real cloud systems (detailed in Section IV-A3). They define the starting time, the size and duration of each VM. 1) Server Model: Simulated servers are based on real servers that we use to run our scientific application. These servers, provided by the Grid'5000 platform [10], are Dell PowerEdge R430 with 16 cores, 64GB of memory and 600GB of disk and equipped with fine-grained wattmeters. They have an average P idle of 75.83W and consume 8.81W when powered down. The originally measured P max was 174.04W. However, IaaS servers are usually of a larger size in order to host an important amount of VMs. To account for this in the simulator, we double the number of CPU cores per server to 32. P max is adapted accordingly to 272.25W. Based on Eq. 1, P CP U is 6.14W. Finally, we also measured energy of on/off state transitions that are respectively 4.92Wh (143 sec) to power up and 0.11Wh (6 sec) to power down. In our evaluation, this server model is used by all IaaS providers.

The number of servers per IaaS, shown in Table II, is defined in order to handle the peak usage of resources defined by their workload.

2) Scientific Application: To load our PaaS cloud with realistic scientific applications, we use Montage [START_REF] Berriman | Montage: A Grid Enabled Engine for Delivering Custom Science-Grade Mosaics on Demand[END_REF] which is CPU-intensive. Montage is a toolkit for assembling images into custom mosaics specialized for astronomical images. By giving a space location and a high/width (in angular degrees), it downloads all the images from its servers that match with the given space area. Then, it detects similarities between all images and finally generates an image of the space area. In a preliminary experiment, we execute Montage in VMs of all possible sizes hosted on a real server. We build a specific Montage VM image optimized to be CPU intensive. This image contains the required input images and has a Montage installation using MPI in order to benefit from all available cores. This is a common PaaS behavior to have several VM images, each of them specialized for specific executions [START_REF] Pierre | ConPaaS: A Platform for Hosting Elastic Cloud Applications[END_REF]. Figure 4 presents the execution time (D[S]) of Montage for each VM size. It also shows the price it would cost to execute it on Amazon EC2 instances in an hourly basis (pay for a full hour, even if reservation is less than an hour) and prorated (only pay for the exact reservation duration). This calculation corresponds to the flavor price given in Table I times the application execution time.

3) Real Datacenter Workloads: We use a real job submission trace as the input workload of each cloud provider. The Parallel Workloads Archive [START_REF]Parallel Workloads Archive[END_REF] provides open source workload logs. The longest utilization log available (2 years long) in the archive is from the MetaCentrum Czech National Grid [START_REF]The MetaCentrum 2 log[END_REF]. In order to find a typical working day (24h long workload) as an input of each cloud provider, a k-mean algorithm is executed to group the days with a similar job submission distribution. An analysis with 6 groups (k-mean with k = 6) exhibits one group with a typical working day trend (load during working hours). Four different days are taken from this group to represent the 4 IaaS workloads used in experiments (listed in Table II). All jobs in workloads are normalized in order to have a size within the available flavor sizes. Also, only the jobs with a duration between 5 min and 12 hours are retained. The duration of each IaaS job is rounded to hours because IaaS clouds, such as Amazon EC2, provide VMs on an hourly basis. For the PaaS layer, another workload is taken from the same k-mean group. This workload contains only 2, 4 and 8 vCPU jobs for easing the simulator task. Indeed, the C3 contract tries to scale up and down the resource size, thus we keep the larger and lower VM sizes for this contract.

B. Experimental Validation of our Approach

We validate our approach by simulation through 4 different experiments. The energy saving is in percent in comparison with the scenario where all PaaS applications execute directly without any modification (C1).

1) Energy Saving According to the Contract Selection: In this first experiment, we show the impact of the contract selection on the IaaS energy consumption and on the total monetary cost of applications. In this context, there is only one IaaS provider (IaaS 1) and the allowed TW is set to 6 hours for all applications. Figure 5 presents the energy consumption (blue bars) and the percentage of energy saving (white curve) when varying the distribution of the contract selection. It also displays the sum of all applications monetary cost (orange bars) as defined in Figure 4. For this experiment, an additional execution contract (C4) is used to define the selection of the less energy consuming contract between the 3 proposed contracts.

Firstly, all applications execute with the same contract (three first bars). In contrast to a normal execution (C1), when all users select C2 or C3, the IaaS provider achieves respectively 2.16% and 5.14% of energy saving. It is important to keep in mind that at datacenter's scale, each percent of energy saved is greater than 1kWh. Moreover, the comparison is made with an IaaS manager that consolidates resources and shuts down unused servers ; a scenario already more energy-efficient than currently deployed datacenter managers. In the case of 100% C2, it is not always possible to find a time slot within T W . Among the 1500 applications, 58 of them, so less than 4% had to run with a C1 contract instead of C2 (i.e. their execution were not delayed). Considering the monetary aspect, costs with C1 and C2 are the same because the resource size does not change. In the case of C3, changing the resource size reduces the total cost by almost 41%. It reveals that resources tend to be reduced in size (less expensive).

Secondly, a percentage of applications executes with the less energy consuming contract among all three (this contract is named C4 for clarity sake) and the remaining applications with contract C1. Increasing the amount of C4 contracts from 0% to 100% shows an increasing percentage of energy saving going up to 5.49%. In the 100% C4 scenario, there is 91.2% of C3 because the high flexibility of C3 shows lower energy consumption in most cases. As C1 and C2 contracts can consume less than C3, 100% C4 saves more energy than when using only C3 contracts. From a monetary point of view, executing all PaaS applications is less expensive when more of them execute with the C4 contract.

Results show that, in this context, users acceptance to delay and change the size of their applications allows to reduce the energy consumed by the IaaS provider and to lower the monetary cost for PaaS users.

2) Impact of the Time Window on Energy Consumption:

In the second experiment, we show the link between the allowed time window T W and the energy savings. Figure 6 presents the results of an IaaS where all applications execute with the C2 contract and T W increases from 1 to 13 hours in 2 hours steps.

It shows that a T W of more than 3 hours is required to reach around 2% of energy savings but remains stable beyond this value. The red line shows the number of C2 contracts that failed to find a delayed time slot. In this context, a T W greater than 7 hours is necessary to have all VMs successfully delayed. This experiment shows that a minimal length of T W is required to save energy, but a too large one does not provide additional savings. II, the T W is fixed to 6 hours, and all applications execute with contract C2. In Figure 7, a cumulative distribution function (CDF) is used to plot the 1500 execution delays for each number of IaaS. The figure shows that whatever the number of IaaS providers, the CDF always reaches 100% before 6h of T W . It means that the system is able to keep execution delays within the time window constraint. It also shows that each time an additional IaaS provider is available, the system is able to execute more applications with a smaller delay. For example, less than 40% of applications have a maximum delay of 1h with 1 IaaS provider whereas 3 IaaS providers cover about 60% of applications for the same delay. In this scenario, 4 IaaS providers allow to have almost all applications with a delay that barely exceeds 3h (50% of the time window limit).

V. DISCUSSIONS Changing the amount of resources is possible as long as the duration table is provided. To obtain this table, users require to execute their applications on all possible VM sizes or to follow the approach presented in [START_REF] Li | Application Execution Time Prediction for Effective CPU Provisioning in Virtualization Environment[END_REF]. A more convenient solution would be to automate this process by learning from previous executions. The PaaS would record all application execution times for each VM size.. With a learning process, the PaaS would then be able to estimate applications execution time depending on resource size. PaaS users wait until completion time to get their results. It corresponds to the sum of the delay and the execution length, both varying depending on the users flexibility. Thus, with our system, users wait for a longer time than usual but our evaluation shows two major benefits. First, it allows users to contribute to the reduction of datacenters energy consumption. And second, flexibility on the resource size tends to lower applications monetary cost.

Each execution contract comes with an energy prediction that can differ from the energy actually consumed after execution. Out of the 1500 applications, only 2.93% of them have a perfect energy prediction with our model. For about 79% of applications, the prediction is over-estimated. Indeed, while waiting for the execution to start, the server might have received other reservations. A higher vCPU reservation reduces the idle energy part attributed to the application's VM (see Eq. 2). The remaining 18% represents applications consuming more than predicted. The reason of this difference is because the continuous VM allocation on servers and the IaaS consolidation algorithm can make some servers to be better for hosting a VM than the one selected during the RA generation. Such a change can make a VM run on a server less used than predicted and thus it accounts from a larger part of the server's idle energy. The difference between predicted and consumed is not an issue because the overall prediction is over-estimated, and as the experiments show, these energy estimations provide significant savings.

The proposed system would not work in reality without a certification of the process. An audit is required in order to avoid IaaS providers from tampering their energy prediction so that they can increase their profit. At the PaaS level, it would force the filtering component to be impartial with contracts from different IaaS providers.

VI. RELATED WORK

To our knowledge, a few studies exist on the topic of IaaS-PaaS co-design for reducing datacenter energy consumption. In their position paper [START_REF] Carpen-Amarie | Towards energy-aware IaaS-PaaS co-design[END_REF], the authors express that an enhanced interaction between the PaaS and the IaaS layers could help achieving energy savings. Dupont et al. [START_REF] Dupont | Energy Efficient Data Centres Within Smart Cities: Iaas and Paas Optimizations[END_REF] present an extension of the traditional SLA that includes a description of the flexibility of applications deployed in cloud environments. This approach allows to have temporarily reduced or increased performance according to the availability of renewable energies. In their approach, the authors make live optimizations, such as resources scaling, while in our work we opt for temporal consolidation of resources that do not change in size over time. Unfortunately, their proposition does not present validation results. In [START_REF] Djemame | PaaS-IaaS Inter-Layer Adaptation in an Energy-Aware Cloud Environment[END_REF], Djemame et al. propose a self-adaptive IaaS-PaaS architecture that allows to run energy-efficient cloud operations. Applications are defined with complex user-defined SLA rules that express constraints on performance, energy and price aspects. The PaaS layer negotiates and compares offers from multiple IaaS providers and selects the best one. These solutions are similar to our proposition as they enhance the relation between IaaS and PaaS layers, and take end-users into account. However, the complex mechanisms proposed to users might prevent them from using such systems. We believe that systems including end-users should be designed to be as simple as possible to use. Additionally, the temporal placement of VMs based on resource usage prediction is not considered in these studies. In a non-virtualized datacenter context, Hsu et al. propose a solution for the waste of resources in datacenters by looking at the waste of power budget [START_REF] Hsu | SmoothOperator: Reducing Power Fragmentation and Improving Power Utilization in Large-scale Datacenters[END_REF]. Their idea is to group jobs with asynchronous peak times under the same power node, thus allowing to reduce the peak power of each server. While their objective is similar to our waste reduction goal, we focus on resource wastage instead of power and we consider two user parameters: resource size and temporal placement.

VII. CONCLUSION

In this work, we propose to benefit from the flexibility in time and size of PaaS applications in order to reduce the energy consumption of underlying IaaS datacenters. Our approach enhances interactions between PaaS and IaaS and includes end-users in the resource allocation decisions. The PaaS layer negotiates execution contracts with multiple IaaS providers for executing users' applications. End-users are then proposed to select between 3 different execution contracts that allow to save more or less energy.

Through our simulation-based evaluation, we validated the gain in energy of our approach using real workload traces and real application execution logs. Allowing PaaS applications to be flexible in terms of execution delay and size allows to save up to 5.49% of energy at datacenters' level in comparison with an already energy-efficient cloud scenario where VMs are consolidated and unused servers are powered down. This energy saving is achievable by only delaying application executions up to a maximum of 6h.

As future work, we would like to enhance the inclusion of PaaS users by providing an additional flexibility in the means of dynamic scaling up and down of virtual resources. Similarly to [START_REF] Dupont | Energy Efficient Data Centres Within Smart Cities: Iaas and Paas Optimizations[END_REF] and the lever we proposed in [START_REF] Guyon | Energy-Efficient User-Oriented Cloud Elasticity for Data-Driven Applications[END_REF], users would be able to set triggers in their application execution to add/remove resources on the fly in order to adapt their size with the current load. These predictable resource updates, motivated with a reward system, may favor improved consolidation and thus increase energy savings.

Figure 1 :

 1 Figure 1: Example of a server's Resource Availability fluctuating according to the arrival and departure of resource reservations.

Figure 2 :

 2 Figure 2: Possible scheduling (execution contracts) of an 8 vCPUs application in an infrastructure with 2 servers and 3 VMs. Starting the application at submission time (C1) requires to turn on server 2. Delaying (C2) or changing the size (C3) can avoid the need of server 2, thus saving energy.

Figure 3 :

 3 Figure 3: Detailed system architecture with the components of both cloud layers and the interactions between them and the end-user.

Figure 4 :

 4 Figure 4: Execution time, hourly price and prorated price of Montage in all possible VM sizes.

Figure 5 :

 5 Figure 5: Energy saving percentage and total monetary cost of all PaaS application according to the distribution of the contract selection made by PaaS users.

Figure 6 :

 6 Figure 6: Percentage of energy saving with C2 contracts according to the allowed time window.

Figure 7 :

 7 Figure 7: CDF of application delays according to the number of IaaS providers when all applications execute with contract C2.

Table I :

 I List of VM flavors considered in the simulator with their Amazon EC2 hourly pricing for the EU Paris region.

	Flavor	vCPU (S)	EC2 cost
	t2.small	1 $0.0264/h
	c5.large	2	$0.101/h
	c5.xlarge	4	$0.202/h
	c5.2xlarge	8	$0.404/h
	c5.4xlarge	16	$0.808/h

Table II :

 II Details of the cloud layers.

	Cloud layer Nb of servers	Workload size
	IaaS 1	60	2858 VMs
	IaaS 2	60	2273 VMs
	IaaS 3	250	3641 VMs
	IaaS 4	40	1429 VMs
	PaaS	N/A 1500 applications

ACKNOWLEDGMENTS

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).