David Cristinacce
email: david.cristinacce@manchester.ac.uk

Tim Cootes

Facial Motion Analysis using Clustered Shortest Path Tree Registration

We describe a method of automatically annotating video sequences, defining a set of corresponding points in every frame. This is an important pre-processing step for many motion analysis systems. Rather than tracking feature points through the sequence, we treat the problem as one of 'groupwise registration', in which we seek to find the correspondence between every image and an automatically computed model reference, ignoring the ordering of frames. The main contribution of this work is to demonstrate a method of clustering the frames and constructing a shortest path tree over the clusters. This tree defines the order in which frames will be registered with an evolving estimate of the mean. This technique is shown to lead to a more accurate final result than if all frames are registered simultaneously. We describe the method in detail, and demonstrate its application to face sequences used in an experiment to assess the degree of facial motion. The resulting ranking is found to correlate well with that produced by human subjects.

Introduction

Our aim is to achieve accurate registration between all frames of a video sequence, by finding corresponding control points in all frames. The corresponding points then define a dense correspondence between frames, which can serve as a basis for motion analysis.

A critical step in the registration process is the initialisation of control points in each frame. A previous approach described by Cootes et al. [START_REF] Cootes | Groupwise construction of appearance models using piece-wise affine deformations[END_REF], initialises a grid of control points at the same location in all images. It then performs a series of translation searches and similarity transforms, followed by more complex warps, to register the whole set of images. Hence each image is treated equally and processed simultaneously.

In our approach we divide the video frames into smaller clusters of similar images, that are easier to register using the groupwise approach. The clustering method also provides an approximate starting point in each frame and the model generated from previous clusters is used to initialise the registration process in later clusters. A shortest path tree structure between the clusters enables an incremental approach to the registration which allows the initial processing to occur on easy data (in our case -frontal faces) and later matching to occur on outlier frames. We find this incremental clustered tree approach to give superior results to the single groupwise method, when tested on the manually annotated FGNET data set [2].

Our tree based registration is applied separately to 59 video sequences of different individuals. The corresponding points are used to build a multi-person statistical shape model [START_REF] Cootes | Active shape models -their training and application[END_REF]. The paths of each individual sequence are analysed in the combined shape space and used to determine the relative amount of facial motion for different individuals.

Background

There are many approaches to vision based motion analysis of human faces.

One method is to track facial features approximately and use machine learning techniques to make deductions about the underlying image, treating small tracking errors as noise. This method is adopted by Pantic et al. [START_REF] Valstar | Combined support vector machines and hidden markov models for modeling facial action temporal dynamics[END_REF] and successfully used to identify facial expressions. Similar machine learning approaches are adopted by face recognition systems [START_REF] Zhao | Face recognition: A literature survey[END_REF] and could be applied to the task of estimating facial motion.

An alternative approach is to compute a dense motion correspondence explicitly from the image data and apply further analysis to tracked features, independent of the original images. To do this successfully we require accurate localisation of corresponding features on each frame of a video sequence. To achieve this we use an extension of the groupwise registration approach due to Cootes et al. [START_REF] Cootes | Groupwise construction of appearance models using piece-wise affine deformations[END_REF], which is similar to the image encoding method described by Baker et al. [START_REF] Baker | Automatic construction of active appearance models as an image coding problem[END_REF]. An overview of general image registration methods is found in Zitova [START_REF] Zitová | Image registration methods: A survey[END_REF].

Image registration is computationally expensive compared to sequential tracking methods, e.g., Simultaneous Modelling and Tracking (SMAT) [START_REF] Dowson | Simultaneus modeling and tracking (smat) of feature sets[END_REF] or Constrained Local Models [START_REF] Cristinacce | Detection and tracking with constrained local models[END_REF], however there are several advantages from processing the whole set of frames together instead of in chronological order. With tracking methods there is the risk of "losing track" and the danger of "drift" accumulating over the sequence. Tracking methods use the previous frame for initialisation, however when considering the whole sequence another frame may be more appropriate. As our processing is offline, we process the whole set of frames simultaneously to obtain a more accurate feature correspondence between frames.

Accurate feature correspondence is essential for motion analysis, therefore the main contribution of this paper is a method of improving image registration across a video sequence. We find that a simple detection and clustering approach combined with a shortest path tree data structure improves the registration significantly.

Given our improved registration and a set of corresponding features, we build a statistical shape model [START_REF] Cootes | Active shape models -their training and application[END_REF] of variation over a set of facial motion videos. We then deduce the average distance travelled across each frame of the sequence to evaluate the magnitude of the motion of each face, which is shown to correlate with subjective scores provided by human subjects.

Methodology

Initialisation and Clustering

The aim of the initialisation step is to find a set of regions r i which encompass the face in each frame F i . Additionally a set of cluster templates {X j } are computed and each region r i is associated with a particular template X j . We therefore aim to simultaneously detect and cluster the face region across a set of video frames.

Our method requires an initial starting region r 0 in the first frame of the sequence F 0 . A template is generated from this region X 0 . This template is then used to search subsequent frames F i and the best matching region r i selected. If the best match of the template falls below a threshold T then a new template X j is created and added to the current list of templates. The current list of templates are ordered according to the their similarity to the current best match and the number of templates used to search any individual frame F i is limited to k max .

Formally the detection and clustering algorithm is as follows:-Algorithm 1 Initialisation and Clustering Algorithm Note in our experiments the matching score m j was computed using normalised correlation and a sliding window search around the previous frame region r i-1 . The normalised correlation threshold limit to create a new template is set to T = 0.90. The maximum number of templates applied to each frame was k max = 10. These settings were able to cluster the FGNET [2] talking face sequence (see Figure 2) into 53 separate clusters. The generated templates {X j } are shown in Figure 1 as nodes of the shortest path tree. The construction of the tree is described in Section 3.2.

Shortest Path Trees

The clustering method described in Section 3.1 produces a set of templates {X j }. We wish to construct a tree from this set of clusters which arranges similar clusters together and also minimises the distance from a given root template.

A suitable structure is the "shortest path" tree which minimises the total distance from the root to each node, over all nodes. Given a selected root node and distance measure between each pair of templates this can be computed using Dijkstra's algorithm [START_REF] Cormen | Introduction to Algorithms, Second Edition-Dijkstra's algorithm[END_REF].

The distance measure we use is the sum of squared pixel errors between the cluster templates, which is found to produce a shortest path tree which split the video sequence into associated templates. For example templates where the head is orientated in a particular direction are grouped together (see Figure 1).

Fig. 1. FGNET talking face template cluster tree diagram

The root template selected to compute the shortest path tree is typically the template associated with the largest number of frames. One frame associated with this template will be labelled manually. Then the registration algorithm will construct a correspondence from this frame to all other frames. The shortest path tree minimises the distance from the labelled frame and therefore reduces the chances of the registration failing.

Groupwise Image Registration

Given a set of unlabelled images, we wish to automatically find corresponding points in all of the images. We use the method described by Cootes et al. [START_REF] Cootes | Groupwise construction of appearance models using piece-wise affine deformations[END_REF] which optimises the total cost of encoding exact copies of the training set, the corresponding control points in each image being synonymous with a compact representation of the data.

Our groupwise approach uses a very simple mean image I m and mean shape S m , which has been shown to be effective in [START_REF] Cootes | Groupwise construction of appearance models using piece-wise affine deformations[END_REF]. The training set images {F i } are then approximated by placing control points {c i } in each image, which represent a mapping from the mean shape S m to each individual frame. The encoding cost of the training set C total is as follows:-

C total = C model + C params + C residuals (1)
Where C model is the cost of encoding the model, C params the cost of encoding the model and warp parameters and C residuals = ∑ i R i is the cost of encoding the residual pixel errors over all frames. In our case C model and C params can be ignored, as the changes in encoding cost of the mean image I m , mean shape S m and warp parameters are negligible compared to the cost of encoding the residuals R i for each frame.

Therefore the only variable to optimise is ∑ i R i , which has two components for each image i.e. R i = S i + T i . Where S i is the cost of storing the residual shape variation between the mean shape S m warped into each frame and the exact control points c i . T i is the cost of encoding the texture error over all pixels for each frame F i given the warp defined by c i . Since both S i and T i depend on the control points c i , the principal task of registration is to optimise the control point location c i in each frame of the video.

