Dima Abi-Abdallah
email: dima.abiabdallah@utc.fr

Agnès Drochon

Vincent Robin

Odette Fokapu

Cardiac and respiratory MRI gating using combined wavelet sub-band decomposition and adaptive filtering

Keywords: ECG, RLS, filter banks, real time synchronization

published or not. The documents may come

Cardiac and respiratory MRI gating using combined wavelet sub-band decomposition and adaptive filtering

Introduction

The constantly evolving MR techniques have resulted in an increasing number of cardiac MRI scans in order to gain important insights into the functional and metabolic bases of heart disease.

However, MRI observations of a moving organ, such as the heart, require synchronization so that the image acquisition might be accurately combined with its phase motion, thus eliminating the movement related blurring effects. It can be assumed that the heart moves independently with two motions: a contractile motion associated with the cardiac cycle, changing primarily the shape of the heart, and a respiratory motion that changes its position and orientation as the diaphragm and chest wall move with each breath [START_REF] Klein | Real-time system for respiratory-cardiac gating in positron tomography[END_REF]. Therefore, efficient synchronization techniques must associate cardiac gating with respiratory blanking, where the heart beat governed acquisitions are only enabled during at-rest periods (between breaths) to overcome the respiratory motion artifacts [START_REF] Cassidy | Assessment of Motion Gated Strategies for Mouse Magnetic Resonance at High Magnetic Fields[END_REF], [START_REF] Maı | Effects of Breathing and Cardiac Motion on Spatial Resolution in the Microscopic Imaging of Rodents[END_REF], [START_REF] Wiesmann | High-resolution MRI with cardiac and respiratory gating allows for accurate in vivo atherosclerotic plaque visualization in the murine aortic arch[END_REF].

Cardiac gating is usually done by detecting the R peaks on the simultaneously recorded electrocardiogram (ECG), which are then used to trigger the consecutive image acquisition sequences. However, accomplishing an accurate R peak detection, hence a correct synchronization is often obstructed by the high corruption levels of the ECG signal due to electromagnetic effects. Especially when high resolution imaging is performed where high static magnetic fields as well as strong and fast switching gradients are needed. In fact the signal S(t), recorded in the NMR environment, is contaminated by many artifact signals :

 A signal due to patient related sensor motions in the magnetic field such as respiration, heart beat and voluntary motion. This artifact becomes more pronounced as the static magnetic field strength increases.

 A high frequency artifact resulting from the radio frequency pulses. This can be easily removed with low pass filtering.

 Noise generated by the temporal variations of the magnetic field gradients. These are spurious signals that often resemble the QRS spikes and cause erroneous detection. In the case of small animal imaging, these artifacts are further aggravated by the exceptionally high heart rate (400-600 beats per minute), the high respiratory rate (30-60 breaths per minute) [START_REF] Cassidy | Assessment of Motion Gated Strategies for Mouse Magnetic Resonance at High Magnetic Fields[END_REF], and the low amplitude of the ECG signal (just a few millivolts).

 An artifact signal induced by the magnetohydrodynamic (MHD) effect. It results mainly in an elevation of the T wave that might reach amplitudes as high as the R wave, consequently leading to misgating.

Various approaches aiming to reduce the ECG signal artifacts have been previously considered in order to produce a synchronization signal. However, with the evolution of high resolution imaging techniques, the contamination levels have greatly increased, thus degrading the efficiency of classical analogical filtering methods [START_REF] Damji | Interference suppression in cardiac synchronization system operating in a high magnetic field NMR imaging system[END_REF], [START_REF] Felblinger | Electrocardiogram Acquisition during MR Examinations for Patient Monitoring and Sequence Triggering[END_REF], [START_REF] Laudon | Minimizing Interference from Magnetic Resonance Imagers During Electrocardiography[END_REF]. Numerical methods were later proposed [START_REF] Abi Abdallah | Reference signal extraction from corrupted ECG using wavelet decomposition for MRI sequence triggering: application to small animals[END_REF], [START_REF] Felblinger | Restoration of Electrophysiological Signals Distorted by Inductive Effects of Magnetic Field Gradients during MR Sequences[END_REF], but could not be directly implemented in real time, yet it's indisputable that real time filtering is mandatory for achieving an efficient cardiac synchronization. Moreover the algorithm must account for the subject's physiological variations during the examination and for the artifact waveform variations, which differ according to the applied imaging sequence, as well as patient movements. Very few real time algorithms have been proposed [START_REF] Fischer | Novel real-time R-wave detection algorithm based on the vectorcardiogram for accurate gated magnetic resonance acquisitions[END_REF] and [START_REF] Abächerli | Suppression of MR gradient artefacts on electrophysiological signals based on an adaptive real-time filter with LMS coefficient updates[END_REF], and the solutions impose various technical constraints. The first algorithm exploits a vectocardiogram (VCG) signal and necessitates the use of multiple sensors. However, positioning a large number of sensors in a NMR environment is not always simple, and might have consequences on the image quality. Not to mention the near impossibility of placing so many sensors on a small animal's chest. The second method, based on adaptive filtering, uses the machine's gradient signals as references. For this, one must have access to the MRI machine's gradient amplifier. In this paper, we propose a real time method, for achieving double cardiac and respiratory gating, that is quite simple and inexpensive in terms of materials.

In order to combine cardiac gating with respiratory blanking, the respiration motion must be detected. This is usually done using motion detectors; however this is not always comfortable for some patients and is cumbersome for very small animals. Nevertheless, the respiration motion information is enclosed within the electrocardiogram as a modulating signal [START_REF] Felblinger | Amplitude Demodulation of the Electrocardiogram Signal for Respiration Monitoring and Compensation during MR Examinations[END_REF]: In fact, the respiratory motion of the thorax, heart and diaphragm changes the orientation of the heart's electric axis, resulting in an amplitude modulation of the ECG, which is itself a vectorial projection of the electrical activity of the heart.

In a previous work [START_REF] Abi Abdallah | Reference signal extraction from corrupted ECG using wavelet decomposition for MRI sequence triggering: application to small animals[END_REF] we developed an off-line algorithm for cardiac triggering based on wavelet sub-band decomposition. By separating the sub-bands where the QRS energy is maximal, a reference signal was extracted from the ECG, and was then subjected to thresholding to produce a trigger signal. The algorithm has proven its efficiency even for very low signal to noise ratio (SNR) signals. In the present work we exploit this algorithm to produce an adapted filter, for real time use, that yields maximum QRS energy using the most appropriate wavelet. Furthermore, by applying an adequate electrode positioning, that maximizes the respiratory dependent modulation, we make use of the respiration information contained in the ECG signal in order to generate, in real time, a synchronization signal that takes into account the respiratory movements as well.

Method

MRI technique and signal acquisition

The MRI experimental setup (figure 1) incorporates a 2 Tesla OXFORD 85/310 horizontal cryomagnet with a 180 mT.m -1 gradient system using a home-made Alderman-Grant coil adapted to the small size of mice. A SMIS MR-operator control console, which has a scan management system with TTL signal input used as a timer for gated sequences.

