Roilhi F Ibarra-Hernández

Nancy Bertin

Miguel A Alonso-Arévalo

Hugo A Guillén-Ramírez

R F Ibarra

A benchmark of heart sound classification systems based on sparse decompositions

Keywords: Heart sounds, Phonocardiogram, Matching Pursuit, Linear Predictive Coding, Cross Validation test, Oversampling, SMOTE

Background: Nowadays, cardiovascular diseases (CVD) remain the main cause of death worldwide. A heart sound signal or phonocardiogram (PCG) is the most simple, economical and non-invasive tool to detect CVDs. Advances in technology and signal processing allow the design of computer-aided systems for heart illnesses detection from PCG signals.

Purpose: The paper proposes a pipeline and benchmark for binary heart sounds classification. The features extraction architecture is focused on the use of Matching Pursuit time-frequency decomposition using Gabor dictionaries and the Linear Predictive Coding method of a residual. We compare seven classifiers with two different approaches: feature averaging and cycle averaging. Methods: We test our proposal on the PhysioNet/CinC challenge 2016 database, which comprises a wide variety of heart sounds recorded from patients with normal and different pathological heart conditions. We conduct a 10-fold stratified cross-validation method to evaluate the performance of different classification algorithms. The feature sets were also tested when using an oversampling method for balancing. Results: The benchmark identified systems showing a satisfying performance in terms of accuracy, sensitivity, and Matthews correlation coefficient. Results can be improved when using feature averaging and an oversampling strategy.

INTRODUCTION

Cardiovascular diseases (CVDs) are one of the main causes of death globally according to the World Health Organization. [START_REF] Who | Cardiovascular diseases, fact sheet[END_REF] A valuable method for CVDs detection is auscultation, i.e. listening to the mechanical valvular activity. When recorded, the resulting sound sequence is called a phonocardiogram (PCG). Auscultation is in fact a primary diagnosis method for initial detection of heart valves malfunctioning. It is the most economical and simple screening test, since the only device needed is a stethoscope. However, this technique has restrictions, since high quality recordings are required for reliable interpretations. Another constraint is the amount of experience and training to master auscultation skills for physicians, who are few in quantity and often located in urban areas. Nowadays, with the recent advances in signal processing and electronic stethoscope technologies, the design of automatic classification schemes from PCG recordings appears as a promising diagnosis method for CVDs detection. This paper aims at defining end-to-end pipelines for automatic classification of heart sounds as healthy or pathological, and at benchmarking instantiations of this pipeline with several signal representations and classifiers, in order to assess their potential in clinical practice.

Heart sound signals

In a healthy or normal state, a PCG recording comprises two main components called fundamental heart sounds (FHS) and denoted S 1 and S 2 . Each FHS can be characterized by a common time length and energy concentration over low frequency regions. For instance the S 1 components are dominant in the region of 10 Hz to 140 Hz while S 2 components usually concentrate their energy around the 10 Hz to 200 Hz band. [START_REF] Abbas | Phonocardiography signal processing[END_REF] Murmurs are sounds stemming from a turbulent blood flow due to valve malfunction, hence denoting a pathological or abnormal state. The energy distribution of murmurs in frequency vary widely and depending on its nature they can be as high as 600 Hz. There exists other sounds called S 3 and S 4 which can represent a pathological state in adults, but are physiological in children. The frequency content of these sounds may overlap with the S 1 and S 2 energy distribution. Figure 1 illustrates the waveform and the time-frequency content representative of a PCG cardiac cycle (union of S 1 and S 2 sounds) in both normal and pathological states. Murmur detection is a key task in the primary diagnosis of pathological states from a PCG recording. However, murmurs are highly nonstationary signals. In addition, PCGs are often corrupted by noises, such as those arising from speech or stethoscope movements. As a result, automatic heart sounds classification is a promising, but challenging task.

Prior art

The 2016 PhysioNet/Computing in Cardiology Challenge (CinC) was open in 2016 to encourage researchers to design and train classification methods for the binary labeling of PCG sounds as in normal or abnormal conditions. [START_REF] Clifford | Classification of normal/abnormal heart sound recordings: The physionet/computing in cardiology challenge 2016[END_REF] It provides the largest public database of PCG recordings registered in a real clinical environment. [START_REF] Liu | An open access database for the evaluation of heart sound algorithms[END_REF] As for any classification task, each approach submitted to the challenge can be mainly characterized by two components: the chosen signal representation (features) and the type of classifier. During the CinC 2016 conference, most of the participants have used time-frequency features, in particular the wavelet decomposition coefficients and Mel-Cepstral Frequency Coefficients (MFCCs). [START_REF] Abdollahpur | Cycle selection and neuro-voting system for classifying heart sound recordings[END_REF][START_REF] Munia | Heart sound classification from wavelet decomposed signal using morphological and statistical features[END_REF][START_REF] Zabihi | Heart sound anomaly and quality detection using ensemble of neural networks without segmentation[END_REF] Linear predictive coding coefficients have been also used. [START_REF] Zabihi | Heart sound anomaly and quality detection using ensemble of neural networks without segmentation[END_REF] In terms of classifiers, competitors mostly used Neural Networks, [START_REF] Abdollahpur | Cycle selection and neuro-voting system for classifying heart sound recordings[END_REF][START_REF] Kay | Dropconnected neural network trained with diverse features for classifying heart sounds[END_REF][START_REF] Potes | Ensemble of feature-based and deep learning-based classifiers for detection of abnormal heart sounds[END_REF] Support Vector Machines [START_REF] Whitaker | Combining sparse coding and time-domain features for heart sound classification[END_REF][START_REF] Bobillo | A tensor approach to heart sound classification[END_REF] and Random Forest techniques. [START_REF] Homsi | Automatic heart sound recording classification using a nested set of ensemble algorithms[END_REF][START_REF] Antink | Nonnegative matrix factorization and random forest for classification of heart sound recordings in the spectral domain[END_REF][START_REF] Singh-Miller | Using spectral acoustic features to identify abnormal heart sounds[END_REF] In parallel, sparse representations of the heart sounds (Matching Pursuit decomposition and Linear Predictive Coding of the residual) have been shown to provide a compact and meaningful representation of PCG signals. [START_REF] Ibarra | A parametric model for heart sounds[END_REF] To our best knowledge, the usefulness of such a representation in a heart sounds classification task has not been assessed so far.

Contributions

This paper proposes a unified but modular pipeline for the development and the assessment of heart sounds classification systems based on sparse representations. A global overview of the pipeline and details on all its modules are given in Sec. 2. We provide a systematic evaluation of the derived systems as well, with a particular emphasis on the comparison between various design choices: cycle averaging vs. feature averaging, optional use of data oversampling for unbalanced classes compensation, and a wide range of classifiers. Experimental results, presented in Sec. 4 and discussed in Sec. 5, are obtained from a 10-fold Cross Validation test outputting accuracy, sensitivity and specificity scores, and Matthews Correlation Coefficients. This exhaustive evaluation allows us to conclude in Sec. 6 on the promising outperformance of a system combining cycle averaging, sparse representation, LPC residual coding and a Random Forest classifier, even in adverse conditions such as unbalanced classes and noise-corrupted recordings.

METHODS

Features extraction

This paper compares two design choices for features extraction in order to represent each PCG: feature averaging and cycle averaging. Each approach gives as output two feature sets: A and B. Both methods share the following steps:

• Cycle segmentation: It consists in the time delineation of heart sound cycles. For this purpose we took the corrected annotations from Physionet/CinC challenge database. [START_REF] Clifford | Classification of normal/abnormal heart sound recordings: The physionet/computing in cardiology challenge 2016[END_REF] Original annotations came as an output from Springer's method [START_REF] Springer | Logistic regression-hsmm-based heart sound segmentation[END_REF] and then corrected by the sponsors.

• Matching pursuit (MP) decomposition: It performs the time-frequency representation of a cycle or a FHS. [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] Features from MP are given by the parameters of the selected atoms from MP algorithm when using a Gabor dictionary (see Sec. 2.1.2).

• Linear predictive coding (LPC): Due to its greedy nature, MP does not provide a perfect reconstruction of the signal. The difference between the MP reconstructed signal and the original PCG is called the residual. LPC models this residual in a compact form (see Sec. 2.1.3).

Fig. 3 describes the methods developed in this work to extract the feature sets A and B. The first method is shown on the left diagram and depicts the feature averaging approach. The output feature set A is actually obtained from the mean value of the MP and LPC parameters for the N segmented PCG cycles. The diagram on the right depicts a second method which uses a cycle averaging approach. MP and LPC parameters are extracted from an averaged cycle. By contrast with the first method (in which parameters are directly extracted from the whole cycle), in this second method, the averaged cycle is split into the two FHS. MP is then applied on each FHS separately, and two frequency parameters F 1 and F 2 are kept for each of them, as in a previous methodology.

Preprocessing

From observations of heart sounds pointed on 1.1, we introduce a band-pass filter between 25 Hz and 600 Hz to remove non informative frequencies in the input PCG signal. The original sampling rate of 4 kHz of recordings is preserved and the filter implemented is a sixth-order Butterworth filter.

Matching Pursuit (MP)

The Matching Pursuit method is a greedy and iterative algorithm which aims to find a linear combination of D elementary waveforms called atoms which approximates a signal x with minimal error. Each atom g d is a elementary signal belonging to a redundant set of all possible predefined atoms called dictionary D. MP iteratively selects atoms g d to provide a decomposition of the form:

x = D d=1 α d • g d + r, (1)
where the coefficient α d is a scalar weighting factor and r is a residual term. When using a dictionary of Gabor functions, the MP decomposition derives an adaptive time-frequency transform, [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] since this dictionary is composed by well concentrated waveforms in the time-frequency plane. In this work we use a pre-defined set of multiscale functions which is a collection D = ∪ J j=1 D j of blocks D j of time-frequency atoms at different scales. The waveform of a Gabor atom in a multiscale dictionary is defined by the modulation, dilation, translation and sampling of a (continuous) window w j (t) as:

g j,n,k (m) = w j (mT s -nT j ) exp 2iπkmT s K j for 1 ≤ m ≤ M, (2)
where the time location or window shift is defined as nT j , the window length or scale L j and is modulated at a frequency k/K j , K j is a predefined number of possible frequencies (FFT size), T s is the sampling period and M the number of time samples. Fig. 4 illustrates the waveform of a Gabor atom, which is actually a cosine modulated Gaussian window. Other waveforms are shown at the right to see the effect of changing the modulation frequency.

