Georges Beurier
email: georges.beurier@ifsttar.fr

Xiaolin Yao

Yoann Lafon

Xuguang Wang
email: xuguang.wang@ifsttar.fr

A Markerless Method for Personalizing a Digital Human Model from a 3D Body Surface Scan

Keywords: Digital Human Models, Template, non-rigid registration, bodyscan

Digital Human Models (DHM) are used for ergonomic design of products. For instance, vehicle ingress/egress motions are simulated for assessing vehicle accessibility. In order to validate simulations, experiments are often needed implying motion capture and motion reconstruction using a DHM. The first step for motion reconstruction is to create a personalized DHM respecting the anthropometric dimensions of the volunteer performing the task. However creating a personalized DHM from external body shape is not straight forward, because the internal skeleton has to be identified from external body shape. Here we propose a four-step method for generating a personalized DHM which matches a 3D scan. The first step is to clean the scan data and to prepare a DHM and a third body surface template. Then, thanks to the use of the third common body template, the correspondence between the DHM and scan surface points is established, making it possible to calculate the transformation parameters by kriging. From estimated position of joint centers, the internal skeleton is scaled and positioned from a known reference posture to the scan position. The third step is then to attach the surface points to their corresponding skeletal segments. The last step is to check and correct the attached skin points around some joints so as to respect the skin to segment structure specific to a DHM. Compared to the method used in the past by manually adjusting a DHM on calibrated photos of several points of views; the proposed method is operator independent and much less time consuming.

Introduction

Digital Human Models (DHM) are now used to simulate the human interacting with an environment for ergonomics assessment of a product such as a vehicle [START_REF] Duffy | Handbook of Digital Human Modeling: Research for Applied Ergonomics and Human Factors Engineering[END_REF]. For example, RAMSIS [START_REF] Seidl | RAMSIS: A New CAD-Tool for Ergonomic Analysis of Vehicles Developed for the German Automotive Industry[END_REF] is a DHM software package largely used for vehicle interior design. We can also mention two other DHMs, SAFEWORK (now Human Builder under CATIA) [START_REF] Fortin | SAFEWORK: a microcomputer-aided workstation design and analysis. New advances and future developments[END_REF] and JACK [START_REF] Badler | The Jack interactive human model[END_REF]. These digital human modeling systems are able to represent future users of different anthropometric characteristics thanks to integrated anthropometric database of different populations. The DHM skin surface is commonly composed of about 2500 triangular faces. The skin points are defined locally relatively to the corresponding skeleton's segment, distributed over a certain number of sections along the segment. They also have a kinematic linkage representing the human skeleton with more than 100 degrees of freedom allowing simulating different postures and movements. In general, joints are simplified with one to three rotation axes without translation. Number of joints and rotation axes (or degrees of freedom) is model specific. For instance, RAMSIS has a simplified spine with only seven joints approximatively approaching the locations of S2/S3 (GLK), L4/L5 (GLL), T12/L1 (GBL), T8/T9 (GBB), T4/T5 (GHB), C4/C5 (GHH), C1/head (GKH), while Human Builder includes a full representation of the thoracic and lumbar spine. In order to simulate complex task oriented motions, different data-based approaches have been proposed in recent years ( [START_REF] Wang | A Motion Simulation Tool for Automotive Interior Design[END_REF], [START_REF] Park | Toward Memory-based Human Motion Simulation: Development and Validation of a Motion Modification Algorithm[END_REF]). For this, real motions performed by volunteers need to be captured, reconstructed and analyzed in order to constitute a task specific motion database. The first step for motion reconstruction is to define a personalized DHM matching to the anthropometric dimensions of the volunteer (See [START_REF] Wang | Motion Capture and Human Motion Reconstruction[END_REF] for a review of existing motion reconstruction methods). Most of existing DHMs have a module for generating manikins from a few external anthropometric dimensions. For Human Builder, up to 70 anthropometric measurements such as stature, sitting height, weight can be used for the DHM personalization. In practice, not all required anthropometric measures are available. The missing parameters are automatically estimated from a chosen anthropometric database.

A similar procedure can be found in RAMSIS, which uses 22 anthropometric dimensions. Though generated manikins could have a good representation of external anthropometric dimensions such as stature, sitting height, etc., more detailed external body shape is not guaranteed compared to the real external shape of the volunteer. An alternative way is to adjust a DHM so as to match it by superimposition with the photos of at least two different views ( [START_REF] Wang | Validation of a model-based motion reconstruction method developed in the Realman project[END_REF]). However, this procedure is tedious, time consuming and operator dependent.

3D body scanning technologies are becoming accessible, and more and more 3D scan data are available. But creating a personalized DHM from a scanned external body shape is not straight forward, because the internal skeleton has to be identified from external body shape. More critically, the local coordinate system associated to each body segment of a DHM need to be defined so as to attach skin points to corresponding body segments. Moreover, attached skin points have to follow the skin-to-segment structure especially around the joint areas. Joint centers can be estimated using external anatomical landmarks (ALs) as proposed by Reed et al [START_REF] Reed | Methods for measuring and representing automobile occupant posture[END_REF]. However it is difficult in locating ALs and requires specific training for palpation. Moreover, joint centers alone are not enough to completely define the local coordinate system associated to each body segment.

Therefore, the main objective of the present paper is to present an operator independent method to personalize a DHM from 3D body scans data without using ALs. To illustrate the proposed approach, RAMSIS is used in the present work as DHM. The first step is to create a DHM reference template (a) and to clean the 3D scan data. The DHM template is composed of a body surface and an articulated internal skeleton. In the present work, a template corresponding to a 50 th percentile male was used. Raw 3D body scan data generally needs to be cleaned to get a usable surface mesh (b). This step is detailed in Section 2.2.

Proposed approach 2.1 Overview