-6 -For cardiac and respiratory signal amplification and transmission outside the Faraday cage, an optical non-magnetic module was used. Signal detection was performed using 3M Red Dot neonatal electrodes, (3M, St Paul, MN). The ECG active electrodes were placed parallel to the heart axis (45°) in order to obtain maximum modulation. The respiration modulated ECG signals were thus recorded during MR scanning according to three high resolution imaging sequences -Gradient Echo (GE), Fast Spin Echo (FSE) and Inversion Recovery Spin Echo (IRSE)-inducing different types of contamination signals. For all three sequences, coronal and transversal planes were imaged. 4 to 6 slices of 1 mm thickness were acquired with an in-plane resolution of 100 µm, a field of view (FOV) of 27mm and a flip angle ( ) of 60°. MR images were reconstructed using a matrix of 256 x 256 pixels.

Repetition time (TR) to echo time (TE) ratios were 385/10 ms for GE, and 2240/50 ms for T2-weighted FSE and IRSE. The system for digital processing of the contaminated modulated ECG, used for sampling, filtering and gating signal generation, consists of a National Instruments' DAQCard 6024E and a real time structure developed using Matlab ® & Simulink ® (The MathWorks Inc., Natick, MA).

Cardiac-respiratory trigger extraction process

The cardiac-respiratory trigger extraction process, illustrated in figure 2, is made up of two simultaneous trigger production operations: the R peaks of the ECG are detected to yield a cardiac trigger, while the respiratory modulation signal is extracted, from that same ECG, to produce the respiration trigger. The final synchronization signal is none other than the combination of both triggers (logical AND) so that the image acquisitions would only be set off by the R waves during the low respiratory phase. that, P and T waves only present a significant spectral density up to 10 Hz, and that most of the QRS power lies in the 3-20 Hz band, beyond which the complex energy decreases gradually [START_REF] Thakor | Estimation of QRS Complex Power Spectra for Design of a QRS Filter[END_REF]. Thus, a strategical approach for detecting heartbeats would be to analyze different sub-bands of the ECG using wavelet based filter banks [START_REF] Afonso | ECG Beat Detection Using Filter Banks[END_REF]. Figure 3 shows an illustration of a previously developed algorithm [START_REF] Abi Abdallah | Reference signal extraction from corrupted ECG using wavelet decomposition for MRI sequence triggering: application to small animals[END_REF]. It consists in decomposing the ECG signal into multiple scales (8 for a signal sampled at f s =1Khz), then reconstructing the details and keeping only those that contain the maximum of the QRS energy in order to compose a reference signal. This signal is then subjected to a double threshold comparator to produce a cardiac trigger. The algorithm was tested on simulated as well as small rodents ECGs that had been contaminated by one of three imaging sequences (Gradient Echo, Fast Spin Echo and Inversion Recovery Spin Echo) and some analyzing wavelets revealed having outstanding performances according to the signal and the sequence type.

2.2.1.Theoretical Basis

But even though the wavelet decomposition algorithm gave proof of an almost infallible QRS detection, it could not be applied in real time for image synchronization, because the decomposition into 8 levels followed by the reconstruction of the 6 th and 7 th detail induce a significant delay. Equation (1) gives the time necessary for the reconstruction of the n th detail using a wavelet with M coefficients

   1 1 2 1   M f n (1) -8 -
For instance, if we use the Db6 wavelet (M=6), to reconstruct the 7th detail of a signal (n=7) sampled at f s =1000Hz we get a 0,635 seconds delay, which is almost equal to one heart cycle.

In order to achieve real time synchronization, while making use of the wavelets capabilities, we conceived an off-line pre-calculation phase, that precedes the actual real time trigger generation. It consists of a coefficients computation phase, where an optimal QRS filter is calculated based on wavelet decomposition, which can afterwards be easily used for real time processing.

The purpose of this first stage is to construct a filter that would yield the same reference signal as the wavelet decomposition method while assuring a minimum delay. A good approach for this filter optimization problem is to minimize the mean square value of the error signal that is defined as the difference between the desired response and the actual filter output. For stationary inputs, the resulting solution is commonly known as the Wiener filter, which is said to be optimum in the mean square sense. However the Wiener filter is inadequate for dealing with situations in which the signal and/or the noise are non stationary (such as the ECG). In such situations the filter has to assume a time varying form. An efficient method for realizing this is to use an adaptive filter.

Adaptive filter

By such a device we mean one that is self designing, in that the adaptive filter relies for its operation on a recursive algorithm which makes it possible for the filter to perform satisfactorily in an environment where complete knowledge of the relevant signal characteristics is not available. In a non stationary environment the algorithm offers a tracking capability whereby it can track time variations in the statistics of the input data provided that the variations are sufficiently slow. As a discrete consequence of the application of a recursive algorithm, where the parameters of an adaptive filter are

Respiratory phases delineation

It is well-known that respiration influences electrocardiographic measurements. During the respiratory cycle, chest movements and alterations in thoracic impedance, as the lungs fill and empty, alter the ECG morphology. In fact, the ECG is a projection of the heart's electrical activity vector along the derivation direction (as perceived by the leads). The respiratory motion of the thorax, heart and diaphragm changes this vector's orientation, resulting in an amplitude modulation of the ECG in the direction of the cardiac mean electrical axis. It has been shown [START_REF] Felblinger | Amplitude Demodulation of the Electrocardiogram Signal for Respiration Monitoring and Compensation during MR Examinations[END_REF] that positioning the ECG electrodes along the heart axis (45°, V2 V3) can maximize the respiratory dependent modulation, thus facilitating the extraction of the respiratory information without however affecting the cardiac signal quality. Note that when electrodes are placed at this position, inhaling would produce a decrease of amplitude (figure 4).

An easy way to retrieve the respiratory modulating signal is to apply a highly selective low pass filter L(n) to the recorded ECG. Using Matlab ® 's Digital Filter Design Tool we calculated a FIR equiripple filter. Note that the equiripple design produces the most -11 -efficient filters, that is, that they just meet the specification with the least number of coefficients. A cut-off frequency of 1.2Hz was deemed acceptable since the rodents respiration frequency is usually less than 1Hz, especially when anesthetized, as is the case during NMR imaging. Naturally, it comes to one's attention that due to the narrow bandwidth of the low pass filter, the extracted respiratory signal would be greatly delayed. In fact, the designed filter L(n) has over 3000 coefficients (for a sampling frequency of 0.5Khz) leading to an approximate 1.5 seconds delay, thus compromising the respiratory synchronization. Nonetheless, given that the respiration is somewhat periodic, the delay can be simply compensated by a prediction of the upcoming signal along a period which value is continuously updated.

2.2.2.Algorithm

i.

Off-line filter computation phase

The off-line phase algorithm proceeds as follows:

 Record a small portion of signal S ori (t) (minimum duration of two or three respiratory periods ) while the imaging sequence gradients are on.

 Process it to produce a reference signal S ref (t) based on the previous method [START_REF] Abi Abdallah | Reference signal extraction from corrupted ECG using wavelet decomposition for MRI sequence triggering: application to small animals[END_REF] using the most appropriate wavelet.

 Calculate the optimal filter coefficients using the RLS algorithm by setting S ori (t) as input and S ref (t) as the desired signal. A number of 64 coefficients for the filter seemed to provide a good compromise for getting low error while having short delay and computation time. (the larger the coefficients number, the smallest the error but with longer delay and computation time)

Thus we obtain a filter F(n) that is adapted to both the ECG signal and the applied imaging sequence.

ii. Real time trigger production

This second phase consists of two simultaneous processes: an R peak detection producing a cardiac trigger, along with a respiratory trigger extraction. The final cardiac-respiratory signal used for triggering the imaging sequences is composed by combining both triggers.

Cardiac trigger

