Tracking of Vehicle Trajectory by Combining a Camera and a Laser Rangefinder

Keywords: Visual tracking, particle filter, sensor fusion

This article presents a probabilistic method for vehicle tracking using a sensor composed of both a camera and a laser rangefinder. Two main contributions will be set forth in this paper. The first involves the definition of an original likelihood function based on the projection of simplified 3D vehicle models. We will also propose an efficient approach to computing this function using a line-based integral image. The second contribution focuses on a sampling algorithm designed to handle several sources. The resulting modified particle filter is capable of naturally merging several observations functions in a straightforward manner. Many trajectories of a vehicle equipped with a kinematic GPS 1 have been measured on actual field sites, with a video system specially developed for the project. This field input has made it

Introduction

We will be presenting a vehicle trajectory tracking method that relies upon a sensor composed of a color camera and a one-dimensional scanning laser rangefinder. The objective of this system is to accurately estimate the trajectory of a vehicle traveling through a curve. The research reported herein lies within the scope of an French ANR-PREDIT project 2 The sensor, installed in a curve, is composed of three cameras placed on a tower approximately 5 m high to cover the beginning, middle and end of the curve, in addition to a scanning laser rangefinder laid out parallel to the ground. Since the cameras offer only limited coverage, their observations do not overlap and we will be 2 Y. Goyat et al. considering in the following discussion that the system can be divided into three subsystems, each composed of a camera-rangefinder pair, with the recalibration between each pair performed by means of rigid transformations, which will not be addressed in the present article. The object-tracking procedure is intended to estimate the state of an object at each moment within a given scene, based on a scene observation sequence. Tracking methods can be broken down into two major categories: the first concerns off-line or non-causal tracking, for which the state estimation at a given point in time uses the entire observation sequence [START_REF] Fleuret | Multicamera people tracking with a probabilistic occupancy map[END_REF]. The second category relates to online or causal tracking, for which the state of the object at a given point in time has been estimated as a function of the record of past and current observations and the record of past states [START_REF] Isard | Bramble: A bayesian multiple-blob tracker[END_REF]. This second category may also encompass the notion of realtime when the period necessary to estimate a state is shorter than the sensor acquisition frequency. The method described herein would be classified as a causal method.

In this article, we are proposing a solution based on a probabilistic formalization of the tracking problem with a stochastic framework (particle filter). The literature, which contains many references in the areas of vision and data merging, proceeds with the recursive time estimation of a state through application of Monte-Carlo methods [START_REF] Arulampalam | A tutorial on particle filters for on-line non-linear/nongaussian bayesian tracking[END_REF].

The rest of the article has been divided into four parts. The first will discuss the tracking principle, on the basis of a probabilistic model. The next part will focus on the likelihood functions proposed for estimating the weights associated with the particle set. The third part will provide a detailed description of our proposed sampling method (called multi-source sampling). The last section will then present the measurement campaign conducted for the purpose of quantifying method accuracy and robustness. A large number of trajectories could be estimated and a kinematic GPS was used to determine the actual field values associated with each trajectory.

Related Work

Recently, traffic video surveillance has become an important topic in the Intelligent Transport Systems (ITS), so vehicle detection, description, and/or recognition have been an active research field. Most of the solutions assume that the camera is at a high angle. Tracking vehicles from a static camera is challenging for several reasons:

-Outdoor vision tracking solutions must handle with variation of illumination and with shadows.

-Generic model based approaches are complex solutions due to the huge appearance variability of the vehicle class.

-Tracking several vehicles simultaneously implies to be able to handle with mutual occlusions.

Tracking objects from a static camera often uses a background/foreground subtraction [START_REF] Magee | Tracking multiple vehicles using foreground, background and motion models[END_REF]. In outdoor environment such method must handle with variations of illumination conditions and shadow. In [START_REF] Stauffer | Learning patterns of activity using real-time tracking[END_REF], Stauffer and al. propose a parametric model with a temporal update.

The resulting solution deals efficiently with small illumination variations. Detecting shadow can be done by color based analysis [START_REF] Ivanov | Fast Lighting Independant Background Subtraction[END_REF] and injecting temporal information into the models [START_REF] Mikik | Moving shadow and object detection in traffic scenes[END_REF].

Many object descriptions can be used in order to obtain a model of a vehicle. Some approaches are based on a generic model, uses to detect, then track the vehicle, from a luminance [START_REF] Avidan | Support vector tracking[END_REF] or Haar Wavelets based model [START_REF] Chateau | Real-time tracking with classifiers[END_REF][START_REF] Schneiderman | A statistical model for 3d object detection applied to faces and cars[END_REF]. 3D wireframes models have been also used but they require to define many models associated to different types of vehicles [START_REF] Zhao | Car detection in low resolution aerial image[END_REF][START_REF] Kim | Fast vehicle detection with probabilistic feature grouping and its application to vehicle tracking[END_REF][START_REF] Koller | Model based object tracking in monocular image sequences of road traffic scenes[END_REF][START_REF] Tan | Efficient image gradient based vehicle localization[END_REF][START_REF] Haag | Combination of edge element and optical flow estimate for 3d-model-based vehicle track-ing in traffic image sequences[END_REF]. In [START_REF] Pece | Tracking with the em contour algorithm[END_REF] a non-rigid 3D model is used into a EM based contour tracking. Recently, Kanhere and al. [START_REF] Kanhere | Vehicle segmentation and tracking from a low-angle off-axis camera[END_REF] have proposed an interest point based method to segment and track vehicles.

The method presented here uses a 3D simplified model of a vehicle, projected into the image, and then compared to the background/foreground subtraction map to provide and efficient observation function.

In [START_REF] Kamijo | Vehicle tracking in low-angle and front view images based on spatiotemporal markov random fields[END_REF], Kamijo and al.

Particle Filter

Vehicle trajectory is estimated recursively using a nonlinear filter, whose resolution entails a widespread stochastic method for vision applications: the particle filter. This choice is dictated by the nonlinear nature of the system.

Particle filtering [START_REF] Arulampalam | A tutorial on particle filters for on-line non-linear/nongaussian bayesian tracking[END_REF][START_REF] Isard | Condensation -conditional density propagation for visual tracking[END_REF] is based on an estimation of the a posteriori probability density p(X t |Z 0:t ) of state X t conditioned by the historical record of measurements Z 0:t , at time t, by a set of N weighted particles {(X n t , π n t )} N n=1 with their associated weights. In the case of an observa-tion stemming from several sources, we are proposing to merge observations intrinsically during the re-sampling phase. The algorithm derived (see Algorithm 1) is a variant of the Condensation algorithm [START_REF] Isard | Condensation -conditional density propagation for visual tracking[END_REF]. A weight vector (composed of weights generated from observations of each source) is to be associated with each particle. This multi-source re-sampling method (called MSS) will be further developed in Section 5, page 10.