{ Figure 4
. Time waveform of a Gabor atom and its defined parameters. The right panel illustrates the resulting waveforms when frequency parameter k/Kj varies.

Linear Predictive Coding (LPC)

As seen in Eq.1, after MP, the initial signal is decomposed into two main parts, a linear combination of Gabor atoms and a residual. The residual term r is expected to be rather uncorrelated with atoms, and thus has to be represented differently to be integrated in the feature sets. Instead of representing the temporal waveform, LPC tries to fit the spectrum of the signal. A number of works have used the LPC parameters [START_REF] Zabihi | Heart sound anomaly and quality detection using ensemble of neural networks without segmentation[END_REF][START_REF] Redlarski | A system for heart sounds classification[END_REF] as features for PCG's classification. Let h k be the set of LPC coefficients. They define a filter called predictor. The principle is that the residual r can be predicted as a linear combination of the previous samples:

r n = - p i=1 h i r n-i + e n , (3)
where n = 0, • • • , N -1 and e n is the final residual. Filter coefficients h i are added to the feature set.

EXPERIMENTAL SETTINGS

Heart sounds database

The set of PCG recordings from the Physionet/Cinc Challenge 2016 4 was employed to conduct the experiments in this work. This set is available online and contains 3,153 recordings from 764 subjects as training set. It has been assembled from six different subsets comprised by PCG sounds recorded under different conditions. Each recording has been labeled as normal or abnormal according to the presence or not of a murmur. Another given label is the quality of the recording as good or bad. However, this dataset is unbalanced since the class distribution of the signals is uneven: 665 recordings were marked as abnormal while 2,488 as normal.

Feature extraction settings

The MP decomposition of PCG signals was configured in the Matching Pursuit Toolkit library (MPTK) designed for MATLAB. [START_REF] Krstulovic | MPTK: Matching pursuit made tractable[END_REF] We used a Gabor dictionary of J = 5 blocks, each block corresponding to a common atom lengths L j of 32, 64, 128, 256 and 512 samples. The selected number of atoms was M = 15 in order to reach almost a 99% of the energy to reconstruct a PCG cycle. [START_REF] Ibarra | A parametric model for heart sounds[END_REF] The atom parameters of frequency, amplitude, length, position (shift) and phase were extracted to be considered as features in features set A while in features set B just two atom frequencies were considered.

For the LPC representation of the residual signal output by MP, we used the MATLAB code from UCLA available on-line. [START_REF] Ozun | Wideband speech coding with linear predictive coding (lpc)[END_REF] The selected order of the filter was p = 15. Both feature sets A and B do not have the same number of samples. Combining the MP and LPC methods, features set A consists on N f eatures = 90 resulting from 75 parameters from MP (15 atoms with 5 parameters each one) and 15 from LPC. This choice is based on our previous work. [START_REF] Ibarra | A parametric model for heart sounds[END_REF] On the other hand, features set B contains N f eatures = 19, among which 4 are provided by MP (frequencies of the first two selected atoms for each FHS). This setting, in addition to the rest of set B extraction scheme (FHS segmentation and feature averaging), is in line with state-of-the-art. [START_REF] Wang | Analysis of the first heart sound using the matching pursuit method[END_REF][START_REF] Wang | Decomposition and analysis of the second heart sound based on the matching pursuit method[END_REF] Thus, set B acts as a comparison baseline for set A, whose performance in heart sound classification has not been assessed so far.

Classifiers settings

We tested the classification state of our proposed pipeline using seven different methods. Table 1 presents a brief description of each one by its name, acronym and the main parameters employed. The classifiers were configured according to the values depicted in Table using the sci-kit learn toolbox of Python. [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] In addition, feature sets were normalized to have zero mean and unit variance when using SVM. For the RF method we changed the number of estimators to 100 as some authors recommend in presence of unbalanced problems. 25, 26

Noisy recordings

Some of the recordings provided by Physionet database do not contain any annotation file to perform segmentation, due to their noise level * . To handle these files, we performed the processing of the signals in frames with time lengths of 900 ms (with respect to the approximate time duration of a cycle) and 200 ms (according to the typical length of a FHS.)

Classification performance

In order to measure the classifiers performance, the following metrics were calculated: accuracy (ACC), Sensitivity (SE), Specificity (SP) and Matthews Correlation Coefficient (MCC), defined as:

ACC = T P + T N T P + F N + T N + F P (4) SE = T P T P + F N × 100 (5) SP = T N T N + F P × 100 (6)
M CC = T P × T N -F P × F N (T P + F P )(T N + F P )(T P + F N )(T N + F P ) , (7)
where the values T P , T N , F P and F N corresponds to the number of true positives, true negatives, false negatives and false positives respectively † . * These were found as low quality, actually. † In our case a true positive (T P ) corresponds to a PCG entry with an abnormal condition correctly predicted.

RESULTS

Classifiers performance test without class balancing

We conducted in sci-kit learn a 10-Fold CV approach. For the first test, data is plural. Table 2 shows the SE, SP , ACC and M CC average and standard deviation output metrics. In terms of SP , the SVM method reached the highest score of 92.29% when using features set A. However, the SE score reached by this combination is the lowest, even the M CC is not the highest presented. The RF method reached the highest values for the remaining metrics when using features set A: SE = 75.80, M CC = 0.55 and ACC = 0.86. Results from Table 2 present also greater standard deviation values for SE. We observe that for all feature sets and all classifiers, SP is considerably higher than SP .

Classifiers performance test with class balancing

For a second experiment, we perform a balancing applying the SMOTE technique when oversampling the minority class. The SMOTE library included in the imbalanced-learn toolbox of Python 34 was used for this purpose. Then, we conducted again the stratified 10-fold CV test to evaluate the classifiers performance. Table 3 shows the average metrics that outcomes from this experiment. Compared to the metrics obtained without applying SMOTE, the SE values showed an increment for all tested classifiers and feature sets. The highest sensitivity SE = 91.60% value was reached for the combination of the features set A and the RF method. This approach also reached the highest accuracy ACC = 92% and Matthews Correlation Coefficient M CC = 0.84. Although KNN classifier with the features set B presented the highest SP , the remaining scores were not the highest values. Fig. 5 plots all the ACC scores obtained. An increase in the mean values is shown, except for the NB method and SVM when using as input the features set B.

DISCUSSION AND FUTURE WORK

Overall, the outcomes of the proposed classification schemes indicate that a combination of MP time-frequency features and LPC features, oversampling, and state-of-the-art classifiers result together in a satisfying performance for heart sound classification.

The obtained performances can be indicatively related to those reported for the 2016 PhysioNet/CinC entrants. [START_REF] Clifford | Classification of normal/abnormal heart sound recordings: The physionet/computing in cardiology challenge 2016[END_REF] In the original challenge, performance was assessed by modified versions of SE, SP , and their arithmetic mean, denoted M acc ‡ . The highest ranked system, 9 based on a convolutional neural network, reached a mean accuracy M acc = 86.02% with SE = 94.24% and SP = 77.81%. In comparison, in our experiments on similar data § , the random forest classifier with feature set A and SMOTE oversampling reached a sensitivity SE = 91.60% and a specificity SP = 92.10%, corresponding to a mean accuracy of M acc = 91.85%. This suggests that some of the approaches presented here could be competitive with state-of-the-art systems submitted to the challenge.

Several novelties of the proposed approaches, emphasized below, contribute to their performance. Their analysis allows us to draw additional insights on heart sound classification strategies and to identify axes for future research.

Features

Several differences between feature sets A and B, and features used as input in state-of-the-art systems, can be stressed out:

• Compared with most of Physionet challenge entries (wavelets, time lengths of FHS and murmurs and MFCC coefficients mainly), the features extracted in both sets came from an MP-based approach using Gabor dictionaries. This proposal has been previously shown to provide with an efficient representation of PCG signals. [START_REF] Zhang | Analysis-synthesis of the phonocardiogram based on the matching pursuit method[END_REF][START_REF] Sava | Application of the matching pursuit method for structural decomposition and averaging of phonocardiographic signals[END_REF] • The LPC extracted features came from a representation of the output residual signal from the MP decomposition, instead of direct modeling of cardiac sounds waveform. In previous work we showed that this representation of the residual preserves some meaningful heart sound signal components. [START_REF] Ibarra-Hernández | Design and evaluation of a parametric model for cardiac sounds[END_REF] • We considered two main approaches to build and extract PCG signal feature sets: feature averaging and cycle averaging. Most of Physionet challenge submissions focused on the first one, while the second one had not been evaluated with these recordings so far. ‡ It must be noted that metrics are slightly different from ours, calling for cautious interpretations. Contrary to the cited paper, we did not weight SP and SE according to the data labels as indicated for the Physionet/Cinc 2016 challenge. Little difference is expected from this. On the other hand, the average M acc is different from the ACC criterion, which puts more importance on the SE/SP trade-off, and they cannot be directly compared.

§ Our experiments were conducted exclusively on data included in the Physionet database. However, due to the unavailability of the two "hidden test sets", they are not strictly identical.

All else being equal, systems using the set A as input systematically outperformed their equivalent with the set B, although the number of features in A is greater. An additional test on cumulative feature importance for RF classification was performed (data not shown). This test revealed the importance of adding the LPC features: in set A, although they represent only about 17% of the input feature vector, 5 of them were ranked in the 10 most important features for the RF classifier. In addition, our results validate the new feature averaging approach against the previous cycle averaging approach.

However, as a 90-dimensional feature space remains demanding and possibly suboptimal (due to correlations between features), future work may include a feature selection step, to better deal with small training datasets and reduce the computational burden.

Oversampling

Unbalance between healthy and pathological classes is a typical situation in medical applications. It is known, in particular, to be very detrimental to sensitivity of classifiers, and to lead to uneven SE -SP trade-offs.

Our results clearly confirm the dramatic impact of using an oversampling strategy to alleviate this problem. We observe an average improvement of 24.4 points in sensivity scores, and a global reduction of the variance on this indicator. In most cases, this is accompanied with a moderate variation of specificity (minor loss or improvement). Hence, the use of SMOTE brings the majority of our best systems to reach higher ACC and M CC scores, which denotes a more balanced trade-off between SE and SP , and a better global performance.

Exceptions are the NB classifiers (with both feature sets), SVM classifier (with set A only) and the LR and LDA classifiers (with set B only) which apparently don't benefit from oversampling, encountering an important drop in specificity. However, none of these classifiers ranked first in any of our experiments.

Classifiers and their tuning

As mentioned above, the RF classifier showed the higher metrics for both experiments, reaching a mean ACC of 92%, SE = 91.60%, SP = 92.10% and M CC = 0.84 when using the feature set A and oversampling. Without oversampling, RF still ranks first with set A and first ex-aequo with set B (together with KNN.)

The achieved SE -SP trade-off varies widely between classifiers. The highest specificity is achieved by KNN with oversampling (97.1% with set A and 95.87% with set B), at the price of a poor sensitivity. This may be linked to the relatively low number of considered neighbors (5) with respect to the feature set size (90). As sensitivity is concerned, RF ranks first again. This plays a large role in its global score, since sensitivity is obviously damaged by unbalanced classes for all classifiers, and the criterion that most of the tested classifiers struggle to fulfill.

It should be noted that all performance scores, but more particularly the SE -SP trade-off, is known to be also largely ruled by the classifier's parameters. In our experiments, their values were chosen from default values and literature, and were not tuned nor optimized. Further work would be needed to explore their impact on performance, defining tuning strategies and identifying the classifiers' various operational regimes. This is all the more relevant than in clinical practice, the smallest SE -SP difference is not always seeked for. As noticed in the original challenge, 3 as false positives and false negatives have very different consequences on patients, physicians may favor one of the criteria over the other, depending on the clinical situation, other existing complementary tests and costs. As PCG classification is meant to act as a first stage screening test, one could expect that a better sensitivity would be preferred, even at the cost of a loss in specificity. This choice, which can be informed by benchmarks like the present one, however belongs to physicians and public health policy experts.

CONCLUSIONS

This paper proposes a methodical benchmark of several phonocardiogram binary classification systems, aiming at detecting pathological heart states. Two sets of time-frequency features were extracted from each heart sound, based on the MP decomposition and the LPC coding technique. [START_REF] Ibarra | A parametric model for heart sounds[END_REF] Seven common state-of-art classification methods were then tested in a 10 fold stratified cross validation technique. An oversampling technique, SMOTE, was also added in order to compensate for unbalanced classes (less recordings with an abnormal condition, which is typical in biological data), and its impact on classifiers performance was assessed. Performance was measured through sensitivity, specificity, accuracy and Matthews correlation coefficient. Among all conducted experiments, the Random Forest classifier combined with SMOTE technique outperformed all other configurations, reaching a competitive score for accuracy (92%), with a good sensitivity-specificity tradeoff. Detailed insights from the experiments also validate the newest of the two used feature sets in the considered classification task. Its originality mainly relies on three novelties: an alternative strategy for feature extraction from a long recording (feature averaging) to the more common heart cycle averaging strategy; it does not require a segmentation of the fundamental heart sounds S1 and S2 within one heart cycle; the LPC coefficients that it includes are not computed on the whole signal, but only on the residual signal left after a sparse decomposition of the PCG in a Gabor dictionary.

As the identified scheme shows promising results for detecting abnormal conditions in PCG signals, future work will target performance and computational efficiency improvements through a refined preprocessing (denoising) and dimensionality reduction approaches (feature selection). Although excluded from our benchmark due to their high training data requirements, neural network based classifiers cannot be ignored nowadays and should also be tested against the proposed feature set.

Answers to reviewers questions

We want to thank the reviewers for the comments made to our manuscript, SIP300-14, A benchmark of heart sound classification systems based on sparse decompositions. We completed the submission successfully taking into account these valuable observations. In this section, we would like to ask a couple of questions formulated during the review.

Question: The authors mention that murmurs are highly nonstationary signals, but aren't S1 and S2 sounds also nonstationary signals?

Our answer: We consider that murmurs are, by nature, highly nonstationary because of their characteristic broad and abrupt frequency changes. On the other hand, S1 and S2 sounds are quasi-stationary due to their waveform shape and limited characteristic bandwidth. However, time location and delineation in these sounds vary according to physiological conditions from each patient.

Question: It could be desirable a discussion about the advantages of the classifiers analyzed over the more recent techniques as deep learning or recurrent neural networks.

Our answer: Deep neural networks usually require much more data than traditional Machine Learning algorithms, at least datasets containing hundreds of thousands if not millions of labeled samples. Unfortunately, we do not have access to such amount of annotated cardiac sounds. Training of deep learning or neural networks is not a trivial task and due to the amounts of data required the implementation of the methods can be hard to parallelize [1].

[1] Chen, X. W., & Lin, X. (2014). Big data deep learning: challenges and perspectives. IEEE access, 2, 514-525.

Figure 1 .

 1 Figure 1. The time waveform and spectrogram representative of a PCG cardiac cycle in normal (top) and abnormal (bottom) conditions.

Fig. 2

 2 Fig.2gives an overview of the proposed pipeline for heart sounds classification. Features extraction step gives two different sets of features: A and B (see Fig.3for details). To address the probable case where the data would be unbalanced (more recordings labeled as normal than as abnormal in the training set), the Synthetic minority oversampling technique (SMOTE)[START_REF] Chawla | Smote: synthetic minority oversampling technique[END_REF] can be employed. After the SMOTE procedure, a 10-fold cross validation (CV) test is performed for different classifiers (summarized in Tab.1.) As described in Fig.2, each fold operates on a different training/test data splitting (one tenth is used for test and the remaining for training).

Figure 2 .

 2 Figure 2. Overview of the proposed pipeline.

Figure 3 .

 3 Fig.3describes the methods developed in this work to extract the feature sets A and B. The first method is shown on the left diagram and depicts the feature averaging approach. The output feature set A is actually obtained from the mean value of the MP and LPC parameters for the N segmented PCG cycles. The diagram on the right depicts a second method which uses a cycle averaging approach. MP and LPC parameters are extracted from an averaged cycle. By contrast with the first method (in which parameters are directly extracted from the whole cycle), in this second method, the averaged cycle is split into the two FHS. MP is then applied on each FHS separately, and two frequency parameters F 1 and F 2 are kept for each of them, as in a previous methodology.[START_REF] Wang | Analysis of the first heart sound using the matching pursuit method[END_REF][START_REF] Wang | Decomposition and analysis of the second heart sound based on the matching pursuit method[END_REF]

Figure 5 .

 5 Figure 5. Accuracy scores for the experiments conducted. At the right chart data was oversampled with SMOTE while in the left is not.

Table 1 .

 1 Brief description of the tested classifiers in this work. Parameters tunning and references for details.

	Acronym Full name	Main parameters	Value of parameters	Reference
			Criterion to measure the
		Classification	quality of each split,	Gini criterion,
	CART	and regression	spitter, strategy to	best split between	Breiman et al 27
		tree	choose the split	trees.
			at each node
			Number of neighbors,
	KNN	K-nearest neighbors	weight function used metric of distance to in prediction,	5 neighbors, Minkowski's distance. uniform weight,	Wenberger et al 28
			use for the tree
			Name of solver to reduce
		Linear	data dimensions, tolerance	Solver: SVD,
	LDA	discriminant	value threshold for	tolerance of	Friedman et al 29
		analysis	rank estimation	0.0001
			in the solver
	LR	Logistic regression	Penalty norm, number of iterations	L2 norm penalty, 100 iterations	Schmidt et al 30
			Prior probabilities	No prior
	NB	Naive Bayes	in the classes to	probabilities	Chan et al 31
			adjust the data	value
			Number of estimators,
	RF	Random forest	criterion to measure the use or not bootstrap to build quality of the split,	100 estimators, bootstrap: true Gini criterion,	Breiman et al 32
			trees.
				C = 2.37 for
	SVM	Support vector machines	Penalty parameter C of the error term, kernel function used in the coefficient gamma algorithm, kernel	unbalanced data and C = 1 for automatic gamma kernel with balanced, radial basis function as	Cortes et al 33
				parameter

Table 2 .

 2 Performance metrics resulting for our heart sounds cross validation test without balancing.

	Model			Average performance
	Classifier Features set	SE	SP	ACC	MCC
	CART	A	50.19 ± 6.33	87.30 ± 1.91	0.79 ± 0.03	0.38 ± 0.08
	KNN	A	48.26 ± 4.85	86.36 ± 1.21	0.78 ± 0.02	0.35 ± 0.06
	LDA	A	58.32 ± 6.23	83.36 ± 0.86	0.81 ± 0.01	0.31 ± 0.05
	LR	A	58.44 ± 6.97	83.36 ± 0.95	0.81 ± 0.02	0.32 ± 0.06
	NB	A	49.56 ± 10.41	80.84 ± 0.68	0.79 ± 0.02	0.18 ± 0.06
	RF	A	75.80 ± 5.32	88.03 ± 1.32	0.86 ± 0.02 0.55 ± 0.06
	SVM	A	46.41 ± 1.99	92.29 ± 0.76	0.76 ± 0.02	0.45 ± 0.03
	CART	B	42.70 ± 3.75	84.76 ± 1.02	0.76 ± 0.02	0.28 ± 0.05
	KNN	B	60.63 ± 3.82	84.54 ± 0.85	0.82 ± 0.01 0.36 ± 0.04
	LDA	B	48.11 ± 10.78	80.65 ± 0.76	0.79 ± 0.01	0.17 ± 0.07
	LR	B	54.31 ± 12.73	80.43 ± 0.88	0.79 ± 0.01	0.17 ± 0.08
	NB	B	58.59 ± 12.43	80.93 ± 0.68	0.80 ± 0.01	0.21 ± 0.07
	RF	B	72.74 ± 8.88	83.17 ± 1.11	0.82 ± 0.02	0.36 ± 0.08
	SVM	B	34.95 ± 4.69	81.91 ± 0.90	0.73 ± 0.02	0.16 ± 0.05

Table 3 .

 3 Performance metrics resulting for our heart sounds cross validation test when using SMOTE balancing.

	Model			Average performance
	Classifier Features set	SE	SP	ACC	MCC
	CART	A	82.61 ± 1.60	84.98 ± 1.85	0.84 ± 0.02	0.68 ± 0.03
	KNN	A	70.74 ± 1.39	97.10 ± 1.63	0.79 ± 0.02	0.62 ± 0.03
	LDA	A	76.95 ± 1.27	82.37 ± 2.14	0.79 ± 0.01	0.59 ± 0.03
	LR	A	78.60 ± 1.60	82.38 ± 2.50	0.80 ± 0.02	0.61 ± 0.04
	NB	A	84.54 ± 4.58	53.29 ± 0.63	0.56 ± 0.01	0.21 ± 0.03
	RF	A	91.60 ± 1.77	92.10 ± 1.70	0.92 ± 0.01 0.84 ± 0.03
	SVM	A	77.20 ± 1.09	78.77 ± 1.99	0.78 ± 0.01	0.56 ± 0.03
	CART	B	80.71 ± 2.43	82.91 ± 1.82	0.82 ± 0.02	0.64 ± 0.03
	KNN	B	79.24 ± 1.56	95.87 ± 0.79	0.86 ± 0.01	0.73 ± 0.02
	LDA	B	70.55 ± 1.52	70.89 ± 2.44	0.71 ± 0.02	0.41 ± 0.04
	LR	B	69.90 ± 1.79	67.97 ± 1.94	0.69 ± 0.02	0.38 ± 0.04
	NB	B	82.80 ± 5.09	52.78 ± 0.66	0.55 ± 0.01	0.19 ± 0.03
	RF	B	89.10 ± 1.46	91.55 ± 1.40	0.90 ± 0.01 0.81 ± 0.02
	SVM	B	67.18 ± 4.03	53.41 ± 0.96	0.56 ± 0.02	0.15 ± 0.04

ACKNOWLEDGMENTS

The first author would like to express his gratitude to the Rennes Métropole department, which funded its stay at Inria Rennes, France research center to support this work. M. Ibarra is a PhD student supported by the Mexican National Council for Science and Technology (CONACYT) through the Graduate Research Fellowship No. 477876. Authors also want to thank to the sponsors of the Physionet/Cinc 2016 Challenge for providing the largest open database of heart sounds to the scientific community.


Corentin Vallée

Pierre Maurel

Isabelle Corouge

Christian Barillot

Resting State ASL : Toward an optimal sequence duration

Resting-state functional Arterial Spin Labeling (rs-fASL) in clinical daily practice and academic research stay discreet compared to restingstate BOLD. However, by giving direct access to cerebral blood flow maps, rs-fASL leads to significant clinical subject scaled application as CBF can be considered as a biomarker in common neuropathology. Our work here focuses on the link between overall quality of rs-fASL and duration of acquisition. To this end, we consider subject self-Default Mode Network (DMN), and assess DMN quality depletion compared to a gold standard DMN depending on the duration of acquisition.

Introduction

In both resting-state functional MRI (rs-fMRI) academic research and clinical applications, Arterial Spin Labeling (ASL) remains rare compared to Blood Oxygenation Level Dependent (BOLD) fMRI. Despite ASL lower estimated signal-to-noise ratio compared to BOLD-fMRI, rs-fASL is a promising valuable technique to be used in clinics since :

It provides cerebral blood flow (CBF), and thus leads to potential significant clinical applications as CBF can be used as a biomarker in various neuropathology (Alzeihmer disease, severe depression...). For over the last decade ASL processing has rapidly improved, and becomes more and more disseminated for clinical usage.

As asked in : "What ASL time series length would be adequate?", our purpose is to provide a first answer to this question. While some indications were given for rs-BOLD , this work has not been done yet for rs-fASL, as far as we know. In this abstract, we will focus on one functional brain network characterizing the resting-state: the Default Mode Network (DMN). Given a duration of acquisition (DOA), we make the assumption that the quality of the DMN defines quite well the quality of the corresponding rs-fASL acquisition. As we want to work on a subject-scaled (first-level) analysis for further applications and as a preliminary work, we will describe the variation of quality of the subjects self-DMN compared to a gold standard reference.

Methods

We considered a dataset of 6 healthy subjects who volunteered for a 24min40s (426 volumes) rs-fASL acquisition. We split each acquisition into 20 parts of various durations (5min to 24min by 1). Each of the 6x20 acquisitions were processed independently with typical preprocessing . In order to work at subject-scaled level (i.e. considering each subject DMN, not group built DMN), we relied on Seed-Based Analysis (SBA). We selected seeds in the medial prefrontal cortex as they give the best results, according to the number of non-noised DMN exhibited among the 120 acquisitions. We built maps of the second order approximation of the FDCR estimator of correlation as we can get rid of the bias of the canonical estimator of correlation considering the given number of volumes (86 for 5min up to 412 for 24min). Then, we tested whether each voxels signal is positively correlated with the seed (risk FWER-corrected set at 0.01). Let's the binary mask obtained be called with the subject . We wanted to compute the Jaccard index between the 120 bmaps with a gold standard reference, let's call it . However, taking a subject's DMN as gold standard leads to an overwhelming dependency between bmaps and the reference. Moreover, we can not build a global this way as would depends on . Hence, we took as gold standard the DMN from the recent atlas of functional area of the brain (MSDL) proposed by Varoquaux et al. . Now, for every subject , we can compute and as MDSL DMN is not subject dependent, we can also compute a global . Finally, we modelled the Jaccard index for each subject and for all subjects with the non-parametric local regression LOESS with two degrees local polynomials and a span equal to 0.8.

Results

Every

and its associated LOESS is shown in Figure 1. Some corresponding overlays of bmaps and MSDL DMN on MNI152 template are shown in Figure 2 to provide qualitative illustrations of Jaccard Index. The global and its associated LOESS is presented in Figure 3.

All

increase before around 16-18min. 1. After a peak around 17 min each except and decrease. 2. For , , , close DOA does not imply close Jaccard Index. 3. Global increases till 17min. 4. After 17min decrease. 5. Considering the LOESS on global , on average, has 95% chance to be higher at 17 minutes than before 12 minutes and after 21 minutes. 6.

Discussion

Surprisingly, the quality of the DMN looks to be not maximal when DOA is maximal. In our dataset, the optimal quality of rs-fASL seems to be reached around 17 minutes. Result (3.) suggests quality may not be stable over time in subject-scaled analysis.

Conclusion

In this preliminary work, we exhibit a subject-scaled and an overall tendency about quality of rs-fASL, depending on DOA (under the assumption of correlation between quality of the DMN and overall quality of the rs-fASL acquisition). For subject-scaled analysis, checking subset of a whole acquisition should be considered as quality seems not to be always stable. In further work, we would like to confirm the tendency by evaluating other usual functional areas of interest like motor, visual or language.

Figures

  Figures

Figure 1 :

 1 Figure 1 : Curves of the computed LOESS for each subject. The grey area is the 95% confidence area related to the LOESS. Framed dots are corresponding to overlays shown in Figure 2. Colour on subjects as well as aspect ratio are matching those in Figure 3.

Figure 2 :

 2 Figure 2 : Some illustrations of values of the Jaccard Index according to the DOA with corresponding maps. All images shown are registered in the MNI 152 template. Red areas correspond to the MSDL DMN. Blue areas correspond to subject voxels whose signal were significantly positively correlated to the chosen seed signal in medial prefrontal cortex.

Figure 3 :

 3 Figure3: Curve of the computed LOESS for all subject. The grey area is the 95% confidence area related to the LOESS. Framed dots are corresponding to overlays shown in Figure2. Colour on subjects as well as aspect ratio are matching those in Figure1. Black lines show the average values of Jaccard index according to the LOESS at 95% chance to be higher than before 12 minutes and after 21 minutes.

Acknowledgements


Christine Ménager
email: menager@ccr.jussieu.fr

Olivier Sandre

Jêrome Mangili

Valérie Cabuil

Jérôme Mangili

Preparation and swelling of hydrophilic magnetic microgels

Keywords: magnetic fluid, microgel, hybrid particles

HAL is a

1.Introduction

The synthesis of hybrid materials by encapsulation of an inorganic material or a liquid in a polymer matrix opens new possibilities in the fabrication of solid particles for materials science. Such composite particles offer potential advantages in applications such as cosmetics, inks and paints coming from the improved compatibility between the filler and the binder [1].

Polymer shells encapsulating inorganic particles are also of great interest in pharmaceutical and biotechnological industries for producing drug release products [2]. However, recent studies have underlined the major difficulty to encapsulate inorganic particles homogeneously, especially in the case of magnetic particles in polymer microgels [3,4].

In the field of drug delivery, many types of colloidal carriers have been tried out, which actively respond to external stimuli such as temperature, pH, electric or magnetic field [5][6][7].

Among the good candidates, one finds microgels that are cross-linked latex particles swollen by a good solvent.

Aqueous magnetic fluids are colloidal dispersions of magnetic nanoparticles in water, exhibiting a giant paramagnetic behavior. They can be used in biosciences, namely as contrast agents for magnetic resonance imaging (MRI) [8], for magnetic guidance of drugs or radioisotopes and for cell sorting processes.

The encapsulation of magnetic particles inside organic swellable particles made of polystyrene or polymethylmethacrylate has been extensively described [9][10]. In this paper, we describe a general method to encapsulate magnetic particles homogeneously inside a hydrophilic gel or hydrogel. The usual way to produce magnetic polymer particles consists in coating the magnetic material by a surfactant and embedding it in the polymer using processes such as suspension, emulsion or precipitation polymerization. Here we use the inverse emulsion process to disperse the aqueous magnetic fluid and the hydrophilic monomer localized together in water droplets just prior to polymerization. As the magnetic particles are themselves hydrophilic, no surface treatment is required before the encapsulation process. This is a key factor since there is neither surfactant nor polymer molecules in the initial magnetic fluid that can interfere with the surfactant used for the emulsification or with the polymer matrix.

Marerials and methods

Materials

Acrylamide 99% (AM) was purchased from Aldrich, N,N'-methylenebisacrylamide (MBA) was obtained from Sigma, ferric chloride (FeCl 3 27% solution), ferrous chloride (FeCl 2, 4H 2 O), ammonia (NH 3 20%), nitric acid (HNO 3 52.5 %), hydrochloric acid (HCl 35%), trisodium citrate (C 6 H 5 Na 3 O 7 ,2H 2 O) and sorbitan monooleate (Span 80) were purchased from Prolabo. V-50 (2-2' azobisamidino propane) was obtained from Fluka.

