Mikaël Le Pendu

Aljosa Smolic

Christine Guillemot

HIGH DYNAMIC RANGE LIGHT FIELDS VIA WEIGHTED LOW RANK APPROXIMATION

Keywords: High Dynamic Range (HDR), Light Fields, Low Rank Matrix Completion, Weighted Low Rank Approximation

In this paper, we propose a method for capturing High Dynamic Range (HDR) light fields with dense viewpoint sampling. Analogously to the traditional HDR acquisition process, several light fields are captured at varying exposures with a plenoptic camera. The RAW data is de-multiplexed to retrieve all light field viewpoints for each exposure and perform a soft detection of saturated pixels. Considering a matrix which concatenates all the vectorized views, we formulate the problem of recovering saturated areas as a Weighted Low Rank Approximation (WLRA) where the weights are defined from the soft saturation detection. We show that our algorithm successfully recovers the parallax in the over-exposed areas while the Truncated Nuclear Norm (TNN) minimization, traditionally used for single view HDR imaging, does not generalize to light fields. Advantages of our weighted approach as well as the simultaneous processing of all the viewpoints are also demonstrated in our experiments.

INTRODUCTION

The emerging technologies in light field capture and HDR imaging allow increased freedom for content creation and manipulation by producing a richer description of a scene than traditional images. For instance, light fields make it possible to change the focus after taking a picture and HDR images can be freely re-exposed without unveiling either noise in dark areas or saturated pixels in the bright parts of the image. A light field is typically represented as an array of viewpoints. Many capture systems have been designed, including camera arrays, single cameras mounted on moving gantries and plenoptic cameras with a micro-lens array placed in front of the sensor. Regarding HDR images, their acquisition generally involves taking several pictures of the scene at different exposures. A large panel of methods have been developed for merging the captured images into an HDR one, including patch-based methods [START_REF] Hu | Hdr deghosting: How to deal with saturation ?[END_REF][START_REF] Sen | Robust Patch-Based HDR Reconstruction of Dynamic Scenes[END_REF][START_REF] Gallo | Locally non-rigid registration for mobile hdr photography[END_REF], Low Rank Matrix Completion (LRMC) [START_REF] Oh | Robust high dynamic range imaging by rank minimization[END_REF][START_REF] Lee | Ghost-free high dynamic range imaging via rank minimization[END_REF][START_REF] Lee | High dynamic range imaging via truncated nuclear norm minimization of low-rank matrix[END_REF][START_REF] Lee | Computationally efficient truncated nuclear norm minimization for high dynamic range imaging[END_REF], and more recently deep learning [START_REF] Kalantari | Deep high dynamic range imaging of dynamic scenes[END_REF]. A comprehensive review of the subject is provided in [START_REF] Gallo | Chapter 3 -stack-based algorithms for HDR capture and reconstruction[END_REF].

While both HDR and light field imaging have received a lot of attention, only a few attempts have been made at combining these technologies. Manakov et al. [START_REF] Manakov | A reconfigurable camera add-on for high dynamic range, multispectral, polarization, and light-field imaging[END_REF] have developed a camera add-on capable of transforming a standard camera into either a HDR or a light field camera. In [START_REF] Georgiev | High dynamic range image capture with plenoptic 2.0 camera[END_REF][START_REF] Georgiev | Using focused plenoptic cameras for rich image capture[END_REF], a focused plenoptic camera with a micro-lens aperture pattern is designed to capture an HDR image in a single shot. However, the authors of [START_REF] Manakov | A reconfigurable camera add-on for high dynamic range, multispectral, polarization, and light-field imaging[END_REF][START_REF] Georgiev | High dynamic range image capture with plenoptic 2.0 camera[END_REF][START_REF] Georgiev | Using focused plenoptic cameras for rich image capture[END_REF] do not discuss the simultaneous HDR and light field capture. Based on a camera array setup, the method in [START_REF] Schedl | Coded exposure hdr light-field video recording[END_REF] is able to capture HDR light field videos by varying the exposure time both sequentially (between successive frames) and spatially (over the different cameras). Note that the sparse view sampling resulting from the large baseline of camera array devices causes strong angular aliasing when generating refocused images from the light field (i.e. sharp structures in the out of focus areas). This issue is observed in the results of [START_REF] Schedl | Coded exposure hdr light-field video recording[END_REF]. Similarly, related methods in multiview HDR capture [START_REF] Orozco | Chapter 3 -HDR multiview image sequence generation: Toward 3d HDR video[END_REF] are not suitable for light field rendering that require dense viewpoint sampling. Finally, Li et al. [START_REF] Li | High dynamic range and all-focus image from light field[END_REF] capture multiple exposures with a plenoptic camera and merge them into a HDR light field by applying directly the method of Debevec and Malik [START_REF] Debevec | Recovering high dynamic range radiance maps from photographs[END_REF] to 4D light fields instead of 2D images. This can be seen as treating each viewpoint independently. However, the method in [START_REF] Debevec | Recovering high dynamic range radiance maps from photographs[END_REF] is known for requiring many different exposures with overlapping well-exposed areas, which complicates the acquisition task. Furthermore, the large amount of data to process increases the computational load, particularly for light fields with densely sampled viewpoints.

In this paper, we propose a new HDR light field acquisition method based on multiple plenoptic captures (typically two or three) with varying exposure as depicted in Fig. 1. By taking the highest exposure as a reference, generating the HDR light field amounts to recovering the saturated areas using the lower exposure captures. For that purpose, soft detection of the saturated pixels is performed from the RAW data. It accounts for the fact that in practice, sensor response is not linear close to saturation. The whole set of views (all sub-aperture images at all exposures) is then arranged in a matrix where each column is a vectorized view. The matrix is completed with Weighted Low Rank Approximation (WLRA) naturally exploiting the redundancies between views. In our approach, the WLRA problem is solved by extending the matrix completion algorithm of [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] to nonbinary weights, thus accounting for the soft saturation model. We show that our non-binary approach better handles the transition between the saturated and non-saturated areas. Furthermore, unlike the the Truncated Nuclear Norm minimization previously used in several HDR imaging methods [START_REF] Oh | Robust high dynamic range imaging by rank minimization[END_REF][START_REF] Lee | Ghost-free high dynamic range imaging via rank minimization[END_REF][START_REF] Lee | High dynamic range imaging via truncated nuclear norm minimization of low-rank matrix[END_REF][START_REF] Lee | Computationally efficient truncated nuclear norm minimization for high dynamic range imaging[END_REF], our rank minimization successfully applies to light fields where the rank is expected to be higher than 1 because of the parallax. We finally show the advantage of processing the different viewpoints simultaneously.

RAW DATA PROCESSING 2.1. De-multiplexing and Soft Saturation Detection

