H Becker

L Albera

P Comon

R Gribonval

F Wendling

I Merlet

A PERFORMANCE STUDY OF VARIOUS BRAIN SOURCE IMAGING APPROACHES

Keywords: EEG, MEG, inverse problem, source localization, priors

The objective of brain source imaging consists in reconstructing the cerebral activity everywhere within the brain based on EEG or MEG measurements recorded on the scalp. This requires solving an ill-posed linear inverse problem. In order to restore identifiability, additional hypotheses need to be imposed on the source distribution, giving rise to an impressive number of brain source imaging algorithms. However, a thorough comparison of different methodologies is still missing in the literature. In this paper, we provide an overview of priors that have been used for brain source imaging and conduct a comparative simulation study with seven representative algorithms corresponding to the classes of minimum norm, sparse, tensor-based, subspace-based, and Bayesian approaches. This permits us to identify new benchmark algorithms and promising directions for future research.

INTRODUCTION

Electroencephalography (EEG) and magnetoencephalography (MEG) are two multi-sensor systems that are routinely used for the diagnosis and management of some diseases such as epilepsy or for the understanding of the brain functions in neuroscience research. The EEG/MEG measurements (X ∈ R N ×T ), acquired with N sensors positioned on the scalp for T time samples, constitute a linear mixture of brain sources (S ∈ R D×T ) in the presence of noise (N ∈ R N ×T ):

X = GS + N. ( 1
)
The mixture is characterized by the lead field matrix (G ∈ R N ×D ), which can be computed given a model of the head H. Becker was supported by Conseil Régional PACA and by CNRS. The work of P. Comon was funded by the FP7 European Research Council Programme, DECODA project, under grant ERC-AdG-2013-320594. The work of R. Gribonval was funded by the FP7 European Research Council Programme, PLEASE project, under grant ERC-StG-2011-277906. We also acknowledge the support of Programme ANR 2010 BLAN 0309 01 (project MULTIMODEL). and a predefined source space that is spanned by D dipoles with fixed positions inside the brain [START_REF] Gramfort | Mapping, timing and tracking cortical activations with MEG and EEG: Methods and application to human vision[END_REF]. We assume in this paper that the source space consists of dipoles that are located on the cortical surface with orientations perpendicular to this surface [START_REF] Dale | Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach[END_REF]. The objective of brain source imaging consists in reconstructing the sources S based on the surface measurements X. However, as the number of dipoles (∼ 10000) generally exceeds the number of sensors (∼ 100), the source imaging problem is ill-posed. In order to restore identifiability, additional hypotheses about the sources have to be made.

After several decades of research, giving rise to an impressive number of algorithms [START_REF] Baillet | Electromagnetic brain mapping[END_REF][START_REF] Michel | EEG source imaging[END_REF][START_REF] Grech | Review on solving the inverse problem in EEG source analysis[END_REF][START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF], the source imaging problem is still one of the major challenges in biomedical engineering. Based on various priors, the proposed techniques pursue different methodological approaches, comprising minimum norm, sparse, tensor-based, Bayesian, and subspacebased methods. However, a thorough comparison including algorithms of all these approaches is still missing. In this paper, we fill this gap by providing a classification of different algorithms based on the exploited priors and conducting a simulation study in which we compare seven representative algorithms of different classes: sLORETA [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sLORETA): technical details[END_REF], MCE [START_REF] Uutela | Visualization of magnetoencephalographic data using minimum current estimates[END_REF], VB-SCCD [START_REF] Ding | Reconstructing cortical current density by exploring sparseness in the transform domain[END_REF], MxNE [START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF], STWV-DA [START_REF] Becker | EEG extended source localization: tensor-based vs. conventional methods[END_REF], Champagne [START_REF] Wipf | Robust Bayesian estimation of the location, orientation, and time course of multiple correlated neural sources using MEG[END_REF], and 4-ExSo-MUSIC [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF]. This permits us to identify new benchmark algorithms and directions for further research.

HYPOTHESES

The hypotheses that are exploited to solve the brain source imaging problem can be distinguished into three categories: hypotheses that apply to the spatial, the temporal, and the spatio-temporal distribution of the sources. In the following, we provide an overview and a short description of the priors which have been used in brain source imaging.

Hypotheses on the spatial distribution S1) Minimum energy

The power of the sources is physiologically limited. A popular approach thus consists in identifying the spatial distribution of minimum energy [START_REF] Hämäläinen | Interpreting magnetic fields of the brain: Estimates of current distributions[END_REF][START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sLORETA): technical details[END_REF]. S2) Minimum energy in a transform domain As the spatial distribution of the sources is unlikely to contain abrupt changes, it can be assumed to be smooth [START_REF] Pascual-Marqui | Low resolution electromagnetic tomograph: A new method for localizing elecrical activity in the brain[END_REF]. In practice, this is achieved by constraining the Laplacian of the source spatial distribution to be of minimum energy.

S3) Sparsity

In practice, it is often reasonable to assume that only a small fraction of the source dipoles contributes to the measured signals of interest in a significant way. The signals of the other source dipoles are thus expected to be approximately zero, leading to the concept of sparsity [START_REF] Uutela | Visualization of magnetoencephalographic data using minimum current estimates[END_REF]. S4) Sparsity in a transformed domain If the number of active dipoles exceeds the number of sensors, the source distribution is not sufficiently sparse for standard methods based on sparsity to yield accurate results, leading to too focused source estimates. In this context, another idea consists in transforming the sources into a domain where their distribution is sparser than in the original source space and imposing sparsity in the transformed domain [START_REF] Ding | Reconstructing cortical current density by exploring sparseness in the transform domain[END_REF]. S5) Separability in space and wavevector of the spacewavevector content For each superficial distributed source, one can assume that the space-wavevector matrix at each time point, which is obtained by computing a local spatial Fourier transform of the measurements, can be factorized into a function that depends on the space variable only and a function that depends on the wavevector variable only [START_REF] Becker | Multi-way space-time-wave-vector analysis for EEG source separation[END_REF]. The space and wavevector variables are thus said to be separable, allowing for the use of tensor decomposition methods.

S6) Gaussian joint probability density function with parameterized spatial covariance

This prior assumes the source signals to be random variables that follow a Gaussian distribution with a spatial covariance matrix that can be described by a linear combination of a certain number of basis covariance functions [START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF][START_REF] Wipf | Robust Bayesian estimation of the location, orientation, and time course of multiple correlated neural sources using MEG[END_REF]. This combination is characterized by so-called hyperparameters, which have to be identified in the source imaging process.

Hypotheses on the temporal distribution T1) Smoothness

Since the autocorrelation function of the sources of interest usually has a full width at half maximum of several samples, the temporal distribution of the sources can be assumed to be smooth [START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF].

T2) Sparsity in a transformed domain

This assumption implies that the source signals admit a sparse representation in a transformed domain, which can, for example, be achieved by a wavelet or Gabor transform applied to the temporal dimension of the data [START_REF] Gramfort | Time-frequency mixed-norm estimates: Sparse M/EEG imaging with non-stationary source activations[END_REF]. The transformed signals can then be modeled using a small number of basis functions or atoms, which are determined by the source imaging algorithm.

T3) Pseudo-periodicity with variations in amplitude If the recorded data comprises recurrent events like a repeated time pattern that can be associated to the sources of interest, such as interictal epileptic spikes, one can exploit the repetitions as an additional diversity [START_REF] Möcks | Decomposing event-related potentials: a new topographic components model[END_REF]. This permits us to employ tensorbased methods. T4) Separability in time and frequency of the timefrequency content Under certain conditions such as oscillatory signals, the time and frequency variables of data transformed into the time-frequency domain (for example, by applying a wavelet transform to the measurements) can be assumed to separate (cf. hypothesis S4)) [START_REF] Miwakeichi | Decomposing eeg data into space-time-frequency components using parallel factor analysis[END_REF]. This hypothesis gives rise to tensor-based approaches for source separation. T5) Non-zero higher-order marginal cumulants Regarding the measurements as realizations of a vector of random variables, this hypothesis is required when resorting to statistics of order higher than two, that offer a better performance and identifiability than approaches based on second order statistics [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF]. It is generally verified in practice, as the signals of interest usually do not follow a Gaussian distribution.

Hypotheses on the spatio-temporal distribution

ST) Synchronous dipoles Contrary to point sources, which can be modeled by a single dipole, distributed sources are composed of a certain number of grid dipoles. These dipoles can be assumed to transmit synchronous signals [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF].

ALGORITHMS

In this section, we classify seven representative algorithms based on the hypotheses that they exploit. Furthermore, we briefly describe the characteristics of each method. sLORETA sLORETA [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sLORETA): technical details[END_REF] belongs to the class of minimum norm estimates (MNE) [START_REF] Hämäläinen | Interpreting magnetic fields of the brain: Estimates of current distributions[END_REF][START_REF] Pascual-Marqui | Review of methods for solving the EEG inverse problem[END_REF], which exploit hypothesis S1) by minimizing the L 2 -norm regularized cost function:

L(s) = ||x -Gs|| 2 2 + λ||s|| 2 2 (2)
where x and s correspond to columns of X and S, respectively. The first term on the right hand side of ( 2) is referred to as the data fit term, while the second term is called the regularization term. The regularization parameter λ manages a balance between both terms. The particularity of sLORETA consists in providing a source estimate that is normalized with respect to the source covariance, leading to an unbiased solution in the case of a single dipole source. MCE The minimum current estimate (MCE) algorithm [START_REF] Uutela | Visualization of magnetoencephalographic data using minimum current estimates[END_REF] has been developed to overcome the problem of blurred source reconstructions encountered with MNE and is based on hypothesis S3). In order to provide focused source estimates, it replaces the regularization term in (2) by the sparsity-inducing L 1 -norm, ||s|| 1 . This leads to a convex optimization problem, for which identifiability conditions and efficient algorithms are available.

VB-SCCD The VB-SCCD algorithm [START_REF] Ding | Reconstructing cortical current density by exploring sparseness in the transform domain[END_REF] imposes sparsity on the variational map, which characterizes variations in amplitude from one dipole source to adjacent dipole sources. This leads to a piecewise constant spatial distribution of the sources. To apply hypothesis S4), a regularization term of the form ||Ts|| 1 is used, where T is a transformation matrix which implements the variation operator. MxNE The mixed norm estimate (MxNE) [START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF] combines the hypotheses S3) and T1), leading to a spatially sparse, but temporally smooth source distribution. To achieve this, a mixed L 1,2 -norm regularization term is employed, leading to the cost function:

L(S) = ||X -GS|| 2 F + λ||S|| 1,2 . (3)
STWV-DA The space-time-wavevector based disk algorithm (STWV-DA) [START_REF] Becker | EEG extended source localization: tensor-based vs. conventional methods[END_REF] belongs to the class of tensor-based approaches and proceeds in two steps: i) the separation of different distributed sources using tensor decomposition based on hypotheses ST) and S5) (alternatively, one could exploit hypothesis T3) or T4)) and ii) the identification of grid dipoles characterizing each distributed source using hypothesis S4).

To construct a third order tensor, a local spatial Fourier transform is applied to the measurements, leading to space-timewavevector data. The tensor is then decomposed using the canonical polyadic decomposition [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF], which permits to identify a vector of spatial characteristics for each distributed source. Subsequently, to identify the dipoles belonging to each distributed source, this vector is compared to the spatial characteristics of elementary (circular) distributed sources whose union should fit the expected distributed source.

Champagne The Champagne algorithm [START_REF] Wipf | Robust Bayesian estimation of the location, orientation, and time course of multiple correlated neural sources using MEG[END_REF] is an empirical Bayesian approach and assumes that the sources at each time point are described by a zero-mean Gaussian random vector of covariance C s corresponding to hypothesis S6). The main step of the algorithm consists in estimating the covariance C s from the data. This is achieved by an optimization process that requires an estimate of the noise covariance matrix. Once the matrix C s is known, it permits to directly derive an estimate of the signals. In the absence of noise and for a source space with sufficiently high resolution, Champagne has been shown to provide a perfect source reconstruction. 4-ExSo-MUSIC The 4-ExSo-MUSIC algorithm [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF] is a representative of the class of subspace-based methods and exploits the 4-th order cumulants of the measurements, thus relying on assumption T5). Exploiting the orthogonality between the 4-th order signal and noise subspaces of the quadricovariance matrix, 4-ExSo-MUSIC identifies the elements of a dictionary of potential distributed sources (assumption ST)) that best describe the measurements based on a metric. The dictionary is constructed in the same way as for STWV-DA, thus also using hypothesis S4). The estimated sources are obtained as a union of a certain number of dictionary elements, determined after thresholding the 4-ExSo-MUSIC metric.

SIMULATIONS

To compare the performances of the different source imaging algorithms, we apply them to realistic simulated data in the context of epileptic EEG activity. We consider a realistic head model with three compartments representing the brain, the skull, and the scalp. The source space is composed of 19626 dipoles located on the cortical surface. The lead field matrix describing the propagation to 91 electrodes is computed numerically using a boundary element method. The epileptic source regions are modeled by patches which consist of adjacent grid dipoles. Highly correlated, physiologically realistic epileptiform spike-like signals are generated for all dipoles belonging to a patch using a coupled neuronal population model [START_REF] Cosandier-Rimélé | A physiologically plausible spatiotemporal model for EEG signals recorded with intracerebral electrodes in human partial epilepsy[END_REF]. To simulate epileptic activity spreading from one brain region to another, for scenarios with two patches, we use the same activities except for a time delay depending on the distance between the patches. All dipoles that do not belong to a patch are considered as generators of background activity. The SNR of epileptic and background activity is adjusted to realistic values. After spatially prewhitening the data based on the noise covariance matrix, we apply the sLORETA, MCE, MxNE, VB-SCCD, STWV-DA, Champagne, and 4-ExSo-MUSIC source imaging algorithms to the data of 10 epileptic spikes, which are concatenated for 4-ExSo-MUSIC and averaged for the other algorithms. To evaluate the source localization results, we employ the distance of localization error (DLE) [START_REF] Cho | Evaluation of algorithms for intracranial EEG (iEEG) source imaging of extended sources: feasibility of using (iEEG) source imaging for localizing epileptogenic zones in secondary generalized epilepsy[END_REF], which is defined as follows:

DLE = 1 2Q k∈I min ∈ Î ||r k -r ||+ 1 2 Q ∈ Î min k∈I ||r k -r ||. (4)
Here, I and Î denote the original and estimated sets of indices of all grid dipoles belonging to an active patch, Q and Q are the numbers of original and estimated active grid dipoles, and r k denotes the position of the k-th source dipole. The DLE is averaged over 50 realizations, obtained with different spikelike signals and varying background activity. First, we study the influence of the patch size on the source localization results. To this end, we vary the size of a patch located in the superior frontal gyrus (SupFr), considering 10, 100, and 400 grid dipoles, which corresponds to 0.5, 5, and 20 cm 2 of the cortical surface. Fig. 1 shows the original patches and a typical example of the reconstructed source configurations obtained by the different source imaging algorithms. sLORETA and Champagne yield similar results, which are blurred for the smallest patch and do not permit to accurately determine the spatial extent of the patches. MCE and MxNE provide very focal source estimates and are not suited to recover the size of the patch either. VB-SCCD overestimates the size of the smallest patch, recovering a much larger source region. However, for medium to large patches, it provides good estimates of the patches. STWV-DA and 4-ExSo-MUSIC yield similar source reconstructions that are very close to the ground truth for the small and medium patch and slightly worse for the large patch. These results are also reflected by the DLE values shown in Table 1. Overall, the best performances are achieved for the patch of medium size.

Second, we examine three scenarios with two patches, 1. While most methods show a tendency to yield somewhat improved source estimates for larger patch distance, the DLE does not decrease systematically with increasing distance for all algorithms and the source imaging performance seems to depend more on the patch position and form than the distance between two active patches. However, these results show that there is a gap between the performance of sLORETA, Champagne, MCE, and MxNE, which achieve a DLE between 2.5 and 3.7 cm for the three scenarios, and the DLE of VB-SCCD, STWV-DA, and 4-ExSo-MUSIC, that is smaller than 1.3 cm. This is due to the fact that VB-SCCD, STWV-DA, and 4-ExSo-MUSIC permit to identify distributed sources, contrary to sLORETA, Champagne, MCE, and MxNE, as has also become apparent in the previous simulation.

CONCLUSIONS

In this paper, we have proposed a classification of different source imaging algorithms based on the hypotheses that are imposed on the temporal and spatial distributions of the sources. To illustrate the use of various hypotheses, we have provided a short presentation of seven representative algorithms, whose performance has been compared in a simulation study. This study has demonstrated the superior performance of VB-SCCD, STWV-DA, and 4-ExSo-MUSIC for the recovery of distributed sources as compared to sLORETA, Champagne, MCE, and MxNE. Among the latter four methods, we have observed that sLORETA generally leads to the best performance. Moreover, we noticed that most algorithms (except for STWV-DA and 4-ExSo-MUSIC) require spatial prewhitening or an estimate of the noise covariance in order to yield accurate results. Therefore, the best overall performance, both in terms of robustness and accuracy, is achieved by STWV-DA. Further studies should be conducted to confirm these results.

To further enhance the performance, one could consider integrating the steps of two-step procedures such as STWV-DA into one single step in order to process all the available information and constraints at the same time. Another promising perspective for future research consists in exploring different combinations of a priori information, for example, by merging successful strategies of different recently established source imaging approaches, such as tensor-based, subspacebased, or Bayesian approaches and sparsity.

Fig. 1 .

 1 Fig. 1. Original patch and source configurations estimated with different source imaging algorithms for the patch SupFr composed of 10, 100, and 400 dipoles.

Table 1 .

 1 DLE (in cm) of source imaging algorithms for different scenarios scenario sLORETA Champagne MCE MxNE VB-SCCD STWV-DA 4-ExSo-MUSIC

	SupFr, small size	3.62	5.15	4.80	4.94	3.30	0.41	0.65
	SupFr, medium size	1.97	2.51	2.95	3.03	0.53	0.43	0.44
	SupFr, large size	5.45	3.38	4.02	4.78	0.89	1.14	1.15
	InfFr+SupFr, small dist.	2.55	3.21	3.71	3.58	1.16	0.73	0.85
	InfFr+InfPa, medium dist.	2.97	3.31	3.51	3.52	1.23	0.59	0.61
	SupFr+SupOcc, large dist.	2.49	2.63	2.95	3.11	0.44	0.77	0.63
	composed of 100 dipoles each and all located on the left hemi-
	sphere: patches InfFr and SupFr (inferior frontal and superior
	frontal gyri) with small distance, patches InfFr and InfPa (in-
	ferior frontal and inferior parietal gyri) with medium distance,
	and patches SupFr and SupOcc (superior frontal and superior
	occipital gyri) with large distance. Due to limited space, we
	only provide the results in terms of DLE, summarized in Table