The second main step is to determine the transformation which deforms the DHM template so as to match the bodyscan mesh without using landmarks. This step is detailed in Section 2.3. The outputs are the DHM surface points fitted to the bodyscan mesh (c) as well as the DHM joint centers (d) in a global coordinate system. The third step (Section 2.4) is to attach the surface points to the underlying skeletal segments. For this, body segment local coordinate systems (f) need to be calculated from the estimated joint centers. This step is detailed in Section 2.3.

The last step (Section 2.5) is to correct the attached surface points for each segment according to the segment-surface points structure specific to each DMH.

At the end, a personalized DHM (g) is generated respecting the skin to skeleton structure specified by each DHM.

Preparation

Preprocessing of scan data

In this study, the scan data in a standing position were obtained (Fig 2 ) using the body scanning system SYMCAD TM II commercialized by Telmat [START_REF] Rennesson | A Full Range of 3D Body Scanning Solutions[END_REF], which provides a textured mesh with a density of about a point every 3 mm. The bodyscan mesh data presents holes at all extremities (hands, feet, and head) as well as under the arms. The meshes can be also locally deformed due to Vicon markers attached on body surface used for motion capture.

Fig. 2. Screen capture of some subjects bodyscan mesh with the Symcad II standard posture -same scale

The scan data were first cleaned using MeshLab [START_REF] Meshlab | [END_REF] by removing duplicate vertices and faces, non-manifoldness… The number of points was also reduced to a density similar to the third body template used for the DHM to scan transformation step. The surface points are distributed in a number of sections along the underlying segment. Each section has 6 to 16 points in a ring shape. When the DHM is set to the zero reference posture (i.e. all joint angles set to zero), most of the sections are parallel. In order to control the shape deformation around a joint area, a few mobile sections are introduced. When the proximal and distal segments are aligned, these mobile sections are superimposed. For instance, a mobile section is defined around the elbow. Figure 4 shows how the mobile sections move when flexing the elbow.

Fig. 4. Elbow shape in different flexion angles by a step of 45° (view from top)

In the present work, a RAMSIS template was created based on an existing model corresponding to a 50 th percentile male. The template was positioned by manually adjusting joint angles approaching the scan standing position.

Fig. 5. RAMSIS template, from the left : zero reference posture, kinematical skeleton, scan posture

Third human body template

In order to calculate the DHM to scan transformation (see 2.3), a third body surface reference template is needed. In the present work, the template mesh corresponds to an average male from MakeHuman [START_REF] Makehuman | makehuman.org[END_REF]. The template mesh was initially posed with Blender [START_REF] Blender | [END_REF] in a posture similar to the scan posture. A particular attention was paid to the forearm axial rotation (pronation/supination) for a good fit to scan. The template mesh was simplified manually using Blender by removing internal parts (tongue, teeth, eyes …) as well as small faces. Moreover, due to the hole at the top of the head in the scan data, the top of the head of the MH template was also removed to avoid large head deformation. At the end, the third template used is composed of 17992 faces and 9020 vertex (Fig 6).

Transformation from DHM template to scan data

In order to avoid the use of body landmarks to fit the DHM template to the scan, we propose to use a third body surface template created from MakeHuman (MH). The MH template is fit to the scan and also to the DHM reference surface using mHBM, a robust markerless non-rigid transformation tool developed by Japan National Institute of Advanced Industrial Science and Technology (AIST) [START_REF] Yamakazi | Markerless landmark localization on body shape scans by non-rigid model fitting[END_REF]. As the two deformed surfaces from a same MH template have the same data structure, a DHM to scan transformation can then be calculated with a kriging operation [START_REF] Dumas | Soft tissue artifact compensation by linear 3D interpolation and approximation methods[END_REF] using all or part of surface points. In the present study, half of the points were used.

Fig. 7. Principle of the determination of the DHM to scan transformation

Attaching the surface points to skeleton

The third step is to attach the surface points to the underlying skeletal segment. For this, the DHM internal skeleton is dimensioned at first and then positioned in a reference posture with the local coordinate system of all skeleton segments being defined. The RAMSIS reference posture with all joint angles equal to zero was used in the present work, corresponding to a standing posture with the two extended arms pointing laterally (e) (figure 5). Then the scaled skeleton is placed in the same position as the scan posture using an inverse kinematic algorithm with the joint centers as targets, which are obtained by applying the DHM to bodyscan transformation. Once the joint angles corresponding to the scan posture are obtained, the local coordinate system of each skeleton segment can be calculated. The surface points are then projected into the corresponding body segment local coordinate systems (f).

Surface points correction around joint centers for respecting DHM surface to skeleton structure

Once surface points are attached to underlying segments, we need to check if the attached points in the reference posture respect the RAMSIS model structure (see 2.2.2). In particular, the points in the mobile sections and next sections around a joint have to be superimposed. In the present work, the points of one mobile section were chosen to impose to other mobile sections. After the correction step, a personalized RAMSIS DHM is created.

Implementation and examples

The proposed approach has been implemented using Matlab and tested under the environment of RPx, a Matlab based digital human modelling software package developed jointly by IFSTTAR and Renault [START_REF] Wang | A Motion Simulation Tool for Automotive Interior Design[END_REF] for human motion reconstruction, analysis and simulation. Figure 8 shows an example of a personalized RAMSIS DHM emulated in RPx. Figure 9 shows the personalized DHM in a car ingress motion.

Discussion and future work

