Romain Saussard
email: romain.saussard@renault.com

Boubker Bouzid
email: boubker.bouzid@renault.com

Marius Vasiliu
email: marius.vasiliu@u-psud.fr

Roger Reynaud
email: roger.reynaud@u-psud.fr

Towards an Automatic Prediction of Image Processing Algorithms Performances on Embedded Heterogeneous Architectures

Keywords: Heterogeneous Architectures, Performance Prediction, Image Processing

Image processing algorithms are widely used in the automotive field for ADAS (Advanced Driver Assistance System) purposes. To embed these algorithms, semiconductor companies offer heterogeneous architectures which are composed of different processing units, often with massively parallel computing unit. However, embedding complex algorithms on these SoCs (System on Chip) remains a difficult task due to heterogeneity, it is not easy to decide how to allocate parts of a given algorithm on processing units of a given SoC. In order to help automotive industry in embedding algorithms on heterogeneous architectures, we propose a novel approach to predict performances of image processing algorithms on different computing units of a given heterogeneous SoC. Our methodology is able to predict a more or less wide interval of execution time with a degree of confidence using only high level description of algorithms to embed, and a few characteristics of computing units.

I. INTRODUCTION

Vehicles increasingly provide ADAS (Advanced Driver Assistance System) capabilities, as passive systems (e.g. a lane departure warning system alerts the driver when the vehicle crosses a lane) or active systems (e.g. a lane centering assist system controls vehicle trajectory). Some of these systems use cameras and image processing to sense the environment and detect potential obstacles, we can cite lane detection [START_REF] Bertozzi | Gold: A parallel real-time stereo vision system for generic obstacle and lane detection[END_REF], obstacle detection [START_REF] Dagan | Forward collision warning with a single camera[END_REF], pedestrian detection [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF], [START_REF] Geronimo | Survey of pedestrian detection for advanced driver assistance systems[END_REF], etc. Image processing algorithms need high computational capabilities, because they process high amount of data. Most of state-of-the-art algorithms for ADAS run on powerful computers or not in real time. However, image processing algorithms can often be parallelized, so they can benefit from hardware accelerator like GPU or multi-core CPU. Automotive industry needs low-power high performance embedded systems to embed image processing applications on vehicles.

Semiconductor companies are moving into the market of embedded systems for ADAS with heterogeneous architectures. These architectures embed several processing units with different capabilities on the same SoC (System on Chip), often with massively parallel computing unit. We can cite the Tegra K1 SoC of Nvidia (which embed ARM, GPU and ISP), the TDA2x SoC [START_REF] Sankaran | TDA2X, a SoC optimized for advanced driver assistance systems[END_REF] of Texas Instrument (ARM, DSP and EVE vectorial processor), or the EyeQ of Mobileye [START_REF] Stein | A computer vision system on a chip: a case study from the automotive domain[END_REF].

The problem of this type of architecture is the complexity to embed algorithms. Indeed, as there are several processing units, it is not easy to find the best mapping between algorithms and processing units. Moreover, adjusting an algorithm to embed it on a single processing unit is quite time-consuming, that is why adjusting an algorithm for all processing units to determine the best association is unachievable. In light of that fact, there is a need for car manufacturers and suppliers to predict performances of algorithms on different processing units to help them to choose the best algorithm-processing unit association.

In this work, we introduce a novel methodology to meet automotive industry needs, to predict performances of image processing algorithms on heterogeneous architectures. Firstly, we describe in section II the Nvidia K1 and TDA2x heterogeneous architectures, their characteristics and capabilities. Then, in section III, we introduce our parallelism classification and the problem of kernel mapping optimization. In section IV, we present state-of-the-art of performance prediction, followed by our novel approach description in section V. We illustrate our methodology with an example extracted from an ADAS application in section VI. Finally, we discuss future works and conclude in section VII.

II. EMBEDDED HETEROGENEOUS ARCHITECTURE

Semiconductor companies such as Nvidia, Texas Instrument and Freescale propose heterogeneous architectures to meet automotive industry needs for embedding image processing algorithms for ADAS. These architectures handle high performance computing with massively parallel computing unit (e.g. GPU or vectorial processor) and lowpower consumption. In this paper, we propose a general approach which can fit with any of these SoCs, and we show some preliminary results obtained on Nvidia Tegra K1 heterogeneous SoC.

A. Nvidia Tegra K1 Architecture

The Nvidia Tegra K1 SoC is composed of quad-core ARM Cortex A15 CPU (1.5 GHz clock rate) providing ARMv7 instruction set and out-of-order speculative issue 3-way superscalar execution pipeline, see Fig. 1, each core has NEON and FPU unit [START_REF] Lanier | Exploring the design of the cortex-A15 processor[END_REF].

The Kepler GPU of the K1 is composed of one streaming multiprocessor (SMX) of 192 cores [START_REF]Nvidia kepler GK110 architecture whitepaper[END_REF], accessible with CUDA [START_REF]Cuda C programming guide[END_REF] and the new standard OpenVX [START_REF] Rainey | Addressing system-level optimization with OpenVX graphs[END_REF]. This framework is a hardware abstraction layer for image processing applications supporting modern hardware architectures such as embedded heterogeneous SoCs. OpenVX is based on the implementation of image processing kernels designed by SoCs manufacturers, beneficing from hardware acceleration of architectures.

The parallelization with CUDA is qualified as SIMT (Single Instruction Multiple Threads). GPU has three memories reachable by all threads:

• Global memory: read and write access.

• Constant memory: low latency, read access only.

• Texture memory: read access only, can interpolate adjacent data value, always handle boundary issues. The global memory of the K1 is the same for GPU and ARM, both can potentially access to the same data, see Fig. 2. Memory area of each processing unit is handled by the OS (L4T or Vibrante, which are both based on Linux kernel).

B. K1 Specific Features

The K1 also has several hardware processing units like Image Signal Processor (ISP) and Image Processing Accelerator providing fast and specific image processing algorithms such as debayering, noise reduction, lens correction etc. These units are very fast but user can only control a limited set of parameters and the chaining order of kernels is restricted.

In addition of basic instructions, ARM and Kepler GPU have specific features which can be used to accelerate image processing algorithms. Thus, ARM processor provides SIMD instructions with NEON units [START_REF] Mitra | Use of SIMD vector operations to accelerate application code performance on low-powered ARM and Intel platforms[END_REF]. In order to handle the issues of concurrent reading and writing, CUDA provides atomic instructions based on hardware design.

Nvidia GPU texture memory provides couple of advantages for reading images data. When accessing to noninteger pixel coordinates (e.g. Im[i-0.5][j-0.5]), texture units handle linear interpolation. This is cost-free because it is computed by hardware design. Moreover, texture units handle access to pixels outside the image (e.g. Im[-1][-1]), also cost-free.

C. Texas Instrument TDA2x SoC

The TDA2x is composed of four different types of programmable units. The architecture is given in Fig. 3. Firstly, it provides a 750 MHz dual core ARM A15 and a dual-Cortex-M4. These computing units do not bring that much computing capability (A15 core on TDA2x is 4 times less powerful than the one in K1), but they can be used for data management, video acquisition control / rendering, high level decision making, etc.