The plenoptic RAW data is first de-multiplexed into light field views by adapting the method in [START_REF] Dansereau | Decoding, calibration and rectification for lenselet-based plenoptic cameras[END_REF] (we use the enhancements in [START_REF] David | White lenslet image guided demosaicing for plenoptic cameras[END_REF]). In order to exploit the full dynamic range of each captured light field, we removed the clipping operations previously used for pixels close to saturation. Instead, soft saturation detection is performed directly on sensor data (i.e. before the devignetting and demosaicing steps of [START_REF] Dansereau | Decoding, calibration and rectification for lenselet-based plenoptic cameras[END_REF]). For a pixel of normalized value x, we compute its saturation s(x) = min((x + (1 -τ )) 12 , 1) as illustrated in Fig. 2. In practice, even if a pixel is not fully saturated (x < 1), neighboring pixels associated to the other color components on the bayer pattern may saturate. Hence, the level of full saturation τ is set to 0.9 to account for unreliable colors of high pixel values. Similarly to the color data, the sensor saturation image is demultiplexed using [START_REF] Dansereau | Decoding, calibration and rectification for lenselet-based plenoptic cameras[END_REF] (without devignetting and demosaicing) to obtain a saturation map of each light field view.

Construction of Matrices for WLRA

In order to cope with possible movements of the camera when capturing the different exposures, a homography alignment is performed. For fast computations, homography parameters are determined only from the central views of the high and low exposure light fields employing homography based low rank approximation (HLRA) [START_REF] Jiang | Light field compression with homography-based low-rank approximation[END_REF]. All the views of the low exposure light fields are then aligned to the high exposure one using the same homography. Prior to the WLRA, we convert RGB data to luminance Y and CIE chromaticity components u v . Independent processing of the Y , u and v components has been shown in [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] to be advantageous both in terms of quality and computing speed for light field completion with low rank matrix approximation. Additional conversion of the saturation s into weights w is also necessary to perform weighted low rank approximation:

w = 1 -s s + 1/wmax . (1)
The non-saturated pixels (s = 0) are thus associated to the maximum weight wmax indicating a high confidence, while the fully saturated pixels (s = 1) have a null weight. In our experiments, we use wmax = 100.

Finally, weights and image data are arranged into respective matrices W and M ∈ R m×n such that each column contains a vectorized light field view (n is the total number of views including all exposures and m is the number of pixels per view). In the case where an object is over-exposed even in the lowest exposure, some rows of the matrix may have only unknown elements (zero weight). This may result in arbitrary completed values in these areas. In order to avoid this situation, the weights of all the views of the lowest exposure are set to wmax.

WEIGHTED LOW RANK APPROXIMATION 3.1. General Formulation

Our WLRA method generalizes the matrix completion algorithm in [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] to the case of non-binary weights. While most methods, relying on the results of [START_REF] Candès | Exact matrix completion via convex optimization[END_REF], solve a simpler convex problem by minimizing the nuclear norm (i.e. sum of singular values), the algorithm in [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] keeps the rank in the objective function and solves the problem using the Alternating Direction Method of Multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. We use the same approach as it was shown to outperform conventional matrix completion for the closely related light field inpainting application. Given the matrix M to approximate, a global noise tolerance parameter , and the element-wise weight matrix W, the problem formulation is:

min X rank(X) s.t. X = Z W • (Z -M ) 2 F ≤ , (2)
where the operator • is the element-wise multiplication. The introduction of the matrix Z and the constraint X = Z makes it possible to use the Alternating Direction Method of Multipliers.

The constraint X = Z is taken into account by defining the augmented Lagrangian function:

L(X, Z, Λ, ρ) = rank(X) + Tr(Λ (X -Z)) + ρ 2 X -Z 2 F , (3)
where Λ is a matrix of lagrangian multipliers and ρ is a positive scalar. At each iteration k, the ADMM algorithm for our WLRA problem then consists in updating X, Z, Λ, and ρ as:

X k = arg min X L(X, Z k-1 , Λ k-1 , ρ k-1 ), (4)
Z k = arg min Z s.t. W •(Z-M ) 2 F ≤ L(X k , Z, Λ k-1 , ρ k-1 ), (5)
Λ k = Λ k-1 + ρ k-1 • (X k -Z k ), (6)
ρ k = t • ρ k-1 (with t > 1), ( 7
)
where t is a parameter controlling the tradeoff between accurate minimization (t ≈ 1) and fast convergence (t 1). For our HDR light field application, we use t = 1.1 and t = 4 respectively for the luminance and the chromaticity components.

In the following subsections we describe how to efficiently solve the X and Z sub-problems introduced in Eqs. ( 4) and ( 5) respectively. For simplicity of notation, the iteration indices k are ignored in the rest of the paper.

X sub-problem

Defining Hτ as the singular value hard thresholding operator with threshold τ , it was shown in [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] that the X sub-problem in Eq.( 4) has a closed form solution given by: arg min

X L(X, Z, Λ, ρ) = H 2 ρ Z - Λ ρ .
(8)

Z sub-problem

The Z sub-problem in Eq.( 5), can be equivalently rewritten:

arg min Z N -Z 2 F s.t. W • (Z -M ) 2 F ≤ , (9)
where N = X + Λ ρ . Let us define the Lagrangian function:

L Z (Z, λ) = N -Z 2 F + λ W • (M -Z) 2 F -. ( 10
)
By deriving the stationary condition (i.e. ∂L Z ∂Z (Z) = 0), we obtain the following expression of Z for all entries Z ij :

Z ij = (N ij + λW 2 ij • M ij )/(1 + λW 2 ij ). (11)
In order to determine the value of the lagrangian multiplier λ such that Z is the solution of Eq.( 9), we use the Karush-Kuhn-Tucker complementary slackness conditions:

λ ≥ 0 λ = 0 or f W • (M -Z) 2 F -= 0. (12) (13)
The expression f in Eq.( 13) can be seen as a function of λ by substituting Z into its expression from Eq.( 11). We obtain:

f (λ) = i,j W ij (M ij -N ij ) 1 + λW 2 ij 2 -= 0. ( 14
)
As Eq.( 14) does not have a closed form solution for λ, we solve it numerically. The value of λ is initialized as 0 and is iteratively updated following Newton's method, by subtracting f (λ) f (λ) , where the derivative f of f is given by:

f (λ) = i,j -2W 2 ij 1 + λW 2 ij • W ij (M ij -N ij ) 1 + λW 2 ij 2 . ( 15
)
Note that f is monotonically decreasing on R + . Thus, if f (0) < 0, there is no positive solution to Eq.( 14). In this case, the value of λ satisfying the conditions ( 12), ( 13) as well as the constraint in Eq.( 9) is 0, and no further iteration is required.