In this paper, a method for generating a personalized digital human model (DHM) from a 3D body scan is proposed. The use of a third body template makes it possible to create the correspondence between DHM and scan surface points and then to calculate the transformation from DHM to scan reference surface points without using body landmarks. By applying the DHM to scan transformation, DHM joint centers are estimated in the scan position. In a recent study [START_REF] Nérot | An assessment of the realism of digital human manikins used for simulation in ergonomics[END_REF], the location of the joint centers of RAMSIS and Human Builder DHMs generated using a few anthropometric dimensions was compared with those obtained from two perpendicular x-ray images. Results showed that both RAMSIS and HB could correctly reproduce external anthropometric dimensions, while the estimation of internal joint centers location presented an average error of 27.6mm for HB and 38.3mm for RAMSIS. As an accurate location of joint centers is required for biomechanical and ergonomics analysis of human activities; the joint centers estimated only by fitting the DHM to external body surface may not be anatomically correct. If the DHM template is too far from the anthropometric characteristics of a target person, a big error in joint centers location may be expected. Work at IFSTTAR in collaboration with LBM of ParisTech Art et Métier is in progress to establish the relationship between the characteristics from external body shape and those of the internal skeleton [START_REF] Nerot | Estimation of hip joint center from the external body shape: a preliminary study[END_REF]. One way is to predict joint centers from external body shape using a statistical shape modeling approach like that used by Yamazaki et al [START_REF] Yamakazi | Markerless landmark localization on body shape scans by non-rigid model fitting[END_REF] or Reed et al [START_REF] Reed | Developing and implementing parametric human body shape models in ergonomics software[END_REF] to predict anatomical landmarks if data containing both external body shape and joint centers are available. Compared to the method by manually adjusting a DHM to the calibrated pictures of a subject, the present method is operator independent and much less time consuming. Using the proposed method, it takes less than 1 minute to generate a personalized DHM instead of 20 minutes by manual adjusting method.

Fig. 1

 1 Fig.1presents an overview of the proposed approach for generating a personalized DHM model from a 3D bodyscan. It consists in four main steps:1) Preparation 2) Transforming DHM template to bodyscan 3) Attaching surface points to the DHM skeleton 4) Correcting the attached surface points to the skeleton respecting the DHM structure.

Fig. 1 .

 1 Fig. 1. Overview of the proposed approach

Fig. 3 .

 3 Fig. 3. Data processing results for the bodyscan post treatment

Fig. 6 .

 6 Fig. 6. Template definition: (a) MakeHuman posed and simplified, (b) with the top of the head removed, (c) in red the rigidified vertices

Fig. 8 .

 8 Fig. 8. From left to right are bodyscan, fitted MH template (grey) superimposed with fitted DHM template (red points), bodyscan superimposed with the personalized DHM skin points (red points), and the personalised DHM.

Fig. 9 .

 9 Fig. 9. Personalized DHM performing an ingress motion


Tecla Bonci

Valentina Camomilla

Raphaël Dumas

Laurence Cheze

Aurelio Cappozzo

When estimating bone

WHEN ESTIMATING BONE POSE USING PROCRUSTES SUPERIMPOSITION, ONLY THE RIGID COMPONENT OF THE SOFT TISSUE ARTIFACT IMPACTS ON END RESULTS

T. Bonci1,2,3,4,5 , V. Camomilla 4,5 , R. Dumas 1,2,3,5 , L. Chèze 1,2,3,5 , A. Cappozzo 4,5 Institutions:

Introduction and Objectives: To reconstruct the pose of a bone in the 3D space, stereophotogrammetry and a skin-marker cluster are used. The movement between skin-markers and the underlying bone is regarded as an artefact (soft tissue artefact, STA) having devastating effects on end results. This STA can be described at marker-cluster level by a series of geometrical transformations, such as rotations and translations (cluster rigid motion: CRM), homotheties and stretches (cluster non-rigid motion: CNRM). Recent studies quantified these STA components [1]- [4], showing that CRM is normally predominant with respect to CNRM. Based on this observation, it is concluded, either explicitly or implicitly, that CNRM has a limited impact on bone pose estimation (BPE) and that STA compensation should concentrate on CRM. This study disputes the message carried by this statement and demonstrates that CNRM does not have a limited effect on BPE accuracy, but, rather, it has no effect whatsoever and that this is the case independently from its magnitude relative to CRM. For this reason, the only STA component to be compensated for is CRM. Methods: The data obtained in [5] and relative to the trajectories of both skin and pin thigh markers recorded during 5 trials of each of 3 running subjects (S1, S2, S3) were used. For each trial and subject, a bone anatomical frame (AF) was defined, based on the pin markers, and the movement of four skin-markers reconstructed in the AF. Relevant displacement vectors were represented, in each k-th sample, as an STA vector field , (k=1:n) [6]. This field was decomposed into modes, by projecting it onto an orthogonal base of unit vectors chosen so that the first six modes represent rotations and translations (CRM), and the further six homotheties and stretches (CNRM) [6]:

The energies associated with CRM and CNRM were calculated as

E 1 and E 1 ,
respectively.

The ratio between the former and the latter energy (R) was determined.

To assess the impact that CNRM has on BPE, a Monte Carlo Simulation was used to generate a set of one thousand STA fields that have the same CRM and amplified CNRMs:

The amplification factor r was randomly generated in the range from 1 to 2√# so that the mean CNRM energy of this set was equal the CRM counterpart. Then, the CRM components were removed from both the measured and the amplified STA fields, generating STA fields affected only by the real, Discussion: Skin-marker trajectories were generated from the reference AF pose and the STA fields available ( , , , ). These trajectories were used to estimate the artefact-affected pose of the AF in the global reference frame with a Procrustes Superimposition (PS) approach. The root mean square difference between the artefact-affected and reference AF pose was calculated and considered as an error (rmse θ , and rmse px , rmse py , rmse pz , for orientation (attitude angle) and position components, respectively). The median (inter-quartile range) values of the R factors were 7(1), 18(4), and 14(15), for S1, S2 and S3, respectively. Obviously, the CNRM amplification caused an increase in the energy percentage of this component with respect to the total energy: mean (±standard deviation) values went from 8±3% (V(k)) to 52±22% (V S (k)) (Fig. a). Before and after the amplification of CNRM, errors in pose estimation were exactly the same (although the total STA energy increased) (Fig. b). In all cases, after removing CRM, not altered throughout the simulation, the error was null. The results empirically showed that, using a PS approach, only CRM has an impact on the accuracy of the BPE, independently of the amplitude of CNRM. Moreover, after removing CRM, the real p was obtained. It must be acknowledged that a reference pose obtained with pin markers was used to compute the modes and to remove CRM from skin-marker trajectories. In the present context, this choice does not constitute a limitation, but does simply show that the CNRM has no effect on BPE. Results suggest that, in a STA compensation perspective, future work should be focused on modelling and removing only CRM.

  Figure -a) Mean values of the CRM and CNRM energy in percentage of the total energy, for the measured (E CRM and E CNRM ) and simulated STA (E S_CRM and E S_CNRM ). Statistics performed over all trials and subjects (S1, S2, S3). b) Boxplots (minimum, lower quartile, median, upper quartile, and maximum) of the BPE errors, for position and orientation, relative to all the STA fields available.