Algorithm 1 Condensation in the multi-source case

Init : particles {(X ′ n 0 , 1/N )} N n=1 according to the initial distribution X 0 for t = 1, ..., T end do Prediction :

generation of {(X n t , 1/N )} N n=1 from p(X t |X t-1 = X ′ n t-1 )
Observation : estimation of the weight vector according to the various sources {(X n t ,

π n t )} N n=1 with π n t ∝ p(Z t |X t = X n t ) Sampling : build {(X ′ n t-1 , 1/N )} N n=1 from {(X n 0 , π n 0 )} N n=1 using Multi Source Sampling (MSS) Estimation : Xt . = 1 N N n=1 X n t
end for Output : The set of estimated states during the video sequence { Xt } t=1,...,T end Implementing a particle filter necessitates defining three models: 1) a state model that serves to define the kinematic characteristics of the object to be tracked;

2) an evolution model that defines the state of an object at a given point in time depending on the state at the previous point in time; and 3) an observation model that defines a measurement between a state hypothesis and the observations. The state and evolution model are presented in the following sections and the observation model is detailed in section 4.

Background Extraction

The majority of methods for tracking objects within a static scene introduce a background extraction step, which consists in a binary case of ascribing each image pixel a "Background/Foreground" class. Most of these assignments are statistical, which forwards the hypothesis that each pixel may be modeled by a random variable capable of assuming either the "Background" or "Shape" state.

Stauffer and Grimson [START_REF] Stauffer | Learning patterns of activity using real-time tracking[END_REF] proposed employing a parametric Gaussian Mixture Model (GMM) in order to depict the probability density associated with each pixel. It is also possible to use a nonparametric probability density model, such as the one in [START_REF] Isard | Bramble: A bayesian multiple-blob tracker[END_REF]. This will be the approach adopted herein.

Let I t be an image acquired at time t, and y t (u) . = {y t (i, u)} i=1,2,3 the function that provide information on pixel u of the CCD sensor according to the three colorimetric components: Red, Green and Blue. Each pixel constitutes a discrete random variable capable of assuming either one of the two following states: 1) "Background" (ω 1 ), and 2) "Foreground" (ω 2 ).

Background Model

We are proposing to discretely model the likelihood p(y t (u))|ω 1 by marginalizing the three color planes, given that:

p(y t (u)|ω 1 ) = 3 i=1 p(y t (i; u)|ω 1 ) (1)
with :

p(y t (i; u)|ω 1 ) ≈ K N j=1 q b t (i, j; u)d(b t (i; u) -j), (2)
where d represents the Kronecker function and K a constant assigned to standardize the term (dependent on i and u).

N j=1 q b t (i, j, u) = 1. q b t (i, j; u) is a weight function associated with the discrete model of the likelihood function p(y t (u)|ω 1 ). The model evolves in sync with this weight function. q b t+1 (i, j; u) is updated with each image by application of the following AR equation:

q b t+1 (i, j; u) = 1 1 + α .[q b t (i, j; u)+α.d(b t+1 (i; u)-j)] (3)

Shape Model

The shape model is defined by the likelihood of assimilation in the shape category p(y t (u)|ω 2 ). Even though statistically speaking, some pixels (e.g. those positioned above the horizon line) display a smaller probability of belonging to this shape, we are considering herein that the likelihood of assimilation to the shape category does not depend on the position of the observed pixel. It thus becomes possible to approximate the likelihood function p(y t (u)|ω 2 ) by :

p(y t (i; u)|ω 2 ) ≈ K N j=1 q f t (i, j)d(b t (i; u) -j), (4)
The weights q f t (i, j) represent a distribution discretization, which has been marginalized a priori from the color of the shape object being sought. When no information is available on the model of objects present in the foreground, an equal probability hypothesis is to be considered. In this particular case, the terms q f t (i, j) stem from a color histogram of the tracked object, as extracted during an initialization phase.

Probabilistic Shape Assimilation Map

A probabilistic shape assimilation map is generated from the likelihood ratios p(y t (i; u)|ω 2 )/p(y t (i; u)|ω 1 ). By expressing this ratio in log-likelihood form, it becomes possible to build a pseudo-image of the log-likelihood ratio in which the value associated with the pixel located at coordinates u is calculated using the following l m,t (u) function :

l m,t (u) . = log (p(y t (i; u)|ω 2 )) -log (p(y t (i; u)|ω 1 )) (5)

State Model

The state model, which identifies the trajectory characteristics to be tracked, must integrate the constraints related to vehicle kinematics. We are proposing to use a bicycle type of model and then recognize and focus on many behavioral and stability properties. The hypotheses inherent in this model are as follows:

no transfer of lateral load; the vehicle is thus compressed onto a single path, These hypotheses all imply that: vehicle acceleration at all times remains below 0.4 g (i.e. linear operating regime for tires); the steering angle, drift angle, etc. are small; and the ground surface is smooth (no suspension displacement). As an initial approach, it has been considered that vehicles are negotiating the turn at low speed; centrifugal forces are thus negligible and the tires must not develop any lateral forces, i.e.:

rolling without sliding or drift, -in order for tires not to slide laterally, the instantaneous center of rotation (ICR) of each tire must lie in the middle of the curve.

Let (x, y) be the selected coordinate system; the kinematic relations governing this model can then be written as follows:

ẋ = v. cos β ẏ = v. sin β β = v L . tan δ (6)
with β representing the vehicle direction and δ the steering angle of the front wheel within the world coordinate system (the extended Lambert II system 3 ) along the vehicle axis.

According to Ackerman's theory at low-speed behavior, the ICR lies at the intersection of the extension to 3 The absolute reference frame used here is the NTF geodesic system in extended Lambert II projection.

the rear wheel axis and the perpendicular to the front wheel plane [START_REF] Gillespie | Fundamentals of vehicle dynamics[END_REF]. The ideal front wheel deviation can be deduced from the construct illustrated in Figure 2 and its angle then written: tan δ = L R , with L denoting the vehicle wheelbase and R the curve radius of the road. By adopting the small angle hypothesis δ = L R , the system state vector is expressed as follows:

X t = (P t , β t , δ t , v t ) t (7)
with P t . = (x t , y t ) representing the vehicle position and v t the vehicle speed within the world coordinate system.

Fig. 2 The bicycle model synthesizes the displacement of a four-wheel vehicle, through the displacement of two wheels whose centers are connected to a rigid axis of length L. Ackerman's theory serves to estimate the steering angle of the front axis of a vehicle traveling at low speed.