Julian Quiroga
email: quiroga.j@javeriana.edu.co

Thomas Brox
email: brox@cs.uni-freiburg.de

Frédéric Devernay

James Crowley

Dense Semi-Rigid Scene Flow Estimation from RGBD images ⋆

Keywords: motion, scene flow, RGBD image

Scene flow is defined as the motion field in 3D space, and can be computed from a single view when using an RGBD sensor. We propose a new scene flow approach that exploits the local and piecewise rigidity of real world scenes. By modeling the motion as a field of twists, our method encourages piecewise smooth solutions of rigid body motions. We give a general formulation to solve for local and global rigid motions by jointly using intensity and depth data. In order to deal efficiently with a moving camera, we model the motion as a rigid component plus a non-rigid residual and propose an alternating solver. The evaluation demonstrates that the proposed method achieves the best results in the most commonly used scene flow benchmark. Through additional experiments we indicate the general applicability of our approach in a variety of different scenarios.

Introduction

The 3D motion field of a scene is useful for several computer vision applications, such as action recognition, interaction, or 3D modeling on nonrigid objects. One group of methods uses a stereo or multi-view camera system, where both scene flow and depth are estimated from the images, while a second group uses RGBD images as input. In the latter case, the depth given by the sensor may be used directly for scene flow estimation. In this paper, we present a new formulation for this second group.

The two main questions for scene flow estimation from RGBD images are: (a) how to fully exploit both sources of data, and (b) which motion model should be used to compute a confident scene flow. When both intensity and depth are provided, there is not a consensus on how to combine these two sources of information most effectively. A straightforward approach would compute the optical flow from RGB images and infer the scene flow from the depth data. It is also possible to generate colored 3D point clouds and try to find out the 3D motion vectors consistent with the input data. However, it is not mandatory to explicitly represents points in 3D to solve for the scene flow: the scene structure can be represented in the image domain, where color and depth data can be coupled using a projective function. This way, depth changes influence the motion in the image domain and consistency constraints can be formulated jointly over the color and depth images. This is the approach we follow in this work. It gives us access to some of the powerful tools developed for 2D motion estimation.

The second question has been addressed less. Scene flow estimation using intensity and depth is an ill-posed problem and regularization is needed. Usually, the 3D motion vector of each point is solved to minimize color and intensity constraints in a data term and the whole 3D motion field is regularized to get spatially smooth solutions while preserving discontinuities. Since depth data is available, a weighted regularization can be used to preserve motion discontinuities along depth edges, where independent motions are probable to appear. However, solving for a piecewise smooth solution of the 3D motion field may not be the best choice for some motions of interest. For example, a 3D rotation of a rigid surface induces a great variety of 3D motions that are hardly well regularized by such an approach. A similar issue occurs when the RGBD sensor is moving with respect to the scene.

In this work, we take advantage of the fact that most real scenes can be well modeled as locally or piecewise rigid, i.e., the scene is composed of 3D independently rigid components. The main contribution of this paper is the definition of an over-parametrized framework for scene flow estimation from RGBD images. We model the scene flow as a vector field of rigid body motions. This representation helps the regularization process which, instead of directly penalizing variations in the 3D motion field, encourages piecewise smooth solutions of rigid motions. Moreover, the proposed rigid body approach can be constrained in the image domain to fully exploit intensity and depth data. This formulation is flexible enough to support different data constraints and regularization strategies, and it can be adapted to more specialized problems. By using the same general framework, it is possible to model the scene flow as a global rigid motion plus a non-rigid residual, which is particularly useful when estimating the motion of deformable objects in conjunction with a moving camera.

Related work

Scene flow was first introduced by Vedula [START_REF] Vedula | Three-dimensional scene flow[END_REF] as the 3D motion field of the scene. Since this seminal work several approaches have been proposed to compute the 3D motion field. If a stereo or multi-view camera system is available, the scene flow can be computed by enforcing consistency with the observed optical flows [START_REF] Vedula | Three-dimensional scene flow[END_REF], by a decoupled [START_REF] Wedel | Efficient dense scene flow from sparse of dense stereo data[END_REF] or joint [START_REF] Basha | Multi-view scene flow estimation: A view centered variational approach[END_REF][START_REF] Huguet | A variational method for scene flow estimation from stereo sequences[END_REF] estimation of structure and motion, or by assuming a local rigidity of the scene [START_REF] Vogel | 3D scene flow estimation with a rigid motion prior[END_REF]. In this work we assume that a depth sensor is available and the estimation of the structure is not needed. The first work using intensity and depth was by Spies et al. [START_REF] Spies | Dense range flow from depth and intensity data[END_REF], where the optical flow formulation by Horn and Schunck [START_REF] Horn | Determining optical flow[END_REF] is extended to include depth data. In this approach, depth data is used simply as an additional channel in the variational formulation of the optical flow. The range flow is estimated using the observed data and enforced to be smooth. However, the scene is assumed to be captured by an orthographic camera and there is no coupling between optical and range flows. The coupling issue can be solved as in [START_REF] Letouzey | Scene flow from depth and color images[END_REF], where the scene flow is directly solved using the depth data to constrain the 3D motion in the image domain. However in that work, the depth data is not fully exploited since there is no range flow constraint. Similar to [START_REF] Spies | Dense range flow from depth and intensity data[END_REF], this method suffers from the early linearization of the constancy constraints and from over-smoothing along motion boundaries because of the L 2 -regularization. Quiroga et al. [START_REF] Quiroga | Scene flow by tracking in intensity and depth data[END_REF] define a 2D warping function to couple image motion and 3D motion, allowing for a joint local constraint of the scene flow on intensity and depth data. Although the method is able to deal with large displacements, it fails on untextured regions and more complex motions, such as rotations.

In order to solve for dense scene flow, a regularization procedure is required. Usually, the 3D motion field is assumed to be piecewise smooth, and total variation (TV) is used as regularizer. The work by Herbst [START_REF] Herbst | RGB-D flow: Dense 3-D motion estimation using color and depth[END_REF] follows this idea, but as [START_REF] Spies | Dense range flow from depth and intensity data[END_REF], it lacks a coupling between optical and range flows, and the regularization is done on the optical flow rather than on the scene flow. In [START_REF] Quiroga | Local/global scene flow estimation[END_REF], a variational extension of [START_REF] Quiroga | Scene flow by tracking in intensity and depth data[END_REF] is presented. A weighted TV is applied on each component of the 3D motion field, aiming to preserve motion discontinuities along depth edges.

All these methods assume spatial smoothness of the scene flow, which is a reasonable assumption for translational motions but not for rotations. Under a rotation, even close scene points present different 3D motions. In case of a moving camera the regularization of the motion field can be a challenge. In this work, we use an over-parametrization of the scene flow, where each scene point is allowed to follow a rigid body motion. This way, the regularization can be done on a field of rigid motions, favoring piecewise solutions, which is a better choice for real scenes. A similar idea is presented in [START_REF] Rosman | Group-valued regularization for analysis of articulated motion[END_REF], where a regularization of a field of rigid motions is proposed. Our approach differs from that work in three ways. First, we use a more compact representation of the rigid-body motion via the 6-parameter twist representation, instead of a R 12 embedding. Second, our approach solves and regularizes the rigid motion field at the same time. Finally, we decouple the regularization of the rotational and translational fields, which simplifies the optimization and allows the use of different TV strategies on each field. Similar to [START_REF] Quiroga | Scene flow by tracking in intensity and depth data[END_REF], we use a depth-based weighting function to avoid penalization along surface discontinuities.

In an alternative approach, Hadfield and Bowden [START_REF] Hadfield | Scene particles: Unregularized particle-based scene flow estimation[END_REF] estimate the scene flow using particle filtering in a 3D colored point cloud representation of the scene. In that approach, a large set of motion hypotheses must be generated and tested for each 3D point, leading to high computational costs. We rather use the 2D parametrization provided by RGBD sensors to formulate an efficient 3D motion exploration. As done in [START_REF] Quiroga | Scene flow by tracking in intensity and depth data[END_REF], we define a warping function to couple the twist motion and the optical flow. A very similar warp is presented in [START_REF] Kerl | Robust odometry estimation for RGB-D cameras[END_REF] to solve for a global rigid motion from RGBD images. In our approach, we use the warping function to locally constrain the rigid motion field in the image domain. The global rigid motion estimation can be seen as a particular case of the proposed formulation. Moreover, unlike [START_REF] Kerl | Robust odometry estimation for RGB-D cameras[END_REF], we define a depth consistency constraint to fully exploit both sources of data. Thus, we can solve for the local twist motion that best explains the observed intensity and depth data, gaining robustness under noise. The local solver, in conjunction with the TV regularization of the twist field, provides an adjustable combination between local and piecewise rigidity.

As we present in the experiments, this framework is flexible enough to estimate a global rigid body motion or to solve for a general 3D motion field in challenging setups. This way, we are able to model the motion of the scene as a global rigid motion and a non-rigid residual.

Scene motion model

We parameterize every visible 3D point into the image domain