Université de Lyon, F-69622, Lyon, France

Université Claude Bernard Lyon 1, Villeurbanne, France

IFSTTAR, UMR_T9406, Laboratoire de Biomécanique et Mécanique des Chocs (LBMC), F-69675, Bron, France

Department of Movement, Human and Health Sciences, Università degli Studi di Roma ''Foro Italico'', Rome, Italy

Interuniversity Centre of Bioengineering of the Human Neuromusculoskeletal System, Università degli Studi di Roma "Foro Italico", Rome, Italy


Georgios Christodoulis

Manuel Selva

François Broquedis

Frédéric Desprez

Olivier Muller

Franc ¸ois Broquedis

An FPGA target for the StarPU heterogeneous runtime system

Keywords: FPGA, StarPU, Vivado HLS, Task programming system, Heterogeneous runtime system

I. INTRODUCTION

Computer architectures are getting more and more complicated, exposing massive parallelism, hierarchically-organized memories, and heterogeneous processing units. The complexity of the platforms above imposes a compromise between portability and performance. To exploit heterogeneous processing units integrated into a single machine, programmers often need to deal with hardware-specific low-level mechanisms explicitly. For example, they must expressly handle memory transfers between the main memory and the private memories of the accelerators.

Solutions have been proposed to tackle the portability versus performance dilemma. Existing parallel programming models and tools evolved towards the support of heterogeneous architectures, broadly adopting the Graphical Processing Units (GPU) as complementary computational resources to the general purpose processors (CPU).

Many of these parallel programming systems rely on the concept of tasks. A task is an autonomous work unit that can be seen as a portion of code with inputs and outputs, possibly depending on other tasks. Within this model, an application is organized in the form of a dynamically created task graph. Tasks are later assigned to processing units by a software layer called the runtime system.

Tasks prevalence is mainly due to their affordable cost and flexibility, enabling the implementation of dynamic load balancing strategies with very low overhead. Unfortunately, none of the well-established task-based parallel programming models for heterogeneous architectures support Field-Programmable Gate Arrays (FPGA).

Unlike CPUs and GPUs, FPGAs can be configured to fit the application needs. They contain arrays of programmable logic blocks that can be wired together to build a circuit specialized to the targeted application. For example, FPGAs can be configured to accelerate portions of code that are known to perform poorly on CPU or GPU [START_REF] Dorrance | A scalable sparse matrixvector multiplication kernel for energy-efficient sparse-blas on fpgas[END_REF]. The energy efficiency of FPGAs is also one of the primary assets of this kind of accelerators compared to GPUs, which encourages the scientific community to consider them as one of the building blocks of large-scale low-power heterogeneous platforms.

Unfortunately, their programming complexity has prevented High-Performance Computing (HPC) developers from using them broadly. When the entry point to hardware design was a Register Transfer Level (RTL) description, with Hardware Description Languages (HDL) like VHDL and Verilog, FPGAs were mostly confined to hardware designers. The capabilities delivered by modern High-Level Synthesis (HLS) tools bridge the gap between software programmers and FPGA. HLS tools generate an RTL description of an algorithm, provided typically in C or C++. Despite their rapid evolution, those tools, like Catapult from Calypto or Vivado HLS from Xilinx, still requires a certain level of expertise in hardware design to benefit from their abilities fully.

Our contribution is a framework, which enables developers to use FPGA as another accelerator within a traditionally heterogeneous HPC platform. This framework relies on existing HLS tools and heterogeneous task programming frameworks. We focused on ease of development and portability without sacrificing performance.

The remainder of the paper is organized as follows. Section II describes related work. Section III introduces the background for this work. Section IV presents the fundamental concepts behind our framework and the way we implemented it. Section V gives an evaluation of the HEAVEN framework, while conclusion and future works lie in Section VI.

II. RELATED WORK

Numerous proposals to ease the programmability of an FPGA enhanced platform have appeared with different features regarding the types of supported processing units, the programming abstraction (from tasks to threads), the way they manage data transfers as well as their scheduling support.

OpenCL [START_REF]Opencl specification[END_REF], [START_REF] Gaster | Heterogeneous Computing with OpenCL: Revised OpenCL 1[END_REF] is a standard framework offering both a device-side language and a host management layer to exploit multiple hardware accelerators from a single application, created initially for GPU and later on extended for FPGA. The device-side language, called OpenCL C, implements both task and data parallelism and provides ways to map data efficiently on specialized processor devices. The application programmer can attach command queues to processing elements and explicitly assign them tasks from the application. This way, OpenCL stands as a great programming environment to experiment with heterogeneous architectures and implement higher-level layers, like runtime systems but may appear too low-level for application programmers, as the assignment of tasks and the data transfers between host and device memories are explicit.

At thread level, hthreads [START_REF] Andrews | Achieving programming model abstractions for reconfigurable computing[END_REF], is a unifying programming model aiming to specify application threads running within a hybrid CPU/FPGA system. A unified multiprocessor abstracts every component of the platform. Application developers can express the concurrency of their computations using threads, arranged in a single program, which can either be compiled to run on a CPU or synthesized to run on an FPGA.