Prediction Model

The bicycle-type kinematic model applied to each particle evolves according to the following nonlinear form:

x t+1 = x t + T.v t . cos (β t ) y t+1 = y t + T.v t . sin (β t ) β t+1 = β t + T. v L . tan δ t δ t+1 = δ t + T.b δ v t+1 = v t + T.b a (8)
with b δ ∼ N (0, σ δ ) and b a ∼ N (0, σ a ). The two terms δ and a are randomly distributed (according to a Gaussian), while the terms σ δ and σ a represent respectively the speed deviation range in the front wheel steering angle and the acceleration range, as performed by a standard vehicle during sampling period T .

Initialization

The filter initialization process consists of ascribing a hypothesis to each filter particle, such that the entire set of particles offers a stochastic representation of the density associated with the particular state, at the initial time.

The vehicle position on the road pavement is initialized based on observations generated from the first image.

The steering angle is initialized from an a priori distribution calculated as a function of characteristics describing the target curve. The heading direction is initialized using an a priori distribution calculated with respect to curve position and characteristics. The speed is calculated by the rangefinder at initial position of the vision process, the two sensors being synchronized temporally.

Determining the a priori distribution on the vehicle position relies upon a method for detecting cluster cen-ters, which associates an assimilation probability at the cluster center with each pixel labeled "Shape".

To proceed with this method we take the hypothesis α << 1 and a pixel becomes "Background" should its value remains stable for an image number k >> α.

Then, equation 2 can be simplified and the set of shape pixels F t are built by means of thresholding the likelihood function associated with the background assimilation :

F t . = u∈It {u|p(y t (u)|ω 1 ) < k.α} (9)
Typically k = 25 and α = 0.01). This probability is obtained from a non-parametric model based on Parzen estimators [START_REF] Duda | Pattern Classification, Second Edition[END_REF]:

π n t ∝ p(Z t |X t = X n t ) ≈ 1 |F t | u∈Ft ϕ(p n t , u) (10)
where ϕ(p n t , u) is a Gaussian kernel defined by:

ϕ(p n t , u) = 1 (2π).|Σ| 1/2 exp   - 1 2 (p n t -u) t Σ -1 (p n t -u)    (11)
The covariance matrix Σ of the kernel depends on the estimated vehicle projection size in the image plane. This dimension, in pixels, is a function of the vehicle position in the image; as the vehicle approaches the camera, its size grows in the image. Figure 3 shows the assimilation probability distribution at the cluster center for the background/shape extraction example.

The complexity of this method lies in O(n 2 ). In [START_REF] Goyat | Vehicle trajectories evaluation by static video sensors[END_REF],

we have proposed a stochastic approach in order to reduce the cost of running this algorithm.

Fig. 3 Illustration of the method for searching the "Shape" point cluster center: Each pixel of the set of "Shape" points displays an assimilation probability at the cluster center. The right-hand side figure is an image of the probability of belonging to the cluster center; as the pixel color becomes redder, the associated weight rises.

The position part (P n t ) of the filter is initialized as a function of Equation [START_REF] Ivanov | Fast Lighting Independant Background Subtraction[END_REF]. Figure 4 displays an example of an initialized particle set. To improve this initialization, the filter is iterated on the first image, by noising just the position and heading as well as by applying the observation model described in the next section.

The shape model (marginalized color histogram q f t ) can then be built from the pixels belonging to the estimated vehicle.

Observation Function

The proposed observation function utilizes the measurements provided by the two available sensors (i.e. color camera and 1D scanning laser rangefinder). A likelihood function must be defined for each sensor.

Vision Likelihood Function

The likelihood function proposed herein relies upon a simplified three-dimensional geometric model of the vehicle, as depicted in Figure 5. This model is composed of two nested parallelepipeds. In a general case, the model may be more complex and contain P M parallelepipeds.

Let M (R0) = {M (R0) i } i=1,..,NM represent the model's set of cube vertices ( N M = 8 × P M ), expressed within a coordinate system associated with model R 0 . This coordinate system is selected such that the 3 axes all lie in the same direction as that of the coordinate system associated with scene R w . For a given particle X n t , the likelihood (weighting) calculation is determined as the product of the shape/background likelihood ratios located inside the vehicle model projection in the image. This computation performed for each particle consumes valuable processing time, and we are proposing herein a fast likelihood calculation method based on an approximation of the 3D model projection in the image through its convex hull. Each model point is projected onto the image via the following equation:

mi ∝ C c . (Rw) T (R0) (X n t ). M(R0) i (12)
with M homogeneous coordinates associated with point M ; C c is the camera projection matrix, and (Rw) T (R0) (X n t ) the homogeneous transformation matrix between the world coordinate system and the system associated with the 3D model (cf. figure 5). This matrix, which depends on X n t , may be simply written as:

(Rw) T (R0) (X n t ) =         cos(θ n t ) -sin(θ n t ) 0 x n t sin(θ n t ) cos(θ n t ) 0 y n t 0 0 1 0 0 0 0 1         (13)
The set M (Ri) = {m i } i=1,..NM is thus built based on the projection of 3D model points within the image.

Let E(M (R0) ; X n t ) = {e i } i=1,..,N (e i = (x e i , y e i ) as coordinates of e i in the image plane) be the list of convex hull points 4 . We will now define E n t . = E(M (R0) ; X n t ) in order to streamline notations. The likelihood calculation may be performed efficiently by use of a line-by-line integral image derived from the log-likelihood ratios calculated in Equation (5) page 5. Since the concept of integral image is often used in computer vision to increase performances, it can not be used directly here because 4 Calculation of the convex hull will not be developed in this article; the calculation procedure is conducted using a classical algorithm with a complexity expressed in O(N. log N ).

the shapes are not combination of rectangles. So, we propose to extend the concept of integral image to line-byline integral image. The resulting image is build using integration along the current raw of the image (each raw is independent) and the likelihood calculation is:

l I,t ((x, y) T ) = x i=1 l m,t ((i, y) T ) ( 14
)
Points e i are categorized by pairs featuring the same y-coordinate values, such that:

E n t = { (x e
1 , y e ), (x e 2 , y e ), (x e 3 , y e + 1), (x e 4 , y e + 1), ... (x e N -1 , y e + N/2), (x e N , y e + N/2)}

Convex hull coding within the set E n t necessitates a shape discretization along the image lines. Moreover, special attention needs to be paid to coding the upper and lower extremities. On the other hand, it is not at all necessary to sort points positioned on the same line. A compliance measurement relative to a convex envelope is calculated in the integral image by application of the following relation:

a(E n t ) = N/2 j=1 [2.(l i,t (e 2j ) -l i,t (e 2j-1 )) -(x e 2j -x e 2j-1 )] ( 16
)
Figure 6 describes the principle behind the likelihood calculation method using the integral image. A line-byline scanning is performed as part of this method.

The vision weight associated with each particle is directly correlated with a(E n t , I I ) by means of the following expression: This projection is approximated by its convex hull (shown in red on the lower image). The likelihood calculation proceeds in a line-by-line integral image of the log-likelihood ratio.

π n v,t ∝ p(Z t |X t = X n t ) . = max(0, a(E n t )) (17)

Telemetric Likelihood Function

The telemetric sensor provides, on a horizontal plane, the distance from the first obstacle with a resolution of 1 degree. The intersection of the laser beam with a vehicle yields a set of points; for each particle, a model is now available for calculating the laser echoes generated from the intersection of the beam with the simplified 3D model of the projected vehicle in the world coordinate system.

The likelihood associated with the telemetric observation can then be calculated from a modified Hausdorff distance (denoted d h ) between the actual and simulated echoes (cf figure 7) using the following expression:

π n l,t ∝ p(Z t |X t = X n t ) . = exp (-λ t d h ) ( 18
)
where λ t is an adjustment parameter (typical value is 20).

Multi-Source Sampling

When the observation function is composed of terms stemming from multiple sources, it becomes necessary to lay out a strategy for combining these terms. A classical approach consists of assembling a weighting function composed of a combination of observations generated from each source. The choice of merge operator often takes place empirically. We propose an alternative to this approach, which calls for performing the combination step as part of the particle filter re-sampling step.

Upon observation, the filter may be represented by a set of N particles with an associated weight vector:

{X (i)
t , π i t } i=1,...,N . The weight vector π i t , of size M in correspondence with the number of sources, is composed of the individual particle weight estimated by each source.

To facilitate comprehension, the notations contained in the remainder of this section will omit the time index t.

Principle

Multi-source sampling entails generating a new particle set, by implementing a three-step approach:

1. M particles are sorted (one for each source) according to an importance sampling strategy associated with each source (importance sampling). The output of this step consists of a set of M candidate particles along with their respective weight vector {X (i) , π (i) } i=1,...,M : 2. A confidence vector, of dimension M , is built using likelihood ratios estimated for each candidate particle (this calculation will be detailed in the following discussion).

3. The chosen particle is derived from a selection performed among candidate particles by applying an importance sampling strategy on the confidence vector. These three steps are then repeated a total of N times in order to obtain the complete set.

Confidence Vector

In the following section, we will describe in detail the second step of this multi-source sampling procedure, which is aimed at building a confidence vector associated with the set of candidate particles. The underlying principle consists of calculating the product of likelihood ratios between weights of the same source, for each pair of candidate particles.

As an example, in the three-sensor case, three candidate particles are drawn; for each particle, a likelihood ratio product is calculated, which for the first candidate particle yields:

r 1 . = π 1 1 π 2 1 . π 1 1 π 3 1 ( 19
)
where π i j represents the jth component of vector π i . In the case of a blind source (i.e. the values returned by the observation function associated with this source are constant), those likelihood ratios in which the source is present equal one and do not exert any influence on the calculation of terms r i . In a general case, it is preferable to introduce log ratios, which makes it possible to compute a vector l r , i.e. the log of r, featuring coefficients r i , as indicated in the following expression:

l r . = M                     1 (1×M )   lπ 1 - 1 M M i=1 l πi    1 (1×M )   lπ 2 - 1 M M i=1 l πi    ... 1 (1×M )   lπ M - 1 M M i=1 l πi                        (20)
where l πi is the log of vector π i and 1 (1×M ) is a matrix composed of a single line and M columns all containing one.

By setting C π . = 1 M M i=1 l πi , l r , the following may be written:

l r = M         1 (1×M ) (l π1 -C π ) 1 (1×M ) (l π2 -C π ) ... 1 (1×M ) (l π M -C π )         (21)
Confidence vector c is obtained by normalizing r through a coefficient C c that results in the sum of its elements equaling one.

c . = C c .exp (l r ) (22)
Algorithm 2 Multi-source sampling Input : Particle set and associated weight vector {X (i) , π i } i=1,...,N , M sources for n = 1 to N do -Choose M candidate particles on the basis of {X (i) , π (i) } i=1,...,N and build {X * (j) , π * (j) } j=1,...,M where X * (j) is derived from an importance sampling drawn on source j weights; -Calculate vector l r based on Equation 21, and then calculate confidence vector c . = C c .exp (l r ) -Select the designated particle X e(n) from among the candidate particles by proceeding with an importance sampling drawing. end for Output : Particle set {X e(i) } i=1,...,N composed of the selected particles. Fig. 8 Illustration of multi-source sampling method operations for 2 different scenarios (one in each column). In the left column, from top to bottom: the first two curves represent the source response (observation). The third shows the result of a weight-based importance sampling, calculated as the product of weights from the two sources. In the fourth curve, this weight-based importance sampling result is generated by summing weights from the two sources. The last curve then displays the result of our proposed multi-source sampling method.

Experimental Validation

This section will present the experimental campaigns carried out in order to validate the trajectory tracking method.

The video system developed as part of the SARI project has made it possible to continuously record three color 640 x 480 video streams (cf. Fig. 9), at a frequency of 30 fps, over a several-day period along with data recorded by a scanning telemetry sensor. The calibration between sensors enables expressing all measurements within a common absolute coordinate system tied to the GPS sensor coordinate system, which provides the actual situation in the field. puted trajectories were then compared with trajectories estimated by a vision-only approach, a rangefinder-only approach and an approach merging the two sources. The Table 1 presents the average error and related standard deviation vs. vehicle speed. It can be noticed that the level of accuracy varies only slightly with respect to the tested speed range. Moreover, the two-sensor merge serves to significantly improve this accuracy, which lies on the order of 15 cm with a standard deviation of approximately 10 cm.

Note that these results come from the right camera, but are extensive to the others cameras. The calibration phase is same for all the system and the reference frame is absolute. Figure 11 shows the three sensors merging with the ground truth.

Conclusion -Prospects

This article has set forth a method for estimating vehicle trajectories by use of a sensor composed of a color camera and a 1D scanning laser rangefinder. Based on a time monitoring approach formalized by a particle filter, the algorithm has output, at every point in time, an estima- The experiments conducted have enabled quantifying the precision associated with the proposed approach, thanks to generating actual field values (using a kinematic GPS accurate to within 1 cm). The accuracy obtained lies on the order of 20 cm. These tests have also served to demonstrate the contribution offered by merging vision with telemetry.

Upcoming research will primarily be focused on improving the background subtraction method, which has the potential to introduce bias due to the presence of a vehicle shadow projected onto the pavement [START_REF] Salvador | Cast shadow segmentation using invariant color features[END_REF].