Ω ⊂ R 2 . The projection π : R 3 → R 2 maps a 3D point X = (X, Y, Z) onto a pixel x = (x, y) on Ω by: π(X) = f x X Z + c x , f y Y Z + c y T , (1)
where f x and f y are the focal lengths of the camera and (c x , c y ) its principal point. The inverse projection π -1 : R 2 × R → R 3 back-projects an image point to the 3D space for a given depth z, as follows:

π -1 (x, z) = z x -c x f x , z y -c y f y , z T . (2)
We will use X or π -1 (x, z) as the 3D representation of image point x.

The scene flow v(x) : Ω → R 3 is defined as the 3D motion field describing the motion of every visible 3D point between two time steps accordingly

X t+1 = X t + v (x t ) .
(3)

Twist motion field

Instead of directly representing the elements of the motion field as 3D vectors, we use an over-parametrized model to describe the motion of every point as a rigid transformation. The group action of a rigid body transformation can be written as T (X) = RX + t, where t ∈ R 3 is a translation, and R ∈ SO(3) is a rotation matrix. Using homogeneous coordinates, a 3D point X = (X, 1) T is transformed into X′ accordingly to

X′ = G X, with G = R t 0 1×3 1 ∈ SE(3). ( 4
)
Since G has only 6 degrees of freedom, a more convenient and compact representation is the 6-parameter twist. Every rigid motion can be described as a rotation around a 3D axis ω = (ω X , ω Y , ω Z ) T and a translation τ = (τ x , τ y , τ z ) T along this axis. Therefore it can be shown that for any arbitrary G ∈ SE(3) there exists an equivalent ξ ∈ R 6 twist representation. A twist ξ = (τ, ω) can be converted into the G representation with the following exponential map:

G = e ξ = I + ξ + ( ξ) 2 2! + + ( ξ) 3 3! + • • • , (5)
where

ξ = ω τ 0 1×3 1 , with ω =   0 -ω Z ω Y ω Z 0 -ω X -ω Y ω X 0   . (6)
Correspondingly, for each G ∈ SE(3) there exist a twist representation given by the logarithmic map via ξ = log(G); see [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF] for more details.The motion of the scene is embedded in a twist motion field ξ(x) : Ω → R 6 , where the motion of every 3D point between two time steps is given by:

Xt+1 = e ξ(xt) Xt . (7)

Twist motion on the image

For every time t the scene is registered as an RGBD image S t (x) = {I t (x), Z t (x)} where images I(x) and Z(x) provide color and depth, respectively, for every point x ∈ Ω. Under the action of a twist ξ between t and t + 1, the motion of a 3D point induces an image flow (u, v). We define the warping function W(x, ξ) : R 2 ×R 6 → R 2 that maps each non-occluded pixel onto its new location after the rigid motion. The warping function is defined as:

W(x, ξ) = π e ξ X , with X = π -1 (x, Z(x)) 1 . (8)
Let S k (x) be a component of the RGBD image (e.g. brightness or depth) evaluated at x ∈ Ω and let ρ S k (x, ξ) be a robust error function between S k t (x) and S k t+1 (W(x, ξ)). Every twist motion vector ξ is solved to minimize one or several consistency functions ρ S k (x, ξ). Optimization problems on manifolds such as SE(3) can be solved by calculating incremental steps in the tangent space to the manifold. Considering that an initial estimate ξ is known, the goal at each optimization step is to find the increment ∆ξ which (approximately) minimizes ρ S k (x, ξ). Since SE(3) is a Lie group (not an Euclidean space) with the composition as operation, the exponential and logarithmic functions are used to update the current estimate ξ with the new increment ∆ξ according to ξ ← log(e ∆ξ e ξ ) [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF]. The iterative solution requires a linearization of the warping function. Given the initial estimate ξ, the image point x becomes x ξ = W(x, ξ) and the warping function satisfies:

W(x, log(e ∆ξ e ξ )) = W(x ξ , ∆ξ) = x ξ + δx(x ξ , ∆ξ), (9)
where δx(x ξ , ∆ξ) is the image flow induced by the increment ∆ξ, which is given by δx(x ξ , ∆ξ) = π(e ∆ξ Xξ )π( Xξ ), with Xξ = π -1 (x ξ , Z(x ξ )), 1 T . Assuming a small rotation increment, the exponential function is approximated as e ∆ξ ≈ I + ∆ξ. Therefore for a small increment ∆ξ the warping function can be well approximated by the following linear version:

W(x, log(e ∆ξ e ξ )) = W(x, ξ) + J(x ξ )∆ξ, (10)
where J(x ξ ) is the Jacobian matrix, given by:

J(x ξ ) =    f x Z(x ξ ) 0 - x ξ Z(x ξ ) - x ξ y ξ fy fx+x 2 ξ fy - y ξ fx fy 0 f y Z(x ξ ) - y ξ Z(x ξ ) - fy+y 2 ξ fx x ξ y ξ fx x ξ fy fx    . ( 11
)

Scene flow formulation

Given two pairs of RGBD images {I 1 , Z 1 } and {I 2 , Z 2 } the goal is to solve for the scene flow field that best explains the observed data. Due to vanishing gradients, the aperture problem and outliers, this motion computation is an illposed problem that cannot be solved independently for each point. Therefore, a smoothness in the motion field must be assumed. In our formulation, we consider only spatial smoothness, but it can be extended to temporal smoothness. In general, we solve for the twist field ξ minimizing the following energy:

E(ξ) = E D (ξ) + αE S (ξ), (12)
where the data term E D (ξ) measures how consistent is the estimated twist-based model with the observed color and depth data, and the smoothness term favors piecewise smooth fields while preserving discontinuities.

Data term

The warping function [START_REF] Huguet | A variational method for scene flow estimation from stereo sequences[END_REF] enables the formulation of consistency constraints for the twist field in the image domain. Using the gray value image we define a brightness constancy assumption:

I 2 (W(x, ξ)) = I 1 (x), (13)
and a gradient constancy assumption:

I g 2 (W(x, ξ)) = I g 1 (x). ( 14
)
where I g (x) = |∇I(x)|. The gradient assumption reduces the effect of brightness changes. We use the gradient magnitude since it is invariant to rotation. The estimated twist motion should satisfy these constraints for most points.

On the other hand, the surface changes under the twist action. This variation must be consistent with the observed depth data. Therefore we define a depth variation constraint given by:

Z 2 (W(x, ξ)) = Z 1 (x) + δ Z (x, ξ), (15)
where δ Z (x, ξ) is the depth variation induced on the 3D point π -1 (x, Z 1 (x)) by the twist ξ, obtained from the third component of the 3D vector:

δ 3D (x, ξ) = e ξ -I 4×4 X with X = π -1 (x, Z 1 (x)) 1 . ( 16
)
This equation enforces the consistency between the motion captured by the depth sensor and the estimated motion. Without regularization, equations ( 13), ( 14) and ( 15) alone are not sufficient to constrain the twist motion for a given point since there is an infinite number of twists that satisfy these constraints. However, real scenes are locally rigid and it is possible to solve for a local twist explaining the observed motion. We formulate the twist motion estimation as a local least-square problem by writing the data term as follows:

E D (ξ) = x x ′ ∈Nx Ψ ρ 2 I (x ′ , ξ (x)) + γρ 2 g (x ′ , ξ (x)) + λΨ ρ 2 Z (x ′ , ξ (x)) , (17)
where N x is an image neighborhood centered on x, and with the brightness, gradient and depth residuals, respectively given by:

ρ I (x, ξ) = I 2 (W(x, ξ)) -I 1 (x) , (18)
ρ g (x, ξ) = I g 2 (W(x, ξ)) -I g 1 (x) , (19)
ρ Z (x, ξ) = Z 2 (W(x, ξ)) -(Z 1 (x) + δ Z (x, ξ)). (20)
Constant γ balances brightness and gradient residuals, while λ weights intensity and depth terms. We use the robust norm Ψ s 2 = √ s 2 + ε 2 , which is a differentiable approximation of the L 1 norm, to cope with outliers due to occlusion and non-rigid motion components in N x .

Smoothness term

A twist motion field ξ can be decomposed into a rotational field ω, with e ω : Ω → SO(3), and a 3D motion field τ : Ω → R 3 . Since they are decorrelated by nature, we regularize each field independently using weighted Total Variation (TV), which allows piecewise smooth solutions while preserving motion discontinuities.

Regularization of the fields ω and τ poses different challenges. Elements of τ lie in the Euclidean space R 3 and the problem corresponds to a vector-valued function regularization. Different TV approaches can be used to regularize τ , as described in [START_REF] Goldluecke | The natural vectorial total variation which arises from geometric measure theory[END_REF]. Particularly the channel-by-channel L 1 norm has been successfully used for optical flow [START_REF] Wedel | An improved algorithm for TV-L1 optical flow[END_REF]. We define the weighted TV of τ as:

TV c (τ ) = x c(x) ∇τ (x) , (21)
where ∇τ

:= |∇τ X | + |∇τ Y | + |∇τ Z | and c(x) = e -β|∇Z1(x)| 2 .
The weighting function c helps preserve motion discontinuities along edges of the 3D surface. Moreover, the L 1 norm can be replaced by the Huber norm [START_REF] Werlberger | Anisotropic Huber-L1 Optical Flow[END_REF] to reduce the staircasing effect. Efficient solvers for both norms are presented in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