The groupwise algorithm proceeds by computing the mean image I m and mean shape S m from a set of starting control points {c i } and frames {F i }. The control points c i can then be independently optimised for each image F i , to lower the total encoding cost C total . When the control points c i have been updated in each image, the mean image I m and shape S m are updated. We refer to the process of computing the mean model and adjusting the control points across all images as a registration "stage".

Typically a full registration consists of several such stages, which perform a image pyramid based coarse-to-fine registration of the whole data set with increasingly complex warps of the control points c i in each image, eventually resulting in the individual control point co-ordinates being optimised at the maximum image resolution.

Incremental Tree Registration

The groupwise registration method described above can be applied incrementally given a clustered tree of templates (see Figure 1).

Each node of the tree actually represents a set of individual frames. Initially the groupwise registration algorithm can be applied to frames associated with the root template in isolation. Then for each node of the tree the grid of control points can be initialised from the closest example in the parent cluster. Also the region r i associated with frame F i is used to warp the grid points from the closest frame. For each node, 16 registration "stages" are applied.

Each time a node registration is completed, the running mean texture and shape are updated, which maintains efficiency as each cluster node can then be processed individually, with the current mean and groupwise registration. The computational complexity of the incremental tree registration is the same as the single groupwise method in the sense that the number of times the control points c i are optimised in each image is the same, whether the stages are applied to all images together or to separate image clusters.

The rational behind the incremental approach is that similar frames will be processed together at the start, which is a relatively simple registration task. Later more difficult clusters will be registered and the mean updated, whilst the registration points in previous clusters remain unchanged.

If all frames are registered together in one single registration process (as described in Section 3.3) there is some risk that the outlier frames will distort the mean shape and texture. The incremental method is designed to reduce the affect of outliers by allowing the registration to concentrate on the largest template clusters in isolation, then attempt to fit to the more difficult images.

Experiments

Incremental Tree vs Single Groupwise Registration

The FGNET talking face data set [2] is used to test the performance of the tree based registration (see Section 3.4) and the simple groupwise registration (see Section 3.3). Example frames from this publicly available data set are shown in Figure 2. The original FGNET data consists of 5000 frames of an individual talking in front of a camera. Each frame has been manually landmarked with 68 feature points. We select every 10 th frame to create a smaller data set of 500 images.

We applied the clustering algorithm described in Section 3.1 to deduce a set of regions r i for each of the 500 frames F i . On the first frame F 0 a regular grid of control points is generated over the region r 0 .

We then applied the following two methods:-1. Initialise the control grid in each frame using regions {r i } computed by the cluster method and run the single groupwise registration (see Section 3.3), for 56 stages. 2. Use the shortest path tree and regions {r i } to run the incremental groupwise registration (see Section 3.4), using 16 stages for each of the 53 nodes of the tree.

When applying single groupwise registration. The evolution of the mean image is shown in Figure 3. Figure 3 shows the mean image becoming sharper as the registration algorithm progresses. The mean image is initially blurred (see Figure 3(a)) but as the grid points in each image move to equivalent locations in each frame the image becomes sharper (see Figure 3(e))

When applying the incremental tree method, the progression of the mean images is different, see Figure 4. shows the mean image after completing the groupwise registration on the root template cluster, which consists of 16 stages on the 30 images in this cluster. The mean image is quite sharp, which reflects the fact that all the images in the first cluster are similar. The later mean images, after completing the registration on the later clusters show similar sharpness, but the mean texture and mean shape are updated as more clusters are added. The final mean image of the incremental registration method (Figure 4(e)) looks similar to the single groupwise mean image (see Figure 3(e)), because at this point all the data has been registered.

Figure 5(a) shows the accuracy of the registration methods relative to the ground truth, for both the incremental and single groupwise approaches. The registration accuracy (mean point error on y axis) is determined at each stage using the following method:-1. Warp all the manually labelled ground truth points from each image into the reference frame (using the computed control points c i in each image) 2. Compute the mean x and y co-ordinates for each point in the reference frame 3. Warp the mean points in the reference frame back to the original images 4. Compute the average point to point error over all points and images This determines the consistency of the registration relative to the ground truth. The average point to point error is normalised by the ground truth interocular eye separation, to give an accuracy measure which is independent of image size. The dashed line on Figure 5(a) shows that when applying the single groupwise registration to the FGNET data the point to point error gradually decreases over time. Therefore with each iteration the registration improves, until the error reduces to 0.0275. Which represents an average point to point distance error of 2.75% (standard error 0.21%) of the eye separation, i.e. a reasonably accurate registration.

When applying the incremental tree registration (see solid line Figure 5(a)) the pattern of point to point errors over time is different. For each individual cluster (i.e. every 16 stages) the point to point error first rises as the new cluster is added, due to the new images being only approximately registered. Within each cluster the error decreases as the registration is applied to the new images. However as more clusters are added the registration task becomes more difficult, which steadily increases the registration error over time, resulting in the jagged pattern shown by the solid line Figure 5(a).

However the final incremental registration error on the whole data set is 2.43% (standard error 0.05%) of the eye pupil separation, which is a 12% improvement over the single groupwise approach. The small standard errors indicate that the incremental method achieves a statistically significant better result for the same computational effort compared to the single groupwise approach.

The cumulative error of the final registration for both methods is plotted in Figure 5(b), which shows that 90% of the images marked up using the incremental tree method (solid line) have < 4% mean point error relative to the ground truth, which is greater than the 86% success rate achieved by the single groupwise approach (dashed line) at the 4% threshold.

Example automatically found feature points from the incremental tree registration are shown in Figure 6. This process simulates the user manually landmarking the mean image then generating the labelled points automatically from the registration. The whole data set (500 images) can be labelled accurately by only manually marking one image (see Figure 6).

Facial Motion Experiment

This experiment was designed to assess the ability of human volunteers to recognise celebrities under difficult conditions (the celebrities being famous people in the country where the experiment was performed, i.e. the UK). The videos are poor quality images collected from television broadcasts. See Figure 7 for two examples from the 59 sequences used.

The psychological experiment tested the hypothesis that human face recognition is more successful when observing celebrities who show a large amount of facial movement. As part of this experiment the volunteer subjects rated the amount of facial motion present in each video (giving a subjective score from 0-9). To corroborate the manual motion ratings, we manually marked up one frame of each sequence and using incremental tree registration (Section 3) to automatically mark the whole set of frames for each video. We then compared the amount of shape variation for each individual with the human rankings.

Figure 8 shows the shortest path trees computed from the videos, by taking the most common template as the root node, for two individuals. In both cases the shortest path tree divides the templates extracted from the 64 frame video sequence in a visually sensible form. Results of automatic annotation for the two sequences are shown in Figure 9, which show that the incremental registration was successful. The method is able to cope with the opening and closing of the subject's mouth and head variation as the celebrity speaks, approximately facing the camera in each video. Given automatically labelled points for each of the 59 sequences. A shape model of the whole set of (64*59=3776) images was built and the euclidean distance travelled in shape space calculated for each of the 59 individuals. This automatic estimate of facial motion was then used to rank the individuals and the ranks compared with the human perceived motion ranks. The top ten ranks are shown in Figure 10. In Figure 10, 7 out of the top 10 individuals are the same. Indicating a strong correlation between the automatic computer estimate of facial motion and the motion perceived by humans. A spearman rank correlation test confirms this with a correlation value of 0.68 which is significant at the p < 0.01 level for a sample size of 59, against the null hypothesis of random permutations.

Conclusions

This paper has described an incremental registration method which is shown to outperform a single groupwise registration on a publicly available data set [2]. The method relies on a simple clustering algorithm and the building of a shortest path tree to cluster the data, see Figure 1.

The tree structure allows the incremental registration of the largest cluster of data first and leaves the registration of outlier data until later. The resulting registration is shown to be 12% more accurate than a simple groupwise approach. Due to the incremental computation of the mean image and shape, the tree registration method has the same complexity as the groupwise method.