To handle heavy image processing tasks, TDA2x provides a mix of Texas Instrument fixed and floating point TMS320C66x DSP (Digital Signal Processor) and up to four EVE (Embedded Vision Engine) cores. TMS320C66x is the most recent DSP from Texas Instrument, it can handle up to 32 multiply-accumulate operations per cycle. Each EVE is a 650 MHz core, optimized for image processing, composed of one specific RISC processor and one 512-bit vector coprocessor.

III. HARDWARE ACCELERATION AND KERNEL MAPPING

In order to increase performance of a given algorithm, we have to identify parts of code which can benefit from hardware accelerations of a given architecture. This implies to separate the algorithm in more or less fine blocks called kernels.

A. Parallelization Level and Classification

Execution time of a given kernel can be accelerated by using different levels of parallelization of a given architecture. First, one can address parallelism at register level, which we usually call SIMD instructions or vectorization (e.g. NEON for ARM). Some compilers can handle automatic vectorization, e.g. GCC [START_REF] Naishlos | Autovectorization in GCC[END_REF], but optimization result can vary depending on compiler.

Secondly, computing units often provide multi-core (ARM provides 4 cores and GPU 192 CUDA cores on K1 SoC). This parallelism can be accessed by using different API, e.g. for multi-core CPU, with pthread (for linux), or at more generic level with OpenMP [START_REF] Dagum | OpenMP: an industry standard API for shared-memory programming[END_REF] and C++11/14 parallel features. Pthread is a low level API, it enables fine-grained control over thread management but may be difficult to implement, sometimes too much complicated (depending on code complexity). On the other hand, OpenMP is a high level API, easier to implement, but may have different performances than pthread.

At higher level, heterogeneous architectures offer multicomputing units with different capabilities. Each of these computing units may execute different kernels concurrently.

Speed increasing brought by parallelism is not the same for all kernels, of course it depends on how kernel can be parallelized. We propose four simple classes to characterize kernels parallelism degree:

• Simple parallelism P 0 : kernels for which there is no dependency between data, each operation can be executed without the result of another one, so concurrently, e.g. convolution. construction. • Parallel reduction P 2 : kernels for which parallel reduction can be used, e.g. computing sum of all values of a vector. • Iterative kernels P 3 : kernels for which each iteration depends on the result of the previous one, e.g. if

u i = f (u i-1
). Due to this dependency, kernels in their original form cannot be parallelized. However, operations executed at each iteration, e.g. f (u), may be parallelized.

Complex algorithms provide a mix of kernels which belong to different parallelism classes.

B. Kernel Concurrency

Multi-core and multi-computing units parallelisms enable to execute multiple kernels concurrently, however execution time of each concurrent kernel may be higher than execution time of one kernel without concurrency. In fact, depending on the architecture, some resources (e.g. global memory or cache) may be shared between cores or computing units, limiting performance of kernels execution.

In the K1 SoC case, we implemented a simple benchmark in order to identify this kind of limits for the ARM processor. This benchmark is a simple kernel which copies an input image to an output image, it processes by using 1, 2, 3, or 4 threads (image is divided by the number of threads, each thread copy one part of the image). As shown in Fig. 4, copy using 1 core is slower than copy using 2 cores, whereas performance is the same when using 2, 3, or 4 cores. Thus cache bandwidth of each core is slower than global memory bandwidth, so executing more than 2 memory transfers concurrently can degrade global execution time.

C. Kernel Mapping Optimization

Embedding a given algorithm on a heterogeneous architecture is a difficult task because one can not easily find how to allocate kernels on the different processing units (kernel mapping) [START_REF] Zhou | Task mapping in heterogeneous embedded systems for fast completion time[END_REF].

Let P be the vector of processing units of a given heterogeneous architecture of size m × 1, K be the vector of kernels of a given algorithm of size n × 1 :

P =       p 1 p 2 . . . p m       ; K =       k 1 k 2 . . . k n       . (1)
The matrix M , of size m × n, constitutes the mapping of K on P , given by

P = M K, (2)
with

∀(i, j) ∈ 1, m × 1, n , M i,j ∈ {0, 1}. Thus M i,j = 1 implies that kernel k j is mapped on processor p i .
Let ϕ be the dependency matrix (dependencies between kernels K), with a n × n size, defined as

∀(i, j) ∈ 1, n × 1, n , ϕ i,j ∈ {0, 1}, ϕ i,j = 1 implies that kernel k i depends on k j results.
Let:

• τ (M ) be the execution function, returning a n×1 vector corresponding to the execution time of each kernel with mapping M ; • δ(M ) be the transfer function, returning a n × 1 vector corresponding to the transfer delay needed for each kernel with mapping M ; • η(M ) be the occupancy function, returning a m × 1

vector corresponding to the occupancy rate of each processor with mapping M ; Then, let f (P, K, ϕ, τ (M ), δ(M ), η(M )) be the cost function, returning the global execution time of the algorithm executed on the heterogeneous architecture with mapping M . The aim of kernel mapping optimization is to find M minimizing f :

arg min M [f (P, K, ϕ, τ (M ), δ(M ), η(M ))] . (3)
Parameters of function f can be measured or predicted for different M . This evaluation needs all kernels implementations on all processing units (very time consuming), but prediction only needs kernels and architectures analysis (very efficient but needs deep SW / HW knowledge). In this paper, we present a methodology to estimate the execution τ (M ) and transfer δ(M ) times with little knowledge of target architecture.

Our approach is based on two levels of accuracy. The first one uses all basic information provided by manufacturers; but this level cannot deal with cache effects, compiler optimizations, memory latency, concurrent memory access, resources starvation, etc. Thus the second level overcomes these limitations by using a generic benchmark-vectors set which automatically extracts parameters (e.g. L1, L2, L3 caches, concurrent access, or compiler auto-vectorization effects on execution and transfer times, etc.). These parameters can be extracted for any architectures and compilers.

IV. RELATED WORK

The aim of performance prediction is to estimate the execution time of a given kernel on a given computing unit. It needs a characterization of computing unit (provided by manufacturer or by benchmark results), and a characterization of kernel (high level description, source code or binary file analysis).

The difficulty for performance prediction techniques is to find the best compromise between complexity of model and precision of predictions. In fact if technique needs a full optimized source code in order to estimate execution time on one target, it does not bring that much gain compared to measure real performance on real embedded system. Whereas, if technique needs only a high level description of kernel in order to predict performance for different computing unit, one can rapidly find what target will bring the best performance for that kernel.

A survey of performance modeling techniques is given in [START_REF] Lopez-Novoa | A survey of performance modeling and simulation techniques for accelerator-based computing[END_REF], according to it three main approaches for performance modeling can be found in literature: analytical modeling, machine learning, and simulation. A performance simulator is able to reproduce a computing unit behavior. It can give a lot of information, identify bottlenecks, predict performance, etc, and generally obtains a fine estimation. Some simulators address hybrid architectures (CPU+GPU), e.g. [START_REF] Ubal | Multi2sim: a simulation framework for CPU-GPU computing[END_REF].

Using an architecture simulator implies to port kernel to simulator, which represents about the same amount of work than embedding kernel on architecture. Moreover, a fine simulation of the system implies a very long execution time.