Synthesis of the aqueous magnetic fluid :

180 g of Fe II salt (0.9 mole) and 715 mL of Fe III solution (1.5 moles) were mixed in acidic medium (HCl). Coprecipitation of magnetite (Fe 3 O 4 ) particles was performed by rapid addition of concentrated ammonia (1L). After acidification with a nitric acid solution (360 mL), the addition of 800 mL of a boiling solution of ferric nitrate (323 g, 1.3 moles) leaded to the chemical oxidation of magnetite to maghemite (γFe 2 O 3 ). After magnetic decantation, the red precipitate was dispersed in water and then 70 g of trisodium citrate salt (0.24 mole) was added. The mixture was heated at 80°C for 1/2h. The magnetic material was precipitated by addition of acetone. The precipitate was isolated and suspended in distilled water leading to a stable colloidal dispersion (pH7).

Emulsion polymerization

Two series of micron-sized and sub-micron sized microgels were synthesized, both of them being based on polyacrylamide cross-linked by N,N'-methylenebisacrylamide (MBA).

Both microgels were synthesized using an inverse (W/O) emulsion polymerization method. Aqueous droplets were stabilized by an emulsifier (Span 80) and dissolved in dodecane. In a first step, the polar monomer and the cross-linker were solubilized directly in the magnetic fluid to be incorporated in the final microgels. In a second step, this aqueous mixture was introduced under stirring into the oily phase. Polymerization was started by the introduction of a free radical initiator (V50). Using either inverse emulsions or miniemulsions, it is possible to investigate two ranges of sizes : the micron or the sub-micron one respectively.

