Keywords: Surgical workflow, surgical microscope, feature extraction, video analysis, surgical process model, DTW, HMM

The need for a better integration of the new generation of Computer-Assisted-Surgical (CAS) systems has been recently emphasized. One necessity to achieve this objective is to retrieve data from the Operating Room (OR) with different sensors, then to derive models from these data. Recently, the use of videos from cameras in the OR has demonstrated its efficiency. In this paper, we propose a framework to assist in the development of systems for the automatic recognition of high level surgical tasks using microscope videos analysis. We validated its use on cataract procedures. The idea is to combine state-of-the-art computer vision techniques with time series analysis. The first step of the framework consisted in the definition of several visual cues for extracting semantic information, therefore characterizing each frame of the video. Five different pieces of image-based classifiers were therefore implemented. A step of pupil segmentation was also applied for dedicated visual cue detection. Time series classification algorithms were then applied to model time-varying data. Dynamic Time Warping (DTW) and Hidden Markov Models (HMM) were tested. This association combined the advantages of all methods for better understanding of the problem. The framework was finally validated through various studies. Six binary visual cues were chosen along with 12 phases to detect, obtaining accuracies of 94%.

I. INTRODUCTION

ver the past few years, the increased availability of sensor devices has transformed the Operating Room (OR) into a rich and complex environment. With this technological emergence, the need of new tools for better resource support and surgical assessment increases [START_REF] Cleary | The operating room of the future[END_REF]. In parallel, better management of all devices, along with improved safety is increasingly necessary. In order to better Manuscript received July 18, 2011. F. Lalys, L. Riffaud, D. Bouget and P. Jannin are with the U746 INSERM Institute and University of Rennes I, Faculté de Médecine, CS 34317, F-35043, Rennes Cedex France (phone: +33223233829; e-mail: florent.lalys@irisa.fr, david.bouget@irisa.fr, pierre.jannin@irisa.fr).

L. Riffaud is with the neurosurgical department of the University Hospital of Rennes, 35000 Rennes, France, (e-mail: laurent.riffaud@chu-rennes.fr).

apprehend theses new challenges in the context of CAS systems and image-guided surgery, recent efforts have been made on the creation of context-aware systems. The idea is to understand the current situation, and automatically adapt the assistance functions accordingly. Being able to automatically extract information from the OR such as events, steps or even adverse events would allow the management of context-aware systems, but also the postoperative generation of reports or the evaluation of surgeons/surgical tools use. The goal of CAS systems based on surgery modelling is to retrieve low-level information from the OR and then to automatically extract high-level tasks from these data [START_REF] Jannin | Surgical models for computer-assisted neurosurgery[END_REF]. Consequently, it would be beneficial for surgeons performing the surgery, allowing better systems management, automatically reporting procedures, evaluating surgeons, increasing surgical efficiency and quality of care in the OR. This is why, in the context of CAS systems, the automatic extraction of high-level tasks in the OR has recently emerged.

To design such models, there are real advantages to automating the data extraction process, mainly because manual work is time-consuming and can be affected by human bias. Automatic data extraction is now easier, thanks to the high number of sensor devices. Among all sensors, teams have recently focused on videos from cameras already employed during the procedure, such as endoscopes or microscopes, which are a rich source of information. Video images provide information that can be processed by image-based analysis techniques, and then fused with other data from different sensors for creating models capturing all levels of surgery granularity. Moreover, computer vision techniques provide complex processing algorithms to transform images and videos into a new representation that can be further used for machine learning techniques through supervised or non-supervised classification. Compared to other data extraction techniques, the use of video not only eliminates the need to install additional material in the OR, but it is also a source of information that does not need to be controlled by humans, thus automating the assistance provided to surgeons, without altering the surgical routine.

A framework for the recognition of high-level surgical tasks from video images for cataract surgeries F. Lalys, L. Riffaud, D. Bouget, and P. Jannin O Current work using surgical videos has made progress in classifying and automating the recognition of high-level tasks in the OR. Four distinct types of input data associated with their related applications have recently been investigated. Firstly, the use of external OR videos has been tested. Bhatia et al. [START_REF] Bhatia | Real-time identification of operating room state from video[END_REF] analysed overall OR view videos. After identifying 4 states of a common surgical procedure, relevant image features were extracted and HMMs were trained to detect OR occupancy. Padoy et al. [START_REF] Padoy | On-line recognition of surgical activity for monitoring in the operating room[END_REF] also used low-level image features through 3D motion flows combined with hierarchical HMMs to recognize online surgical phases. Secondly, the use of endoscope videos in Minimally Invasive Surgery (MIS) has been investigated. The main constraints in MIS range from the lack of 3D vision to the limited feedback. However, studies on the subject have recently shown that the use of videos in this context was relevant. Speidel et al. [START_REF] Speidel | Situation modeling and situation recognition for a context-aware augmented reality system[END_REF] focused on surgical assistance for the construction of context-aware systems. Their analysis was based on augmented reality and computer vision techniques. They identified two scenarios within the recognition process: one for recognizing risk situations and one for selecting adequate images for the visualisation system. Lo et al. [START_REF] Lo | Episode Classification for the Analysis of Tissue-Instrument Interaction with Multiple Visual Cues[END_REF] used vision and particularly visual cues to segment the surgical episode. They used colour segmentation, shape-fromshading techniques and optical flows for instrumenttracking. These features, combined with other low-level visual cues, were integrated into a Bayesian framework for classification. Klank et al. [START_REF] Klank | Automatic feature generation in endoscopic images[END_REF] extracted image features for further scene analysis and frame classification. A crossover combination was used for selecting features, while Support Vector Machines (SVMs) were used for the supervised classification process. Also in the context of endoscopic interventions, Blum et al. [START_REF] Blum | Modeling and Segmentation of surgical workflow from laparoscopic video[END_REF] automatically segmented the surgery into phases. A Canonical Correlation Analysis was applied based on tool usage to reduce the feature space, and resulting feature vectors were modelled using Dynamic Time Warping (DTW) and Hidden Markov Model (HMM). Thirdly, also based on videos but in the context of robotic assistance, with the Da Vinci robot, Voros and Hager [START_REF] Voros | Towards "real-time" tool-tissue interaction detection in robotically assisted laparoscopy[END_REF] used kinematic and visual features to classify tool/tissue interactions in real-time. Similarly, Reiley and Hager [START_REF] Reiley | Decomposition of robotic surgical tasks: an analysis of subtasks and their correlation to skill[END_REF] focused on the detection of subtasks for surgical skill assessment. Finally, our first work [START_REF] Lalys | Surgical phases detection from microscope videos by combining SVM and HMM[END_REF] proposed the extraction of surgical phases with microscope videos and validated it in the context of neurosurgical procedures. In this paper, we extend this approach by proposing aframework that can be adjusted (subject to further researches) to any type of surgery. The objective is to automatically detect surgical phases from microscope videos. Surgical procedures can be decomposed into four main levels of granularity: phases, steps, tasks and motions. Surgical phases are thus defined as sequences of steps performed by the surgeon at a high level of granularity. The idea of the framework is first to manually defined visual cues that can be helpful for discriminating high-level tasks. The visual cues are automatically detected by state-of-the-art image-based classifiers, obtaining a semantic signature for each frame. These time series are then aligned with a reference surgery using the DTW algorithm to recognize surgical phases. Compared to traditional video understanding algorithms, this framework extracts generic application-dependant visual cues. The combination of image-based analysis and time series classification enables high recognition rates to be achieved. We evaluated our framework with a dataset of cataract surgery videos through various cross-validation studies, and compared results of the DTW approach to the HMM classification.

II. MATHERIALS AND METHODS

A. Application-dependant visual cues

The proposed framework (Fig. 1) was created to be adapted, if needed, to other types of surgery. Therefore, five sub-systems based on different image processing tools were implemented. Each of these sub-systems is related to one type of visual cue: visual cues recognizable through colour were detected with simple histogram intersection. For shape-oriented visual cues such as object recognition, a Haar classifier was trained. For texture-oriented visual cues, we used a bag-of-words approach using local descriptors, and finally for all other visual cues we used a conventional image classification approach including a feature extraction process, a feature selection process and a supervised classification with SVM. In all cases, the features were considered to be representative of the appearance of the cues to be recognized.

Preliminary pupil segmentation step