The Heterogeneous System Architecture Foundation [START_REF] Foundation | The hsa foundation specifications version 1.0[END_REF] (HSA) is a consortium of industrial and academic partners discussing ways of standardizing a whole ecosystem for heterogeneous programming. In particular, HSA provides a unified view of fundamental computing elements. HSA allows writing applications that seamlessly integrate CPUs (called latency compute units) with GPUs (called throughput compute units) while benefiting from the best attributes of each. The essence of the HSA strategy is to create a single unified programming platform providing a strong foundation for the development of languages, frameworks, and applications that exploit parallelism. The HSA platform is designed to support high-level parallel programming languages and models, including C++, OpenCL, OpenMP, Java, and Python, to name a few. HSAaware tools generate program binaries that can execute on HSA-enabled systems supporting multiple instruction sets and also can run on existing architectures without HSA support. An HSA implementation is a system consisting of a heterogeneous hardware platform that integrates both CPUs and GPUs, which operate coherently in shared memory, a software compilation stack, a user-space runtime system, and some kernel-space system components.

The OpenMP accelerator support provides explicit ways of offloading portions of applications to accelerator devices. Relying on this OpenMP support, recent work [START_REF] Sommer | Openmp device offloading to fpga accelerators[END_REF] have been proposed to handle FPGA. This work focuses on making data transfers automatic based on dependencies provided at the OpenMP level and does not discuss heterogeneous scheduling questions. Relying on StarPU, our proposal will benefit from smart scheduling heuristics based on cost models for both computation and data transfers times.

OmpSs is a programming environment developed at BSC that enables the execution of task-based applications on heterogeneous platforms. It can be seen as a fork of OpenMP with original extensions. Like OpenMP, OmpSs provides compiler directives to turn serial C, C++ and Fortran applications into parallel applications initially able to offload computations on GPU. The programming environment comes with a dedicated compiler and a runtime system responsible for distributing tasks over the CPUs and the accelerator devices, managing data transfers between host and devices memories in a transparent way. OmpSs has then been extended to support offloading on FPGA [START_REF] Filgueras | Ompss@zynq allprogrammable soc ecosystem[END_REF]. Compared to our proposal, this extension focuses on embedded heterogeneous architectures such as the Zynq from Xilinx integrating a CPU and an FPGA on the same die. Blaze [START_REF] Huang | Programming and runtime support to blaze fpga accelerator deployment at datacenter scale[END_REF] provides programming and runtime support that enables rapid and efficient deployment of FPGA accelerators at warehouse scale. It builds upon Apache Spark, a widely used framework for writing Big Data processing applications.

Among many industrial attempts to use FPGA as an accelerator in a heterogeneous environment, Amazon EC2 F1 instances can be considered to deliver high flexibility to developers [START_REF]Amazon ec2 f1[END_REF]. Instances come in two sizes, with up to eight 16 nm Xilinx UltraScale Plus FPGA devices connected to each one of them through PCIe. They provide all the necessary tools to develop an AFI (Amazon FPGA Image) that can later be used to configure one of the devices of the instance. The runtime decisions are yet to be managed by the developer though.

III. BACKGROUND

We now introduce the background required for our work. In this paper, we use the term host to refer to the general purpose processor of the machine along with its main RAM, and with all its software environment including its operating system.

A. Vivado HLS and Vivado

Vivado HLS is one of the broadly used HLS tools. From a C++ description of a kernel, Vivado HLS automatically generates an RTL description. Vivado HLS delivers performance by allowing the developer to monitor and tune the behavior of its kernel by using some pragmas to exploit the fundamentally concurrent parts of it. From an HPC developer, a relatively efficient hardware description of a C/C++ application should arrive after a basic understanding of the design principles through the rich documentation and code examples the tool provides. Once an RTL description has been generated by Vivado HLS, it needs first to be instantiated along with all the other required components for the final design. Then, the complete RTL description of the final design is synthesized to logic gates by Vivado.

B. RIFFA:

The link between the host and the FPGA RIFFA (Reusable Integration Framework for FPGA Accelerators) is a simple framework for communicating data from a host to a FPGA via a PCIe bus [START_REF] Jacobsen | Riffa 2.1: A reusable integration framework for fpga accelerators[END_REF]. It is then made of a host part and a FPGA part. RIFFA is available on Linux and for several series of boards, from different vendors such as AVNet, Xilinx and Altera.

On the host side there is a clean interface to exchange data between the host and the device, that relies on a custom Linux kernel driver. The communication is managed using direct memory accesses transfers and interrupt signaling.

On the FPGA side, RIFFA provides a clean interface for accessing up to twelve signal controlled completely independent bidirectional channels. Each channel has two sets of signals; one for receiving and one for sending data. For a transaction, the sender has to notify the receiver for the upcoming operation, providing supplementary information regarding the data size, a potential offset, and whether other transactions are expected or not. Once they receive the acknowledgment from the receiver, they synchronize themselves using a two signals handshaking protocol for the consumption of the data.

C. StarPU

StarPU [START_REF] Augonnet | Starpu: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] is a task-based runtime system targeting heterogeneous machines comprising several kinds of processing units, from traditional multicore CPUs to GPU accelerators and coprocessors like the Intel Xeon Phi. StarPU comes with an API to express task parallelism from a C/C++ or Fortran application. StarPU is data-flow based, meaning the programmer has to specify for each task the input and output data, thus creating dependencies between tasks. This way, a StarPU application generates a graph of tasks at runtime, with edges representing the dependencies between them. Unlike most task-based parallel programming libraries, each StarPUtask can embed multiple implementations, one for each kind of processing unit for example. When a task becomes ready, the StarPU tasks scheduler dynamically assigns it to what it considers to be the best processing unit at the given time. To do so, StarPU takes into account the availability of targetspecific implementations for the current task, the availability of the hardware resource at the given time and some modelbased performance estimation of the task on the different processors/accelerators. On top of that, StarPU also handles data transfers automatically and transparently, and can make smart decisions on whether some data should stay in the accelerator memory or not, depending on the input data of the next tasks to be executed.