Polymerization in emulsion:

We investigate here two weight ratios MBA/AM of crosslinker with respect to the monomer (5 wt% and 2.5 wt%) at a constant monomer concentration (28 wt% AM in the aqueous magnetic phase). The hydrophilic monomer acrylamide (AM, 2g) and the cross-linker N,N'-methylenebisacrylamide (MBA, 0.1g and 0.05g respectively for the two investigated weight fractions were solubilized in 5 mL of the pH7 magnetic fluid. Both aqueous phase and organic solution of the emulsifier (Span 80 in 15mL dodecane, m Span80 /m dod = 2.25%) were placed under N 2 for ½ h and then added together and stirred. Polymerization was initiated by adding 20 mg of V-50 solubilized in 0.5 mL of water (1wt% /AM) and raising temperature to 50°C in a water bath. The reaction mixture was stirred one hour at 50°C under N 2 . At the end of the reaction, the mixture had become heterogeneous: the magnetic microgels were settling at the bottom of the flask, the biggest ones being often stuck together. However, the magnetic microgels were easily separated from the organic solvent under a magnetic field gradient (a strong Na-Fe-B magnet) and washed several times with ethanol. They were then redispersed in water but due to their size, such suspensions of magnetic microgels in water were not stable.

Polymerization in miniemulsion