The recorded NMR contaminated signal is filtered by the calculated minimum delay F(n) to obtain an artifact free ECG, where the QRS peaks constitute the most prominent features. This signal can then undergo a simple threshold comparison to produce the cardiac trigger. Note that the threshold values are adjustable via a graphical interface.

Respiratory trigger

The collected ECG is filtered using the selective low pass filter L(n), to produce the respiration related modulating signal. Then, the respiratory gating trigger is produced by subjecting the filtered signal to a simple threshold at the mean value.

The algorithm described above was implemented on Simulink ® 's Real Time Windows Target platform. The conceived Simulink ® model (outlined in figure 5) was compiled using Matlab ® 's

Real Time Workshop, to create a real time software application. The application was then attached to a user friendly interface (figure 6) that allows the experimenter to launch the different procedures, and to easily pick and set the algorithm's parameters: At the beginning of the session, after recording a chunk of signal, he can pick a wavelet, test its performances and adjust the R peaks detection thresholds. Once the wavelet is selected, he can run the off-line calculation phase, and then start the online trigger production.

Results

The cardiac-respiratory trigger extraction algorithm was tested for mice cardiac MRI synchronization, using three imaging sequences GE, FSE and IRSE. The animal study was conducted with approval from the Institutional Animal Care and Use Committee of the University of Lyon 1. The mice were anesthetized by intraperitoneal injections of pentobarbital sodium (50 mg/kg, Sanofi, France).

According to the applied sequence the algorithm was run using the most appropriate wavelet [START_REF] Abi Abdallah | Reference signal extraction from corrupted ECG using wavelet decomposition for MRI sequence triggering: application to small animals[END_REF]:

Coif3 for GE, Sym8 for FSE and Sym4 for IRSE. The real time process was able to generate effective online cardiac-respiratory triggering for all 3 sequences. An illustration is given in figure 7. The QRS peak detections, achieved in real time by the adapted filter, were later compared to the triggers that the wavelets would have produced off-line using the prior algorithm. For each sequence, the detection error rates were as low as those given by the best performing wavelets [START_REF] Abi Abdallah | Reference signal extraction from corrupted ECG using wavelet decomposition for MRI sequence triggering: application to small animals[END_REF]. Furthermore, we did not perceive any significant influence of slice orientation on the gating process. The error rates were similar whether the signals were acquired during coronal or transversal slice acquisition. In all cases we attained an almost perfect QRS detection.

The evaluation of the respiration detection algorithm was carried out on simulated signals. In fact, it was difficult to delineate the respiration phases manually because the sequence noise alters the R peak amplitude thus preventing correct visual detection of the signal modulation.

Previously recorded contaminated ECGs [START_REF] Abi Abdallah | Reference signal extraction from corrupted ECG using wavelet decomposition for MRI sequence triggering: application to small animals[END_REF] were modulated by a sine wave that comprises variations of amplitude (up to 150%) and frequency (0.2Hz -1Hz). The algorithm performed just as expected: the modulation was extracted successfully, however when a change in frequency and/or amplitude occurs the algorithm takes around one respiration period or two to adapt to the changes. During that time, a high respiration phase might be mistaken for a low one and taken into account for synchronization. Once the algorithm recovers, it will once again correctly retrieve the modulating signal (figure 8). Note that the modulation retrieval was independent of the sequence type; we observed the same behavior with all three sequences (the noise is eliminated by the low pass filter).

The global evaluation of the algorithm was done in terms of image quality. The images, acquired with synchronization, revealed a very clear enhancement compared to cardiac MRIs taken without cardiac or respiratory triggering (figures 9, 10 and 11). Without any synchronization the image is quite blurry especially in the parts that move along with each heart beat. With cardiac triggering, clearer structures appear on the image close to the heart. When applying both cardiac and respiratory synchronization, we observe a better contrast in the whole image, and the structures are more emphasized.

Discussion

Cardiac-respiratory synchronization is essential for obtaining good quality cardiac MR images.

However, the perturbation of the ECG signals acquired in the NMR environment, especially during small animal imaging, greatly obstructs gating tasks.

In this work, a method for triggering MR imaging sequences on cardiac and respiration rate is

presented. An adapted filtering, based on wavelet sub-band decomposition, is applied to the contaminated simultaneously recorded modulated ECG to rid it of the NMR environment artifacts, while emphasizing the QRS peaks. This allows the extraction of a cardiac trigger that is completely synchronized with the heart beats. The presence of an off-line pre-calculation phase is very convenient, in the sense that we can add any kind of operation to the wavelet process to improve the reference signal, and then regroup all of it in a single filter with minimum delay.

In addition, based on the fact that the amplitude of the ECG varies according to the respiration motions, a demodulation of this "cardiac respiratory" signal yields a respiration curve which is used to build the respiratory triggering signal. A final trigger is then obtained by combining both cardiac and respiratory triggers and used for MRI gating. The algorithm has proven efficiency for small rodents' cardiac imaging. The trigger extraction was almost flawless and the image quality was greatly enhanced.

Our method could constitute a good alternative for the solutions proposed by Fischer et. al. [START_REF] Fischer | Novel real-time R-wave detection algorithm based on the vectorcardiogram for accurate gated magnetic resonance acquisitions[END_REF] and Abächerli et. al. [START_REF] Abächerli | Suppression of MR gradient artefacts on electrophysiological signals based on an adaptive real-time filter with LMS coefficient updates[END_REF]. Unlike these solutions, the technique presented in this paper does not have any technical constraints. The use of one ECG sensor is quite practical and can be easily positioned even for small animals. Furthermore, the extracted reference signal is independent of the MRI machine.

As for the algorithm, the following remarks can be discussed:

Adaptive improvement: For images which acquisitions take a relatively long time to complete, significant variations may occur in the signal morphology. Therefore, better results would be obtained if the filter coefficients were recalculated during the synchronization operation and updated every few minutes, in order to adapt the processing to the signal changes. Nevertheless, since no major signal variations occurred in all ECGs recorded throughout our experiments, using the same filter during the whole imaging session was just as good as using updated filters.

That is why calculating the filter only in an off-line pre-processing phase was satisfactory enough. However, in a future work we will incorporate continuous filter estimations and coefficients updating during the real time process. We will also give the user the possibility to adjust the R peak detection thresholds during the online phase.

Off-line phase duration:

The length of the signal to be used for the pre-calculation phase can be defined by the experimenter as he sees adequate. However a minimum duration of two or three respiratory periods is required so that the algorithm could function properly. Even though the RLS convergence is a lot faster than that, additional time is needed to make an initial estimation of some parameters that are necessary to get the respiration trigger production started.

Respiration extraction: the implemented respiration extraction algorithm could be evaluated more effectively by using a traditional respiration detector, and comparing the extracted signal with that of the detector. However placing any additional detector on the mice was nearly impossible given the small thorax area. We could also use a ventilator system and impose breathing frequency changes; nevertheless, given the narrow dimensions of the MRI tunnel and of the antenna, which have been maximally reduced to enhance the image quality, placing such system would be very difficult. It is conceivable, though, to test it on human subjects in a later stage. For the time being we have only evaluated the algorithm using simulated signals. The modulation variations imposed on these signals imitated, as close as possible, the variations that would most likely occur in real signals: smooth frequency and amplitude modifications. In that case, errors would only occur during the respiration period that follows the change. Even so, this did not seem to impair the image enhancement. Note that in order for the method to function properly, there has to be steady stages where the signal is not altered. For example, if the signal has a frequency that is constantly changing, say in a steep slope linear manner, the detection error rate suffers a great deal. In fact, the algorithm needs to predict the upcoming period, so if the latter keeps changing the prediction would always be faulty.