An analytical model is a set of equations which represents characteristics of the system (kernel and computing unit). Machine learning techniques extract some characteristics by using a set of code and hardware features, then it performs by using feature selection, clustering and regressions techniques to estimate execution times of a kernel, e.g. [START_REF] Kerr | Modeling GPU-CPU workloads and systems[END_REF] which addresses CPU and GPU.

Most of analytical model handle performance prediction for only one type of architecture. In [START_REF] Hoste | Performance prediction based on inherent program similarity[END_REF], authors propose a model based on 47 architecture independent characteristics. Then performances are predicted by using programs in a benchmark suite and measuring similarity between benchmark programs and application, using the 47 architecture independent characteristics. Authors show results on different CPU architectures, but do not address neither GPU nor embedded architectures like DSP.

The famous model of Hong and Kim [START_REF] Hong | An analytical model for a GPU architecture with memory-level and thread-level parallelism awareness[END_REF] addresses performance prediction for GPU with two metrics: Memory Warp Parallelism (MWP) and Computation Warp Parallelism (CWP). The aim of this two metrics is that if M W P ≤ CW P then performance is limited by memory bandwidth and latency, but if CW P > M W P then memory latency is hidden by computing operations. Authors report a 13.3% mean error on execution time estimation.

The roofline model [START_REF] Williams | Roofline: an insightful visual performance model for multicore architectures[END_REF] uses approximatively the same approach by studying the arithmetic intensity of application and memory / computational bandwidths of architecture. However, it is not used for performance prediction, only for bottleneck highlighting and code optimization. The boat hull model [START_REF] Nugteren | The boat hull model: enabling performance prediction for parallel computing prior to code development[END_REF] adapts the roofline model to provide performance prediction. It is based on a set of primitives, kernel complexity, and memory / computational bandwidths of architecture. It shows good results, 3 and 8% of error for 2 applications. However, the model is limited to kernels which fit in classification given in [START_REF] Nugteren | A modular and parameterisable classification of algorithms[END_REF], and cannot handle arbitrary code.

V. PERFORMANCE PREDICTION USING COMPUTING PROFILE

We propose a novel model for performance prediction which:

• can address multiple architectures,

• is not limited to a set of primitives,

• can be used without specific optimized source code,

• does not need a deep knowledge of architectures. Our methodology is based on a kernel descriptor, the computing profile.

A. Basic Level of Classification

A kernel can be defined as a set of well-defined instructions, and each instruction can be classified in two classes: the computing instructions, and the memory instructions. Secondly, we can classify each computing instruction. The key is that one class of computing instructions is executed by one type of unit on the computing unit, e.g. additions will be executed by the ALU, and floating point operations by the FPU.

Classes of instruction are:

• S Int: simple operations on integer, e.g. add, sub, cmp, etc. • M Int: operations with multiplications on integers, this includes multiply accumulate operations. • Float: floating points operations, instructions computed by the FPU. • Specific: specific operations which often encapsulate several instructions, e.g. div, sqrt, etc. • Branch: branching instructions, e.g. for loop, if, etc.

• Address: addressing operations, e.g. accessing to the ith member of an array. • Memory: load and store instructions. The arithmetic intensity [START_REF] Harris | Mapping computational concepts to GPUs[END_REF] is defined as the number of operations per memory instruction (load and store). Given N M , the number of memory instructions and N C the number of computing instructions, the arithmetic intensity Ia is defined as:

Ia = N C N M . (4)
This metric can be used to identify the way to optimize the algorithm, for example with the roofline model [START_REF] Williams | Roofline: an insightful visual performance model for multicore architectures[END_REF]. Thus, performance of algorithm with small arithmetic intensity is limited by memory bandwidth and latency, whereas performance is limited by computing capability for algorithm with high arithmetic intensity. In fact, memory access latency can be hidden by multiple computing operations.

On state-of-the-art architectures, it is not easy to predict the number of global memory access and the delay of each one because of sophisticated cache systems. Indeed there are different levels of cache with different latencies and bandwidths.

B. Computing Profile and Throughput

The computing profile of a kernel is an illustration of resources needed by this kernel. The aim is not to get the complexity (number of operations), but to know which classes of instructions are used.

The profile is the ratio of each class, this is the number of instructions of one class divided by the number of instructions of all computing classes. Let C be the set of computing instructions (no memory instructions):

C = {S int, M int, F loat, Specif ic, Branch, Address} (5)
r c = N c i∈C N i . (6)
The computing profile shows us which class of instructions is the most used (the maximum of the r c for all c ∈ C). This information can be used to choose the best architecture for the algorithm. For example an algorithm with a lot of branching operations will have poor performance if embedded on GPU, but could have good performance on ARM, with speculative execution.

Moreover, if we know the throughputs of different classes of instructions for a given architecture, we are able to estimate computation time for each class of instructions on this architecture. Let p c,a be the throughput (in operations per cycle) of architecture a for class of instructions c, computation time (in cycles) T c,a is:

T c,a = N c p c,a . (7)
If architecture a is scalar, different instructions of different classes cannot be executed simultaneously, so total computation time is the sum of all T c,a :

t a = i∈C T i,a . (8)
If architecture a is superscalar, and in the optimal case (when instructions of different classes follow each other with good timing in order to get the maximum benefit from the superscalar capability of the architecture), computation time is given by the maximum of all T c,a :

t a = max {i∈C} T i,a . (9)
State-of-the-art architectures are all superscalar, thus we can estimate an execution time interval for the kernel associated with architecture a. In the best case (when the superscalar capability is fully exploited), the computation time is given by the maximum of all T c,a . In the worst case (when the superscalar capability is not exploited), the total computation time is the sum of all T c,a :

t min,a = max {i∈C} T i,a t max,a = i∈C T i,a . (10)
Given a kernel and a computing unit, we are able to estimate a computation time interval. In our case, that can be helpful to estimate the best computing unit for a given kernel. Indeed, our method is applicable to different kind of computing units, and results can be compared by knowing clock frequency of each computing unit.

C. Prediction for Consecutive Kernels

With computing profile prediction, we can obtain an interval of predicted execution time for a given kernel and different computing units. However a complete algorithm is composed of consecutive kernels, so we have to associate predictions of multiple kernels in order to predict performance of a complete algorithm.

First, we can simply use interval arithmetic, to sum intervals of the different kernels in order to obtain another interval. This implies a loss of precision because of a larger interval. Let an algorithm be composed of two consecutive kernels k 1 and k 2 , let x 1 , x 2 be the real execution time of

k 1 and k 2 on a given computing unit. Let [a 1 , b 1 ], [a 2 , b 2 ]
be the predicted interval obtained with computing profile for k 1 and k 2 , such as

x 1 ∈ [a 1 , b 1 ] and x 2 ∈ [a 2 , b 2 ].
This approach returns another interval for the predicted total execution time:

(x 1 + x 2 ) ∈ [a 1 + a 2 , b 1 + b 2 ].
The exact value of execution time is unknown (until we measure it in real condition), so this value belongs to the predicted interval [t min , t max ]. Without any other specific information, we can modelize the possible value for the execution time as a uniform distribution over the predicted interval [t min , t max ]. As predictions of two different kernels are independent (as long as they are not executed concurrently), summing two predicted execution times implies summing two independent random variables.

Let us now consider two independent random variables X 1 and X 2 representing the estimated execution times, and two uniform distributions f X1 (p) and f X2 (p), describing the relative likelihood that X 1 = x 1 and X 2 = x 2 , defined such as:

f Xi (p) = 1 bi-a1 for a i ≤ p ≤ b i , 0 for p < a i or p > b i , i ∈ {1, 2}.
(11) Thus:

P r[p 1 ≤ X i ≤ p 2 ] = p2 p1 f Xi (p)dp. (12)
As X 1 and X 2 are independent, the probability density function of their sum, f X1+X2 , is given by the convolution of the two density functions:

f X1+X2 (p) = ∞ -∞ f X1 (p -u)f X2 (u)du. ( 13
)
The density function

f X1+X2 (p) is non-zero for p ∈ (a 1 + a 2 , b 1 + b 2 ) and is maximum for p ∈ [min(a 1 + b 2 , a 2 + b 1 ), max(a 1 + b 2 , a 2 + b 1 )].
For both approaches, predicted interval become larger each time we add a kernel. However the probabilistic approach enables to reduce the estimated interval with a degree of confidence, as shown in ( 12), thus we obtain a better precision for performance prediction.

VI. CASE OF STUDY: A LANE DETECTION ALGORITHM EMBEDDED ON THE K1 SOC

To illustrate our method for performance prediction, we take an automotive use-case: lane detection application. The image processing part of our application takes about 95% of computational needs. It is divided into four levels:

1) Gradient computation: the input color image is converted to grayscale, the gradient is computed with a horizontal Sobel operator, and finally a threshold is applied in order to get a binary image. 2) A Bottom-Hat filter is applied to reduce the number of points for the next step.

3) The lane detection is performed by a Hough transform. 4) High level kernel for decision-making. The gradient computation and the Bottom-Hat filter outputs are illustrated in Fig. 5.

A. Throughputs of the K1 SoC

As seen in previous section, the model is based on throughputs of architectures, p c,a . Manufacturers may give information about throughput for each class, e.g. for Nvidia this can be found in [START_REF]Cuda C programming guide[END_REF].

Moreover, we establish throughputs for different computing units by using a benchmark, which measures the time to execute a known number of instructions. Time is measured with performance counter registers of each architecture. Throughputs for ARM and GPU of the K1 SoC are given in Tab. I. As Specific operations on ARM imply multiple instructions, there is not a unique throughput for this class. Throughputs are given for one core, or one SMX for the GPU. Speed increasing provided by using 4 cores on ARM vary depending on parallelism degree of kernel and API used. For a P 0 kernel with high Ai, we assume that using the best acceleration technique provides a speed increasing equal to the number of threads (e.g. using 4 cores reduces execution time by 4).

The throughputs of NEON operations on ARM are the same than corresponding regular operations, e.g. NEON unit can execute 2 SIMD additions per cycle. However using NEON implies that data have to be loaded into NEON registers from ARM registers, and have to be stored from NEON registers to ARM registers, throughput for loading and storing NEON instructions can vary depending on data alignment.

B. Example of Performance Prediction with Horizontal Gradient

In order to illustrate our methodology for performance prediction, we present an example with the Gradient algorithm of our lane detection application. The algorithm processes in two steps, first a RGB to gray conversion and then a horizontal Sobel filter. In order to achieve high Ai, we choose to store temporary values in a small buffer, benefiting from low latency of L1 cache. Thus complete algorithm performs only 3 loads (one for each color component) and 1 store from / to global memory for each pixel. RGB to gray conversion is given by listing 1, with pIn pointer on input image, and pOut pointer on output image. For each pixel, it performs 3 M Int, 5 S Int plus 2 if we consider that multiplication-accumulation is not used, and 1 S Int plus 1 Branch for the for loop. We choose to separate the Sobel filter into two filters. The first one is given by listing 2, for each pixel it performs 5 S Int, and 1 S Int plus 1 Branch for the for loop. The second one is given by listing 3, the Step variable is the distance between vertically consecutive pixels. For each pixel it performs 7 S Int (multiplication by 2 is counted as S Int, because corresponds to an arithmetic left shift), 1 S Int plus 1 Branch for the if condition, 1 S Int for the abs operator , and 1 S Int plus 1 Branch for the for loop. 1) Prediction on ARM: First, we apply prediction model for ARM using one core and no vectorization. We assume Ai is high enough to consider performance is limited by computing capability of the architecture, not by memory bandwidth. Predicted execution time, t, is given in number of clock ticks per pixel:

• RGB to Gray: t ∈ [START_REF] Geronimo | Survey of pedestrian detection for advanced driver assistance systems[END_REF][START_REF]Nvidia kepler GK110 architecture whitepaper[END_REF] • First Sobel Filter: t ∈ [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF][START_REF] Geronimo | Survey of pedestrian detection for advanced driver assistance systems[END_REF] • Second Sobel Filter: t ∈ [START_REF] Sankaran | TDA2X, a SoC optimized for advanced driver assistance systems[END_REF][START_REF] Lanier | Exploring the design of the cortex-A15 processor[END_REF] The global execution time, obtained by summing the intervals, is t ∈ [START_REF] Naishlos | Autovectorization in GCC[END_REF][START_REF] Hong | An analytical model for a GPU architecture with memory-level and thread-level parallelism awareness[END_REF], or t = 15.5 ± 22.6% clock ticks per pixel.

The probability density distribution f (p), obtained by convolving probability density distribution of each kernel, is: [START_REF] Lopez-Novoa | A survey of performance modeling and simulation techniques for accelerator-based computing[END_REF][START_REF] Ubal | Multi2sim: a simulation framework for CPU-GPU computing[END_REF]. According to the density function, the probability that execution time t ∈ [START_REF] Zhou | Task mapping in heterogeneous embedded systems for fast completion time[END_REF][START_REF] Kerr | Modeling GPU-CPU workloads and systems[END_REF] (or t = 15, 5 ± 9.7%) is 70%.

f (p) =                             (p -12) 2 /
We implement Gradient algorithm on ARM, and apply it on a 1280 × 720 image. Measured execution time is 9.750 ms, 15.87 clock ticks per pixel which is in the interval with 70% confidence, and represents a difference of 2.5% compared with the average of the interval.

To benefit from multi-core and decrease execution time, image is divided into four parts, each thread processes one part of the image. In order to avoid side effect, each thread processes their parts plus one line, in our case this represents blocks of 1280×181 pixels. Based on prediction with 1 core, we can estimate execution time using 4 cores. As all kernels used for Gradient algorithm are P 0 , and pthread is used, we assume there is no waste of time.

One thread processes on 1280×181 pixels, thus predicted time is t ∈ ( [START_REF] Naishlos | Autovectorization in GCC[END_REF][START_REF] Hong | An analytical model for a GPU architecture with memory-level and thread-level parallelism awareness[END_REF] × 1280 × 181) clock ticks. If we divide by the total number of pixels, then t ∈ [3.017, 4.776] clock ticks per pixel. Measured value is 2.560 ms, 4.17 clock ticks per pixel, it represents a difference of 7% compared with the average of the interval.