In this work, we focus on using the C/C++ API of StarPU. In this case, once the different tasks making the application have been identified, the programmer has to define task types, called codelets in the StarPU terminology. A codelet establishes the number of inputs and outputs of a task along with all its implementation for one or more different types of processing units. Once codelets have been defined, the application developer has to create the main StarPU program. This program first registers in the StarPU runtime, the inputs and outputs data that will be read and produced by tasks. Then the program creates tasks along with their codelets and their inputs and outputs and submits them to the StarPU runtime system. The task graph represents the execution of a StarPU application.As a consequence, it is fully dynamic and can, for example, depend on user-specified inputs.

When a task is submitted, the StarPU runtime first checks for input dependencies. If these dependencies are satisfied, the scheduler decides on which processing units the task will execute, among all the possible ones specified in the task codelet. At this time, the runtime pushes the task into the ready queue of the chosen processing unit. If some input dependencies are not satisfied, the runtime pushes the task into a global waiting queue instead.

At initialization time, StarPU creates a worker for each processing unit of the underlying machine. A worker is a thread, running on the host, responsible for popping tasks the runtime has scheduled in its ready queue and for executing them. The worker is also responsible for handling data transfers required for the task execution.

Regarding the scheduling decisions, StarPU comes with many policies which can be chosen by the application. Many of these strategies rely on cost models for task execution and data transfers. These models can either be initialized from a previous run of the application or from within the runtime system itself, by profiling during the execution.

IV. THE HEAVEN FRAMEWORK

This section describes in details the HEAVEN framework. We first present the components involved in the framework with a focus on our contributions. Then, we present the HEAVEN framework from the user point of view, showing how to use it.

A. Architecture

Figure 1 shows all the components involved in the HEAVEN framework. Because we target heterogeneous machines including FPGA accelerators, the framework is split between the host part shown on the top of the figure, and the FPGA part displayed on the bottom. The red components of the figure are the contributions of this work. Because we base our framework on StarPU, a complete figure would also include the alreadysupported GPU side. For clarity, only the FPGA side which is our contribution is shown here.

On the host's side, we then had to extend StarPU to integrate the capability to use FPGA accelerators. To ease the communication between this extended version of StarPU and the FPGA, we developed Conor, a library providing abstractions on top of RIFFA.

On the FPGA's side, we developed a connector allowing to send inputs and get outputs to and from hardware tasks generated from Vivado HLS. This connector enables the extended StarPU runtime running on the host side to communicate with the hardware tasks through PCIe without any intervention of the application developer. In the current version of the framework, the configuration of the FPGA must be done once, before the execution of the StarPU application. The framework supports up to five FPGA devices, while the number of tasks StarPU can execute on each one of them is constrained by the number of available communication channels as well as the spatial constraints. In the current version of the framework, if several hardware tasks are present in the system, either on several FPGA devices and a single FPGA device, they must all have the same type. Said differently, they must expect the same inputs and provide the same outputs.

We now describe the internal of the components at the heart of this work. The components provided by the application developer, shown with dashes in Figure 1 will be described in detail in Section IV-B.

1) Extended StarPU: To include FPGA support in StarPU, we created a new type of worker dedicated to FPGA devices. At initialization time, a worker is created for each hardware task present on the FPGA devices of the machine as shown in Figure 2. As any other StarPU's worker, each one of these FPGA workers has its own tasks ready queue which is filled by the StarPU's scheduler.

The role of the worker is then to pop tasks from this queue and to launch their execution on the FPGA side. When the worker is idle, it pops the ready task on the top of the queue and first send its input data to the FPGA through Conor. This sending operation is currently done by a function provided by the programmer as shown in Section IV-B. Nevertheless, the extended StarPU runtime has all the required information so that it can automatically perform this operation. Sending the inputs data triggers the computation on the FPGA side automatically, as soon as all data have been transmitted. Then, the function provided by the programmer must wait for the output data to come back from the FPGA. To completely decouple the StarPU runtime and the StarPU application from low level concerns about communication with the FPGA, we developed a communication library called Conor. All the FPGA communications happening inside StarPU go through Conor using its simple API for sending and receiving data.

2) Conor: The Conor FPGA communication library is built on top of RIFFA as shown in Figure 1. RIFFA provides a primitive infrastructure to send and receive non-structured data to the device, where low-level information should still be managed for a clean adaptation.

Conor first provides functions for the transmission of structured data, such as arrays and vectors, over the PCIe in an optimal way. These functions abstract away all the low-level hardware concerns, and just let the caller say "send this buffer which size is that to the FPGA". Conor also provides functions to handle the control signals regarding the state of the device. These functions allow controlling the coherency of the state of the board.

3) Connector: On the FPGA side, as shown on the bottom half of Figure 1, our framework includes a connector allowing the communication between hardware tasks created with Vivado HLS and RIFFA. Figure 3 shows the internals of this connector. On the hardware tasks side, the framework relies on the ap_hs protocol of Vivado HLS to handle the communication. The data transfers are controlled by a two signals handshake protocol. On the other hand, a RIFFA channel has two sets of signals, one for receiving and one for sending data. The transfers are also controlled by a handshake protocol on two phases, one to inform the parties about the type of the upcoming transaction, and another to perform the transaction itself. The connector ensures the consistency of the communication between these two different protocols using information provided by the hardware tasks through the control stream.

Fig. 3: The connector allows communication between RIFFA and the hardware tasks.

B. Developing an application