The effect of cardiac pathologies on the performances of the extraction method still needs to be examined more thoroughly. It has come to our attention that some diseases, might introduce very low frequency components into the ECG signal. These components could get through the low pass filter, and distort the extracted modulating signal. The filtered signal would not only contain the respiration but it will also include some additional disruptions.

Conclusion

The method proposed in this paper for synchronizing image sequences on the subject's physiological rhythm is simple and original. Using one ECG sensor to achieve double cardiac -17 -and respiratory synchronization is a potentially useful technique, especially in the NMR environment where the number of measurement devices must be as limited as possible. The efficiency of the method relies mainly on the idea of combining the wavelet transform with adaptive filtering, thus solving, with this association, most frequently encountered difficulties in such applications. The developed user interface allows the experimenter to optimize the denoising parameters according to the MRI sequence type, and enhance the trigger signal quality throughout the imaging session. The algorithm has been successfully applied to modulated ECG signals recorded on healthy small animals. It now needs to be tested on signals recorded in pathological cases. The signal is subjected to an adapted filtering to enable QRS detection, as well as a low pass filtering to extract the respiration signal. The cardiac and respiratory triggers are then combined to produce a synchronization signal that controls the image sequence triggering The final trigger reproduces the occurrences of the R peaks only during the exhalation phase i.e.

when the respiratory modulation signal is higher than its mean value. Around the 28 th second, the wave's frequency is doubled and its amplitude is reduced by 40%.