In the context of cataract surgery, some visual cues are identifiable only inside the pupil. The regions around the pupil may therefore bias the detection and a preliminary segmentation step is needed to ensure good detection. Detecting this Region Of Interest (ROI) will allow the retrieval of more specific visual cues, consequently improving surgical process analysis and modelling. The pupil segmentation procedure can be divided into three steps and is based on the colour difference between the pupil and the remaining eye regions. The first step allows the creation of an outline mask from the input image (Fig. 2-a) transformed into the YUV colour space. Within this first step, smoothing, thresholding and morphological operations were performed to the input image, obtaining a mask (Fig. 2-b). Using the mask, the second step consists in determining circles through the image using the Hough [START_REF] Swain | Color indexing[END_REF] transform (Fig. 2-c). The third step can be considered as a verification step. Incomplete circle outlines in the image mask may occur, leading to Hough circle detection failure. In order to tackle this problem, an iterative search is performed on the mask to identify the most probable circular zone, based on both the counting of corresponding pixels and circle diameters. Following this procedure, the ROI around the patient pupil can be retrieved (Fig. 2-d

Colour-oriented visual cues

Colour is one of the primary visual features used to represent and compare visual content [START_REF] Hough | Machine Analysis of Bubble Chamber Pictures[END_REF]. In particular, colour histograms have a long history as a method for image description, and can also be used for identifying colour shade through images. Here we used the principle of histogram intersection to extract colour-oriented visual cues, by creating a training image database composed of positive and negative images. Two complementary colour spaces [START_REF] Smeulders | Contentbased image retrieval at the end of the early years[END_REF] were extracted: RGB space (3 x 16 bins) along with Hue (32 bins) and Saturation (32 bins) from HSV space. For classifying visual cues, we used a KNN classifier with the correlation distance to compare histograms composed of feature vectors.

Texture-oriented visual cues

For whole-image categorization tasks, bag-of-visualwords (BVW) representations, which represent an image as an orderless collection of local features, have recently demonstrated impressive levels of performance along with relative simplicity of use. The idea of BVW is to treat images as loose collections of independent patches, sampling a representative set of patches from the image, evaluating a visual descriptor vector for each patch independently, and using the resulting distribution of samples in descriptor space as a characterization of the image. A bag of keypoints is expressed as a histogram recounting the number of occurrences of each particular pattern in an image. Given the occurrence histograms of positive and negative regions of a training database, a classifier can be trained. Considering the objective of getting binary visual cues, we used a SVM classifier with a Gaussian kernel.

In order to find these keypoints, we tested 4 keypoints detection methods: SIFT [START_REF] Lowe | Object recognition from scale-invariant features[END_REF], SURF [START_REF] Bay | SURF: Speeded Up Robust Features[END_REF], Harris [START_REF] Harris | A combined corner and edge detector[END_REF] and STAR [START_REF] Agrawal | Center surround extremas for realtime feature detection and matching[END_REF], providing access to local image information. All of them provided a similar result, which is a sample of keypoints, though they differed radically in the methods used to obtain them and by the nature of the keypoints found. After detection, a keypoint is then described as a local, rectangular or circular, patch of the image and is represented in a formal way. They are thus represented by a descriptor vector whose length is variable and highly correlated to the chosen descriptor method. Two main descriptors are generally used: SIFT and SURF descriptors. In this study, we focused on SURF descriptors, for computational reasons. Indeed, the vector space dimension is reduced by a half (from 128 to 64) when switching from SIFT to SURF descriptors.

Shape-oriented visual cues

The presence of instruments in the surgical layout is a vital piece of information to access a lower granularitylevel in surgical process analysis. The main limitation is that instruments frequently have similar shapes and are therefore very difficult to recognize through image-based analysis only. Two methods were thus implemented: one for recognizing and categorizing a specific instrument (i.e. "specific instrument categorization"), and one for detecting the presence of a surgical instrument without being able to categorize it (i.e. "Detection of other instruments"). Specific instrument categorization: For this we used a Viola-Jones object detection framework [START_REF] Viola | Rapid real-time face detection[END_REF]. It is mainly used to detect specific objects in an image, such as human faces. We chose this approach for computational reasons, getting a robust method that minimizes computation time while achieving high detection accuracy. The basic idea is to create a feature-based classifier based on features selected by AdaBoost [START_REF] Freund | A decision-theoretic generalization of on-line learning and an application to boosting[END_REF]. AdaBoost is a method of combining weighted weak learners to generate a strong learner. Here, weak learners of the algorithm are based on the Haar-like rectangular features [START_REF] Papageorgiou | A general framework for object detection[END_REF]. It is based on comparing the sum of intensities in adjacent regions inside a detection window. Strong learners are then arranged in a classifier cascade tree in order of complexity. The cascade classifier is therefore composed of stages each containing a strong learner. The strong learners, computed during the learning stage, are optimized based on negative images composed of background, and positive images of the object. During the detection phase, a window looks through the image with different scales and positions for detecting the object. The idea is to determine, at each stage of the cascade, whether the given sub-window is not or may be the searched object. This enables a high number of subwindows to be very rapidly rejected. As a result, the false positive rate and the detection rate are the product of each rate of the different stages. This algorithm is known to work well for rigid objects. We applied it in the case of surgical instrument categorization.

Detection of other instruments: This method was meant to be able to automatically detect any instrument appearing in the microscope's field of view. For this detection we used a pixel-based approach composed of 3 steps: segmentation, description and classification. For the segmentation step, the goal was to create as many ROIs as instruments in the image. The better the ROIs around the instruments, the better the identification will be. Pre-processing operations were first performed to create a mask composed of all strong outlines. This approach is based upon the fact that there is a distinct colour difference between instruments and the background. A Gaussian blur transform, a Laplacian transform, along with threshold and dilatation operations were applied to the input image (Fig. 3-a) to create the mask. By applying a connected component method to the mask, we were able to detect and remove all small connected components which were assumed to be noise. The final mask is then clean, and only contains strong outlines (Fig. 3-b). Lastly, we retrieved the two largest remaining connected components respectively (no more than 2 instruments can be present at a same time within the surgical scene), and created a mask for each one. At this stage, these selected connected components are very likely to be the instruments. By applying these masks to the input image, we obtained two different images, each one with only a ROI of the input image (Fig. 3-c, Fig. 3-d For the description and the segmentation step, the aim was to provide a robust and reproducible method for describing the ROIs that have been isolated by the segmentation step, and then to identify the ROI as an instrument or not using the descriptors. Given the fact that the goal was to detect instruments, we were willing to extract scale and rotation invariant features. We used the same approach as that described in subsection "Texture-oriented visual cues" in II.A, with local description of features using a BVW approach, and finally a supervised classification using SVM. Fig. 4-a and Fig. 4-b show an example of the extraction of local features for the input image of Fig. 3-a. If the same instrument appears with various scales and orientations, we will be able to extract the same feature points with the same descriptors. Using this approach, ROIs can be classified as belonging to an instrument or not.

Alternative method

This approach can be seen as conventional image classification. Each frame is represented by a signature composed of global spatial features. The RGB and the HSV spaces, the co-occurrence matrix along with Haralick descriptors [START_REF] Haralick | Textural features for image classification[END_REF], the spatial moments [START_REF] Hu | Visual pattern recognition by moment invariants[END_REF], and the Discrete Cosine Transform (DCT) [START_REF] Ahmed | Discrete Cosine Transform[END_REF] coefficients were all computed. Each signature was finally composed of 185 complementary features, that had to be reduced by feature selection. We combined a filter and a wrapper approach [START_REF] Duda | Pattern classification and scene analysis[END_REF] using the method described by Mak and Kung [START_REF] Mak | Fusion of feature selection methods for pairwise scoring SVM[END_REF] using the union of both results for improving the final selection. The Recursive Feature Elimination (RFE) SVM [START_REF] Guyon | Gene selection for cancer classification using support vector machine[END_REF] was chosen for the wrapper method and the mutual information (MI) [START_REF] Hamming | Coding and Information Theory[END_REF] was chosen for the filter method. Finally, 40 features were kept, and a SVM was applied to extract the desired visual cue. This particular procedure has been presented and validated in Lalys et al. [START_REF] Lalys | Automatic phases recognition in pituitary surgeries by microscope images classification[END_REF].

Once every visual cue of a particular surgery has been defined and detected, creating a semantic signature composed of binary values for each frame, the sequences of frame signatures (i.e. the time series) must be classified using appropriate methods. Two different approaches were tested here: the HMM modelling and the classification by DTW alignment.

Hidden Markov Model

HMMs [START_REF] Rabiner | A tutorial on Hidden Markov Models and selected applications in speech recognition[END_REF]

Dynamic Time Warping

We used the DTW algorithm [START_REF] Keogh | An enhanced representation of time series which allows fast and accurate classification, clustering and relevance feedback[END_REF] as a method to classify the image sequences in a supervised manner. DTW is a well-known algorithm used in many areas (e.g. handwriting and online signature matching, gesture recognition, data mining, time series clustering and signal processing). The aim of DTW is to compare two sequences X := (x1 , x2 , . . . , xN) of length N and Y := (y1 , y2 , . . . , yM) of length M. These sequences may be discrete signals (time-series) or, more generally, feature sequences sampled at equidistant points in time. To compare two different features, one needs a local cost measure, sometimes referred to as local distance measure, which is defined as a function. Frequently, it is simply defined by the Euclidean distance. In our case however, considering that we are using binary vectors, we will use the Hamming distance, well-adapted for this type of signature. In other words, the DTW algorithm finds an optimal match between two sequences of feature vectors which allows for stretched and compressed sections of the sequence. To compare each surgery, we created an average surgery based on the learning dataset using the method described by Wang and Gasser [START_REF] Wang | Alignment of curves by dynamic time warping[END_REF]. Each "query" surgery is first processed in order to extract semantic information, and then the sequence of image signature is introduced in the DTW algorithm to be compared to the average surgery. Once warped, the phases of the average surgery are transposed to the unknown surgery in a supervised way. Additionally, global constraints (also known as windowing functions) can be added to the conventional algorithm in order to constrain the indices of the warping path. With this method, the path is not allowed to fall within the constraints window. For surgery modelling, we used the Itakura parallelogram [START_REF] Niennattrakkul | Learning DTW global constraint for time series classification[END_REF]. This prevents the warping path from straying too far away from the diagonal path.

B. Data Video data-set

Our framework was evaluated on cataract surgeries (eye surgeries). The principle of cataract surgery is to remove the natural lens of the eye and insert an artificial one (referred as an IntraOcular Lens, or IOL) in order to restore the lens' transparency. In this project, twenty cataract surgeries from the University Hospital of Munich were included (with mean surgical time of 15 min). Three different surgeons performed these procedures. All videos were recorded using the OPMI Lumera surgical microscope (Carl Zeiss) with an initial resolution of 720 x 576 at 25 fps. Considering the goal of recognizing only high-level tasks surgical tasks, each video was downsampled to 1 fps. Original frames were also spatially downsampled by a factor 4 with a 5-by-5 Gaussian kernel. Twelve surgical phases were defined (Fig. 6.).

Definition of the visual cues

Six pieces of binary visual cues were chosen for discriminating the surgical phases. The pupil colour range, defined as being orange or black, was extracted using a preliminary segmentation of the pupil along with a colour histogram analysis. Also analysing only the pupil after the segmentation step, the global aspect of the cataract (defined as parcelled out or not) was recognized using the BVW approach with local spatial descriptors. The presence of antiseptic, recognizable by virtue of its specific colour, was detected using colour histogram analysis, but on the entire image. Concerning the detection of surgical instruments, only one had a characteristic shape, the knife. We trained a Haar classifier using 2000 negative images and 500 positive images for its detection. All other instruments have very similar shapes and are very difficult to categorize. For this reason, we chose to detect the presence of an instrument as a particular visual cue. Lastly, the IOL instrument was not readily identifiable through only colour or shape analysis and we chose a classical approach using many spatial features along with a SVM classification to detect this visual cue.

C. Validation studies

Initial indexing was performed by surgeons for each video, by defining the phase's transitions, along with all visual cues for each video. With this labelled video database, we were able to evaluate both aspects of our framework, i.e. detection of the different visual cues and the global recognition rate of the entire framework. From each video, we randomly extracted 100 frames in order to create the image database, finally composed of 2000 labelled images, each being associated to one video. This image database was used for the assessment of visual cue detection, whereas the videos, along with their corresponding frames from the image database, were used to assess the entire framework.

One step of our procedure didn't require any training stage: the preliminary step of pupil segmentation. This step was simply validated over the entire video database by testing each frame of each video. During this validation, a pupil was considered correctly segmented if and only if the circle found by the Hough transform precisely matched the pupil. A percentage was then obtained corresponding to the accuracy of the segmentation.

The second aspect of our framework that was assessed was the recognition of all visual cues. Independently, for knife recognition, the training stage was performed using manually selected positive and negative images for better object training. For the validation, 1000 test images were used and global recognition accuracy was computed. The four other visual cues classifiers were the assessed through 10-fold cross-validation studies. The image database was therefore divided into 10 subsets, randomly selected from the 20 videos. Nine were used for training while the prediction was made on the 10 th subset. One test subset was consequently composed of frames from two videos, while the validation sets was composed of the rest of the frames from the 18 other videos. This procedure was repeated 10 times and results were averaged. After evaluating each classifier, we validated their use compared to a traditional image-based classifier, i.e. compared to feature extraction, selection and classification as performed by the conventional classifier. Additionally, we also compared the four local keypoint detectors presented in subsection "Texture-oriented visual cues" in II.A, along with the optimal number of visual words, for both the textureoriented classifier and the detection of instrument presence. This validation, performed under the same conditions as the detection of the visual cues, was necessary to optimize the detection of texture and shape-oriented information by the BVW approach.

Lastly, we evaluated the global framework, including visual cue recognition and the time series analysis with the same type of cross-validation. Similarly to the visual cues recognition, at each stage, 18 videos (and their corresponding frames from the image database) were used for training and recognitions were made on the 2 others. For this assessment, the criterion chosen was the Frequency Recognition Rate (FRR), defined as the percentage of frames correctly recognized over a video by the recognition framework.

III. RESULTS

Taking all frames from each video (at 1 fps), the pupil was correctly extracted with an accuracy of 95% (Tab. 1.). The worse video was very difficult to segment, with 78% of all frames correctly segmented. The best video, on the other hand, had almost its entire frame correctly segmented (99%). Fig. 7-a,b. shows the BVW study for choosing the best parameters for both the detection of instrument presence and the texture-oriented classifier respectively. Surprisingly, for both figures, the number of visual words did not appear to be a major parameter to be enhanced. Indeed, the accuracy didn't vary significantly from 1 to 60 visual words, and this result was true for the 4 keypoint detectors and the two BVW studies. For both studies, the best accuracy was still obtained for a number of visual words equal to 12. On the contrary, the influence of keypoint detectors was significant. For the detection of instruments presence (Fig. 7-a), the SURF keypoint detector showed best recognition accuracies (with 12 visual words: 86%), whereas for the detection of the cataract aspect, the SIFT keypoint detectors shows best results (with 12 visual words: 83%). Tab. 2. gives all parameter values used in the 6 detection methods of visual cues.

The results of the cross-validation study for the recognition of all visual cues (Tab. 3.) showed that very good detection accuracies were obtained, which outperform the standard classifier in all cases. The best recognition was obtained for the presence of the Knife (Haar classifier), achieving a recognition rate of 96.7%, whereas the lower rate was obtained for the recognition of the instrument presence. Similarly, the other visual cue detected using a BVW approach (the aspect of the cataract) has not a very high accuracy (87.2). Histogram approaches shows good results (96.2% for the pupil colour range detection and 96.1% for the antiseptic detection), whereas the IOL instrument has also a good recognition rate of 94.6%, even detected with the conventional classifier. From Tab. 4., the time series study showed better results using the DTW approach than with HMM classification. With HMM, a mean FRR of 91.4% was obtained with a high maximum of 99.7%, and a low minimum of 80.5%. The DTW approach shows a mean FRR of 94.4% with quite high minimum (90.6%) and particularly low maximum (96.4) compared to its mean FRR. Fig. 8. shows an example of video recognized by the system (by HMM classification and with the DTW approach respectively) along with the warping path from the DTW approach. On this particular recognized video, DTW shows the best recognition compared to the HMM classification.

IV. DISCUSSION

In this paper, we proposed a framework that automatically recognizes surgical phases of cataract surgery. The recognition relies on data provided by microscope videos, which is a novel way of addressing situation recognition issues. Twelve phases were defined by surgeons, and their recognition was based on the detection of 6 visual cues, along with time series analysis. This combined approach allowed a high degree of automatic recognition system accuracy, achieving accuracies of approximately 94% using DTW modelling.

A. Adaptation to other surgical procedures

The detection framework is based on the recognition of visual cues within the surgical scene. The recognition of such information allows the modelling of the surgical procedure and final surgical phase recognition. Due to the varying facilities between departments, the numbers of phases, as well as the colours, tools and shapes will differ. Consequently, considering that surgical environments are different in each hospital, one recognition system should be created for each department. The adopted solution was to create a framework using specific image-based sub-systems in order to be as generic as possible, and to provide as many tools as possible for being exhaustive. This way, our method addresses the issue of system adaptability.

Even though the framework was created to be adaptable, each kind of surgical environment has its own particularities and characteristics. This is why preliminary pre-processing steps may be mandatory in order to tune the recognition system according to the type of surgery. For instance, in the case of low-resolution video images, the purpose would be to improve image quality for further processing. In the context of cataract surgery, the microscope field of view is precisely delineated, thus enabling the use of a preliminary step of segmentation to restrict the search for a specific visual cue within a ROI defined by the pupil outlines. This is the only step that is specific to cataract surgery. For the adaptation to other surgical procedures, this segmentation step could be either adapted or even removed. Taking as example neurosurgical procedures and specifically hypophyse surgeries that we used in our previous work [START_REF] Lalys | Surgical phases detection from microscope videos by combining SVM and HMM[END_REF], segmentation would not be necessary as the field-of-view is already zoomed and adapted for image-based analysis and visual cues extraction. Pre-processing steps for image quality enhancement would also not be required because of the high-resolution of neurosurgical microscopes, neither intensity corrections nor specular reflection removal. This example would be true for this specific type of easy and reproducible surgical procedures. However, other adaptations could be conducted. Dealing with more complex surgeries would involves further researches on the pre-processing step, on the segmentation of surgical tools before their categorizations and possibly on the definition of other sub-systems for the detection of visual cues. We proposed in this paper a complete framework that we tested on one surgical procedure, but the evolution to other surgical procedures should be experimented.

Once the framework has been tuned a dedicated surgical procedure, its use is fully automatic and will work with any microscope video of this type of surgery in its environment. Similarly, other variability factors may affect recognition, such as the manner in which surgeons operate. With this system, a training stage is necessary for each surgical department, assuming that the surgeons within the department use identical materials and follow the same sequence of phases during the procedure, making image features invariant to task distortion.

B. Pupil segmentation

Using an adapted method composed of image-based analysis, the segmentation of the pupil returns highly accurate results. For 95% of the frames, the ROI correctly contains the entire pupil. Moreover, to avoid distorting any further detection that could be done within the pupil we decided to define a constant circumference value. Thus, each time a region of interest is detected, the centre is kept and the circumference value is reset to the default value. Due to its high accuracy over the entire video database, it allows all potential colour-associated noise to be removed from around the pupil for further recognition. The very low accuracy obtained for one video can be explained by the presence of the retractors, rendering the field of view very narrow. Automatic segmentation turns out to be difficult when the retractors, or even the surgical instruments, around the eye, occupy too much space within the field of view.

As a drawback, our approach, which always returns a ROI, is not always perfectly centred on the middle of the pupil. We can explain this issue by the fact that the pupil is not always completely inside the microscope's field of view. Sometimes the pupil outlines are too distorted due to surgical tools or the surgery itself. Sometimes the retractor is as wide as the pupil and sometimes the surgeon's fingers are in the field of view. In that case, it's difficult to extract the exact position of the pupil and its outlines and to adjust an intensity threshold accordingly. If the surgical microscope had a permanent position, or if we could precisely estimate the position of the pupil in each image, it would be possible to adjust a threshold for the segmentation automatically.

C. Application-dependant visual cues

Before the visual cue recognition training stage, the user will need to choose the visual cues and the associated image-based recognition classifier. In image classification problems, users usually do not think in terms of low-level features, resulting in poor recognition of the high-level semantic content of the images. Here, during the stage of visual cue definition, the colour, texture and shape behaviour of the visual cues are often intuitively known, allowing the most effective classifiers to be chosen. When visual cue is unknown or undocumented, the solution proposed is to choose the generic approach, integrating a large number of image features. This approach, combining global spatial features and SVM, may therefore be adapted to the recognition of any type of cue. The feature selection step allows the user to select discriminatory features and remove unsuitable ones, which is the intended objective. To improve recognition, however, the three other specific classifiers seem to be well-adapted when the behaviour of the visual cue is well perceived.

Generally, the main drawback of a global colour histogram representation is that information concerning object shape, and texture is discarded. In our case, however, it is only used for colour-related visual cue detection. Similarly, the main drawback of shape-based approaches is the lack of well-defined outlines. The Haar classifier was used in our framework for specific objects only, e.g. the knife, which is significantly different from all others instruments used in cataract surgery. The use of this approach to categorize other instruments, such as the cannula, was tested, but gave very poor results due to the narrow field of view and the difficulty in discriminating that specific instrument from the others. For this reason, we chose to use a second approach for object recognition, allowing the system to gain an information concerning object presence, without categorizing it. This type of information was still relevant for phase detection and allowed complete image signatures to be achieved using other information with a different level of granularity. The use of a BVW approach combined with local descriptors was also validated. Local descriptor comparisons (Fig. 5.) enabled selection of the most appropriate features, and application with the recognition of the global aspect of the cataract gave very promising results.

With the exception of the Haar classifier, the three other classifiers are all based on a training image database. The power of discrimination of the image database is thus vital. We can easily imagine that accuracy may decrease sharply if the images do not efficiently represent all phases or all scene possibilities within the phases. Additionally, the training stage is time-consuming and requires human efforts. In our particular case, the best method, used here in our validation studies, was to annotate surgical videos before randomising the initial sample. The randomisation process is thus no longer performed on all frames, but on each video independently, extracting the same number of frames per video.

D. Time series analysis

Combined with state-of-the-art computer vision techniques, time series analysis displayed very good performance, opening the way for further promising work on high-level task recognition in surgery. Without this step of time series analysis, results are less accurate (~80%). This can be explained because some visual cues don't appear only during one particular phase, and the information of sequentiality is needed. For instance, the knife always appears twice during cataract surgery: once during phase n°4 (principal corneal incision), and once during phase n°10 (expansion of the principal incision). All other visual cues are not present during these two phases. The discrimination of both phases appears to be possible with an information of time only that the HMM or the DTW can bring.

In particular, DTW captures the sequentiality of surgical phases and is well-adapted to this type of detection. The cost function between 2 surgical procedures with the same sequence of phases, but with phase time differences, will be very low. The advantage is that it can accurately synchronize two surgical procedures by maximally reducing time differences. The main limitation concerning the use of DTW, however, can be seen phase sequence differences appear between two surgeries. The warping path would not correctly synchronize the phases and errors would occur. In the context of cataract surgery, the procedure is standardized and reproducible, justifying the very good results of the recognition. But we can imagine that for surgeries that are not completely standardized, DTW would not be adapted. In this case, HMM could be used by adding bridges between the different states of the model (the transition matrix should be adapted in that case), allowing the sequence to be resumed and perform the same phase multiple times. The state machine would not be a left-right structure but would include more complex possibilities with many bridges between states. As a drawback, the complexity of such HMMs could completely affect the recognition accuracy. For each surgical procedure, the HMM structure should be created by minimizing the possibilities of transitions from states to states not to affect the classification phase. In the particular case of cataract surgery, the results showed that the DTW algorithm was, not surprisingly, quite better than HMM. Compared to HMM, the other limitation of the DTW algorithm is that it cannot be used on-line, as the entire procedure is required to determine the optimum path. For on-line applications, the HMM classification should be used.

E. Use in clinical routine

The automatic recognition of surgical phases is very useful for context-awareness applications. For instance, there is a need for an analysis methodology specifying which kind of information needs to be displayed for the surgeon's current task. The use of videos, such as endoscope or microscope videos allows automating the surgeons' assistance without altering the surgical routine. Moreover, it can also support intra-operative decision making by comparing situations with previously recorded or known situations. This would result in a better sequence of activities and improved anticipation of possible adverse events, which would, on the one hand optimize the surgery, and on the other hand improve patient safety. Because of the time-series methods used in the framework that are not fully on-line algorithms, the on-line use of the system remains a long-term objective. In its present form, the computation time of the recognition process for one frame (visual cues detection + DTW/HMM classification) was evaluated to around 3s on a standard 2-Ghz computer. It could be introduced into clinical routine for post-operative video indexation and creation of pre-filled reports. Surgical videos are increasingly used for learning and teaching purposes, but surgeons often don't use them because of the huge amount of surgical videos. One can imagine a labelled database of videos with full and rapid access to all surgical phases for easy browsing. The video database created would contain the relevant surgical phases of each procedure for easy browsing. We could also imagine the creation of pre-filled reports that would need to be completed by surgeons.

V. CONCLUSION

The goal of context-aware assistance is to collect information from the OR and automatically derive a model that can be used for advancing CAS systems. To provide this type of assistance, the current situation needs to be recognized at specific granularity levels. A new kind of input data is more and more used for detecting such highlevel tasks: the videos provided by existing hardware in the OR. In this paper, our analysis is based exclusively on microscope video data. We proposed a recognition system based on application-dependant image-based classifiers and time series analysis, using either an HMM or DTW algorithm. Using this framework, we are now able to recognize the major surgical phases of every new procedure. We have validated this framework with cataract surgeries, where twelve phases were defined by an expert, as long as six visual cues, achieving a global recognition rate of around 94%. This recognition process, using a new type of input data, appears to be a non-negligible progression towards the construction of context-aware surgical systems. In future work, lower-level information, such as the surgeon's gestures, will need to be detected in order to create more robust multi-layer architectures.

Fig. 1 .

 1 Fig. 1. Framework of the recognition system.

Fig. 2 .

 2 Fig. 2. Different steps of the pupil segmentation. a) input image, b) 1 st step: creation of the mask, c) 2 nd step: Hough transform computation, d) 3 rd step: final segmentation of the pupil.