To illustrate how to use the HEAVEN framework, we will focus on a basic application consisting in computing the product of two matrices by blocks, as shown in Figure 4. In this example, each square block represents a sub-part of a matrix with a specific size, let's say 64 by 64. Then, the matrix A is of size 192x64, B is of size 64x192, and C is of size 192x192. To compute the content of the matrix C, we can perform in parallel the computation of the nine blocks making it. Each one of these sub-computations will require a sub-part of matrix A and a sub-part of matrix B as inputs. To describe this application using the HEAVEN framework, the application developer must write both the StarPU application along with the hardware implementation of its tasks as shown by dashed boxes in Figure 1. In the case of the blocked matrix multiplication we use in this section, tasks are completely independent. Nevertheless, because the HEAVEN framework is an extension to StarPU, it also supports applications with task dependencies.

1) Writing the StarPU application: First, the programmer must describe the codelet for the tasks computing one block in the result matrix. Figure 5 shows how to do that using StarPU's C API. This codelet provides both a CPU implementation and an FPGA one. It also specifies that the task has two inputs and one output through the nbuffers and modes attributes. Then, the fpga_mult function referenced in the codelet must be defined has shown in Figure 6. This function is mainly divided in four steps:

• Ask to Conor to reserve an hardware task;

• Ask to StarPU where are inputs and outputs and on which part of them the task must operate; • Ask to Conor for the result and then copy back this result into the C matrix. In this function, the order used to send the inputs is crucial. It must be aligned with the expected inputs order in the hardware task implementations as explained in Section IV-B2.

Again, this function is not yet automatically generated by the extended StarPU runtime but all the information required to do so is available to it. Said differently, in the next version of the framework, this function will no more be written by the programmer. The runtime will automatically send and receive inputs and outputs using the information associated with the task.

Then the programmer must write the main part of the application responsible for allocating matrices, register them with StarPU, launch tasks and wait for their completion as shown in Figure 7.

The procedure for launching the task is completely uniform; it does not know anything about task types as shown in Figure 8. It consists of creating the tasks, associating them with the codelet defined in Figure 5, setting their dependencies, and finally submitting them.

2) Writing hardware tasks: On the FPGA side, the application developer must provide the implementation of hardware tasks. The HEAVEN framework relies on Vivado HLS for this task. Using this tool, the application programmer describes the implementation of its hardware tasks in C++. This description consists in:

• Receive the block from A and B;

• Compute the block of C as the product of the blocks from A and B; • Send this C computed block. For receiving inputs and sending outputs, the application developer is provided two high-level objects which type is Vivado HLS hls::stream. The first one is used to receive inputs and the second one to send outputs. The function also has a third hls::stream parameter used for control. The interface between these three high-level objects with the PCIe bus is handled transparently by the framework through the connector and RIFFA as described in Section IV-A.

For the receiving part of the function, the programmer must read the blocks for the A and B input matrices in a given order and save them in local variables which are automatically translated to BRAM blocks by Vivado HLS. This reading order must be respected in StarPU when inputs are sent to the FPGA as presented in Section IV-B1.

For the computation step, the application developer has first to declare the C matrix as a local variable that will also be allocated into BRAM blocks by Vivado HLS. Then the for (uint32_t x = 0; x < 9; x++) { for (uint32_t y = 0; y < 9; y++) { spu_task * task = spu_task_create(); task->cl = &cl; / * Get handlers for each block * / task->handles[0] = spu_get_sub_data( A_handle, 1, y); task->handles [START_REF] Dorrance | A scalable sparse matrixvector multiplication kernel for energy-efficient sparse-blas on fpgas[END_REF] = spu_get_sub_data( B_handle, 1, x); task->handles [START_REF]Opencl specification[END_REF] = spu_get_sub_data( C_handle, 2, x, y); starpu_task_submit(task); } } Fig. 8: Submit tasks to the StarPU runtime.

effective computation of C is done using three nested for loops just as in a software implementation.

Finally, the C matrix must be sent back to host side using the second hls::stream object provided to the function. In case of tasks with several outputs, as for the inputs, the sending order must be respected on the StarPU side when receiving the data.

3) Generating the bitstream: Once the RTL description of the hardware tasks has been generated with Vivado HLS, the application developer must generate the final bitstream to be used to configure the FPGA. For that task, the HEAVEN framework provides a parametric wrapper written in Verilog. This wrapper is only configured by the application developer to specify the number of hardware tasks to be instantiated. The wrapper combines the specified number of hardware tasks along with the hardware part of RIFFA, and the connector between RIFFA and the hardware tasks. The application developer then uses Vivado to synthesize the complete design to the final bitstream.

V. PRELIMINARY EVALUATION

A. Environment of experiments

We executed the experiments in a hybrid machine, where the FPGA was connected to the host via PCIe. The host CPU was a 64-bit Intel Xeon W3530, a two-way hyperthreaded quad-core, for a total of 8 virtual cores, running at the speed of 2.80GHz, with a smart cache of 8MB, 256KB of L2 and 8192KB of L3, a thermal design reference power of 130W, and a semiconductor size of 45nm. It was coupled with 12GB of DRAM at 1333MHz.

The FPGA we used for our experiments was Xilinx Virtex-7 VC709, a board using the XC7VX690T chip, with 693,120 logic cells, 3,600 DSP slices, 52,920kb of BRAM memory, up to 4GB of RAM at 1866Mbps, and an 8-lane PCIe edge connector.

B. Application set-up

Within the context of this study, we opted to evaluate the performance we could obtain using our platform on the blocking version of matrix multiplication introduced in Section IV-B. For the given analysis, we did not consider hybrid scenarios, where the scheduler could use within the same execution CPU and FPGA implementations for the tasks. Thus, for every experiment under test, we went through two executions, one using only FPGA tasks and one using only CPU tasks. We expected that the communication overhead would significantly downgrade the performance of the FPGA version. As a consequence, to study this effect, we focused on two constant task sizes for the entire range of experiments. A light-weight task operates on blocks of 64x64 integers while a heavy-weight task operates on 256x256 blocks.