Otherwise, there is a unique positive solution. In practice, for our problem, less than 10 iterations are usually required to find λ with the stopping criterion: |f (λ)| < 10 -4 . In order to accelerate each iteration of the λ computation, the elements of the input weight matrix W are quantized to a limited set of values {wq} q∈ 1,qmax . Then, instead of computing f and f by summing terms over all matrix elements i, j, they are determined as a sum of only qmax terms:

f (λ) = qmax q=1 w 2 q (1 + λw 2 q ) 2 • Rq -, (16)
f (λ) = qmax q=1 -2w 2 q 1 + λw 2 q • w 2 q (1 + λw 2 q ) 2 • Rq , (17)
where Rq does not depend on λ and therefore only needs to be computed once as:

Rq = i,j (M ij -N ij ) 2 • δ Wij wq , (18)
with δ the Kronecker delta (equal to 1 if wq = W ij and 0 otherwise). Using this method with a number of weight values qmax = 256, the computing cost of λ is negligible with respect to the rest of the algorithm. Hence, the per-iteration cost of our WLRA method is not increased compared to the algorithm [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] restricted to binary weights. Additional discussion on computing cost is given in the next section.

EXPERIMENTAL RESULTS

For all our experiments, light fields are captured using a Lytro Illum camera. We first validate the choice of the rank minimization formulation by comparing our results with those obtained by replacing the rank in Eq.( 2) by the Truncated Nuclear Norm X r (sum of the singular values excluding the r highest ones). This problem is solved by replacing the hard thresholding in Eq.( 8) by the Partial Singular Value Thresholding operator P r, 1 ρ defined in [START_REF] Oh | Partial sum minimization of singular values in rpca for low-level vision[END_REF]. For both rank and TNN minimizations, Z is initialized as Z 0 = M and we use the parameters = 0.003

• W • M 2 F
for the luminance component, and = 0.04

• W • M 2
F for the u v components. For our rank minimization, the penalty parameter ρ is initialized according to [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] as ρ 0 = 8/(σ 1 + σ 2 ) 2 (with σ 1 and σ 2 the two largest singular values of M ). For TNN minimization, we use ρ 0 = 1.25/σ 1 as suggested in [START_REF] Lee | Computationally efficient truncated nuclear norm minimization for high dynamic range imaging[END_REF]. In single view HDR imaging, existing low rank based methods [START_REF] Oh | Robust high dynamic range imaging by rank minimization[END_REF][START_REF] Lee | Ghost-free high dynamic range imaging via rank minimization[END_REF][START_REF] Lee | High dynamic range imaging via truncated nuclear norm minimization of low-rank matrix[END_REF][START_REF] Lee | Computationally efficient truncated nuclear norm minimization for high dynamic range imaging[END_REF] minimize the TNN with r = 1 as the matrix rank can be assumed to be 1. However, this assumption does not hold for light fields with disparities between views, and produces blur in the completed regions as shown in Fig. 3(c),(g). Using higher values of r keeps the details but not the low spatial frequencies (Fig. 3(a),(b),(e),(f)). Our rank minimization successfully recovers both low and high frequencies (Fig. 3(d),(h)).

Fig. 4 provides further justification of our approach compared to independently applying the WLRA to each viewpoint. By considering all the viewpoints simultaneously our method is able to reduce the noise and thus address the challenging case of large exposure gaps where some saturated areas of the high exposure capture are under-exposed and noisy in the low exposure one. Note that it also possibility to extend to correction of artifacts that are not consistent between views (e.g. lens flare) as long as it can be detected and used to reduce the weights.

We finally illustrate in Fig. 5 the advantage of using soft saturation detection along with our WLRA method instead of the matrix completion of [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] that requires binary weights. Note that in this example, the matrix size is 267030 × 98 (2 exposures for 7 × 7× views of size 430 × 621). With our Matlab implementation, the processing time was 54.4s for WLRA (respectively 103, 8 and 8 iterations for Y,u',v'), and 54.8s for the low rank matrix completion (respectively 104, 8 and 8 iterations for Y,u',v'). Hence, for our application, the proposed extension of matrix completion to non-binary weights does not affect the convergence speed and has similar computing cost. Low Rank Matrix Completion of [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] (binary weights), (b) our WLRA (non-binary weigths). In this example, the white parts of the label were saturated in the high exposure capture. Traditional matrix completion produces artifacts at the boundary of these regions.

CONCLUSION

We have presented a new method for the acquisition of high dynamic range, densely sampled light fields, using successive plenoptic captures at varying exposures. In order to recover the saturated regions of the high exposure light field, we have defined a new weighted low rank approximation algorithm. It exploits simultaneously the information of all the viewpoints at all exposures, as well as soft saturation maps previously extracted from the RAW data. We have shown that, unlike the low rank methods used in single view HDR acquisition, our algorithm remains suitable for processing dense light fields with parallax. Our weighted approach also produces more natural results than traditional low rank matrix completion. One can note that we do not address the case of dynamic scenes if moving objects are partly saturated in the high exposure capture. However, we have shown that the proposed algorithm remains effective with large gaps between exposures, thus simplifying the acquisition task.

As a final note, the simultaneous processing of all the views allows for potential extension to plenoptic devices generating differently exposed or filtered viewpoints in a single capture. Although existing lenslet-based plenoptic cameras already present such characteristics (darker external views caused by vignetting of the micro-lenses), this is not sufficient to significantly increase the dynamic range. Hence, further research is necessary to achieve instantaneous capture of HDR light fields.

ACKNOWLEDGEMENT

This project has been supported in part by the Science Foundation Ireland (SFI) under the Grant Number 15/RP/2776 and in part by the EU H2020 Research and Innovation Programme under grant agreement No 694122 (ERC advanced grant CLIM).

Fig. 1 .

 1 Fig. 1. HDR Light Field acquisition pipeline.

Fig. 2 .

 2 Fig. 2. Soft saturation detection model.

Fig. 3 .

 3 Fig. 3. 13x13 HDR Light field reconstructed from 3 exposures: (a)-(d) show the central view in false color for each method. (e),(f),(g),(h) are the epipolar images corresponding to the blue segment in (a),(b),(c),(d) respectively (vertical axis shows vertical viewpoint change). The epipolar image in (g) shows that the structures completed using TNN minimization with r=1 are blurred and identical in all the views (parallax is not preserved).

  (a) Viewpoints processed independently (b) All viewpoints processed simultaneously

Fig. 4 .

 4 Fig. 4. One view of a HDR light field captured from two significantly different exposures (images are gamma encoded for display purpose). Some over-exposed areas in the high exposure capture are under-exposed in the low exposure one. Independent processing of the viewpoints (a) highlights the presence of noise due to underexposition. Using all the viewpoints in the WLRA (b) naturally exploits redundancies to reduce the noise.

Fig. 5 .

 5 Fig. 5. Recovered saturated areas from two exposures using: (a)


Hamon Arnaud
email: arnaud.hamon@irccyn.ec-nantes.fr