Fig. 3 .

 3 Fig. 3. Different stages of the segmentation step for the detection of instruments a) the input image, b) the clean mask, c) the region of interest corresponding to the first connected component, d) the ROI corresponding to the second connected component.

Fig. 4 .

 4 Fig. 4. SURF features detected on image from Fig. 3-a., and shown as blue circles. a) SURF points on the first connected component, b) SURF points on the second connected component.

1 .

 1 Here, outputs of the supervised classification were treated as observations for the HMM. States were represented by the surgical phases that generated a left-right HMM (Fig.5.). Transition probabilities from one state to its consecutive state were computed for each training video, and then averaged. If we set one probability to α , the probability of remaining in the same state is then α Output probabilities were computed as the probability of making an observation in a specific state. Training videos were applied to the supervised classification for extracting binary cues and output probabilities were obtained by manually counting the number of occurrences for each state. Then, given the observations and the HMM structure, the Viterbi algorithm[START_REF] Viterbi | Errors bounds for convolutional codes[END_REF] identifies the most likely sequence of states. HMM training, like feature selection process, must be performed once only for each learning database.

Fig. 5 .

 5 Fig. 5. Left-right HMM, where each state corresponds to one surgical phase

Fig. 6 .

 6 Fig. 6. Typical digital microscope frames for the 12 surgical phases: 1-preparation, 2-betadine injection, 3-lateral corneal incision, 4principal corneal incision, 5-viscoelastic injection, 6-capsulorhexis, 7-phacoemulsification, 8-cortical aspiration of the big pieces of the lens,