During the adaptation period the algorithm misses 3 peaks and detects 2 false ones. Otherwise all detected peaks are correctly detected within the low respiration phase. Note that the represented modulating signal is shifted for visibility's sake. -23 -

  characterized by a cyclic occurrence of patterns (QRS complexes, P and T waves) with different frequency contents. Power spectral analysis of the ECG show

  6. increment n by 1 and go to step 2 where n : current iteration number u(n): input at iteration n P(n):inverse of the correlation matrix k(n): gain vector ) (n h  : filter coefficient vector η(n) : estimation error d(n) : desired output i.

Figure 1 -

 1 Figure 1 -MRI gating experimental setup

Figure 2 -

 2 Figure 2 -Real time cardiac-respiratory trigger extraction:

Figure 3 -

 3 Figure 3 -Overview of the off-line cardiac trigger extraction algorithm

Figure 4 -

 4 Figure 4 -Respiratory modulation of the ECG signal

Figure 5 -

 5 Figure 5 -Simulink® bloc model

Figure 6 -

 6 Figure 6 -The synchronization software's graphical user interface

Figure 7 -

 7 Figure 7 -Synchronization signal extraction:

Figure 8 -

 8 Figure 8 -Trigger extraction from a simulated ECG:

Figure 9 -

 9 Figure 9 -Mouse cardiac MRIs acquired with a GE sequence:

Figure 10 -

 10 Figure 10 -Mouse cardiac MRIs acquired with a FSE sequence:

Figure 11 -

 11 Figure 11 -Mouse cardiac MRIs acquired with an IRSE sequence:

Acknowledgments

The study was supported by the Ministry of Industry of France (RNTS 2003).

The authors would like to thank the "Laboratoire de RMN, UMR CNRS 5012",University Claude Bernard Lyon 1 for the use of their MR machine and the support by their MR staff.

Another thank you is addressed to the J. Felblinger of the IADI laboratory University Hospital Nancy,Brabois at Vandoeuvre-les-Nancy for his technical help.

updated from one iteration to the next, the parameters become data dependent. The operation of an adaptive filtering algorithm involves two basic processes: a filtering process designed to produce an output in response to a sequence of input data, and an adaptive process, the purpose of which is to provide a mechanism for the adaptive control of an adjustable set of parameters used in the filtering process. These 2 processes work interactively with each other.

There is no unique approach to the development of an adaptive filter; there is a variety of recursive algorithms each of which offers desirable features of its own. The most popular solutions are the LMS (least mean square) and the RLS (recursive least square). In this study we preferred the RLS approach because it is more robust than the LMS, and it has a faster convergence rate: Contrarily to the LMS, where the correction, applied in updating the old estimate of the coefficient vector, is based on the instantaneous sample value of the tap-input vector and the error signal, the RLS utilizes all the information contained in the input data from the start of the adaptation up to the present [START_REF] Haykin | Standard recursive least-squares estimation in adaptive filter theory[END_REF]. Since it converges fast, the RLS filter computation allows us to record a smaller chunk of signal to be used in this stage. The superior performance of the RLS algorithm compared to the LMS, however, is attained at the expense of a large increase in computational complexity. But this is of no inconvenience to us since we operate in an off-line computation process, and we only have a small number of coefficients (short filter delay) to calculate.

Summary of the RLS algorithm:

1. let n=1 2. compute the gain factor : 3. compute the true estimation error : 4. update the estimate of the coefficient vector: 5. update the error correlation matrix:


M Bui

V Frémont

D Boukerroui

P Letort

Multi-sensors people detection system for heavy machines

Keywords: Heavy machines, sensor fusion, pedestrian detection, deformable part model, fisheye, histogram of oriented gradients

In this paper, we propose a multi-sensors system to detect people in the context of construction sites using heavy machines. The system includes a LIght Detection And Ranging (Lidar) sensor and a fisheye camera. We present an effective method to determine regions of interest (ROI) on fisheye images using Lidar data, which can be used with different sensors configurations. A Deformable Part Model (DPM) approach is adapted and used as the main people detector on image. We also present a specific dataset built using the multi-sensor system mounted on a heavy machine for evaluation.

I. INTRODUCTION

Heavy machines in construction sites constitute a high risk of serious injuries for people working in their surroundings. Due to their large different functionalities and shapes, accidents caused by heavy machines are various. An investigations of accidents in the period starting from 1997 to 2008 for France was conducted and published in the EPICEA database 1 [START_REF] Marsot | Collisions engins-piétons, analyse des récits daccidents de la base epicea[END_REF]. From 2157 accident cases of all kinds, there are about 15% caused by collisions between machines and people. In view of these worrying statistics, efforts have been made to improve the safety for people working around heavy machines. To better prevent than cure, research has moved toward active safety intelligent systems which are able to predict dangerous situations and anticipate the accidents. They are referred as advanced driver assistance systems (ADAS), in the sense that they help the driver by providing warnings, assistance to take decisions and even taking automatic evasive actions in extreme cases. Notwithstanding many years of progress the problem of detecting obstacles and recognizing people from other ordinary objects is still a tough issue.

A. Related works

To solve the problem of people detection, different sensors have been used and combined. The most popular type of sensor for people detection is the camera. It offers the best trade-off as a low-cost and polyvalent sensor. Recently, following the Euro New Car Assessment Program (NCAP) recommendations, people detection systems in automobile using cameras known important progresses [START_REF] Dollár | The fastest pedestrian detector in the west[END_REF], [START_REF] Enzweiler | A multi-level mixture-of-experts framework for pedestrian classification[END_REF], [START_REF] Felzenszwalb | Object detection with discriminatively trained part based models[END_REF]. Although the objective is similar in both automobiles and heavy machines fields, we can clearly distinguish between the two. Indeed, in the automobile context, cars need to stop if there is an obstacle, no matter if it is a pedestrian or not. The task of recognizing people is more important for heavy machines where the main requirement is human's safety and since heavy machine can pass through specific obstacles (rocks, holes, woods, etc.). Besides, cars often operate at a higher speed and on straight ways. While it is important for the system on automobile to be able to detect people at far distances, heavy machines need preferably a larger field of view (FOV) to cover the nearby area that constitutes the main dangerous zone. Moreover, construction machines often have a complicated shape and a larger size, which can also benefit from the large FOV. These features reveal the necessity to use wide angle sensors such as fisheye cameras as the main sensor for people detection and recognition [START_REF] Bui | People detection in heavy machines applications[END_REF].

However, objects detection in fisheye images is a challenging task. Indeed, unlike in perspective images, the objects appearances are strongly distorted. Wrapping the image into a local perspective image is the direct way to avoid nonperspective deformations but, besides adding computational load, this approach also creates undesirable effects [START_REF] Bülow | Spherical diffusion for 3d surface smoothing[END_REF], [START_REF] Daniilidis | Image processing in catadioptric planes: Spatiotemporal derivatives and optical flow computation[END_REF]. Recently, the Deformable Parts Model (DPM) [START_REF] Felzenszwalb | Object detection with discriminatively trained part based models[END_REF] and its variants have gained a lot of attentions in object detection and recognition. The DPM represents an object model as different parts floating around their reference locations and finds the associated optimal part-configuration at every root position. This elegant way of representing complex shapes of objects, brings a lot of benefit in detecting people in distorted appearances. Still, the DPM approach has few drawbacks that can be stated here:

• The high performance of the DPM approach comes with a high computation cost and complexity. • The DPM approach using the sliding windows detection method needs to scan all over the image to detect the person. Moreover, it is not invariant to the rotation of the person's appearance on image. This drawback is critical in fisheye images. For example, it might not be suitable to detect people in images from a camera that looks down from high positions (see Fig. 6).

Since cameras are low-cost sensors and since they exhibit high computation load, they can be combined with other sensors. For example, range sensors, like radar, Light Detection And Ranging (Lidar) and ultrasonic [START_REF] Arras | Using boosted features for the detection of people in 2d range data[END_REF], which have good performance in detecting obstacles. Heavy machines often work in complicated terrains with a lot of nearby objects.

Sometimes they even need to crush these obstacles. In these situations, range sensors will trigger a permanent alarm, which is useless and annoying for the drivers. Range sensors are still very valuable in our context as they provide very precise and robust measurements. Therefore, such sensors are very good options in combination with cameras [START_REF] Pradeep K Atrey | Multimodal fusion for multimedia analysis: a survey[END_REF], [START_REF] García | Context aided fusion procedure for road safety application[END_REF], [START_REF] Oliveira | Context-aware pedestrian detection using lidar[END_REF], [START_REF] Premebida | Fusing lidar, camera and semantic information: A context-based approach for pedestrian detection[END_REF].

B. Article's contributions and structure

The contributions of this paper are two folds. The first one lies in the combination of Lidar data to bypass the drawbacks of the DPM approach in people detection. The specific difficulties to process fisheye images in the context of heavy machines are taken into consideration. The proposed fusion system has various advantages:

• It helps to reduce the detection area on images resulting in a faster detection. • It helps to give a lower false detection rate without sacrificing the detection performance. • It helps to adjust the rotation of objects' appearances caused by image distortions in different camera view angles. The second contribution concerns the evaluation of our approach on real data acquired in the context of heavy machine.

The paper is organized as follows. First, two adaptive sensors configurations on heavy machines and details about our dataset are given in section II. Then, details on the data fusion algorithm and the projection process of the regions of interest (ROIs) on fisheye images are explained in section III. The evaluation protocol and the experimental results are discussed in section IV. Finally, section V presents the conclusions and the perspectives of our work.

II. ACQUISITION SYSTEM AND DATASET

Datasets take a very important role in the process of the development of an obstacle detection algorithm. Indeed, a well-defined dataset is not only useful for evaluating the approach but it also takes part in the training process to improve the performance of the detector. To the best of our knowledge, there are no others available datasets which provides at the same time synchronized fisheye images and Lidar data in the context of heavy machines environments. Moreover, this context has some special features, such as the outdoor changing light conditions, the strong vibrations coming from the engine, and the brutal shocks that might make the detection process much harder. There are essential needs for a new dataset in order to identify conditions under which current detectors fail. We will focus our attention on these difficult cases.

A. System Configuration

The danger is highly dependent on the action taken by the machine and its direction. Based on the survey of accidents caused by heavy machines (Tab. I), we can conclude that accidents rarely happen on the sides of the machines, except for "rotating-base machines", such as excavators. The back and the front of a machine in motion are the most dangerous parts. Fisheye cameras mounted at the back and in front of the machine seem to be a good option in these cases.

We defined two positions for the fisheye camera denoted by index 1 and 2 in Fig. 1 which led to two different sensors configurations of the system:

• Configuration 1: The two sensors are at position 1. In this configuration, the sensors are kept at low position (height h = 110cm) and are parallel to the ground plane. This is the commonly used configuration of many people detection systems because the appearance of the person on the image is relatively the same at any relative position in the FOV of the camera. This advantage is not conserved in fisheye image since the object appearance distortions depends on its distance and its angle to the camera. A quantitative analysis of this distortion phenomenon is presented in [START_REF] Bui | People detection in heavy machines applications[END_REF]. • Configuration 2: The fisheye camera is mounted at position 2 (height h = 210cm, at the forks level), looking down with an angle of 30°. The Lidar is at the same position as in configuration 1. The high position of the camera is commonly used in the context of heavy machines. The reasons are to get a better coverage of blind angles around the machine and to avoid collisions that can damage the sensors. In our case, it gives an additional advantage in observing the whole people appearance at a very close range, even when the person touches the machine. The main drawback is the complicated transformation of the person's appearance obtained on fisheye images.

Let assume that the Lidar is always parallel to the ground.

In the targeted application, this constraint will be broken whenever the heavy machine is moving on a non-flat ground, or simply due to large vehicle vibrations with respect to the ground. In practice, these artifacts are minor when we use fisheye camera. A way to avoid this constraint is to use a multi-layers Lidar, along with a method to track the vehicle pitch.

B. The Acquisition System

The proposed dataset consists of images captured from one fisheye camera (Point Grey Firefly MV USB2.0), one conventional camera (Sony PlayStation Eye for PS3) and one range-sensor (LIDAR Hokuyo UTM-30LX-EW). The sampling frequency is 10Hz for both cameras and 40Hz for the Lidar. The middleware Pacpus2 is used to manage the data acquisition process. Pacpus can handle data coming from different sensors in real time and save them to the hard disk driver (HDD) for replay and post-processing. These data are timestamped and synchronized so they can be used later to simulate real-time processing (Fig. 2).

These data are taken on board of a telescopic-forklift, namely a Bobcat-TL470, as shown in Fig. 3 (a professional driver has been recruited to operate the machine during the experiments). The experiments are divided into 4 scenarios for each configuration. All of the scenarios are predefined and took place under security control. In most of the data sequences, the machine moves on a clear terrain. The scenarios defined within the experiments aim to simulate frequently meet situations on a construction site. People wore Table II: The testing dataset statistics.

different kind of clothes, including helmet, reflective vests and civil clothes. Different situations of occlusions were also simulated.

The fisheye images are partially annotated and used for testing purpose. The annotations for the ground truth of these image sequences are done using the labeling tool of Dollár et al [START_REF] Dollár | Pedestrian detection: An evaluation of the state of the art[END_REF]. This tool requires drawing the bounding box around objects in some key-frames and provides linear interpolation to infer the bounding boxes of the same object in intermediate frames. The objects can be labeled, in our case as: "person", "person_sitting" and "occluded". In the evaluation, only "person" label are considered for instance. Tab.II summarizes mains characteristic of image sequences used in our tests.

III. FROM LIDAR POINTS TO IMAGES REGION OF INTEREST

The sliding-window is the simplest candidate generation method [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] which makes the bridge between feature extraction and classification modules by scanning the entire spatial support of the image at different scales. One of the key problems of this popular method is the huge amount of negative candidates which are potentially false alarms. Rather than applying sliding windows detection techniques over the whole image, we propose to use a range sensor (here a Lidar) to roughly localize all potential obstacles inside the field-of-view (FOV) of the camera. The amount of negative candidates are significantly reduced without missing the positive candidate. Fig. 4 shows the Lidar-based sensor fusion architecture used within our system. In this section, we present the whole process starting from the Lidar segmentation, through obstacles definition, coordinate transformation and ROIs localization. Preprocessing is the decisive phase since failures during this stage will strongly affect all the subsequent modules.

A. Lidar Data Segmentation

The Lidar data segmentation process is the primary stage for obstacle detection. An entity of interest that we call here an obstacle, can be anything: buildings, vehicles, people, etc. The goal of our system is to localize an obstacle and to classify it as a "person" or a "non-person". The localization is done by performing a segmentation on the Lidar data.

The raw Lidar data is a 3D cloud of points. A single 2D scan is a sequence of N S laser measurements S = {(θ l , d l )|l = 1, . . . , N S } where (θ l , d l ) denotes the polar coordinates of the l th scan point. Each obstacle is characterized by a cluster of Lidar points. The cluster S k is obtained by a convolution operation of the points in S with a kernel mask [-1, 1] as proposed in [START_REF] Premebida | A lidar and vision-based approach for pedestrian and vehicle detection and tracking[END_REF]. Those points whose distances to neighbor points are smaller than a given threshold Θ, belong to the same cluster. A cluster S k can be expressed as

B. Lidar Points Projection In the Camera Frame

The rigid body transformation between the Lidar L and the camera C frames is denoted as C [R, t] L where C R L is the rotation matrix and C t L is the translation vector. A point L p in the Lidar frame is projected into the camera frame as:

C p = C R L • L p + C t L (1)
The extrinsic parameter C [R, t] L can be obtained through the calibration process presented in [START_REF] Fremont | Circular targets for 3d alignment of video and lidar sensors[END_REF] or deduced from direct measurements.

C. Fisheye Camera Model

Wide-angle cameras show noticeable geometric distortions. In order to follow the pinhole projection model, it is generally desirable to remove them in many applications in computer vision. The geometric distortions include two major components: radial and tangential components. Radial distortions cause image points to be translated by a proportional amount to their radial distance from the optical center. Tangential distortions (or decentering distortions) are generally less significant than radial distortions and are mainly originated from the misalignment of the lens elements. with f x and f y the focal lengths of the camera optic. In the case of a fisheye camera, the position of the distorted point on the image is given by [START_REF] Heikkila | A four-step camera calibration procedure with implicit image correction[END_REF]:

p = ũ ṽ = p + δp + p 0 ( 2
)
where

δp = δu δv = δu (r) + δu (t) δv (r) + δv (t)
is the approximated distortions and p 0 = (u 0 , v 0 ) is the principal point of the camera. The variables δu (r) , δv (r) represent the radial distortions and δu (t) , δv (t) are the tangential distortions along the two image axes. Among different distortion models, the standard polynomial model [START_REF] Hughes | Wide-angle camera technology for automotive applications: a review[END_REF] is commonly used. In this paper, we use a third degree polynomial model:

δu (r) δv (r) = ũ(k 1 r 2 d + k 2 r 4 d + k 3 r 6 d ) ṽ(k 1 r 2 d + k 2 r 4 d + k 3 r 6 d ) (3)
δu (t) δv (t) = 2p 1 ũṽ + p 2 (r 2 d + 2ũ 2 ) p 1 (r 2 d + 2ṽ 2 ) + 2p 2 ũṽ (4)
The optimal values of (f x , f y , u 0 , v 0 ) and (k 1 , k 2 , k 3 , p 1 , p 2 ) are estimated through a calibration process [START_REF] Bouguet | Camera calibration toolbox for matlab[END_REF], [START_REF] Mei | Calibration between a central catadioptric camera and a laser range finder for robotic applications[END_REF].

D. Adaptive ROI Projection

For each detected obstacle from the Lidar data segmentation, we project three characteristic points ( L C k , L A, L A ′ ) as ( I C k , I Ã, I Ã′ ) on the fisheye image. The entities L A and L A ′ are two 3D real points that create a horizontal line through the obstacle following the projection process described in sections III-B and III-C. We denote (∆x, ∆y) = I Ã -I Ã′ . As depicted in fig. 5, the angle between the obstacle appearance and the horizontal line in the fisheye image is therefore given by µ = arctan ∆y ∆x . In the proposed approach, we do not use the geometric characteristic of the cluster as a cue to recognize the obstacle. Only the cluster size helps to determine the size of the ROI in the fisheye image: For any configuration of the Lidar/Camera sensors, one Lidar cluster defines one ROI = (centroid point, ROI width, µ). This is very convenient for an image-based people detection algorithm and the requirement of a rotation invariant method is not necessary anymore.

ROI width = D k × F rame width × fy d C k × Sensor width (5)

IV. EXPERIMENTAL RESULTS

The dataset and the evaluation protocol have been presented in section II and IV-A. Our evaluation focus on the role of the Lidar in the fusion system with two different sensors configurations as stated in section II-A.

A. Evaluation protocol

The detection system takes an image and returns bounding boxes with corresponding scores or confidence indicators. A detected bounding box A and a ground truth bounding box B form a match if they have a sufficient overlap area. In the PASCAL challenge [START_REF] Everingham | The pascal visual object classes (voc) challenge[END_REF], the overlap criterion between the two bounding box A and B is defined as t = A∩B A∪B > t 0 where t 0 is a given threshold. A value of t 0 = 0.5 is considered as reasonable and is commonly used. The protocol of evaluation is adapted from the tool used in [START_REF] Dollár | Pedestrian detection: An evaluation of the state of the art[END_REF]. As the context of heavy machines requires reducing the false detection rate, the results are presented in miss-rate against false positive per image (FPPI). Only bounding boxes with a height more than 50 pixels are considered in the evaluation. Each detected bounding box may be matched once with the ground truth and redundant detections are considered as false positives (FP). In Sect. IV-B, we plot the miss-rate versus FPPI (lower curves indicate better performance) by varying the threshold of the detector. Decreasing the threshold level results in reducing the miss-rate but will increase false positive per image. At a given FPPI, it is easy to compare the performance of different detectors. The log-average miss-rate is used to summarize detector performance. The log-average miss-rate is computed by averaging miss-rate at nine FPPI rates evenly spaced in log-space in the range 10 -2 to 10 0 (for curves that end before reaching a given FPPI rate, the minimum achieved miss-rate is used). When curves are somewhat linear in this range, the log-average miss-rate is similar to the performance at 10 -1 FPPI [START_REF] Dollár | Pedestrian detection: An evaluation of the state of the art[END_REF], [START_REF] Hussein | A comprehensive evaluation framework and a comparative study for human detectors[END_REF]. The displayed legend entries are ordered by log-average miss-rate from the worst to the best case.

B. Results

1) Configuration 1:

In terms of computational cost as expected, the combination of a Lidar with a fisheye camera improves significantly the speed of the people detection system. While the complexity of the obstacles detection in Lidar data and in the ROIs projection states are always the same, the speed of the people detection by DPM approach depends on the area of the ROIs. As a result, we noticed a speed factor of about 20 in detection time, measured using 300 images at VGA resolution on a desktop computer without any specific hardware acceleration. It is worth noting that the motions of heavy machines are normally much slower than automobiles and a detection frequency of 10Hz can be considered as enough.

The results shown in Fig. 7a highlights the performance of 3 detectors: the standard HOG of [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF], the DPM approach [START_REF] Felzenszwalb | Object detection with discriminatively trained part based models[END_REF] and the proposed DPM approach used in combination with Lidar data segmentation. The DPM approach has a much better performance in people detection on fisheye images than the HOG detector, thanks to the flexible model of the deformable parts. Moreover, introducing a Lidar in a multi-sensor system apparently get an even better result by eliminating false detections in regions where it is impossible to have an obstacle. The measurements given by the Lidar are generally precise and robust, so the risk that the Lidar misses an obstacle is minor.

2) Configuration 2: Fig. 7b shows the performance of the Lidar-Camera system compare to the Camera-only approach in configuration 2. People staying on the sides of the sensors are at the same time strongly distorted and rotated. Visualbased person detection methods which are not designed to take into account these rotations, have very bad performance. The system's miss detection rate reach 95% over the range of FPPI between 10 -2 and 1, which could be considered as totally failed. With the help of the Lidar, the ROIs corresponding to obstacle are rotated (detailed in section III-D) before launching the visual recognition. In this favorable case, the DPM approach performs much better. Most of the miss detections are due to the low resolution of the appearance on image when the person is staying at a long distance (over 5m).