Gradient algorithm can also be accelerated by using NEON. Performance with auto-vectorization depends on compiler capabilities, so our methodology using computing profile can not work. Moreover, in this example memory delay is hidden by computing time. In order to address prediction where performance depends on other parameters than computing profile and throughput, we are using the second level of our estimation approach (a generic benchmark-vectors set injected on this specific architecture). The extracted parameters can be acceleration provided with auto-vectorization, memory access delays, transfers delays between computing units, etc. Then, both information from computing profile and benchmarks are used to obtain predicted performances (cost function f ) and the best mapping for a given algorithm and heterogeneous SoC.

2) Prediction on GPU: We apply the same methodology on the same algorithm for the GPU of K1. However, in order to maximize occupancy [START_REF] Harris | Optimizing cuda[END_REF], each thread computes four pixels, and Sobel filter is not separated. Moreover RGB to Gray conversion uses floating point operations. Temporary values are stored in shared memory in order to minimize the number of global memory access.

With CUDA API, image indexes are obtained with threads IDs, there is no need of for loop. Thus in RGB to Gray kernel, each thread performs 12 Float, 32 S Int and 1 Branch (to handle borders of blocks). In Sobel kernel, each thread performs 29 S Int, 4 Branch, and 4 Specific (absolute value operator belongs to Specific class for the GPU). Thus prediction, in clock ticks per pixel, for each kernel are: Measured execution time on a 1280 × 720 image is 0.606 clock ticks per pixel, 990 Mp/s, which is in the interval with 63% of confidence, and represents a difference of 9.66% compared with the average of the interval. Fig. 8 shows results of predictions and measured execution times for both ARM 4 cores and GPU with different image sizes.

C. Kernel Mapping

Our application is divided into four parts and K1 SoC has two programmable computing units, thus more than 128 different configurations can be explored (different mapping, of core to use, execution pipeline, etc). On TDA2x, there are four programmable computing units, so more than 1024 different configurations can be studied. Our methodology aims to automatically find one of the best mapping without porting any kernel. On K1, predictions show that computing time is slower on GPU than on ARM for all kernels except for the high level decision making. Acquisition is handle by the ARM, thus executing a kernel on GPU implies copying data to GPU memory. According to our benchmarks (see Fig. 9), copying a RGB image and running the kernel on GPU is more time consuming than running it on ARM. The Hough transform kernel only need few pixels coordinates as input data, returns lines coordinates, and prediction shows it has good performance on GPU. Thus one of the best mapping for our application is to embed Hough transform on GPU and everything else on ARM with 4 cores. With this mapping, the application runs in 10.4 ms on a 1920×1080 HD image.

On TDA2x, we are working to embed lane detection application. Our methodology gives us that one of the best configuration is to run the two first kernels on DSP, then Hough transform on EVE, and finally high level decision making on ARM.

VII. CONCLUSION AND FUTURE WORK

In this work, we have introduced a novel approach for performance prediction addressing multi-architectures, and operating without specific and optimized source code of application for each architecture. Our methodology is based on a high level description of kernels, the computing profile, and computing throughputs of architectures. We have shown, with an example extracted from an ADAS application, our approach is able to predict a more or less wide interval of execution time with a degree of confidence.

We are working to improve our methodology by taking into account memory delay to predict performances of ker-nels with small arithmetic intensity. Thus, we are adapting the approach of boat hull model to our methodology, classifying type of memory access found in image processing algorithms and studying there behaviors. Moreover, we are applying our method to TDA2x SoC to confirm the general approach for computing units like vectorial processor.

In a future work, we will propose results on kernels which belong to other parallelism classes than P 0 . Finally we will show performance prediction using both computing profile and parameters extracted with our set of benchmarks, in order to predict all terms of (3). Thus, we aim to deal with performance prediction of a complex application for different heterogeneous SoC, to help with kernel mapping optimization, and to choose the best suited SoC for a given application.

Figure 1 .

 1 Figure 1. ARM A15 3-way superscalar instruction pipeline.

Figure 2 .

 2 Figure 2. Tegra K1 heterogeneous architecture, global memory is the same for CPU, GPU and ISP.

Figure 3 .

 3 Figure 3. Block diagram of TDA2x heterogeneous architecture.

Figure 4 .

 4 Figure 4. Execution time for an image copy kernel for different image size, using 1, 2, 3, and 4 cores of the ARM.

  Let N c be the number of instructions associated with class c, with c ∈ C. The ratio for the class c, r c is:

Figure 5 .

 5 Figure 5. Successive steps of the lane detection application, and their computing profile.

Listing 1 .

 1 RGB to Gray conversion * ( pOut ++) = ( * ( p I n ) * 28 + * ( p I n + 1 ) * 151 + * ( p I n + 2 ) * 7 7 ) >> 8 ; p I n +=3;

Listing 2 .

 2 First Sobel filter * ( pOut ++) = * ( pIn -1) - * ( p I n + 1 ) ; p I n ++;

Listing 3 .

 3 Second Sobel filter v a r = a b s ( * ( pIn-S t e p ) + * ( p I n ) * 2 + * ( p I n + S t e p ) ) i f ( v a r > 2 5 5 ) * ( pOut ++) = 2 5 5 ; e l s e * ( pOut ++) = v a r ; p I n ++;

Figure 6 .

 6 Figure 6. Probability density function of predicted execution time for gradient algorithm using 1 core and no vectorization on ARM. Measured execution time is 15.87 clock ticks per pixel.

Figure 8 .Figure 9 .

 89 Figure 8. Measured and predicted performances for the Gradient computation on ARM using 4 cores (top) and GPU (bottom). Both axis are in log scale.

Table I THROUGHPUTS

 I FOR ARM AND GPU OF THE K1 ARCHITECTURE. VALUES ARE GIVEN IN INSTRUCTIONS PER CYCLE.

	Class of instructions	ARM A15	GPU Kepler
	S Int	2	160
	M Int	1	32
	Float	1	192
	Specific	*	32
	Branch	1	32
	Address	1	32
	NEON load & store	0.5	/

•

  RGB to Gray: t ∈ [0.2, 0.29275] • Sobel: t ∈ [0.18125, 0.43125] Predicted global execution time is t ∈ [0.38125, 0.724] clock ticks per pixel, or t ∈ [829, 1574] Mp/s (clock rate of K1 GPU is 600 MHz). According to the probability density function, given in Fig. 7, t ∈ [0.474, 0.63125] clock ticks per pixel with 63% of confidence.

  Figure 7. Probability density function of predicted execution time for gradient algorithm using GPU. Measured execution time is 0.606 clock ticks per pixel.

						GPU
	6					Probability density function f(p)
	5
	4
	3
	2
	1
	0 0.35	0.4	0.45	0.5	0.55	0.6	0.65	0.7	0.75
				Execution time in clock ticks per pixel
				Execution Time for the Gradient Computation on ARM
		10000		Reality t min t max
	Execution Time in us	1000
		100 300k	400k		600k	800k	1M
						Number of Pixels
				Execution Time for the Gradient Computation on GPU
		10000		Reality t min t max
	Execution Time in us	1000
		100 300k	400k		600k	800k	1M
						Number of Pixels


Romain Saussard

Boubker Bouzid
email: boubker.bouzid@renault.commarius.vasiliu

Marius Vasiliu

Roger Reynaud
email: roger.reynaud@u-psud.fr

THE EMBEDDABILITY OF LANE DETECTION ALGORITHMS ON HETEROGENEOUS ARCHITECTURES