: the magnetic miniemulsions were prepared under the same conditions as the emulsions except for mixing, which differs by the use of ultra-sons.

After degassing, the organic phase (0.3 g Span 80 and 11 g dodecane) was mixed with the aqueous solution (1.4 g AM, 0.07 g MBA in 1.9 mL of the aqueous magnetic fluid). The radical initiator (20 mg V50 in 0.5 mL of water) was introduced at the end (after nitrogen bubbling). After vortex stirring (5 minutes) the mixture was placed in an ultrasonic bath (Elma Transsonic T 460H, 285 Watts, 35 kHz). The reaction mixture was sonicated at 50°C for 15 minutes. After reaction, we always noted the presence of a coagulum and a supernatant.

The microgels in the supernatant were characterized in their solvent of synthesis and also after several washings in ethanol and finally in water.

Determination of the magnetic load

The mass fraction of magnetic material inside the composite microgels was obtained from the ratio of the iron oxide weight divided by the total dry weight of a given sample. In the case of the stable sub-micron microgels, magnetic measurements could be performed with a vibrating sample magnetometer because the suspension did not settle down under the strong applied magnetic field (up to 1 tesla).

Concerning the larger micron-sized microgels, the molar content of iron was titrated by atomic absorption spectroscopy with a Pelkin Elmer Analyst 100 apparatus after degradation of the particles in boiling HCl (6 mol/L).

Optical microscopy

The samples were observed with an inverted microscope (Leica DMIL, 40x magnification).

The orange color of the microgels is due to the trapped magnetic nanoparticles that strongly absorb in the visible spectrum. The pictures recorded with a color video camera are digitized with a PC graphic card (ATI All-In-Wonder Pro).

Transmission Electron Microscopy (TEM)

TEM was performed using a JEOL 100CXII top entry UHR electron microscope. The diluted solutions were deposited on a carbon-coated copper grid and left to dry in the air. No further contrasting technique was applied. In the case of miniemulsions only, the microgels were deposited in their solvent and then washed directly on the grid with ethanol.

Multi-angle Photon Correlation Spectroscopy (PCS)

PCS was performed on a Macrotron spectrometer (Amtec) equipped with a He/Ne laser (λ=632.8 nm) and a Malvern 7132 correlator card with 192 time channels. Samples (about 1 mL in cylindrical glass vials) were immersed in a filtered toluene bath, which temperature was regulated (±0.2°C) by a water jacket. The autocorrelation functions of the scattered intensity were recorded using a 0.5 µs sampling time and 30 sub-runs of 20 s averaged together for a total duration of 10 min. The scattering angle θ was varied from 45° up to 155° every 5°. Precise measurements of the hydrodynamic diameter by PCS necessitate the exact knowledge of the solvent viscosity (η solvent ) and of the refraction (n solvent ). Therefore PCS could not be performed directly in the synthesis emulsion which contains many species (polymer, surfactant, …) which modify those two parameters. Indeed the particles were isolated with a strong permanent magnet, washed several times and resuspended in filtered absolute ethanol (η EtOH =1.159 cPoise at 22°C and n EtOH =1.3611) or freshly distilled water (η water =0.957 cPoise at 22°C and n water =1.3326). Those values of dynamic viscosities were calculated from the solvent densities (ρ water =1, ρ EtOH =0.7937) and experimental values of the kinematic viscosities (ν water =0.957 cStokes, ν EtOH =1.461 cStokes at 22°C) measured with an Ubbelohde viscosimeter (capillary No. I, Schott-Geräte, Germany). Tabulated values were used for the indexes of refraction. The autocorrelation functions were analyzed using the standard program CONTIN [11][12], in order to extract the dominant decay time τ max . Then the decay rate Γ max = (τ max ) -1 was used in the primary equation of PCS theory:

solvent 2 B max 3 η π q d T k Γ H =
where the scattering vector writes:

      = 2 sin 4 solvent θ λ π n q
The hydrodynamic diameters d H were computed from the slope of a linear regression of Γ max as a function of q 2 / η solvent extrapolated in the zero q limit. The reliability of this method was tested on monodisperse beads of (Stöber) silica with d TEM =278 nm, for which the linear regression yields d H =288 nm (the discrepancy of 10 nm can be attributed to a hydration layer).

Analyzing the auto-correlation functions by the method of cumulants rather with the CONTIN program decreases the quality of the linear regressions of Γ max vs. q 2 / η solvent , presumably due to the influence of a secondary characteristic time (see inset in Figure 5a). This mode approximately ten times faster than the main one corresponds to a hydrodynamic diameter about 15 nm, therefore attributed to Brownian diffusion of traces of magnetic nanoparticles released from the gel into the solution.

This mode is ten times faster than the main one and therefore attributed to Brownian diffusion of traces of magnetic nanoparticles released from the gel into the solution.

Results and discussion

All related studies revealed that a successful encapsulation first requires a good dispersion of the inorganic particles in the medium prior to carrying out the polymerization [13]. In fact, this is a prerequisite to their use in a variety of applications. The magnetic load must be uniformly dispersed in the gel. The clustering of the magnetic nanoparticles in a localized spot inside the composite colloids has often been observed in the literature [3].

Compatibility between the magnetic fluid and the polymer

The encapsulated magnetic nanoparticles are maghemite (γFe 2 O 3 ) grains, which were synthesized using a method described above and elsewhere [14]. The diameters of particle range between 5 and 14 nm. This rather broad size distribution can be narrowed through a post synthesis treatment [15,16]. In order to avoid aggregation of the ferric oxide particles in the polymer network, we used a magnetic fluid that is stabilized in water at pH7. The trapping of the same kind of particles in a polymer network has already been described in the case of bulk (centimeter large) hydrogels [17,18]. These magnetic nanoparticles are covered by citrate ligands and are negatively charged at pH7. The electrostatic repulsions ensure the stability of the colloidal dispersion. The control of the stability of the inorganic particles dispersions during emulsification and polymerization appears to be the major issue [19]. The volume fraction of magnetic particles set to 2.3% and a concentration of unadsorbed citrate species in the aqueous phase equal to 2.5 10 -2 mol.L -1 are appropriate to avoid any microscopic destabilization of the colloidal dispersion after mixing with the monomer. The volume fraction of magnetic particles was determined by magnetic measurements. The magnetization curve was measured with a vibrating sample magnetometer [20]. A fit of this curve allows the determination of the volume fraction of maghemite nanoparticles (φ=2.3% from the magnetization at saturation and of the size distribution modeled by a Log-normal law (d 0 =7nm, σ=0.4).

Depending on technical applications, it is also necessary to control the size of the microgels.

For this purpose, two kinds of synthesis are described here, leading to final microgels either in the micrometer or in the sub-micrometer range of size.

Micron-sized microgels

Magnetic behavior and magnetic content

Due to their size, the microgels prepared in emulsions and dispersed in water can be observed using optical microscopy (Figure 1a). The microgels are spherical just like the emulsion droplets and their diameters vary between 10 and 60 µm. When a magnetic field of low intensity (200 Gauss) was applied, we observed the chaining of the magnetic microgels with perfect reversibility (Figure 1b). Please note that there was no visible deformation of the sphere (magneto-striction) at least for this low field intensity.

The determination of the magnetic load inside those microgels yielded 64 wt% for the larger cross-linking density (5%wt MBA/AM) and 34 wt% for the lower one (2.5%wt MBA/AM).

Swelling degree and dynamics

Microgels swell in suitable solvents and their swelling degree depends on the concentration of cross-linker. We followed the swelling of microgels (5 and 2.5 wt% MBA/AM) individually under the microscope. The microgels were deposited from an alcoholic suspension onto a glass slide; after evaporation of ethanol, the beads of collapsed gel were allowed to hydrate by adding a water droplet on them. The swelling process was recorded and shown on Figure 2.

Thanks to the coloration by iron oxide, we observe the flow of solvent towards the center of the microgel concomitant with the increase of its external diameter. The flow of solvent (water) inside the microgel is concomitant with the increase of its diameter, as seen on the plot of the diffusion front (closed markers) and of the external diameter (open markers) as a function of time (Figure 2). The degree of cross-linking affects the swelling ratio significantly: qualitatively a 2.5 wt% cross-linked microgel absorbs about twice the volume of water of a 5 wt% MBA microgel for the same total monomer concentration.

These experiments to follow kinetics of minigel swelling are reminiscent of those reported for polyacrylamide spherical macrogels with diameters of the order of the millimeter [21].

Interestingly, the characteristic time of swelling τ is about 10 4 -10 5 shorter in the case of minigels, in accordance with the theory of Tanaka et al. that shows that τ scales with the square of the bead diameter. On these pictures, we also observe the transient rippling instability of hydrogel surface, which was already reported for centimeter large hydrogels with a surface constrained to adhere on a substrate and a free interface with the swelling solvent [22].

Please note that there is no detectable release of the magnetic content, this reversible swelling resulting only from the diffusion of water and ions from the outer solution towards the porous medium constituting the hydrogel.

TEM analysis

Due to their inorganic load, the magnetic microgels can be analyzed by TEM without any coloration, the electronic contrast being already provided by the magnetic particles (Figure 3a). The magnetic particles are well dispersed in the gel, the darker color at the center compared to periphery being due to a larger thickness. Figure 3b shows a magnification of the microgel sphere at the periphery. Although the density of the film is high, we can still visualize individual particles, which means that they are not distributed as islands. In addition, the saw-tooth profile of the gel wedge is attributed to a wetting instability particular to elastic microgels that will be developed in a forth-coming paper.

Sub-micron sized microgels

Ugelstad et al [23] first published results on miniemulsions polymerization leading to polystyrene polymer particles. The principles of miniemulsion polymerization in non-aqueous systems (inverse miniemulsions) have been described to produce hydrophilic polymer latexes [10]. It has been shown that the addition of salt has a positive effect on the stability of the inverse miniemulsion by acting as osmotic pressure agent preventing Ostwald ripening. Here the salt was added with the magnetic fluid (tri-sodium citrate salt).

Magnetic content

The magnetic content has been calculated from the magnetic measurement made on stable colloidal dispersion of magnetic microgels in water and from the dried mass of polymer. For a sample with an overall solid fraction (dry weight /total weight) of 0.77 wt%, we measured a magnetization at saturation M s =330 A/m (Figure 4). Using the specific magnetization of colloidal maghemite m s =3.1 10 5 A/m which is experimentally about 75 % of the value for bulk maghemite [24], we obtain the volume fraction of magnetic particles φ=M s /m s =0.106% and thus their mass fraction 0.53 wt% (ρ = 5 g.cm -3 ). The magnetic mass fraction is thus 0.53 / 0.77 = 69 wt% in these microgels.

Swelling degree

Contrary to samples polymerized in conventional emulsions, the miniemulsions were stable (remained dispersed) several weeks after polymerization in their solvent (dodecane). We did not check all the criterions to ensure that the presented system could really be qualified as a polymerization in miniemulsion, but it corresponds to miniemulsion in terms of stability and of particles size. Polymerization in miniemulsions with 5 wt% MBA/AM and 42 wt% AM monomer in the aqueous magnetic phase allowed to get magnetic microgels that exhibit colloidal stability.

By performing multi-angle photon correlation spectroscopy (PCS) on dispersions of microgels in bad and good solvent respectively, we were able to demonstrate a change in size in response to a change in solvency. By plotting the dominant correlation rate Γ max versus the square of the scattering vector divided by the solvent viscosity, q 2 / η solvent , we obtain the hydrodynamic diameter from the inverse of the slope near the origin: d H EtOH =158 nm, d H water =195 nm and d H =288 nm for the test beads (Figure 5a). As a result, multi-angle PCS measurements and extrapolation to the low-regime enabled to bypass the issue of polydispersity, thereby reporting a significant increase (+23.4%) of the hydrodynamic diameters (Figure 5a). The microgels thus undergo (on average) a (1.234) 3 ~ 1.9-fold increase in volume when they are transferred from a bad solvent (ethanol) into a good solvent (water) [25]. Please note that this measurement of the hydrodynamic diameter represents the (volume) averaged value among a large distribution of diameters. Unlike the straight line over the whole range of q vectors obtained with monodisperse silica beads, the curves on Figure 4a for minigels both in ethanol and water exhibit a positive deviation from linearity at large q. This upward curvature is attributed to a large polydispersity of the samples, estimated as a Lognormal distribution with σ=0.45 from the width at half-height of the CONTIN peak.

TEM analysis

Figure 5b is a TEM picture of microgels synthesized in dodecane. Most of them appear spherical and the measurement of diameters (d TEM ) over 82 minigels shows a Log-normal law of peak value d 0 =65 nm and standard deviation σ=0.25 (no unit).

Interestingly, all the microgels are loaded with magnetic nanoparticles, which appear homogeneously distributed within the polymer.

Conclusion

Rapid and reversible response to a weak magnetic field, solvent swelling without release of the magnetic content and homogeneous distribution of the magnetic colloid inside the matrix, these are the major properties of the hydrophilic magnetic microgels described here. All these properties concomitant with a high magnetic oxide content are essential for biotechnological applications. Here we have described a simple method of homogeneous encapsulation of inorganic particles in hydrogel microbeads. The monomer was just added to the colloidal solution made of inorganic particles in water avoiding any destabilization of the colloidal solution by the careful adjustment of the ionic strength and the pH. However the homogenization of the miniemulsion could be improved by means of turbidity and surface tension measurements as described by Landfester [9].

The process described in this paper is widely applicable to a large variety of water-soluble polymers (polyNIPAM for instance) and to inorganic materials that can be dispersed as colloidal solutions. Moreover, it allows the co-encapsulation of the magnetic fluid with water soluble drugs or markers for biomedical applications.

Figure 1 :

 1 Figure 1: Optical micrographs of magnetic microgels (5%wt MBA/AM) prepared in water in

Figure 2 :

 2 Figure 2: Optical micrographs of the swelling dynamics of magnetic microgels (5% wt

Figure 3 :

 3 Figure 3: a) TEM image of a microgel (5%wt MBA/AM) prepared in emulsion. The microgel

Figure 4 :

 4 Figure 4: Magnetization curve of sub-micron sized microgels measured with a vibrating

Figure 5 :

 5 Figure 5: a) PCS measurements on magnetic minigels prepared in a water/dodecane