The incremental method was successfully applied to low quality videos collected from television broadcasts and the resulting motion analysis was able to achieve a high correlation with human facial motion perception. The approach is general and could easily be applied to other types of video motion analysis.

Fig. 2 .

 2 Fig. 2. Example of frames from the FGNET talking face sequence

Fig. 3 .

 3 Fig. 3. Single groupwise mean texture on FGNET talking face data, after each stage of the registration

Fig. 4 .

 4 Fig. 4. Incremental tree registration mean texture on FGNET talking face data, after registration of each tree node

Fig. 5 .

 5 Fig. 5. Single Groupwise vs Incremental Tree Registration on FGNET data

Fig. 6 .

 6 Fig. 6. Example of automatic labelling on FG-NET talking face sequence

Fig. 7 .

 7 Fig. 7. Example frames from two video sequences of 'celebrities'

Fig. 8 .

 8 Fig. 8. Example shortest path trees built from 2 of the 59 celebrity videos sequences

Fig. 9 .

 9 Fig. 9. Example of automatic labelling for two video sequences of 'celebrities'

Fig. 10 .

 10 Fig. 10. First 10 individual faces ranked by human perceived motion ratings (column 2) and automatic facial motion estimate using distance travelled in the shape space (column 3)


Lior Wolf

Tal Hassner

Yaniv Taigman

Descriptor Based Methods in the Wild

Recent methods for learning similarity between images have presented impressive results in the problem of pair matching (same/notsame classification) of face images. In this paper we explore how well this performance carries over to the related task of multi-option face identification, specifically on the Labeled Faces in the Wild (LFW) image set.

In addition, we seek to compare the performance of similarity learning methods to descriptor based methods. We present the following results:

(1) Descriptor-Based approaches that efficiently encode the appearance of each face image as a vector outperform the leading similarity based method in the task of multi-option face identification. (2) Straightforward use of Euclidean distance on the descriptor vectors performs somewhat worse than the similarity learning methods on the task of pair matching.

(3) Adding a learning stage, the performance of descriptor based methods matches and exceeds that of similarity methods on the pair matching task. (4) A novel patch based descriptor we propose is able to improve the performance of the successful Local Binary Pattern (LBP) descriptor in both multi-option identification and same/not-same classification.

Introduction

The Labeled Faces in the Wild (LFW) database [START_REF] Huang | Labeled faces in the wild: A database for studying face recognition in unconstrained environments[END_REF] offers a unique collection of annotated faces captured from news articles on the web. The dataset is published with a specific benchmark, which focuses on the face recognition task of pair matching. In this task, given two face images, the goal is to decide whether the two pictures are of the same individual. This is a binary classification problem, in which the two possible outcomes are "same" or "not-same".

The simple binary structure of Same-Not-Same classification simplifies the design of benchmark experiments. However, in many face recognition applications the task is quite different, and can be defined as follows: given a gallery containing labeled face images of several individuals (one or more face images for each person), classify a new set of probe images. The classification label can be "unknown person" or that of one of the individuals in the gallery.

The existence of the "unknown person" label is crucial at the application level. However, it makes benchmarking and experimental design much more challenging. First, the design should address the distribution of unknown faces. This can be done for example by providing a training set containing images of individuals not present in the gallery. Second, the likelihood of encountering an 1 unknown individual should be defined and may greatly affect the results. For these reasons, in this work, as in many other face recognition reports, a more limited task is studied: the multi-option identification task, where the "unknown person" is excluded.

Methods developed for each of the three identification tasks discussed above can be adopted for any of the other tasks. For example, one can compute scores between all gallery and all probe images using pair matching techniques and then use them for identification based on a winner takes all approach.

Since the LFW benchmark focuses on the pair matching problem, it is important to understand whether the reported success of algorithms on this benchmark carries over to other identification tasks. In addition, we test how descriptor based techniques, which are often studied in multi-option recognition tests, perform in the pair matching task.

Finally, we develop novel image descriptors that are able to improve performance both in the multi-option task and in the pair matching task. These descriptors are based on patch statistics, and we suggest using them in combination with other features.

Existing methods

In the previous section we have focused on the importance of the specific face recognition task. We now turn our attention to relevant face recognition algorithms and techniques.

Modern image similarity learning techniques Recently, some effort has been devoted to the estimation of visual similarities between two unseen images [START_REF] Jain | Discriminative training of hyper-feature models for object identification[END_REF][START_REF] Ferencz | Building a classification cascade for visual identification from one example[END_REF]4,[START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF], and such methods have been applied to determine whether two images belong to the same person. One method [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF] that has shown good results for uncontrolled imaging conditions uses Randomized Decision Trees [START_REF] Geurts | Extremely randomized trees[END_REF] and Support Vector Machines. In the first image of the pair, image patches (fragments of the image) are selected at random locations. For each patch the most similar patch in the second image is searched at a nearby image location. A decision tree is trained to distinguish between pairs arising from matching images and those arising from non-matching images. Given a pair of unseen images, a Support Vector Machine (SVM) classifier is used to determine if they match by aggregating the Decision Tree output of many image patches. More specifically, the leaf (terminal classification node) of the decision tree is recorded for each pair of patches, and the SVM input consists of a binary vector that indicates if patches arising from a particular pair of images have reached a certain leaf.

Descriptor based methods for face recognition Face Images can be most readily described by statistics derived from their intensities. Intensities have thus served in many template-based methods. The intensities were sometimes normalized and sometimes replaced by edge responses [START_REF] Brunelli | Face recognition: Features versus templates[END_REF]. More recently [START_REF] Wiskott | Face recognition by elastic bunch graph matching[END_REF][START_REF] Meyers | Using biologically inspired features for face processing[END_REF][START_REF] Tan | Fusing gabor and lbp feature sets for kernel-based face recognition[END_REF], Gabor wavelets have been used to describe the image appearance. A texture descriptor called Local Binary Patterns (LBP) [START_REF] Ojala | A comparative-study of texture measures with classification based on feature distributions[END_REF][START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF][START_REF] Ojala | A generalized local binary pattern operator for multiresolution gray scale and rotation invariant texture classification[END_REF]has been shown to be extremely effective for face recognition [14]. The most simple form of LBP is created at a particular pixel location by threshholding the 3×3 neighborhood surrounding the pixel with the central pixel's intensity value, and treating the subsequent pattern of 8 bits as a binary number (Fig. 1). A histogram of these binary numbers in a predefined region is then used to encode the appearance of that region. Typically, a distinction is made between uniform binary patterns, which are those binary patterns that have at most 2 transition from 0 to 1, and the rest of the patterns. For example, 1000111 is a uniform binary pattern while 1001010 is not. The frequency of all uniform LBPs is estimated, while all non-uniform LBPs, which are typically around 10% of the patterns in an image, are treated as equivalent and given only one histogram bin. The LBP representation of a given face image is generated by dividing the image into a grid of windows and computing histograms of the LBP values within each window. The concatenation of all these histograms constitutes the image's signature.

In this work we propose a patch-based descriptor that has some similarities to a variant of LBP called Center-Symmetric LBP (CSLBP) [START_REF] Heikkilä | Description of interest regions with center-symmetric local binary patterns[END_REF]. In CSLBP, eight intensities around a central point are measured. These intensities are spread evenly at a circle every 45 degrees starting at 12 o'clock. The binary vector encoding the local appearence at the central point, consists of four bits which contain the comparison of intensities to intensities on the symmetric position (180 degrees/ 6 hours difference).

Multi-block LBP [START_REF] Zhang | Face detection based on multiblock lbp representation[END_REF] is an LBP variant that replaces intensity values in the computation of LBP with the mean intensity value of image blocks. Despite the similarity in terms, this method is very much different from our own. Multiblock LBP is shown to be effective for face detection, and in our initial set of experiments does not perform well for face recognition.

Patch-based approaches in recognition As mentioned above, a patch based approach [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF] provides state of the art capabilities in similarity learning of faces and of general images. Other successful object recognition systems based on patches include the hierarchical system of [START_REF] Ullman | Object classification using a fragment-based representation[END_REF].