Aoustin Yannick
email: yannick.aoustin]@irccyn.ec-nantes.fr

Walking gait of a planar bipedal robot with four-bar knees

The design of a knee joint is a key issue in robotics and biomechanics to increase the compatibility between prostheses and human movements and to improve the performances of the bipedal robot. We propose a novel design for the knee joint of a planar bipedal robot, based on a four-bar linkage. The advantage of this structure is to produce a translation of the Instantaneous Center of Rotation (ICR) of the knee joint like in the human case. We propose to study the sthenic criterion of a bipedal robot using four-bar knee joint during a walking gait. This walking gait is a succession of finite time double support phases, single support phases and impacts. During the double support phase, both feet rotate. This phase is ended by an impact on the ground of the toe of one foot, the other foot taking off. The single support phase is ended by an impact of the swing foot heel, the other foot keeping contact with the ground through its toe. A parametric optimization problem is presented for the determination of the parameters corresponding to the optimal cyclic walking gait. The main contribution of this paper is to show the influence of double support phases for cyclic stable walking gaits with this novel bipedal robot.

Introduction

The researchers in biomechanics have improved a lot the understanding of the human lower limb and especially the knee joint [START_REF] Wilson | Ligaments and articular contact guide passive knee flexion[END_REF] and the ankle joint [START_REF] Leardini | A geometric model of the human ankle joint[END_REF]. Indeed, these two joints have a complex architecture formed by non symmetric surfaces. This common architecture of the knee has not really changed the last 300 million years despite an important diversity of the functional need [START_REF] Dye | An evolutionary perspective of the knee[END_REF]. The motions of the femur with respect to the tibia are limited due to many ligaments and the patella. In addition to the flexion in the sagittal plane, there is an internal rotation with a displacement of the Instantaneous Center of Rotation (ICR) of the knee joint and a posterior translation of the femur on the tibia. The mobility of the knee joint consists of three important rotations and two translations, which are the flexion-extension, adduction-abduction, and external-internal rotations, and translation in the sagittal plane, see [START_REF] Farhat | musculo-skeketic model for knee joint forces estimation in sport activities[END_REF]. These motions are guided by the cruciate ligaments and the articular contacts, see [START_REF] Frank | Anatomy of the cruciate ligaments and their function in extension and flexion of the human knee joint[END_REF] and [START_REF] Wilson | Ligaments and articular contact guide passive knee flexion[END_REF]. Then they cannot be only represented by one or two revolute joints. Different studies have confirmed these results by an observation of the motions of the human knee in the 3D space [START_REF] Landjerit | Cinématique spatiale de l'articulation fémoro-tibiale du genou humain: caractérisation expérimentale et implications chirurgicales[END_REF].

To reproduce the movement of the knee several models exist. Freudenstein et al [START_REF] Freudenstein | Kinematics of human knee joint[END_REF] proposed a mathematical description, based on the Euler-Savary equation and relaters concepts, of the tibiofemoral relative planar motion. The femoral condyle was approximated in [START_REF] Kurosawa | Geometry and motion of the knee for implant and orthotic design[END_REF] by a spherical joint, from measurement carried out on the posterior portion. Delp et al [START_REF] Delp | An interactive graphicsbased model of the lower extremity to study orthopaedic surgical procedures[END_REF] chose an ellipse in the sagittal plane to approximate the femoral condyle. In 1974, A. Menschik [START_REF] Menschik | Mechanics of the knee-joint. part i[END_REF] proposed to represent the knee joint by a crossed four-bar linkage. A classical polycentric knee is the crossed four-bar linkage [START_REF] Gard | The influence of four-bar linkage knees on prosthetic swing-phase floor clearance[END_REF], which is used in most prosthetic knees. This mechanism is the dual solution of the classical four-bar knee but does not have any singularity in the range of motion typically used. As a matter of fact, the kinematic singularities of the mechanism may limit its range of motion, mainly when an important knee flexion is required. The dimensions of the crossed four-bar linkage can be chosen by measuring, on a real subject: (i) the length of the anterior and posterior ligaments; (ii) the position of the crossed ligament attachments on the tibia and the femur, projected on the sagittal plane in the maximum extension position [START_REF] Feikes | A constraint-based approach to modelling the mobility of the human knee joint[END_REF]. Authors in [START_REF] Farhata | A procedure for estimating the relevant forces in the human knee using a four-bar mechanism[END_REF] developed a dynamic model for the lower extremity capable of estimating forces generated in the knee. The proposed model considers on a four-bar mechanism in the knee, a spherical joint in pelvis and a revolute one in the ankle. Authors proved experimentally through a force platform and high-speed cameras that the crossed four-bar mechanism could reproduce the knee flexion-extension and the two main translations. As a consequence for the knee the adduction-abduction and external-internal rotations are not modeled. However, these angles are relatively small when compared with the flexion-extension angle.

Roboticists have come up with new and better bipedal robots recently. For instance, HRP-2LR [START_REF] Kajita | A running controller of humanoid biped hrp-2lr[END_REF] and RABBIT [START_REF] Chevallereau | Rabbit: a testbed for advanced control theory[END_REF], which are able to run and and to be dynamically stable. A. Grishin et al. [START_REF] Grishin | Dynamic walking of a vehicle with two telescopic legs controlled by two drives[END_REF] focused on the design of a bipedal robot with telescopic legs. T. Yang et al. [START_REF] Yang | Design and control of a planar bipedal robot ernie with parallel knee compliance[END_REF] used a compliant parallel knee to improve the walking motion. Some authors also dealt with the walking and running gaits using the toe rotation [START_REF] Kajita | Zmp-based biped running enhanced by toe springs[END_REF], [START_REF] Tajima | Fast running experiments involving a humanoid robot[END_REF], [START_REF] Romero | Design of a walking cyclic gait with single support phases and impacts for the locomotor system of a thirteen-link 3d biped using the parametric optimization[END_REF]. Our objective is to improve the bipedal robot performance thanks to a new design of the knee joint. Several papers deal with the bipedal robots equipped with complex knees, like G. Gini et al. [START_REF] Gini | Human-oriented biped robot design : insights into the development of a truly antropomophic leg[END_REF] which used knee joints based on the human knee surfaces. F. Wang et al. [START_REF] Wang | Design and implementation of coordinated control strategy for biped robot with heterogeneous legs[END_REF] developed a bipedal robot with two different joints, namely, a revolute joint and a four-bar linkage. We proved in [START_REF] Hamon | Cross four-bar linkage for the knees of a planar bipedal robot[END_REF] that a knee based on a four-bar linkage is better than a knee designed with a revolute joint in terms of energy consumption. With a good choice of the length of the links of the four-bar and its cross configuration the flexion of the knee joint is usually not too limited with the kinematic singularities.