Figures

  Figures

Figure 1

 1 Figure 1

Figure

  Figure 2

  Figure 2

Figure 3 Figure 4

 34 Figure 3

Figure 5 18

 5 Figure 5

Acknowledgments

The authors thank Chamsoudine Youssouf for preliminary experiments. The electron microscopy was performed by Fabienne Warmont of the "Service de microscopie électronique de l'Université Paris 6".


Tele-MAGMaS: an Aerial-Ground Co-Manipulator System Nicolas Staub 1 , Mostafa Mohammadi 2,3 , Davide Bicego 1 , Quentin Delamare 4 , Hyunsoo Yang 6 , Domenico Prattichizzo 2,3 , Paolo Robuffo Giordano 5 , Dongjun Lee 6 and Antonio Franchi 1

I. INTRODUCTION

The manipulation of large size objects by robotic systems is a challenge for applications in the field of construction industry, industrial decommissioning, and Urban Search and Rescue (USAR). These are associated to dangerous environments and thus motivate the need of devising robotic solutions for replacing human presence. Furthermore, they often require manipulation of long objects, such as pipes, bars, beams and metal frameworks, with limited access to their center of mass (CoM).

Robotic object handling has been tackled using (mobile) ground manipulators (GM). The rich literature on GMs proposes solutions with single or multiple robots for, e.g., cooperative transportation of large objects [START_REF] Machado | Multi-constrained joint transportation tasks by teams of autonomous mobile robots using a dynamical systems approach[END_REF], offshore robotic sensing and manipulation [START_REF] Pfeiffer | Offshore robotics -survey, implementation, outlook[END_REF], or for cooperative assembly [START_REF] Knepper | Ikeabot: An autonomous multi-robot coordinated furniture assembly system[END_REF]. The use of GMs has two major drawbacks: first, typical small industrial manipulators have limited joint torques, resulting in poor maximal admissible Cartesian torques at the end-effector (EE) and second, GMs have a rather small workspace around their base, reducing their manipulation capabilities. These limitations can be particularly problematic for long objects if the GM cannot grasp them by their CoM, as that demands high torques at the EE and large workspace.

An emerging approach is the use of aerial manipulators (AM) for construction and large load handling. This has been made possible by the recent developments in Unmanned Aerial Vehicles (UAVs) control, especially the latest progress in the field of physical interaction, see e.g., [START_REF] Ruggiero | Aerial manipulation: A literature review[END_REF] for a panorama of the field. Their use has been demonstrated, in cooperative load transportation using cables [START_REF] Manubens | Motion planning for 6-D manipulation with aerial towed-cable systems[END_REF], [START_REF] Sreenath | Dynamics, control and planning for cooperative manipulation of payloads suspended by cables from multiple quadrotor robots[END_REF] (splitting the overall payload among members), in multi-robot assembly [START_REF] Augugliaro | The flight assembled architecture installation: Cooperative construction with flying machines[END_REF], [START_REF] Lindsey | Construction with quadrotor teams[END_REF], and aerial manipulation [START_REF] Kim | Aerial manipulation using a quadrotor with a two DOF robotic arm[END_REF]- [START_REF] Mohammadi | Cooperative aerial tele-manipulation with haptic feedback[END_REF]. A major drawback of these platforms is their limited payload.

In this paper, we present a novel class of heterogenous systems which tackles the problem of manipulating long objects that cannot be grasped close to their CoM. Such systems go beyond the limitations of the previous approaches by leveraging the advantages of both AMs and GMs together. The small payloads of the AMs are compensated by the strength of the GM, while the limited workspace and poor Cartesian torque at EE of the GM is balanced by the virtually unlimited workspace and the favorable lever provided by AMs. The AM act as Flying Assistants, as they assist the GM from the air. The proposed class is called Tele-MAGMaS, where MAGMaS stands for 'Multi Aerial Ground Manipulator System' and Tele reflects the tele-operation capabilities, in fact, such systems also allow for remote human operation at different autonomy levels. The presence of a human operator, who owns superior intelligence and cognitive capabilities, is necessary in most of the realistic applications to cope with unknown/partially known environments with possibly unpredicted events. In the proposed Tele-MAGMaS scheme the human operators are also provided with haptic feedback to enrich their tele-presence and improve their performances.

Loosely related to the Tele-MAGMaS concept, a theoretical scheme for cooperative transportation between an AM and a mobile cart has been presented for the planar case in [START_REF] Nguyen | Control of a UAV and a UGV cooperating to manipulate an object[END_REF]. Another scheme with GM and tethered UAVs is presented in simulation in [START_REF] Lippi | Cooperative object transportation by multiple ground and aerial vehicles: modeling and planning[END_REF]. The Tele-MAGMaS concept considered in this work consists, instead, of two real world manipulators -aerial and ground, respectively, and with any number of degrees of freedom (DoF) -which have to perform a real manipulation task.

This paper highlights the key components of Tele-MAGMaS and demonstrates practically the implementation of the described algorithm in software and hardware for aerial-ground co-manipulation. Theoretical groundings are discussed in our previous works [START_REF] Staub | Towards robotic MAGMaS: Multiple aerial-ground manipulator systems[END_REF], [START_REF] Staub | Towards a flying assistant paradigm: the OTHex[END_REF]. To the best of our knowledge, this is the first time a Flying Assistant has been implemented. In particular, for the first time a cooperative manipulation task between a ground industrial manipulator and an aerial manipulator has been robustly demonstrated.

The design, architectural and experimental aspects raise many scientific and technological challenges. Indeed, a cooperative manipulation implies the exchanges of forces which need to be handled carefully to ensure system stability. Additionally experimental validation asks for robustness to model imperfections and over non-linearities. Based on our experience we felt the necessity to develop our own fully-actuated AM, able to resist external forces disturbances without changing its orientation and with an extra aperture for object manipulation. The proposed system is over-actuated w.r.t the task and the leveraging of redundancy in the system is not straightforward. In particular, attention to all actuation constraints is necessary during the cooperative manipulation. Also, the integration of the human operator needs to be handled with care, so that no provided commands can push the system beyond its limits and in unstable configurations. This requires careful analysis of the system behavior and the construction of safety margins around the operator behavior. The high complexity of the Tele-MAGMaS necessitated to develop software in a modular way in order to allow easier development, this required the careful design of the full software architecture used in this work, which implements a control framework allowing for the different modalities of i) full autonomy, ii) tele-operation and iii) shared control of the system, thus allowing the system to cope with the different complexity levels of various environments by also leveraging (when needed) the cognitive abilities of a human operator.