We also vary the size of the input matrices A and B. We evaluated nine different sizes as shown in Figure 9. For each one of these nine scenarios, the number of tasks that can be deduced from the way the input matrices were sliced. Each scenario is further divided in two by using either light-weight of heavy-weight tasks. For example, the scenario in Figure 9e corresponds to a setup for the 64x64 task size and to another setup for the 256x256 one. In this scenario, for both 64x64 and 256x256 setups the A and B matrices are split into five slices leading to a total of 25 tasks. In the 64x64 setup the size of A is then 320x64 and the size of B is 64x320. In the 256x256 setup the size of A is then 1280x256 and the size of B is 256x1280.

In all our experiments, the FPGA was configured with four instances of the hardware task computing a single block of the result matrix. Hence for both architectures, CPU and FPGA we had an equal amount of independent computational entities because in the CPU case, StarPU allocates by default a number of workers equals to the number of physical cores in the host (4 in our case).

C. Results

We first evaluated the per task execution time for each architecture. This corresponds to the scenario presented in Figure 9a. In that case, the computation of the block corresponded to the entire application. The per-task performance we obtained is presented in I. The bigger task shows a significant performance difference when executed on an FPGA compared to when executed on a CPU. In the smaller task, the ratio between the computation and the communication is small, resulting in comparable behavior for both architectures. With the increase of the task size, the cost of data transfer (linear) has a minor impact compared to the cost of the computations (cubic).

Then we evaluated how the extended StarPU runtime behaves when the number of tasks increases. In this case, the expected total execution time should theoretically, follow the formula:

T ime =    tpt + α ,#tasks < #conT asks #tasks #conT asks × tpt + α ,#tasks ≥ #conT asks (1
) In this formula, T ime refers to the total execution time of the application, tpt to the execution time of a single task which value is shown in Table I, #tasks to the number of tasks, #conT asks to the number of concurrent tasks and α to a constant representing the initialization overhead. Since within our experiments we could execute up to four tasks concurrently in every scenario, the expected total execution time should follow the:

T ime =    tp + α ,#tasks < 4 #tasks 4 × tpt + α ,#tasks ≥ 4 (2)
Figure 10 shows how the total execution time evolves with the number of tasks both for the CPU only version and the FPGA only version. For both light-weight and heavy-weight tasks, Figures 10a and Figure 10b, we observe a linear evolution of the total execution time validating the theoretically expected one. We can then conclude that for the sizes studied in our experiments, there is no additional overhead coming from task management.

VI. CONCLUSION AND PERSPECTIVES

We have presented our ongoing effort to create a task programming systems targeting heterogeneous architectures including FPGA. Leveraging the StarPU programming system and its runtime along with HLS tools, the HEAVEN framework provides an easy and portable way for writing and running applications on top of heterogeneous architectures including FPGA.

The perspective of this work is numerous. First, as mentioned in Section IV-B1 we must make the communication with the FPGA completely transparent for the application developer. This is possible because the StarPU runtime already knows what the inputs and the outputs of each task are. The only thing the programmer will have to specify is the order in which inputs must be sent and outputs must be received. We also plan to extend the framework so that different types of hardware tasks could be used at the same time. This new feature will impact mainly the Conor library that must be aware about the type and the location of each hardware task in the FPGA devices. In this paper we have presented results for executions using only the cores of the host or only the hardware tasks in the FPGA. Because our goal is to leverage heterogeneous architectures, now that we have an operational framework, an immediate perspective is to execute applications combining CPU tasks, GPU tasks and FPGA tasks. We already started to study these executions and we believe that it may lead to interesting research directions regarding the scheduling heuristics of the StarPU runtime. As presented in Section IV-B2, the framework is currently using the FPGA BRAM memory to store inputs and outputs for the hardware tasks. Because the amount of available BRAM is usually in the order of megabytes, we are planning to investigate how the framework should be extended so that we can benefit from the biggest main memory coming along with the FPGA in boards such as the VC709 we used in our experiments. In this case, the application developer will be provided pointers into this main memory as the entry point for the implementation of hardware tasks. The StarPU runtime will then be responsible for managing the memory allocations in the board's main memory with the help of an additional hardware component in the FPGA. Finally, we want to study how dynamic reconfiguration could be integrated into the framework. For sure, this would be viable only if reconfiguration time is low compared to the task granularity.

The StarPU scheduler will then need to take into account the reconfiguration time.

Fig. 1 :

 1 Fig. 1: Overview of the different components involved in the HEAVEN framework. The upper part of the diagram shows the components on the host's side, while the lower part shows the components on the FPGA's side. The red components are the contributions of this work and the dashed components are the ones to be provided by the application developer.

Fig. 2 :

 2 Fig. 2: One FPGA worker is created for each hardware task present in the machine. Workers communicate with the FPGA through Conor only.

Fig. 4 :

 4 Fig. 4: Blocked matrix multiplication example. A is of size 192x64, B is of size 64x192 and C is of size 192x192.

Fig. 5 :

 5 Fig.5: Task codelet for a task computing one block of the result matrix. The task has two inputs and one output. The codelet provides both a CPU implementation and an FPGA one.

Fig. 6 :

 6 Fig. 6: Definition of the fpga mult function.

Fig. 7 :

 7 Fig. 7: Main part of the StarPU application.

(a) 1 Fig. 9 :

 19 Fig. 9: The different inputs matrices size evaluated. Both the A and B matrices are split into N slices, with N from 1 to 9.

  (a) Light-weight tasks. (b) Heavy-weight tasks.

Fig. 10 :

 10 Fig. 10: Performance results for block sizes of 64x64 and 256x256 and both architectures (CPU, FPGA).

TABLE I :

 I Performance of a single task for each architecture and task size.

ACKNOWLEDGEMENT This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program Investissement davenir