Elements of ω are rotations in the Lie group SO(3) embedded in R 3 through the exponential map, and the regularization has to be done on this manifold. In order to apply a TV regularization, a notion of variation should be used for elements of SO(3). Given two points e ω1 , e ω2 ∈ SO(3) the residual rotation can be defined as e -ω1 e ω2 . The product in logarithmic coordinates can be expressed as e ω1 e ω1 = e µ(ω1,ω2) , where the mapping µ can be expanded in a Taylor series around the identity, using the Baker-Campbell-Hausdorff formula:

µ(ω 1 , ω2 ) = ω1 + ω2 + 1 2 [ω 1 , ω2 ] + O(|(ω 1 , ω2 )| 3 ), (22)
where [•, •] is the Lie bracket in so(3). Close to the identity, ( 22) is well approximated by its first-order terms, so that for small rotations the variation measure can be defined as the matrix subtraction in so(3), or equivalently, as a vector difference for the embedding in R 3 . Accordingly, the derivative matrix Dω := (∇ω X , ∇ω Z , ∇ω Z ) T : Ω → R 3×2 approximates the horizontal and vertical point-wise variations of ω on the image. Following [START_REF] Goldluecke | The natural vectorial total variation which arises from geometric measure theory[END_REF], we define the TV as the sum over the largest singular value σ 1 of the derivative matrix:

TV σ (ω) = x c(x) σ 1 (Dω(x)) . (23)
This TV approach supports a common edge direction for three components, which is a desirable properties for the regularization of the field of rotations. Moreover, deviations are less penalized with respect to other measures (e.g. Frobenius norm [START_REF] Rosman | Group-valued regularization for analysis of articulated motion[END_REF]) and efficient solvers are available. However, this TV definition approximates the real structure of the manifold yielding to a biased measures far from the identity. The more the rotation is away from the identity, the more its variations in SO(3) are penalized as is shown hereinafter. Given two rotations

ω 1 = θ 1 ← - ω 1 and ω 2 = θ 2 ← - ω 2
, with θ the angle and ←ω the unitary axis vector of the rotation, and writing

θ 2 = θ 1 + δθ, leads to ω 2 -ω 1 = θ 1 ( ← - ω 2 -← - ω 1 ) -δθ ← - ω 2
. This linearly dependent penalization usually is not a problem, since large rotations imply larger motion on the image. Thus, a stronger regularization can be reasonable. Moreover, large rotation caused by a global motion of the scene or the camera can be optimized separately and compensated, as we show in Sec. 4.4. Optionally, this over-penalization can be removed by expressing each rotation as ω = θ ←ω and applying a vectorial TV on ←ω and a scalar TV on θ. In our approach, the full smoothing term is given by:

E S (ξ) = TV c (τ ) + TV σ (ω).
(24)

Optimization

The proposed energy ( 12) is minimized by decomposing the optimization into two simpler problems. We use the variable splitting method [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF] with auxiliary variable χ, and the minimization problem becomes:

min ξ,χ E D (ξ) + 1 2κ x |ξ(x) -χ(x)| 2 + αE S (χ), ( 25
)
where κ is a small numerical variable. Note that the linking term between ξ and χ is the distance on the tangent space at the identity in SE(3). The solution of (25) converges to that of ( 12) as κ → 0. Minimization of this energy is performed by alternating the two following optimization problems:

i. For fixed χ, estimate ξ that minimizes (25). This optimization problem can be solved point-wise by minimizing:

1 2κ |ξ -χ| 2 + x ′ ∈Nx Ψ ρ 2 I + γρ 2 g + λΨ ρ 2 Z , (26)
where the parameters (x ′ , ξ) are considered implicit. This energy can be linearized around an initial estimate ξ, using a first-order Taylor series expansion:

1 2κ log(e ∆ ξ e ξ ) -χ 2 + x ′ ∈Nx Ψ |ρ I + I x J∆ξ| 2 + γ |ρ G + I g x J∆ξ| 2 + λΨ |ρ Z + (Z x J -K) ∆ξ| 2 . ( 27
)
where I x , I g x and Z x are the row vector gradients of I 2 (x), I g 2 (x) and Z 2 (x), respectively, and J is the Jacobian [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF] of the warp, all evaluated at x ξ = W(x, ξ). The 1 × 6 vector K is defined as K = D ([X ξ ] × |I 3×3 ) with D = (0, 0, 1) isolating the third component and [•] × the cross product matrix. Finding the minimum of (27) requires an iterative approach. Taking the partial derivative with respect to ∆ξ and setting it to zero, the increment ∆ξ can be computed as:

∆ξ = -H -1 x ′ ∈Nx Ψ ′ ρ 2 I + γρ 2 G (I x J) T ρ I + γ(I g x J) T ρ G +λ Ψ ′ ρ 2 Z (Z x J -K) T ρ Z + 1 κ log (e ξ ) -1 e χ , ( 28
)
where Ψ ′ is the derivative of the robust norm, which is evaluated at the current estimate ξ. The 6 × 6 matrix H is the Gauss-Newton approximation of the Hessian matrix and is given by:

H = x ′ ∈Nx Ψ ′ ρ 2 I + γρ 2 G (I x J) T (I x J) + γ (I g x J) T (I g x J) + λ Ψ ′ ρ 2 Z (Z x J -K) T (Z x J -K) + 1 2κ I 6×6 . ( 29
)
ii. For fixed ξ = (ω, τ ), compute χ = (̟, π) that minimizes:

TV c (π) + η 2 x |π(x) -τ (x)| 2 + TV σ (̟) + η 2 x |̟(x) -ω(x)| 2 (30)
where η = (κα) -1 . Each side of equation ( 30) corresponds to a vectorial image denoising problem with a TV-L 2 model (ROF model). Efficient primal-dual algorithms exist for solving both problems. The left side is solved component-wise using the first-order primal-dual algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. For the right side we use the vectorial approach [START_REF] Goldluecke | The natural vectorial total variation which arises from geometric measure theory[END_REF] which allows an optimal coupling between components. A small modification is necessary in both cases to include the weighting function.

Scene flow with camera motion estimation

In many applications, the sensor itself moves relative to the observed scene and causes a dominant global motion in the overall motion field. In this situation, compensating for the motion of the camera can simplify the estimation and regularization of the scene flow. Moreover, for 3D reconstruction of deformable objects, the camera motion is needed to register partial 3D reconstructions. Therefore, we consider splitting the motion of the scene into a globally rigid component ξ R = (τ R , ω R ) ∈ R 6 , capturing the camera motion relative to the dominant object/background, plus a non-rigid residual field ξ = (τ, ω). We assume that a large part of the scene follows the same rigid motion. Accordingly, the scene flow is defined as the composition χ = log(e ξ + e ξR -I 4×4 ), and the estimation problem is formulated as:

min χ E Rig (χ) + E Res (χ) , (31)
with E Rig (χ) and E Res (χ) the rigid and non-rigid energies, respectively. It is worth noting that the separation of the camera motion is in addition to the framework presented above, i.e., the non-rigid part can still deal with local motion.

Rigid energy. The camera motion can be estimated using the data term [START_REF] Sturm | A benchmark for the evaluation of RGB-D slam systems[END_REF], by considering every pixel (or a subset of Ω) to solve for a unique twist ξ R . Accordingly, the rigid component of the energy is defined as:

E Rig (χ) = x Ψ ρ 2 I (x, χ) + γρ 2 g (x, χ) + λΨ ρ 2 Z (x, χ) . (32)
Non-rigid residual energy. The residual motion can be computed following [START_REF] Quiroga | Scene flow by tracking in intensity and depth data[END_REF], with the non-rigid energy given by:

E Res (χ) = E D (χ) + αE S (χ). (33)
We minimize (31) by an iterative, alternating estimation of ξ R and ξ:

a. Given a fixed ξ, solve for ξ R that minimizes E Rig (χ). This is done by iteratively applying (28) and ( 29), with a zero auxiliary flow.

b. Given a fixed ξ R , solve for ξ that minimizes E Res (χ). This is done by iterating steps i and ii in Sec. 4.3.

Experiments

Implementation details

In order to compute the scene flow, the proposed method assumes that a pair of calibrated RGBD images is provided. Regardless of the depth sensor, depth data is always processed in cm and RGB color images are transformed to intensity images and normalized. For all the experiments α = 10, β = 1, γ = 0.1 and λ = 0.1. For each scene a depth range is defined and only pixels having a valid depth measure inside the range are taken into account for the data term and final measurements. However, all pixels are considered in the regularization.

We use a multi-scale strategy in order to deal with larger motions. We construct an image pyramid with a downsampling factor of 2. We apply a Gaussian anti-aliasing filter to the intensity image and the pyramid is built using bicubic downsampling. For the depth image, a 5 × 5 median filter is used and the pyramid is constructed by averaging pixels in non-overlapped neighborhoods of 2 × 2, where only pixels with a valid depth measure are used. Having a pyramid with levels l = {0, 1, ..., L}, with 0 the original resolution, the computation is started at level L an the estimated twist field is directly propagated to the next lower level. The camera matrix is scaled at each level by the factor 2 l . The neighborhood N x is defined as a N × N centering window. At each level we perform M loops consisting of M GN = 5 iterations of the Gauss-Newton procedure followed by M TV = 50 iterations of the TV solver. The constant κ is styled at each scale.

Middlebury datasets