This paper aims to study the performance of a planar bipedal robot, equipped with knees based on four-bar linkages, during a complex walking gait. This gait is composed of finite time double support phases, single support phases and impacts. During the double support phase, both feet rotate. This phase ends with a contact of one toe with the ground while the other foot is taking off. The single support phase ends with an impact of the swing foot heel, the toe of the other foot being in contact with the ground. We also present the dynamic model of the planar bipedal robot whose knees are composed of a four-bar linkage. We developed a parametric optimization method to define a set of optimal reference trajectories for both cyclic gaits. We use this optimization process to study the effect of the double support for cyclic walking gaits.

This paper is organized as follows. Section 2 presents the novel planar bipedal robot whose knees are based on four-bar linkages and its dynamic models. Section 3 deals with the trajectory planning. Section 4 proposes numerical results about the evaluation of the energy consumption as a function of the bipedal robot velocity. Finally, Section 5 presents some conclusions and future works.

Presentation of the bipedal robot with knees composed of a four-bar linkage

Let us introduce the bipedal robot, which is depicted in Figure 1. The dimensions of the four-bar linkage are chosen with respect to the human characteristics measured by J. Bradley et al. through radiography in [START_REF] Bradley | Orientation of the cruciate ligament in the sagittal plane[END_REF] and in accordance with the hydroid geometrical parameters to avoid any singularity in the knee joints during the walking gaits.

Figures 1; 2(a) and 2(b) depict the bipedal robot under study and its four-bar knee linkage. Figure 2(a) represents the four-bar knee linkage. The angular variable α 1 is the actuated variable of the four-bar linkage.

The bipedal robot is equipped with two closed-loop knees. Let us introduce the constraint equations solving the dynamic model [START_REF] Khalil | Modeling, identification and control of robots[END_REF]. Equations for knee joints 1 and 2 are similar. For a sake of clarity we consider knee joint 1 only. The equations of the closed-loop geometric constraints are defined as follows:

a cos q 1 -b sin q g11 + c cos q 2 + d sin q g12 = 0 a sin q 1 + b cos q g11 -c sin q 2 -d cos q g12 = 0 (1) Their first and second time derivatives are -a q1 sin q 1 -b qg11 cos q g11 -c q2 sin q 2 + d qg12 cos q g12 = 0 a q1 cos q 1 -b qg11 sin q g11 -c q2 cos q 2 + d qg12 sin q g12 = 0

q 5 q 3 q 2 q 1 q 4 q p1 q p2 Γ p1 Γ 3 Γ p2 Γ 2 s 5 s 1 s 3 s 4 s 2 Γ 4 Γ 1
(2) and -aq 1 sin q 1 -bq g11 cos q g11 -cq 2 sin q 2 + dq g12 cos q g12 -a q2

1 cos q 1 + b q2 g11 sin q g11 -c q2 2 cos q 2 -d q2 g12 sin q g12 = 0 aq 1 cos q 1 -bq g11 sin q g11 -cq 2 cos q 2 + dq g12 sin q g12 -a q2

1 sin q 1 -b q2 g11 cos q g11 + c q2 2 sin q 2 + d q2 g12 cos q g12 = 0.