There are also extreme cases where the fisheye camera gives a close to top-view figure of the person (Fig. 8). The ROIs given by the Lidar is correct but the DPM detector cannot recognize the person. These specific problems can be solve by defining a mixture-model system where dedicated DPM model is trained in function of the position of the person.

V. CONCLUSIONS

Given that the DPM approach shows very good performance in detecting and recognizing objects in perspective images, we built a fisheye image dataset and evaluated the DPM approach in the context of people detection using fisheye images for heavy machines applications. It turned out that the deformable models can handle very well the strong fisheye distortions, however requires high computation cost and is not invariant to rotation of the person's appearance. Therefore, we proposed a multi-sensor system that consists of a fisheye camera and a Lidar in order to bypass these limitations. The results are promising both in term of speed and performances. Moreover, the fusion architecture and the Lidar data segmentation algorithm used in our system are sequential and very simple. In the near future, we are planning to integrate an object tracking algorithm to improve the quality of the Lidar segmentation and to solve the occlusion problems during the people detection phase.

Figure 1 :

 1 Figure 1: Sensor configurations map. Sensor heights are h = 110cm and H = 210cm.

Figure 2 :

 2 Figure 2: Architecture of the acquisition system.

Figure 3 :

 3 Figure 3: Real image of the acquisition system setup in the configuration 1. Testing dataset Number of sequences Positive images Positive samples Configuration 1

Figure 4 :

 4 Figure 4: The Lidar-based sensor data fusion architecture.

  S k = {(θ m , d m ) | m = b k , ..., b k + n k } where n k and b k are respectively the number of points and the initial scan point of the cluster S k . Now, an obstacle is represented as a cluster S k in the Lidar frame. In practice, we are interested only in the centroid point L C k = (θ C k , d C k ) and the size D k of the cluster k. In the rest of the paper, we assume that the size of a cluster is represented as the Euclidean distance between its endpoints.

Figure 5 :

 5 Figure 5: Computation of the ROI's angle on image

Figure 6 :

 6 Figure 6: An example of the Lidar-camera people detection. (a) The Lidar data presented in Top-view angle. (b) Result of the Lidar segmentation. (c) Lidar points projected on image. The maximum range is limited to 10 m. (d) Obstacles' ROI are represented by a blue box on image. The result of the person recognition algorithm are represented with a red box inside the ROIs.

Figure 7 :

 7 Figure 7: Comparing the detection performance of the visionbased approach versus the fusion system Lidar-fisheye camera using ROC curve (a) Configuration 1.(b) Configuration 2.

Figure 8 :

 8 Figure 8: Extreme case where the DPM approach fails.

Table I :

 I Repartition of accidents cause by different heavy machine operations.

	Operation	Repartition of accidents (%)
	Static	20
	Move backward	42
	Move forward	27
	Not specified	11

EPICEA is a data base on working accidents, collaboration between CNAMTS (Caisse nationale de l'assurance maladie), CRAM (les Caisses régionales d'assurance maladie) and INRS (Institut national de recherche et de sécurité) since 1988.

Pacpus is the experimental platform on Intelligent Vehicles (IV) owned by Laboratory Heudiasyc. See http://www.hds.utc.fr/pacpus

The authors are with 1 Université de Technologie de Compiègne (UTC), France. CNRS Heudiasyc UMR 7253, France. 2 Technical center for the Mechanical Industry (CETIM), France


Rémy Belmonte
email: remybelmonte@gmail.com

Archontia C Giannopoulou
email: archontia.giannopoulou@gmail.com

Daniel Lokshtanov
email: daniello@ii.uib.no

Dimitrios M Thilikos

Structure of W_4 -immersion free graphs

whether they are published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

All grahs in this paper are undirected may have multiple edges. See the monograph by Diestel [START_REF] Diestel | Graph theory[END_REF] for undefined notation and terminology. A recurrent theme in structural graph theory is the study of graphs which exclude a fixed pattern. The notion of appearing as a pattern gives rise to various graph containment relations. Maybe the most famous example is the minor relation that has been widely studied, in particular with the fundamental results of Kuratowski and Wagner who proved that planar graphs are exactly those graphs that contain neither K 5 nor K 3,3 as a (topological) minor. Another famous such example is Wagner's theorem, which describes the structure of K 5 -minor free graphs.

The structure of graphs that exclude a fixed graph as a minor was studied in the seminal Graph Minor series of papers by Robertson and Seymour [START_REF] Robertson | Graph Minors. XXI. Excluding a non-planar graph[END_REF]. However, while the structure of graphs that exclude a fixed graph H as a minor has been extensively studied, the structure of graphs excluding a fixed graph H as a topological minor or as an immersion has not received as much attention. While a general structure theorem for topological minor free graphs was very recently provided by Grohe and Marx [START_REF] Grohe | Structure theorem and isomorphism test for graphs with excluded topological subgraphs[END_REF], finding an exact characterization of the graphs that exclude K 5 as a topological minor remains a notorious open problem. Similarly, Wollan gave a structure theorem for graphs excluding complete graphs as immersions [START_REF] Wollan | The structure of graphs not admitting a fixed immersion[END_REF]. A graph G contains a graph H as a immersion if H can be obtained from G by a sequence of vertex deletions, edge deletions and replacing edge-disjoint paths with single edges. Observe that if a graph G contains a graph H as a topological minor, then G also contains H as an immersion, as vertex-disjoint paths are also edge-disjoint. In 2011, DeVos et al. [START_REF] Devos | Minimum degree condition forcing complete graph immersion[END_REF] proved that if the minimum degree of a graph G is at least 200t then G contains the complete graph on t vertices as an immersion. In [START_REF] Ferrara | On H-immersions[END_REF] Ferrara et al. provided a lower bound on the minimum degree of any graph G in order to ensure that a given graph H is contained in G as an immersion.

In this paper, we prove a structural characterization of the graphs that exclude W 4 as an immersion and show that they can be constructed from graphs that are either subcubic or have treewidth bounded by a constant. We denote by W 4 the wheel with 4 spokes, i.e., the graph obtained from a cycle on 4 vertices by adding a universal vertex. The structure of graphs that exclude W 4 and W 5 as a topological minor has been studied by Farr [START_REF] Farr | The subgraph homeomorphism problem for small wheels[END_REF]. For the case of W 4 , he proved that these graphs can be constructed via clique-sums of order at most 3 from subcubic graphs. However, this characterization only applies to simple graphs. More recently, Robinson and Farr studied the structure of graphs that do not contain W 5 [START_REF] Robinson | Structure and recognition of graphs with no 6-wheel subdivision[END_REF] and W 6 [START_REF] Robinson | Graphs with no 7-wheel subdivision[END_REF] as a topological minor.

We conclude this section with some elementary definitions that will be required later.

Definition 1 An immersion of H in G is a function α with domain V (H) ∪ E(H), such that: • α(v) ∈ V (G) for all v ∈ V (H), and α(u) = α(v) for all distinct u, v ∈ V (H);
• for each edge e of H, α(e) is a path of G with ends α(u), α(v);

• for all distinct e, f ∈ E(H), E(α(e) ∩ α(f )) = ∅.
Definition 2 Let G, G 1 , and G 2 be graphs. Let t ≥ 1 be a positive integer. The graph G is a t-edge-sum of G 1 and G 2 if the following holds. There exist vertices

v i ∈ V (G i ) such that |E G i (v i )| = t for i ∈ [2]
and a bijection π :

E G 1 (v 1 ) → E G 2 (v 2 ) such that G is obtained from (G 1 -v 1 ) ∪ (G 2 -v 2 )
by adding an edge xy for every pair of edges e 1 and e 2 such that e 1 = xv 1 , e 2 = yv 2 , and v 2 = π(v 1 ). We say that the edge-sum is internal if both G 1 and G 2 contain at least 2 vertices and denote the internal t-edge-sum of G 1 and G 2 by G 1 ⊕t G 2 .

The (elementary) wall of height r is the graph W r with vertex set

V (W r ) = {(i, j) | i ∈ [r + 1], j ∈ [2r + 2]
} in which we make two vertices (i, j) and (i ′ , j ′ ) are adjacent if and only if either i = i ′ and j ′ ∈ {j -1, j + 1} or j ′ = j and i ′ = i + (-1) i+j , and then remove all vertices of degree 1. The vertices of this vertex set are called original vertices of the wall. A subdivided wall of height r is the graph obtained from W r after replacing some of its vertices by internally vertex-disjoint paths.

It is well known that large treewidth ensures the existence of a large wall as a topological minor. Recently, Chekuri and Chuzhoy proved that polynomial treewidth suffices to ensure the existence of a wall as a topological minor, solving a long standing open problem.

Theorem 1 [START_REF] Chekuri | Polynomial Bounds for the Grid-Minor Theorem[END_REF] Let G be a graph and r a positive integer. If G does not contain W r as a topological minor, then G has treewidth O(r 98 ).

Structure of graphs excluding W 4 as an immersion

In this section, we prove the main result of our paper, namely we provide a structure theorem for graphs that exclude W 4 as an immersion. We first show that the property of containing W 4 as an immersion is closed under 3-edge sums.

Theorem 2 Let G, G 1 , and G 2 be graphs such that G = G 1 ⊕t G 2 , with t ∈ [START_REF] Devos | Minimum degree condition forcing complete graph immersion[END_REF]. Then, G contains W 4 as an immersion if and only if G 1 or G 2 does as well.

We now provide a technical lemma that will be crucial for the proof of Theorem 4.

Lemma 1 There exists a function f such that for every integer r ≥ 60000 and every graph G that does not contain W 4 as an immersion, has no internal 3 edge-cut, and has a vertex u with d(u) ≥ 4, if tw(G) ≥ f (r), then there exist sets Z = {z 1 , . . . , z r }, S 1 , . . . , S r , and X, that satisfy the following properties:

(i) z i ∈ S i , ∀i ∈ {1, . . . , r}; (ii) z i ∈ S j , ∀i = j ∈ {1, . . . , r}; (iii) u ∈ i∈{1,...,r} S i ; (iv) ∂(S i ) ≤ 6; (v) G[S i ] is connected, ∀i ∈ {1, . . . , r}; (vi) X ∩ S i = ∅, ∀i ∈ {1, . . . , r}; (vii) For every Z ′ ⊆ Z such that |Z ′ | ≥ 7, there is a 7-flow from Z ′ to X;
Lemma 1 essentially states that large treewidth yields a large number of vertex disjoint cycles that are highly connected to each other, and an additional disjoint set that is highly connected to these cycles. However, this, together with the assumption that W 4 does not immerse in G, implies that there cannot be a large flow between a vertex of degree at least 4 and one of the cycles. We will combine this fact with the notion of important separators to obtain Lemma 2. Please refer to e.g., [START_REF] Chen | An improved parameterized algorithm for the minimum node multiway cut problem[END_REF][START_REF] Marx | Parameterized graph separation problems[END_REF] for a definition of (X, Y )-important separator.

Theorem 3 [START_REF] Chen | An improved parameterized algorithm for the minimum node multiway cut problem[END_REF][START_REF] Marx | Parameterized graph separation problems[END_REF] Let X, Y ⊆ V (G) be two sets of vertices in graph G, let k ≥ 0 be an integer, and let S k be the set of all (X, Y )-important separators of size at most k.

Then |S k | ≤ 4 k and S k can be constructed in time |S k | • n O(1) .
Theorem 3 states that the number of important separators of a certain size is bounded. The next lemma uses this fact together with Lemma 1.

Lemma 2 Let G be a graph such that G does not contain W 4 as an immersion, has no internal 3 edge-cut and has a vertex u with d(u) ≥ 4. Then the treewidth of G is upper bounded by a constant.

We are now ready to prove the main theorem of our paper.

Theorem 4 Let G be a graph that does not contain W 4 as an immersion. Then the prime graphs of a decomposition of G via i-edge-sums, i ∈ [START_REF] Devos | Minimum degree condition forcing complete graph immersion[END_REF], are either subcubic graphs, or have treewidth bounded from above by a constant.

We conclude this section by noting that Theorem 4 is that it is in a sense tight: indeed, both the fact that we decompose along edge-sums of order at most 3 and the requirement that a unique vertex of degree at least 4 is sufficient to enforce small treewidth are necessary. The fact that decomposing along internal 3 edge sums is necessary can be seen from the fact that there are internally 3 connected graphs that have vertices of degree at least 4 and yet do not contain W 4 as an immersion, e.g. a cycle where every edge is doubled.

* The first author was supported by the NFR SCOPE and the MEXT ELC projects. The research of the second author leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 267959. The fourth author of this paper was co-financed by the European Union (European Social Fund ESF) and Greek national funds through the Operational Program "Education and Lifelong Learning" of the National Strategic Reference Framework (NSRF) -Research Funding Program: ARISTEIA II.