The Middlebury stereo dataset [START_REF] Scharstein | High-accuracy stereo depth maps using structured light[END_REF] is commonly used as benchmark to compare scene flow methods [START_REF] Basha | Multi-view scene flow estimation: A view centered variational approach[END_REF][START_REF] Huguet | A variational method for scene flow estimation from stereo sequences[END_REF][START_REF] Hadfield | Scene particles: Unregularized particle-based scene flow estimation[END_REF][START_REF] Quiroga | Local/global scene flow estimation[END_REF]. Using images of one of these datasets is equivalent to a fixed camera observing an object moving in X direction. As in [START_REF] Hadfield | Scene particles: Unregularized particle-based scene flow estimation[END_REF][START_REF] Quiroga | Local/global scene flow estimation[END_REF], we take images 2 and 6 as the first and second RGBD image, respectively, and use the ground truth disparity map of each image as depth channel. Stereobased methods [START_REF] Basha | Multi-view scene flow estimation: A view centered variational approach[END_REF][START_REF] Huguet | A variational method for scene flow estimation from stereo sequences[END_REF] do not assume RGBD images and simultaneously estimate the optical flow and disparity maps by considering images 2, 4, 6 and 8 of each dataset. The ground truth for the scene corresponds to the camera motion along the X axis, while the optical flow is given by the disparity map. The scene flow error is measured in the image domain using root mean squared error (RMS) and the average angle error (AAE) of the optical flow. In order to compare with optical flow methods, we include results for the scene flow inferred using LDOF [START_REF] Brox | Large displacement optical flow: Descriptor matching in variational motion estimation[END_REF] and the depth data, as is described in [START_REF] Hadfield | Scene particles: Unregularized particle-based scene flow estimation[END_REF]. Results and comparison for Teddy Table 1. Middlebury dataset: errors on the optical flow extracted from the scene flow (except for [START_REF] Brox | Large displacement optical flow: Descriptor matching in variational motion estimation[END_REF], which is an optical flow method). See Sec. 5.2 for details. 1, where stereo methods are denoted with 2 views. In this experiment, we used a 5-level pyramid, with M = 5, N = 3 and for each level l, we set κ = 10 4 10 -l . A non-optimized implementation of our method processes each dataset in about 60 The proposed approach outperforms previous methods. Because the 3D motion field resulting from camera translation is constant, Middlebury datasets are not well suited to fully evaluate the performance of scene flow methods.

Teddy

Scene flow from RGBD data

We performed further experiments on more complex scenes using two RGBD sensors: the Microsoft Kinect for Xbox and the Asus Xtion Pro Live. We consider three different setups: i) a fixed camera, ii) a moving camera observing a rigid scene and iii) a moving camera capturing deformable objects. In each case, we show the input images, the optical flow, and one or more components of the scene flow. To give an idea of the motion, the average of the two input RGB images is shown. The optical flow is visualized using the Middlebury color code [START_REF] Scharstein | High-accuracy stereo depth maps using structured light[END_REF]. For the scene flow, we show each component using a cold-to-warm code, where green color is zero motion, and warmer and colder colors represent positive and negative velocities, respectively. In the first experiment, the Kinect sensor is fixed to compute the scene flow from two sequences (Fig. 1). The top row shows the deformation of a poster, which produces a non-uniform deformation. The proposed method is able to capture the poster deformation, thus it is possible to accurately estimate the changes in depth when the poster is folded. The gradient constancy constraint plays an important role here, since the sensor applies automatic white balancing. The bottom images show a motion performed with arms and hands. While hands are rotating inwards, the elbows lift, and a region of both arms remains almost still. This composite motion generates a discontinuous optical flow, which is well estimated by our method. Moreover, it can be seen that small rotations and articulated motions are well described for the proposed motion model.

The second experiment considers a static scene observed by a moving camera. We use Teddy images of the RGBD dataset [START_REF] Sturm | A benchmark for the evaluation of RGB-D slam systems[END_REF]. Unlike the Middlebury dataset, [START_REF] Quiroga | Local/global scene flow estimation[END_REF]. Results of the proposed method are clearly more accurate than those of [START_REF] Quiroga | Local/global scene flow estimation[END_REF].

this scene presents a changing 3D motion field due to the translation and rotation of the sensor. We first estimate a single rigid motion as is presented in Sec. 4.4. We estimate the parameters using every second pixel with a range between 50 cm and 150 cm from the sensor. At each level of the pyramid we run 100 iterations. This result is taken as baseline and compared in Figure 2 with resulting optical flows of our method and the dense estimation provided by [START_REF] Quiroga | Local/global scene flow estimation[END_REF]. The proposed method produces a close estimation of the camera motion without assuming a global rigidity of the scene. In contrast, the direct regularization on the scene flow components, as is performed in [START_REF] Quiroga | Local/global scene flow estimation[END_REF], fails to capture the diversity of 3D motions introduced by the camera rotation.

In the last experiment, a moving Asus Xtion camera observes a non-rigid scene. To estimate the rigid and non-rigid components we perform 3 rounds of alternation at each level of the pyramid. As is shown in Figure 3, the proposed approach allows the joint estimation of both components. Particularly, it is possible to compute the motion of both hands and the rotation of the face while the camera is turning. Capturing the motion of thin objects, as the fingers of the left hand, is a challenge since depth data is incomplete and very noisy in this area.

The size of the local neighborhood was kept fixed for all the scales and for every position on the image. Results could be improved by adjusting the size of the window using the depth data, in order to get a constant resolution.

Summary

We have presented a new method to compute dense scene flow from RGBD images by modeling the motion as a field of rigid motions. This allows for piecewise smooth solutions using TV regularization on the parametrization. We have decoupled the regularization procedure for the rotational and translational part and proposed some approximations to simplify the optimization. Future advances on manifold regularization may provide even more accurate and faster solvers that can be used with our parameterization. In order to fully exploit both intensity and depth data, we constrain the rigid body motion in the image domain. This way, we can solve for the local rigid motion as an iteratively reweighted least squares problem. The proposed approach provides an adjustable combination between local and piecewise rigidity, which, in conjunction with a global rigid estimation, is able to capture the motion in real world scenes.

Fig. 1 .

 1 Fig. 1. Scene flow estimation with a fixed camera. a) Left: Input images. b) Middle: Optical flow c) Right: Z-component of the scene flow.

Fig. 2 .

 2 Fig. 2. Camera motion. From left to the right: a) Input images. b) Optical flow of the rigid motion estimation. c) Optical flow using the proposed method. d) Optical flow by the scene flow method[START_REF] Quiroga | Local/global scene flow estimation[END_REF]. Results of the proposed method are clearly more accurate than those of[START_REF] Quiroga | Local/global scene flow estimation[END_REF].

Fig. 3 .

 3 Fig. 3. Motion estimation with a moving camera. Top, from left to right: a) Input images. b) Optical flow of the global rigid motion. c) Optical flow of the non-rigid residual. Bottom, from left to right: (X, Y, Z) components of the non-rigid 3D motion.

⋆ This work was supported by a collaborative research program between Inria Grenoble and University of Freiburg. We gratefully acknowledge partial funding by CMIRA 2013 (Region Rhône-Alpes), GDR ISIS (CNRS), and COLCIENCIAS (Colombia).


M K Smail

J Waeytens

F Bourquin

THE HEALTH MONITORING OF A PRESTRESSED CONCRETE BEAM USING INVERSE MODELING TECHNIQUE AND MEASURED DYNAMIC RESPONSE

Keywords: Inverse Problem, Adjoint Problem, Beam Parameters, Dynamic Response

This paper presents an inverse modeling technique using dynamic strain responses to identify the prestressed concrete beam parameters: the pretension in the cable and the prestress in the concrete. The direct problem corresponds to 2D elastodynamics equations. To get the beam parameters, we minimize a classical data misfit function using a gradient-like algorithm. A low-cost computation of the functional gradient is performed using the adjoint problem.

INTRODUCTION

Prestressed beams are widely employed in civil engineering structures. Between 1945 and 1965, 1000 prestressed concrete bridges were constructed, 550 were of the VIPP type (independent span viaducts with prestressed beams) [START_REF] Chatelainl | Post-tensioned concrete bridges[END_REF]. Thus, the maintenance of such structures is a central issue. To reduce maintenance costs, the damaged beams in the structure have to be identified early. To achieve this purpose, the beams should be monitored and non-destructive control methods should be used. Inverse techniques using dynamic responses have been developed in [START_REF] Lu | Identification of Prestress force from measured structrural responses[END_REF] to identify the prestress force in a beam, considering a one-dimensional (1D) beam modelisation. In the present article, we focus on inverse problems based on two-dimensional (2D) beam modelisation. As a matter of fact, by means of 2D modelisation, we can specify the 2D position of the sensors and we can precisely describe the geometry and the position of cables. Previous works [START_REF] Waeytens | Elastodynamics Model Updating for the Monitoring of Reinforced Concrete Beam: Methodology and Numerical Implementation[END_REF] using 2D models and dynamic responses deal with the determinations of the steel bar cross section and the concrete Young Modulus in potential damage area. Using dynamic strain responses, we are interested in identifying the pretension in the cable and the prestressed in the concrete. For that, an inverse modeling technique based on the optimal control theory [START_REF] Lions | Contrôle Optimal des Systèmes Gouvernés par des Equations aux Dérivées Partielles[END_REF] is proposed.

METHODOLOGY TO IDENTIFY THE TENSION IN THE STEEL BAR AND THE PRESTRESS IN THE CONCRETE IN A 2D CONCRETE BEAM