The ability to detect local texture properties by examining the cross correlation between a central patch and nearby patches on both sides have been demonstrated in the texture segmentation system of [START_REF] Wolf | Patch-based texture edges and segmentation[END_REF]. In [START_REF] Shechtman | Matching local self-similarities across images and videos[END_REF] a central patch was compared to surrounding patches to create a descriptor which is an extension of the shape-context [START_REF] Belongie | Shape context: A new descriptor for shape matching and object recognition[END_REF] descriptor. The resulting descriptor has been shown to be highly invariant to image style and local appearance.

Improving descriptors by learning A large body of literature exists on the proper way of learning classifiers and distances for face recognition. In this work, we build upon basic classification methods, such as Linear Discriminant Analysis and One-Vs-All Support Vector Machine [START_REF] Allwein | Reducing multiclass to binary: a unifying approach for margin classifiers[END_REF], which is the simplest way to construct multi-class classifiers out of binary SVM.

Novel patch based LBPs

The LBP descriptor and its variants use short binary strings to encode properties of the local micro-texture around each pixel. CSLBP [START_REF] Heikkilä | Description of interest regions with center-symmetric local binary patterns[END_REF], for example, encodes in each pixel the gradient signs at the pixel in four different angles. Here we propose a family of related descriptors each designed to encode additional types of local texture information. The design of these descriptors is inspired by the Self-Similarity descriptor of [START_REF] Shechtman | Matching local self-similarities across images and videos[END_REF]. Specifically, we explore different ways of using bit strings to encode the similarities between neighboring patches of pixels, possibly capturing information which is complementary to that of pixel-based descriptors. Thus, employing patch based descriptors and pixel based ones in concert may improve the over-all accuracy of a classification system.

Three-Patch LBP Codes

As its name implies, the Three-Patch LBP (TPLBP) code is produced by comparing the values of three patches to produce a single bit value in the code assigned to each pixel. For each pixel in the image, we consider a w × w patch centered on the pixel, and S additional patches distributed uniformly in a ring of radius r around it (Fig. 2). For a parameter α, we take pairs of patches, αpatches apart along the circle, and compare their values with those of the central patch. The value of a single bit is set according to which of the two patches is more similar to the central patch. The resulting code has S bits per pixel. Specifically, we produce the Three-Patch LBP by applying the following formula to each pixel:

TPLBP r,S,w,α (p) = S i f (d(C i , C p ) -d(C i+α mod S , C p ))2 i
(1) Where C i and C i+α mod S are two patches along the ring and C p is the central patch. The function d(•, •) is any distance function between two patches (e.g., L 2 norm of their gray level differences) and f is defined as:

C 3 C 1 C P C 0 C 4 C 6 C 2 C 5 C 7 r α w w TPLBPr,8,3,2(p) = f (d(C0, Cp) -d(C2, Cp))2 0 + f (d(C1, Cp) -d(C3, Cp))2 1 + f (d(C2, Cp) -d(C4, Cp))2 2 + f (d(C3, Cp) -d(C5, Cp))2 3 + f (d(C4, Cp) -d(C6, Cp))2 4 + f (d(C5, Cp) -d(C7, Cp))2 5 + f (d(C6, Cp) -d(C0, Cp))2 6 + f (d(C7, Cp) -d(C1, Cp))2 7 (a) (b) (c)
f (x) = 1 if x ≥ τ 0 if x < τ (2)
We use a value τ slightly larger than zero (e.g., τ = 0.01) to provide some stability in uniform regions, similarly to [START_REF] Heikkilä | Description of interest regions with center-symmetric local binary patterns[END_REF]. In practice, we use nearest neighbor sampling to obtain the patches instead of interpolating their values, as this speeds up processing with little or no effect on performance. Once encoded, an image's signature is produced similarly to that of the CSLBP descriptor [START_REF] Heikkilä | Description of interest regions with center-symmetric local binary patterns[END_REF]. The image is divided into a grid of none-overlapping regions and a histogram measuring the frequency of each binary code is computed for each such region. Each of these histograms are normalized to unit length, their values truncated at 0.2, and then once again normalized to unit length. An image is represented by these histograms concatenated to a single vector.

Four-Patch LBP Codes

For every pixel in the image, we look at two rings of radii r 1 and r 2 centered on the pixel, and S patches of size w × w spread out evenly on each ring (Fig. 3). To produce the Four-Patch LBP (FPLBP) codes we compare two center symmetric patches in the inner ring with two center symmetric patches in the outer ring positioned α patches away along the circle (say, clockwise). One bit in each pixel's code is set according to which of the two pairs being compared is more similar. Thus, for S patches along each circle we have S/2 center symmetric pairs which is the length of the binary codes produced. The formal definition of the FPLBP code is as follows:

(p) = f (d(C10, C21) -d(C14, C25))2 0 + f (d(C11, C22) -d(C15, C26))2 1 + f (d(C12, C23) -d(C16, C27))2 2 + f (d(C13, C24) -d(C17, C28))2 3 (a) (b) (c)
FPLBP r1,r2,S,w,α (p) = S/2 i f ( d(C 1i , C 2,i+α mod S )- d(C 1,i+S/2 , C 2,i+S/2+α mod S ))2 i (3)
The final image signature is produced by using the same two-step normalization procedure described in Section 3.1.

Face same-not-same classification