Fig. 7 .

 7 Fig. 7. BVW validation studies comparison of accuracies with different number of visual words and different keypoints detectors: a) Detection of the instruments presence, b) Recognition of the cataract aspect.

Fig. 8 .

 8 Fig. 8. Distance map of two surgeries and dedicated warping path using the Itakura constraint (up), along with transposition of the surgical phases (middle), and the visual cues detected by the system (down).

Table 1 .

 1 Mean accuracy, minimum and maximum of the segmentation of the pupil over the entire video database.

		Accuracy (Std)	Minimum	Maximum
	Detection	95,00 (6)%	78,00%	99,00%

Table 2 .

 2 Parameters of the classification algorithms used for extracting visual cues.

	Type of visual cues	Algorithm	Parameters
			Type of color space: RGB,
	Color-oriented	Color histogram intersection	HSV Classifier: KNN
			Distance: correlation
			Classifier: SVM with Gaussian
			kernel
	Texture-oriented	BVW approach	Interest points detectors: SIFT
			Feature representation: SURF
			Codebook generation: KNN
	Instrument categorization	Viola-Jones approach	Features: Haar-like rectangular Negative images: 2000 Positive images: 500
			Classifier: SVM with Gaussian
			kernel
	Detection of other instruments	BVW approach	Interest points detectors: SURF
			Feature representation: SURF
			Codebook generation: KNN
			Spatial features: RGB, HSV
			spaces, Haralick descriptors,
	Alternative method	Global features classification	DCT, spatial moments Wrapper method: RFE-SVM Filter method: MI
			Classifier: SVM with Gaussian
			kernel