(

Through the virtual work principle, these constraint equations can be expressed in the dynamic model by adding the Lagrange multipliers J t 1e λ. Here J 1e is the 2 × 13 Jacobian matrix such as equations ( 2) and ( 3) can be rewritten under the compact forms

J 1e ẋ = 0 (4)
and

J 1e ẍ + J1e ẋ = 0. ( 5
)
and vector λ = f c1 = [f x1 , f y1 ] t defines the two Lagrange multipliers vector associated to the constraint of the closed-loop linkage. These multipliers represent the exerted forces on point A of the knee joint 1 (Fig. 2(a)). The generalized vector x is such as

x = [q p1 , q p2 , q 1 , q 2 , q 3 , q 4 , q 5 , q g11 , q g12 , q g21 , q g22 , x h , y h ] t .

We apply the same principle for the knee joint 2 to obtain the dynamic model of the bipedal robot with the four-bar knees:

A e (x)ẍ + h e (x, ẋ) = D Γe J t 1e J t 2e Γ m f c + J t r1 r 1 m 1z + J t r2 r 2 m 2z (6)
with the constraint equations:

J ri ẍ + Jri ẋ = 0 for i = 1 to 2 J 1e J 2e ẍ + J1e J2e ẋ = 0 (7)
Γ m = [Γ p1 , Γ p2 , Γ 1 , Γ 2 , Γ 3 , Γ 4 ]
t is the vector of the applied joint torques, r i m iz t , with i = 1 to 2, are the resultant wrenchs of the contact efforts with the ground reaction in both feet, and

f c = [f t c1 , f t c2 ] t .
Torques Γ 1 and Γ 4 respectively actuate joint D of knee 1, see (Fig. 2(a)), and the similar joint of knee 2. where x h and y h are the hip coordinates. J r1 and J r2 are the 3 × 13 jacobian matrix of the constraint equations in position and orientation for the two feet, respectively. A e (x) is the 13×13 symmetric positive definite inertia matrix, h e (x, ẋ) is the 13 × 1 vector, which groups the centrifugal, Coriolis effects, and the gravity forces. D Γ is a 13 × 6 matrix composed of the 0 and ±1 given by the principle of virtual works [START_REF] Appell | Dynamique des Systèmes -Mécanique Analytique[END_REF].

(a) q 2 q g11 α 1 -q g12 ICR Δ -q 1 P c A B C D (b)
This dynamic model ( 6) with constraints ( 7) is valid in single and double support phases. During a single support phase the ground reaction force is zero on the swing foot.

Reduced model in single support

The aim is to propose a dynamic model with an implicit liaison of the stance foot, here assumed to be foot 1, with the ground to calculate the torques during the optimization process, with the knowledge of the reference trajectories for the generalized coordinates. This reduced dynamic model is only valid if the stance foot does not take off and there is no sliding during the swing phase.

Then, during the single support phase, the stance foot is assumed to remain in flat contact on the ground, i.e., there is no sliding motion, no take-off, no rotation (q p1 = 0). We can use a new generalized vector q = [q p2 , q 1 , q 2 , q 3 , q 4 , q 5 , q g11 , q g12 , q g21 , q g22 ] t . The reduced dynamic model does not depend on the ground reaction force, which is applied in the stance flat-foot. The dynamic model in single support phase for the biped equipped with the four-bar knees is given by the simplification of the dynamic model ( 6) and (7):

A(q)q + h(q, q) = D Γ J t 1 J t 2 Γ m f c . ( 8
)
with the constraint matrix equation,

J 1 J 2 q + J1 J2 q = 0. ( 9
)
Here vector f c is such as

f c = [f t c1 , f t c2 ] t .
To define the constraints about the ground reaction, the no take-off of the stance flat-foot during the optimization process, we recall the calculation of position of the Zero Momentum Point. The resultant force r 1 of the ground reaction can be calculated by applying the second Newton law at the center of mass of the biped:

mγ = r 1 + mg, ( 10
)
where m is the global mass of the biped, γ = [ẍ g , ÿg ] t are the horizontal and vertical components of the acceleration for its center of mass in the world frame. g = [0, -g] t is the vector of the acceleration of the gravity. This equation allows directly to get r during the single support. Assuming the center of mass G f of the foot has for coordinates (s px , s py ), see Fig. 3. Let m f be the mass of the foot. For the planar biped the coordinate of the Zero Moment Point ( ZM P ), which is merged with the centre of pressure to keep the flat contact, can be obtained through the calculation of the global equilibrium of the bipedal robot around axis z, which gives (Vukobratovic [START_REF] Vukobratovic | On the stability of anthropomorphic systems[END_REF]):

l ZMP = Γ p1 + s px m f g -H p r x r y , (11)
To calculate the applied joint torques, which are used to obtain the sthenic criterion of the bipedal robot during the optimization process, we use the dynamic model [START_REF] Delp | An interactive graphicsbased model of the lower extremity to study orthopaedic surgical procedures[END_REF].

We assume the friction effects due to the four-bar mechanism are negligible with respect to those in the gearbox of the actuators. Then only the performances of actuators are considered to study the eof volition sthenic criterion for the biped. No friction terms are included in the model.

Impact model

During the biped's gait, impacts occur, when the sole, the heel or the toe of the swing foot touches the ground. Let T be the instant of an impact. We assume that the impact is absolutely inelastic and that the foot does not slip. Given these conditions, the ground reactions at the instant of an impact can be considered as impulsive forces and defined by Dirac delta-functions r j = i j δ(t -T ) (j = 1, 2). Here i j = [i jx , i jy ] t is the magnitudes vector of the impulsive reaction in foot j (see [START_REF] Formalskii | Locomotion of Anthropomorphic Mechanisms[END_REF]). Impact equations can be obtained through integration of the matrix motion equation ( 6) for the infinitesimal time from T -to T + . The torques provided by the actuators in the joints, Coriolis and gravity forces have finite values. Thus they do not influence the impact. Consequently the impact equations can be written in the following matrix form:

A e (x(T ))( ẋ+ -ẋ-) = J t 1e J t 2e i fc + J t r1 i 1 + J t r2 i 2 (12)
Here x(T ) denotes the configuration of the biped at instant t = T , (this configuration does not change at the instant of the impact), ẋ-and ẋ+ are respectively the velocity vectors just before and just after an inelastic impact. Vectors i 1 , i 2 and i fc = [i t fc1 , i t fc2 ] t are respectively the impulsive wrenchs in both feet and the impulsive forces in both knees. To take into account of the closed-loop of the four-bar knee linkage we have to complete [START_REF] Feikes | A constraint-based approach to modelling the mobility of the human knee joint[END_REF] with:

J 1e J 2e ẋ+ = 0 (13)
The velocity of the contact part of the stance foot (j = 1) before an impact is null.

J r1 ẋ-= 0 (14)
The swing foot (j = 2) after the impact becomes a stance foot. Therefore, the velocity of its contact part with the ground becomes zero after the impact,

J r2 ẋ+ = 0 (15)
Generally speaking, two results are possible after the impact, if we assume that there is no slipping of the stance feet. The stance foot lifts off the ground or both feet remain on the ground. In the first case, the vertical component of the velocity of the taking-off foot just after the impact must be directed upwards. Also there is no interaction (no friction, no sticking) between the taking-off foot and the ground. The ground reaction in this taking-off leg tip must be null. In the second case, the stance foot velocity has to be zero just after the impact. The ground produces impulsive reactions (generally, i j = 0, j = 1, 2) and the vertical components of the impulsive ground reactions in both feet are directed upwards. For the second case, the passive impact equation ( 12) must be completed by one matrix equations.

J r1 ẋ+ = 0 (16)
In general, the result of an impact depends on two factors: the biped's configuration at the instant of an impact and the direction of the swing foot velocity just before impact [START_REF] Formalskii | Locomotion of Anthropomorphic Mechanisms[END_REF]. After an impact for a biped, there are two possible phases: a single support or a finite time double support. The resolution of the system composed ( 12), ( 13), [START_REF] Frank | Anatomy of the cruciate ligaments and their function in extension and flexion of the human knee joint[END_REF] and eventually [START_REF] Gard | The influence of four-bar linkage knees on prosthetic swing-phase floor clearance[END_REF] gives the velocity vector ẋ+ just after the impact, the impulsive reaction wrenchs i 1 , i 2 and the impulsive forces i fc = [i t fc1 , i t fc2 ] t relatively to the velocity vector ẋ-just before the impact.

To calculate the position of the ZM P at the impact with flat-foot, we have to take into account of the impulsive ground reaction in the global equilibrium of the stance foot and the result is:

l ZMP = - H p i jx i jy . ( 17
)
3 Gait optimization for the cyclic walking

Principle

The studied trajectories are composed of double support phases and single support phases. During the double support phase, we denote a rotation of the forward foot around the heel and a rotation of the backward foot around the toe. The double support phase is ended by the flat contact of the forward foot on the ground. The single support phase is ended by the contact of the swing foot on the toe. The biped with a stance flat-foot is driven by six torques. Its configuration is given in double support phase with five generalized coordinates only because the positions of the front heel and the rear toe are constant. These five independent coordinates are the orientation of feet q p1 and q p2 , the orientation of one of the leg 1 q 1 , q 2 and the inclination of the torso q 5 . For these coordinates, cubic spline functions θ i (t), i = 1, ..., n j = 5 are defined with an initial value in position and velocity and a final value in position and velocity, with two selected knots, n = 2. Thus, for each joint we need to define four parameters of position and velocity to design the trajectories. The final configuration of the biped is with one flat-foot on the ground. Thus only four independent variables are necessary to define it. We use q p2 , q 1 , q 2 and the inclination of the torso q 5 . Only five components of the joint velocities are independent during this phase, thus the final velocity can be defined with five variables at maximum. At the end of this phase from the impact, which occurs on the front toe, six initial angular velocities can be determined for the next single phase. When functions θ i (t), i = 1, ..., 5 are chosen, the joint velocity and the joint acceleration can be deduced by the derivation of the polynomial function. For a given motion, the torques required could be calculated solving the inverse dynamic model [START_REF] Chevallereau | Optimal reference trajectories for walking and running of a biped robot[END_REF]. However, there are 14 unknown variables (six torques, four components for the internal forces in both knees and four components for both ground reactions), for 13 scalar equations only. Consequently, we have to prescribe one unknown variable as parameter. To highlight this point let us consider the global equilibrium in translation and rotation of the biped (18), with Fig. 5. The center of mass G of the biped is located with the coordinates x g and y g , r 1 and r 2 are the ground reaction acting in the front hell and the rear toe, δ g is the dynamic momentum of the biped with respect to its center of mass. We have four unknown variables, r 1x , r 2x , r 1y and r 2y for three equations only.

   y g (r 1x + r 2x ) + (dx g )r 2y + x g r 1y = δ g r 1x + r 2x = mẍ g r 1y + r 2y -mg = mÿ g [START_REF] Gini | Human-oriented biped robot design : insights into the development of a truly antropomophic leg[END_REF] If the sum r 1x + r 2x is unique, there are an infinity of solutions for r 1x or r 2x which satisfy [START_REF] Gini | Human-oriented biped robot design : insights into the development of a truly antropomophic leg[END_REF]. As a consequence r 2x can be defined as a second order polynomial function in time:

r 2x (t) = d 0 + d 1 t + d 2 t 2 (19)
The number of parameter is equal to 29. P = q 1 (T DS ), q 2 (T DS ), q 5 (T DS ), q p2 (T DS ), q p2 ( TSS 2 ), q 1 ( TSS 2 ), q 2 ( TSS 2 ), q 3 ( TSS 2 ), q 4 ( TSS 2 ), q 5 ( TSS 2 ), qp1 (T DS ), qp2 (T DS ), q1 (T DS ), q2 (T DS ), q5 (T DS ), x h (T SS ), y h (T SS ),

q 5 (T SS ), d, qp2 (T SS ), q1 (T SS ), q2 (T SS ), q3 (T SS ), q4 (T SS ), q5 (T SS ), d 0 , d 1 , d 2 , T DS ] (20)
The variables T SS and T DS are the duration of the single support and the double support respectively for the cyclic walking gait 2.

Parametric optimization problem

Basis functions as polynomial functions, see [START_REF] Channon | Derivation of optimal walking motions for a bipedal walking robot[END_REF], truncated Fourier series, see [START_REF] Cabodevilla | Near optimal gait for a biped robot[END_REF] or cubic spline function, see [START_REF] Romero | Design of a walking cyclic gait with single support phases and impacts for the locomotor system of a thirteen-link 3d biped using the parametric optimization[END_REF] can be chosen for the parametric optimization algorithm. A basis function of high order can lead to oscillatory phenomenon during the optimization process, see see [START_REF] Chevallereau | Optimal reference trajectories for walking and running of a biped robot[END_REF] or [START_REF] Hu | Biped gait optimization using spline function based probability model[END_REF]. By parameterizing the joint motion in terms of cubic spline functions, the optimization problem is reduced to a constrained parameter optimization problem of the form:

Minimize C W (P) subject to g j (P) ≤ 0 for j = 1, 2, ...l [START_REF] Hu | Biped gait optimization using spline function based probability model[END_REF] where P is the set of optimization variables. C W (P) is the criterion to minimize with l inequality constraints g j (P) ≤ 0 to satisfy. Cubic function of third-order are chosen for the single and double support phase with two selected knots. Thus the cubic spline are uniquely defined by specifying an initial configuration, an initial velocity, a final configuration and a final velocity. The criterion and constraints are given in the following sections. We used the SQP method (Sequential Quadratic Programming) [START_REF] Gill | Practical optimization[END_REF], [START_REF] Powell | Variable metric methods for constrained optimization[END_REF] with the fmincon function of Matlab R to solve this problem.

The criterion

Many criteria can be used to produce an optimal trajectory. A sthenic criterion is chosen to obtain optimal trajectories:

C W = 1 d T 0 6 i=1 Γ t mi Γ mi dt, ( 22
)
where T is the step duration and Γ mi is the torque applied in joint i. During an optimization process the step length d is an optimization variable and the walking speed v is fixed, such as the step duration is directly given through the relation

T = d/v.
The resulting optimal control is continuous and cancels the risks of a jerky functioning [START_REF] Chevallereau | Bipedal Robots[END_REF]. This smoothness property also guarantees a better numerical efficiency for the algorithm used for the optimization problem-solving

The constraints

Two types of constraints are used to obtain a realistic gait.

• The necessary constraints, which ensure a valid walking gait.

The first constraint ensures the supporting leg tip does not take off or slide on the ground. So, the ground reaction force is inside a friction cone, defined with the coefficient of friction f :

max(-f r iy -r ix ) ≤ 0 max(-f r iy + r ix ) ≤ 0 (23)
j = 1 or 2. r x and r y are the normal and tangential component of the reaction force. Moreover, we can introduce a constraint on the ground reaction at the impact:

(-f i 1y -i 1x ) ≤ 0 (-f i 2y + i 2x ) ≤ 0 (24)
To ensure the no rotation of the supporting foot we introduce a constraint on the ZM P during the single support phase and at the instant of the impact:

(l p -L p ) ≤ l ZMP ≤ l p (25)
Here L p is the length of the foot and l p is the distance between the heel and the ankle along the horizontal axis, see Fig. 3.

Just after the impact, the velocity of the taking-off foot should be directed upward. In consequence, the positivity of the vertical component of the velocities for the heel and the toes is added to the set of constraints.

The last constraint allows to ensure the non penetration of the swinging foot in the ground.

• The unnecessary constraints, which ensure a technological realistic gait.

We introduced mechanical stops on the joint variables. Moreover, we limited the torques with a constraint, which sets a template of the maximum torque of the motor relatively to the velocity [START_REF] Chevallereau | Rabbit: a testbed for advanced control theory[END_REF].

Results

With this optimization process, a set of optimal trajectories is developed for different walking velocities for a biped equipped of four-bar knees. We use for the SQP algorithm the fmincon function of Matlab R Ṫhe SQP algorithm solves a local optimization problem. However we tried success different initial conditions, which satisfy the constraints defined in the previous section 3.2.2. We can see Fig. 6.a a comparison of the sthenic criterion according to the walking velocities for the biped with the four-bar knees and with the revolute joint knees. Despite the local optimization problem we can see the the four-bar knees offer a benefice with respect to the sthenic criterion: for all the walking velocities the sthenic criterion is bigger with the revolute joint knees than with the four-bar knees. For the biped with the four-bar knees we can see Fig. 7.a the ratio between the part of the sthenic criterion during the double support phase and the sthenic criterion and between the part of the sthenic criterion during the single support phase and the sthenic criterion. For the slow walking gaits the most part of the sthenic criterion is consumed during the double support phase. For the fast walking gaits the ratio of the part of the sthenic criterion between the single support phase and the sthenic criterion increases. It is logical from a physical point of view. The parts of the sthenic criterion consumed by the biped during the double support and single support phases are quasi equal for around a velocity of 4 Km/h. For walking velocities between 1.5 Km/h and 2.5 Km/h the motion of the transfer leg with a four-bar knee is quasi passive, Fig. 7.b. Its consumed part of the sthenic criterion is very small with respect to the sthenic criterion. It is one of the remarkable characteristics of the walking gaits for human [START_REF] Formalskii | Locomotion of Anthropomorphic Mechanisms[END_REF]. We do not observe this property with a walking gait, which is only composed of single support phases and impacts, see [START_REF] Hamon | Cross four-bar linkage for the knees of a planar bipedal robot[END_REF]. Then it is interesting to insert a double support phase to come closer to the anthropomorphic walking gait. and9.c show the vertical sides of the heel and big toe of the left foot. The balancing phase takes place between 0.6 s and 0.9 s. For this cycle both feet are grounded for 0-15% of the cycle. the leg, which is the swing leg during the single support. Figure 10.a shows the flexion and the extension movements according to the time measurement of the four-bar knee in the sagittal plane during the step of a walking gait. Figures 9.b and 9.c show the vertical sides of the heel and big toe of the foot. The balancing phase takes place between 0.4 s and 0.6 s. Of course there are many differences between the results from the numerical test and the human data, specially for the knee joint. However despite that the size of the walking man is different from those of the bipedal robot the simulated walking gait leads to similar scales of size with respect to human.

Conclusions

We have studied a closed-loop linkage for the knees of a planar bipedal robot. We have produced with a parametric optimization method a set of optimal reference walking trajectories for this bipedal robot. The walking motions are composed of single support and double support phases separated by impulsive impact. The benefice on the sthenic criterion of a four-bar knee is real with respect to the revolution joint because the needed torques for a similar walking gaits are smaller with the four-bar knees. Then in practice, it is possible to reduce the size of the actuators and to decreased the energy consumption. The study of the sthenic criterion of this biped has shown interesting results in term of energy distribution, between the single support phase and the double support phase, as a function of the walking velocity, which can be compared with the characteristics of the human walking gait. The characteristics of the walking gait of the bipedal robot equipped with four-bar knees are still far from those of human. But the human walking gait is more complex with a rolling sole on the ground of the stance foot than our defined walking gait. Furthermore a great attention about the optimal dimensions of the four-bar mechanism is necessary.

The perspective of this study is to extend this work for a 3D bipedal robots and to test this new structure with a experimental bipedal robots. Moreover, a measure of the evolution of the center of rotation of the human knee will be utilized to choose with a better accuracy the dimensions of the fourbar mechanism. Finally, a single support phase with rotation of the support foot can be introduced to come closer to the results of an anthropomorphic walking motion.

Figure 1 :

 1 Figure 1: Schematic of a planar bipedal robot. Absolute angular variables and torques.

Figure 2 :

 2 Figure 2: Details of the four-bar joint and position of the Instantaneous Center of Rotation (ICR)

Figure 3 :

 3 Figure 3: Details of the foot.

Figure 4 :

 4 Figure 4: Presentation of the studied walking gait.

2 GFigure 5 :

 25 Figure 5: Ground reactions in double support phase and the center of mass of the biped.

Figure 6 :

 6 Figure 6: a): Evolution of the sthenic criterion of bipedal robot according to the walking velocity, (o) revolute joint knees, (+) four-bar knees.

Figure 7 :

 7 Figure7: a): Evolution of the sthenic criterion of bipedal robot according to the walking velocity, (o) revolute joint knees, (+) four-bar knees. a) Ratio between the part of the sthenic criterion during the double support phase and the sthenic criterion (plus marker) and between the part of the sthenic criterion during the single support phase and the sthenic criterion (round marker) as a function of the walking velocity. b) Evolution of the part of the sthenic criterion for the transfer leg.

Figure 8 .

 8 Figure 8.a presents the evolution of the double support and single support duration during a walking gait as a function of the walking velocity.Figure 8.b gives the evolution of the step length as a function of the walking velocity. From Fig. 8.a and Fig. 8.b we observed that, when the walking velocity increases, the duration of the single support phase decreases. The durations of the double support phase and the single support phase tend to be equal when the walking velocity increases.

Figure 8 .

 8 Figure 8.a presents the evolution of the double support and single support duration during a walking gait as a function of the walking velocity.Figure 8.b gives the evolution of the step length as a function of the walking velocity. From Fig. 8.a and Fig. 8.b we observed that, when the walking velocity increases, the duration of the single support phase decreases. The durations of the double support phase and the single support phase tend to be equal when the walking velocity increases.

Figure 8 :

 8 Figure 8: a): Evolution of the double support phase duration (plus marker) and the single support phase duration (round marker) as a function of the walking velocity. b) Evolution of the step length as a function of the walking velocity.

Figure 9 .

 9 Figure 9.a 1 shows the flexion and the extension movements according to the time measurement of the left knee in the sagittal plane during the same cycle of a walking male of a height of 1.7 m. Figures 9.band 9.c show the vertical sides of the heel and big toe of the left foot. The balancing phase takes place between 0.6 s and 0.9 s. For this cycle both feet are grounded for 0-15% of the cycle.

Figure 9 :

 9 Figure 9: a): Flexion and extension of the left knee according to the time measurement of the left knee. b) Vertical side of the heel according to the time measurement. c) Vertical side of the big toe according to the time measurement.

Figures 10 are devoted

  Figures 10 are devoted to the leg, which is the swing leg during the single support. Figure10.a shows the flexion and the extension movements according to the time measurement of the four-bar knee in the sagittal plane during the step of a walking gait. Figures 9.b and 9.c show the vertical sides of the heel and big toe of the foot. The balancing phase takes place between 0.4 s and 0.6 s. Of course there are many differences between the results from the numerical test and the human data, specially for the knee joint. However despite that the size of the walking man is different from those of the bipedal robot the simulated walking gait leads to similar scales of size with respect to human.

  Figures 10 are devoted to the leg, which is the swing leg during the single support. Figure10.a shows the flexion and the extension movements according to the time measurement of the four-bar knee in the sagittal plane during the step of a walking gait. Figures 9.b and 9.c show the vertical sides of the heel and big toe of the foot. The balancing phase takes place between 0.4 s and 0.6 s. Of course there are many differences between the results from the numerical test and the human data, specially for the knee joint. However despite that the size of the walking man is different from those of the bipedal robot the simulated walking gait leads to similar scales of size with respect to human.

Figure 10 :

 10 Figure 10: a): Flexion and the extension of the four-bar knee according to the time measurement. b) Vertical side of the heel according to the time measurement. c) Vertical side of the big toe according to the time measurement.

Table 1

 1

	gathers the physical data of

Table 1 :

 1 Physical parameters of the bipedal robot.

The following measurements were communicated to us by Franck Multon, who is part of the Physiology and Muscular Exercise Laboratory of the University of Rennes

and the ENS of Cachan