The method discussed in this article is currently operating within the scope of an ANR-PREDIT project. For the purpose of covering the entire curve, the system is composed of three color cameras with very little overlap plus a scanning laser rangefinder, which has successfully analyzed the observations recorded under actual traffic conditions over several-day periods. Fig. 10 Tracking example drawn from the testing campaign performed as part of this research (for the right camera). The model corresponding to vehicle localization is re-projected onto the current image. For the left-hand column, the method is confined solely to data stemming from the vision sensor. For the middle column, the method uses only the telemetric data, with the green crosses corresponding to simulated laser firings re-projected onto the image and red crosses indicating actual re-projected laser data. The right-hand column depicts method application by combining the two sensors.

  propose a probabilistic model to track multiple vehicles with spatio-temporal Markov random fields. Recently, stochastic methods [13,22] have been presented to handle with realtime multi object tracking. In [25] Yu and al. propose a Monte Carlo Markov Chain method to estimate, from a video, global trajectories of the vehicles. In this paper, we propose a probabilistic framework to solve the problem of online vehicle trajectory estimation. The trajectory is modeled by a random state vector and the distribution associated to this random vector is approximated in a sequential scheme with a particle filter. Since the trajectory is highly driven by the kinematic model of the vehicle, we propose to inject this model within the dynamics associated to the filter. Moreover, an original data fusion sampling algorithm is proposed to handle with several observation functions.3 Principle of the MethodThis section will set forth the principle behind the trajectory tracking method.The core of this proposed method (see Fig.1) consists of a recursive filter that has been formalized stochastically (using a particle filter). The vehicle state is repre-sented by its corresponding bicycle model. The prediction function adopts the hypothesis of a constant driving angle acceleration and speed. A 3D geometric model of the vehicle, projected onto the image, is then compared with image data described by a probabilistic shape assimilation map. The observation function also comprises a likelihood function that reflects the consistency of telemetric data with respect to the hypothesis. The filter is able to produce, at each iteration, an estimation of the vehicle state (position, heading direction, speed, steering angle).

Fig. 1

 1 Fig. 1 Block diagram of the proposed method

-

  no longitudinal transfer, -no roll or pitch motion, -tires in a linear configuration, -constant forward speed V, -no aerodynamic effects, -position control, and no effect from suspension and chassis flexibility.

Fig. 4

 4 Fig. 4 Example of particle distribution during the initialization phase (shown in green).

Fig. 5

 5 Fig.5Example of the simple three-dimensional geometric model used for a vehicle: it is composed of two cubes. The coordinate system associated with the cube and the other system associated with the scene are related according to pure translation. The plane (Oxy) of the world coordinate system and component axes are merged with the GPS coordinate system (Lambert II).

Fig. 6

 6 Fig. 6 Illustration of the vision likelihood computation. The 3D model of the vehicle (shown in green) is re-projected onto the image generated from the background-shape extraction.This projection is approximated by its convex hull (shown in red on the lower image). The likelihood calculation proceeds in a line-by-line integral image of the log-likelihood ratio.

Fig. 7

 7 Fig. 7 Illustration of the distance computation between simulated echoes on the 3D model of the cube (shown in green) and the actual echoes (red).

Figure 8

 8 Figure 8 demonstrates how the multi-source sampling method works for two distinct scenarios, by comparing method behavior with that displayed when sampling from both a weight function composed of the weight product from each source and another weight function composed of the sum of weights from each source. In

Fig. 9

 9 Fig.9The video system is composed of three color cameras, capable of being easily deployed onsite.

  error was quantified as the average distance between each estimated vehicle position and the straight line passing through the two closest GPS points. For each test, at least five vehicle passes were carried out, which enabled deriving a very rough statistic on the recorded measurements. For the tests actually conducted, the vehicle has been tracked in a curve over a distance of approximately 100 m. The vehicle model used for these purposes comprises a single cube. Finally, according to the empirical results of the

Figure 10

 10 Figure 10 shows a tracking example from the test campaign carried out as part of this research. The model corresponding to vehicle localization has been re-projected onto the current image. For the left-hand column, the method makes exclusive use of data stemming from the vision sensor. The middle column reflects results based solely on telemetric data. The green crosses correspond

Fig. 11

 11 Fig. 11 representation of the three trajectories calculated by each camera and the ground truth associated.

Table 1

 1 Precision for the right camera (error and standard deviation of error between the estimation and actual field value output by a kinematic GPS) of the proposed method, for several travel speeds. It may be remarked that merging the two sensors serves to considerably improve estimation precision. Actual field values have been furnished by a kinematic GPS type of sensor.

	Speed	Vision	Rangefinder Sensor merge
	km/hr	ave/std	ave/std	ave/std
	40	0.25/0.18	0.65/0.54	0.17/0.10
	60	0.19/0.16	0.72/0.67	0.09/0.06
	80	0.18/0.15	0.33/0.22	0.14/0.10

Table 2

 2

	particules number	50	100 150 200 300
	average	0.61 0.22 0.13 0.14 0.14
	standard deviation 0.61 0.24 0.14 0.14 0.14

, for every tracking, we decide to use 150 particules.

Table 2

 2 Precision for the right camera (error and standard deviation of error between the sensor merge estimation and GPS) of the proposed method, for several particules number.