Table 3 .

 3 Mean accuracy (standard deviation) for the recognition of the 6 binary visual cues.

		Pupil colour	Presence	Presence	Presence	Cataract aspect	Presence
		range	antiseptic	Knife	IOL instrument		instrument
	Specific image-based classifier (%)	96.2 (3.6)	96.1 (0.7)	96.7 (3.4)	94.6 (1.1)	87.2 (5.4)	84.1 (8.6)
	Classical approach (%)	94.1 (4.6)	95.6 (0.4)	88.5 (4.3)	X	54.1 (3.6)	58.7 (6.1)

Table 4 .

 4 Mean, minimum and maximum FER of the DTW

		Accuracy (Std)	Minimum	Maximum
	HMM (%)	91.4 (6.4)	80.5	99.7
	DTW (%)	94.4 (3.1)	90.6	96.4


Angelo Iollo

Damiano Lombardi

Advection Modes by Optimal Mass Transfer

Classical model reduction techniques approximate the solution of a physical model by a limited number of global modes. These modes are usually determined by variants of principal component analysis. Global modes can lead to reduced models that perform well in terms of stability and accuracy. However, when the physics of the model is mainly characterized by advection, the non-local representation of the solution by global modes essentially reduces to a Fourier expansion. In this paper we describe a method to determine a low-order representation of advection. This method is based on the solution of Monge-Kantorovich mass transfer problems. Examples of application to point vortex scattering, Korteweg-de Vries equation, Von Kármán wake and hurricane Dean advection are discussed.

Introduction

Principal Component Analysis (PCA) is the main tool behind many techniques to perform model reduction of systems of partial differential equations (PDEs). The objective of model reduction is to obtain a simpler model that retains the main dynamical features of the full model, typically for optimization and control purposes. The principle of PCA is to find the base of a small dimensional subspace in such a way that the solution of a given PDE is accurately represented in this subspace. This is the main idea behind Proper Orthogonal Decomposition (POD) [START_REF] Holmes | Turbulence, coherent structures, dynamical systems and symmetry[END_REF][START_REF] Sirovich | Low dimensional description of complicated phenomena[END_REF]: for a given database of model solutions, POD * angelo.iollo@math.u-bordeaux1.fr, Univ. Bordeaux, IMB, UMR 5251, F-33400 Talence, France. CNRS, IMB, UMR 5251, F-33400 Talence, France. INRIA, F-33400 Talence, France.

† damiano.lombardi@inria.fr, INRIA, Paris Rocquencourt, B.P.105, 78153 Le Chesnay

1 X ξ x Y Ω 0 Ω
Figure 1: Lagrangian description of transport: the reference configuration is Ω 0 , points ξ ∈ Ω 0 are transported by the direct mapping in X(ξ, t). Given the actual configuration Ω, a point x ∈ Ω is sent back to its counterimage in the reference configuration by backward characteristics, i.e., the inverse mapping Y (x, t).

builds a basis that minimizes the L 2 average distance between the reduced representation and the solution database. POD approximations are usually satisfactory for those problems where the solution has a global behavior, but tends to perform poorly for systems characterized by concentrated structures that are advected.

Alternatives to POD taking implicitly into account the notion of transport in the definition of a global reduced basis are provided by Koopman modes [START_REF] Mezić | Spectral properties of dynamical systems, model reduction and decompositions[END_REF][START_REF] Rowley | Spectral analysis of nonlinear flows[END_REF] and Dynamic Mode Decomposition (DMD) [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]. These approaches assume the existence of a linear propagator (relative to some possibly non-linear map) whose spectrum provides a frequency based mode decomposition of the considered evolution. Starting from snapshots of some observables of the physical system, an Arnoldi-type algorithm is employed to estimate the linear propagator.

Our objective here is to introduce the notion of transport modes to represent by a reduced basis the evolution of coherent structures when advection is the leading phenomenon. To plainly present the main ideas let us proceed in two steps.

Transport. In Fig. 1 a conceptual description of transport is shown. Given a point ξ ∈ Ω 0 , where Ω 0 ⊆ Ê d is a reference configuration, transport at time t is described by a mapping X(ξ, t). The point x = X(ξ, t) belongs to the actual physical configuration Ω ⊆ Ê d . Let us consider a point x in the actual physical configuration. The inverse mapping, denoted by Y (x, t) (called otherwise backward characteristics), identifies the point in the reference configuration that has been transported by the direct map in x at time t. The following relations hold:

x = X(ξ, t), ξ = Y (x, t), Y = X -1 , [∇ ξ X][∇ x Y ] = I, (1)
where [∇ ξ X] is the jacobian of the transformation X(ξ, t) and [∇ x Y ] its inverse, i.e., the jacobian of the inverse mapping. Also, we have:

∂ t Y + v • ∇ x Y = 0, Y (x, 0) = x v(x, t) = ∂ t X, X(ξ, 0) = ξ, (2)
where v is the velocity field.

Let us consider, as an example, the inviscid Burgers equation:

∂ t v + v • ∇ x v = 0. (3)
This equation describes a pressureless Euler flow. Since no force is acting on the medium, each component of the velocity field is purely advected. In lagrangian coordinates we have:

∂ 2 t X(ξ, t) = 0 =⇒ X(ξ, t) = ξ + v(ξ, 0)t. ( 4
)
The solution consists of particles moving on straight lines (no acceleration).

Advection modes. An optimal subspace of given size to represent a solution database is usually sought by minimizing the distance (to be defined) between the solution database and its projection in the subspace. Let u(x, t) represent the solution database, where x and t stand here for discrete time and space. With POD, the distance to be minimized is defined by the norm of an appropriate Hilbert space, and the projected database is represented by

û(x, t) = a i (t)φ i (x), i = 1, ..., N P OD , (5)
where here and in the following we use the summation convention. The function φ i is called the i-th POD mode, and it is computed as a linear combination of the elements of the database, i.e., the snapshots. The solution u(x, t) is approximated by global fields (defined everywhere in the actual physical configuration) which are typically oscillating in time.

Let now ũ be another approximation of the database such that:

ũ(x, t) = [u 0 (Y (x, t)) + R(Y (x, t), t)] det(∇ x Y (x, t)), (6)
where u 0 (ξ) is a reference mode, Y (x, t) is a suitable backward mapping and R(ξ, t) is a residual to be defined in the following. This relation states that the database elements can be represented as the sum of a mapped constant reference configuration u 0 (ξ), and the mapping of a time-dependent residual field R(ξ, t) in the reference configuration.

Three elements have to be defined in order to obtain a reduced representation: a reference mode u 0 (ξ), a reduced-order expansions for the inverse characteristics Y (x, t)

and the residual R(ξ, t).

The following ansätze are introduced:

R(ξ, t) = α k (t)φ k (ξ), k = 0, ..., N r Y (x, t) = β j (t)W j (x), j = 0, ..., N y , (7)
where φ k are POD modes for the residual R and W j are the advection modes to be found.

The rest of this paper explains how to extract the reference mode, the advection modes as well as residual modes starting from a snapshots database u i = u(x, t i ), i = 1, ..., N s .