In a 2D concrete beam we aim at identifying the pretension N 0 in the steel bar and the prestress σ 0 in the concrete from the knowledge of measured strain. Note that the dynamic loading applied to the concrete beam is considered known. To determine the beam mechanical parameters, an inverse modeling technique is employed. As usual in the optimal control theory, the pretension N 0 in the steel bar and the prestress σ in the concrete are determined by minimizing a data misfit function. The minimization process is performed in an iterative way using the steepest decent direction (gradient method). At each iteration, the mains steps are the followings:

• solve the direct problem (elastodynamics equations forward in time) considering the beam parameters N 0 and σ 0 of the previous iteration; • solve the adjoint problem (elastodynamics equations backward in time);

• compute the functional gradient using the direct and the adjoint states;

• update the pretension N 0 applied in the steel bar and the prestress σ 0 in the concrete.

DIRECT PROBLEM

In the present study, the shear and the bending work in the cable are neglected. We also neglect the quantity of acceleration in rotation in the cable. We suppose that the cable and the concrete are perfectly adherent. The beam is instrumented with strain sensors. We seek the displacement field

u ∈ U 0 = {u * ∈ H 1 (Ω)\u * = 0 on ∂ Ω 1 , u * .y = 0 on ∂ Ω 2 } such that : Ω ρ c ü.u * dΩ + Γ ρ b S b ü.u * dΓ + Ω ε(u) : K c ε(u * )dΩ + Γ E b S b ∂ u x ∂ x ∂ u * x ∂ x dΓ + Ω σ 0 : (∇u.∇u * )dΩ + Γ N 0 ∂ u x ∂ x ∂ u * x ∂ x dΓ - ∂ Ω f F d .u * ∂ Ω = 0, ∀u * ∈ U 0 (1)
The initial conditions are supposed to vanish. where

σ 0 = σ 0 0 0 0 and σ 0 = - N 0 h.p . 1 - 12.e.y h 2 (2)
• σ 0 : Prestressed concrete;

• N 0 : Tension on the cable;

• h : The height of the beam;

• p : The depth of the beam;

• e : The distance between the middle of the beam and the cable;

For solving the direct problem, the finite element method and Newmark time integration schema are used.

INVERSE PROBLEM

In order to determine the beam parameters, i.e the pretension N 0 in the cable and the prestress σ 0 in the concrete, an inverse modeling technique is employed. As usual in the optimal control theory, we seek N 0 and σ 0 such that it minimizes a data misfit function J. The minimization process is performed in an iterative way using the steepest descent direction (gradient method). The data misfit at sensors i is defined as:

δ i (t) = Ω ε sim xx (x,t)ψ i (x -x i )dΩ -(ε mes xx ) i (t) (3)
We seek the normalized pretension N0 applied in the cable, minimizing the data misfit functional defined as:

J( N0 ) = 1 2 n s ∑ i=1 T 0 [δ i (t)] 2 dt (4)
where :

• n s is the number of strain sensors;

• (ε mes xx ) i (t) and (ε sim xx ) i (t) are the xx measured and simulated strain respectively, by sensor i located at x i ; • ψ i (xx i ) is spatial weight function associated to sensor i.

• N 0 = N0 .N 0 ud , with: N0 ∈ [0, 1] and N 0 ud is the undamaged pretension in the cable.

To get the functional gradient at a low computational cost, we use the adjoint state. We seek the displacement ũ

∈ U 0 = {u * ∈ H 1 (Ω)\u * = 0 on ∂ Ω 1 , u * .y = 0 on ∂ Ω 2 } satisfying: Ω ρ c ü.u * dΩ + Γ ρ b S b ü.u * dΓ + Ω ε( ũ) : K c ε(u * )dΩ + Γ E b S b ∂ ũx ∂ x ∂ u * x ∂ x dΓ + Ω σ 0 : (∇u * .∇ ũ)dΩ + Γ N 0 ∂ ũx ∂ x ∂ u * x ∂ x dΓ - n s ∑ i=1 Ω ( σ0 xx ) i ∂ u * x ∂ x dΩdt = 0, ∀u * ∈ U 0 (5)
where ( σ0 xx ) i is the loading of the adjoint problem which corresponds to a xx internal stress in the vicinity of the i th sensor location, it is defined as:

( σ0 xx ) i = Ω ε xx (x,t)ψ(x -x i )dΩ -(ε mes xx ) i (t) ψ i (x -x i ) (6)
Note that the problem (5) is solved using classical finite element method. At each iteration, we solve the direct problem (elastodynamics equations forward in time), we solve the adjoint problem (elastodynamics equations backward in time) and we compute the functional gradient using the direct and adjoint states. When considering the analytical format of the prestress depending on N 0 , the functional gradient is given by:

∂ J ∂ N0 = T 0 Ω 1 h.p . 1 - 12.e.y h 2 N ud 0 .(∇u.∇ ũ) xx dΩdt - T 0 Γ N ud 0 ∂ u x ∂ x ∂ ũx ∂ x dΓdt (7)
Once we get the functional gradient, several solves of the direct problem are performed to determine the optimal descent step and the new set of beam parameters are deduced.

NUMERICAL EXAMPLE

A 2D concrete beam with a single horizontal pretension cable is considered, Figure 1. A dynamic load F d (t) is applied to the top of the concrete beam, it is defined as;