Acknowledgements This research work was conducted as part of the PREDIT project entitled "SARI" (French acronym for Automated Road Monitoring to Inform drivers and facility managers). This project is intended to help lead to significant reductions in the number of accidents due to loss of contact with the roadway or loss of vehicle control, by means of better informing drivers of the upcoming difficulties facing them on the road (http://www.sari.prd.fr).


Santiago Arroyave-Tobon

Kalenia Marquez-Florez

Paul Heymann

Jean-Marc Linares

Generative design of joint contact surfaces inspired by biological morphogenesis

Keywords: Biologically inspired design, Geometry, Stress

Recent knowledge about biological joint morphogenesis opens new perspectives in mechanics to automate joint contact surface design. The present work evaluates the feasibility of a generative design method inspired by the joint morphogenesis process to develop contact surfaces. A finite element model of joint morphogenesis reported in literature was implemented. This morphogenesis process was adapted and implemented for mechanical applications. The results show that the bioinspired joint shaping process adds matter in the zones next to the contact zone decreasing the contact pressure (up to 57%). The results demonstrate the feasibility of implementing biological growing rules in generative design.

Introduction

Improving the engineering design process is a subject of interest for science and industry. Hence, several theories and methodologies have been developed [START_REF] Tomiyama | Design methodologies: Industrial and educational applications[END_REF], most of them as decision support tools. Computer-automated tools have been proposed to help engineers in the different phases of the design process. For the detailed design of mechanical parts, several automated approaches have been proposed; they can be classified into topology [START_REF] Rozvany | A critical review of established methods of structural topology optimization[END_REF], shape [START_REF] Sysaykeo | Clearance and design optimization of bio-inspired bearings under off-center load[END_REF] and size optimization [START_REF] Krish | A practical generative design method[END_REF]. These methodologies are specialized in structural design and are less applicable to mechanical joints. Currently in mechanical design, geometry of contact surfaces of a joint is imposed by the designer, usually resulting from technological background. The literature about mechanical joint design has been mainly focused on the optimisation of roller bearings contact profiles [START_REF] Mermoz | New methodology to define roller geometry on power bearings[END_REF] and gears lead contact profiles [START_REF] Bergseth | Logarithmical crowning for spur gears[END_REF] under normal loads. Although the proposed solutions have contributed to improving the load transmission performance of mechanical joints, they cannot be easily generalized to other loading conditions. Additionally, the manufacturing costs have limited the development of new mechanical joints made up of free-form contact surfaces, even if it allows to increase the durability of the systems [START_REF] Sysaykeo | Clearance and design optimization of bio-inspired bearings under off-center load[END_REF]. However, current sustainability trends suggest reviewing our design practices [START_REF] Roucoules | Coevolution of digitalisation, organisations and Product Development Cycle[END_REF].

For its part, vertebrates have already been confronted with the mechanical problems mentioned above: optimization of load transmission, durability, energy efficiency and frugality. With many biological solutions offered to us by nature, bioinspiration is an emerging field in engineering design [START_REF] Byrne | Biologicalisation: Biological transformation in manufacturing[END_REF]. Some structural design methodologies have been inspired by natural objects or processes. Recent knowledge developed by biologists on the natural process of joint shaping (joint morphogenesis) [START_REF] Carter | The role of mechanical loading histories in the development of diarthrodial joints[END_REF] opens new perspectives. The knowledge behind the process of joint morphogenesis could be used in mechanical design to automatically shape functional contact surfaces. However, the application of the joint morphogenesis rules cannot directly be used in a mechanical context because of the differences in the operating conditions (material properties, contact conditions, and self-repairing capacities). Based on this, the research question investigated in this work is: how to adapt the joint morphogenesis process to design of contact surfaces? In order to answer this question, the present work evaluates the feasibility of a generative design method inspired by the joint morphogenesis process to design functional contact surfaces.

Synovial joint morphogenesis

Joint morphogenesis theory

In synovial joint morphogenesis, the shaping of contact surfaces is due to a smart growth of the cartilage tissue. Studies have shown that cartilage growth might be regulated by mechanical stimuli. According to Carter's theory [START_REF] Carter | The role of mechanical loading histories in the development of diarthrodial joints[END_REF], the behaviour of the chondrocytes (cartilage tissue cells) is highly dependent on the shear octahedral stress (or von Mises stress) (σVM) and the hydrostatic stress (σhyd). Von Mises stress promotes the ossification of the cartilage through the stimulation of the chondrocyte hypertrophy, whereas the compressive hydrostatic stress inhibits the ossification by keeping the cells in a proliferative state [START_REF] Carter | The role of mechanical loading histories in the development of diarthrodial joints[END_REF]. Experimental studies have shown that static compression stress inhibits cartilage growth, while dynamic compressive stress, conversely, promotes it [START_REF] Guilak | Mechanical and biochemical changes in the superficial zone of articular cartilage in canine experimental osteoarthritis[END_REF]. In addition to mechanical stimuli, cartilage growth is also affected by some biochemical factors.

Morphogenesis theory for mechanical applications

Synovial joint morphogenesis is a frugal shaping process in which nature disposes of matter in a very selective way. In contrast, the design of mechanical joints is usually driven (or even constrained) by the manufacturing costs. Hence, for ages, mechanical parts have been made up of simple geometries (planes and cylinders). As described by [START_REF] Carter | The role of mechanical loading histories in the development of diarthrodial joints[END_REF], growth rules of the biological joint morphogenesis process are physically analogue with the multiaxial fatigue failure criteria; both are affected by cyclic von Mises and hydrostatic stresses. Hence, there exists some knowledge behind this biological process that has not been used yet in mechanical design to automatically shape functional contact surfaces.

Computational models of joint morphogenesis

During joint morphogenesis, local tissue growth, directly related to chondrocytes proliferation, is defined as locally isotropic. That means that the growth tensor G can be assimilated to a scaled identity tensor I, as G = γ •I; where γ is the growth magnitude. This means that the shaping of biological joints is dependent on the magnitude of γ so that some areas grow faster than others, locally and iteratively adapting the geometry. Representing the mechanical and biological growth factors, the growth tensor G can be divided into two parts:

G = (a1•Cb + a2•Cm) I (1)
with Cb and Cm biological and mechanical contributions, respectively. The constants a1 and a2 correspond to the strength of the mechanical or biological stimuli regarding the overall growth. Cb is related to the chondrocytes' density [START_REF] Heegaard | Mechanically modulated cartilage growth may regulate joint surface morphogenesis[END_REF]. Cm is defined by σVM and σhyd. From this base, several research works have proposed computational models about the joint morphogenesis process. The first one was developed by [START_REF] Heegaard | Mechanically modulated cartilage growth may regulate joint surface morphogenesis[END_REF] who considered that only the positive σhyd promoted cartilage growth. Later, [START_REF] Giorgi | Mechanobiological simulations of prenatal joint morphogenesis[END_REF] proposes a model in which the cartilage growth was regulated by the negative σh, and they were the first to include the synovial capsule (as a spherical uncompressible domain surrounding the cartilaginous bone rudiments). This model was enhanced by [START_REF] Carrera-Pinzón | Computational model of a synovial joint morphogenesis[END_REF] including the effect of the von Mises stress (σVM) as a growth inhibitor in the regulation of the articular cartilage growth. All models found in the literature follow an iterative loop, in which the tissue grows on each step until a stop criterion is reached.

Implementation of a joint morphogenesis computational model

The starting point of this work was to replicate a computational model of a joint morphogenesis existing in the literature. The implemented model was the one presented in [START_REF] Carrera-Pinzón | Computational model of a synovial joint morphogenesis[END_REF]. From that work, we simulated the case of the morphogenesis of an interphalangeal joint developing under palsy conditions. The initial geometry, material properties and boundary and loading conditions were reproduced as closely as possible. A mechanical simulation was performed to calculate the stress state of the bone rudiments. Then, σVM and σh stress fields resulting from this simulation were extracted. These values were used to calculate the field of Cm (see Eq. ( 1)). Cb was calculated following the same procedure of [START_REF] Heegaard | Mechanically modulated cartilage growth may regulate joint surface morphogenesis[END_REF]. A ratio of a1/a2 =0.1 was considered to weigh Cb and Cm. The growth magnitude γ was considered as a temperature in a thermal expansion simulation to imitate the tissue growth (as it has been done in [START_REF] Giorgi | Mechanobiological simulations of prenatal joint morphogenesis[END_REF]). Thermal and mechanical problems were solved with finite element methods and implemented in the open-source (finite element analyses) FEA software tools Salome-Meca and Code-Aster. At the end of each iteration, the computed deformed geometry was extracted in a stress-free condition so the loading conditions could be applied again for the next step of growth.

The results of the reference study and those of our simulation can be seen in Fig. 1(a) and (b), respectively. The simulation lasted 7 steps as in the reference study. Contact surfaces from both models converged to the same shape. In both simulations, the bottom surface of the distal phalanx evolved into a concave shape, whereas the upper surface of the proximal bone converged more to a convex one. Although the obtained geometry is like the biological reference, the load transmission performance was not evaluated. Then, this implementation of joint morphogenesis was adapted to a mechanical context to assess its applicability for the generative design of contact surfaces of mechanical parts.

Bioinspired design of contact surfaces of mechanical joints

In order to adapt the concepts of synovial joint morphogenesis to the design of contact surfaces of mechanical joints, two contexts were defined: one related to the application case and another related to the morphogenesis process (see Fig. 2). The application context (Fig. 2.a) is the one in which the materials (Young's modulus E and Poisson's ratio υ) and the load (F) imposed by the design case are considered. In the application context, the contact pressure and stress distributions are calculated to evaluate the load transmission performance during the simulation. Load transmission between parts is considered by means of direct dry contact. The morphogenesis context (Fig. 2.b) is an abstract context in which the shaping of the "morphogenesis pieces" takes place. This context simulates the cartilaginous state of the bones during the first stages of the morphogenesis, with a load Fm. In the morphogenesis context, load transmission between parts is done through a capsule (equivalent to a synovial capsule) represented in light blue in Fig. 2b. The contact zone susceptible to growth is called hereafter malleable zone and represents the epiphysis in a biological context (region represented in yellow in Fig. 2b). The proposed bioinspired algorithm (presented in Fig. 3.) requires as inputs an initial geometry of the mechanical parts, the loading conditions and the material properties of the application case. The mechanical parts adapt their shape according to the mechanical stimuli following bioinspired growing rules. The obtained geometry is evaluated in the application context at each iteration by means of the peak contact pressure, hydrostatic stress and von Mises' stress. Hereafter, each step of the algorithm is explained. Create initial mesh. The first stage of the biological joint morphogenesis process is the cavitation, where the bone rudiments and the synovial capsule develop. In the bioinspired design algorithm, this stage corresponds to the creation of the initial geometry of the parts in contact. Using this geometry, the capsule is then created with Boolean operations. The capsule is considered as an incompressible domain surrounding the joint [START_REF] Giorgi | Mechanobiological simulations of prenatal joint morphogenesis[END_REF][START_REF] Carrera-Pinzón | Computational model of a synovial joint morphogenesis[END_REF]. Afterwards, the contact parts and the capsule are meshed. The Young's modulus of the capsule was defined in function on that of the mechanical parts: Ec=E/1000, as in other studies [START_REF] Carrera-Pinzón | Computational model of a synovial joint morphogenesis[END_REF]. Structural FEA (morphogenesis context). Cartilage growth is mainly regulated by the mechanical stimuli generated by the interaction of the bone rudiments. In the bioinspired design algorithm, this mechanical interaction is simulated by means of a structural FEA. This FEA is performed with the meshed contact elements and the capsule. The loading and boundary conditions applied to the soft parts are the same as within the application case defined by the user (Fm = F). The transmission of the load between the contact parts is done through the capsule.

Get stresses σhyd and σVM. During joint morphogenesis, the hydrostatic and the von Mises stresses are the ones that influence the most on the cartilage growth [START_REF] Carter | The role of mechanical loading histories in the development of diarthrodial joints[END_REF]. Therefore, it is necessary to obtain the morphogenesis parts' stress distribution in the bioinspired methodology. From the results of the structural FEA, the σhyd and σVM stress fields are calculated.

Set growth potential β. The processes related to joint morphogenesis occur mostly at the epiphyses of the bones (distal and proximal portions of the long bones which joints with the adjacent bones). Nature sets the epiphyses' size based mainly on bio-molecular factors, which inhibit or promote certain cellular processes, depending also on the stage and type of cell. Within the bioinspired algorithm, this aspect is reproduced by the parameter named as growth potential β. By means of this parameter, the region where the morphogenesis parts are allowed to grow is defined. This region is where more matter can be added (malleable zone in Fig. 2). The growth potential can take value from 0 to 1, being 0 at the non-malleable zone and 1 at the malleable zone:

𝛽 = { 1 if abs(𝑦) ≤ 𝜆 • 𝑎 0 if abs(𝑦) > 𝜆 • 𝑎 (2)
The size of the malleable zone is defined λ times proportional to a geometrical characteristic of the initial geometry: the halfwidth a (Figure 2). The depth along the symmetrical axis of the parts is named y. The transition of the growth potential between the malleable zone and the non-malleable zone was smoothed.

Set growth factor Gf. In nature, cellular processes are regulated by growth factors, which are substances that act as signalling between cells. These signals (molecules) are produced by the cells in order to communicate to their neighbour peers the stimuli they are sensing. When a cell receives these signals, a chain of events is triggered, including the amplifying of the signal and cell processes, among which, cell proliferation. Through the biological joint morphogenesis, cartilage cells express growth factors to signal the mechanical (hydrostatic and von Mises stresses) and bio-molecular stimuli. This process regulates how cartilage cells should proliferate. This aspect was included in the algorithm by calculating a term named growth factor, Gf. This factor depends on the hydrostatic stress, σhyd; the von Mises stress, σVM; the stress limit, σlim; and the growth potential, β. The stress limit, σlim, which inhibits cartilage growth, was used as a switch in the algorithm. A limit value of σlim = 0.5•σmax was defined, where σmax corresponds to the maximum von Mises stress within the malleable zone. The growth factor, Gf, was calculated as follows:

𝐺𝑓 = { 𝛽 • 𝜎 ℎ𝑦𝑑 if 𝜎 ℎ𝑦𝑑 < 0 & 𝜎 𝑉𝑀 < 𝜎 𝑙𝑖𝑚 0 in other cases (3)
Set growth coefficient αi. Nature regulates tissue growth either by cell hypertrophy, cell proliferation or a combination of both. During joint morphogenesis, proliferation is the predominant mechanism. Cell proliferation rate is the factor that influences the tissue proportions (as well as the orientation of the cell division). The proliferation rate depends on the type of cell and its current stage. This effect within the bioinspired algorithm is represented by a growth coefficient, αi. This term was defined as evolving linearly from a high growth coefficient (highly proliferative), αini, to a null growth coefficient (non-proliferative), with each step:

𝛼 𝑖 = 𝛼 𝑖𝑛𝑖 ( 1-𝑖 𝑁-1 + 1) (4)
where i is the current growth step, and N is the total number of growing steps. Isotropic growth FEA (addition of new matter). The growth can be modelled through the isotropic expansion of elements [START_REF] Giorgi | Mechanobiological simulations of prenatal joint morphogenesis[END_REF]. This effect can be simulated as a thermal expansion using FEA; being the thermal coefficient equivalent to the growth coefficient, αi, and the temperature to the growth factor, Gf. Due to the thermal expansion, the morphogenesis parts get deformed. The volume variation is considered as addition of new matter, and the deformed state is set as stress-free (stresses are reset). With the obtained geometry a new capsule is generated with Boolean operations.

Contact performance assessment FEA (application context). With the new geometry, a contact FEA is carried out in order to assess the evolution of contact performance. Load transmission between parts is considered by means of direct dry contact. The peak contact pressure is evaluated as a performance variable.

Convergence test. At this stage, it is evaluated if the number of growth steps is reached.

Results

The proposed bioinspired methodology for the generative design of contact surfaces was tested in a design case of a unilateral contact joint. A case of two steel (E=200 GPa and υ=0.3) cylindrical parts in contact was considered. The application case was treated as an axisymmetric model. A static force F of magnitude 160 N was applied perpendicular to the upper surface of the top part, as represented in Fig. 2.a. The top part was only allowed to move vertically. The nodes of the lower surface of the bottom part were only allowed to move laterally, while all displacements of the centre node of the same surface were constrained. The diameter of the parts (2a) was defined as the same for both and equal to 5 mm. The height of both parts was defined high enough to prevent the boundary conditions from affecting the stress distribution on the contact surfaces. The height of the malleable zone was defined by λ=1. The initial geometry of the contact surface of the top part was defined by an ellipsoidal face (major radius: 2.50 mm, minor radius: 0.75 mm). The initial geometry of the contact surface of the bottom part was defined by a horizontal planar face. The peak contact pressure was chosen as a monitoring variable to evaluate the load transmission performance during the bioinspired shaping process. A convergence analysis allowed to define an optimal mesh size of 9E-3 mm. The results obtained in the first iteration were successfully compared with the theoretical Hertzian solution. The bioinspired algorithm was implemented using the programming language Python and the FEA software Salome-Meca and Code-Aster. Each iteration (growth step) lasted around 40 min on a standard laptop (64-bit i7 processor @ 2.7 GHz, 16 GB RAM, Ubuntu 18.04 operating system). The number of iterations was set up manually here to 20. For that number of iterations, a stabilization of the monitored variable (peak contact pressure) was ensured.

The results of the simulation are presented in Fig. 4. The evolution of the peak contact pressure is plotted in function of the growth step. In Fig. 4, the shaped contact geometry is presented in blue colour after (a) 4 growth steps, (b) 8 growth steps, (c) 11 growth steps and (d) 20 growth steps. In those figures, the initial geometry is represented by a grey shaded region to show clearer the new matter added by the algorithm. The adjustment of the contact geometry improves the load distribution which results in a 57% reduction of the maximum contact pressure compared to its initial value. It is important to remark that the bioinspired joint shaping algorithm adds matter around the contact zone, increasing the contact area and gradually decreasing the contact pressure. At the end of the bioinspired morphogenesis process, the geometry is more conformal than that at the beginning. The von Mises and hydrostatic stress distributions are presented in Fig. 5 (a) and (b) respectively, for the initial and the final geometry. A reduction of about 80.4% and 81.4% of the maximum von Mises and hydrostatic stresses, respectively, was achieved. These results demonstrate that the shaping of the geometry, thanks to the bioinspired growing rules, allows a better stress distribution.

Conclusion

Knowledge about biological joint morphogenesis opens new perspectives in mechanics to automate joint contact surface design. A mechanical transcription of cartilage growth rules was proposed and implemented to automatically shape contact surfaces. This work is a proof-of-concept demonstration about a generative design method inspired by biological morphogenesis. It demonstrates that the automated design of contact surfaces is feasible by mimicking cartilage growth and ossification.

The results obtained in this work are encouraging and can be of interest in mechanical design. A better distribution of the contact stresses (lower peak contact pressure) in mechanical joints means lower potential wear and therefore more durable systems. Another potential application of this work is custom prothesis contact surfaces shaping, in which the real joint forces of the subject are considered during the design process.

It must be noted that the algorithm, as presented, is sensitive to the initial geometry, which now is user defined. For future work, the starting point of the algorithm might be defined as a function of the design constraints (loads and material properties). Additionally, the algorithm simplifies some biological processes such as the ossification. Nevertheless, future work will enable local changes in material properties by mimicking local ossification (as in bone development). This would allow also to explore alternative ways to define the convergence of the algorithm.

The work proposed in this article focuses on the genesis of joint surfaces. Prospects of development of this method would be to consider complementary biological processes, such as bone remodelling and self-repairing properties. This work is a step toward a computed-automated methodology for shaping functional mechanical parts.

Figure 1 .

 1 Figure 1. Results of morphogenesis of an interphalangeal joint developing under palsy conditions from (a) a biology study from literature [13] and (b) this study.

Figure 2 .

 2 Figure 2. Equivalences between: (a) the application and (b) the morphogenesis contexts for the bioinspired design of contact surfaces.

Figure 3 .

 3 Figure 3. Flowchart of the proposed methodology for the generative design of contact surfaces inspired by biological joint morphogenesis.

Figure 4 .

 4 Figure 4. Evolution of the peak contact pressure and evolution of the generative geometry after: (a) initial configuration, (b) 4 growth steps, (c) 8 growth steps, (d) 11 growth steps and (e) 20 growth steps. For (a), (b), (c), (d) and (e) the aspect ratio vertical/horizontal was set equal to 4.

Figure 5 .

 5 Figure 5. (a) von Mises stress and (b) hydrostatic stress distribution for the initial geometry and the final geometry for the case of a unilateral contact problem. Aspect ratio vertical/horizontal was set equal to 4.

Acknowledgements

This research was supported by French Research National Agency (ANR) Grant No. ANR-20-CE10-0008.