Advection mode expansion

In order to determine the advection modes, we define a suitable optimal transport problem. Let us associate a scalar density function ̺(u) ≥ 0 to the solution u(x, t) , in such a way that:

Ω ̺(x, t) dx = 1, ∀t ∈ Ê + (8)
so that the non-negative density is normalized to 1 for all times. The choice of the density function is for the moment arbitrary. If u is a non-negative scalar and satisfies this normalization, it may be directly used as a density function.

Let ̺ i , i = 1, ..., N s be the snapshots of the density function. The Wasserstein distance (denoted by W) relative to a density pair is defined as:

W 2 (̺ i , ̺ j ) = inf X Ω ̺ i (ξ)| X(ξ) -ξ| 2 dξ , subject to ̺ i (ξ) = ̺ j ( X(ξ)) det(∇ ξ X). (9)
The Wasserstein distance is the distance induced by the optimal mapping X * that minimizes the cost of the L 2 transport (Monge) problem, among all the change of coordinates X(ξ) locally keeping constant mass between the densities i and j. The solution to this problem exists and is unique (see, for a comprehensive review [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Villani | Optimal Transport, old and new[END_REF][START_REF] Evans | Partial differential equations and Monge-Kantorovich mass transfer[END_REF]).

The squared Wasserstein distance can be computed for all i, j = 1, ..., N s : 1 2 N s (N s -1) Monge problems are solved. Then, the following matrix is defined:

D ij = W 2 (̺ i , ̺ j ), (10)
that is the matrix of the squared distances between the densities. This matrix is symmetric

(W(̺ i , ̺ j ) = W(̺ j , ̺ i )
) and all the elements on the diagonal are zero (W(̺ i , ̺ i ) = 0).

In an Euclidean vector space, one can uniquely transform a matrix of canonical squared distances relative to a set of points, in a positive semi-definite matrix whose entries are the corresponding scalar products between the position vectors χ i of those points. In other words

D ij = χ i -χ j 2 = B ii + B jj -2B ij (11)
where B ij = χ i • χ j are the entries of a matrix B ∈ Ê AE××AE× . Let χ ∈ E ⊆ Ê AE× and let ½ ∈ Ê Ns be a column vector whose components are all 1. We assume that the vectors representing the points are centered in the origin of the vector space, i.e., that B ½ = 0.

Then Eq. ( 11) can be inverted and we have

B = - 1 2 JDJ (12)
and

J = I - 1 N s ½½ T , (13)
where I is the identity matrix. We assumed that E is embedded in an Euclidean space. In this case it can be shown that Eq. 12 leads to a matrix B that is positive semidefinite. For a Wasserstein distance matrix this is not necessary the case. However, one can look for an euclidean set of vectors that gives a squared distance matrix which is an approximation of the Wasserstein square distance matrix.

In order to do this, we adopt the same technique presented in [START_REF] Muskulus | Wasserstein distances in the analysis of time series anddynamical systems[END_REF]. Matrix B is real and symmetric and hence it can be decomposed as:

B = USU T , ( 14
)
where U is a unitary matrix and S is the diagonal matrix whose entries are the eigenvalues As stated, the embedding is performed with respect to the origin of the euclidean space, which corresponds to a density distribution whose properties will be investigated in the following sections. The ith component of the kth eigenvector of B + represents the weight of the ith mapping (that transports the origin in the ith snapshot) to build the optimal transport corresponding to the kth eigenvalue. In the following we define the kth transport mode as the kth optimal transport, once a barycentric density is given. Conversely, when a base of mappings is taken in the space {v 1 , ..., v k }, the mapping that transports the barycenter in the ith snapshot is derived by summing the base mappings multiplied by the coordinates of the point representing the ith snapshot.

In POD the snapshots reduced space is generated by the snapshots themselves, while for transport modes the representation of the snapshots is provided by a set of optimal transports that map a barycentric density into the snapshots. This barycentric density may be used as u 0 (ξ). The mappings obtained in this way are the basis of transport mappings. Let us remark that two series of optimal mappings are obtained: the direct mappings and their corresponding inverses.

Residual representation

In this section the representation of the function residual function R is discussed. The transport modes are defined according to the technique presented in the section above.

The entire set of transport modes account for an exact representation of the densities.

A subset of N T transport modes is chosen, in order to provide a reduced representation of the direct and inverse mapping, respectively X(ξ, t) and Y (x, t):

X(ξ, t i ) = α j (t i )Z j (ξ), Y (x, t i ) = β j (t i )W j (x), j = 0, ..., N T ≤ N s , (16)
where the coefficient α j (t i ) is known from the components of the j-th eigenvector of the embedding matrix B and the coefficient β j may be derived by considering the definition of Y as the inverse mapping.

Then we map each snapshot to the reference domain and compute the residual with respect to the initial condition

R i (ξ) = R i (X(ξ, t i )) = u(X(ξ, t i ), t i ) det(∇ ξ X(ξ, t i )) -u 0 (ξ), (17)
where u(x, t i ) = u(X(ξ, t i ), t i ) are the solution snapshots. Then the classical POD technique is applied:

A ij = Ω 0 R i (X(ξ, t i ))R j (X(ξ, t j )) dξ, (18)
and the POD residual modes are defined through the eigenvalues and eigenvectors of the autocorrelation matrix:

φ k (ξ) = h b (k) h R h (ξ) λ 1/2 k , (19) where b (k)
h is the h-th component of the k-th eigenvector of A and λ k its corresponding eigenvector.

We finally obtained a representation in a reference configuration of a time dependent field, decomposed in a mean field (the barycentric density or equivalently the initial condition) and a POD expansion. This time dependent field is mapped to an actual configuration at time t i by the mapping Y (x, t i ), determined through a decomposition based on the Wasserstein distance. In order to actually compute the Wasserstein distance matrix in this work we refer to the technique presented in [START_REF] Iollo | A lagrangian scheme for the solution of the optimal mass transfer problem[END_REF], where the optimal mass transfer problem is solved by a lagrangian method based on the fact that characteristics reduces to straight lines.

Numerical experiments

In this section we present four illustrations of increasing complexity ranging from an exact idealized case to experimental data.

Ideal vortex scattering

We consider two couples of counter-rotating ideal point vortices in the plane. The flow is incompressible and the vorticity is represented by four Dirac masses located at the vortex centers, so that the flow is irrotational almost everywhere. Under these hypotheses the flow is potential and the trajectories of the point vortices are obtained by the solution of an hamiltonian dynamical system. The detailed derivation of the governing equations is found in [START_REF] Meleshko | Vortex Dynamics and chaotic phenomena[END_REF].

The flow domain is Ê 2 and the coordinates of the vortex cores are:

x a = r 1 cos(θ 1 ) y a = r 1 sin(θ 1 ), x b = r 2 cos(θ 2 ) y b = r 2 sin(θ 2 ), x c = r 1 cos(θ 1 + π) y c = r 1 sin(θ 1 + π), x d = r 2 cos(θ 2 + π) y d = r 2 sin(θ 2 + π), (20)
where r 1 , r 2 , θ 1 , and θ 2 are the only variables necessary to describe the interaction. They are initialized as follows:

r 1 (0) = l 2 + (1 + β)f 2 1 2 , r 2 (0) = l 2 + (1 -β)f 2 1 2 , θ 1 (0) = arctan (1 + β)f l , θ 2 (0) = arctan (1 -β)f l , (21)
where l, β and f are three parameters that determine the distance of the dipoles, the relative distance of the counter rotating vortices and the offset of the dipole axis. The initial geometry of the system determines the nature of the scattering occurring between the dipoles. The ODEs describing the evolution are:

ṙ1 = - 2 sin(2(θ 1 -θ 2 ))r 1 r 2 2 π (r 4 1 -2 cos(2(θ 1 -θ 2 ))r 2 1 r 2 2 + r 4 2 ) , ṙ2 = - 2 sin(2(θ 1 -θ 2 ))r 2 r 2 1 π (r 4 1 -2 cos(2(θ 1 -θ 2 ))r 2 1 r 2 2 + r 4 2 ) , θ1 = 3r 4 1 -2 cos(2(θ 1 -θ 2 ))r 2 1 r 2 2 -r 4 2 2πr 2 1 (r 4 1 -2 cos(2(θ 1 -θ 2 ))r 2 1 r 2 2 + r 4 2 ) , θ2 = r 4 1 + 2 cos(2(θ 1 -θ 2 ))r 2 1 r 2 2 -3r 4 2 2πr 2 2 (r 4 1 -2 cos(2(θ 1 -θ 2 ))r 2 1 r 2 2 + r 4 2 ) . ( 22
)
The equations of motion are integrated via an adaptive-step fourth-order Runge-Kutta scheme in the time interval [0, 2.5]. In Fig. 2 We consider the L 2 norm of vorticity (enstrophy) as the transported density ̺. In this motion enstrophy is constant, so that we analyse the dynamics of four unitary Dirac masses. The Wasserstein distance between the time snapshots was computed by means of an exact combinatorial algorithm. For all the following cases 50 time frames were taken between initial and final time.

The embedding technique presented in the previous section was adopted. The spectrum of matrix B includes a few negative eigenvalues due to the fact that the distance is not euclidean. They are small in modulus so that B ≈ B + . In Fig. 3 the eigenvalues of the embedding matrix are represented for the first case described (a). Only two eigenvalues are relevant in the approximation of the phenomenon. The corresponding eigenvectors are represented in a phase-plane plot. The circles represent the components of the eigenvectors and can be associated to the time frames. Two directions that represent the optimal transports occurring before and after the interaction can be observed. The points which are not aligned represents the snapshots of the enstrophy configurations occurring during the interaction.

For case (b) (see Fig. 4) the spectrum of the embedding matrix is similar to that obtained for the first case: two eigenvalues emerge. The phase plot of the first two eigenvectors shows that the vortex interaction is quite different, but again two main transport direction corresponding to the dynamics before and after the interaction are present.

The third case (see Fig. 5) is different from the others. The interaction is weak and the resulting motion is practically an optimal transport. This third case may be considered as a perturbation of a single optimal transport. Indeed, in Fig. 5.a) the plot of the eigenvalues confirms that only one eigenvalue is important. The plot of the corresponding eigenvector in Fig. 5.b) show that most of the snapshots are aligned, that is, they may be obtained by non-linear interpolation (i.e., by transport) of the barycentric density via a unique optimal transport. In this simple example the comparison with standard POD is straight forward.