Keywords: ADAS, Lane Detection, Embedded Processing

Lane detection plays a crucial role for Advanced Driver Assitance System (ADAS) or autonomous driving applications. Literature shows a lot of lane detection algorithms can work in real time with good results. However, they require much computer processing and cannot be embedded in a vehicle ECU without deep software optimizations. In this paper, we discuss the embeddability of lane detection algorithms by comparing state-of-the-art algorithms in terms of functional performance and computational timing. We identify what essential parts of lane detection are time consuming, and show these parts can be computed in real time on embedded systems.

INTRODUCTION

On-vehicle lane detection system is an important component for ADAS application. In fact it can be use for warning systems (e.g. alerting the driver in case of dangerous situation: Lane Departure Warning) or active systems (Lane Keeping Assist or Lane Centering Assist).

Lane detection has been deeply studied this past twenty years [START_REF] Bertozzi | Gold: A parallel real-time stereo vision system for generic obstacle and lane detection[END_REF][START_REF] Mccall | Video-based lane estimation and tracking for driver assistance: survey, system, and evaluation[END_REF][START_REF] Kim | Robust lane detection and tracking in challenging scenarios[END_REF][START_REF] Sivaraman | Integrated lane and vehicle detection, localization, and tracking: A synergistic approach[END_REF]. However, embedded lane detection is not much studied despite the fact car manufacturers and suppliers are interested in this kind of works. Some automotive providers already offer embedded solutions for lane detection. They produce intelligent cameras which perform embedded real time lane detection, the algorithm is computed by FPGA/ASIC in camera. However, state-of-the-art algorithms show better results.

Semiconductor companies are moving into the market of embedded systems for ADAS with heterogeneous architectures. Heterogeneous architectures embed several processing units with different capabilities on the same System on Chip (SoC), often with massively parallel computing unit. We can cite the Tegra K1 SoC of Nvidia (which embed ARM, GPU and an ISP), the TDA2x SoC [START_REF] Sankaran | TDA2X, a SoC optimized for advanced driver assistance systems[END_REF] of Texas Instrument (ARM, DSP and EVE vectorial processor). This type of SoC offers a good solution for embedding image processing algorithms in the automotive field because it can handle high performance computing with low-power consumption.

In this paper, we discuss the embeddability of lane detection algorithms on heterogeneous architectures. First we will study the different approaches for lane detection, their results and identify their bottlenecks. Then we will discuss how they can be adapted for embedded processing and what performances can be expected with the different approaches.

LANE DETECTION ALGORITHMS

The aim of lane detection algorithms is to find positions of lanes on the roads and to characterize them. Based on [START_REF] Hillel | Recent progress in road and lane detection: a survey[END_REF], a lane detection algorithm can be separated into four parts:

• Image cleaning: This step enhances the input image in order to get better results with features extraction. It can be an algorithm dealing with abrupt illumination change (e.g. when vehicle goes under a tunnel), color space conversion, ROIs extraction, etc.

• Features extraction: Low level features are extracted from the image. Algorithms used for this step can be gradient computation, color segmentation, thresholding, Inverse Perspective Mapping (IPM). The IPM (also called bird-eyes view) consists in applying a homography on the image in order to obtain parallel lines.

• Lane model fitting: A geometric model of lane (2D or 3D) is used on extracted features in order to find potential lanes. Hough transform (straight lines or curves) is often used in this step.

• Temporal integration: Results of previous frames are used in order to detect lane in current frame. It can be used to reduce the search area in the image and the computational cost. Temporal integration is not always used, and often implies data fusion with other sensors, e.g. IMU, GPS, etc.

In spite of some datasets can be found for pedestrian detection [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF], to our knowledge, no free dataset with ground truth can be found for lane detection algorithms. That is why it is hard to compare results of state-of-the-art lane detection algorithms.

However, some papers give results and performance of their algorithm. Thus in [START_REF] Lu | An improved approach for vision-based lane marking detection and tracking[END_REF], where both IPM and Hough transform are used, algorithm has been tested on KITTI database [START_REF] Andreas | Vision meets robotics: The KITTI dataset[END_REF] (1392 × 512 color images) with an average of 68.18 ms processing time per frame. The average error between marked points and lines located with GPS + map is 15 cm.

The algorithm in [START_REF] Aly | Real time detection of lane markers in urban streets[END_REF] performs an IPM transformation with integration on Y axis for features extraction and uses RANSAC for lane fitting. It achieves 50 frames per second rate with 640 × 480 input images. Their provided correct detection rate is 90.89% in urban conditions.

In [START_REF] Wang | Lane detection and tracking using B-Snake[END_REF], the lane fitting is based on active contours with initialization performed by Canny & Hough lines detection. It shows good results, 95% of correct lane detection, but poor performance, only 2 frames per seconds with a resolution of 240 × 256. We can notice it has been tested on old CPU, executing the algorithm on a recent one may show better performance.

Most state-of-the-art algorithms use IPM and/or Hough transformation. These two kernels are essential to have good results for lane detection, but they have a high computational cost. That is why we choose to focus on these two kernels, to implement and optimize them on an embedded SoC.

HETEROGENEOUS ARCHITECTURES

Semiconductors companies such as Nvidia, Texas Instrument and Freescale propose heterogeneous architectures to meet automotive industry needs for embedding image processing algorithms for ADAS. These architectures handle high performance computing with massively parallel computing units (e.g. GPU or vectorial processor) and low-power consumption. In this paper, we propose a general approach which can fit with any of these SoCs, and we show some preliminary results obtained on Nvidia Tegra K1 heterogeneous SoC.

Nvidia Tegra K1 SoC

The Nvidia Tegra K1 SoC is composed of a quad-core ARM Cortex A15 CPU (1.5 GHz clock rate) providing ARMv7 instruction set and an out-of-order speculative issue 3-way superscalar execution pipeline, each core has a NEON and FPU unit [START_REF] Lanier | Exploring the design of the cortex-A15 processor[END_REF].

The Kepler GPU of the K1 is composed of one streaming multiprocessors of 192 cores [START_REF] Nvidia | Nvidia kepler GK110 architecture whitepaper[END_REF], accessible with CUDA [START_REF]Cuda C programming guide[END_REF] and the new standard OpenVX [START_REF] Rainey | Addressing system-level optimization with OpenVX graphs[END_REF]. The parallelization with CUDA is qualified as SIMT (Single Instruction Multiple Threads). The global memory of the K1 is the same for GPU and ARM, both can potentially access to the same data. Memory area of each processing unit is handled by the OS.

K1 Specific Features

The K1 also has several hardware processing units like Image Signal Processor (ISP) and Image Processing Accelerator providing fast and specific image processing algorithms such as debayering, noise reduction, lens correction etc. These units are very fast but the user can only control a limited set of parameters and the chaining order of kernels is restricted.

In addition of the basic instructions, ARM and Kepler GPU have specific features which can be used to accelerate image processing algorithms. Thus, ARM processor provides SIMD instructions with NEON units [START_REF] Mitra | Use of SIMD vector operations to accelerate application code performance on low-powered ARM and Intel platforms[END_REF]. In order to handle the issues of reading and writing accesses of different threads, CUDA provides atomic instructions based on hardware design. This feature is very useful for histogram construction and the voting process of Hough transform.