II. ROBOTIC SUB-SYSTEMS

In this section, we present the three robotic systems that compose the Tele-MAGMaS, the complete system depiction can be found in the Technical Report attached in the multimedia material. The GM is a state-of-art industrial manipulator, the AM is a custom platform, called OTHex, specifically designed for bar lifting. Lastly, the haptic interface is a commercial delta manipulator, with force feedback capabilities. The three robotic systems have all been chosen following the principles exposed in the Technical Report attached in the multimedia material.

A. Ground Manipulator -LBR-iiwa

As GM we used the 7-DoF LBR-iiwa 14 R820, allowing a 14kg payload and 820 mm reach (see Tab. I for more data), this satisfies the task requirements for the GM in terms of payload and workspace, see the attached Technical Report. The embedded joint torque sensors allow to retrieve the external wrench without the need to integrate a F/T sensor at the EE. Despite the significant payload, its maximum admissible joint torques for the last joins is 40 N and its rated torque of 23.3 N m corresponds to holding horizontal a 1.42 kg bar of length 3 m, which may not reflect typical bar-weight for several of the considered applications.

The software architecture of the LBR-iiwa is dictated by the proprietary solution. The control cabinet contains two computers, a Windows machine for the high-level interface and a Vx-Work machine for real-time execution of the applications. In our setup, we use a KUKA library, called FRI for Fast Robot Interface, which uses UDP to achieve real-time control from a remote computer. The application on the control cabinet runs a server, while also running the low-level safety features, to which the client sends commands at high frequency (up to 1kHz) and receives system information. In this modality, the remote computer has access to various quantities in order to close the control loop and sends joint commands to a low-level joint control loop that runs on the control cabinet. This allows real-time interaction between the operator and the manipulator. Any other manipulator would have worked as long as hard real-time remote control is supported.

During our implementation only the joint controllers (position or impedance) were available in the FRI library. Hence we built our control architecture on top of the joint position controller, by implementing both an admittance filter and an inverse kinematics layer. For compliance, admittance filter is preferred over impedance control despite being more tedious to tune, because of the availability of low-level position controller, see [START_REF] Ott | Unified impedance and admittance control[END_REF]. That means that the new, compliant, reference trajectory for the EE pose is decoded into joint trajectory based on closed-loop inverse kinematics [START_REF] Antonelli | Stability analysis for prioritized closed-loop inverse kinematic algorithms for redundant robotic systems[END_REF].

As for the gripper, we integrated a Schunk PGN-plus-E 80-1 gripper on the flange of the LBR-iiwa, with custom designed claws. The choice of this gripper was made considering the significant torques that could be involved at the EE during our experiments (as pointed out in the Technical Report).

B. Aerial Manipulator -OTHex

The Open Tilted Hexarotor (OTHex) is an aerial manipulator specifically developed for this project and tailored to perform physical interaction tasks with the environment, i.e., meeting the requirements presented in the Technical Report.

Description

Its main features and peculiarities are presented here, more details can be found in [START_REF] Staub | Towards a flying assistant paradigm: the OTHex[END_REF], while its key specifications are listed in Table II. The first feature is the tilted propeller actuation. In classical multi-rotor design all the propellers are collinear, which leads to an underactuated system dynamics, i.e., position and orientation cannot be controlled independently. The OTHex design enforces non-collinear propeller orientation, as illustrated in Fig. 2-top. This design feature allows the system to have multidirectional thrust ability, i.e., the total thrust can be exerted in a polytope attached to the body frame instead of just along a single upward direction as in the underactuated case. This means that the robot can exert lateral forces without the need for re-orienting itself thus being able to track a decoupled reference trajectory in position and orientation, within the physical limits of the actuators. Such a design is gaining popularity in the literature [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally-bounded input force[END_REF]- [START_REF] Nguyen | A novel robotic platform for aerial manipulation using quadrotors as rotating thrust generators[END_REF] mainly because the multidirectional thrust aerial vehicles are, by design, particularly well suited for physical interaction tasks. Indeed, they can exert a decoupled set of forces and torques on their environment, independently of the contact point position. In particular, this means that they can resist external disturbances while completing a manipulation task without the need to change their orientation, which could jeopardize the manipulation task itself.

The second feature of the OTHex consists of an aperture of 85 between two bars of its skeleton (see bottom Fig. 2). This configuration has been preferred over the regular hexagonal positioning of the actuation units to facilitate the manipulation of a long object, e.g., a beam. Thanks to this unique design the object can pass through the aperture, allowing for a wider variety of beam manipulation tasks.

The extra payload of the OTHex AM can be used to assist the GM in the cooperative manipulation task. For example only 7 N at the end of the bar (0.25% of its extra payload) are sufficient to cooperatively transport a 3 m bar of weight that fits the considered applications.

Lastly we designed a lightweight mechanical system consisting of a passive revolute joint with two grippers, which endows the OTHex with grasping and manipulation capabilities required to perform physical interaction. A passive joint has been preferred over an actuated one to reduce the complexity and weight of the system; this solution is viable thanks to the full-actuation of the OTHex. Indeed, in the studied use case, the OTHex acts as a Flying Assistant following the GM motion, hence a 1-DoF passive joint is sufficient to accomplish the task.

The OTHex control architecture is articulated over three main components: low-level controller, wrench estimator and admittance filter, whose details can be found in our previous work [START_REF] Staub | Towards a flying assistant paradigm: the OTHex[END_REF] and therein.

To guarantee robustness of aerial physical interaction, software compliance was necessary. We preferred admittance over impedance framework, for two main reasons: i) the availability of a precise pose controller for the OTHex and ii) the ease to enable/disable the software compliance during pose-based tasks. In practice, we did not experience the tuning drawbacks of admittance framework w.r.t. inconsistencies for steady state commands with different virtual inertia.

C. Haptic Interface

Human operators can interact with robotic systems via different interfaces, such as touch displays, game-pads and haptic devices. A haptic device (or haptic interface) is a robot that works in the master side of a tele-operation system. The human operator can move the master robot by applying force to its handle and reversely the haptic interface can apply forces to the human operator, called "haptic cues" or "haptic feedback." Haptic feedback is used to provide the human operator with information about the state of the remote system.

In our testbed, we used an Omega.6 haptic device, which has six DoFs, the three translations are actuated by independently controlled DC motors and the three rotations are passive. The device communicates through USB 2.0 with the main PC and can be controlled at up to 4 kHz (the faster the control loop the better the force rendering).

III. SYSTEM AUTONOMY FRAMEWORK

Following the discussion on system autonomy in the Technical Report, the Tele-MAGMaS has three distinct operation modes: i) fully autonomous, ii) tele-operated and iii) sharedcontrol. In the following we review these approaches and we show how the human operator can assist or intervene in each approach by means of a haptic device.

A. Fully autonomous MAGMaS

The MAGMaS fully autonomous operation mode is targeted for when the system evolves in a well-structured fully known environment. It relies on two key components: a trajectory planner and finite state machine (FSM); which together compose the task planner, depicted in Fig. 3-top. The FSM is detailed in lower part of Fig. 3 and defines the policy used by the task planner for generating the robot motion trajectories and triggering the grippers actions, based on the robots and environment information and the operator. The autonomous cooperative manipulation phase (S 3 ) can be implemented in a centralized or decentralized way. In centralized cooperative manipulation, both robots are commanded based on the manipulated object position. A detailed description of this approach can be found in [START_REF] Staub | Towards robotic MAGMaS: Multiple aerial-ground manipulator systems[END_REF]. In this case, even if the proposed path planner relies on model inversion, the sensibility to parameter uncertainty and disturbance is mitigated via a disturbance observer, the compliant control of both manipulators and mechanical compliance. This made the full system reasonably robust to parameters uncertainties (see [START_REF] Staub | Towards a flying assistant paradigm: the OTHex[END_REF] for details on the OTHex side). In the decentralized modality, such as leaderfollower approach, during the cooperative manipulation the GM manipulates the object and the flying manipulator assists the GM following its lead and producing an additional upward force on the bar.

B. Tele-operated MAGMaS

At the other hand of the autonomy spectrum, the bilateral tele-operation approach can be used to cope with unknown environment, uncertainties and also to facilitate completion of complex tasks. This is useful in scenarios like USAR, where MAGMaS must work in a partially or even completely unknown environment. In the bilateral tele-operation approach, skilled human operator drive the system in a precise and safe way while provided with force feedback in order to increase their situational awareness on the remote side.

In fact, in the bilateral tele-operation approach, the human intelligence performs the task planning with the help of visual and haptic feedback. If direct visual feedback is not possible, cameras can be mounted on the robots. In the non-cooperative parts of the task (S 1 , S 2 and S 5 ) the human operator drives the robots, while in the cooperative part (S 3 ) the operator command the bar while the robots cooperatively manipulate The operator is provided with a haptic feedback, f hap (see Fig. 4), which depends both on the inertia of the whole system, through f M and on a repulsive viscoelastic virtual force, f o , generated with the purpose of letting the operator feel the obstacles in the environment, this yields

f hap = f M + f o f M = M x x x ẍ x x f o = K P d min + K D ḋmin
where M x x x is the inertia of the system projected at the objects point considered for the tele-operation, x x x and d min and ḋmin denote the shortest distance, components wise, toward an obstacle and its derivative and K P , K D are two positive definite matrices chosen to generate the viscoelastic force. The total rendered force is saturated to satisfy the actuation constraints of the haptic device.

C. Shared-Control MAGMaS

The third possible autonomy level for the system is called shared control, which is the combination of the two previous approaches. An automatic task planner plans the motion commands and changes of logic states and the human operator can modify the planned trajectories. Human operator can locally change the trajectories or change the time law of the trajectories for reacting to changes in the environment. This autonomy level is well suited for semi-structured environment, i.e., partially mapped environment, or fully mapped but with the possible occurrence of unpredictable events, e.g., a construction site or an industrial environment where humans, robots and other machines work together.

IV. SOFTWARE ARCHITECTURE

The Tele-MAGMaS is composed by 3 main robotic parts, i) the LBR-iiwa arm, ii) the OTHex aerial manipulator and iii) the Omega.6 haptic interface, commanded by a centralized controller, which is interfaced with a simulator/visualizer.

In order to develop the experimental framework we relied on the GenoM3 1 abstraction layer, which allows defining middleware-independent software components for robotics (the middleware, e.g., ROS, can then be chosen at compilation time). GenoM3 components are written in C/C++ and can be controlled via tcl-shell, Matlab command line, Matlab-Simulink or middleware specific means. This allows high flexibility in the development and use of the components.