F d (t) =    -F max sin 2π T c t y, t ∈ [0, T c /2[ 0, t ∈ [T c /2, T ] (8)
In this study, the section of the steel bar S b is assumed constant throughout the beam, the same for the Young Modulus of the concrete E c , we take S b = 0.05 and E c = 27GPa . The Pretension N 0 in the cable and the prestrees σ 0 in the concrete are given by in the equation [START_REF] Lu | Identification of Prestress force from measured structrural responses[END_REF]. We consider seven strain sensors in the concrete beam, there are located at : S 1 : (3L/16, h/8); S 2 : (5L/16, h/8); S 3 (7L/16, h/8); S 4 (L/2, h/8); S 5 (9L/16, h/8); S 6 : (11L/16, h/8); S 7 (13L/16, h/8).

(

) 9
where L is the length of the beam (L = 8m )and h is the height (h = 0.60m).

Solution of the inverse modeling technique

Figure2 illustrates the convexity of the functional J described in equation 4. At the fourth iteration, Figure 3, we note that the error between the the reference pretension (N 0 ) re f and the reconstructed pretension N 0 is about 1% .

CONCLUSION

The identification technique based on 2D elastodynamics model and dynamic strain responses is proposed. On the 2D concrete beam with single cable, the pretension in the cable and the prestress in the concrete are identified after 4 iterations with an error less than 1%. In the futur work, real data outputs will be considered and more representative 2D beam model will be employed in the inverse technique.

Figure 1 :Figure 2 :

 12 Figure 1 : Concrete beam with a horizontal steel bar -Beam instrumented with 7 strains sensors.

Figure 3 :

 3 Figure 3 : The Tension N 0 at each iteration

July 8-11, 2014. La Cité, Nantes, France Copyright © Inria (2014)

EWSHM 2014 -Nantes, France

EWSHM 2014 -Nantes,

ACKNOWLEDGMENTS

This research was supported by the region Ile de France. This work was conducted as part of the project "SIPRIS", which involves several partner: Advitam (leader), ASF (Highways of the South of


S F Zedek

C Codreanu

A Belisario

J-Y Fourniols

C Escriba

T Camps

R Francois

ADAPTABLE SYNCHRONOUS DETECTION FOR THE SHM OF COMPOSITE COMPONENT: AUTONOMOUS FPAA ARCHITECTURE ON POLYIMIDE FILM

Keywords: FPAA1, Kapton2, SHM3, Strain Gauge, Synchronous detection, CFRP6

Structural Health Management strategies usually involve monitoring complex components, during their working life. Carbon fiber reinforced polymer components are becoming more and more used, throughout the structural engineering fields, because of their superior mechanical properties. Moreover, due to their construction, these materials open the path for new "smart" components, with integrated sensors and monitoring capabilities. Flexible circuits are essential in this context. In order to adapt to complex geometries, sensors and their connectivity have to be flexible and as unintrusive as possible. Following previous work, aimed at inserting sensors with flexible connectivity in between the sheets of carbon fibers, we conceived and implemented a reconfigurable acquisition system that can be easily adapted to perform strain measurements in various EM5 noise conditions. This paper presents a comparison between three circuits dedicated to strain measurement using a commercial strain gauge and flexible connectivity on Kapton2. The first circuit consists of a classic Wheatstone bridge under continuous voltage. The second uses the synchronous detection principle with classic circuits and functions. The third innovates upon the synchronous detection principle by implementing several functions in a reprogrammable analog chip called FPAA (Filed Programmable Analog Array). The analysis of these circuits shows the added value of using reconfigurable devices with flexible connectivity and sensors in order to provide a universal solution for the monitoring of complex components.

INTRODUCTION

tructural Health Monitoring usually refers to the process of implementing a damage detection strategy for aerospace, civil or mechanical engineering [START_REF] Charles | Structural Health Monitoring a Machine Learning Perspective[END_REF]. The human being is surrounded by buildings, bridges and complex structures which have to be safe for use. In order to comply with safety regulations several techniques have been developed for monitoring these structures throughout their lifetime.

Almost all private and government industries want to detect damage in their products as well as in their manufacturing infrastructure at the earliest stage possible. Such detection requires these industries to perform some form of SHM and is motivated by the potential life gain, safety and economic impact of this technology S [START_REF] Charles | Structural Health Monitoring a Machine Learning Perspective[END_REF].

Strain gauges are widely used for SHM in civil engineering because of their low cost, good sensitivity and simple usage. They exist in different sizes, forms and configurations, allowing the measure of strain and torsion. They are usually implemented in thin polymer films and are un-intrusive.

Previous work [START_REF] Torres | Valeur ajoutée d'une encapsulation de capteurs silicium pour l'instrumentation à coeur des structures composites[END_REF] has shown the feasibility of inserting sensors with Kapton connectivity in between the carbon fiber layers of a CFRP component. The circuits fared well during the high temperature molding process and the sensors were exploitable. The Kapton connectivity also withstands folding, allowing the creation of complex paths through a component. Figure 1 shows a pre-mold setup, with sensors and flexible connections placed on the second layer of carbon fibers, before they are covered with 8 more layers and molded at high temperature. Using the gauges to monitor strain applied to the board requires a dedicated circuit capable of transforming a small change in electrical resistance into a voltage level. The starting block of this circuit is the Wheatstone bridge under continuous voltage. This circuit is the first one we created and analysed in order to have a reference set of characteristics. An improvement over the continuous voltage setup involves working at a high frequency (100kHz). This shifts the working spectrum away from 0 Hz and into an area where the sensor will be less sensitive to noise. The synchronous detection was implemented with classic components and with a reprogrammable analog chip. These two circuits were tested and benchmarked against each-other and against the Wheatstone bridge under continuous voltage.

For testing purposes we constructed an experimental bench using a strain gauge glued onto a carbon fibre rod. The rod is fixed at one end onto a stable support, as shown in figure (2) The following recordings of free oscillations were made during a step release test using this bench.

EWSHM 2014 -Nantes, France 1 Implementation of the proposed circuits

The Wheatstone bridge under continuous voltage

This type of bridge is widely used for the measurement of sensors [START_REF]Dynamically programmable analog signal processing[END_REF]. As shown in the figure below the strain gauge is placed in one of the four branches of the bridge. Its nominal resistance is 121Ω and a trimmer was used in parallel on the left hand side in order to precisely balance the bridge (not shown in the figure). Once the bridge is properly balanced, VA = VB where:

) /( ) /( 4 3 4 1 R R R E V R R R E V B SG SG A      
A small change RSG caused by mechanical strain produces an offset:

SG SG B A M R R E R V V V       1
This signal is extracted using an INA128 instrumentation amplifier. It is amplified by a factor of 500 using low-noise amplifiers. The final stage is a low-pass second order filter with a cut-off frequency of 2.5 kHz.

Figure 4 shows a recording of the free oscillations of the component under test during a step release test using this circuit.

The synchronous detection circuit

The synchronous detection is a signal processing technique for extracting a weak signal flooded in noise. The following block diagram details its main components. The useful signal is shifted to a frequency, imposed by the wave generator Fm, where the noise level is lower. The modulated signal is amplified; filtered and then demodulated. After the multiplier, the low-pass filter removes the 2Fm component. In order to extract the sensor's signal, we use a modified version of the Wheatstone bridge, the Maxwell Bridge. Capacities were introduced in parallel to each resistance on the left hand side of the bridge (not shown on figure 6) in order to balance the complex impedances of the components. Again, the bridge's offset is directly proportional to the change in gauge resistance. Due to the complex equivalent impedance, the signal VM = VA -VB changes both its amplitude and phase. Again, we used an INA128 instrumentation amplifier to extract the VM signal. There is a pre-amplification step by a factor of 10 before filtering. Following filtering the signal is again amplified by a factor of 50 and filtered in order to remove frequencies above 2.5 kHz.

Since our objective is to compare the discrete architecture with the reprogrammable architecture, we sized the filters such that the circuits are similar in terms of bandwidth. The following figure shows the band pass of the filter created for this circuit and the one programmed in the following circuit :

Reprogrammable synchronous detection circuit

As seen previously, one of the key parameters of the synchronous detection system is Fm, the work frequency. This is usually chosen depending on the strain gauge's noise profile but external interferences must also be taken into account. Changing this frequency in the previous architecture implies modifying a big part of the design since it requires changing the wave generator and the band pass filter. The re-programmable analog chip FPAA offers an excellent solution to this problem. FPAA stands for Field Programmable Analog Array. It can incorporate both the oscillator and the band pass filter and both can be changed programmatically without modifying the circuit board. Furthermore, since the analog multiplier is one of the most expensive components on the board, it's worthwhile to remove it. The FPAA is capable of incorporating a multiplier CAM with equivalent performances. CAM stands for Configurable Analog Module. FPAA is an analog reprogrammable device based on operational amplifiers and switch capacitors distributed in four Configurable Analog Blocs, CABs. Figure 9 is a diagram of the internal architecture of the AN231E04. It can be used for low power circuits in many electronic applications for filtering and adaptable design. In our application, we used it for signal conditioning. Anadigm provides the AnadigmDesigner®2, a software package for the users to configure, simulate and program their FPAA boards. We used the AN231E04 and we configured it with the following schematic: For the band pass filter we used the included "Filter Designer" and we chose the following parameters:

Experimentally we measured the magnitude graph shown in figure 7. This band pass filter occupied two operational amplifiers, eight switched capacitors in one CAB [START_REF]Dynamically programmable analog signal processing[END_REF] and it requires 55±17 mW of power. This leaves enough room in the CABs for two additional filter CAMs. The programmed FPAA is inserted into the synchronous detection circuit as shown on the following schematic:

Noise analysis

Strain gauges are mainly used for long term surveillance of quasi-static low-amplitude deformations. Thus they must be stable in time and have as little low frequency noise as possible. One of the main advantages of the synchronous detection system is the improved signal/noise ratio. In order to quantify this we measured the noise power spectrum. Figure 13 shows the results for the three circuits under test: As expected, the noise level for the synchronous detection circuit is lower than the one obtained with the bridge under constant voltage. On the other hand, using the reprogrammable system introduces more noise and cancels the improvements achieved by the synchronous detection. At low frequencies ( < 1 Hz ), the noise traces have the same shapes but using the FPAA removes the stray frequencies observed for the two previous circuits. The noise level could be further reduced by directly integrating the FPAA onto the printed circuit board instead of linking with the evaluation board.

Power ratings

The following table shows the measured consumptions of the three circuits under test:

Passing from the Wheatstone bridge to the synchronous detection circuit the power consumption increases slightly due to the addition of several active components. As shown by the noise spectrum, this is a good compromise. Using the reprogrammable circuit requires three times more power. This is due to the addition of several active components required by the multiple interfaces with the FPAA, but the biggest power consumers are the auxiliary systems present on the development board ( USB controller, RS232 programmer, EEPROM, PIC16F876A ) which are not mandatory for our application.

Conclusions

In our work we demonstrated the feasibility of implementing synchronous detection using a reprogrammable analog chip. This brings many advantages in terms of design adaptability: the working frequency can be easily changed, the bandwidth of the filters can be adjusted to comply with a specific need, additional functions can be programmed into the FPAA. The main drawbacks of using the FPAA are the noise level and the power consumption. The noise level can be improved by directly integrating the FPAA onto the circuit board. Also power consumption can be reduced by eliminating several auxiliary systems. Furthermore, changing the design so that it works on an asymmetric 0 -5V supply and exploiting the fully differential in-out of the FPAA, could help reduce power consumption by a factor of 2. Flexible circuits will change structural engineering by allowing sensors and circuitry to be directly integrated into the building blocks of mechanical systems. Reconfigurable circuits will play a key role in this technological advance by smoothing the path from the sensors to the data processing and the monitoring strategies which are required to ensure the proper operation of increasingly complex systems.

Figure 1 :

 1 Figure 1: Flexible circuits inserted in between the layers of carbon fibers before molding. Also shown in the bottom left are the external connection pads. In the bottom right is shown the Kapton folding.

Figure 2 :

 2 Figure 2: Experimental setup with the sensor shown in red.

Figure 3 :

 3 Figure 3 : Schematic of a Wheatstone bridge under continuous voltage.

Figure 4 :

 4 Figure 4 : Step response using the Wheatstone bridge under constant voltage.

Figure 5 :

 5 Figure 5 : Schematic of the synchronous detection system

Figure 6 :

 6 Figure 6 : Circuit schematic for the synchronous detection

Figure 7 :

 7 Figure 7 : Band pass filters magnitudes

Figure 8 :

 8 Figure 8 : Step response using the synchronous detection circuit

Figure 9 :

 9 Figure9: Architectural overview of the AN231E04 device with dynamic re-configurability[START_REF]Structural health Monitoring System based on Strain gauge enabled wireless Sensor Nodes[END_REF]

Figure 10 :

 10 Figure 10 : Programmed circuit using the AnadigmDesigner suite.

Figure 11 :

 11 Figure 11 : Synchronous detection integrating the reprogrammable circuit.

  The step response of the experimental setup is shown in the following figure:

Figure 12 :

 12 Figure 12 : Step response using the reprogrammable synchronous detection circuit.

Figure 13 :

 13 Figure 13 : Noise power density spectrum

TABLE I PASS

 I BAND FILTER PARAMETERS

July 8-11, 2014. La Cité, Nantes, France Copyright © Inria (2014)

EWSHM 2014 -Nantes, France

Acknowledgment

We take this opportunity to thank we would like to thank Mr. P. Marcoul, working with us here at LAAS-CNRS, for routing the circuits and fabricating the printed circuit boards.