The Nvidia GPU texture memory provides couple of advantages for reading images data. When accessing to noninteger pixel coordinates (e.g. Im[i-0.5][j-0.5]), texture units handle linear interpolation, this is cost-free because computed by hardware design. Moreover texture units handle access to pixels outside the image (e.g. Im[-1][-1]), also cost-free. This feature can highly decrease the IPM computational timing.

KERNEL MAPPING AND EMBEDDABILITY

According to section 2, IPM and Hough transform are two kernels used by state-of-the-art algorithms with best results, but have a high computational cost. Embedding a given algorithm on a heterogeneous architecture is a difficult task because one can not easily find how to allocate kernels on the different processing units (kernel mapping) [START_REF] Zhou | Task mapping in heterogeneous embedded systems for fast completion time[END_REF].

Kernel Mapping Optimization

Let P be the vector of the processing units of a given heterogeneous architecture, K be the vector of kernels of a given algorithm, the matrix M constitutes the mapping of K on P , given by P = M.K. Let ϕ be the dependency matrix (dependencies between kernels K), τ (M ) be the execution function (returning the execution time for each kernel), δ(M ) be the transfer function (returning the transfer delay needed for each kernel), η(M ) be the occupancy function (returning the occupancy for each computing unit) and f (P, K, ϕ, τ (M ), δ(M ), η(M )) be the cost function (returning the global execution time of the algorithm). The aim of kernel mapping optimization is to find M minimizing f :

arg min M [f (P, K, ϕ, τ (M ), δ(M ), η(M ))]
Parameters of function f can be measured or predicted for different M . Measure needs all kernels implementation on all processing units (very time consuming), but prediction needs kernels and architectures analysis (very efficient but needs deep SW / HW knowledge). Actually, our approach mixed both techniques, starting with some specific representative kernels implementations, and finishing with performance prediction analysis to achieve the best kernel mapping [START_REF] Saussard | Predicting ADAS algorithms performances on K1 architecture[END_REF].

Gradient Computation

Gradient computation is often performed with convolution operators (e.g. Sobel filter). First, it can easily be accelerated by using separable filter. Moreover convolution operators are easily parallelizable because the result of one pixel is not dependent on other pixels result. Thus it can be accelerated with GPU or by using NEON instructions and multicore on CPU ARM.

According to our results, an algorithm performing horizontal and vertical Sobel filter is 14 times faster by using 4 cores and 8×int16 NEON vectors than using only one core. Our implementation of this algorithm is 20 times faster on the GPU of K1 than the implementation on ARM with one core. However, using GPU implies extras delays for transferring data between CPU and GPU (about 2ms for a 2MB image).

Inverse Perspective Mapping (IPM)

Many state-of-the-art lane detection algorithms use IPM, it results an image with vertical and parallel lines which reduces the complexity of lane marking. Camera intrinsic (focal length and optical center) and camera extrinsic (pitch angle, yaw angle and height above the ground) parameters are needed to construct the homography matrix [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. Pitch angle and yaw angle may not be constant (e.g. vehicle vibrations), that is why camera is often associated with an IMU [START_REF] Lu | An improved approach for vision-based lane marking detection and tracking[END_REF].

A pixel I IP M of the bird-eyes view image with coordinate (u, v) corresponds to the pixel I in with coordinate (x, y) of the input image, so I IP M (u, v) = I in (x, y). Coordinates are linked by: s.x s.y s T = H. u v 1 T , with H the homography matrix. However, x and y will be non-integer values, so interpolation should be performed in order to compute I in (x, y). Despite the complexity of the IPM is linear, this transformation is associated with some computational cost and is often apply on a ROI. Linear interpolation can be handled very quickly with texture memory of the GPU which highly reduces the computing time. Indeed only the matrix product has to be computed for each pixel, so it can be parallelized.

Our IPM implementation on GPU provides 5 ms on a 1920 × 1080 HD image (with a transfer delay 8.5 ms). These results on HD image are here to illustrate reachable performance, real implementation uses a ROI of a smaller image.

Hough Transform

The computational complexity of a two-dimensional Hough transform is given by: O(n.s), where n is the number of input points and s the number of iterations needed. Based on [START_REF] Van Den Braak | Fast hough transform on GPUs: Exploration of algorithm trade-offs[END_REF], Hough transform computation can be accelerated on GPU.

In fact reducing the number of input points and parallelizing each iteration decreases the computational timing. We implemented a simple algorithm which performs a Sobel filter, a thresholding, a morphological filter and apply the Hough transform. According to our results, the best mapping is to embed Hough transform on GPU, and other kernels on ARM. Our algorithm has been tested with a 1920 × 1080 HD image, about 0.2% of total pixels are used for Hough transform. The ARM processes in 5.4 ms, and the Hough transform on GPU in 4.5 ms. Thus, by adding the memory delay and latency for GPU kernel execution, the full lane detection algorithm is executed in 10.4 ms. Results for tested algorithms are given in table 1.

Our next hardware target for studying ADAS algorithms embeddability is the TDA2x platform [START_REF] Sankaran | TDA2X, a SoC optimized for advanced driver assistance systems[END_REF]. Gradient will probably be computed by the DSP, and IPM and Hough transform by the EVE vectorial processors.

Conclusion

The embeddability of image processing algorithms is a big issue for automotive industry, it is not a simple task to predict if a robust known algorithm can be embedded on a given SoC. In this paper we discussed the embeddability of lane detection algorithms on a heterogeneous architecture, the K1 SoC. It has been shown that costly operations such as IPM or Hough transform can be optimized and executed in real time on an embedded SoC. However, doing such studies for multiple algorithms and SoCs represents a lot of work. The most effective solution for dealing with embeddability of image processing algorithms is performance prediction. Thus, we are working on performance prediction for algorithms embedded on different SoCs, in order to help car manufacturers and suppliers to choose what algorithm -SoC association presents the best performance and efficiency.

Table 1 .

 1 Best performances for tested algorithms, execution time is given in pixels per second. As Hough transform performance depends on the number of input points, which can vary with a fixed image size, the value is given in points per second.

	Algorithm	Performance Computing Unit
	Horizontal & vertical	450 Mp/s	ARM
	gradient
	IPM transformation	415 Mp/s	GPU
	Hough transform	960 kPts/s	GPU
	Simple lane detection	200 Mp/s	ARM + GPU


Giuseppe Iurato
email: giuseppe.iurato@unipa.it

Giuseppe Ruta
email: giuseppe.ruta@uniroma1.it

The role of virtual work in Levi-Civita's parallel transport

According to current history of science, Levi-Civita introduced parallel transport solely to give a geometrical interpretation to the covariant derivative of absolute differential calculus. Levi-Civita, however, searched a simple computation of the curvature of a Riemannian manifold, basing on notions of the Italian school of mathematical physics of his time: holonomic constraints, virtual displacements and work, which so have a remarkable, if not dominant, role in the origin of parallel transport.

Introduction

Tullio Levi-Civita (1873-1941)'s teachers gave him a strong training in mathematical physics and mechanics, and developed absolute differential calculus, an inheritance of Beltrami's investigations on Riemannian manifolds [START_REF] Capecchi | Strength of Materials and Theory of Elasticity in 19th Century Italy[END_REF][START_REF] Tazzioli | [END_REF]. In the mid-1910s, Levi-Civita investigated relativity by absolute differential calculus: he studied the curvature of four-dimensional Riemannian manifolds, modelling space-time, through the parallel transport of vectors over these manifolds.

We usually find that Levi-Civita's parallel transport is motivated by the search for a geometrical interpretation to the covariant derivative of absolute differential calculus. However, in the introduction to [START_REF] Levi-Civita | Opere matematiche. Memorie e note[END_REF] Levi-Civita states that he wishes to simplify the computation of the curvature of a Riemannian manifold V n with dimension n ≥ 2, and introduces parallelism on it. Going further, the principle of virtual work seems to be one of the basis for parallel transport according to Levi-Civita.

Virtual work and parallel transport

Balance exists if and only if the total power (work) on any admissible, or virtual, velocity (firstorder displacement) field vanishes. This law is basically the same since Lagrange's formulation, where inertia is a force to be added to the active ones, and dynamics is reduced to statics. 1 The reactions of perfect constraints spend no work on admissible displacements; if δP i is the first-order displacement of the point where active forces plus inertia F i are applied, the law of virtual work is ( [START_REF] Lagrange | Mécanique analytique[END_REF], p. 14)

δL = i F i • δP i = 0, (1)
Levi-Civita in [START_REF] Levi-Civita | Opere matematiche. Memorie e note[END_REF] embeds V n in a Euclidean space S N with dimension N ≤ n(n + 1)/2, and considers two unit vectors α, α ′ emerging from nearby points P, P ′ ∈ V n . In S N , α, α ′ are parallel if, for any auxiliary direction f emerging from P , angle( α, f ) = angle( α ′ , f ). Parallelism in V n asks this relation to hold for any direction f of the tangent plane T S N P (V n ), which is an intrinsic definition, since it depends solely on metrics in V n , ds 2 = n i,k=1 a ik dx i dx k . The manifold V n may be described by the system ( [3], eq. ( 1), p. 4)

y ν = y ν (x 1 , . . . , x n ), ν = 1, 2, . . . , N (2)
with y ν ∈ S N , x n ∈ V n . Eq. ( 2) describes an n d.o.f. system subjected to N smooth holonomic bilateral constraints. The point P varies on a smooth curve C ∈ V n , parameterized by the abscissa s, thus the components α ν = α ν (s), so that

C ≡ y ν (s) = y ν (x 1 (s), . . . , x n (s)), ν = 1, . . . , N ⇒ y ′ ν = n i=1 ∂y ν ∂x i x ′ i ν = 1, 2, ..., N, (3)
In the analog constrained system, s is an evolution parameter (e.g., time), and C is a trajectory in the manifold of admissible configurations. Derivation with respect to s yields the direction cosines y ′ ν ∈ S N , while x ′ i are the direction cosines of the same unit direction in V n . Levi-Civita poses the direction cosines of α to be ξ (i) , i = 1, 2, . . . , n with respect to V n , and α ν , ν = 1, . . . , N with respect to S N . Then, from eq. ( 3) we have ( [START_REF] Levi-Civita | Opere matematiche. Memorie e note[END_REF], eq. ( 7), p. 6)

α ν = n l=1 ∂y ν ∂x l ξ (l) ν = 1, 2, . . . , N. (4)
Levi-Civita considers an arbitrarily fixed direction f of S N , with direction cosines f ν , thus

cos f , α = N ν=1 α ν f ν ⇒ d cos f , α = ds N ν=1 α ′ ν (s)f ν . (5)
Levi-Civita's intrinsic condition of parallelism requires the variation in eq. ( 5) to vanish only for the directions tangent to V n at P ∈ C, which are those compatible with the constraints (2) ( [3], p. 7). Thus, by replacing the f ν with quantities proportional to them, Levi-Civita's intrinsic definition of parallelism implies ( [3], eq. (I), p. 7)

N ν=1 α ′ ν (s)δy ν = 0 (6)
for any variation δy ν , that is, for any admissible virtual displacement compatible with the constraints in eq. ( 2). If the α ′ ν (s) are some mechanical actions in S N , eq. ( 6) is a formulation of the virtual work law in S N related to the smooth bilateral holonomic system defined by eq. ( 2), hence related to a Riemannian manifold. From eq. ( 2) it follows that ( [3], p. 7)

δy ν = n k=1 ∂y ν ∂x k δx k ν = 1, 2, . . . , N, ⇒ N ν=1 α ′ ν (s) ∂y ν ∂x k = 0 (k = 1, 2, . . . , n), (7)
for the arbitrariness of the δx k and (6). Eq. ( 7) are the parallelism conditions for the directions α moving along C. To have an intrinsic expression, Levi-Civita replaces the direction cosines α ν with their expression eq. ( 4), to include the intrinsic direction cosines ξ (i) , and, if Γ jl i are Christoffel symbols of second kind, he deduces ([3], eq. (I a ), p. 8)

dξ (i) ds + n j,l=1 Γ jl i x ′ j ξ (l) = 0 (i = 1, 2, . . . , n), (8)
Levi-Civita shows that eq. ( 8) may be expressed by the covariant quantities associated with the ξ (i)2 , thus the intrinsic condition of parallelism has the same form of Lagrange equations of motion on a Riemannian manifold ( [START_REF] Levi-Civita | Opere matematiche. Memorie e note[END_REF], eq. (I c ), p. 12). It seems apparent that one key guide to Levi-Civita is mathematical physics, not only pure geometry, and, in particular, the principle of virtual work for a mechanical system under smooth holonomic bilateral constraints.

Final remarks

We claim that the virtual work principle played a key role in the origin of Levi-Civita's parallel transport in a Riemannian manifold, eq. ( 6). Levi-Civita's jargon clearly refers to this, e.g., with the locutions 'constraint' and 'displacements compatible with constraint'. Furthermore, he specifies that eq. ( 6) is obtained "for all the displacements δy ν compatible with the constraints (2)", highlighting in italic this phrase in [START_REF] Levi-Civita | Opere matematiche. Memorie e note[END_REF]. Indeed, for the analog constrained mechanical system, the unit directions emerging from the points of V n assume the role of admissible displacements, and the dual forms on them are reactions provided by the geometrical links between the elements of the mechanical system. ( [START_REF] Sommerfeld | Vorlesungen über Theoretische Physik[END_REF], ch. II, sect. 8).

In his monograph on absolute differential calculus [START_REF] Levi-Civita | Lezioni di calcolo differenziale assoluto (A. Stock[END_REF] (ch. V, (b), sects. 10-15), Levi-Civita highlights the role of analytical mechanics in developing differential geometry notions, and says that his symbolic equation of parallelism formally recalls the principle of virtual work. Heuses the same principle in discussing the geodetic principle for the dynamics of a material particle moving in a four-dimensional space-time manifold ( [START_REF] Levi-Civita | Lezioni di calcolo differenziale assoluto (A. Stock[END_REF], ch. XI, sect. 12).

It seems apparent to us that the key notions of the Italian school of mathematical physics were a strong basis for Levi-Civita's non-Euclidean geometrical definitions.

Levi-Civita adopts the locution 'moment', traditional in the Italian school of mathematical physics of his time, defining a mechanical action dual to a Lagrangean parameter of admissible (virtual) displacements.