The global software architecture is depicted in Fig. 5. The high-level control of the full system is realized in Matlab-Simulink, linked to the hardware via GenoM3 components (OTHex, Joypad, MoCap, LBR-iiwa) or Matlab S-function drivers (Haptic device, V-REP), which are C++ function natively interacting with Matlab. We chose this approach because we realized that development and test of controller in Matlab-Simulink is considerably faster than in pure C/C++. However, since Matlab-Simulink is not meant for real-time execution, the hardware is commanded via GenoM3 components. These GenoM3 components are essentially drivers for the hardware to satisfy real-time constraints as we chose to keep most of the algorithmic part in Matlab-Simulink. In order to satisfy the hard real-time constraints of the communication with the LBRiiwa, inverse kinematics and other related utilities are performed in the GenoM3 component. In our architecture Matlab-Simulink runs at 500 Hz, the task/path planner (see Sec. III), the human input interpreter, the system state estimator and the OTHex controller all run inside Matlab-Simulink.

This component-based architecture allows easy repartitioning of the load between processes and machines. In our experiment, we chose to use ROS as middleware, which provides sufficient 'real-timeness' for the intended purpose. The component-based design also allows seamless change of the operator inputs, perception components or of the MAG-MaS hardware, as each of those are separated from the main 1 https://git.openrobots.org/projects/genom3/wiki algorithmic part and provide standard interfaces which are not hardware specific.

V. SIMULATION ENVIRONMENT

Simulation has to trade-off accuracy of the simulated behavior and computational load, weighted by the real-time constraint imposed by the human-in-the-loop control scheme. The chosen compromise is a multi rigid-body dynamics integration with collisions processing and reconfigurable joints. The physics simulation is handled by the software V-REP2 , a robot simulator providing a design environment and incorporating the simulation library Bullet. As the overall controller of the system relies on Matlab-Simulink, a bridge enabling bilateral communication between V-REP and the controller has been implemented, exploiting the the Remote API functions, which are integrated within S-functions in Matlab-Simulink, thus providing data transceiver blocks.

Additionally, a graphical user interface was implemented for providing some telemetry to the operator, both in simulations and during real experiments. The displayed information is a visual feedback, via emulated camera mounted on the OTHex, the interaction wrenches and general metadata about the subsystems state. Considering the large amount of possible data to be shown, we have chosen to set up a hierarchical user interface, which allows the user to choose the level of detail. The operator is informed about the status of the system and sub-systems by means of color-changing indicators. Each detailed sub-window can be folded to reduce the visual load on the screen, see Fig. 6.

VI. EXPERIMENTS

The Tele-MAGMaS first experiments were showcased at the Hanover Fair 2017, as part of the finals of the KUKA 2017 Innovation Award, see Fig. 7. The demonstration presented the main features of the Tele-MAGMaS system and a proof-ofconcept application. During the fair week, the demonstration was running every hour (or more), thus demonstrating the high reliability of the proposed system and control architecture. Videos highlighting the key features of the demonstration, i.e., cooperative aerial-ground manipulation with human telepresence, can be found online 3,4 .

Later on, we conducted a set of experiments with a successful co-manipulation of a 2.5 m long bar. The goal is to validate cooperative manipulation with a Tele-MAGMaS for both horizontal displacements of the object and lifting. These basic motions are considered to be representative of the possible construction/decommissioning scenario. Additionally the tele-presence framework is validated for tele-operation and full autonomy. The experiment sequence is depicted in Fig. 8 and consists in the following: at first the OTHex is manually flown to grasp the bar from one of its ends, while the ground robot autonomously grasps the other end. Once both manipulators are attached to the bar the co-manipulation is fully autonomous: they lift the bar from its supports, move it twice along a line in the horizontal plane (blue part) and then synchronously lift the bar up to 30 (green part). Then they bring the bar back to its starting position. This experiment highlights both the vibration stabilization induced by the OTHex and the feasibility of MAGMaS. Indeed we realized that the OTHex was subject to larger external forces than expected, i.e., predicted by the model, due to vibrations of the long object. This is also illustrated in the video contained in the multimedia material. Key quantities of the system are displayed in Figs. 9-10. In particular, Fig. 10 illustrates that a passive joint is sufficient to complete the task, as the bar motion is guided by the GM. Note that the small oscillation would have required some active suppression with an actuated joint, leading to more external wrench on the AM side. This implies that the passive joint as a stabilizing effect on the system. In Fig. 9 the Cartesian wrench at the GM end-effector is depicted (top), which transcribes the wrench exchanged between the manipulated object and the GM. Note that during the initial lifting of the bar from the stand there is a transient with oscillations which disappears during motion, this arises most likely because the model does not take into account the discontinuity of the object breaking contact with the stand, nevertheless, the system responds in a satisfactory manner. Also, note that the torque generated by the object's weight (t y ) diminishes during the lifting of the object, as it is proportional to the cosine of the object inclination. The bottom part of Fig. 9 depicts the joint torque of the GM, it is clear that joint torques start to be noisy when the breaks are disengaged (⇡30 s). Given the particular configuration of the GM most efforts are furnished by the joint A2 (shoulder) and the three joints not aligned with the motion are not solicited (A1, A3, A5). All joint torques are well within their limits, see Tab. I. The AM desired trajectory from the admittance filter is presented in Fig. 10, the good tracking performances validate our approach. Indeed thanks to the multi-directional thrust the error in orientation does not impact the position during the cooperative manipulation.

VII. CONCLUSION AND FUTURE WORK

We introduced a new kind of heterogeneous multi-robot system with a tele-operation capability, the Tele-MAGMaS. This new category of systems leverages the advantages of both mobile ground manipulators and aerial manipulators, in order to perform tasks which would be impossible for homogeneous systems. A scenario of object manipulation (pose tracking), out of the working conditions that a single robot could achieve, is enabled through the tele-operation of the MAGMaS. In order to demonstrate the benefit of aerial manipulators, we designed a fully-actuated platform tailored for bar lifting, the OTHex. Integration and implementation of the framework have been demonstrated in a set of experiments (including a week of live demos during the 2017 Hanover Fair), thus proving the reliability of our approach. Moreover, the software architecture is component-based, hence adding higher-level task planner or other perception modules to operate outdoor could be easily obtained.

Aside from advantages provided by its architecture, the Tele-MAGMaS opens a new field of study, that may eventually bring improvements to the presented method. In particular, although the idea of Tele-MAGMaS is generic, further work is necessary to reach the point where three arbitrary robotic components can be strapped together as a viable Tele-MAGMaS for ad-hoc cooperative manipulation. Currently, the system relies on the Motion Capture in-formation for its localization and general state estimation: our goal is to replace the motion capture system with truly embedded perception in order to run experiments in any environment, especially outdoor. Our future work will also include the extension of the system and framework to a ground manipulator mounted on a mobile platform, which is of paramount importance in USAR or construction missions. It is also important to address the estimation of the load parameters and active rejection of vibrations in the system. Another envisioned plan to improve the work is to robustify the tele-operation and shared control schemes against time delay, to make the system usable in Internet-based tele-operation scenarios in which considerable time delays can occur during internet traffic jams.

Fig. 1 :

 1 Fig. 1: Experimental setup with an LBR-iiwa (3) mounted on a FlewFellow (4). The smartPad (5) is carried by the FlexFellow (4) and the gripper (2) grasps the bar via custom-made claws (1).

Fig. 2 :

 2 Fig.2: Top, OTHex (hovering) with main features highlighted: one of the tilted motors (1), the electronics case (2), the passive joint (3) and the two coordinated grippers (4). Bottom, Comparison of the spatial location of the propellers for a standard hexarotor configuration and the OTHex. The OTHex configuration is such that an aperture of 85 is left to facilitate the manipulation of a long object (green area).

Fig. 3 :

 3 Fig. 3: Top, task planner inputs and outputs for a cooperative manipulation task performed by a Tele-MAGMaS. Bottom, the simple finite state machine (FSM) implemented in our task planner.

Fig. 4 :

 4 Fig.4: Illustration of the bilateral tele-operation scheme, with operator commanding the position p and receiving a force feedback f hap rendering the system inertia and proximity to obstacles.

Fig. 5 :

 5 Fig. 5: Software architecture used in the Tele-MAGMaS project. In green Matlab-Simulink links, in blue C S-function links, in orange GenoM3 links and in black low-level links.

Fig. 6 :

 6 Fig.6: Display of the 3D scene and graphical user interface within VREP, in red the sub-system status indicators (that can be folded), in blue the emulated FVP camera on the OTHex, in green the visualization of the interaction forces and torques.

Fig. 7 :

 7 Fig. 7: A snapshot of the cooperative manipulation state of Tele-MAGMaS using bilateral tele-operation approach performed during the KUKA 2017 Innovation Award at the Hanover Fair.

Fig. 8 :

 8 Fig.8: Time-lapse of a MAGMaS cooperative manipulation task. Both robots are at their initial position (1), approach to the bar (2), grasping the bar (3), cooperative lifting (4), cooperative lateral motion (5), cooperative lifting up to 30 (6-7-8) and release of the bar[START_REF] Kim | Aerial manipulation using a quadrotor with a two DOF robotic arm[END_REF]. See attached multimedia for the corresponding video of the experiment.

Fig. 9 :

 9 Fig.9: On top and middle, external wrench as sensed from the joint sensors and projected on Cartesian space. Bottom, joint torques for each articulation of the LBR-iiwa. The three instants highlighted are LBR-iiwa grasping (1), cooperative lifting (2) and LBR-iiwa un-grasping (3). The blue part highlights the horizontal motion and the green part the bar tilting.

Fig. 10 :

 10 Fig. 10: On top and middle, position and orientation of the OTHex aerial manipulator. Bottom, evolution of the OTHex passive joint angle, during a typical task with free-flight, horizontal motion and object tilting.The four instants highlighted are take-off (a), OTHex grasping (b), cooperative lifting (c) and un-grasping(d). The blue part highlights the horizontal motion and the green part the bar tilting.

CNRS, LAAS-CNRS, Université de Toulouse, Toulouse, France,

http://www.coppeliarobotics.com/

https://vimeo.com/217252361

https://youtu.be/GRnGSvJGUKk

VIII. ACKNOWLEDGMENTS Partially funded by: EU Horizon 2020 grant agreement No 644271 AEROARMS; ANR, Project ANR-17-CE33-0007 MuRoPhen; EU FP7/2007-2013 grant agreement n. 601165 WEARHAP; Basic Science Research Program (2015R1A2A1A15055616) of the NRF funded by MSIP, Korea. We thank KUKA and its Innovation Department, Schunk France and Mathworks for their support in this research project. Other collaborators, who helped during KUKA 2017 Innovation Award, are acknowledged in the Tech Report for their contribution