For the second and the third cases (Fig. 2 b) and c)) the autocorrelation matrix (i.e., the matrix of scalar products of the snapshots) is the identity matrix. Hence, no reduction can be provided by using POD. For the first case (Fig. 2.a)), the trajectories intersect, so that the autocorrelation matrix may not be diagonal for some specific samplings. However, only few extradiagonal elements appear, so that, even in this case, no signal reduction is possible.

One-dimensional Korteweg-de-Vries equation

The one-dimensional Korteveg-de-Vries equation is an interesting test case because transport and dispersion can be modulated. The model equation reads:

∂ t u + u∂ x u + γ∂ xxx u = 0, ( 23
)
where u(x, t) is the solution at time t ∈ [0, 2.5] and position x ∈ [0, 2π], γ = 1e -3 is the dispersion parameter. Periodic boundary conditions are imposed. The initial condition is u(x, 0) = u 0 (x) and for the present case is taken as: This equation was numerically integrated using a spectral discretization in space (256 Fourier modes), and a second order Cranck-Nicholson method with 10 3 steps in time.

u 0 (x) = 0.1 + exp -2(x - π 2 ) 2 . ( 24
)
The preserved measure is in this case the energy of the signal. Let e = u 2 , then it can be proved that ∀t,

E(t) = 2π 0 e(x, t) dx = E(0) = 2π 0 u 0 (x) 2 dx. (25)
The solutions for this model, shown in Fig. 6, show transport and dispersion. In particular, as γ is small, in the first part of the evolution the system evolves as an inviscid Burgers model (Fig. 6.b)). When a shock tends to form, the third derivative norm becomes large and dispersion makes the solution brake in several waves, with different characteristic velocities. For this case, a comparison between POD and the proposed technique is carried out. We take set of 24 snapshots of the solution energy, normalized by its integral (which is constant) between the time t = 0 and the time t = T = 2.5.

Let us compute the advection modes defined by means of Wasserstein distance. We computed the optimal transport between all the possible couples of snapshots and hence we formed the embedding matrix as described in the previous sections. The eigenvalues This advection mode may be defined using the straight line connecting the first snapshot with the last one (see Fig. 7.b)), as it is a good interpolant of all the snapshots. This is by definition the optimal mapping that transports the first snapshot, e(x, 0), into the last one e(x, T ). Let us denote by X 1 (ξ, t) this mapping.

The second step of the procedure consists in computing the residuals of this representation. The snapshots e(x, t i ) are mapped back to the initial configuration defined by X 1 (ξ, t). The differences between the solution snapshots pushed back and e(x, 0) = e 0 (ξ) is computed by

R i = R(ξ, t i ) = e(X 1 (ξ, t i ), t i )∂ ξ X 1 -e 0 (ξ), (26)
where e(X 1 (ξ, t i ), t i )∂ ξ X 1 are the snapshots e i (x, t i ) pushed back by the transport mode.

These residuals are then decomposed in the reference configuration by means of a standard POD analysis. Let us compare the results of this analysis with that of the POD applied Normalized eigenvalues of the autocorrelation matrix to the original set of snapshots. In Fig. 8 the normalized spectrum of the autocorrelation matrix is shown for the standard POD applied to the snapshots of the energy, and for the POD applied to the residuals pushed back. The eigenvalues cascade for the latter is steeper, i.e., a smaller number of modes is necessary for given representation error.

In Fig. 9.a) the transport mode U 1 is represented as function of X 1 (ξ, 0) = ξ. The lagrangian coordinate used to pushback the residuals is: X 1 (ξ, t) = ξ +tU 1 (ξ). In Fig. for POD and advection mode decomposition. In all three cases the advection mode decomposition shows an error that is significantly smaller with respect to that of POD. Also, it tends to diminish faster as the number of modes is increased. This could be anticipated from the spectrum decay. The reconstructions in the physical space for both methods are compared to the snapshot for which the representation given by advection mode decomposition is worst, that is n = 23.

When few modes are used, like in Fig. 12.a), POD is not able to reproduce the peaks that characterize the solution. In the cases shown in Fig. 12.b-c) POD provides a nonphysical solution due to its essentially oscillatory nature: the reconstructed kinetic energy is negative in a few regions of the domain. The advection mode decomposition in Fig. 12.b) is such that all the key features of the solution are well reproduced. In Fig. 12.c) the agreement is even improved.

Two-dimensional vortex shedding

In this section the vortex shedding past a confined circular cylinder is analyzed. We consider here the kinetic energy of the flow as the transported density. For each snapshot this quantity has been normalized in order to fulfill the constant mass constraint. In this case the normalization does not considerably affect the energy distribution since for this flow its integral is practically constant. Half a period of vortex shedding is considered.

The Reynolds number is Re = 200. The simulations are performed by a standard state of the art simulation code.

We took 10 snapshots of the normalized kinetic energy of the flow and computed the Wasserstein squared distance matrix. To do so, the space resolution adopted was 200 × 100, resulting in 2 • 10 4 collocation points. We used a multilevel scheme with 4 grids and linear interpolation kernels to compute the Wassestein distance. The core of this method is similar to what presented in [START_REF] Iollo | A lagrangian scheme for the solution of the optimal mass transfer problem[END_REF] In Fig. 13.a) the first 9 normalized eigenvalues (in absolute value) of the embedding matrix associated to the problem are shown. The cascade is less steep with respect to the one observed for the vortex scattering and the Korteweg-de-Vries equation. However, even in this case a small number of modes may be retained since there are two dominant eigenvalues, see Fig. 13.a). In Fig. 13.b) the phase plot of the first two eigenvectors is shown, revealing a remarkable structure: the points are located on a circle. This means that, given two orthogonal mappings φ 1 , φ 2 the flow is is the time corresponding to one period. This analysis suggests that three snapshots are sufficient to compute the barycentric density and two orthogonal base mappings. Three snapshots are taken at t = 0, t = 0.25, t = 0.5, equally distributed on half a period. In the following ̺ 0 is the kinetic energy distribution at the very beginning, ̺ 1 the kinetic energy at a quarter of period and ̺ 2 that of at half a period. The optimal transport is computed by means of the multilevel algorithm, then, the obtained mapping is used to transport ̺ 0 into the barycentric one, according to:

X G = ξ + 1 2 ∇ ξ φ 2 , (27)
where φ 2 is the optimal transport from ̺ 0 to ̺ 2 and the factor 1/2 means that the collocation points are moved by half the displacement that allows to map ̺ 0 into ̺ 2 .

Hence:

̺ 0 (ξ) = ̺ G (X G ) det(∇ ξ X G ). (28)
In Fig. 14 isocontours, 30 lines between the maximum and the minimum between the kinetic energy distributions ̺ 0 and ̺ 2 . It is not perfectly symmetrical with respect to the x axis, and this is due to the fact that the considered snapshots have a slight asymmetry too. Let us remark that the barycentric density is not a configuration happening in the physical evolution of the system and it is not the snapshot average. However, the average position of the structures and the characteristic distance of the vortex may be observed. In Fig. 15 shown corresponding to the kinetic energy at t = 1/8, which is the farthest from the snapshots used to build the reduced space. Some errors appear concerning the position and the intensity of the structures. However, these errors are localized in the wake. The structures near the body are well captured. This test case shows that a sufficiently good approximation of the flow patterns can be inferred by using little of information about the flow. The representation adopted is the minimal one, consisting in two advection modes computed using three snapshots.

As previously done for the Korteweg-de-Vries equation, the residual pushed back to the reference configuration (in this case the barycentric density) has been computed and the results are compared to classical POD applied to the snapshots sequence. In Fig. 17 to describe the periodic oscillation for all times.

Hurricane Dean