The Labeled Faces in the Wild (LFW) dataset has two versions: the original version and the funneled version, in which images are automatically aligned using the method of [START_REF] Huang | Unsupervised joint alignment of complex images[END_REF]. In all of our experiments we use the funneled version only. We plan to add the results on the original images in the future. There are also two proposed pair matching benchmarks. We report results on the benchmark protocol called "image restricted training", for which public results are available for the algorithm of [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF] on the LFW web-site (http://vis-www.cs.umass.edu/ lfw/results.html).

The image restricted pair matching benchmark is a challenging one. It consists of 6000 pairs, half of matching subjects and half not, which are divided into 10 equally sized sets. The benchmark experiment is repeated 10 times. In each repetition one set is used for testing and nine others are used for training. The goal of the tested method is to predict which of the testing pairs are matching, using only the training data (the decision is done one pair at a time, without using information from the other testing pairs).

We test the performance of descriptor based methods on this benchmark, and focus on two questions: (1) How well do these methods perform compared to the average recognition rate of 0.7333 ± 0.006 obtained by the binaries of [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF]; (2) How can learning be applied to descriptor based methods in the pair matching (same-not-same binary classification) setting.

Distance thresholding for pair matching

The most straightforward approach for pair matching using image descriptors is to consider the distance between the vectors which encode the appearance of the two images. i.e., given two face images I 1 and I 2 which are encoded using some image descriptor g as g(I 1 ), g(I 2 ), the pairs are considered to match if d(g(I 1 ), g(I 2 )) < T , where d is a distance function and T is threshold.

The distance d can vary. We use the Euclidean distance and the Euclidean distance applied to the square roots of the values in g(I 1 ). The motivation for the second distance is that our descriptor vectors consist mostly of histograms, and applying square root prior to the distance calculations corresponds to the Hellinger distance between probabilities [START_REF] Pollard | A user's guide to measure theoretic probability[END_REF].

In order to learn the threshold one can choose the threshold that gives the highest recognition score on the 5400 examples of the training set. A method that gives similar performance, and can be generalized to more than one distance score per pair is to employ a binary Linear SVM. In the single distance case, we train a SVM classifier on the 5, 400 one-dimensional vectors each containing the distance between the two images of a pair. Then, this classifier is used to predict whether the 600 test pairs are matching or not using similar 1D input vectors. This experiment is repeated for the 10 train/test splits, and we record mean recognition rate as well as the standard deviation of it.

A simple way to combine multiple image descriptors and multiple distances is to create a vector of distances and run Linear SVM on this vector. Indeed, we find that combining multiple distances together improves results considerably. In Table 1 below, we report the recognition rate for the two distances (Euclidean, Hellinger) and for each of the four descriptors: LBP, Gabor (C1), Three Patch LBP (TPLBP) and Four Patch LBP (FPLBP). The Gabor descriptor is the C1 descriptor of [START_REF] Riesenhuber | Hierarchical models of object recognition in cortex[END_REF] used for face recognition in [START_REF] Meyers | Using biologically inspired features for face processing[END_REF]. The parameters of the other descriptors are given in the appendix. Note that other descriptors such as CSLBP [START_REF] Heikkilä | Description of interest regions with center-symmetric local binary patterns[END_REF] and direct use of gray values did not produce good results and are omitted. Table 1 also contains the combined result obtained by classifying the vector of 8 distances at once using a SVM classifier.

As can be seen, a direct application of the distance function perform somewhat worse than the 0.7333 ± 0.0060 recognition rate achieved by the much slower random-tree based similarity prediction [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF]. The LBP based descriptions perform similarly to one another, and the choice of distance does not change the results significantly. The combination of all descriptors and distances improves results beyond the best combination of a descriptor and a distance function.

Model learning for pair matching

In the distance based experiments, the training data was used in an improvised manner to classify very small vectors after the distances were taken. How can we use the training data to learn a matching score between pairs of descriptor vectors? There is a vast literature on supervised learning of similarity or of distance functions, e.g., [START_REF] Xing | Distance metric learning with application to clustering with side-information[END_REF][START_REF] Bar-Hillel | Learning distance functions using equivalence relations[END_REF][START_REF] Weinberger | Distance metric learning for large margin nearest neighbor classification[END_REF]. We tried several of these methods using the DistLearnKit Matlab toolbox (http://www.cse.msu.edu/ ∼ yangliu1/ distlearn.htm) and were unable to improve the classification accuracy.

An alternative learning framework is to consider the pair matching problem as related to the problem of learning from one example ("one-shot learning"; see [START_REF] Huang | Labeled faces in the wild: A database for studying face recognition in unconstrained environments[END_REF] for a comparison of this problem and the unseen pair matching problem). We learn a model of the person in the first image of the pair and try to classify the second image. We then replace the roles of the images and repeat. The average prediction of the two classifications is taken as the matching score.

Our procedure required two training sets: a set A containing face images that is used for the negative examples during the one-shot learning, and a second set B of matching and non-matching pairs, which is used to learn the decision threshold as before. Given a pair of images, we vectorize them using an image descriptor. We then train two classifiers by using either one of the two images as the positive example, and the set A as the negative examples. Afterwards, we apply each classifier on the other image of the pair and obtain a score. These two scores are averaged to obtain a similarity score.

This process is repeated to all training image pairs of set B, and a threshold on the obtained scores is learned using a SVM. Given a new test pair, we build a one-shot classifier to each of the images using the same procedure and average the two classification scores. We then apply the SVM threshold on the average to obtain a prediction.

Since two classifiers need to be trained per training pair in set B, and since the number of "negative" images (set A) may be large, using SVM to learn the underlying one-shot classifiers may be computationally demanding if proper care is not taken. Instead, we suggest using a Linear Discriminant Analysis (LDA) classifier [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. Note that due to the special structure of the problem, the LDA computation can be performed very efficiently. The within-class covariance matrix is constant, depends only on the set A, and can be inverted once. Moreover, the direction of the LDA projection can be obtained by extracting the leading eigenvector of a 4 × 4 matrix regardless of the dimensionality of the image descriptors.

For the LFW image matching benchmark, we use in each repetition, out of the nine training sets, one to construct the set A and eight for B. The single negative image partition contains a total of 1,200 images. Note that none of the subjects in the 1,200 images appears in the test set, since the LFW benchmark is constructed to have the persons in the training and test splits mutually exclusive [START_REF] Huang | Labeled faces in the wild: A database for studying face recognition in unconstrained environments[END_REF].

The results of the experiments are described in Table 2 Below. The descriptors used are the same as in the direct distance learning experiments. Here again we use either the original descriptor vectors, or their square root. In the latter case, instead of using the vector g(I) we use g(I).

The 8 descriptor/mode scores in the table are obtained by training SVM on 4, 800 (8 sets) 1D vectors containing the average of the two: the LDA projection of the second image obtained using one-shot model learning of the first image and the LDA projection of the first image obtained from the second model. The "Combined" classification is based on learning and classifying the 8D vectors which are the concatenations of the eight 1D vectors.

Table 2. Mean (± standard error) recognition rates on the funneled pair matching benchmark of LFW (Image-Restricted Training, "view 1") using the per image model learning method. Each row corresponds to one image descriptor, and the columns represent the use of the original descriptor or of its square root. The last row ("Combined") corresponds to training a SVM on the vector containing the 8 prediction scores. The results of LBP alone, using the model learning framework are similar to the results obtained by the method of [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF], which are the best results on the benchmark known to us. The results of combining all 8 scores outperform previous results significantly.

Image

Hybrid method

A last experiment was done by combining the direct distance method and the one-shot model-learning method. This was done by concatenating the 8 distances Fig. 4. ROC curves averaged over 10 folds of View 2 of the LFW data set. Each point on the curve represents the average over the 10 folds of (false positive rate, true positive rate) for a fixed threshold. The propose hybrid method is compare to the method of [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF].

of the first method with the 8 scores obtained in the second method. The second method provides predictions for only 4800 of the training pairs, and the SVM classifier was trained on 4800 training 16-D vectors . A recognition rate of 0.7847 ± 0.0051 is obtained, which is significantly higher than each of the other methods, indicating that the distance base method and the model learning method employ different aspects of the data. The ROC plot of the hybrid method is depicted in Figure 4.

Face identification

Next, we evaluate the performance of the descriptor based methods and the similarity approach of Nowak and Jurie [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF] on the task of image classification. To this end we use the LFW dataset, choosing only those subjects having enough images to contribute for both "probe" and "gallery" sets. In our experiments we use two images per person as probes and two as gallery. Thus, we employ a subset of the LFW image set which consists of the 610 subjects having at least four images. This subset contains a total of 6733 images.

The performance of the various methods as a function of the number of subjects (N ) was compared. We perform 20 repetitions per experiment. In each, we select N random subjects and choose two random gallery images and a disjoint set of two random probes from each. Note that both N and the number of repetitions were limited by the computational demands of [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF].

In order to learn the similarity function of [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF], matching and non-matching pairs are required. We are not able to use the training pairs given with the LFW benchmark (either "view 1" or "view 2") since the images in those pairs and the images in our train/test splits overlap. Therefore, in order to create independent training and testing sets we considered all subjects in the LFW database for which there are no more than three images per subject. There are 5139 such subjects (a total of 6500 images). From these images we extract the maximal number of matching pairs -1652 pairs, and 2000 randomly selected non matching pairs. This training set is somewhat smaller than the data set of 5400 pairs (half matching, half non-matching) used in the same/not-same experiments. Once a similarity function is learned, given a probe image we compute its similarity to each of the gallery images. The probe image is then classified by a winner-take-all approach.

We test the following descriptors: LBP [START_REF] Ojala | A comparative-study of texture measures with classification based on feature distributions[END_REF][START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF][START_REF] Ojala | A generalized local binary pattern operator for multiresolution gray scale and rotation invariant texture classification[END_REF], CSLBP [START_REF] Heikkilä | Description of interest regions with center-symmetric local binary patterns[END_REF], C1 Gabor descriptors [START_REF] Riesenhuber | Hierarchical models of object recognition in cortex[END_REF][START_REF] Meyers | Using biologically inspired features for face processing[END_REF], Image intensities, TPLBP, and FPLBP. In some experiments, descriptors are combined by concatenating several descriptors into a single vector. Classification in the description based tests is performed by training a One-Vs-All multiclass Linear SVM on the gallery images and using it to classify the probe images one by one.

Fig. 5 (a) presents the performance of the three LBP variants (LBP, TPLBP, FPLBP), as well as the performance of the combination of the three. Also shown is the performance of the Gabor (C1) descriptor, and that of the method of [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF]. To remove clutter, we omit the CS-LBP and the direct image intensities descriptors since they consistently under-perform the other descriptors. Adding Gabor to the LBP descriptors did not improve performance and those results are also omitted.

As can be seen, the similarity method of [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF] outperforms the descriptor based methods and their combinations. This comes at a price, however, as the testing stage of the similarity based approach is more than an order of a magnitude slower than that required for descriptor based classification. Moreover, the similarity based approach uses a large data set for training, which was not utilized by the descriptor based methods. A natural question now arising is how to modify the learning stage of the descriptor based methods in order to make use of this training set?

The learning algorithm we use is One-Vs-All SVM, in which for each subject a SVM is trained using the subject's gallery images as positive examples and all other gallery images as negative examples. We modify this method such that an extra set of negative training images is used for each binary classifier learned. These additional negative training images are randomly selected form the set of images used to train the similarity based method (see above). Hence the descriptor based classifiers use no more training information than it is used for the method of [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF].

Figure 5(b) and Figure 5(c) present classification results with the additional negative examples. As can be seen, the performance of the descriptor based methods improves and several of the descriptor methods outperform the similarity based approach. Classification time, however, did not change and remains significantly lower than that required by the similarity based approach.

In order to estimate the effect of varying the number of gallery images per person on the recognition rates, we have varied this parameter. Note that the

Conclusions and future work

We evaluate the performance of a similarity-learning method in comparison to descriptor based methods. The similarity based method performs well on both pair matching and multiple option identification. Descriptor based methods, while performing worse than the similarity learning method when applied directly, can be combined with appropriate learning techniques in order to make better use of the training set and outperform the similarity based method.

There is still much work to be done. Some of the experiments are partial due to the computational complexity of the similarity based method; for the same reason, we were unable to conduct experiments on more than 100 classes. We also wish to study the significance of image alignment to the descriptor based methods. Finally, we are very interested in studying possible combinations of similarity based techniques with descriptor based techniques.

Fig. 1 .

 1 Fig. 1. (a) The LBP image-texture descriptor is computed locally at each pixel location. It considers a small neighborhood of a pixel, and thresholds all values by the central pixel's value. The bits which represent the comparison results are then transformed into a binary number. The histogram of these numbers is used as a signature describing the texture of the image. (b-c) Present an example image from the LFW data set, and its LBP encoding (different intensities representing different codes.)

Fig. 2 .

 2 Fig. 2. (a) The Three-Patch LBP code with α = 2 and S = 8. (b) The TPLBP code computed with parameters S = 8, w = 3, and α = 2. (c) Code image produced from the image in Fig. 1(b).

Fig. 3 .

 3 Fig. 3. (a) The Four-Patch LBP code. Four patches involved in computing a single bit value with parameter α = 1 are highlighted. (b) The FPLBP code computed with parameters S = 8, w = 3, and α = 1. (c) Code image produced from the image in Fig. 1(b).

Fig. 5 .

 5 Fig. 5. Classification results in the multiple-option identification test. The X axis shows the number of classes (a-c) or the number of traninig examples per class (d). The Y axis shows the recognition rate. (a) Training of the descriptor based methods is done with no additional negative examples. (b) 100 additional negative examples are used; (c) 1000 additional negative examples are added to the gallery images during training. (d) performance, with no additional negative training images, as the number of gallery images per class increases (100 classes). The method of [5] is omitted from (d) since the running experiments did not finish in the time of submission. To avoid clutter we refrained from adding error bars. The standard deviations are of the magnitude of 0.05 in all cases. For example, the recognition rate of LBP for 50 classes and no extra negative examples 0.294 ± 0.0098

Table 1 .

 1 Mean (± standard error) recognition rates on the funneled pair matching benchmark of LFW (Image-Restricted Training, "view 1"). Each column represents a different distance measure. Each row corresponds to an image descriptor. The last row ("Combined") corresponds to training a SVM on the vector containing the 8 distances.

	Image Descriptor Euclidian Distance Hellinger Distance
	LBP	0.6767 ± 0.0068	0.6782 ± 0.0063
	Gabor (C1)	0.6293 ± 0.0047	0.6287 ± 0.0046
	TPLBP	0.6875 ± 0.0044	0.6890 ± 0.0040
	FPLBP	0.6865 ± 0.0056	0.6820 ± 0.0055
	Combined	0.7062 ± 0.0057

Acknowledgments

This research is supported by the Israel Science Foundation (grants No. 1440/06, 1214/06), the Colton Foundation, and a Raymond and Beverly Sackler Career Development Chair.

Appendix A: The parameters used in our experiments

Preprocessing For descriptor based methods, all LFW-funneled images used in our tests were cropped to 110 × 115 pixels around their center. Following [START_REF] Heikkilä | Description of interest regions with center-symmetric local binary patterns[END_REF] we further applied an adaptive noise-removal filter (Matlab's weiner2 function) and normalized the images to saturate 1% of values at the low and high intensities. The similarity based method of [START_REF] Nowak | Learning visual similarity measures for comparing never seen objects[END_REF] does not seem to benefit from the preprocessing stage, and we employ it on the original images.

Descriptor parameters Some parameter tuning was done on "view 1" of the LFW dataset, which is intended for such tests. Naive gray-level descriptors are produced by rescaling the cropped and normalized images to half their original size, and sampling all pixels. The image descriptors for all LBP variants are constructed by concatenating histograms produced for 35 non-overlapping blocks of up to 23 × 18 codes. To produce the LBP descriptors we use the MATLAB source code available from [29]. Results are obtained with "uniform" LBP of radius 3 and considering eight samples. The parameters of the patch based LBP descriptors are r1 = 2, S = 8, w = 5 for TPLBP, and r1 = 4, r2 = 5, S = 3, w = 3 for FPLBP.


Giacomo Sorelli

Vyaceslav N Shatokhin

Andreas Buchleitner

Photonic orbital angular momentum in turbulence: vortex splitting and adaptive optics

Keywords: orbital angular momentum, vortex splitting, adaptive optics

Recent works revealed that transmission of light beams carrying orbital-angular-momentum (OAM) through turbulence causes the optical vortex defining these beams to split into multiple vortices with unit topological charge. Here, we consider the numerical propagation of orbital-angular-momentum (OAM) modes through a horizontal atmospheric channel. By analysing the beam's phase front after transmission through turbulence, we confirm the occurence of vortex splitting, but we also witness the emergence of vortex-antivortex pairs. Moreover, by performing performing a decomposition of the transmitted wave into OAM modes, we show that while adaptive optics cannot cancel vortex splitting, it still is pretty efficient in diminishing the turbulence-induced crosstalk between different OAM modes.

INTRODUCTION

High-dimensional states of light bear the possibility of increasing the information capacity and the security of free-space quantum and classical communication. [START_REF] Cerf | Security of quantum key distribution using d-level systems[END_REF] Spanning an unbounded space, the orbital angular momentum (OAM) [START_REF] Allen | Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes[END_REF] of light can be used to realize such high dimensional states. The phase front of a light beam carrying an OAM h (twisted light) is characterised by an optical vortex with topological charge on the beam axis. [START_REF] Dennis | Singular optics: optical vortices and polarization singularities[END_REF] Such a peculiar phase front is however very fragile with respect to turbulence-induced distortions. Hence, it is mandatory to resort to compensation techniques -such as adaptive optics (AO) [START_REF] Leonhard | Protecting the entanglement of twisted photons by adaptive optics[END_REF][START_REF] Sorelli | Entanglement protection of high-dimensional states by adaptive optics[END_REF] -in order to mitigate these effects.

In a recent experiment on the propagation of OAM-carrying beams over a 1.5-km-long intra-city free space link, [START_REF] Lavery | Free-space propagation of high-dimensional structured optical fields in an urban environment[END_REF] and in subsequent theoretical work, [START_REF] Lavery | Vortex instability in turbulent free-space propagation[END_REF] it was suggested that the most detrimental turbulence-induced distortion on twisted light is vortex splitting. By vortex splitting, we intend the separation of an optical vortex of charge , as the one characterizing the phase profile of an optical mode with OAM h , into vortices with unit topological charge. [START_REF] Gbur | Vortex beam propagation through atmospheric turbulence and topological charge conservation[END_REF][START_REF] Ge | Evolution of phase singularities of vortex beams propagating in atmospheric turbulence[END_REF] It was also claimed that traditional AO cannot compensate for vortex splitting, and therefore it is of little help for OAM based free-space communication. The aim of this contribution is to ponder the potential usefulness of AO, its shortcomings notwithstanding.

For this goal, we numerically propagate OAM modes through a horizontal atmospheric channel which is similar to the one considered previously. [START_REF] Lavery | Free-space propagation of high-dimensional structured optical fields in an urban environment[END_REF] By analysing the beams' phase front after transmission through turbulence, we observe that vortex splitting not only occurs but also is often accompanied by the emergence of vortex-antivortex pairs. Finally, we consider AO corrections from the fundamental point of view. Thereby, we exclude from our analysis practical aspects of retrieving the phase of an optical field containing vortices. [START_REF] Fried | Branch point problem in adaptive optics[END_REF]

By

Send correspondence to giacomo.sorelli@lkb.upmc.fr 1 means of an OAM-mode decomposition of the transmitted wave, we show that while AO cannot cancel vortex splitting, it can be pretty efficient in diminishing turbulence-induced crosstalk between different OAM modes. Given that the latter is the most detrimental effect for quantum and classical communication, we conclude that despite vortex splitting, we should not dismiss the benefit of AO in free-space OAM-based communication.

METHODS

We consider a light beam with wavelength λ = 809 nm, whose spatial profile is given by a Laguerre-Gaussian (LG) mode [START_REF] Siegman | Lasers[END_REF] with azimuthal index l 0 = 3, radial index p = 0 and beam waist w 0 = 1.5 cm. For comparison with the experimental results, [START_REF] Lavery | Free-space propagation of high-dimensional structured optical fields in an urban environment[END_REF] we propagate the above mentioned LG mode through a turbulent channel of length z = 1500 m and analyse the received light through a circular aperture of diameter 15 cm. The transverse length-scale of phase fluctuations in a turbulent channel can be characterized by the Fried parameter r 0 = (0.423k 2 C 2 n z) 3/5 , with k = 2π/λ the wave number, and C 2 n the turbulence structure constant. [START_REF] Ishimaru | Wave Propagation and Scattering in Random Media[END_REF][START_REF] Andrews | Laser Beam Propagation Through Random Media[END_REF] When this length scale is much smaller than the typical size of the beam, then propagation will significantly distort the light's spatial profile. Conversely, when this length is much larger than the beam size, then the beam will experience only weak distortions. Accordingly, as a measure of the turbulence strength in our channel, we consider the ratio between the size of our beam and the Fried parameter: W = w 0 /r 0 . Another important parameter is the Rytov variance σ 2 R = 1.23C 2 n k 7/6 z 11/6 , which quantifies the amplitude of turbulence-induced intensity fluctuations, also known as scintillation. [START_REF] Ishimaru | Wave Propagation and Scattering in Random Media[END_REF][START_REF] Andrews | Laser Beam Propagation Through Random Media[END_REF] In particular, the Rytov variance allows to distinguish between the weak (σ 2 R < 1) and strong (σ 2 R > 1) scintillation regimes. 12, 13 For sufficiently short propagation distances, we are always in the weak scintillation regime, where turbulence distortions are dominated by phase fluctuations. [START_REF] Ishimaru | Wave Propagation and Scattering in Random Media[END_REF][START_REF] Andrews | Laser Beam Propagation Through Random Media[END_REF] We exploit this fact to simulate light propagation through turbulence with a multiple-phase-screen approach. [START_REF] Schmidt | Numerical Simulation of Optical Wave Propagation with examples in MATLAB[END_REF][START_REF] Lukin | Adaptive beaming and imaging in the turbulent atmosphere[END_REF] This method consists in splitting a turbulent channel into shorter segments where turbulence-induced phase distortions are introduced via random phase screens, while the evolution between two successive screens is described as free-space propagation in vacuum. In particular, we consider turbulence strengths up to W = 1.5. The weak scintillation at each partial propagation step is verified by splitting our 1.5-km-long channel into 24 segments. We generate the random phase screens with Kolmogorov statistics, defined by the phase power spectral density Φ(κ) = 0.49r -5/3 0 κ -11/3 , using a sub-harmonic method with seven sub-harmonics. [START_REF] Schmidt | Numerical Simulation of Optical Wave Propagation with examples in MATLAB[END_REF][START_REF] Lane | Simulation of a kolmogorov phase screen[END_REF][START_REF] Johansson | Simulation of stellar speckle imaging[END_REF] To assess the efficiency of AO in presence of vortex splitting, we considered two different AO models (see also [START_REF] Leonhard | Protecting the entanglement of twisted photons by adaptive optics[END_REF][START_REF] Sorelli | Entanglement protection of high-dimensional states by adaptive optics[END_REF] ), where phase distortions of a Gaussian beacon are measured and used to correct the LG mode. In the first model, we assume that the phase of the propagated beacon is fully reconstructed with spatial resolution limited only by our numerical grid. This corresponds to an ideal AO system, whose performance is experimentally unachievable. On the contrary, the second AO system only stabilizes the direction of the incoming beam with a tip-tilt mirror. Such a system only compensates for the two lowest-order distortions (tip and tilt) of atmospheric turbulence, which are however dominant in weak turbulence. [START_REF] Tyson | Principles of Adaptive Optics[END_REF] By considering the simplest AO system, as well as its ideal counterpart, we can put a lower and an upper bound on the performance of any realistic AO system.

RESULTS

Let us start by observing how the phase of an OAM-carrying beam evolves upon propagation through the described above turbulent channel for a turbulence strength W = 0.5. In Fig. 1, we see that already in the early stage of the propagation the original optical vortex, with topological charge Q 0 = 3, splits into three spatially separated vortices with charges Q i = 1. When the beam propagates further away from the sender, the vortices wander off the beam axis. Figure 1 also shows us that in addition to vortex splitting, turbulence induces the creation and annihilation of vortex-antivortex pairs, i.e., couples of vortices with opposite topological charge (in general, Q = ±1). The reader should also notice that the fine phase modulations visible in Fig. 1, especially in the first part of the propagations, are numerical artefacts which, however, occur in regions of zero intensity, and do not affect the successive propagation.

Higher-order vortex instability is well known in singular optics [START_REF] Dennis | Singular optics: optical vortices and polarization singularities[END_REF] and it can be understood as follows. Vortices, as all field singularities, can only occur at points of vanishing intensity. Around these points, a general solution of the wave equation featuring a vortex of charge s can be written as

3 ψ(x, y) = a(x + iy) |s| + b(x -iy) |s| , (1)
with a and b complex numbers, and |a| > |b| for s > 0 and |a| < |b| for s < 0. Equation ( 1) tells us that at the vortex position, the field ψ and all its mixed partial derivatives up to order |s| -1 must vanish, while all mixed partial derivatives of order s must be different from zero. Accordingly, |s|(|s| + 1) conditions must be exactly satisfied for an optical vortex of order s to appear, and every perturbation of the field that breaks at least one of this conditions will result in vortex splitting. On the contrary, the sum of the topological charges of all vortices in an optical field,

Q = Nv i=1 Q i , (2)
is a conserved quantity in optical propagation. [START_REF] Dennis | Singular optics: optical vortices and polarization singularities[END_REF][START_REF] Gbur | Vortex beam propagation through atmospheric turbulence and topological charge conservation[END_REF][START_REF] Ge | Evolution of phase singularities of vortex beams propagating in atmospheric turbulence[END_REF] The combination of these two effects, higher-order vortex instability and conservation of Q, explain the behaviour observed in Fig. 1.

To quantify the amount of turbulence-induced vortex splitting, we use the vortex-splitting ratio V, namely the average distance of the vortices from the beam axis divided by the beam radius. Such a quantity was introduced by Lavery. [START_REF] Lavery | Vortex instability in turbulent free-space propagation[END_REF] We modified this definition by weighting the vortices' positions with their topological charges to take into account the creation and annihilation of vortex-antivortex pairs. Formally, we can write the vortex-splitting ratio as

V = 1 N v w 0 Nv i=1 Q i ρ (i) v , (3)
where for each vortex, identified by the index i, ρ

v is the distance from the beam axis, Q i is the topological charge, N v is the total number of vortices in the phase profile, and w 0 is the beam waist.

The left panel of Fig. 2 shows the evolution of the vortex splitting ratio V with the propagation distance z for different turbulence strengths. These results show a monotonic growth of V with a rate that increases with the turbulence strength. The apparent slow-down of the vortex-splitting observed for W = 1.5 close to the receiving plane (purple diamond in the left panel of Fig. 2) is due to the degrading performance of our vortex search routine at large turbulence strengths. [START_REF] Sorelli | Quantum state transfer in diffractive and refractive media[END_REF] Our results are in qualitative agreement with those of Lavery. [START_REF] Lavery | Vortex instability in turbulent free-space propagation[END_REF] However, in the latter case 7 the vortex-splitting ratio grows slower with the propagation distance. This stems from the different numerical simulation strategy, [START_REF] Lavery | Vortex instability in turbulent free-space propagation[END_REF] which underestimate the turbulence strength in the channel. In the right panel of Fig. 2, we present the sum Q (see ( 2)) of all topological charges in the phase profile. As discussed above, the total charge Q must be a conserved quantity. We see that the total charge is correctly predicted in most cases, however, it deviates from its expected value Q = 3 for short propagation distances, as well as high turbulence strengths.

The underestimation of Q that we observe for z = 125 m is a simple spatial resolution issue. In fact, so close to the sender, even though the vortex splitting already happened, the vortices are still not well-separated (see the top left panel in Fig. 1), and our detection algorithm is unable to distinguish them. However, these unresolved vortices do not critically influence the estimation of the small vortex-splitting ratio close to the sender.

For larger turbulence strengths, there are two effects that may lead to an incorrect determination of the total topological charge. First, to reduce the sensitivity of our results with respect to numerical noise at the boundary of the aperture (where the amplitude should vanish, but may be garnished by finite values of the order of the numerical precision), we search for vortices only in a finite region around the beam axis. For large turbulence strengths, it may happen that a vortex wanders away from this region and therefore remains undetected. This physical effect causes underestimation of the vortex-splitting ratio and of the total topological charge. Second, at larger turbulence strengths, the phase acquires finer spatial structures which are difficult to resolve numerically and which lead to failures of the vortex detection routine. [START_REF] Sorelli | Quantum state transfer in diffractive and refractive media[END_REF] Numerical errors may result in both under-and overestimation of the total topological charge (see Fig. 2).

The above discussion proved that our numerical routine is able to correctly describe vortex splitting within a given range of turbulence conditions. We now aim at analysing how this phenomenon affects the capability of AO to protect light beams carrying OAM. A previous work 7 considered this problem and concluded that AO cannot be efficient in presence of vortex splitting. However, we believe that the poor AO-performance observed in this work [START_REF] Lavery | Vortex instability in turbulent free-space propagation[END_REF] is mainly due to the unfavorable diffraction conditions that he assumed and not by vortex splitting. To support this point, we slightly modified our simulation parameters to mitigate diffraction effects. In particular, we increased the beam waist to w 0 = 0.03 m, and the receiver aperture to 0.2 m, leaving all other parameters unchanged.

Our results, for a turbulence strength W = 1, are plotted in Fig. 3. The first row in Fig. 3 shows the stability of the intensity ring of the LG modes that we discussed in previous work [START_REF] Sorelli | Entanglement protection of high-dimensional states by adaptive optics[END_REF] and which was observed experimentally. [START_REF] Krenn | Twisted light transmission over 143 km[END_REF] Let us have a closer look at the phase distributions. In the uncorrected case (panel (d)), we observe a highly distorted phase profile where the three vortices originating from the splitting of the original Q = 3 vortex are accompanied by two vortex-antivortex pairs. Tip-tilt corrections (panel (e)) do not significantly reduce the phase distortions, and leave the vortices unaltered. In panel (f), we see that, as expected, also the ideal correction does not modify the positions of the vortices, but it significantly modifies the phase distribution so that it closely resembles the one of a mode which propagated in the absence of perturbations. The observation that for an ideally corrected beam both intensity and phase resemble those of an unperturbed LG mode suggests that AO can still provide some benefits, despite vortex splitting.

To make this conclusion more quantitative, we decomposed the propagated field in the basis of LG modes with p = 0. How the intensity of the input l 0 = 3 LG mode spreads over this basis is plotted in panels (g-i) of Fig. 3, displaying the probability distribution P(l 0 → l) of finding the output field in the LG mode with the azimuthal index l. Comparing the different panels, we see that already tip-tilt corrections (panel (h)) mitigate the spreading of the OAM distribution by depleting the populations of those modes which are further away from the input modes. When using ideal AO (panel (i)) the populations of modes different from the input one, are substantially negligible. We also see that the average OAM of the distribution is significantly lower than the one of the input mode (l 0 = 3), while it is close to this value in the tip-tilt corrected case and is exactly l = 3 in the ideally corrected one. We point out that this shift of the mean OAM towards a lower value, which cannot be reproduced by weak turbulence models (e.g. single phase screen models [START_REF] Paterson | Atmospheric turbulence and orbital angular momentum of single photons for optical communication[END_REF][START_REF] Tyler | Influence of atmospheric turbulence on the propagation of quantum states of light carrying orbital angular momentum[END_REF][START_REF] Leonhard | Universal entanglement decay of photonic-orbitalangular-momentum qubit states in atmospheric turbulence[END_REF] ), was the experimental observation 6 that triggered the investigation of vortex splitting. [START_REF] Lavery | Vortex instability in turbulent free-space propagation[END_REF] Finally, in panels (g -i) of Fig. 3, we also included the average vortex splitting ratio, and we see that it is substantially unaffected by the corrections * . This confirms our idea that AO can be very efficient in correcting OAM modes even in the presence of vortex splitting.

CONCLUSIONS

We investigated the impact of vortex splitting on the efficiency of AO protection of OAM modes propagated through atmospheric turbulence. We proved that our numerical simulation routine can accurately reproduce the monotonic increase of vortex splitting with the propagation distance and the turbulence strength. [START_REF] Lavery | Vortex instability in turbulent free-space propagation[END_REF] It was also suggested [START_REF] Lavery | Vortex instability in turbulent free-space propagation[END_REF] that given its monotonicity, the vortex splitting ratio (3) could be an appropriate measure to quantify the turbulence-induced distortions of OAM beams. Although this idea deserves further investigations, our results show that vortex splitting is accompanied by the generation of vortex-antivortex pairs, which wander through the phase profile and are not distinguishable from the vortices originating from the splitting of a higher-order vortex. These pairs of vortices with opposite topological charge may affect the monotonicity of the vortex-splitting under strong turbulence. Ultimately, we observed that -even in presence of vortex splitting -AO remains capable of mitigating turbulence-induced crosstalk between different OAM modes, which is the most detrimental effect for quantum and classical communication applications.

Figure 1 .

 1 Figure 1. Phase profile of an l0 = 3 LG mode transmitted through a turbulent channel (turbulence strength W = 0.5) recorded every 125 m between the sender (at z = 0) and the receiver (at z = 1500 m). The dots indicate the location of optical vortices, with colour-encoded topological charge: black for Q = +1 and red for Q = -1.

Figure 2 .

 2 Figure 2. (Left) Vortex-splitting ratio V (see (3)) as a function of the propagation distance z for different turbulence strengths W . (Right) Total topological charge (see (2)) plotted against the propagation distance z for the same values of the turbulence strengths W as in the left plot. Both panels result from the propagation of an LG mode (l0 = 3, p = 0) averaged over 30 random realizations of the turbulent channel. Error bars are standard deviations from the mean.

Figure 3 .

 3 Figure 3. (a-c) Intensity and (d-f) phase distributions of an LG mode with azimuthal index l0 = 3 and radial index p = 0 transmitted through a turbulent channel with turbulence strength W = 1. The dots in the phase panels represent the positions of the optical vortices with color-coded topological charge: black for Q = +1 and red for Q = -1. Panels (g-i) show a decomposition of the received field into p = 0 LG modes averaged over 100 realizations of the turbulence channel together with the average OAM (in units of h) l of the distribution, and the average vortex splitting ratio V. The different columns correspond to different degrees of AO corrections: panels (a), (d), (g) no correction, (b), (e), (h) tip-tilt correction, (c), (f), (i) ideal AO.

* The slightly different value in the ideally corrected case can be imputed to numerical noise.

Acknowledgements

VS acknowledges support by the Strategiefonds der Albert-Ludwigs-Universitaet Freiburg.