We consider the satellite images of the Caribbean see between the 17th and 22nd August 2007 showing the trajectory of hurricane Dean. The images are based on data from NOAA archive (http://www.aoc.noaa.gov/archive.htm) to get location and timeline of the storm. The images were generated by a NASA Goddard Space Center application (http://disc.sci.gsfc.nasa.gov/daac-bin/). Again, only 10 snapshots are considered. In this case, the density is simply defined as the normalized grey scale image. The time sampling was 6 hours, so that the time interval between the first snapshot and the last one is T = 54 hours. In Fig. 18 Here, we concentrate on the possibility of using advection mode decomposition and POD to estimate a subsequent image not included in the initial database. The most energetic POD mode obtained from the 10 snapshots database is shown in Fig. 19.a): it is the average of the snapshots. The average position of the hurricane may be inferred from As for advection mode decomposition, a set of N s -1 mapping is consider (N s = 10 is the number of snapshots). This allows to map the snapshot ̺ n into ̺ n+1 . This set may be used to define, by numerical integration, a lagrangian coordinate for the system. Let φ n be the mapping that transports ̺ n into ̺ n + 1. The following holds:

X(ξ, t n ) = tn 0 v dτ = tn 0 ∇φ dτ = Q(ξ, φ n ), ( 29
)
where X is the lagrangian coordinate and Q is a quadrature formula. In this case a second order Adams-Bashforth scheme was used to integrate the mappings. Let us remark that X is not the lagrangian coordinate of the physical system, but it is a mapping that allows to recover its state starting from the initial condition thanks to the Wasserstein map between each image. In Fig. 20 the lagrangian coordinate is represented at times t = 0, T /2, T . The deformation of the space is more intense around the hurricane (see for instance Fig. 20.b-c) ). This is due to the fact that its motion is more coherent than that of all the other structures. Optimal transport allows to highlight this feature in a straightforward manner.

We computed the transport modes and the first one is shown in Fig. ??. These modes represent potentials that suitably combined allow to map the first image of the sequence to all the others. With three transport modes, we can recover all the images with an error which is about 9% in norm L 2 norm of the gray scale.

We considered a subsequent hurricane image at time T * = 60h not included in the database to compute the POD modes and the advection modes. For POD representation, a simple problem of approximation is investigated, i.e., how accurately the new snapshot is represented using the POD modes computed using the database whose last snapshot is at T = 54h. The error of the reconstruction of this last image based on POD modes is about 16% using all the POD modes and does not vary significantly with the number of the modes used.

The approach based on Wasserstein distance allows to extrapolate the regularized lagrangian coordinate in a natural way. The discrete points (computed as explained above) X(ξ, t n ), n = 1, . . . , 10, are used to estimate X(ξ, t n+1 ) via a standard polynomial extrapolation. The error between the extrapolated lagrangian coordinate and that computed using of the novel datum at T = 60h differs by 0.0108 in norm L 2 . This means that the position of the hurricane is extrapolated with an error that is at most of 1%. In terms of gray scale image reconstruction, the L 2 error of the whole field is about 9% and it is due essentially to the small random structures that are not well captured. The extrapolated image via advection modes is still better reconstructed compared to standard POD, for which we used the best possible estimate, i.e., the actual image.

Conclusion

In this paper we have proposed an advection based modal decomposition. It exploits the Wassestein distance to define a transport mode hierarchy that describes the main features of advection. Additional features are included by using a POD mode expansion of the residual in a reference domain. When dealing with systems in which transport is the leading phenomenon, this method provides an efficient and meaningful low-order representation. This has been shown by contrasting the results of this approach to standard POD analysis for different cases ranging from idealized vortex motion, to actual hurricane data.

Appendix: Normalization of the embedding

Let us suppose that the densities ̺ i are taken by uniformly sampling in time an optimal transport between ̺ 0 and ̺ Ns-1 . In this case, X i , mapping ̺ 0 into ̺ i is written by interpolation (see [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF]):

X i = ξ + i∆t∇ ξ Φ(ξ), ( 30
)
where ∆t is the sampling time, and Φ is a (almost everywhere) convex potential. In this particular case, all the densities are aligned on a one-dimensional subspace of the Wasserstein space, since they belong to the same optimal transport. Hence, we expect that only one eigenvalue of matrix B is different from 0.

The squared Wassertstein distance between the ith and the jth sample is:

W 2 (̺ i , ̺ j ) = Ω i ̺ i (η)|X ij (η) -η| 2 dη, (31)
where X ij (η) is the optimal mapping between ̺ i and ̺ j . For these particular mappings the elements of the matrix have the form:

D ij = W 2 (̺ i , ̺ j ) = C N 2 s (i -j) 2 , ( 32
)
where C is a constant, representing the squared Wasserstein distance (i.e. twice the kinetic energy) of the unique mapping linking all the snapshots. The time at which the last snapshot is taken is taken to be T = 1.

Let us consider the matrix D = (ij) 2 , D ∈ Ê n×n and prove that the associated embedding matrix B has only one zero eigenvalue and that this value is λ = n(n+1)(n-1)

12

.

The elements of B are computed using two standard results in finite series:

By performing all the matrix vector products and exploiting the projector properties we have:

-

of B .

 . Let us take the positive part of B B + = US + U T where S + = S + |S| 2 (15) so that D ≈ B + , D -B + = S -|S| 2 and B + is of course positive semidefinite. It is now possible to determine matrix X whose columns are the position vectors χ as X = U √ S + . A normalization technique, which is useful when a large number of snapshots is taken is discussed in Appendix.

  scheme in the time interval [0, 2.5]. In Fig. 2, three different situations are represented: a) a scattering where the vortices keep their partner (the parameters used are: l = 1.5, β = 0.5, f = 0.25); b) a case where the vortices exchange their partner and escape with the counter rotating vortex belonging to the other dipole (l = 1.0, β = 0.75, f = 0.15); c) a weak interaction in which the dipoles simply move on (almost) straight lines (l = 2.0, β = 0.15, f = 0.30).

Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 2345 Figure 2: Three different scattering trajectories of vortex cores for: a) l = 1.5, β = 0.5, f = 0.25 b) l = 1.0, β = 0.75, f = 0.15 c) l = 2.0, β = 0.15, f = 0.30. 10

Figure 6 :

 6 Figure 6: Solution of the KdV model: a) contour of x-t diagram, 20 values between minimum and maximum b) solution at different times, where T=2.5.

Figure 7 :

 7 Figure 7: Embedding matrix analysis: a) first 10 normalized eigenvalues in modulus; b) components of the first eigenvector.

Figure 8 :

 8 Figure 8: Comparison between standard POD and the POD of the residuals pushed back.

  9.b) and c) the first and second residual mode is represented as function of ξ. Let us remark that these modes have significant variations in the portion of the domain where the initial energy distribution varies. Let us compare these modes with the classical POD modes for the energy distribution, represented in Fig. 10. The POD modes are quite different from the residual modes: they are global and their support approximately extends to the whole domain. We compare now the reconstruction of the snapshots using the advection mode decompoition (AMD) and POD. The same number of modes is used for both the techniques. Hence, for the POD reconstruction we will use one additional mode compared to the number of residual modes. This allows a fair comparison because for the advection mode decomposition we also have (in this case) one transport. Three different reconstructions are considered using 4, 6 and 10 POD modes respectively. In Fig. 11 the L 2 normalized errors are plotted as function of the snapshot number top (b) middle (c) bottom

Figure 9 :

 9 Figure 9: Advection modes expansion for KdV model: a) Transport mode b) First mode of the residuals c) Second mode of the residuals.

Figure 10 :Figure 11 :

 1011 Figure 10: POD modes for the normalized kinetic energy: a) First mode b) Second mode; c) Third mode.

Figure 12 :Figure 13 :

 1213 Figure 12: Reconstruction of the solution. Black line is the solution, blue line is the advection mode decomposition, red line is POD representation when: a) 4 modes, b) 6 modes, c) 10 modes are used 20

Figure 14 :

 14 Figure 14: Barycentric density for the kinetic energy of the flow around a circular cylinder:

Figure 15 :

 15 Figure 15: Isocontours of the base mappings: 30 lines between the maximum (1.25e -3) and the minimum -1.25e -3.

Figure 16 :Figure 17 :

 1617 Figure 16: Isocontours (30 lines between the maximum and the minimum) of the reconstruction (upper line) and the simulation (lower line) for a) t = 0 and for t = 1/8

  three images are shown, at time t = 0, T /2, T . The resolution is 512 × 256.

  Figure 18: Images of the hurricane Dean at a )t=0, b) t=27 hours and c) 54 hours

  Figure 19: POD modes: a) first mode b) second mode c) third mode 28

(a) top (b) middle (c) bottom Figure 20: Lagrangian Coordinate at time a) t=0 b) t=T/2 c) t=T

The relation written above states that B is the tensor product of a unique vector, whose components are y i = (ki). This is sufficient to prove the first point. Let us now explicitly compute the only non-zero eigenvalue λ. Again, the results on finite series are used, leading to:

The normalization condition for the generic B ∈ Ê Ns×Ns is then: