Jean Gilles

Yukiko Denis

Laurent Hugues

Nicolas Marie-Pierre Marc

Applications of Digital Topology For Real-Time Markerless Motion Capture

Keywords: markerless, motion capture, digital topology, thinning, alignment, homeomorphism, betweenness

This manuscript deals with the problem of markerless motion capture. An approach to this problem is model-based and is divided into two steps: an initialization step in which the initial pose is estimated, and a tracking which computes the current pose of the subject using information of previous ones. Classically, the initialization step is done manually, forbidding the possibility to be used online, or requires constraining actions of the subject.

We propose an automatic real-time markerless initialization step, that relies on topological information provided by skeletonization of a 3D reconstruction of the subject. This topological information is then represented as a tree, which is matched with another tree used as model description, in order to identify the different parts of the subject. In order to provide such a method, we propose some contributions in both digital topology and graph theory research fields.

As our method requires real-time computation, we first focus on the speed optimization of skeletonization methods, and on the design of new fast skeletonization schemes providing good results.

In order to efficiently match the tree representing the topological information with the tree describing the model, we propose new matching definitions and associated algorithms.

Finally, we study how to improve the robustness of our method by the use of innovative constraints in the model. This manuscript ends by a study of the application of our method on several data sets, demonstrating its interesting properties: fast computation, robustness, and adaptability to any kind of subjects.

Titre:

Applications de la topologie discrète pour la captation de mouvement en temps réel et sans marqueurs.

Résumé

Durant cette thèse, nous nous sommes intéressés à la problématique de la captation de mouvement sans marqueurs. Une approche classique est basée sur l'utilisation d'un modèle prédéfini du sujet, et est divisée en deux phases: celle d'initialisation, où la pose initiale du sujet est estimée, et celle de suivi, où la pose actuelle du sujet est estimée à partir des précédentes. Souvent, la phase d'initialisation est faite manuellement, rendant impossible l'utilisation en direct, ou nécessite des actions spécifiques du sujet.

Nous proposons une phase d'initialisation automatique et temps-réel, utilisant l'information topologique extraite par squelettisation d'une reconstruction 3D du sujet. Cette information est représentée sous forme d'arbre (arbre de données), qui est mis en correspondance avec un arbre utilisé comme modèle, afin d'identifier les différentes parties du sujet. Pour obtenir une telle méthode, nous apportons des contributions dans les domaines de la topologie discrète et de la théorie des graphes.

Comme notre méthode requiert le temps réel, nous nous intéressons d'abord à l'optimisation du temps de calcul des méthodes de squelettisation, ainsi qu'à l'élaboration de nouveaux algorithmes rapides fournissant de bons résultats.

Nous nous intéressons ensuite à la définition d'une mise en correspondance efficace entre l'arbre de données et celui décrivant le modèle.

Enfin, nous améliorons la robustesse de notre méthode en ajoutant des contraintes novatrices au modèle.

Nous terminons par l'application de notre méthode sur différents jeux de données, démontrant ses propriétés: rapidité, robustesse et adaptabilité à différents types de sujet.

Mots Clés

captation de mouvement sans marqueurs; topologie discrète; squelettisation; alignement; homeomorphisme; interposition. Global Introduction

List of Figures

"Notre nature est dans le mouvement, le repos entier est la mort."

("Our nature consists in motion; complete rest is death.")

Blaise Pascal.

Background

Motion capture consists in automatic estimation of the pose (i.e. relative position and orientation of each part of the subject), motion and actions of a subject, usually a full human being. It is an highly active research field since more than twenty years, notably due to its industrial applications in computer-animated feature films and Human-Computer Interfaces (HCI). For the latest, some specific properties are required:

• Real-time computation, in order to provide good interaction.

• The subject must not need to wear specific equipment (markers).

Usual markerless motion capture systems are composed by a set of cameras, providing the input information of the subject. Numerous markerless methods of motion capture use an a priori model describing the structure and shape of the subject whom motion has to be captured. Such methods are composed of at least two steps: the initialization step, where the initial pose and the parameters of the model are detected, and the tracking step, consisting of finding the current pose in function of the previous ones.

In many of such methods, the initialization step is done manually by the user, forbidding the possibility to be used for HCI purpose, or requires constraining actions or poses of the subject.

Contribution

In this thesis, we propose an automatic real-time markerless initialization step using a very simple a priori model, described by a tree (called model tree), which can be easily designed in order to use the method for any kind of subject, e.g. full human being, hand, or dog.

Our method uses topological information (the skeleton) extracted from a 3D reconstruction of the subject, obtained by using a multi-camera system. This topological information is then represented by a tree (called data tree), which will be matched with the model, in order to find the position of specific parts, e.g. head, torso, crotch, arms and legs in the case of full human being.

In order to reach our goal, we have done some contributions in different fields of research:

Digital Topology. As our method requires real-time computation, we have worked on optimizations of the speed of skeletonization algorithms and on design of new fast and efficient skeletonization schemes.

Graph Theory. The matching between the model tree and the data tree has to take into consideration some specific noises on the data tree, due to its acquisition. We have proposed a definition of matching that is adapted to our needs, and associated computationally efficient algorithms.

Motion Capture. Merging the results of our previous works, we provided an efficient automatic real-time markerless initialization step. We also introduced novel intuitive constraint definitions, which can be optionally added to the model definition, in order to improve the robustness of the method.

The structure of this manuscript is in four parts: the first one introduces the background and the backbone of our method, and the three others deal with our contributions and results in the three fields of research enumerated above, respectively.

Chapter 2

Motion Capture Overview

Automatic capture and analysis of articulated 3D subjects (e.g. humans, hands or animals) motion is an highly active research area since more than twenty years.

Motion capture presents issues in several kinds of applications:

3D models animation: motion capture is widely used for movies and video games characters animation. Usage of real actors movements to animate virtual avatars is faster, more accurate and realistic than doing it manually.

Human-computer interfaces (HCI): several usages can be done of motion capture in the field of HCIs, including domotic (control of TV or HIFI by hand moves), sign language understanding, or virtual reality system interactions.

Surveillance: motion capture can be used in surveillance field in order to detect suspect behaviors, e.g. fights or theft, or crowd analysis.

Medical analysis: motion capture can be used for medical purpose, in order to analyze gait problems or high level sportsmen performances.

Of course, all these applications do not require the same properties for motion capture. Animation and medical analysis request a very accurate motion capture and analysis in order to be usable, but it can be done offline (i.e. after the performance of the subject). On the other hand, HCI and surveillance applications need a real-time motion capture, in order to allow an immediate response for some given situations but the accuracy is less essential, and only requires to allow the understanding of the action or the approximate location of the subject parts.

There exists numerous different methods to perform motion capture. In the sequel of this chapter, we propose a non exhaustive state of the art on motion capture, focusing on the class of methods used in the thesis, and on their specificities. For exhaustive studies, the reader can refer to [START_REF] Moeslund | A survey of computer vision-based human motion capture[END_REF][START_REF] Moeslund | A survey of advances in vision-based human motion capture and analysis[END_REF][START_REF] Poppe | Vision-based human motion analysis: An overview[END_REF].

Overview of Motion Captures Systems

Acquisition Systems Using Markers

Several acquisition systems estimate motion of a subject by measuring the position and/or the orientation of physical objects (called markers), which are fixed on the subject at specific locations. Acquisition systems with markers are widely proposed as commercial solutions. Different classes of such system can be defined, in regard of the method used to measure marker positions:

• Mechanical systems (e.g. Gypsy TM system, see Figure 2.1a) estimate the motion by direct measurement of articulation angles, via potentiometers fixed on an exoskeleton. The main problem of such systems is the fragility of the skeleton, which can be easily broken.

• Magnetic systems (e.g. MotionStar TM and Fastrak TM systems) consist in generating an electromagnetic field then finding the position and orientation of markers (coils of electric wire) by their perturbation of the field. The main problem of such system is the sensitivity to interferences from metal objects of the environment.

• Inertial systems (e.g. Colibri TM and Xsens TM systems, see Figure 2.1b) use markers composed by accelerometers and gyroscopes. These sensors provide motion informations, as speed and rotation of the marker. Then, softwares are used in order to find the pose of the subject from these informations. A very famous inertial motion capture system is the Wiimote TM . The main problem of such system is the lack of accuracy, in regard of other methods.

• Optical Systems can be classified in two categories:

-Passive optical systems (e.g. OptiTrack TM , see Figure 2.1c, and Qualisys TM systems) use cameras in order to track markers coated with a retroreflective material reflecting light generated near the cameras lens.

-Active optical systems (e.g. PhaseSpace TM system) use markers powered to emit their own light (via a LED).

The main problem of such systems is the sensibility to self occlusions: all the markers are not detected by all cameras, some of them being hidden by parts of the body.

These acquisition systems using markers have in common to be very expensive, due to the requested specific hardware.

Marker Free Optical Acquisition Systems

Contrary to other acquisition systems, marker free optical acquisition systems only require cameras. We can classify marker free methods in regard of the requested number of cameras: Monocular Methods. Pose estimation from only one video stream is a very complex problem, notably due to self occlusions of subject parts. This kind of method is mainly used for surveillance purpose, crowd study and hand gesture recognition.

Multi-view Methods. During the last twenty years, the increase of computer power has allowed the management of multiple synchronized cameras and fast computation of 3D representation of the subject. It led to the development of numerous multi-view methods, freed from the problem of self occlusions. A majority of these methods requires a static background and a prior camera calibration.

According to Menier et al. [START_REF] Menier | 3d skeleton-based body pose recovery[END_REF], most of the markerless methods fall into three categories:

• Learning-based methods [START_REF] Agarwal | 3D human pose from silhouettes by relevance vector regression[END_REF][START_REF] Alisi | Motion Capture Based on Color Error Maps in a Distributed Collaborative Environment[END_REF]9,[START_REF] Grauman | Inferring 3d structure with a statistical image-based shape model[END_REF][START_REF] Mori | Estimating human body configurations using shape context matching[END_REF][START_REF] Shakhnarovich | Fast pose estimation with parameter-sensitive hashing[END_REF][START_REF] Stenger | Filtering using a tree-based estimator[END_REF][START_REF] Tanaka | Human pose estimation from volume data and topological graph database[END_REF][START_REF] Urtasun | Temporal motion models for monocular and multiview 3D human body tracking[END_REF] store a set of training examples whose poses are known and estimate pose by searching for training images similar to the input images. This category of methods requires larger storage memory space than others.

• Model-free methods [START_REF] Brostow | Novel skeletal representation for articulated creatures[END_REF][START_REF] Chu | Markerless kinematic model and motion capture from volume sequences[END_REF][START_REF] Cuzzolin | Robust spectral 3D-bodypart segmentation along time[END_REF][START_REF] Cuzzolin | Coherent Laplacian 3-D protrusion segmentation[END_REF][START_REF] Drouin | Incremental discovery of object parts in video sequences[END_REF][START_REF] Mateus | Articulated shape matching using locally linear embedding and orthogonal alignment[END_REF][START_REF] Mateus | Articulated shape matching using laplacian eigenfunctions and unsupervised point registration[END_REF][START_REF] Starck | Spherical matching for temporal correspondence of non-rigid surfaces[END_REF][START_REF] Varanasi | Temporal surface tracking using mesh evolution[END_REF] do not require any a priori knowledge and automatically recover articulated structures.

• Model-based methods use an a priori model which is tracked using image information.

Moeslund et al. [START_REF] Moeslund | A survey of computer vision-based human motion capture[END_REF][START_REF] Moeslund | A survey of advances in vision-based human motion capture and analysis[END_REF] defined a general structure of such methods in four steps:

1. Initialization. Ensuring that the system starts motion capture with a correct interpretation of the scene: initial subject pose, background images, model parameters are some of the possible informations acquired during this step.

2.

Tracking. Using temporal information in order to find the current parameters of the model (e.g. position and orientation of the different parts) from the previous ones.

3. Pose estimation. Estimating the pose of the subject in function of the model parameters.

4.

Recognition. Recognizing the actions, activities and behaviors of the subject.

Classically, tracking and pose estimation steps are merged in one step. Here, we consider that recognition step is a problem distinct from motion capture.

Since in this thesis, we exclusively deal with markerless multi-view model-based methods, we will focus on the main characteristics of this kind of method.

A priori Model Definitions

Numerous model definitions have been proposed in the literature. A wide majority of them includes a kinematic structure (also called kinematic skeleton) defined by:

• a set of segments with specified lengths, representing limbs,

• a set of joints linking segments, with specified degrees-of-freedom, representing articulations.

Although some methods only require a kinematic structure [START_REF] Cheung | Shape-from-silhouette of articulated objects and its use for human body kinematics estimation and motion capture[END_REF][START_REF] Menier | 3d skeleton-based body pose recovery[END_REF], most of them require to link it with shape information.

A usual way to provide shape information is the use of geometric primitives, like ellipsoids [START_REF] Caillette | Real-time markerless human body tracking with multi-view 3d voxel reconstruction[END_REF],

truncated elliptical cones [START_REF] Knossow | Human motion tracking with a kinematic parameterization of extremal contours[END_REF][START_REF] Knossow | Tracking with the kinematics of extremal contours[END_REF][START_REF] Knossow | Articulated-Body Tracking Through Anisotropic Edge Detection[END_REF], super-quadric [START_REF] Sundaresan | Markerless motion capture using multiple cameras[END_REF][START_REF] Sundaresan | Multi-camera tracking of articulated human motion using motion and shape cues[END_REF], or combinations of several kinds [START_REF] Delamarre | 3D Articulated Models and Multiview Tracking with Physical Forces* 1[END_REF][START_REF] Michoud | Real-time and marker-free 3D motion capture for home entertainment oriented applications[END_REF][START_REF] Michoud | Real-Time Marker-free Motion Capture from multiple cameras[END_REF][START_REF] Mikić | Articulated body posture estimation from multi-camera voxel data[END_REF][START_REF] Mikić | Human body model acquisition and tracking using voxel data[END_REF].

Another usual way to provide shape information is to use a surface, which can be generated from geometric primitives [START_REF] Horaud | Human motion tracking by registering an articulated surface to 3D points and normals[END_REF][START_REF] Plänkers | Articulated soft objects for video-based body modeling[END_REF] (implicit surface), or an accurate mesh acquired in a prior step by a specific device [START_REF] Gall | Motion capture using joint skeleton tracking and surface estimation[END_REF][START_REF] Kehl | Full body tracking from multiple views using stochastic sampling[END_REF]. Notice that in some cases [START_REF] Aguiar | Performance capture from sparse multi-view video[END_REF][START_REF] Gall | Optimization and filtering for human motion capture[END_REF], such meshes are used without kinematic skeletons.

In addition to kinematic and shape information, color information can be used, by texturing the mesh [11], or giving some additional color priors [START_REF] Michoud | Real-time and marker-free 3D motion capture for home entertainment oriented applications[END_REF] (e.g. noticing that hands and head are the only parts to have the same color).

Adaptability of the Models

The adaptability of the method, i.e. the simplicity to modify the method in order to use it for another subject than the one initially targeted, obviously depends on the model definition. We define two levels of adaptability:

• inner-class adaptability is the simplicity to modify the model in order to use it for another subject of the same kind (e.g. two different human beings)

• trans-class adaptability is the simplicity to modify the model in order to use it for a subject of different kind (e.g. from full human body to hand motion capture) Models using only kinematic structure, geometric primitives and implicit surfaces have a high inner-classes adaptability, model parameters being adjusted during initialization step, contrary to meshes, which have to be rebuild in a prior step, using specific devices and methods.

Concerning the trans-classes adaptability, almost all the methods can be used for any kind of subject, if an adapted model is provided (such methods are called generic). The exception is the case of color priors [START_REF] Michoud | Real-time and marker-free 3D motion capture for home entertainment oriented applications[END_REF] and methods requesting a specific detection for each part of the subject [START_REF] Michoud | Real-time and marker-free 3D motion capture for home entertainment oriented applications[END_REF][START_REF] Mikić | Articulated body posture estimation from multi-camera voxel data[END_REF][START_REF] Mikić | Human body model acquisition and tracking using voxel data[END_REF]. For all the others, the trans-classes adaptability is inversely proportional to the complexity of the model.

On the Use of Multi Camera Systems

Motion capture methods use information provided by multi-camera systems in many different ways, and can be classified into two categories: methods using 3D reconstruction of the subject and methods performing projection of the model in the images.

Projection in Images

Projection of the model in images can be performed by several methods, in order to use specific information which is easier and faster to extract in 2D than in 3D. Examples of information are:

• silhouette [START_REF] Gall | Optimization and filtering for human motion capture[END_REF] and contours of the subject [START_REF] Delamarre | 3D Articulated Models and Multiview Tracking with Physical Forces* 1[END_REF][START_REF] Knossow | Human motion tracking with a kinematic parameterization of extremal contours[END_REF][START_REF] Knossow | Tracking with the kinematics of extremal contours[END_REF][START_REF] Knossow | Articulated-Body Tracking Through Anisotropic Edge Detection[END_REF], methods searching for a optimal fitting of the model in each of them.

• connected components of a same color [START_REF] Aghajan | Distributed Vision Networks for Human Pose Analysis[END_REF], methods using them to facilitate the fitting.

• optical flow of contours [START_REF] Sundaresan | Multi-camera tracking of articulated human motion using motion and shape cues[END_REF], specific projected points [11] or shape descriptors [START_REF] Aguiar | Performance capture from sparse multi-view video[END_REF][START_REF] Gall | Motion capture using joint skeleton tracking and surface estimation[END_REF] (e.g.SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF]) can be performed in order to track the motion in all the images independently, then reconstruct the 3D motion by triangulation.

3D Reconstruction of the Subject

Another approach consists of directly using 3D clues, considering 3D data resulting from multiview modeling methods.

Stereo-Correlation

A famous binocular modeling method is the stereo-correlation, used in some methods [START_REF] Delamarre | 3D Articulated Models and Multiview Tracking with Physical Forces* 1[END_REF][START_REF] Plänkers | Articulated soft objects for video-based body modeling[END_REF]:

when cameras have a small baseline and are pointing in about the same direction, a part of the surface of the subject is visible in both images. It is thus possible to estimate the depth of this surface part by studying the disparity between corresponding points in the two images.

Visual Hulls

The most used multi-view modeling method is the visual hull, introduced by Laurentini [START_REF] Laurentini | The visual hull concept for silhouette-based image understanding[END_REF]:

intuitively, an object lies inside the volume generated by back-projecting its silhouette through the camera center (called silhouette cone). The intersection of the silhouette cones of all the cameras results in a volume called the visual hull (see Figure 2.3 for some illustrations), which is guaranteed to contain the object. Two kind of visual hulls can be distinguished:

• surface visual hull, used in [START_REF] Menier | 3d skeleton-based body pose recovery[END_REF], consists in only representing the surface of the visual hull by a mesh. Some methods, e.g. [START_REF] Horaud | Human motion tracking by registering an articulated surface to 3D points and normals[END_REF], use in addition the information from normals to surface elements.

• volumetric visual hull, used in [START_REF] Mikić | Articulated body posture estimation from multi-camera voxel data[END_REF][START_REF] Mikić | Human body model acquisition and tracking using voxel data[END_REF], consists in representing the visual hull volume by a subset of a voxel grid positioned at the intersection of all images cones: a voxel belongs to the visual hull if and only if the projection of it center in each image belong the corresponding silhouette. The voxels can be colorized [START_REF] Caillette | Real-time markerless human body tracking with multi-view 3d voxel reconstruction[END_REF][START_REF] Cheung | Shape-from-silhouette of articulated objects and its use for human body kinematics estimation and motion capture[END_REF][START_REF] Kehl | Full body tracking from multiple views using stochastic sampling[END_REF][START_REF] Michoud | Real-time and marker-free 3D motion capture for home entertainment oriented applications[END_REF][START_REF] Michoud | Real-Time Marker-free Motion Capture from multiple cameras[END_REF] using ray tracing from the cameras, providing additional information. Volumetric visual hull can be computed in real-time, using GPU [START_REF] Hasenfratz | Real-time capture, reconstruction and insertion into virtual world of human actors[END_REF], or hierarchical carving [START_REF] Caillette | Real-time markerless human body tracking with multi-view 3d voxel reconstruction[END_REF].

Notice that some methods [START_REF] Sundaresan | Markerless motion capture using multiple cameras[END_REF] use both visual hull and projection in images.

Initialization Step of Model-Based Methods

A wide part of model-based motion capture methods found in the literature focus on the tracking and pose estimation step, and often omit the initialization step.

The initialization step is though essential, as providing the parameters of the model (limb lengths, dimensions of geometric primitives, color information, initial position and orientation of the different parts) which are primordial for an efficient tracking.

Different solutions have been proposed in the literature:

1. manual initialization [START_REF] Menier | 3d skeleton-based body pose recovery[END_REF] consists of manually setting the initial values of the model.

2. limb by limb identification [START_REF] Cheung | Shape-from-silhouette of articulated objects and its use for human body kinematics estimation and motion capture[END_REF] consists of asking to the subject to successively move only each of its limbs, in a specific order, in order to easily identify which part of subject information (in the case of Cheung et al. [START_REF] Cheung | Shape-from-silhouette of articulated objects and its use for human body kinematics estimation and motion capture[END_REF], a visual hull) is associated to each part of the model.

3. some methods require a specific initial pose [START_REF] Caillette | Real-time markerless human body tracking with multi-view 3d voxel reconstruction[END_REF][START_REF] Cheung | Shape-from-silhouette of articulated objects and its use for human body kinematics estimation and motion capture[END_REF][START_REF] Kehl | Full body tracking from multiple views using stochastic sampling[END_REF], most often an "X"-pose, where the subject have spread legs and arms. It allows a fast identification of the model parameters, as the limb positions are already approximately set.

4. the initialization can be automatically done, using a hierarchic fitting [START_REF] Michoud | Real-time and marker-free 3D motion capture for home entertainment oriented applications[END_REF][START_REF] Michoud | Real-Time Marker-free Motion Capture from multiple cameras[END_REF][START_REF] Mikić | Articulated body posture estimation from multi-camera voxel data[END_REF]: due to its unique shape and size, the head is the easiest part to find and is identified first. Then, the torso is located using head position information, and finally limbs are detected using torso position and orientation information.

The method 1 is difficult to use online, for HCI purpose for example. It is even more a shame, as some methods using it (e.g [START_REF] Menier | 3d skeleton-based body pose recovery[END_REF]) are able to perform real-time tracking.

On the other hand, methods 2 and 3 can be used online, but such methods involve some constraints for the subject, which have to perform special actions in order to be tracked.

Finally, method 4 can be used online and does not require any specific action or pose from the subject, but is not generic: if the subject does not contain a part which is easy to detect, the initialization cannot be performed.

For our best knowledge, no real-time generic initialization step method has been proposed in the literature.

Chapter 3

Motivation

As seen in the previous chapter, numerous methods of marker-free 3D motion capture require a manual initialization of the subject pose or a constraining one, e.g. needing a special pose of the subject.

The aim of our work is to provide a pose initialization method (i.e. resulting only in the positions of specific parts of the subject) with several properties. First, our method has to be generic, in the sense it must be adaptable for a large set of subjects (full human body, hands, upper part of the human body, animals...) only by creating a very simple a priori model. Our method has to be real-time, in order to be usable with real-time tracking for online applications.

Furthermore, we request our method to be robust, i.e. to work for a large set of subject poses, avoiding constraints for the subject. Finally, we consider that the input of the method is only a volumetric visual hull of the subject.

This chapter is composed of two sections. In the first one, we discuss the definition of the a priori model required by our method. Then, in the second section, we propose a brief overview of the pipeline of our method.

Model Definition

Constraints of Descriptors

The choice of the descriptors used for our a priori model definition is linked to different observations, according to the nature of the application and to our philosophy.

As we consider only the visual hull in input, we do not use any color information, contained in the camera images. No color prior is used in our method. Furthermore, some models have homogeneous color, e.g. hands, thus such information is not enough generic.

13

The quality of the 3D shape depends of the number of views for the visual hull reconstruction, and the respective positions of the cameras. In a multi-view system without enough cameras, or not well positioned ones, the visual hull can present important variations of limbs thickness in regard of the subject. Due to this fact, we do not use accurate shape representations, like geometric features.

As we want to create a method that should be usable with any kind of subject, we have to propose a representation allowing to easily describe the subject. We thus avoid the use of shape primitive description and laser scanned models.

Our method has to be able to find the pose of the subject in a maximum number of situations, in order to be able to perform the initialization step as often as possible. By this way, we can use our method even when the subject does not perform a specific action or pose (like it is usually the case), and we can check the tracking at closer intervals, increasing the accuracy.

The weakest constraint we found is that the topology of the shape and its salient geometric features have to be the same as for the model.

Model Description

Taking into consideration all the constraints enumerated in the previous section, we propose an a priori model providing a description of the topology and of the geometric features of the subject, using a tree representation.

The model is built as follows:

• Each vertex represents a characteristic part of the model. For example:

-for the full human body: hands, feet, head, torso and crotch.

-for the hand: the arm, the palm (in two parts), the end of the fingers.

• Two vertices are linked if the associated parts are adjacent in the model.

• To each edge is associated a weight, corresponding to the approximative distance (for a given distance unit used for all the weights in the graph) between the two linked parts.

Two kinds of models can be considered (see Figure 3.1 for examples):

• the complete models, for subjects fully contained in the 3D acquisition space, as in the case of full human body motion capture, for example.

• the incomplete models, which are part of a biggest shape, as in the case of hand pose estimation. In this case, a part of the shape intersects the border of the 3D acquisition space (e.g. the arm, in case of hand motion capture), and we represent this part in the model as an infinite limb, as we do not know it length.

: H-T=4 T-C=10 T-H1=12 T-H2=12 C-F1=16 C-F2=16 H T C A 1 A 2 F 1 F 2 A P F T I M R L Vertices: H T C A1 A2 F1 F2
Vertices:

Figure 3.1: Description of two different models. On the Left, the hand model (incomplete model): notice that the edge between the arm and the palm is weighted by "infinity", as we do not know its length. On the right, the full human body model (complete model).

Notice that the use of a tree, instead of a graph, for the model description, forbid the representation of models containing holes (represented by a loop in a graph), but allow faster algorithms for the matching with data.

Initialization Method Overview

Now we have defined the inputs of our method: a visual hull representation of the subject and the a priori model.

Our method consists in extracting a descriptor from the visual hull (called the data descriptor ) which is similar to the one used in the a priori model (called the model descriptor ). The main difference is that the data descriptor has to contain some spatial information, where as that the model descriptor contains informations about the nature of the parts of the model.

Finding a matching between the two descriptors results in finding spatial information for each part of the model and by this way, the pose of the subject.

Extraction of Data Tree

First we have to find a way to extract topological and geometrical information from the visual hull, and to represent it by a tree. We solve this problem in two steps: first, we consider the skeleton of the visual hull. The skeleton is an efficient shape descriptor, representing in a compact form the topology and the main geometric features of a shape. It can be obtained by iteratively removing voxels without changing the topology (see Figure 3.2 for an example).

Then, we "translate" the skeleton into a tree: • Each vertex represents a characteristic point of the skeleton: ending point or intersection point.

• Two vertices are linked if the associated points are adjacent to a same curve in the skeleton.

• To each edge is associated a weight, which is the number of voxels in the corresponding curve.

In the case of incomplete model, we have to consider that a branch ended by a point touching the border of acquisition space have an infinite length. See

Matching Data Tree with Model Tree

Once we have a tree description of the subject, we have to find a way to match it with the model tree. In the literature, numerous methods are proposed to find a matching between two trees. However, in our case we have to take into consideration several kind of noises in the data tree. Furthermore, in the literature, proposed methods are designed for special kind of trees, usually rooted, oriented, and with information on vertices instead of on edges, as in our case.

We have to define a new way to match our data tree with our model tree, corresponding to our constraints.

Pipeline and Main Difficulties of Our Method

The pipeline of our method is summarized in Figure 3.4. We have now defined the different steps of our method, however, two main difficulties remain.

Visual Hull

First, the skeletonization process has to be both speed efficient and robust to visual hull surface noise. To our best knowledge, no skeletonization algorithm satisfies these two properties. In the second part of this manuscript, we propose some solutions to speed up existent algorithms and a new algorithm providing good skeletons in real time.

The second problem is the matching with the model tree. As explained in Section 3.2.2, no method found in the literature is adapted to our special case of unrooted edge weighted trees, with special noises on the data tree. In the third part of this manuscript, we propose new kind of matching specially designed for our problem, and efficient algorithms for their computation.

Finally, after solving these two main difficulties, we propose in the last part of the manuscript to improve the robustness of the result by addition of intuitive constraints, and we study the results and the speed of our method.

Part II

Skeletonization Optimizations

Chapter 4

State of the Art of Skeletonization

The skeleton was originally defined by Blum [START_REF] Blum | An associative machine for dealing with the visual field and some of its biological implications[END_REF] based on a " grass fire" analogy. Imagine a shape as a field covered by dry grass; if you set on fire the contour of the field, then the meeting points of the flame fronts would constitute the skeleton of the shape. In the continuous framework, this definition is equivalent to saying that the skeleton is the set of points which are centers of maximal balls (balls included in the object, and not strictly included in any other such ball) ( [START_REF] Calabi | A study of the skeleton of plane figures[END_REF]). In 1969, Hilditch gave four properties that a skeleton in a bi-dimensional space should possess [START_REF] Hilditch | Linear skeletons from square cupboards[END_REF]. Adapted to the general case of n-dimensional skeletons, these properties are:

1. a skeleton should be homotopic to the original object, In the continuous framework, the set of centers of maximal balls, called the medial axis, satisfies these properties [START_REF] Lieutier | Any open bounded subset of n has the same homotopy type as its medial axis[END_REF][START_REF] Matheron | Examples of topological properties of skeletons[END_REF]. In the discrete framework Z n , the discrete medial axis does not satisfy two of these properties: it is not always homotopic to the original object, and it is not always thin (See Figure 4.2 for an example). Various methods have now been developed for performing skeletonization of a discrete object.

According to Palágyi [START_REF] Palágyi | A 3d fully parallel surface-thinning algorithm[END_REF], discrete skeletons can be computed using four types of methods:

Voronoi-based transformations [START_REF] Brandt | Continuous skeleton computation by voronoi diagram[END_REF][START_REF] Näf | 3d voronoi skeletons and their usage for the characterization and recognition of 3d organ shape[END_REF], distance-based transformations [START_REF] Borgefors | Computing skeletons in three dimensions[END_REF][START_REF] Toriwaki | Distance transformation and skeletonization of 3d pictures and their applications to medical images[END_REF], general-field methods [START_REF] Ahuja | Shape representation using a generalized potential field model[END_REF][START_REF] Rumpf | A continuous skeletonization method based on level sets[END_REF] and thinning. In this thesis, we will focus only on thinning methods, as they are the most studied and well defined.

Thinning Theory

Thinning consists of iteratively removing points of the object without changing its topology.

These points have specific characteristics and are called simple points. In this section, we will first provide usual definitions and notations of digital topology. Using these definitions, we will in a second time describe the characteristics of simple points.

Definitions and Notations

First we have to introduce some fundamental definitions and notations of digital topology.

In digital topology, the framework is the discrete grid Z n . In our case, we will only consider the cases n = 2 (bi-dimensional) and n = 3 (tri-dimensional). The object is represented by X ⊂ Z n , and its complementary Z n \ X is denoted by X.

A point p ∈ Z n is defined by (p 1 , ..., p n ), with p i ∈ Z.

Neighborhood

The notion of neighborhood is central for digital topology. In the 2D case, two neighborhoods (see Figure 4.3) are considered:

• the 4-neighborhood of p is the set N 4 (p) = {q ∈ Z 2 ; |q 1 -p 1 | + |q 2 -p 2 | ≤ 1}.
• the 8-neighborhood of p is the set N 8 (p) = {q ∈ Z 2 ; max(|q In 3D case, three neighborhoods (see Figure 4.4) are considered:

• the 6-neighborhood of p is the set N 6 (p) = {q ∈ Z 3 ; |q 1 -p 1 | + |q 2 -p 2 | + |q 3 -p 3 | ≤ 1}.
• the 26-neighborhood of p is the set

N 26 (p) = {q ∈ Z 3 ; max(|q 1 -p 1 |, |q 2 -p 2 |, |q 3 -p 3 |) ≤ 1}.
• the 18-neighborhood of p is the set For a k-neighborhood, we define N * k (p) = N k (p) \ {p}.

N 18 (p) = {q ∈ N 26 (p); |q 1 -p 1 |+|q 2 -p 2 |+|q 3 -p 3 | ≤ 2}.

p p p

Connectivity

From the definition of neighborhood we can propose the definition of connectivity.

Let p be an element of X, we define the k-connected component of X containing p, denoted by C k (p, X), as the maximal subset of X containing p, such that for all points a ∈ C k (p, X), there exists a sequence of points of p 0 , ..., p n such:

• p 0 = p and p n = a,

• ∀i ∈ {0, ..., n}, p i ∈ C k (p, X),

• ∀i ∈ {1, ..., n}, p i-1 is in the k-neighborhood of p i .

The set of all k-connected components of X is denoted by C k (X). We can notice that the union of all the elements of C k (X) is equal to X, and the intersection of two different elements of

C k (X) is always the empty set. A subset Y of Z n is k-adjacent to a point p ∈ Z n if Y ∩ N * k (p) = ∅.
The set of all k-connected components of X which are k-adjacent to a point p is denoted by C p k (X).

When we consider the k-connectivity of X, it is mandatory to consider a different k-connectivity for X [START_REF] Kong | Digital topology: introduction and survey[END_REF]. In bi-dimensional space, if k = 4, then k = 8, and inversely. In tri-dimensional space, if k = 6, then k = 26 and inversely. This rule is necessary in order to retrieve some important topological properties such as the Jordan theorem.

Connectivity Numbers

In order to provide local description of simple points and other characteristic points, we have to introduce the notion of connectivity numbers.

In the 2D case, let X ⊆ Z 2 and p ∈ Z 2 . For k ∈ {4, 8} the connectivity number T k (p, X) is defined by:

T k (p, X) = |C p k (N * 8 (p) ∩ X)|
In the 3D case, the definition of connectivity numbers lies on the notion of geodesic neighborhood.

Let X ⊆ Z 3 and p ∈ Z 3 . The t-order k-geodesic neighborhood of p in X is the set N t k (p, X) recursively defined by:

• N 1 k (p, X) = N * k (p) ∩ X • N t k (p, X) = {N k (q) ∩ N * 26 (p) ∩ X, q ∈ N t-1 k (p, X)}
The geodesic neighborhoods G k (p, X) are defined by: G

6 (p, X) = N 2 6 (p, X) and G 26 (p, X) = N 1 26 (p, X).
We can now define the connectivity numbers in 3D, for k ∈ {6, 26} as:

T k (p, X) = |C k (G k (p, X))| .

Simple points

Simple Points in 2D

Intuitively, a point is simple if it can be removed from an object without changing its topology.

In the digital topology framework, the topology of an object depends on the chosen connectivity;

for this reason, when considering a k-connected object, we will talk about k-simple points. The notion of simple point is central for homotopic thinning in the digital framework: a skeleton is obtained by iteratively removing simple points from an object.

According to [START_REF] Lam | Thinning methodologies-a comprehensive survey[END_REF], in the 60s, 2D simple points were characterized based on connectivity: a point p is k-simple for an object X if the removal of p does not change the number of k-connected components of X nor the number of k-connected components of X [START_REF] Brain | Graphical data processing research study and experimental investigation[END_REF][START_REF] Golay | Hexagonal parallel pattern transformations[END_REF]. This definition does not lead to efficient algorithms : indeed, in order to test if a single point is simple, it requires to scan the whole object in order to enumerate its connected components. Fortunately, local characterization of deletable points in 2D began to appear in the mid 60s [START_REF] Hilditch | Linear skeletons from square cupboards[END_REF][START_REF] Rosenfeld | Connectivity in digital pictures[END_REF][START_REF] Rutovitz | Pattern recognition[END_REF][START_REF] Yokoi | An analysis of topological properties of digitized binary pictures using local features[END_REF].

All these works established that, in order to decide whether a point is deletable or not, it is only necessary to look at the configuration of the point's neighbourhood (no need to count the number of connected components of the whole object). Consequently, in 2D, deciding if a point is simple can be done in constant time.

Proposition 1. Let X ⊂ Z 2 , and p ∈ X. If T k (x, X) = 1 and Tk(x, X) = 1, then p is k-simple for X.

Simple Points in 3D

In 3D, the removal of a point may not only change the number of connected components of the object, but may also change the tunnels of the object (see Figure 4.5 for an example). As in the 2D case, 3D simple points can be locally characterized [START_REF] Morgenthaler | Three-dimensional simple points: serial erosion, parallel thinning and skeletonization[END_REF]. Further work on 3D simple points have established that only connectivity of X and X is sufficient in order to characterize 3D simple points [START_REF] Bertrand | A boolean characterization of three-dimensional simple points[END_REF][START_REF] Bertrand | A new characterization of three-dimensional simple points[END_REF][START_REF] Malandain | Fast characterization of 3D simple points[END_REF][START_REF] Saha | Detection of 3-d simple points for topology preserving transformations with application to thinning[END_REF][START_REF] Saha | Topology preservation in 3d digital space[END_REF]. As in 2D, deciding if a point is simple can be done in constant time in 3D. Bertrand and Malandain propose a local definition of simple points using connectivity numbers T 6 and T 26 :

Proposition 2. [19] Let X ⊆ Z 3 and x ∈ X. The point x is k-simple for X iff T k (x, X) = 1
and Tk(x, X) = 1.

(a) (b) (c)

Simple Points in Higher Dimensions

Studies of simple points in 4-dimensions have also been achieved, leading once more to a local characterization of such points [START_REF] Couprie | New characterizations of simple points in 2d, 3d, and 4d discrete spaces[END_REF][START_REF] Kong | Topology-preserving deletion of 1's from 2-, 3-and 4-dimensional binary images[END_REF]. Thanks to these works, characterization of simple points in 4D can be done again in constant time.

Different Kind of Skeletons

In this section, we will review the different kinds of results which can be expected from a thinning process. The design and the constraints of the thinning algorithm will obviously depend on the desired result.

Ultimate Skeleton

Let X ∈ Z n and Y ∈ Z n be two objects. We say that X and Y are homotopic if we can transform X in Y by iteratively adding and removing simple points. We say that Y is lower homotopic to X if Y can be obtained from X by iteratively removing simple points.

Let X ∈ Z n be an object. A subset Y ⊆ X is an ultimate skeleton of X if:

• Y is lower homotopic to X

• for any point p of Y , p is not simple for Y .

Ultimate skeleton can be used to study the topology of the object, but possibly important geometric features informations are lost (see Fig.

Curvilinear and Surface Skeletons

Ultimate skeleton, providing only topological information about the shape, is not efficient as a shape descriptor. Some geometric information, like salient parts, have to be kept in the result of the thinning algorithm.

In 2D case, the skeleton has to be 1D or 0D, according to the Hilditch's criteria. A 1-dimensional skeleton, regardless of the dimension of the original object, is called a curvilinear skeleton (see In the 3D case, always according to the Hilditch's criteria, the skeleton can be 0D, 1D or 2D. A 2-dimensional skeleton is called a surface skeleton(see Fig. 4.6(e)).

These skeletons provide both topological and geometrical information about the initial shape, which summarize the general form of the object.

The choice of the skeleton dimension depends on the application. We can notice that a surface skeleton corresponds to the analogy of grass fire in 3D, and provides a more accurate descriptor of the initial object. On the other hand, the curvilinear skeleton is smaller than the surface skeleton, and the provided geometric information can be enough for some applications.

In order to obtain curvilinear or surface skeletons instead of ultimate skeleton, several approaches exist: a first one consists of forcing the preservation of medial axis points, another one consists of preserving some locally characterized points, as extremity points or isthmuses.

Extremity Points

Extremity points are simple points which have to be kept in the skeleton in order to preserve elongated geometric features.

In the case of curvilinear skeletons, the considered points are those which end a curve, and are called curve extremities. These points can be locally described:

• Let X ⊆ Z n , p ∈ X, p is a k-curve extremity iff |N * k ∩ X| = 1.
In the case of surface skeletons, the considered points are those that can be found on "surface borders". To our knowledge, there is no consensus on the definition of surface border in digital topology.

In addition to this problem, the other one is the sensitivity of this approach to small variations.

It results in skeletons containing more branches than needed to represent the geometric features.

These additional uninteresting branches are called spurious branches.

Isthmus

In order to solve the problems mentioned above, another strategy has been recently developed by Bertrand and Couprie [START_REF] Bertrand | Transformations topologiques discrètes[END_REF]. Instead of searching for simple points which have to be kept in order to preserve geometric features, we can search for non simple points having certain topological characteristics. These points are called isthmuses.

For a curvilinear skeleton, the considered isthmuses are the 1Disthmuses (see Figure 4.7 for some examples). In the case of surface skeleton, the considered isthmuses are the 2D-isthmuses.

The main advantage in regard of surface borders is that 2D-isthmus are locally well defined.

In 2D, let X ⊆ Z 2 , p ∈ X, p is 1D-isthmus iff T k (p, X) ≥ 2. In 3D, let X ⊆ Z 3 , p ∈ X, p is 1D-isthmus iff T k (p, X) ≥ 2. The point p is 2D-isthmus iff Tk(p, X) ≥ 2.

Symmetric versus Asymmetric Skeletons

Another distinction has to be done between the different kind of skeletons: symmetric and asymmetric ones. This notion of symmetry is due to interpretation in the digital framework of two of the four properties of the skeleton proposed by Blum: the centering and the thinness properties.

Assume that we have an efficient thinning process which thins an object layer by layer. It can lead to sets of simple points which cannot be all removed together without breaking the topology.

We have then two choices: we can keep all these points or we can arbitrarily remove a subset of them in order to have a thin skeleton.

The first choice leads to a symmetric skeleton: it is well centered but it is not thin. A symmetric skeleton is closer to the medial axis as it preserves more maximal ball centers.

The second choice leads to an asymmetric skeleton: it is thin but not well centered, due to arbitrary choices.

Thinning Algorithms

Homotopic thinning in the digital framework consists of removing simple points from an object, until either no more simple point can be found (resulting in an ultimate skeleton), or a satisfactory subset of voxels has been reached (resulting in a surface or curvilinear skeleton).

Two main strategies are possible for removing simple points: sequential removal and parallel removal.

Sequential algorithms

Sequential removal of simple points can be achieved by detecting simple points in an object, and removing them one after the other, until no more simple point can be found. After removing a simple point, the new set of simple points of the object must be computed. Such basic strategy does not guarantee the result, which is an ultimate skeleton, to be centered in the original object.

It is important, when designing a sequential thinning algorithm, to decide of a removal order of simple points, and of a strategy for preserving interesting visual features of the object (as explained in Section 4.2).

In order to obtain a centered skeleton, one must define a precise order of removal of simple points. Usually, in order to get a centered skeleton, it is necessary to delete simple points "layer by layer", from the outer layer to the inner one. Many strategies have been proposed in 2D for choosing of a removal order (a very exhaustive survey of thinning methods in 2D before 1992 can be found in [START_REF] Lam | Thinning methodologies-a comprehensive survey[END_REF]) : for example, in the 80s, it was proposed to follow an object's contour in order to find and remove simple points "layer by layer" [START_REF] Arcelli | Pattern thinning by contour tracing[END_REF][START_REF] Pavlidis | A thinning algorithm for discrete binary images[END_REF]. However, this thinning scheme (as many others) hardly generalizes to 3D.

A widely used strategy to obtain a centered skeleton with a sequential thinning process consists of computing a priority function on the object and removing the simple points of X according to the value of this function [START_REF] Davies | Thinning algorithms: A critique and a new methodology[END_REF] : at each step, the simple point that is removed is one with the lowest possible value. The Euclidean distance map (storing distances to the complementary) is widely used as priority function [START_REF] Couprie | Discrete bisector function and euclidean skeleton in 2d and 3d[END_REF][START_REF] Malandain | Euclidean skeletons[END_REF][START_REF] Talbot | Euclidean skeletons and conditional bisectors[END_REF]. Other works use discrete distances, such as the chamfer distance [START_REF] Pudney | Distance-ordered homotopic thinning: A skeletonization algorithm for 3d digital images[END_REF], to decide of a removal order. inhibitor set, is a set of points of the input object which must be in the resulting skeleton, and D is the priority function used to decide an order of points removal (here, the lower priority means faster removal, so it is possible to use a distance map as priority function).

An inhibitor set allows to choose "anchor points" for the skeleton, and therefore preserve the visual aspect of the original object in the skeleton. However, when performing a thinning guided by an Euclidean distance map, the points of the inhibitor set and the directions of thinning followed by the algorithm are not always "compatible". In [START_REF] Talbot | Euclidean skeletons and conditional bisectors[END_REF], the authors use a thinning algorithm where the slope of the priority function is used to dynamically add points to the constraint set. In [START_REF] Couprie | Discrete bisector function and euclidean skeleton in 2d and 3d[END_REF], the authors propose to merge the slope calculation into the priority The main problem of parallel thinning is that removing simple points simultaneously from an object usually "breaks" the topology. Thus, additional conditions must be introduced in order to solve this problem. In order to describe the different strategies, we introduce some vocabulary found in [START_REF] Palágyi | A 3d fully parallel surface-thinning algorithm[END_REF].

A thinning operator transforms an object X only by removing some points, without changing the topology of the shape. A parallel thinning operator simultaneously deletes all points of the object X satisfying some conditions (called the deletion conditions). The support of a thinning operator is the minimal set of points whose values must be examined to verify the deletion conditions for a point p ∈ X.

According to Hall [START_REF] Hall | Parallel connectivity-preserving thinning algorithms[END_REF], parallel thinning algorithms can be divided into three categories:

• In directional algorithms, the main loop is divided into sub-iterations, and the deletion conditions are changed from one sub-iteration to another.

• In subfield-based algorithms, the points of the object are decomposed into subsets, and at a given iteration of the algorithm, only points in a given subset are studied.

• In fully parallel algorithms, no sub-iteration takes place : the same thinning operator is used at each iteration of the main loop. However, if for the other strategies the usual support is a 3 × 3 × 3 neighborhood (i.e. N 26 ), fully parallel operators require an extended support (usually, 5 × 5 × 5).

Directional Algorithms

A directional algorithm is a parallel algorithm where each iteration is divided in sub-iterations using different thinning operators. Each thinning operator can only delete points belonging to a specific border orientation. This border orientation classification is performed regarding the belonging to complementary of points in a given neighborhood.

The first directional thinning algorithm has been proposed in 1971 by Stefanelli and Rosenfeld [START_REF] Stefanelli | Some parallel thinning algorithms for digital pictures[END_REF] in the Z 2 framework. The authors used the 4-neighborhood to classify the borders points (see Figure 4.8), leading to a 4 sub-iteration algorithm. In the Z 3 framework, the most usual kind of directional algorithms use a classification in regard of the 6-neighborhood, leading to algorithms with 6 sub-iterations [13,[START_REF] Gong | A simple parallel 3D thinning algorithm[END_REF][START_REF] Lee | Building skeleton models via 3-D medial surface/axis thinning algorithms[END_REF][START_REF] Lohou | A 3d 6-subiteration curve thinning algorithm based on p-simple points[END_REF][START_REF] Ma | Parallel thinning algorithms on 3D (18, 6) binary images[END_REF][START_REF] Mukherjee | On connectivity issues of ESPTA[END_REF][START_REF] Palágyi | A 3D 6-subiteration thinning algorithm for extracting medial lines[END_REF][START_REF] Xie | A topology-preserving parallel 3D thinning algorithm for extracting the curve skeleton[END_REF] or 3 sub-iterations [START_REF] Palágyi | A 3-subiteration 3D thinning algorithm for extracting medial surfaces[END_REF][START_REF] Palágyi | A 3-Subiteration Surface-Thinning Algorithm[END_REF][START_REF] Palágyi | A subiteration-based surface-thinning algorithm with a period of three[END_REF], considering in one sub-iteration a class and its opposite.

Directional algorithms with 8 sub-iterations [START_REF] Palágyi | Directional 3D thinning using 8 subiterations[END_REF] (classification in regard of N 26 \ N 18 neighborhood) and 12 sub-iterations [START_REF] Lohou | A 3d 12-subiteration thinning algorithm based on p-simple points[END_REF][START_REF] Palágyi | A parallel 3D 12-subiteration thinning algorithm[END_REF] (classification in regard of N 18 \ N 6 neighborhood) have also been proposed.

Generally, a directional algorithm is designed on the model of Algorithm 3, using classification of the directions in the order proposed in Figure 4.9 for 3D, or Figure 4.8 for 2D.

Algorithm 3: Directional Thinning(X) Data: A connected shape X Result: A skeleton of X repeat foreach direction i do S = {x ∈ X|x satisfies deletion conditions of the ith thinning operator}

X = X\S until stability ; return X 4.3.2.2 Subfield-Based Algorithms In subfield-based algorithms, Z n is divided in m ≥ 2 disjoint subfields Z n 0 , ..., Z n m-1 , such m-1 i=0 Z n i = Z n .
An iteration of subfield-based algorithm is divided in sub iterations corresponding to the activation of each subfield. A point can be deleted only when the subfield it belongs to is activated.

The first subfield-based algorithm, at our knowledge, has been proposed in 1984 by Groen and

Foster [START_REF] Groen | A fast algorithm for cellular logic operations on sequential machines[END_REF] in the Z 2 framework.

In the Z 3 , algorithms using 2 subfields [START_REF] Kardos | Topology preserving 2-subfield 3d thinning algorithms[END_REF][START_REF] Ma | A medial-surface oriented 3-d two-subfield thinning algorithm[END_REF][START_REF] Ma | Extracting medial curves on 3D images[END_REF], 4 subfields [START_REF] Ma | Three-dimensional topology preserving reduction on the 4subfields[END_REF][START_REF] Németh | Topology Preserving 3D Thinning Algorithms Using Four and Eight Subfields[END_REF] and 8 subfields [START_REF] Németh | Topology Preserving 3D Thinning Algorithms Using Four and Eight Subfields[END_REF] have been proposed.

The subfields are designed in sort that if 2 points p, q belong to the same subfield, q is not in a given neighborhood of p (see Figure 4.10):

• for m = 2, ∀0 ≤ i < m, ∀(p, q) ∈ Z 3 i × Z 3 i , q / ∈ N 6 (p) • for m = 4, ∀0 ≤ i < m, ∀(p, q) ∈ Z 3 i × Z 3 i , q / ∈ N 18 (p) • for m = 8, ∀0 ≤ i < m, ∀(p, q) ∈ Z 3 i × Z 3 i , q / ∈ N 26 (p)
It is this property which allow the preservation of the topology, the deletion of a point having a controlled impact on its neighborhood.

Generally, a subfield-based algorithm is designed on the model of Algorithm 4, using subfields proposed in Figure 4.10.

Notice that this strategy is known to usually provide bad quality skeletons. Algorithm 4: Subfield Based Thinning(X) Data:

A connected shape X Result: A skeleton of X repeat foreach subfield F i do S = {x ∈ X ∩ F i |x satisfies deletion conditions } X = X\S until stability ; return X 4.3.2.

Fully Parallel Algorithms

Fully parallel algorithms use only one thinning operator in the full space of the same iteration.

In order to do it without breaking the topology, the support has to be bigger than for other classes of thinning algorithms.

Rutovitz was the first to propose a (fully) parallel thinning algorithm, in 1966 ( [START_REF] Rutovitz | Pattern recognition[END_REF]). However, it has been proved that Rutovitz's algorithm does not always preserve topology and "patches" exist in order to correct it ( [START_REF] Couprie | Note on fifteen 2d parallel thinning algorithms[END_REF]). It is in 1981 that Pavlidis published the first fully parallel homotopic thinning algorithm ( [START_REF] Pavlidis | A flexible parallel thinning algorithm[END_REF]), in 2D for 8-connected objects, that was later proved to preserve topology ( [START_REF] Couprie | Note on fifteen 2d parallel thinning algorithms[END_REF]).

In the 3D digital topology framework, the problem is very complex, and several theories have been developed in order to help to design correct fully parallel thinning algorithms.

Minimal non-deletable sets Minimal non-deletable sets were introduced by C. Ronse ([146]) in order to characterize under which conditions simple points could not be removed simultaneously from a 2D object without "breaking" the topology.

Minimal non-deletable sets define "forbidden features" that should not appear in a set of pixels in order for it to be simple. Minimal non-deletable sets were designed in order to prove that 2D parallel thinning algorithms were topology preserving (and therefore valid) by testing only a small number of configurations of points. A computer-based implementation of these tests was later proposed in [START_REF] Hall | Tests for connectivity preservation for parallel reduction operators[END_REF], a 3D implementation of these sets was proposed in [START_REF] Kong | Problem of determining whether a parallel reduction operator for ndimensional binary images always preserves topology[END_REF], [START_REF] Ma | On topology preservation in 3d thinning[END_REF] and [START_REF] Kong | On topology preservation in 2-d and 3-d thinning[END_REF], and a 4D implementation was proposed in [START_REF] Gau | Minimal non-simple sets in 4d binary images[END_REF] and [START_REF] Kong | Minimal non-simple sets in 4-dimensional binary images with (8,80)-adjacency[END_REF].

P -simple points In 1995, Bertrand introduced the P -simple points in order to characterize, in 3D, which simple points could be removed simultaneously ( [13]).

As with the minimal non-deletable sets, the P -simple points allow to check if existing parallel 3D thinning algorithms work: indeed, in [START_REF] Bertrand | Sufficient conditions for 3D parallel thinning algorithms[END_REF], the author gives a method for checking, based on the P -simple points framework, the topological validity of thinning algorithms. Moreover, P -simple points were widely used in order to propose new parallel 3D thinning algorithms (an example of a new algorithm is given in [START_REF] Bertrand | Sufficient conditions for 3D parallel thinning algorithms[END_REF]).

Critical kernels

The recent critical kernels, introduced in [START_REF] Bertrand | On critical kernels[END_REF], present a new framework for performing parallel thinning in 2D, 3D and 4D (see [START_REF] Bertrand | A new 3d parallel thinning scheme based on critical kernels[END_REF] and [11]). In relation to critical kernels, a new definition and new characterizations of simple points in 2D, 3D and 4D (see [START_REF] Couprie | New characterizations of simple points in 2d, 3d, and 4d discrete spaces[END_REF]) have been proposed, and links between this framework, P -simple points and minimal non-deletable sets were established in [START_REF] Couprie | New characterizations of simple points, minimal non-simple sets and p-simple points in 2d, 3d and 4d discrete spaces[END_REF]. More precisely, critical kernel theory is a generalization of the two other theories. Critical kernels were also used to prove that some thinning algorithms were valid, while others were not correct ( [START_REF] Couprie | Note on fifteen 2d parallel thinning algorithms[END_REF]).

Chapter 5

Faster Implementation of Thinning

Schemes

In our context of real time motion capture, the considered 3d objects contain between 50 000 and 1 000 000 voxels, and have to be skeletonized in less than 40 milliseconds (remember that skeletonization is only a part of the process which has to be "real time"). No direct implementation of existing thinning algorithms running on an average computer (with CPU frequency lower than 4GHz) can reach this time for these sizes of object. In this chapter we will first discuss speed optimizations schemes found in the literature, then propose some new ones.

Classical Optimizations

In the literature, the thinning algorithms are most usually described using pseudo code, only describing the main lines of the real implementation, and skipping details about, e.g. how to reach low execution time. However, thinning algorithms used in full applications (like medical imaging or material studies) are implemented using some common optimizations which can be called "classical". In this section, we will review the most usual of them.

Restriction of Points to Check

The first and the most common optimization method consists of a restriction of points to test for deletion conditions at each iteration. If we implement directly the algorithms described in pseudo code in the majority of the literature, we have to check all the points in the image grid at each iteration, leading to a time complexity in O(|G| × I), with I being the number of iterations and |G| being the size of the image grid.

37

First, as it is a deletion condition test, we can test only the points which can be deleted (i.e. the points belonging to the object). Furthermore, in order to be deletable in a thinning process, a point has to be simple. Considering this fact, we can test only the points at the border of the object at each iteration. This is naturally done by sequential thinning using distance transform as a priority function [START_REF] Couprie | Discrete bisector function and euclidean skeleton in 2d and 3d[END_REF][START_REF] Malandain | Euclidean skeletons[END_REF][START_REF] Pudney | Distance-ordered homotopic thinning: A skeletonization algorithm for 3d digital images[END_REF][START_REF] Talbot | Euclidean skeletons and conditional bisectors[END_REF], by sorting the points in regard of their distance to the border.

In the case of parallel algorithms, a list can be used in order to store the border points [START_REF] Palágyi | A 3d fully parallel surface-thinning algorithm[END_REF][START_REF] Palágyi | Fully parallel 3D thinning algorithms based on sufficient conditions for topology preservation[END_REF].

Before the iterative part of the algorithm, the list is initialized by adding all the border points of the image. During the iterative process, only points of the list are tested: if a point has to be deleted, we remove it from the list and we add its neighbors belonging to the object which are not already in the list.

Using this optimization, the time complexity is linear with regard to the size of the image.

Optimization of Deletion Conditions Tests

Depending on the considered algorithm, the tests of deletion conditions can be very lengthy.

Assuming that deletion tests take into consideration only a determined finite set of neighbors (the support of the thinning operator), pre-computation methods can be used in order to obtain quicker and constant time.

Usage of Look-Up Tables

In thinning algorithms, a widely used tool(explicitly in [START_REF] Palágyi | A 3-subiteration 3D thinning algorithm for extracting medial surfaces[END_REF][START_REF] Palágyi | A 3-Subiteration Surface-Thinning Algorithm[END_REF][START_REF] Palágyi | A subiteration-based surface-thinning algorithm with a period of three[END_REF][START_REF] Palágyi | A 3d fully parallel surface-thinning algorithm[END_REF], for example) to reduce computation time is the look-up table (usually shortened in LUT ).

A look-up table is an array containing the precomputed results of some time-consuming problem P for a given set of input values V . In order to find the index in the look-up table storing the result P (v), v ∈ V , we have to define a bijective function K : V → {0, ..., S -1}, with S being the size of the look-up table.

In the case of thinning algorithms, the problem P is the deletion condition test returning a boolean (true if the point can be deleted, false otherwise) and the set of input values is the set of all possible configurations of points in the support.

For a 3 × 3 × 3 support, we have to consider that each point in the 26 neighborhood may or may not belong to the object, leading to a look-up table size of 2 26 = 67 108 864 values. Since each result is a boolean, the look-up table uses 4 megabytes in memory.

The bijective function K is defined using some properties of the possible configurations of points in the support. Assuming that each point in the support has a binary status (1 if the point belongs to the object, 0 otherwise), each configuration can be represented by a vector v = (v 0 , ..., v 25 ) ∈ {0, 1} 26 , were each v i represents the status of the ith point for a given navigation order of the support. We define the bijective function providing an index to each configuration by:

K(v) = 26 i=0 2 i * v i
The navigation order of the support depends on the algorithm: if the algorithm needs only one look-up table, the order will be the same as the order of the image representation in memory, in order to minimize transfer between RAM and CPU cache. If the algorithm needs several lookup tables based on different supports for the same point, the order will place first the points of the 6-neighborhood, then those of the 18-neighborhood and finally those of 26-neighborhood.

Usage of Binary Decision Diagrams

Look up tables are not the only way to store and access pre-computed results of a function. For functions taking as input a vector of booleans, as in the case of deletion condition tests, specific data structures have been designed: the binary decision trees, and their compact representation, the binary decision diagrams [START_REF] Bryant | Graph-Based Algorithms for Boolean Function Manipulation[END_REF] (usually shortened in BDD). The use of binary decision diagrams is proposed in some thinning literature papers [START_REF] Lohou | A 3d 12-subiteration thinning algorithm based on p-simple points[END_REF][START_REF] Lohou | A 3d 6-subiteration curve thinning algorithm based on p-simple points[END_REF].

A binary decision tree is a binary tree with several specificities:

• each leaf represents a result;

• if the input vector contains n values, the tree has n levels of interior nodes;

• all the interior nodes of the ith level of the tree represent the ith value of the input vector;

• the leaves of the left (respectively, right) subtree of an interior node at level i contains results obtained if the ith value equals 0 (respectively 1).

Finding the result for a given boolean vector is done using Algorithm 5. A binary decision diagram is obtained from a binary decision tree by applying two operations:

0 b 1 b 1 b 2 b 2 b 2 b 2 b 3 b 3 b 3 b 3 b 3 b 3 b 3 b 3 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 level0 level1 level2 level3
1. all subtrees which are equal are merged;

2. remove nodes with two subtrees which are equal.

Look Up Tables versus Binary Decision Diagrams

A good point for the binary decision diagram is the potentially lower space needed to store the results. Concerning the speed-up, contrary to the LUT which provides an access in constant time, access with BDD depends both on the considered deletion conditions and the order in which the values of the configuration are considered. Furthermore, even if the number of points to consider is reduced, navigation in the binary decision diagram can be more time consuming than the direct access to an array as required by LUT methods.

Look Up Tables for Critical Kernel Based Algorithms

Look up tables can not be directly used for algorithms using extended 5 × 5 × 5 support: it implies more than 2 × 10 37 values, and storing a such look-up table is not possible, even using all the hard drives existing in the world. Furthermore, the computation of values of the whole look-up table with the fastest super-computer would take more time than the age of the universe.

In 2009, Palágyi and Németh [START_REF] Palágyi | Fully parallel 3D thinning algorithms based on sufficient conditions for topology preservation[END_REF] proposed fully parallel algorithms using a thinning operator with 5 × 5 × 5 support, the deletion conditions taking into account not only the configuration of the 26-neighborhood, but also the simplicity of the points in this neighborhood. However, the authors proposed algorithms using tests only in the 3×3×3 neighborhoods, using a decomposition of the tests, and allowing for the usage of look-up tables.

In 2006, Bertrand and Couprie [START_REF] Bertrand | A new 3d parallel thinning scheme based on critical kernels[END_REF] proposed a new fully parallel thinning scheme based on critical kernels, allowing to easily design thinning algorithms for any kind of skeletons (homotopic kernel, curvilinear or surfacic, asymmetric as well as symmetric). Even if the support of the thinning operator is 5 × 5 × 5, we will propose in this section a way to decompose the tests as evoked above and to use look-up tables for critical kernel based algorithms, in order to improve their speed.

Critical Kernel Based Algorithms

First we briefly recall the scheme of critical kernel based algorithms. These algorithms use the notion of crucial cliques, which are set of points which can not be all removed in order to preserve topology. The crucial points are classified in several subtypes : 3-crucial, 2-crucial, 1-crucial and 0-crucial. An iteration of the algorithm consist in keeping only points of the shape which are crucial.

The belonging to a crucial clique is detected using a set of templates (see Figure 5.4). Definition 3.

A B C 1 C 2 C 3 C 4 C 5 C 6 C 7 C 0 D 1 D 0 D 2 D 3 D 4 D 5 D 6 D 7 M 2 mask R Q 0 Q 1 Q 2 Q 3 Q 4 Q 5 Q 6 Q 7 A B C D E F G H J K L I A B C D E F G H M ′ 2 mask M 1 mask M 0 mask
[17] Let X ⊆ Z 3 the object, and let M be a set of voxels of X.

1) The set M matches the mask M 2 if the three following conditions hold: i) the voxels in M are simple for X.

ii

) M = {A, B}. iii) the 2D configuration M ′ 2 obtained by setting R ∈ M ′ 2 and setting Q i ∈ M ′ 2 if {C i , D i } ∩ X = ∅, with i ∈ [0, ..., 7], is such that R is non-simple in the 2D sense.
2) The set M matches the mask M 1 if the four following conditions hold: i) the voxels in M are simple and not 2-crucial for X.

ii) M = {E, F, G, H} ∩ X. iii) the set {E, G} or the set {F, H} (or both) is included in M . iv) we have either [U ∩ X = ∅ and V ∩ X = ∅] or [U ∩ X = ∅ and V ∩ X = ∅], with U = {A, B, C, D} and V = {I, J, K, L}.
3) The set M matches the mask M 0 if the three following conditions hold: i) the voxels in M are simple and neither 2-crucial nor 1-crucial for X.

ii) M = {A, B, C, D, E, F, G, H}.

iii) at least one of the sets {A, G}, {B, H}, {C, E}, {D, F } is a subset of M . Definition 4. [START_REF] Bertrand | A new 3d parallel thinning scheme based on critical kernels[END_REF] Let X ⊆ Z 3 and let M be a set of voxels of X. i) M is a 2-crucial clique for X if and only if M matches the mask M 2 or a mask obtained from it by π 2 rotations.

ii) M is a 1-crucial clique for X if and only if M matches the mask M 1 or a mask obtained from it by π 2 rotations.

iii) M is a 0-crucial clique for X if and only if M matches the mask M 0 .

Using these definitions, Bertrand and Couprie proposed several generic fully parallel thinning algorithms, like Algorithm 6 for symmetric skeletons and Algorithm 7 for asymmetric ones.

Note that these algorithms ensure the preservation of some points defined by the set K and the use of any kind of extremity point or isthmus, defined by the function Ψ.

Function Ψ will return the subset of the input set S ∈ Z 3 which has to be preserved in the skeleton (for an homotopic kernel, Ψ always returns ∅).

Algorithm 6: CK Thinning Data:

X ⊆ Z 3 , K ⊆ X, Ψ : X → Ψ(X) Result: A skeleton of X repeat K = K ∪ Ψ(X)
R = set of non-simple points for X or which are in K foreach possible position of mask M 2 (and any of its

π 2 rotations) do if M 2 matches then R = R ∪ {A, B}
foreach possible position of mask M 1 (and any of its

π 2 rotations) do if M 1 matches then R = R ∪ [{E, F, G, H} ∩ X] foreach possible position of mask M 0 do if M 0 matches then R = R ∪ [{A, B, C, D, E, F, G, H} ∩ X] X = R until stability ; return X Algorithm 7: ACK Thinning Data: X ⊆ Z 3 , K ⊆ X, Ψ : X → Ψ(X) Result: A skeleton of X repeat K = K ∪ Ψ(X)
R = set of non-simple points for X or which are in K foreach possible position of mask M 2 (and any of its

π 2 rotations) do if M 2 matches and {A, B} ∩ R = ∅ then R = R ∪ {A}
foreach possible position of mask M 1 (and any of its

π 2 rotations) do if M 1 matches and {E, F, G, H} ∩ R = ∅ then P = first voxel in {E, F, G, H} ∩ X R = R ∪ {P } foreach possible position of mask M 0 do if M 0 matches and {A, B, C, D, E, F, G, H} ∩ R = ∅ then P = first voxel in {A, B, C, D, E, F, G, H} ∩ X R = R ∪ {P } X = R until stability ; return X

Reformulation of templates

In order to be compatible with LUT usage, a mask has to take into consideration configurations of points with binary state.

Some observations can be done on the used templates:

• In all cases two different informations are taken into consideration:

1. The configuration of the points belonging to the object (we denote this configuration by C X ).

2. The configuration of the points belonging to the object and holding the same set of additional properties, e.g. being non-simple (we denote this configuration by C S ).

Notice that these two configurations are contained in a 3 × 3 × 3 support.

• Due to both symmetries of the templates and their use in the scheme, it is possible to consider only deletion conditions on one voxel, instead of a set, the other voxels of the set being tested for the same neighborhood at a different step of the same iteration.

Using these two observations, we can redefine the masks M i , i ∈ {0, 1, 2}, using for each of them two templates: one providing conditions for C X and the other one for C S . This separation in two masks allows for the use of look-up tables.

In the sequel of this section we say that the configuration C S is compatible with the mask M iS , i ∈ {0, 1, 2}, if and only if the points represented by black disks in the mask belong to the set S and the points represented by white disks in the mask do not belong to the set S. A detection of 2-crucial points can be done using masks shown in Figure 5.5.

5.2.2.1 Reformulation of M 2 P C 1 C 2 C 3 C 4 C 5 C 6 C 7 C 0 D 1 D 0 D 2 D 3 D 4 D 5 D 6 D 7 M 2X mask R Q 0 Q 1 Q 2 Q 3 Q 4 Q 5 Q 6 Q 7 M ′ 2 mask P M 2S mask
Definition 5. Let X ⊆ Z 3 . A) A point P matches the mask M 2X if the 2D configuration M ′ 2 obtained by setting R ∈ M ′ 2 and setting Q i ∈ M ′ 2 if {C i , D i } = ∅, with i ∈ [0, ..., 7], is such that R is non-simple in the 2D sense (condition on C X ).
B) A point P matches the mask M 2S if the configuration of the set of simple points is compatible with M 2S . C) A simple point P is 2-crucial if and only if it matches the masks M 2X and M 2S , or the masks obtained by a same π or π 2 rotation of these two masks.

Reformulation of M 1

A detection of 1-crucial points can be done using masks shown in Figure 5.6.

Definition 6. Let X ⊆ Z 3 .
A) A point P matches the mask M 1X if the two following conditions hold: i) the set {E, G} or the set {F } (or both) is included in X. ii) we have either

[U ∩ X = ∅ and V ∩ X = ∅] or [U ∩ X = ∅ and V ∩ X = ∅], with U = {A, B, C, D} and V = {I, J, K, L}.
B) A point P matches the mask M 1S if the configuration of the set of 2-crucial or not simple points belonging to X is compatible with M 1S .

C) A simple and non 2-crucial point P is 1-crucial if and only if it matches the masks M 1X and M 1S , or the masks obtained by the same combination of π 2 rotations of these two masks.

Reformulation of M 0

A detection of 0-crucial points can be done using masks shown in Figure 5.7.

Definition 7. Let X ⊆ Z 3 .
A) A point P matches the mask M 0X if at least one of the sets {A, G}, {B}, {C, E}, {D, F } is a subset of X.

B) A point P matches the mask M 0S if the configuration of the set of 2-crucial or 1-crucial or non simple points belonging to X is compatible with M 0S .

C) A simple and neither 2-crucial nor 1-crucial point P is 0-crucial if and only if it matches the masks M 0X and M 0S , or the masks obtained by the same combination of π 2 rotations of these two masks.

Look Up Tables for New Masks

Now we have masks taking into consideration the configuration of points sharing the same properties, so we can use look-up tables.

We can first observe that masks M 2S and M 1S do not take into consideration all the points in N 26 , but only points in N 6 and N 18 , respectively. For these masks we can use look-up tables with size reduced to 2 6 and 2 18 , respectively.

As each couple of masks has several orientations and as the properties checking requires the matching of the two masks for a same orientation, we have to store separately the result of the matching of each mask for each orientation.

Due to this fact, in each of the six look-up tables (one per mask) the result for a configuration

C will be a vector of bits v = v 0 , ..., v D-1 , such v i , i ∈ [0, ...D -1] indicate if C matches (v i = 1)
or not (v i = 0) the mask in its ith orientation.

The properties checking is then obtained as follows. Let LU T iS and LU T iX , i ∈ {0, 1, 2}, the look-up tables storing the results of matching masks M iS and M iX , respectively. Let index of C be the index of the configuration C computed as explained in Section 5. and LU T iX (index of C X (p)) have a common bit to 1 in their vector.

This verification can be done using the bitwise AN D operator: p will match M iS and M iX if

and only if LU T iS (index of C S (p)) AN D LU T iX (index of C X (p)) = 0.
To sum up, we use six look-up tables for checking of crucial points, as summarized in the

Reformulation of the Thinning Schemes

Using the look-up tables described above, we can now reformulate the thinning schemes.

Algorithm 8: Optimized ACK Thinning Data:

X ⊆ Z 3 , K ⊆ X, Ψ : X → Ψ(X) Result: A skeleton of X repeat K = K ∪ Ψ(X) S = set of simple points for X which are not in K foreach point p ∈ S do if LU T 2X (index 26 (p, X)) AND LU T 2S (index 6 (p, S)) = 0 then S = S \ {p} foreach point p ∈ S do if LU T 1X (index 26 (p, X)) AND LU T 1S (index 18 (p, X \ S)) = 0 then S = S \ {p} foreach possible position of mask M 0 do if LU T 0X (index 26 (p, X)) AND LU T 0S (index 26 (p, X \ S)) = 0 then S = S \ {p} X = X \ S until stability ; return X
We denote by index n (p, S) the index in the look-up table of the configuration of the set S in the n-neighborhood of a point p. Algorithms 9 and 8 are the LUT-optimized reformulations of the thinning Algorithms 6 and 7, respectively.

Algorithm 9: Optimized CK Thinning Data:

X ⊆ Z 3 , K ⊆ X, Ψ : X → Ψ(X) Result: A skeleton of X repeat K = K ∪ Ψ(X) S = set of simple points for X which are not in K R = ∅ foreach point p ∈ S do if LU T 2X (index 26 (p, X)) AND LU T 2S (index 6 (p, S)) = 0 then R = R ∪ {p} S = S \ R R = ∅ foreach point p ∈ S do if LU T 1X (index 26 (p, X)) AND LU T 1S (index 18 (p, X \ S)) = 0 then R = R ∪ {p} S = S \ R R = ∅ foreach possible position of mask M 0 do if LU T 0X (index 26 (p, X)) AND LU T 0S (index 26 (p, X \ S)) = 0 then R = R ∪ {p} S = S \ R X = X \ S
until stability ; return X

Speed up

The use of look-up tables on critical kernel based thinning (already using border list defined in Section 5.1.1) proposed in this section allows a speed-up yielding a division by ten of the computation time for the thinning on some objects. Please refer to the Section 5.4 for more details.

Configurations Computation Speed-Up

As explained in Section 5.1.2.1, in order to access to the results stored in a look-up table, we have first to compute the index corresponding to the configuration of the point. Computing this index is time consuming, due to the fact that values of the neighborhood are spaced in the memory, whatever the used representation of the image, involving cache misses (transfers between CPU cache and main memory).

Furthermore, for the same point, several conditions must be checked, depending on the thinning method: in directional thinning, the same point must be checked for all directions; in critical kernel based thinning, points are checked for simplicity, then for the different crucial properties.

As each check is done in a different step of the iteration, the configuration index is recomputed for each one, involving a waste of time.

We can notice that in a thinning algorithm, the configuration of a point changes only when one of its neighbors is deleted, thus the configuration changes a limited number of times. In this section, we propose a method for avoiding multiple useless computations of the configuration index of the same point.

At the beginning of the algorithm, we build a map M with the same dimensions as the image I, such that for each point p ∈ I, M (p) contains the index of the configuration of p in I.

During the iterative process, when a point p is deleted, we update the configuration index of its neighbors in M , by putting to 0 the bit representing p in their configuration.

Using this method, the configurations are computed only when a point is deleted, instead of each time a point is checked.

We call this optimization the Configuration Computation Speed-Up, abbreviated as CCSU.

Speed up

The use of this configuration computation speed-up divides by two (in average) the computation time of a thinning. Please refer to Section 5.4 for more details.

Benchmarks

In this section, we study the speed gain of the proposed optimizations.

Methodology

For this purpose, we implemented several versions of the same algorithm with different levels of optimization: a first version without the use of LUT, a second version using LUT, and a third one using both LUT and CCSU. All the versions are optimized using border lists, as explained in Section 5.1.1. Then, we applied the different algorithms on a set of various images, and we compared their computation speed.

Implemented Algorithms

For the evaluation of LUT usage on critical kernel based methods, we implemented the asymmetric scheme (see Algorithm 7), using a constraint function Ψ returning the set of 1D isthmus of the shape. We denote this algorithm by ACK3.

In order to observe the impact of CCSU (in addition to LUT) on different thinning strategies, we also implemented the 6-directional thinning algorithm proposed by Palágyi and Kuba in 1998 [START_REF] Palágyi | A 3D 6-subiteration thinning algorithm for extracting medial lines[END_REF]. We denote this algorithm by PKD6.

Tested Images

The algorithms are tested on a set of images of different sizes, containing various objects. Our image set is represented in Figure 5.8.

Results

The results are obtained by averaging times of ten thinnings for each image. All the experiments have been done on a computer with a processor Intel(R) Q8200 at 2.33 GHz and 3 gigabytes of RAM.

Figures 5.9 and 5.10 show the computation speed results of ACK3 and PKD6, respectively.

Discussion

CCSU Speed Gain

An observation of Figures 5.11 and 5.12 shows some interesting characteristics about CCSU.

First we can notice that the speed gain is higher for PKD6 than for ACK3: it is due to the larger amount of deletion condition tests for the same point in this directional algorithm than in the fully parallel algorithm ACK3.

We also observe the quasi inexistent impact of CCSU optimization on ACK3 for some shapes (9,13,[START_REF] Blum | An associative machine for dealing with the visual field and some of its biological implications[END_REF][START_REF] Borgefors | Computing skeletons in three dimensions[END_REF]. We can see that the LUT optimization is inefficient too for these shapes. We can observe that these shapes are thinner than the other ones, and thus require less iterations for the same image size.

In the case of shape 13, the dancer, we can notice that CCSU slows down PKD6. It is due to the higher time consuming pre-iteration step, which cannot be offset by the gain during the iterative process, due to the low amount of iterations.

To sum up, CCSU is an optimization that is easy to set up, and which provides a sensible speed-up in a majority of cases. However, in the case of shapes close to their skeleton, it is inefficient.

LUT Speed Gain for Critical Kernel Based Thinning

Observing Figure 5.11, we can notice pretty spectacular speed-up by using LUT on ACK3 (more than 12 times faster in the case of shape 7). However, we can notice that the effect is inexistent for some shapes (already mentioned in the previous section).

Combining both LUT and CCSU highly improve the speed of critical based algorithms, but the speed gain highly depends to the shape in input. For each frame of each data set, we compute the visual hull at different voxel grid resolutions and we measure the time to skeletonize the obtained object using our two optimized algorithms.

Application in Motion Capture Framework

The average computation time for each resolution on each data set are reported in Figure 5.13.

We can observe that despite of our efforts for speed-up critical kernel based thinning, ACK3 is still five times slower than directional thinning, and is not fast enough for its use in real time applications. However, we can notice that directional algorithms like PKD6 are very fast and reach the real time even for high resolutions of the voxel grid.

Chapter 6

Isthmus Based Directional Thinning

In our context of motion capture, we have to consider both the speed of the thinning process (see Chapter 5), and the quality of the produced skeletons. The amount of noise in the resulting skeleton has an impact on the good matching with the model. Due to the use of 1D-isthmus, ACK3 provides curvilinear skeleton with very low amount of spurious branches. However, as seen in the end of the previous chapter, this algorithm is not the fastest one.

On the other hand, faster algorithms, like directional thinning or subfield based algorithms, provide noisier skeletons. It is mainly due to the consideration of curve end points as constraint points, which are more likely to appear during thinning than isthmuses.

In this chapter, we design a directional thinning algorithm considering isthmuses as constraints.

Furthermore, we propose a constraint check at the iteration level, instead of sub-iteration level, in order to consider less constraint points than usually in this kind of strategy.

Background

Our method is based on two different works. We use constraints that were originally introduced by Bertrand and Couprie [START_REF] Bertrand | Transformations topologiques discrètes[END_REF] in the context of sequential thinning, and the design of the masks used in our parallel directional algorithm is based on those defined by Palágyi and Kuba [START_REF] Palágyi | A 3D 6-subiteration thinning algorithm for extracting medial lines[END_REF].

Bertrand and Couprie Isthmus Based Thinning

In 2007, Bertrand and Couprie [START_REF] Bertrand | Transformations topologiques discrètes[END_REF] The consideration of isthmuses instead of extremities points is interesting for two reasons: isthmuses can be easily locally defined and detected (it is not the case of surface end points), and appear less often than extremities during thinning process, leading to less noisy skeletons. However, the sequential approach proposed is less adapted than parallel ones for this purpose, and provides skeletons of lower quality. The algorithm is summarized in Algorithm 10.

Algorithm 10: PKD6 [START_REF] Palágyi | A 3D 6-subiteration thinning algorithm for extracting medial lines[END_REF] Data:

X ⊆ Z 3 Result: A skeleton of X Y = X repeat Y = Y \ {p ∈ Y ; p matching one mask of M U } Y = Y \ {p ∈ Y ; p matching one mask of M D } Y = Y \ {p ∈ Y ; p matching one mask of M N } Y = Y \ {p ∈ Y ; p matching one mask of M S } Y = Y \ {p ∈ Y ; p matching one mask of M E } Y = Y \ {p ∈ Y ; p matching one mask of M W } until stability ; return Y
Palágyi and Kuba proved that their algorithm preserves topology, using the following theorem proved by Ma in [START_REF] Ma | On topology preservation in 3d thinning[END_REF].

Definition 8. [95] Let X ⊂ Z 3 be an object. The set D = {d 1 , ..., d k } ⊆ X is called a simple set of X if D can be arranged in a sequence d i 1 , ..., d i k in which each d i j is simple for X \ {d i 1 , ..., d i j-1 } for j = 1, ..., k.
(By definition, let the empty set be simple.) 1. Only simple points are deleted.

M 1 M 2 M 3 M 4 M 5 M 6
2. If two corners belonging to the object, p and q, of a 2 × 2 square in Z 3 are deleted, then the set {p, q} is simple.

3. If three corners belonging to the object, p, q and r, of a 2 × 2 square in Z 3 are deleted, then the set {p, q, r} is simple.

4. If four corners belonging to the object, p,q,r and s, of a 2 × 2 square in Z 3 are deleted, then the set {p, q, r, s} is simple.

5.

No component of the object contained in a 2 × 2 × 2 cube in Z 3 can be deleted completely.

Theorem 10. [START_REF] Palágyi | A 3D 6-subiteration thinning algorithm for extracting medial lines[END_REF] The thinning algorithm PKD6 preserves topology, in the sense of Theorem 9.

New Method

We propose a 6-directional thinning algorithm based on isthmuses. Our algorithm is designed in two steps, separating the thinning process (using a 6-directional approach) and the isthmuses detection. Each iteration consists of:

1. updating the constraint set, by adding points of the object detected as isthmuses, 2. removing deletable points of the object which are not in the constraint set (see Algorithm 11).

Notice that this design is inspired by the one of critical kernel based methods (see Algorithms 6 and 7) reported in the previous chapter.

Algorithm 11: D6I1D Data:

X ⊆ Z 3 , K ⊆ X Result: A skeleton of X repeat K = K ∪ Ψ I1D (X) foreach direction d in (U, D, N, S, E, W ) do X = τ d (X, K) until stability ; return X
This design allows for the preservation of some points defined a priori by the initial set K (which may be empty).

The function Ψ I1D (X) returns the set of all the 1D-isthmuses of X (see Section 4.2.2.2).

The function τ d (X, K) returns the set X \ S, S being the set of all points from X \ K respecting deletion conditions in X for direction d. This function is defined in the sequel of this section.

Design of Deletion Condition Masks

Contrary to deletion conditions used in directional algorithms found in the literature, those of τ d permit the deletion of curve end points, the constraint points being considered separately.

Deletion conditions can be represented using a set of masks: a point respects deletion conditions if and only if it matches at least one of the masks. In this section, we present the masks used by τ d , obtained by modifying those used by PKD6.

In order to separate geometrical and topological conditions, we have to allow in τ d the deletion of ending points (i.e. points with only one neighbor in the object). For this purpose, we add three masks to M d , representing ending points which have to be removed for the direction d.

These masks are shown in Figure 6.2 for the case d = U .

Proving that our algorithms preserve topology (Proposition 11 and Proposition 12) can be done using the framework of critical kernels, introduced by Bertrand [START_REF] Bertrand | On critical kernels[END_REF]. This framework is indeed the most powerful one for both proving and designing n-dimensional homotopic thinning algorithms. However, in this particular case it is simpler to build on the proof that was given by Palágyi and Kuba for PKD6, as our algorithm only slightly differs from the latter.

M α M β M γ
Proposition 11. Algorithm PKD6 with mask sets M d extended by adding convenient rotations of M α , M β and M γ preserves topology.

Proof. We prove that all conditions of Theorem 9 hold.

Let X ⊆ Z 3 and p ∈ X. If p matches a mask in M U , then we know that it is simple from Theorem 10. If p matches a mask in {M α , M β , M γ }, then p is obviously simple.

Let Z be a set of object points contained in a 2 × 2 square in Z 3 , or being an object component

contained in a 2 × 2 × 2 cube in Z 3 , with |Z| ≥ 2.
If all points of Z match masks of M U , then we know that all conditions of Theorem 9 hold, from Theorem 10.

Otherwise, if one point p of Z matches a mask in {M α , M β , M γ }, then clearly |Z| = 2. Let q be the other point in Z. By examination of all configurations, it can be seen that any mask in M U ∪ {M α , M β , M γ } positioned on q does not match, thus Z cannot be deleted.

Mask Set Reduction

The set of masks M U ∪ {M α , M β , M γ } can be compacted as proposed in Figure 6.3. We can observe that mask M ′ 1 matches exactly the same configurations as those matched by masks M 1 and M α . In the same way, mask M ′ 5 matches exactly the same configurations as those matched by masks M Algorithm 12: τ d Data:

X ⊆ Z 3 , K ⊆ X Result: A subset of X R = set of points matching a mask belonging to M ′ d R = R \ K return X \ R M ′ 1 M ′ 2 M ′ 3 M ′ 4 M ′ 5 M ′ 6 M ′ 7 Figure 6.3: Final set of masks M ′ U for direction UP.
Proposition 12. Algorithm D6I1D preserves topology.

Proof. From Proposition 11, algorithm PKD6 with mask sets M d extended by adding convenient rotations of M α , M β and M γ preserves topology. We observe that:

• D6I1D differs only from this algorithm by constraining some points to be preserved;

• if Theorem 9 holds for an operation A, then it also holds for an operation B that removes only a subset of the points removed by A.

These observations lead to the conclusion that D6I1D preserves topology.

Comparative Results

We implemented D6I1D with all the optimizations proposed in the previous chapter: border list, LUT and CCSU. We apply the same procedure as proposed in Section 5.5 to observe the usability of our algorithm in a real context of motion capture. The results are shown in Figure 6.4.

We observe that the computation speed of D6I1D is very close to the one of PKD6, and can be used in real time. Concerning the quality of the resulting skeleton, as no general quality measurement exists (to our knowledge), we propose to observe the results of the three studied thinning, presented in Figures 6.5 and 6.6. It can be observed that our algorithm provides results with less spurious branches than PKD6, and in certain cases less than ACK3 (see images 20 and 23, for example).

However, ACK3 provide better results for flat shape parts; like the chair (image 18) or the turbine (image 14).

Other Kinds of Skeletons

In fact, our algorithm offers great flexibility, due to the fact that constraint set detection (function Ψ) can be changed to implement other conditions. We provide two examples of such variations.

Ultimate Skeletons

In order to obtain ultimate skeletons instead of curvilinear ones, we propose a very simple variation of D6I1D, consisting in replacing Ψ I1D by the function returning empty set (K is unchanged). We call this variation D6U. Figure 6.7 shows some results obtained by using this algorithm.

Surface Skeletons

In order to obtain surface skeletons instead of curvilinear ones, we propose a very simple variation of D6I1D, consisting in replacing Ψ I1D by a function Ψ I2D returning the set of 2D-isthmuses (see Section 4.2.2.2), which can be locally detected. We call this variation D6I2D. Figure 6.8

shows some results obtained by using this algorithm. Problematics and Usual Definitions

ACK3 D6I1D P KD6

In our motion capture method, we propose to transpose the 3d shape of the subject in a graph.

Then, we will perform a matching between this graph and the model graph in order to identify the 3d position of each part described in the model.

We have two problems to solve:

1. find a way to transpose the 3d shape to a usable graph 2. find an appropriate method to match this graph with the model graph

In this chapter, we first provide usual definitions and notations about graph theory which will be used in the following. Then, we explain how we extract a data tree from the shape and what are the possible noises occurring on the data tree which have to be taken into account during the matching.

Representation of the Shape

In order to perform a fast and efficient matching between two shapes, it is very important to describe them in a concise way. Numerous shape descriptors can be found in the literature. In this section, we propose a quick overview of them. 69

Global Descriptors

A common class of descriptors is the one of global descriptors. A global descriptor represents the whole object, often by a vector of numbers. Global descriptors include:

• classical measurements: perimeter (in 2D), area, volume (in 3D), area and volume of the convex hull

• ratios between measurements

• Fourrier transform

• moments

Global descriptors are very useful for a quick similarity measurement between two shapes, but not suitable for matching.

Moments

Moments represent the shape by a numerical sequence (finite or not). The longer the sequence, the more precise the description. If the Cartesian moments are not used as a shape descriptor because non-orthogonal, it is not the case of orthogonal ones, like Legendre moments, used by Mukundan and Ramakrishnan in [START_REF] Mukundan | Moment functions in image analysis: theory and applications[END_REF] and Zernike moments used by Novotni and Klein in [START_REF] Novotni | Shape retrieval using 3D Zernike descriptors[END_REF].

Features Distribution

An alternative to global descriptors is the consideration of the distribution of simple features of a shape. Shape distributions have been introduced by Osada et al. in [START_REF] Osada | Shape distributions[END_REF]. The authors propose to compute the distribution of simple features of random surface points of the shape:

distance between two random points, angles and triangle area for three random points, etc.... The similarity between two shapes is then obtained by computing an appropriate distance between their feature distribution. According to the authors, the feature distribution based on distance between two random points (called "D2 shape distribution") is the one providing experimentally the best results.

Ip et al. proposed in [START_REF] Ip | Using shape distributions to compare solid models[END_REF] to improve the D2 shape distribution, by separating the distributions of distances for two points, depending on whether the segment between them is fully in the shape ("IN shape distribution"), fully out of it ("OUT shape distribution"), or a mix ("MIXED shape distribution"). The similarity between two shapes is here provided by a weighted combination of the similarities of the four distributions (D2, IN, OUT et MIXED).

Contours

Contours constitute a very intuitive way to describe a 2D shape. In the literature, we can distinguish two different ways to use the contours as a descriptor. Describing the contours as a curve (Bunke and Bühler [START_REF] Bunke | Applications of approximate string matching to 2D shape recognition[END_REF], Gdalyahu and Weinshall [START_REF] Gdalyahu | Flexible syntactic matching of curves and its application to automatic hierarchical classification of silhouettes[END_REF], Kaygin and Bulut [START_REF] Kaygin | Shape recognition using attributed string matching with polygon vertices as the primitives[END_REF]), it is possible to measure the similarity of two shapes by searching the minimal transformation changing one curve to the other. Another way to use the contours is to code it by coefficients in an appropriate base and the similarity is computed by comparing the coefficients (Brun and Pruvot [START_REF] Brun | Hierarchical Matching Using Combinatorial Pyramid Framework[END_REF], Shen

and Shen [START_REF] Shen | Orthogonal Legendre moments and their calculation[END_REF]). Notice that obviously contours description is only usable in 2D.

Spatial Maps

Spatial maps are a way to represent shapes by their spatial distribution. Ankerst et al. [START_REF] Ankerst | 3D shape histograms for similarity search and classification in spatial databases[END_REF] propose to segment the space containing the shape in sectors and to compute the histogram of the distribution of the shape in these sectors (a bin representing the filling of the corresponding sector by the shape). The similarity measure is obtained by similarity computation between the histograms representing the shapes, taking into account the distance between the sectors represented by the bins. Such a method is not affine invariant and requires a pose normalization. Kazhdan et al. [START_REF] Kazhdan | Rotation invariant spherical harmonic representation of 3D shape descriptors[END_REF] propose a decomposition in spherical harmonics, providing a rotation invariant descriptor.

Skeleton Based Approaches

Skeleton is an interesting descriptor for several reasons. First, it describes efficiently and concisely the shape and preserves its topology. Secondly, it can be simplified, e.g. by direct graph representation preserving only the shape structure information, or on the opposite can be enriched by adding information on branches, e.g. distance to the complementary. Finally, it yields a lot of different approaches, providing similarity measurement, parts matching, or both.

Shock Graphs

From our best knowledge, shock graphs are the most famous skeleton-based approach for shape matching. Siddiqi et Kimia introduced the shock graphs in [START_REF] Siddiqi | A shock grammar for recognition[END_REF] in 1996: considering the skeleton of a 2D shape, the authors propose to segment its branches in regard to the evolution of radius of the maximal balls along them and for this purpose define four types of evolution:

• type 1 : continuous and monotone evolution of the radius

• type 2 : radius is locally minimal

• type 3 : radius is constant

• type 4 : radius is locally maximal Then, a graph is built from these segments, considering types 1 and 3 segments as edges and types 2 and 4 as vertices. Such a graph can be reduced to a tree. The matching and similarity measurement between the shapes are computed from the trees representing them.

Skeleton Graph

In 2004, Di Ruberto [START_REF] Ruberto | Recognition of shapes by attributed skeletal graphs[END_REF] proposed another representation of the skeleton by a graph: the skeletal graph. Each vertex of a skeletal graph corresponds to an ending point (point with one and only one neighbor) of the skeleton, or to an intersection of branches (point with more than two neighbors). Each edge corresponds to the branch linking such points. Different informations are linked to both vertices (radius of maximal ball) and edges (curvature variation of the corresponding branch, length, orientation, ratio between the size of the branch and the size of the skeleton, etc...). Once again, the matching and the similarity measurement are done via the graphs, using a gradual assignment.

Skeletal Shape-Measure

Torsello and Hancock propose in [START_REF] Torsello | A skeletal measure of 2D shape similarity[END_REF] to assign to each branch of the skeleton two segments of the contours, bounded by the points common to the bi-tangents circles centered at the extremities of the branch. Several informations are extracted from these contour segments: curvature evolution, area of the shape part associated to the branch. A measure is computed using these informations, with interesting properties:

• a branch representing a small part of the contour will have a measure smaller than the one of an "important" branch

• the measure homogeneously variates in function of the deformations of the object

• the measure is bending invariant, thus interesting for the matching of articulated shapes

The authors use this measure to improve shock graphs. Notice that, in reason of contours consideration, such measure cannot be directly used for 3D shapes.

Skeletal Segments

In 2008, Goh [START_REF] Goh | Strategies for shape matching using skeletons[END_REF] proposed a new way to consider the skeleton, by dividing it into "primary skeletons" in regard of different criteria. These primary skeletons can then be merged depending to their visual continuity into "secondary skeletons". Several measures (local and global) are linked to these two kinds of skeletons and their combination provides them a weight. The matching of the skeleton is then done by a variation of the Hungarian algorithm on the matrix of these weights.

Skeleton Paths

Bai and Latecki propose in [10] a new way to find a matching between two skeletons, considering all the skeleton paths between ending points. Each path is segmented in M equidistant points, and for each point, the radius of the maximal ball centered in this point is considered. Each path is thus represented by a sequence of M radii.

The authors provide a matching approach using this path representation. A distance is computed for each couple of paths (one in each skeleton). For each couple of ending points (one in each skeleton), a distance can be computed, using a new method (called "Optimal Subsequence

Bijection") on the set of the distances between all the paths starting by one of the considered points. Finally, the distance between the skeletons is computed using the Hungarian algorithm on the matrix of the distances between each couple of ending points. The Hungarian algorithm provides optimal correspondences for each point.

Choice of the Description

In the case of our application, the problem is slightly different: the matching and similarity computation are not between two shapes, but between a manually described model and a shape.

Furthermore, the model has to be simple to express in order to allow for an easy use of our method. Finally, the acquisition of the shape implies different kinds of noise which have to be taken into consideration, and prevent the use of some informations.

Skeletons are interesting in our case, as they allow for both a similarity measurement and a matching between two shapes. Considering the approaches previously described, the skeletal graphs proposed by Di Ruberto [START_REF] Ruberto | Recognition of shapes by attributed skeletal graphs[END_REF] are the most adapted, as their definition is intuitive and a similar model is easy to define by the user. Concerning the stored informations on the edges, the length of the branches (and then, the ratio of their length to the size of the skeleton) is the only information which is not noised by the acquisition. We propose to use such description for our method, and explain in the sequel of this chapter how to define it and what are the possible noises encountered which have to be taken into account.

Usual Definitions on Graphs

As our method will be based on graphs, we propose in this section a brief recall on the basic graph definitions and the notations we will use in this part.

Undirected graphs definitions

Undirected graphs

An undirected graph is a pair (V, E), where V is a finite set, and E a subset of {{x, y}, x ∈ V, y ∈

V, x = y}. An element of E is called an edge, an element of V is called a vertex.

If {x, y} ∈ E, then x and y are said to be adjacent or neighbours. The set of all neighbours of

x is denoted by N (x). The number of vertices adjacent to a vertex v is called the degree of v, and is denoted by deg(v). For example, in

Paths and cycles

Let G = (V, E) be an undirected graph, and let x, y be in V , a walk from x to y in G is a

sequence of vertices v 0 , ..., v k such that x = v 0 , y = v k and {v i-1 , v i } ∈ E, 1 ≤ i ≤ k.
The number k is called the length of the walk. If k = 0 the walk is called a trivial walk. The walk

is closed if x = y.
A path is a walk in which no vertex occurs more than once in the sequence of vertices of the path (except possibly x = y). A non-trivial closed path in which all edges are distinct is called a cycle.

Trees and forests

A graph is connected if for all {x, y} ⊂ V , a path from x to y exists in G. A tree is a connected graph with no cycles. A simple path from x to y in a tree is unique and is denoted by π(x, y).

A graph with no cycles is called a forest, each of its connected components being a tree.

Directed graphs definitions

Directed graphs

A directed graph is a pair (V, A), where V is a finite set, and A is a subset of V × V . An element of A is called an arc, an element of V is called a vertex. Let G = (V, A) be a directed graph, and let x, y be in V , a path from x to y in G is a sequence

of vertices s 0 , ..., s k such that x = v 0 , y = v k and (v i-1 , v i ) ∈ A, 1 ≤ i ≤ k.

Associated undirected graphs

The undirected graph associated to G is the undirected graph G ′ = (V, E), such that {x, y} ∈ E

Rooted trees and rooted forests

Let G be a directed graph. A vertex r ∈ V is a root of G if for all x ∈ V \ {r}, a path from r to x in G exists. The directed graph G is antisymmetric if for all (x, y) ∈ A such that x = y,

(y, x) / ∈ A.
The directed graph G is a rooted tree (with root r) if r is the unique root of G, G is antisymmetric and if the undirected graph associated to G is a tree.

A directed graph, where each of its connected components is a rooted tree, is called a rooted forest.

Remark 13. in order to avoid possible confusion, trees and forests for undirected case will be respectively called unrooted trees and unrooted forest.

Let G = (V, A) be a rooted tree. If (y, x) ∈ A, we say that y is the parent of x (denoted by par(x)), and that x is a child of y. The set of all children of y is denoted by C(y). A vertex with no children is called a leaf.

A rooted forest is ordered if for all of its vertices, the children are ordered.

The length of a path between the root and a vertex v is called the depth of v. The maximum depth of all the vertices is called the height of the tree.

The vertices on the path from the root to a vertex x are called the ancestors of x. We denote the set of the ancestors of x by anc(x).

Common definitions

Unless otherwise indicated, all the other definitions and notations are similar for both directed and undirected graphs. We present these definitions for directed graphs, the versions for undirected graphs can be obtained by replacing arcs by edges. 7.6 for some examples). From the top to the bottom: an example of isomorphism between undirected graphs (dotted lines represent the bijection), an example of isomorphism between directed graphs, and an example of non isomorphism between directed trees (bold vertex cannot find any correspondence respecting the definition).

Isomorphism

Two graphs G = (V G , A G ) and G ′ = (V G ′ , A G ′ ) are said to be isomorphic if there exists a bijection f : V G → V G ′ , such as for any pair (x, y) ∈ V G × V G , (x, y) ∈ A G if and only if (f (x), f (y)) ∈ A G ′ (see Figure

Attributed graphs

An attributed graph is a quadruplet (V, A, υ V , υ A ), where V is a finite set, A a subset of V × V , and υ V (respectively υ A ) a mapping which associate each element of V (respectively A) to an attribute.

An attribute can be a label (an element of a discrete finite set), a real value, or any other information.

A vertex-attributed graph is a triplet (V, A, υ V ) and an edge-attributed graph is a triplet (V, A, υ A ), in regard of the above notations.

Labeled graphs

A labeled graph is a vertex-attributed graph where all attributes are labels (see Figure 7.7 for an example). In order to avoid confusion with vertex identifier, the label associated to each vertex is write between brackets. In the case where the identifier must be represented, we use the representation shown on the left graph: the identifier with the label between brackets in exponent. In other case, we just write the label between brackets, like in the right graph.

a [u] b [v] c [c] d [λ] [u] [v] [c] [λ]

Weighted graphs

A weighted graph is an edge-attributed graph where all attributes are non negative reals (see In a weighted tree, the weight of the unique path from x to y, denoted by ω(x, y), is the sum of the weights of all arcs traversed in the path.

Data Tree Extraction

In this section, we will first present the required properties of the graph representation of the shape (called data tree). Then, we propose a method to construct the data tree from the skeleton of the shape. Finally we enumerate the different types of noise occurring in the data tree, which have to be taken into consideration for an efficient matching.

Data Tree Specificities

From the definitions presented in the above parts, we can now define formally the graph used in a priori model definition (called the model tree or pattern tree) as an attributed unrooted tree

M = (V M , E M , Σ M , ω M )
, with Σ M a mapping from vertices to labels representing parts of the model and ω M a mapping from edges to non negative reals representing the distances between the parts.

We can make two observations:

• To match the model tree, the data tree must have approximately the same topology than the model tree, and edge or vertex similar information.

• As the final aim is to obtain spatial positions of the model parts, the data tree have to store spatial informations in the same way that the parts are stored in the model.

The above observations involve that the data tree have to be an attributed unrooted tree D = (V D , E D , Ξ D , ω D ), with Ξ D a mapping from vertices to 3d positions of parts of the shape and ω D a mapping from edges to non negative reals representing the distances between the 3d positions.

Then, the matching can be performed between the weighted unrooted tree (V M , E M , ω M ) associated to the model tree and the weighted unrooted tree (V D , E D , ω D ) associated to the data tree.

The obtained matching give, for any parts represented by vertices in V M , a spatial position represented by vertices in V D .

Data Tree Construction Design

Now we have defined the specificities of the data tree, we have to find a way to build it from the shape. We consider two steps in the design of the building:

1. find a method which automatically provides relevant 3d positions in the shape.

2. find a method which automatically provides both the links between these 3d positions and their appropriate distance.

Figure 7.9 shows areas where we have to find usable 3d positions in the case of hand model (at left) and full human body model (at right). We can observe there exists in all areas at least one intersection or end of branches of curvilinear skeleton of the shape (represented in orange).

Furthermore, we can see that areas are linked in a similar way than in the model by branches of the skeleton.

From these observations, we can deduce an efficient method to build the data tree, following the two steps described above:

From Skeleton to Data Tree

The skeleton points can be classified into three classes, in regard of their number of neighbors included in the skeleton:

1. zero or one neighbor → E-points class (ending points)

2. exactly two neighbors → B-points class (branch points)

3. strictly more than two neighbors → I-points class (intersection points) In order to construct the data tree of the shape, we will use its asymmetric skeleton, in reason of the above observation. The construction can be done using Algorithm 13, based on detection of connected components of points of a same class (see Figure 7.11 for an example). Empirically, we can considered that data trees obtained by this method have a maximal degree bounded by 6. It is interesting and important to notice it, as a bounded maximal degree allows a better complexity for some matching algorithms.

Data Tree Noises

Due to the acquisition method, several kinds of noise and deformities can appear in the data tree.

Ghost limbs and Spurious branches

Due to the reconstruction from silhouettes, parts of space cannot be carved, resulting in "limbs" of the object which do not exist on the real model. Due to the skeletonization algorithm and to the amount of noise of the shape surface, spurious branches without important topological signification can appear on skeleton. These noises are illustrated in Figure 7.13.

Useless 2-degree vertices

Vertices with exactly two neighbors are not useful to describe the topology of a shape, and then uselessly split an edge (and its weight) into two parts, making difficult a good matching. This kind of vertices can appear when removing spurious branches or ghosts limbs. See

Splitted vertices

Vertices with more than three neighbors in the model tree can correspond to a cluster of vertices linked by weakly weighted edges in the data tree, due to the skeletonization algorithm. See Chapter 8

State of the Art of Edit-Based Matching

Graphs are among the most common and well-studied combinatorial structures in computer science. In particular, the problems of matching graphs and measure the similarity between them occurs in several diverse areas such as compiler design [START_REF] Aho | Compilers: principles, techniques, and tools[END_REF], image processing [START_REF] Samet | Distance Transform for Images Represented by Quadtrees[END_REF], signal processing [START_REF] Cheng | Waveform correlation by tree matching[END_REF], molecular biology [START_REF] Shapiro | Comparing multiple RNA secondary structures using tree comparisons[END_REF], chemistry [START_REF] Stobaugh | Chemical substructure searching[END_REF], and many others.

Similarity and Matching between Graphs

In the case of shape recognition, two graph problems are considered:

• measuring the similarity between two graphs

• matching two graphs

For our application, we request both matching (to find associated positions for each part of the subject) and similarity measurement (to evaluate the significance of the matching). In our context, the graphs to match are acyclic, it is interesting as specific solutions with better complexities have been proposed for such problems on rooted and unrooted trees.

Graph Matching

Graph matching is a very difficult problem and the algorithms providing them are very time consuming. In the field of shape matching, the representations are most often acyclic graphs or rooted trees, which are easier to match.

85

Tree matching problems can be separate in different classes of problems:

• exact matching problems consist in finding exact patterns in a tree, or finding common subparts of two trees. Examples of exact matching problems are:

-subtree isomorphism problem [START_REF] Matula | An algorithm for subtree identification[END_REF][START_REF] Matula | Subtree Isomorphism in 0(n 5/2 )[END_REF][START_REF] Shamir | Faster subtree isomorphism[END_REF]: given a pattern tree P and a data tree D, find a subtree of D which is isomorphic to P or decide that there is no such tree.

-subtree homeomorphism problem [START_REF] Chung | O(n 2.5 ) time algorithms for the subgraph homeomorphism problem on trees[END_REF][START_REF] Schlieder | Approximate tree embedding for querying XML data[END_REF]: a variant of the subtree isomorphism problem, where 2-degree vertices can be deleted from the data tree.

-maximal subtree isomorphism [START_REF] Pelillo | Matching hierarchical structures using association graphs[END_REF]: given a pattern tree P and a data tree D, find the biggest subtree of D which is isomorphic to a subtree of P .

• inexact matching problems consist in finding the best approximation of a given pattern in a data tree, or finding the mapping from a tree to another which is the closest to the isomorphism. Different classes of inexact matching exist:

-one-to-one matching requires that for two trees, each vertex in both trees match with at most one vertex of the other tree.

-one-to-many and many-to-one matchings allow multiple vertices of only one tree to match with a same vertex of the other tree.

-many-to-many matching allow clusters of vertices in each tree to match together.

Inexact graph and trees matching are well suited for shape matching and several approaches have been developed to compute them efficiently.

Association Graph

Association graphs [12,[START_REF] Pelillo | Matching free trees, maximal cliques, and monotone game dynamics[END_REF][START_REF] Pelillo | Many-to-many matching of attributed trees using association graphs and game dynamics[END_REF][START_REF] Pelillo | Matching hierarchical structures using association graphs[END_REF] are an interesting tool for both exact and inexact graphs matching. Let consider two graphs G = (V, E) and G ′ = (V ′ , E ′ ), an association graph of G and G ′ is a graph where the vertex set is V × V ′ and an edge link two vertices (v 1 , v ′ 1 ) and (v 2 , v ′ 2 ) if and only if a same relation exists between v 1 and v 2 and between v ′ 1 and v ′ 2 .

The matching between G and G ′ is then computed by searching the maximal clique in their association graph. The kind of matching will depend of the relation considered for the linking in the association graph.

Graph Eigenspace

Another way to obtain an inexact graph matching is to use a graph eigenspace [START_REF] Kosinov | Inexact Multisubgraph matching using Graph Eigenspace and Clustering Models[END_REF], where each vertex of both graphs is represented by an eigenvector, which can be seen as a point in n-dimensions. Then, the computation of the clusters of eigenvectors in this eigenspace provide the correspondences between the two sets of vertices.

Similarity Measurement

From our best knowledge, the only methods providing a similarity measurement without providing a matching are those based on graph kernels. Graphs kernels are very useful, as they provide an easy way to classify and compare graphs. For two graphs, computing their kernel consists in projecting them in a Hilbert space and to compute the scalar product of these projections.

Different projection functions have been proposed, based on random walks [START_REF] Gärtner | Exponential and geometric kernels for graphs[END_REF][START_REF] Kashima | Marginalized kernels between labeled graphs[END_REF][START_REF] Mahé | Extensions of marginalized graph kernels[END_REF][START_REF] Vishwanathan | Graph kernels[END_REF],

tree patterns [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF][START_REF] Ramon | Expressivity versus efficiency of graph kernels[END_REF], or bags of paths [START_REF] Borgwardt | Shortest-path kernels on graphs[END_REF][START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF]. Some kernels have also been proposed for trees [START_REF] Vishwanathan | Fast kernels for string and tree matching[END_REF].

Methods Providing both Graph Similarity and Matching

In the trees context, an approach widely used for the matching and similarity measurement is the one using edit-based distances. Edit-based distances are a very powerful and highly studied tool to solve inexact matching problems. It consists in finding the cheapest sequence of operations which transforms the pattern tree into the data tree. In order to solve our specific problem of tree matching, we have developed a new edit-based distance and its associated matching.

In this chapter, we propose a state of art of edit-based distances, in order to position our method.

Edit-based Distance Basics

In this section, we introduce the important notions necessary to easily understand the following state of art of edit-based distances.

Edit Operations

The most fundamental notion is obviously the notion of edit operation. An edit operation is an atomic treatment which will modify the tree on which it is applied. Depending from both the method and the kind of trees, the defined edit operations will differ. We propose to express only the classical ones on labeled rooted trees, widely used in the literature, and those on weighted trees, as our method has to use it.

Labeled rooted tree edit operations For a labeled rooted tree G = (V, A, Σ), classical edit operations are:

Relabel: Change the label of a vertex v ∈ V .

Delete: Delete a non-root vertex v ∈ V with parent v ′ , making the children of v become the children of v ′ .

Insert: The dual of delete. Insert a new vertex v as a child of v ′ , making v the parent of a subset of the children of v ′ .

Examples of these operations are shown in Figure 8.1.

[b] [a] [a] [b] [b] [c] [a] [b] [b] [a] [b] [c] [c] [b] [a] [c] [c]

Cost Function

A cost function is associated to the edit operations. In the literature about labeled trees, it is usual to consider that a relabeling with a different label, the deletion and the insertion have a cost of 1, and the relabeling with the same label have a cost of 0.

For our weighted graph purpose, we use a function γ : R + × R + → R + . For a given operation, its cost is γ(w, w ′ ), where w (respectively w ′ ) is the total weight of the arcs involved in the operation before (respectively, after) its application.

As a consequence, the cost of a resizement can be denoted by γ(w, w ′ ), where w is the former weight of the resized arc and w ′ is the new weight, the cost of a deletion can be denoted by γ(w, 0), where w is the weight of the deleted arc, and the cost of an insertion, γ(0, w), where w is the weight of the created arc.

We choose γ to be a translation invariant metric, and is required to satisfy the following constraints:

• γ(w, w ′ ) ≥ 0 (non-negativity constraint)

• γ(w, w ′ ) = 0 if and only if w = w ′ (identity of indiscernible constraint)

• γ(w, w ′ ) = γ(w ′ , w) (symmetry constraint)

• γ(w, w ′ ) ≤ γ(w, x) + γ(x, w ′ ) (triangle inequality constraint)

• γ(w, w ′ ) = γ(w + x, w ′ + x) (translation invariance constraint)

Tree Edit Distance

An edit script S between two graphs G 1 and G 2 is a sequence of edit operations turning G 1 into G 2 . The cost of S is the sum of the costs of the operations in S. An optimal edit script between G 1 and G 2 is an edit script between G 1 and G 2 of minimum cost and this cost is the tree edit distance.

Mapping

An edit distance mapping (also just called a mapping) between two graphs G 1 and G 2 is a representation of the edit script.

Mapping for labeled rooted trees In the literature, a mapping between two labeled ordered rooted trees

T 1 = (V 1 , A 1 , Σ 1 ) and T 2 = (V 2 , A 2 , Σ 2
) can be formally defined as the triplet (M, T 1 , T 2 ), with M ⊆ V 1 × V 2 and for any pair (v 1 , w 1 ), (v 2 , w 2 ) ∈ M :

• v 1 = v 2 if and only if w 1 = w 2 (one-to-one condition)

• v 1 ∈ anc(v 2 ) if and only if w 1 ∈ anc(w 2 ) (ancestor condition)

• v 1 is to the left of v 2 if and only if w 1 is to the left of w 2 (sibling condition)

An element of M is called a mapping line.

The mapping for labeled unordered rooted tree is defined in the same way, but without the sibling condition.

An example of mapping between labeled rooted trees is shown in Figure 8.3.

[a]

[b] [b] [c] [a] [b] [b] [c] G 1 G ′ 1 [a] [b] [b] [c] [a] [b] [b] [c] G 2 G ′ 2 Figure 8.3: Left: A mapping from G 1 to G ′ 1 .
The dotted lines represent the mapping. Vertices in G 1 not touched by a mapping line have to be deleted. Vertices in G ′

1 not touched by a mapping line have to be inserted. The mapping shows a way to transform G 1 to G ′ 1 . Right: An invalidate mapping from G 2 to G ′ 2 : the ancestor condition is not respected.

Remark 15. In the sequel, we will present some mappings with no mapping line between the roots. As it is in contradiction with the definition of delete operation, which occurs only on non-rooted vertices, we have to consider that the deletion of the root r 1 of T 1 , which split T 1 in the forest of all subtrees rooted in a vertex ∈ C(r 1 ), is possible in two cases:

1. the root r 2 of T 2 is deleted too, then the mapping occurs on the forests of both trees.

2. the root r 2 of T 2 is mapped to a vertex v in T 1 , such that v is an ancestor of all vertices mapped with vertices of T 2 (due to ancestor condition).

Mapping for weighted trees In the case of weighted trees, we can adapt the above definitions of mapping.

A mapping between two weighted rooted trees T 1 = (V 1 , A 1 , ω 1 ) and T 2 = (V 2 , A 2 , ω 2 ) can be formally define as the triplet (M, T 1 , T 2 ), with M ⊆ A 1 × A 2 and for any pair (a

1 = (v 1 , v ′ 1 ), b 1 = (w 1 , w ′ 1 )), (a 2 = (v 2 , v ′ 2 ), b 2 = (w 2 , w ′ 2 )) ∈ M : • a 1 = a 2 if and only if b 1 = b 2 (one-to-one condition) • v ′ 1 ∈ anc(v ′ 2 ) if and only if w ′ 1 ∈ anc(w ′ 2 ) (ancestor condition)
An example of invalidate mapping between weighted trees is shown in Figure 8. c is an ancestor of d, contrary to g which is not an ancestor of h).

Mapping Cost By definition, the cost of a mapping is equal to the cost of the associated edit script.

Complexities notations

In the sequel of the chapter, we will give the complexity of algorithms computing edit distances between two trees T 1 and T 2 , mainly in function of their size, their number of leaves, their height and their maximal degree. For i ∈ {1, 2}:

• the size of T i is denoted by S i

• the number of leaves in T i is denoted by L i

• the height of T i is denoted by H i

• the maximum degree of T i is denoted by D i

General Edit Distance

The edit distance problem between labeled ordered rooted trees was introduced by Tai in [START_REF] Tai | The tree-to-tree correction problem[END_REF].

It consists in finding the optimal edit script and the associated mapping (called the optimal mapping) between labeled ordered rooted trees. Notice that the optimal mapping is not necessary unique. In [START_REF] Zhang | Simple fast algorithms for the editing distance between trees and related problems[END_REF], Zhang and Shasha proposed an algorithm with an O(S 1 S 2 min(L 1 , H 1 )min(L 2 , H 2 ))

time complexity. This algorithm was modified by Klein in [START_REF] Klein | Computing the edit-distance between unrooted ordered trees[END_REF] to get a better worst case time complexity in O(S 2 1 S 2 logS 2 ).

In 2001, Chen [START_REF] Chen | New algorithm for ordered tree-to-tree correction problem[END_REF] proposed an algorithm with O(S 1 S 2 +L 2 1 S 2 +L 2.5 1 L 2 ), which improves previous bounds for certain kinds of trees (for example, if L i < √ S i , i = 1, 2). More recently, Demaine et al. [START_REF] Demaine | An O (nˆ3)-Time Algorithm for Tree Edit Distance[END_REF] proposed an algorithm in O((S 1 + S 2 ) 3 ) time complexity.

For the edit distance problem between labeled unordered rooted trees, Zhang et al. [START_REF] Zhang | On the editing distance between unordered labeled trees[END_REF] shown that the problem was NP-complete, and proposed an algorithm with an O(S

1 S 2 + L 1 !2 L 1 (L 3 1 + H 2
2 )S 2 ) time complexity.

Other Edit Distances

Numerous other edit distances have been proposed in order to reduce the complexity of the general edit distance problem or to provide special characteristics to the associated mapping.

These edit distances are derived from the general edit distance by adding constraints on the mapping, on the edit operations or on the edit script. We will here propose a quick overview of these edit distances.

Alignment Distance

Jiang et al. proposed in [START_REF] Jiang | Alignment of trees -an alternative to tree edit[END_REF] the notion of alignment of trees.

Let T 1 = (V 1 , A 1 , Σ 1 )
and

T 2 = (V 2 , A 2 , Σ 2
) be two labeled rooted trees. We denote by λ a space, which is a label not occurring in T 1 and T 2 .

An alignment of T 1 and T 2 is obtained by first inserting vertices labeled by λ in T 1 and T 2 such that the two resulting trees T ′ 1 and T ′ 2 are isomorphic, then drawing a mapping line from each vertex in T ′ 1 to its corresponding vertex in T ′ 2 . See left part of Figure 8.5 for an example.

An alignment mapping can be provided between T 1 and T 2 by drawing a mapping line between In 2001, Jansson and Lingas [START_REF] Jansson | A Fast Algorithm for Optimal Alignment between Similar Ordered Trees[END_REF] proposed an algorithm to compute the optimal alignment between similar trees. Two trees are similar if there is an optimal alignment between them using at most s spaces. Their algorithm had a O(S 1 + S 2 )log(S 1 + S 2 )(D 1 + D 2 ) 4 s 2 ) time complexity.

each pair (v 1 , w 1 ) ∈ V 1 × V 2 ,
Remark 16. In his thesis [START_REF] Höchsmann | The tree alignment model: algorithms, implementations and applications for the analysis of RNA secondary structures[END_REF], Höchsmann shown that an alignment distance is not a distance, due to non respect of triangle inequality constraint.

[a]

[a]

[b] [b] [c] [c] [d] [d] [e] [e] [f ] [f ] [g] [g] [a] [a] [b] [b] [c] [c] [d] [d] [e] [e] [f ] [f ] [g] [g] [λ] [λ] [λ] [λ] [λ] [λ]

Constrained/Isolated Subtree Edit Distance

According to Wang and Zhang [START_REF] Wang | Finding similar consensus between trees: an algorithm and a distance hierarchy[END_REF], the isolated subtree edit distance (also called the constrained edit distance in the litterature [START_REF] Bille | A survey on tree edit distance and related problems[END_REF]) was first proposed in 1988 by Tanaka and Tanaka [START_REF] Tanaka | The tree-to-tree editing problem[END_REF] and involves that disjoint subtrees are mapped with disjoint subtrees.

A constrained mapping can be obtained by adding a condition on mapping. In a rooted tree T = (V, A), we denote by nca(x, y) the nearest common ancestor of x and y in T , which is the vertex in anc(x) ∩ anc(y) with maximal height. The additional condition for the mapping

(M, T 1 , T 2 ) is, for any triplet (v 1 , w 1 ), (v 2 , w 2 ), (v 3 , w 3 ) ∈ M : • nca(v 1 , v 2 ) ∈ anc(v 3 ) \ {v 3 } if and only if nca(w 1 , w 2 ) ∈ anc(w 3 ) \ {w 3 }.
An example of constrained mapping is shown in Figure 8.6.

For the application between labeled ordered rooted trees, Zhang [START_REF] Zhang | Algorithms for the constrained editing distance between ordered labeled trees and related problems* 1[END_REF] gives an algorithm with a time complexity in O(S 1 S 2 ). In [START_REF] Zhang | A constrained edit distance between unordered labeled trees[END_REF], he proposed an algorithm for the unordered case with In 2003, Ferraro and Godin [START_REF] Ferraro | Optimal mappings with minimum number of connected components in tree-to-tree comparison problems[END_REF] proposed an algorithm to compute the optimal constrained edit mapping with a minimal number of connected components in both trees, with the same complexity than Zhang. A constrained mapping between two labeled unrooted trees

a O(S 1 S 2 (D 1 + D 2 )log(D 1 + D 2 )) time complexity. v 1 v 2 v 3 v 4 v 5 v 1 v 2 v 3 v 4 v 5 w 5 w 4 w 3 w 1 w 2 w 1 w 2 w 3 w 4

Constrained

T 1 = (V 1 , E 1 , Σ 1 ) and T 2 = (V 2 , E 2 , Σ 2 )
can be formally defined as the triplet (M, T 1 , T 2 ), with M ⊆ V 1 × V 2 and for any triplet (v 1 , w 1 ), (v 2 , w 2 ), (v 3 , w 3 ) ∈ M :

• v 1 = v 2 if and only if w 1 = w 2 (one-to-one condition) • (c, c ′ ) ∈ M , with c = center(v 1 , v 2 , v 3 ) and c ′ = center(w 1 , w 2 , w 3 ) (center relationship preservation condition)
An example of invalidate constrained mapping is shown in Figure 8.8.

[a] [b] [a] [a]

[b]

[a]

[a]

Less Constrained Edit Distance

In 2001, Lu et al. [START_REF] Lu | A new measure of edit distance between labeled trees[END_REF] introduced the less-constrained edit distance, which relaxes the constrained mapping. Here, the additional condition for the mapping (M, T 1 , T 2 ) is, for any triplet (v 1 , w 1 ), (v 2 , w 2 ), (v 3 , w 3 ) ∈ M , such none of v 1 , v2 and v 3 is an ancestor of the others:

• depth of nca(v 1 , v 2 ) ≥ depth of nca(v 1 , v 3 ) and also nca(v 1 , v 3 ) = nca(v 2 , v 3 ) if

and only

if depth of nca(w 1 , w 2 ) ≥ depth of nca(w 1 , w 3 ) and also nca(w 1 , w 3 ) = nca(w 2 , w 3 )

The authors proposed an algorithm for this edit distance between labeled ordered rooted trees

with a O(S 1 S 2 D 3 1 D 3 2 (D 1 + D 2 )
) time complexity. In the unordered case, the problem is NPcomplete. An example of less constrained mapping is shown in Figure 8.9.

v 1 v 2 v 3 v 4 v 5 v 1 v 2 v 3 v 4 v 5 w 5 w 4 w 3 w 1 w 2 w 1 w 2 w 3 w 4 Figure 8
.9: Left: a mapping which is not constrained but less constrained. Right: a mapping which is not less constrained.

1-Degree/Top-Down Edit Distance

In [START_REF] Selkow | The tree-to-tree editing problem[END_REF], Selkow proposed an edit distance with for additional constraint that insertions and deletions are restricted to leaves of the tree. This edit distance is called 1-degree edit distance or top-down edit distance, due to the translation of this constraint in the mapping.

Let

T 1 = (V 1 , A 1 , Σ 1 ) and T 2 = (V 2 , A 2 , Σ 2
) be two labeled trees rooted respectively in r 1 and r 2 .

A mapping (M, T 1 , T 2 ) is a top-down mapping if for all couple (v 1 , w 1 )

∈ (V 1 \ {r 1 }) × (V 2 \ {r 2 }):
• (v 1 , w 1 ) ∈ M if and only if (par(v 1 ), par(w 1 )) ∈ M An example of top-down mapping is shown in Figure 8.10.

Yang proposed in [START_REF] Yang | Identifying syntactic differences between two programs[END_REF] an algorithm which computes the top-down edit distance between labeled ordered rooted trees with a O(S 1 S 2 ) time complexity. For our best knowledge, there does not exist algorithm for the unordered case.

v 1 v 2 v 3 v 4 v 5 w 1 w 2 w 4 w 3 w 5 v 1 v 2 v 3 v 4 v 5

Bottom-Up Edit Distance

Valiente proposed in [START_REF] Valiente | An efficient bottom-up distance between trees[END_REF] a bottom up mapping defined by adding the following constraint.

Let

T 1 = (V 1 , A 1 , Σ 1 ) and T 2 = (V 2 , A 2 , Σ 2
) be two labeled trees rooted respectively in r 1 and r 2 . A mapping (M, T 1 , T 2 ) is a bottom up mapping if for all couple (v 1 , w 1 ) ∈ M :

• (v 1 , w 1 ) ∈ M then ∀c i ∈ C(v 1 ), ∃c j ∈ C(w 1 ) such (c i , c j ) ∈ M and ∀c j ∈ C(w 1 ), ∃c i ∈ C(v 1 ) such (c i , c j ) ∈ M .
An example of bottom-up mapping is shown in Figure 8.11.

The algorithm provided by Valiente to compute the bottom up distance have a complexity in O(S 1 + S 2 ) for both ordered and unordered cases.

v 1 v 2 v 3 v 4 v 5 w 1 w 2 w 4 w 3 w 5 v 1 v 2 v 3 v 4 v 5

Inclusion of Mappings

Mappings based on general edit mapping can be hierarchically classified in regard to their constraints. For example, we can assume that all mappings derived from the general edit mapping by adding restriction are general edit mappings. Some of these relations are directly derived from the mappings definitions:

• all constrained (or isolated subtree) mappings are less constrained mappings.

• all 1-degree (or top down) mappings are 2-degree mappings.

• all 2-degree mappings are isolated subtree mappings.

• all bottom-up mappings are isolated subtree mappings.

In 2001, Wang and Zhang [START_REF] Wang | Finding similar consensus between trees: an algorithm and a distance hierarchy[END_REF] demonstrated that all isolated subtree mappings are alignment mappings. From all these relations, we can established a hierarchy on all the mappings described above. This hierarchy is shown on Figure 8.12.

General Edit Mapping

Edit Distances with Extended Set of Operations

Edit distances using additional operations have also been proposed, in order to treat specific problems, as shape matching for example.

Edge Merging and Edge Pruning

Klein et al. [START_REF] Klein | A tree-edit-distance algorithm for comparing simple, closed shapes[END_REF] proposed a modified edit distance on edge-attributed ordered rooted trees using two new operations: edge merging Merges two edges sharing a 2-degree vertex. The resulting attribute is derived from the attributes of involved edges.

edge pruning Removes an edge ended by a 1-degree vertex at one side and a 3-degree vertex v at the other side. Then, the two edges sharing v are merged.

Cut Operation

In [START_REF] Wang | Finding approximate patterns in undirected acyclic graphs[END_REF], Wang et al. proposed to add the cut operation, which consist in removing complete subtrees in one of the two trees for free. It allows to use edit distance for inexact subtree matching problem.

Horizontal/Vertical Merge and Split

In 2003, Bille [START_REF] Bille | Ordered Tree Edit Distance with Merge and Split Operations[END_REF] proposed to extend the general edit distance for labeled ordered rooted trees using merge and split operations: Bille proposed two new edit distances using these operations. The k-way horizontal edit distance consists in allowing horizontal merges and splits in the edit script, with the restriction to merge less than k vertices together. The author provides an algorithm that reaches O(S 2k+2 ) time complexity.

The vertical edit distance consists in allowing vertical merges and splits in the edit script. An algorithm is proposed with a O(S 7 ) time complexity.

We can notice that the use of merge operations involve to break the one-to-one constraint of edit mapping: this kind of mapping is many-to-many.

Discussion for our Purpose

In our purpose, none edit mapping is well adapted to take into consideration the noises, especially the 2-degree useless vertex noises and the spurious branches. However, we can consider that these two problems can be solved by adding respectively merge operation of Klein et al. [START_REF] Klein | A tree-edit-distance algorithm for comparing simple, closed shapes[END_REF] and the cut operation of Wang et al. [START_REF] Wang | Finding approximate patterns in undirected acyclic graphs[END_REF].

The last remaining problem, the splitted vertex problem, can be taken into consideration by edit mappings, using the edge deletion. Among the edit mappings, we will use the one which well preserves the topology and which has a good complexity. Isolated subtree mapping and its sub cases (regarding the mapping hierarchy previously established in Figure 8.12) are too restrictive. The alignment mapping seems a good compromise between topology preservation and computation speed.

Alignment Distance for Weighted Trees

In [START_REF] Jiang | Alignment of trees -an alternative to tree edit[END_REF], Jiang et al. proposed the alignment between vertex-labeled trees, that we transpose here for edge-weighted graphs.

Let G 1 = (V 1 , A 1 , ω 1 ) and G 2 = (V 2 , A 2 , ω 2 ) be two weighted graphs. Let G ′ 1 = (V ′ 1 , A ′ 1 , ω ′ 1 ) and G ′ 2 = (V ′ 2 , A ′ 2 , ω ′
2 ) be weighted graphs obtained by inserting arcs weighted by 0 in G 1 and G 2 , such that there exists an isomorphism

I between (V ′ 1 , A ′ 1 ) and (V ′ 2 , A ′ 2 ). The set of all couples of arcs A = {(a 1 , a 2 ); a 1 ∈ A ′ 1 , a 2 ∈ A ′ 2 , a 2 = I(a 1 )} is called an alignment of G 1 and G 2 . The cost C A of A is given by C A = (a 1 ,a 2 )∈A γ(ω ′ 1 (a 1 ), ω ′ 2 (a 2 )) . (8.1)
The minimal cost of all alignments from G 1 and G 2 , called the alignment distance, is denoted by α(G 1 , G 2 ).

An example of alignment and alignment mapping between weighted trees is shown in Figure 8.13.

G 1 G 1 G 2 G 2 G ′ 2 G ′ 1

Chapter 9 Homeomorphic Alignment

As shown in the previous chapter, there is no matching method fitted for our need in the literature, as far as we know. However, we have seen that solutions exist for the three kind of noises, considered separately:

• alignment mapping seems a good compromise between time complexity and topology preservation, and solve the problem of splitted vertices

• merging operations on edges can be a solution for useless 2-degree vertices

• cut operation is an interesting tool for spurious branches removing

In this chapter we propose a new kind of alignment, specially designed to take into account the problem of useless 2-degree vertices, by extending the set of edit operations. First, we give a formal definition of this alignment. Then, we show how to compute it efficiently for both rooted trees and unrooted trees, using bottom-up based algorithms. Finally, we show how to integrate the cut operation in this alignment, in order to solve the problem of spurious branches.

Preliminar Definitions

Merging Operation

The merging is an operation that can be applied only on arcs sharing a 2-degree vertex (for directed graphs, the vertex have to be the end of one arc and the start of the other). The merging of two arcs (u, v) and (v, w) (we say that we merge on vertex v) in a weighted graph G = (V, A, ω) consists of removing v in V , replacing (u, v) and (v, w) by (u, w) in A, weighted by ω((u, w)) = ω((u, v)) + ω((v, w)). Notice that the cost of a merging is equal to 0, as no weight variation occurs. Examples of merging for directed and undirected graphs are shown in

Homeomorphism

Two weighted graphs

G = (V G , A G , ω G ) and G ′ = (V G ′ , A G ′ , ω G ′ )

Merging Kernel

Considering that a merging on a vertex v on the graph G = (V, A, ω) does not affect the degree of any vertex in V \ {v} (by definition of merging operation) and therefore the possibility of merging this vertex, the number of possible mergings decreases by one after each merging.

In consequence, the maximal size of a sequence of merging operations, transforming G into another graph G ′ = (V ′ , A ′ , ω ′ ) is equal to the initial number of possible mergings in G.

It can be remarked that any sequence of merging operations of maximal size yields the same result. The graph resulting of such a sequence on G is called the merging kernel of G, and is denoted by MK(G).

The following definition is straightforward:

Definition 17. Two graphs G 1 = (V 1 , A 1 , ω 1 ) and G 2 = (V 2 , G 2 , ω 2 ) are homeomorphic iff MK(G 1 ) and MK(G 2 ) are isomorphic.
Examples of merging kernels of directed and undirected graphs are shown in Figures 9.5

Homeomorphic Alignment Definition

Let G 1 = (V 1 , A 1 , ω 1 ) and G 2 = (V 2 , A 2 , ω 2 ) be two weighted graphs. Let G ′ 1 = (V ′ 1 , A ′ 1 , ω ′ 1 ) and G ′ 2 = (V ′ 2 , A ′ 2 , ω ′
2 ) be weighted graphs obtained by deleting arcs in G 1 and G 2 , such that there exists an homeomorphism between G ′ 1 and G ′ 2 (not necessarily unique). Let

G ′′ 1 = (V ′′ 1 , A ′′ 1 , ω ′′ 1 ) and G ′′ 2 = (V ′′ 2 , A ′′ 2 , ω ′′
2 ) be the merging kernel of G ′ 1 and G ′ 2 , respectively. From definition 17, there exists an isomorphism I between G ′′ 1 and G ′′ 2 .

The set of all couples of arcs H = {(a, a ′ );

a ∈ A ′′ 1 , a ′ ∈ A ′′ 2 , a ′ = I(a)} is called an homeomorphic alignment of G 1 with G 2 (see Figure 9.6). The graph G ′′ 1 is called the left graph of H. The graph G ′′ 2 is called the right graph of H.
The cost C H of H is the sum of the costs of all operations used to homeomorphically align G 1 and G 2 : the deletion of arcs in G 1 and G 2 , to obtain G ′ 1 and G ′ 2 respectively, and the resizement for each arc a 1 ∈ A ′′ 1 to the weight of H(a 1 ). More formally :

C H = (a,a ′ )∈H γ(ω ′′ 1 (a), ω ′′ 2 (a ′ )) + a d ∈A 1 \A ′ 1 γ(ω 1 (a d ), 0) + a ′ d ∈A 2 \A ′ 2 γ(0, ω 2 (a ′ d )) .
This minimal cost of all homeomorphic alignments between G 1 and G 2 , called the homeomorphic

alignment distance, is denoted by η(G 1 , G 2 ).

Algorithm for rooted trees

In order to compute homeomorphic alignment distance with good complexity, we will use a bottom-up approach. This approach, widely used in the literature on edit-based distances, is based on a reformulation of the distance between trees in terms of both the distance between their subtrees and the matching of their roots. Then, using this reformulation, the computation starts by the distances between the simplest subtrees (i.e the leafs), and goes up to the root. The distances between subtrees are kept in memory, avoiding redundancy of computation. However, there are two main differences to take into account in the case of homeomorphic alignment distance:

1. Information is on arcs, rather than on vertices. It involves to take into consideration several arc weights in the reformulation (one for each arc between the root and its children), instead of the unique value linked to the root, as in the literature.

2. An arc can be matched with a set of arcs merged together (due to merging kernel application). It involves to take into account not only the subtrees of the tree, but also those which can be generated by some merging.

Our approach to solve these two problems is to use a special kind of tree, the root of which has only one child. Then, we make a reformulation of the homeomorphic alignment distance between this kind of trees, in function of subtrees with the same property. By this way, there is only one new arc to take into consideration at each step, the one between the root and its child. Moreover, we consider that this arc can be the result of a merging. In combination with the recursive nature of the bottom-up approach, it will lead to consider subtrees with all possible mergings. Finally, we also need a reformulation of the homeomorphic alignment distance between "classic" trees in function of the homeomorphic alignment distances between such subtrees.

Definitions and notations

First, we need to define more formally the particular kind of tree described above. Let T = (V, A, ω) be a weighted tree rooted in r T . We denote by T (v), v ∈ V , the subtree of T rooted in v (see Figure 9.7). Let v a be an ancestor of v, we denote by T (v, v a ) the subgraph of T obtained from T (v a ) by removing all complete subtrees which do not contain vertices of T (v), and by merging on each vertex n ∈ anc(v) \ {v a , v}. We say that T (v, v a ) is the subtree of T rooted in v a pruned in v, and we call this kind of trees a pruned tree (see Figure 9.7 for an example).

We denote by F(T, v) the forest, the connected components of which are the trees T (p, v), for all p ∈ C(v) (see Figure 9.7 for an example). By abuse of notation we also denote by F(T, v)

the set of all connected components of this forest.

Reformulations

In order to use a bottom-up approach, we have to express the homeomorphic alignment distance:

• between trees in function of the distances between their forests (Proposition 18). • between pruned trees in function of the distances between their pruned subtrees (Proposition 20).

• between forests in function of the distances between their subforests (Proposition 21).

In addition, we have to reformulate these distances in the special cases where at least one of the trees is empty (Proposition 19).

In the sequel we consider two weighted trees P = (V P , A P , ω P ) and D = (V D , A D , ω D ), rooted respectively in r P and r D .

Proposition 18. We have :

η(P, D) = η(F(P, r P ), F(D, r D )).

Proof. Since the root of a tree cannot be eliminated by any operation (merging, deletion) involved in the definition of η, it may be seen that any homeomorphic alignment H of P with D has a corresponding alignment H ′ of F(P, r P ) with F(D, r D ) of same cost, and the converse also holds.

Proposition 19. Let i ∈ V P \ {r P }, j ∈ V D \ {r D }, i a ∈ anc(i), j a ∈ anc(j), η(∅, ∅) = 0 η(P (i, i a ), ∅) = η(F(P, i), ∅) + γ(ω P (i a , i), 0) η(F(P, i), ∅) = i ′ ∈C(i) η(P (i ′ , i), ∅) η(∅, D(j, j a )) = η(∅, F(D, j)) + γ(0, ω D (j a , j)) η(∅, F(D, j)) = j ′ ∈C(j)
η(∅, D(j ′ , j)).

Proof. Straightforward.

Proposition 20. Let i ∈ V P \ {p}, j ∈ V D \ {d}, i a ∈ anc(i), j a ∈ anc(j). η(P (i, i a ), D(j, j a )) = min             
η(P (i, i a ), ∅) + η(∅, D(j, j a )), γ(ω P (i a , i), ω D (j a , j)) + η(F(P, i), F(D, j)), min jc∈C(j) {η(P (i, i a ), D(j c , j a )) + η(∅, F(D, j) \ D(j c , j))}, min ic∈C(i) {η(P (i c , i a ), D(j, j a )) + η(F(P, i) \ P (i c , i), ∅)}.

Proof. Let H be an homeomorphic alignment of P (i, i a ) with D(j, j a ), and let

H L = (V L , A L , ω L )
and

H R = (V R , A R , ω R
) the left and right graphs of H, respectively.

There are seven possible cases:

1. H is an empty set (it is cheaper to remove both P (i, i a ) and D(j, j a ) than to align them).

2. {(i a , i), (j a , j)} ∈ H.

3. ∃f ∈ A R , f being obtained by merging (j a , j) with other arcs. In this case, there is one and only one child j c of j, such (j, j c ) is merged with (j a , j) in f , and then all D(j ′ c , j), j ′ c ∈ C(j) \ {j c } are deleted. 4. ∃f ∈ A L , f being obtained by merging (i a , i) with other arcs. In this case, there is one and only one child i c of i, such (i, i c ) is merged with (i a , i) in f , and then all

P (i ′ c , i), i ′ c ∈ C(i) \ {i c } are deleted.
Cases 1,2,3,4 justify, respectively, the lines 1,2,3,4 of the expression of η(P (i, i a ), D(j, j a )) in the proposition.

The three last cases cannot lead to a better homeomorphic alignment: 5. The deletion of (i a , i) and (j a , j) cannot be preferred to the resizement (possible case 2), because γ(ω P (i a , i), 0) + γ(0, ω D (j a , j)) ≥ γ(ω P (i a , i), ω D (j a , j)).

6. If (i a , i) was deleted, and not (j a , j), then only one P (i c , i), i c ∈ C(i) is aligned with D(j, j a ), the other being removed. It is less expensive to merge (i a , i) with (i, i c ) (possible case 3), because γ(ω P (i a , i), 0) + γ(ω P (i, i c ), ω D (j a , j)) ≥ γ(ω P (i a , i c ), ω D (j a , j)) (Prop. 14).

7. The deletion of (j a , j) is more expensive than the merging of (j a , j) (possible case 4), for the same reasons as above.

Proposition 21. ∀A ⊆ F(P, i), B ⊆ F(D, j), η(A, B) = min                                  min D(j ′ ,j)∈B {η(A, B \ {D(j ′ , j)}) + η(∅, D(j ′ , j))}, min P (i ′ ,i)∈A {η(A \ {P (i ′ , i)}, B) + η(P (i ′ , i), ∅)}, min P (i ′ ,i)∈A,D(j ′ ,j)∈B {η(A \ {P (i ′ , i)}, B \ {D(j ′ , j)}) +η(P (i ′ , i), D(j ′ , j))}, min P (i ′ ,i)∈A,B ′ ⊆B {η(A \ {P (i ′ , i)}, B \ B ′ ), +η(F(P, i ′ ), B ′ ) + γ(ω P (i, i ′ ), 0)}, min A ′ ⊆A,D(j ′ ,j)∈B {η(A \ A ′ , B \ {D(j ′ , j)})+ η(A ′ , F(D, j ′ )j) + γ(0, ω D (j, j ′ ))}.
Proof. Let H be an homeomorphic alignment of A ⊆ F(P, i) with B ⊆ F(D, j), and let P (i ′ , i) ∈

A and D(j ′ , j) ∈ B. There are five possible cases:

1. D(j ′ , j) is not aligned with element of A,

2. P (i ′ , i) is not aligned with element of B, 3. P (i ′ , i) is aligned with D(j ′ , j), 4
. disjoint subparts of P (i ′ , i) are aligned with elements of B, 5. disjoint subparts of D(j ′ , j) are aligned with elements of A.

Cases 1, 2, 3, 4, 5 justify, respectively, the lines 1, 2, 3, 4, 5 of the expression of η(A, B) in the proposition.

Algorithm

Using the reformulations defined above, we can now design the Algorithm 14 following a bottomup approach. For a tree T = (V, A, ω), we define a suffix order for the vertices in V , such that for all v ∈ V , all the children of v are before v.

Complexity

We denote by N the maximum degree of a vertex in the pattern tree or in the data tree, H the maximal height of both trees and S the maximal size (number of vertices) of both trees. We recall that a set of size n has 2 n subsets.

The time complexity of computing each reformulation used in the algorithm is:

• η(P (i, i a ), ∅) and η(∅, D(j, j a )) → O(1).

• η(A, ∅) and η(∅, B) → O(N ).

• η(P (i, i a ), D(j, j a )) → O(N ).

• η(A, B) → O(N 2 N ).
Combining these complexities with the different loops of the algorithm, we obtain that the total

computation is in O(S 2 (N 2 3N + H 2 N )) time complexity.
If the maximal degree is bounded, the total computation is in O(S 2 H 2 ) time complexity.

Algorithm for unrooted trees

First, let us give an expression of the homeomorphic alignment distance between unrooted trees, in function of the distances between all their possible rooted versions.

Let G = (V, E, ω) be a weighted tree, let r ∈ V , we denote by G r , the directed weighted tree rooted in r, such that G is the undirected graph associated to G r (see Figure 9.8). Proposition 22. Let P = (V P , E P , ω P ) and D = (V D , E D , ω D ) be two weighted trees. We have:

η(P, D) = min i∈V P ,j∈V D {η(P i , D j )}.
Proof. Let G = (V, E, ω) be a graph, and r ∈ V a vertex of G. Notice that:

• a merging occurring on a vertex v ∈ V \ {r} in G can occur in G r ,
• deletion, insertion, resizement, and division occurring in G can occur in G r .

On the other hand, for each optimal homeomorphic alignment H of P in D, it is easy to see that there exists p ∈ V P and d ∈ V D , such that p and d are not affected by a merging. For example, if

D ′ = (V ′ D , E ′ D , ω ′ D )
is a subgraph such as η(P, D) = η(P, D ′ ), p and d can be chosen as 1-degree vertices of V P and V ′ D , respectively. As a result, η(P, D) = η(P p , D d ). Since the homeomorphic alignment is more constrained in the case of rooted trees than in the case of unrooted trees, we can assure than η(P, D) ≤ η(P a , D b ), a ∈ V P , b ∈ V D . To sum up, knowing that η(P, D) ≤ η(P a , D b ), a ∈ V P , b ∈ V D , and that there exists p ∈ V P and d ∈ V D , such that η(P, D) = η(P p , D d ), we conclude that η(P, D) = min i∈V P ,j∈V D η(P i , D j ).

Naive algorithm

Let P = (V P , E P , ω P ) and D = (V D , E D , ω D ) be two weighted trees. From Proposition 22, we propose a first, naive algorithm to compute η(P, D), consisting of computing the homeomorphic alignment distance for all couples of weighted rooted trees we can obtain from P and D, and keeping the minimum reached.

Complexity

Let P = (V P , E P , ω P ) and D = (V D , E D , ω D ) be two weighted trees. As the number of rooted trees we can obtain from an undirected tree is equal to the number of vertices of this graph, the number of couples of weighted rooted trees we can obtain from P and D is equal to

|V P | * |V D |.
The total computation is then in O(S 4 (N 2 3N + Ĥ2 N )) time complexity, with Ĥ representing the maximal height of a rooted tree obtained from P or D (N and S are already defined in Section 9.3.4).

If the maximal degree is bounded, the total computation is in O(S 4 Ĥ2 ) time complexity.

Optimized algorithm

It is easy to see that the above algorithm computes the alignment of subparts of P and D more than one time. Using dynamic programming and an adapted order of navigation in the tree, we can avoid useless computation.

Let P = (V P , E P , ω P ) and D = (V D , E D , ω D ) be two undirected weighted trees.

We denote by F(P, a, b), a, b ∈ V P the set of rooted trees P r , r ∈ V P , such that b is an ancestor of a in P r . We denote by anc(P, a, b), a, b ∈ V P the set of vertices x ∈ V P such that x is an ancestor of a in at least one rooted tree in F(P, a, b).

We denote by C(P, a, b), a, b ∈ V P the set of of vertices x ∈ V P such that x is a child of a in at least one rooted tree in F(P, a, b).

For computing η(P p (i, i a ), D d (j, j a )) and η(F(P p , i), F(D d , j)), we need to know η(P p (i c , i), D d (j c , j)) for all i c ∈ C(P, i, p), j c ∈ C(D, j, d).

We can start by computing η(P p (i, i a ), D d (j, j a )) and η(F(P p , i), F(D d , j)), for all i (respectively, j) being a leaf of P p (respectively, D d ), which have no child, by definition, and continue iteratively with all vertices which have all their children already computed.

Adapted order of navigation

An adapted order of navigation for a tree T = (V, E, ω) can be obtained by the Algorithm 15.

Algorithm 15: Order of navigation computation algorithm (computeOrder) Data: unrooted tree T Combining the time complexities of computing the reformulations (given in the Section 9.3.4)

Result: list of couples of vertices L begin list of couples of vertices L ← ∅ FIFO queue of vertices Q ← ∅ foreach v ∈ V do if deg(v) = 1 then push v in Q while Q = ∅ do pop v of Q foreach w ∈ N (v) do if (v, w) / ∈ L then if ∀x ∈ N (v) \ {w}, (x, v) ∈ L then push w in Q add (v,
with the loops of the algorithm, and taking into account the above observation, we obtain that the total computation time of this algorithm is in O(S 2 N (2 3N + S 2 )) complexity, with N and S defined as in Section 9.3.4.

If the maximal degree is bounded, the total computation is in O(S 4 ) time complexity.

Algorithm for rooted tree with unrooted tree

As seen in the previous section, the optimal homeomorphic alignment between two unrooted trees P and D is equal to the minimal optimal homeomorphic alignment between all possible rooted versions of P and D.

From this fact, we can define the optimal homeomorphic alignment between a tree P r rooted in r and an unrooted tree D as the minimal optimal homeomorphic alignment between P r and all possible rooted versions of D.

From Proposition 22 we can propose Proposition 23, which is straightforward.

Proposition 23. Let P r = (V P , A P , ω P ) a weighted tree rooted in r ∈ V P and

D = (V D , E D , ω D )
an unrooted weighted tree. We have:

η(P r , D) = min j∈V D {η(P r , D j )}.
The only difference between matching a rooted tree P r with an unrooted tree D and matching two unrooted trees P and D is that no merging can occur on the root r of P r , by definition. It is interesting if we want to be sure that a vertex of P is preserved by the matching.

Algorithm

The algorithm to compute the homeomorphic alignment between rooted tree and unrooted tree (see Algorithm 17) can be obtained by crossing Algorithm 14 and Algorithm 16. We denote by N the maximal degree of both trees, S their maximum size, and H the height of the rooted tree. The total computation time of this algorithm is in O(S 2 N (2 3N + SH))

complexity.

If the maximal degree is bounded, the total computation is in O(HS 3 ) time complexity.

Usage of Cut Operation

For the purpose of removing spurious branches without any cost, we propose to integrate the cut operation in our alignment.

In [START_REF] Wang | Finding approximate patterns in undirected acyclic graphs[END_REF], Wang et al. propose a new operation allowing to consider only a part of a tree. Let G = (V, A, ω) be a weighted tree. Cutting G at an arc a ∈ A, means removing a, thus dividing

G into two subtrees G 1 and G 2 .
The cut operation consists of cutting G at an arc a ∈ A, then considering only one of the two subtrees. Let K a subset of A. We use Cut(G, K, v) to denote the subtree of G containing v and resulting from cutting G at all arcs in K. In the case of a rooted tree, we consider that the root r G of G cannot be removed by the cut operation, and then we use the notation Cut(G, K) = Cut(G, K, r G ). In the case of a rooted forest, we consider that the root of each rooted tree composing the rooted forest cannot be removed by the cut operation, and then we use the same notation than above: Cut(G, K).

See Figure 9.9 for some examples of cut operation. Our main problem can be stated as follows: Given a weighted tree P = (V P , A P , ω P ) (the pattern tree) and a weighted tree G D = (V D , A D , ω D ) (the data tree), find

η cut (P, D) = min K⊆A D ,v∈V D {η(P, Cut(D, K, v)}
and the associated homeomorphic alignment. In the case of rooted trees and rooted forests, η cut (P, D) = min K⊆A D {η(P, Cut(D, K))}.

Integration of cut operation in our algorithm

It can be seen that the cut operation can be integrated in the homeomorphic alignment by replacing the deletion of complete subtrees by their cut.

In order to obtain the algorithms for the computation of η cut , we just need to replace:

• η by η cut in the reformulations and algorithms of the Sections 9.3 and 9.4

• Proposition 19 by the following:

Proposition 24. Let i ∈ V P \ {p}, j ∈ V D \ {d}, i a ∈ anc(i), j a ∈ anc(j), η cut (∅, ∅) = 0 η cut (P (i, i a ), ∅) = η cut (F(P, i), ∅) + γ(ω P (i a , i), 0) η cut (F(P, i a ), ∅) = i ′ ∈C(ia) η cut (P (i ′ , i a ), ∅) η cut (∅, D(j, j a )) = 0 η cut (∅, F(D, j a )) = 0.
It is interesting to remark that the complexity is not modified.

Limitations

The homeomorphic alignment is an efficient tool for our purpose of motion capture: it affords to match our model tree with data tree, taking into account the possible noises occurring in our context. However, we can observe some limitations.

The merging operation have to be more restricted. For our application, merging in the model tree are useless, the problem of useless 2-degree vertices occurring only in the data tree. Furthermore, merging operations in model tree can lead to bad matchings, as shown in Figure 9.10.

Restricting the merging operation on the data tree could lead to both optimize the computation speed and increase the accuracy of results. Example of bad homeomorphic alignment matching due to merging on model tree: the "thumb" edge is removed for a cost of 6, as the merging of "arm" edge and "palm" edge in both trees reduce the matching cost by 8.

Chapter 10

Asymmetric Homeomorphic

Alignment

In the previous chapter, we proposed a first solution to our problem of matching by the way of homeomorphic alignment. We have closed this chapter by showing the limitations of this method, due to non restrictions on merging operation.

The alignment and the homeomorphic alignment allow to find the best way to modify two trees in order to reach a specific relation between them (isomorphism in the first case, homeomorphism in the second). Isomorphism is a sub case of homeomorphism (i.e. two trees isomorphic are homeomorphic) and by extension, an alignment between two trees is an homeomorphic alignment as well.

In this chapter, we first define an asymmetric relation being sub case of homeomorphism and for which isomorphism is a sub case. Then, we propose new alignment allowing to find the best way to modify two trees in order to reach this relations. Finally, we design algorithms in order to compute it efficiently. • a R a (reflexivity)

• if a R b then b R a (symmetry) • if a R b and b R c then a R b (transitivity)
The equivalence class of an element s of S by R is defined by R(s) = {x ∈ S|xRs}.

The binary relation R is an order relation if and only if for all a, b and c in S:

• a R a (reflexivity)

• if a

Asymmetric Homeomorphic Alignment

Let G 1 = (V 1 , A 1 , ω 1 ) and G 2 = (V 2 , A 2 , ω 2 ) be two weighted graphs. Let G ′ 1 = (V ′ 1 , A ′ 1 , ω ′ 1 ) and G ′ 2 = (V ′ 2 , A ′ 2 , ω ′
2 ) be weighted graphs obtained by deleting arcs in G 1 and G 2 , such that

G ′ 2 H ≥ G ′ 1 . Let G ′′ 2 = (V ′′ 2 , A ′′ 2 , ω ′′
2 ) be a graph obtained from G ′ 2 by a sequence of merging, such there exists an isomorphism I between G ′ 1 and G ′′ 2 .

The set of all couples of arcs H = {(a, a ′ );

a ∈ A ′ 1 , a ′ ∈ A ′′ 2 , a ′ = I(a)} is called an asymmetric homeomorphic alignment of G 1 with G 2 .
The cost C H of H is the sum of the costs of all operations used to asymmetric homeomorphically align G 1 and G 2 : the deletion of arcs in G 1 and G 2 , to obtain G ′ 1 and G ′ 2 respectively, and the resizement for each arc a 1 ∈ A ′ 1 to the weight of H(a 1 ). More formally :

C H = (a,a ′ )∈ H γ(ω ′ 1 (a), ω ′′ 2 (a ′ )) + a d ∈A 1 \A ′ 1 γ(ω 1 (a d ), 0) + a ′ d ∈A 2 \A ′ 2 γ(0, ω 2 (a ′ d )) .
This minimal cost of all homeomorphic alignments between G 1 and G 2 , called the asymmetric homeomorphic alignment distance, is denoted by η(G 1 , G 2 ).

Algorithm for Rooted Trees

The algorithms for compute the asymmetric homeomorphic alignment distance is designed in a similar way that those we described for compute the homeomorphic alignment distance: we use again a bottom-up approach, but here we consider only mergings on data tree. Due to this fact, we consider pruned trees only on the data tree side.

Reformulations

In the sequel we consider two weighted trees P = (V P , A P , ω P ) and D = (V D , A D , ω D ), rooted respectively in r P and r D .

As merging is not used in the reformulations 25, 26 and 28, the proofs are the same as for Propositions 18, 19 and 21, respectively.

Proposition 25. We have :

η(P, D) = η(F(P, r P ), F(D, r D )). Proposition 26. Let i ∈ V P \ {r P }, j ∈ V D \ {r D }, i a = par(i), j a ∈ anc(j), η(∅, ∅) = 0 η(P (i, i a ), ∅) = η(F(P, i), ∅) + γ(ω P (i a , i), 0) η(F(P, i a ), ∅) = i ′ ∈C(ia) η(P (i ′ , i a ), ∅) η(∅, D(j, j a )) = η(∅, F(D, j)) + γ(0, ω D (j a , j)) η(∅, F(D, j a )) = j ′ ∈C(ja)
η(∅, D(j ′ , j a )).

Proposition 27. Let i ∈ V P \ {r P }, j ∈ V D \ {r D }, i a = par(i), j a ∈ anc(j).

η(P (i, i a ), D(j, j a )) = min

            
η(P (i, i a ), ∅) + η(∅, D(j, j a )), η(γ(ω P (i a , i), ω D (j a , j)) + η(F(P, i), F(D, j)), min jc∈C(j) { η(P (i, i a ), D(j c , j a )) + η(∅, F(D, j) \ D(j c , j))}, η(γ(ω P (i a , i), 0) + min ic∈C(i) { η(P (i c , i), D(j, j a )) + η(F(P, i) \ P (i c , i), ∅)}.

Proof. Let H be an asymmetric homeomorphic alignment of P (i, i a ) with D(j, j a ), and let

H L = (V L , A L , ω L ) and H R = (V R , A R , ω R
) the left and right graphs of H, respectively.

There are six possible cases:

1. H is an empty set (it is cheaper to remove both P (i, i a ) and D(j, j a ) than to align them).

2. {(i a , i), (j a , j)} ∈ H.

10.4 Algorithm for Unrooted Trees and for Rooted Tree with

Unrooted Tree

Algorithms for unrooted trees and for rooted tree with unrooted tree can be obtained using the same method that for homeomorphic alignment. In case of asymmetric homeomorphic alignment, the algorithm for unrooted trees is in O(S 2 * (N * 2 3 * N + S D * N )) time complexity (O(S 3 ) if the maximal degree is bounded).

The algorithm for rooted tree with unrooted tree is in O(S 2 * (N * 2 Chapter 11

Alignments Comparison

All the experiments in this section have been done on a computer with a processor Intel(R) Q8200 at 2.33 GHz and 3 Gio of ram.

Complexities Comparison

In our method, we assume that both pattern tree and data tree have a bounded maximal degree. We denote by S P and S D the size of the pattern tree and the size of the data tree, respectively.

We denote by H P and H D the height of the pattern tree and the size of the data tree, respectively.

Accuracy

In this section we propose to study the robustness of the different alignments in regard of the noises occurring in our motion capture framework.

Protocol

Each experiment consists of the following:

1. We randomly generate a pattern tree P .

2. We generate a data tree D from P in two steps:

weight variation: we randomly alter the edge weights by a given amount of variation. 3. We compute the different alignments between P and D.

4. The precision is given by the percentage of good matchings, that is the percentage of pattern tree vertices that match with their equivalent in the data tree (which is known, according to the above protocol to generate the data tree).

Results

The results, obtained by averaging precisions over 1000 experiments, are shown on Figures 11.3,11.4 and 11.5.

Discussion

The first observation is that the homeomorphic alignment and asymmetric homeomorphic alignment provide a perfect matching if no weight variation occurs, even for large amount of structural noise.

On the other side, for high weight variation and few structural noises, we can observe that the percentage of good matchings for homeomorphic alignment is lower than those for classic alignment. It is due to the use of mergings on both trees, as shown on Figure 11.6: as the vertex d disappear by merging in case of homeomorphic alignment, the percentage of good matchings decrease for this alignment. The asymmetric homeomorphic alignment allow the merging only on data tree, this special case is treated in the same way than for classic alignment.

Concerning the cut operation, we can observe that using it for small weight variations increase the percentage of good matching, but for large weight variations, the results are worst than without cut operation. Then the choice of using cut operation or not will depend to the a priori knowledge of amount of weight variation and structural noise in the data tree and to the presence or not of others mechanisms increasing the robustness. Finally, it is interesting to remark than in all the cases, the asymmetric homeomorphic alignment provide more good matchings than the other alignments, and then is the most adapted in our context of motion capture.

Computation Speed

In this section we propose to study computation speed of the different alignments in a general purpose.

Protocol

Each experiment consists of the following:

1. We randomly generate n trees T 0 , T 1 , ..., T n-1 2. For each tested algorithm:

• we record the starting time t s

• we apply the algorithm on each possible couple (T i , T j )

• we record the ending time t e

• the average speed of the tested algorithm is equal to (t et s ) n 2

Results

The results, obtained for n = 32, are shown on Figure 11.7.

Discussion

The speed difference of the algorithms is well correlated with their complexities. However, we can observe a variation of speed between the classic alignment for rooted trees and the classic alignment for unrooted trees, even if they have the same complexity. This variation is due to the maximal degree observed: in the rooted tree, the degree considered is in fact the number of children. Assuming that the maximal number of children is equal to the maximal degree decreased by one, the constant linked to the computation of forests alignments is divide per two for the alignment between rooted trees.

The constant C linked to the computation of forests alignments (due to bounded maximal degree) have to be take into consideration when we observe the computation speed. Table 11

Frequencies in Motion Capture Context

In this section we propose to study frequencies (number of iterations per second) of the alignments applied in our motion capture context. We first give some statistics on trees for different subjects and different sizes of grid, then, we show the frequencies results in regard of these statistics.

Statistics on data trees

We propose to study two models: the full human body model and the hand model. Full human body model tree have seven vertices and a maximal degree equal to four. Hand model tree have heights vertices and a maximal degree equal to five.

For the data, we use here the dancer data set (height views) from 4Drepository from INRIA for full human body application, and a home made data set with four views for hand application.

We use these sets on voxels grids with side 40,80,120. We show the results of distributions of tree sizes in Figures 11.8,11.9,11.10 and 11.11

Choice of algorithm

First, we have to determine which alignment is the most adapted: as we want preserve all vertices in the model tree, and take into consideration useless 2-degree vertices in the data tree, we can use the asymmetric homeomorphic alignment in our application.

In a second time, we have to specify if both the model tree and the data tree can be rooted or not. If the both trees are rooted, it involves that the roots will be matched together. If only one tree is rooted, it involves that no merging can occurs on the root, then the root will already be matched.

In the case of full human body, and for all complete models, the only assumption we can do is that all vertices have to be matched in the model, but their is no particular vertex in the data tree. We can use a rooted version of the model tree, with an unrooted version of the data tree.

The choice of the root have consequences on the speed results: if the root of full human body model tree is the torso, the rooted tree have a max number of children equals to four, and only three if the root is an other vertex.

In the case of hand, and for all incomplete models, we can consider that the 1-degree vertex which is an extremity of the infinite weighted edge in the model tree will always match with the vertex associated with the grid border point (supposed unique) in the data tree. Considering this fact, both trees can be rooted on the vertices described above.

Protocol

Each experiment consists of the following:

1. We randomly generate n data trees 2. we record the starting time t s 3. we compute the asymmetric homeomorphic alignment between the model tree and each data tree 4. we record the ending time t e 5. the frequency is equal to n (t et s )

Results

The results, obtained for n = 10000, are shown on Figure 11.12 for hand application, and on

Discussion

The average frequency can be done by crossing the graphs statistics and the frequencies using the following formula: Then, in the average case:

• for the dancer data set, we can perform more than 708 alignments per second

• for the hand data set, we can perform more than 548 alignments per second

We can consider that the matching step of our method will be performed very fast.

In the worst case:

• for the dancer data set (data tree with 17 vertices and a degree max equal to 5), we can perform more than 89 alignments per second

• for the hand data set (data tree with 24 vertices and a degree max equal to 5), we can perform more than 46 alignments per second It involve that even if the worst case occurs all the time, the matching with model tree can be performed in real time.

Of course, as the size and the maximal degree of data trees is correlated with the size of the grid, we can choose this size in order to obtain a good compromise between accuracy of positions and computation speed.

Part IV

Motion Capture Applications

Chapter 12

Generic Pose Initialization Step

In the two previous parts of this manuscript, we proposed solutions for both fast skeletonization and adapted matching.

In this chapter, we first describe in details the pipeline of our method, using the works described in the previous parts. Then, in order to improve the robustness of our method, we propose the optional addition of intuitive constraints in the model definition.

Summary of our Method

Now, as we have found efficient solutions for skeletonization and thinning, we can integrate them in our pipeline. For the skeletonization, the use of our D6I1D thinning algorithm (described in Chapter 6) is interesting, in regard of its speed and it resulting skeleton, which has few spurious branches. In the case of incomplete model, we constrain the thinning process to preserve the centroid of the set of points on the boundaries of the grid, belonging to the shape. For the data tree extraction from the skeleton, we use the Algorithm 13 (defined in Section 7.3).

Use of Asymmetric Homeomorphic Alignment

The matching of the model tree with the data tree is done by asymmetric homeomorphic alignment (described in Chapter 10), specially designed for this purpose. Three different algorithms have been proposed for the computation of asymmetric homeomorphic alignment, depending on the kind of trees to align: a) two unrooted trees: this algorithm can be used for matching model tree and data tree without constraints. 141 b) one rooted tree and one unrooted tree: applying this algorithm involves that the root of the rooted tree matches with a vertex of the unrooted tree, and cannot be removed by a merging. c) two rooted trees: applying this algorithm involves that the two roots match together.

When applied to the same data, algorithm c) is faster than algorithm b), which is itself faster than algorithm a). In order to use the full potential of the proposed algorithms, we have to consider different cases, in regard of the class of model.

In case of incomplete model, it is easy to see that infinity weighted edges in both model tree and data tree have to match. Rooting the trees in the two unique 1-degree vertices bounding these edges involves matching them together, and by this way, to match together the unique edges linked to them.

In case of complete model, as we are mostly interested in matchings that preserve a large part of the model tree vertices, we can root the model in one of the vertices which have less chances to be removed, a high degree vertex for example. As we have no information about the importance of data tree vertices, this tree is considered as unrooted.

For the purpose of scale invariance, we start by normalizing the weights of the two trees, so that the sum of all non-infinite weights in a tree is equal to 1. Then, we compute the asymmetric homeomorphic alignment (and the associated distance) between the trees. If the distance is greater than a given threshold T d , which is defined by the user, we assume that the data tree is not enough similar to the model tree to provide a good matching, and we stop the pipeline for this frame. Otherwise, we use the optimal asymmetric homeomorphic alignment to match the labels associated with the model vertices, with the 3D positions associated to the data vertices.

The number of non-aborted matchings, and their quality, obviously depend on the choice of the threshold value T d : as both trees are normalized, T d = 2 means that no matching will be aborted (both trees can be deleted), but the resulting matching can be poor. On the other hand, a lower value yields a greater amount of aborted matchings, but with better probability of good matching.

Pipeline

A detailed view of the pipeline of our method is shown in Figure 12.1.

The second special case occurs when different limbs of a model are represented by edges with similar lengths in data tree. In the case of the hand for example, we can consider the case of four fingers: index, middle finger, ring finger and little finger.

In the reality, little finger is always (or at most in a wide majority of case) smaller than the three others, and the middle finger is longer than the other fingers. It seems natural to represent them in this way in the model, by giving a higher value to the edge representing the middle finger, and a smaller to the little finger, in order to differentiate them.

However, the branches of the skeleton representing these fingers can have very variating lengths, due to the visual hull reconstruction, for example. The branch representing little finger can be longer than the one representing the middle finger, as shown in Figure 12.3. To solve this problem, we have to represent the four fingers by edges with same length, and to find an other way to differentiate them. In order to solve this kind of problem, we propose to add optional constraints in our method.

These constraints have to be intuitive and easy to describe, in order to preserve the philosophy of our method ("keep it simple"). We propose two different kinds of constraints, both requesting a specific positioning of a subject part relative to others.

A coordinate constraint corresponds to constraining the position of a subject part in regard of another linked part, for a given axis (for example, constrain the head to be above the torso).

A betweenness constraint corresponds to constraining the position of a subject part in regard of two others, requesting to be between them (for example, constrain the middle finger to be between the index and the ring finger).

In the sequel of this section, we propose a way to integrate these optional constraints in the model definition and in the pipeline of the method.

Coordinate Constraints

A coordinate constraint consists of forcing a part of the model to be positioned in a subspace defined by one of its linked parts. For example, in the case of full human body model, we can constrain the head to be above the torso: let P H = (x H , y H , z H ) be the 3d position of the head and P T = (x T , y T , z T ) the one of the torso, the constraint is

P H ∈ {P = (x, y, z) ∈ R 3 ; y > y T }.
Considering any rooted version of both the model tree and the data tree, and using the property that coordinate constraint is defined on linked parts only, we can translate coordinate constraint into a constraint on the tree matching.

Let v m and w m be the two vertices associated to the constraint in the model tree. For example, for the head-torso constraint, v m (respectively, w m ) represents the torso (respectively, the head).

By definition of the constraint (the two parts have to be linked), there exists an arc a m = (v m , w m ) belonging to the model tree. Let v be any vertex of the data tree, we denote by

P (v) = (P x (v), P y (v), P z (v))

Betweenness Constraints

In case of similar limbs, the matching with the model tree can generate multiple solutions. In the case of our model of hand, as the descriptions of index, middle finger, ring finger and little finger are identical, the matching will give a set of four possible positions for each finger. To solve this problem, we introduce the betweenness constraints. It consists of constraining the position of some limb to lie between two other ones.

An usual definition of a betweenness ternary relation is based on collinearity [START_REF] Greenburg | Euclidean and non-Euclidean geometries: Development and history[END_REF] : a point P 3 is said to be between two other points P 1 and P 2 (denoted by R c (P 1 , P 2 , P 3 )) if it belongs to the segment P 1 P 2 . However, this definition is too restrictive.

On the other hand, we could say that a point P 3 is between P 1 and P 2 (denoted by R p (P 1 , P 2 , P 3 ))

if its orthogonal projection on the line (P 1 P 2 ) is between P 1 and P 2 . In this case, the problem is that there exist triplets (P 1 , P 2 , P 3 ) such that each point is between the two other ones (for example, the three vertices of an equilateral triangle).

These two relations are illustrated in Figure 12.4.

Design of Intuitive Betweenness Relation

We propose a definition which is not too restrictive, and which gives at most one possibility of "betweenness" for three points P 1 , P 2 , P 3 : consider B the unique ball with diameter P 1 P 2 minus the unique sphere with diameter P 1 P 2 . The point P 3 is between P 1 and P 2 (denoted by R b (P 1 , P 2 , P 3 )) iff P 3 ∈ B (see Figure 12.5 for an illustration).

An other formulation is that P 3 is between two points P 1 and P 2 iff the angle between ---→ P 3 P 1 and ---→ P 3 P 2 is greater than π 2 . This definition is easy to use in our context, and provides a unique result for a wide set of configurations (all the configurations where the 3 points are the vertices of obtuse or right triangles). However, in order to cover all possible configurations, we have to consider other definitions.

Design of Mathematically Efficient Betweenness Relation

In order to design a betweenness relation on 3D points adapted to our purpose, we first need to define its required properties. A well adapted set of properties in our case can be based on the postulates of betweenness proposed by Huntington and Kline in 1917 [START_REF] Huntington | Sets of independent postulates for betweenness[END_REF] and refined by Huntington [START_REF] Huntington | A new set of postulates for betweenness, with proof of complete independence[END_REF].

relation is only satisfied on distinct points. Finally, Postulate 9 is required to derive an order on S from the ternary relation. In our case, it is not necessary.

These postulates have been proposed in a framework considering only 1D sets, like numbers or points of a line. In higher dimensions, some configurations are unsolvable: if three points are equidistant, it is impossible to find a relation such that one and only one of the points is between the others. We propose to modify the postulates to take this into consideration. To sum up, we have to design a ternary relation B such the system (R 3 , B) is a near-betweenness system satisfying the postulates A,B',C and D'.

However, we have to notice that these postulates are necessary but not sufficient for our context.

In order to show it, we propose a betweenness relation satisfying the four postulates but which is unusable in our case. However, even if this relation satisfies the four postulates, it is not usable in our case, due to the fact that it is not intuitive for the common conception of betweenness (see Figure 12.6 for some examples). Due to this fact, it is very difficult to propose betweenness constraints using this relation.

Unusable Betweenness Relation

As our model description has to be simple, we have to find a betweenness relation which satisfies the four postulates and which is intuitive enough for an easy usage. We can notice that the green cross is between the black ones, in the common conception of betweenness, but is not in the sense of R u . In the same way, the red cross is not between the black ones, in the common conception of betweenness, but is so in the sense of R u . We propose now a definition of intuitive betweenness relation, in the sense of the above definition.

More Intuitive Betweenness Relation

The non-equilateral triangles can be classified in two categories: isosceles triangles and nonisosceles triangles.

In the case of isosceles triangles, one of the three edges has a length different to those of the two others, by definition. The vertex at the opposite of this edge can be distinguished from the two others by this way. It is interesting to notice that the orthogonal projection of this vertex on its opposite edge is always between its boundaries.

In the case of non isosceles triangles, one of the three edges is greater than the two others. This way, the vertex at the opposite of this edge can be distinguished from the two others. It is interesting to notice that the orthogonal projection of this vertex on its opposite edge is always between its boundaries.

Application in our Method

In our model, a betweenness constraint is represented by a sequence of at least three model tree vertices. For the example of hand fingers, the constraint is formulated as "betweenness: thumb, index, middleFinger, ringFinger, littleFinger".

The application of this constraint is done on the output of the tree matching step, which provides multiple matching possibilities for some vertices. Then, if there exists a betweenness constraint on model vertices, defined by a sequence m 0 , ..., m n , the corresponding data vertices d 0 , ..., d n , with associated positions P 0 , ..., P n must be chosen such as, for a given betweenness relation R,

for each i ∈ [1, n -1], R(P i-1 , P i , P i+1 ).
The choice of the best betweenness relation will be empirically discussed in the next chapter.

Complete Method and Model Definition

Now we have solutions for problematic cases, thanks to the new optional constraints. The final pipeline of our method, including the management of these constraints, is shown in Figure 12.8.

In both cases, the program computes our pose initialization method in the following loop that consists of four steps:

1. loading of images and silhouette extraction (or direct silhouette loading when available)

2. visual hull computation 3. pose initialization step application

display of results

Several options can be set using command line arguments, e.g. allowing to only perform the two or three first steps.

Used Data Sets

For our tests, we used several data sets provided by INRIA Perception Group 1 :

• for the application to full human body, we used the dancer data set (see Figure 13.1a),

showing a professional dancer wearing a dress doing some moves. This data set is composed of the synchronized images from 8 calibrated cameras and also provides the white on black silhouette of the subject for each image. This data set is interesting for several reasons:

-as the subject wears a dress, it allows us to observe the impact of length variation of observed limbs;

-the subject does different poses, not conventional for usual initialization steps, but with all limbs separated from the others, allowing the study of the robustness of our method;

-the dress is attached to a floating piece of fabric, which is difficult to differentiate from a limb in visual hull. It is interesting to observe the robustness of our method in regard of this additional "noise".

• also for the application to full human body, we used the children data set (see Figure 13.1b),

showing a young boy performing some cartwheels. This data set is composed of the synchronized images (and their silhouettes) from 16 calibrated cameras (only 8 are used in our program). This data set is interesting as it forbids the use of coordinate constraints (the torso is above the head in a large amount of frames).

• for the application to animals, we used the dog data set (see Figure 13.1c), showing a dog trying to catch a ball held by a man. This data set is composed of synchronized images from 16 calibrated cameras (only 8 are used in our program). This data set is interesting for showing the trans classes adaptability of our method. However, the automatic silhouette extraction is very difficult on this data set, due to the fact that dog legs have the same color as the shadows. For our needs, we did manually extract silhouettes for a small set of frames.

In addition, we used a home made hand data set (see Figure 13.1d) for hand pose estimation.

This data set is composed of the synchronized images from 4 calibrated cameras. As camera positions and silhouette extractions are not performed with a good quality, the visual hull is very noisy: fingers are sometimes cut. However, this data set is interesting for the observation of the robustness in case of bad quality 3D reconstruction. The visual hull is computed for a cubic voxel grid with parameterizable resolution positioned in such a way that the whole set of actions performed by the subject is covered as tightly as possible.

13.2 Results

Speed

In order to measure the computation speed of our method, we use the following methodology:

for each data set (dancer, hand and children) and for each resolution of voxel grid, we measure the computation times of the application performing only its two first steps (image loading and visual hull computation) and its three first steps (the third being initialization step application).

We call these two computation times CT 2 and CT 3, respectively. Let N b f rames be the number of frames of the current data set. The computation time of our method is thus obtained using the following formula: CT 3 -CT 2 N b f rames . The average frame rate of our method can be obtained by inversing the previous formula. interesting to remark that method has a linear time complexity in regard of the voxel grid size (see Figure 13.2 top right plot).

Proportion of False Positives

A critical point in initialization step methods is the proportion of false positives, i.e. the proportion of invalid output poses. In order to appreciate the quality of an output pose for our method, as no ground truth data is available, we propose to classify them in regard of their usability, which can be visually measured:

1. an output pose is a good matching (see Figure 13.3a We propose to study the quality of the output pose in regard of the Asymmetric Homeomorphic Alignment (AHA) distance from the model tree to the data tree, in order to highlight the correlation between them.

For this purpose, we measure both the AHA distance and the quality of output pose (manually estimated) for each frame of different data sets, and for different grid resolutions. From these measures, we compute, for each AHA distance threshold, three values:

• the ratio of the number of true positives over the number of frames of the data set;

• the ratio of the number of false positives over the number of frames of the data set;

• the ratio of the number of false positives over the number of matchings (true positives plus false positives) .

The study of these data provides a lot of interesting information, as the detection of an efficient AHA distance threshold, the impact of grid resolution on the validity of the results, or the choice of the most adapted betweenness relation.

Study of Full Human Body Matchings

Results for dancer data set and children data set are shown in Figure 13.4 and Figure 13.5, respectively. We remind that the main difference between these data sets is that no coordinate constraint can be set in the case of children data set, due to the performed action (cartwheel).

Due to this fact, we obviously expect better results for the dancer data set. First, we observe from measures for both data sets that a higher grid resolution provides better quality (more true positives, missing limbs matching being replaced by good matchings). However, this improvement of the quality is more sensible from 50 3 resolution to 100 3 resolution than from 100 3 resolution to 150 3 resolution. Thus, we can assume that for all grid resolutions precise enough for a good visual hull definition, where each limb is clearly separated from others, and allowing a good thinning, the output quality will be about the same. In the same way, we can observe that false positives are fewer for higher grid resolutions.

Concerning the dancer data set, the study of ratios shows that very few false positives are associated to AHA distance lower than 0.4. The use of our method with AHA distance threshold set at this value provides a good matching for between 70 and 90 percent of the frames, with a fail rate lower than five percent.

Concerning the children data set, the study of ratios shows that for 50 3 grid resolution, no AHA distance threshold can be found in order to provide acceptable amount of matching with low rate of false positives. However, for higher resolutions, setting threshold to 0.2 provides a good matching for between 15 and 20 percent of the frames, with a fail rate lower than five percent.

The impossibility to use our method for 50 3 grid resolution for children data set is partially due to the fact that grid volume is higher for children data set than for dancer data set, as the subject is more mobile. A solution should be to refocus the grid on the subject at each frame.

Study of Hand Matchings

The main goal of our study of false positive for hand pose estimation is to find the most adapted Concerning the model definition, we encourage the use of constraints whenever possible, in order to obtain higher rate of true positives.

Possible Usages of our Method

Our method can be used in many different ways. In this section, we propose some of the most natural of them.

Combination with Tracking

The more natural usage is obviously to associate it with a tracking method, in order to obtain a full motion capture application.

Combination with Model Free Tracking

An easy way to use our method with tracking method is to associate it with model-free tracking, using for example set of points matching, like those proposed by Mateus et al. [START_REF] Mateus | Articulated shape matching using locally linear embedding and orthogonal alignment[END_REF][START_REF] Mateus | Articulated shape matching using laplacian eigenfunctions and unsupervised point registration[END_REF]. As those methods do not require any kind of a priori model, the goal of our method will be to perform a body part segmentation of the visual hull.

A simple method to obtain an efficient part labeling of each voxel from our method outputs (i.e.

skeleton with matched specific points) consists of two steps:

• color labeling each part of the skeleton bounded by matched points (e.g. voxels of the skeleton branch linking head and torso are colored in red, those of the branch linking the torso and the crotch are colored in cyan, etc).

• flooding the visual hull simultaneously with all the color labeled voxels.

Results obtained with this method are shown in Figure 13.10.

Combination with Model-Based Tracking

Model-based tracking methods found in literature requires specific models, e.g. kinematic models or geometric shapes, more complex than ours. However, we can use information provided by our method, in order to initialize them easily.

In the case of kinematic model, our method's output can be used as a guide to provide both position and scale for the kinematic model. A way to do it is to set positions of kinematic model extremities to those of our output, then to perform a minimization step in order to find the pose as close as possible to the skeleton. In the case of geometric shape-based model, we can use the above method, then use a fitting process in order to find parameters for each geometric shape, being the most adapted to the subject.

We can also use our method in order to provide color information on the different parts of the subject, e.g. using the segmentation method proposed in previous section, and looking for the color associated to these parts by projecting them in the input images. This color information can be used to easily track the motion of the subject, as proposed in [START_REF] Caillette | Real-time markerless human body tracking with multi-view 3d voxel reconstruction[END_REF][START_REF] Cheung | Shape-from-silhouette of articulated objects and its use for human body kinematics estimation and motion capture[END_REF][START_REF] Kehl | Full body tracking from multiple views using stochastic sampling[END_REF].

Specific Pose Detection and Pose Clustering

Another possible way to use our method is to detect some specific poses, or to cluster similar poses.

Specific Pose Detection

For a pose detection purpose, we can describe some poses combining the syntax used to describe constraints and boolean operators, e.g.:

• for a "hands up" pose, it will consist in searching if both hands are above the head, and the description is "(zgreatest A 1 -H) AN D (zgreatest A 2 -H)".

• for an "ask a question" pose, it will consist in searching if one and only one hand is above the head, and the description is "(zgreatest A 1 -H) XOR (zgreatest A 2 -H)".

• for a "wide spread arms" pose, it will consist in searching if head, torso and crotch are between the two arms, and the description is "(order A 1 -H -A 2 ) AN D (order A 1 -T -A 2 ) AN D (order A 1 -C -A 2 )".

• for an "acrobatic" pose, it will consist in searching if at least one of the feet is above the head, and the description is "(zgreatest F 1 -H) OR (zgreatest F 2 -H)".

Thus, the current pose is detected as correct if the positions of the parts fit with the description.

The good point of this method is the description syntax, which allow the detection of very specific poses, as well as pose of only some parts of the subject.

Pose Clustering

For a pose clustering purpose, consisting in regrouping frames where subject is in similar poses, our method can be used as follows: for each frame where an output pose is provided, we compute the set of all betweenness relations between three parts of the subject. Two poses are considered similar if their betweenness relation sets are similar enough (it can be measured using the cardinality of their intersection over the the cardinality of their union).

The two applications described above are interesting as they propose a high level description of the required pose, via understandable syntax (in case of specific pose detection) or an example (in case of pose clustering).

Chapter 14

Conclusion

In this thesis, we have provided an efficient real-time generic pose estimation step for multi-view marker-free motion capture, with multiple possible applications.

Contributions

Our contributions have been done in different fields of research:

Digital Topology. In this field, our contributions are:

• An efficient method to optimize the speed of critical kernel-based thinning algorithms, by providing a reformulation of the detection of crucial points, allowing the use of look-up tables.

• An easy-to-implement method to optimize the speed of any thinning algorithm using look-up tables, by reducing the time requested to compute the configuration of the neighborhood of a point.

• A new 6-directional thinning scheme based on isthmuses, speed efficient and providing good quality results.

We submitted a paper about this method at the 16th IAPR International Conference on Discrete Geometry for Computer Imagery.

Graph Theory. In this field, our contributions are:

• The Homeomorphic Alignment, a new alignment between weighted graphs, based on homeomorphism notion, instead of isomorphism one. We also proposed efficient algorithms computing it for weighted trees, rooted as well as unrooted.

We presented a paper about Homeomorphic Alignment at the 7th IAPR Workshop on Graph-based Representations in Pattern Recognition [START_REF] Raynal | Homeomorphic alignment of edge-weighted trees[END_REF], and an extended version has been published in Pattern Recognition [START_REF] Raynal | Homeomorphic alignment of weighted trees[END_REF].

169

• The Asymmetric Homeomorphism, a new binary relation on homeomorphic graph sets.

• The Asymmetric Homeomorphic Alignment, a new alignment between weighted graphs, based on our new binary relation. We also proposed efficient algorithms computing it for weighted trees, rooted as well as unrooted.

Motion Capture. In this field, our contributions are:

• The introduction of new constraints in a priori model definitions, increasing the robustness of motion capture methods.

• The formalization of betweenness for 3D points, and the proposition and evaluation of several definitions of betweenness relations.

• The definition of a new method for pose estimation step for multi-view marker-free motion capture, usable in real-time context, adaptable to different kinds of subject, and not constraining for the user.

We presented a paper about this method at the International Conference on Image Analysis and Recognition 2010 [START_REF] Raynal | Generic initialization for motion capture from 3d shape[END_REF].

Papers about our unpublished contributions are in redaction.

Future Works

The contributions of this thesis can be extended in different ways.

First, concerning motion capture applications, we want now focus on tracking methods, by providing new solutions using our knowledge in digital topology and mathematical morphology.

It can also be interesting to develop new methods for pose detection and pose clustering using our method.

Concerning the digital topology field of research, our aim is now to provide new thinning schemes adapted to CPU-parallelization, and to search for new parameterizable constraints more powerful than isthmuses.

Finally, concerning our graph theory contributions, we will try to design efficient algorithms for the computation of both Homeomorphic Alignment and Asymmetric Homeomorphic Alignment between any kinds of graph.

2. 1

 1 Examples of systems using markers . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Examples of model definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Illustration of visual hull concept . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Description of two different subjects . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Example of skeletonization process . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Examples of data tree extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Pipeline of our method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Example of skeleton in continuous framework . . . . . . . . . . . . . . . . . . . . 4.2 Example of medial axis in discrete framework . . . . . . . . . . . . . . . . . . . . 4.3 Illustration of 2D neighborhood . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Illustration of 3D neighborhood . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Illustration of a 3D hole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Examples of the different kinds of skeleton . . . . . . . . . . . . . . . . . . . . . . 4.7 Examples of 1D isthmuses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Illustration of 2D directions used in directional thinning . . . . . . . . . . . . . . 4.9 Illustration of 3D directions used in directional thinning . . . . . . . . . . . . . . 4.10 Illustration of subfields used in subfield-based thinning . . . . . . . . . . . . . . . 5.1 Illustration of neighbors navigation order used for LUT index computation . . . 5.2 Example of binary decision tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Example of binary decision diagram . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Masks used for crucial cliques detection. . . . . . . . . . . . . . . . . . . . . . . . 5.5 New masks for 2-crucial points detection . . . . . . . . . . . . . . . . . . . . . . . 5.6 New masks for 1-crucial points detection . . . . . . . . . . . . . . . . . . . . . . . 5.7 New masks for 0-crucial points detection . . . . . . . . . . . . . . . . . . . . . . . 5.8 Images used for the benchmark. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Computation speed results for ACK3 implemented with different levels of optimization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Computation speed results for PKD6 implemented with different levels of optimization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11 Speed gain obtained by different optimizations on ACK3. . . . . . . . . . . . . . 5.12 Speed gain obtained by different optimizations on PKD6. . . . . . . . . . . . . . 5.13 Thinning speed of fully optimized ACK3 and PKD6 on two data sets. . . . . . . 6.1 Masks used in Palágyi and Kuba 6-directional thinning . . . . . . . . . . . . . . 6.2 New end points masks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Masks used in our isthmus-based 6-directional thinning . . . . . . . . . . . . . . 6.4 Thinning speed of fully optimized ACK3, PKD6 and D6I1D on two data sets. . . xiii List of Figures xiv 6.5 Skeletons resulting of different algorithms for several shapes . . . . . . . . . . . . 6.6 Skeletons resulting of different algorithms for several shapes . . . . . . . . . . . . 6.7 Examples of homotopic kernels obtained using D6U. . . . . . . . . . . . . . . . . 6.8 Surface skeletons obtained using D6I2D. . . . . . . . . . . . . . . . . . . . . . . . 7.1 Example of undirected graph representation . . . . . . . . . . . . . . . . . . . . . 7.2 Illustrations of undirected graph properties . . . . . . . . . . . . . . . . . . . . . 7.3 Example of directed graph representation . . . . . . . . . . . . . . . . . . . . . . 7.4 Example of associated undirected graph . . . . . . . . . . . . . . . . . . . . . . . 7.5 Illustrations of directed graph properties . . . . . . . . . . . . . . . . . . . . . . . 7.6 Illustration of isomorphism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.7 Example of labeled graph representation . . . . . . . . . . . . . . . . . . . . . . . 7.8 Example of weighted graph representation . . . . . . . . . . . . . . . . . . . . . . 7.9 Areas in the shape where we have to find usable spatial positions to represent in the data tree. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.10 Classes of skeleton points in both asymmetric and symmetric skeletons . . . . . . 7.11 Example of connected components of classes of skeleton points . . . . . . . . . . 7.12 Example of data tree construction for incomplete model . . . . . . . . . . . . . . 7.13 Examples of spurious branches and ghost limbs . . . . . . . . . . . . . . . . . . . 7.14 Illustration of useless vertices problem . . . . . . . . . . . . . . . . . . . . . . . . 7.15 Illustration of splitted vertices problem . . . . . . . . . . . . . . . . . . . . . . . . 8.1 Examples of edit operations for labeled rooted trees. . . . . . . . . . . . . . . . . 8.2 Examples of edit operations for weighted trees. . . . . . . . . . . . . . . . . . . . 8.3 Example of mapping for labeled rooted trees . . . . . . . . . . . . . . . . . . . . 8.4 Example of invalidate mapping for weighted rooted trees . . . . . . . . . . . . . . 8.5 Example of alignment and alignment mapping . . . . . . . . . . . . . . . . . . . . 8.6 Example of constrained mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.7 Illustration of the center of three vertices . . . . . . . . . . . . . . . . . . . . . . 8.8 Example of invalidate constrained mapping . . . . . . . . . . . . . . . . . . . . . 8.9 Illustration of less constrained mapping . . . . . . . . . . . . . . . . . . . . . . . 8.10 Example of top down mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.11 Example of bottom-up mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.12 Hierarchy of mappings. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.13 Example of alignment for weighted trees . . . . . . . . . . . . . . . . . . . . . . . 9.1 Example of merging on directed graph. . . . . . . . . . . . . . . . . . . . . . . . . 9.2 Examples of merging on undirected graph. . . . . . . . . . . . . . . . . . . . . . . 9.3 Examples of homeomorphism. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.4 Example of merging kernel of an undirected graph. . . . . . . . . . . . . . . . . . 9.5 Example of merging kernel of a directed graph. . . . . . . . . . . . . . . . . . . . 9.6 Example of homeomorphic alignment. . . . . . . . . . . . . . . . . . . . . . . . . 9.7 Example of pruned tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.8 Illustration of possible rooting of unrooted tree . . . . . . . . . . . . . . . . . . . 9.9 Examples of cut operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.10 Example of bad homeomorphic alignment matching due to merging on model tree 10.1 Example of asymmetric homeomorphic relations on a set of homeomorphic trees. List of Tables 5.1 Look Up Tables for Crucial Point Detection . . . . . . . . . . . . . . . . . . . . . 5.2 Summarize of Optimizations Speed Gain (results expressed as "average [min, max]") 11.1 Time complexities of the different alignments. . . . . . . . . . . . . . . . . . . . . 11.2 Time complexities of the different alignments, taking into account the bounded maximal degree. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii List of Algorithms 1 Basic Thinning(X, W , D, k) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Parallel Thinning(X) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Directional Thinning(X) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Subfield Based Thinning(X) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Result Retrieval in Binary Decision Tree . . . . . . . . . . . . . . . . . . . . . . . 6 CK Thinning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 ACK Thinning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Optimized ACK Thinning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 . 1 :

 21 Figure 2.1: Examples of systems using markers. (a) exoskeleton of Gypsy TM mechanical system. (b) markers of Xsens TM inertial system. (c) markers of OptiTrack TM passive optical system.

Figure 2 .

 2 Figure 2.2 shows examples of the different model definitions.

Figure 2 . 2 :

 22 Figure 2.2: Examples of model definitions. (a) Kinematic model used by Menier et al. [111]; (b) Super-quadratic representation used by Sundaresan and Chellappa [166]; (c) Ellipsoid representation used by Horaud et al. [66] and (c') its associated implicit surface; (d) A surface scan of an actress and (d') corresponding mesh used by De Aguiar et al. [3]; (e) Textured mesh and (e') its underlying skeleton used by Ballan and Cortelazzo [11].

Figure 2 . 3 :

 23 Figure 2.3: (a) Illustration of silhouette cones, reprinted from Matusick et al. paper [109]. (b) Example of input camera image and associated silhouette. (c) and (d) Illustrations of surface and volumetric visual hulls reconstructions, respectively: object is represented in blue, its silhouettes in green, silhouette cones by dashed lines, and result by bold lines. (e) and (f) surface and volumetric visual hulls of a same object, respectively. Figures (b), (e) and (f) are reprinted from Menier thesis [110].

Figure 3 . 2 :

 32 Figure 3.2: Example of skeletonization process. Left: a visual hull. Middle: different steps of skeletonization process. Right: resulting skeleton (in red).

Figure 3 .

 3 3 for some examples of data tree extractions. More details about data tree extraction are given in Chapter 7.

Figure 3 . 3 :

 33 Figure 3.3: Examples of data tree extraction. Left: for a full human body subject. Right: for a hand subject. The black disk represents a point of the skeleton touching the border of the acquisition space.

Figure 4 .

 4 [START_REF] Agarwal | 3D human pose from silhouettes by relevance vector regression[END_REF] shows an example of skeleton in continuous framework.

Figure 4 . 1 :

 41 Figure 4.1: Example of skeleton in continuous framework. The shape (a rectangle) is represented in grey. The skeleton is represented by bold lines. Some maximal balls are represented by dashed lines, with their center represented by black points.

2 .

 2 a skeleton should be thin (should have lower dimension than the object), 3. a skeleton should be centered in the original object, 4. skeletonizing a skeleton should not change anything.

Figure 4 . 2 :

 42 Figure 4.2: Example of medial axis in discrete framework. The shape is represented in grey. The medial axis is represented by crosses. Some maximal balls are represented by dashed lines, with their centers represented by bold crosses.

Figure 4 . 4 :

 44 Figure 4.4: From the left to the right: 6-neighborhood, 18-neighborhood and 26-neighborhood of point p, represented by black disks and p.

Figure 4 . 5 :

 45 Figure 4.5: (a) A shape with a tunnel. (b) Deletion of the red point. (c) The shape has no tunnel.

  4.6,(a) and (c)).

Figure 4 . 6 :

 46 Figure 4.6: Examples of the different kinds of skeleton. Top, 2D skeletons: (a) ultimate skeleton, (b) curvilinear skeleton. Bottom, 3D skeletons: (c) ultimate skeleton, (d) curvilinear skeleton, (e) surface skeleton.

Fig. 4 . 6 ,

 46 Fig. 4.6,(b) and (d)).

Figure 4 . 7 :

 47 Figure 4.7: Examples of 1D isthmuses.

Algorithm 1 :

 1 Basic Thinning(X, W , D, k) Data: A k-connected shape X, a priority function D and a subset W ⊆ X Result: A skeleton of X while there exists a k-simple point in X\W do A = {y ∈ X\W |y is k-simple for X} B = {x ∈ A| for all y ∈ A, D(x) ≤ D(y)} Let z ∈ B X = X\{z} return X Algorithm 1 shows the basic thinning scheme based on a priority function. The set W , called

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Each number i represents a pixel which has to be out of the object to allow the deletion of p in the ith sub-iteration.

Figure 4 . 10 :

 410 Figure 4.10: Top, from the left to the right: Divisions of Z 3 into 2, 4 and 8 subfields. For an algorithm using k subfields, in a ith sub iteration, only points marked i will be considered. Bottom, from the left to the right: black points represent considered points in the first iteration of a k-subfield based algorithm, for k = 2, 4, 8.

Figure 5 .

 5 Figure 5.1 shows two examples of these orders.

Figure 5 . 1 :

 51 Figure 5.1: Left: an order optimized for reducing amount of cache miss. Right: an order useful for algorithms using several look-up tables for different neighborhoods.

Algorithm 5 :

 5 Result Retrieval in Binary Decision Tree Data: Binary Decision Tree T , input vector v = (b 0 , ..., b n ) Result: Value of function for v in input Node cur = root of T while cur is not a leaf do if value represented by cur equals 0 then cur = left child of cur else cur = right child of cur return value contained in cur Example Let consider the function f : {0, 1} 4 → {0, 1} such f (b 0 , b 1 , b 2 , b 3 ) = (b 0 AN D b 1 ) OR (b 2 AN D b 3 ). The results of f are represented by the binary decision tree in Figure 5.2. If we are searching for the value of f (0, 0, 1, 0), from the root, we go down twice in the left child (b 0 = 0 and b 1 = 0), then in the right child (b 2 = 1) and finally in the left child again (b 3 = 1): the result (0) is contained in the leaf where we are.

  b

Figure 5 . 2 :

 52 Figure 5.2: A binary decision tree containing the results of the function f .

Figure 5 .

 5 Figure 5.3 shows an example of transformation of a binary decision tree into a binary decision diagram.

Figure 5 . 3 :

 53 Figure 5.3: Transformation of a binary decision tree into a binary decision diagram.

Figure 5 . 4 :

 54 Figure 5.4: Masks used for crucial cliques detection.

Associated orientation Orientation 0 Orientation 5 Orientation 1 4 Figure 5 . 5 :

 51455 Figure 5.5: New masks for 2-crucial points detection, and all possible orientations of these masks.

10 Figure 5 . 6 :

 1056 Figure 5.6: New masks for 1-crucial points detection, and all possible orientations of these masks.

1 2 3 Orientation 4 5 6 7 Figure 5 . 7 :

 123456757 Figure 5.7: New masks for 0-crucial points detection, and all possible orientations of these masks.

Figure 5 . 8 :

 58 Figure 5.8: Images used for the benchmark. Above each image, the id number, the name, the size of the object and the size of the image.

Figure 5 . 9 :

 59 Figure 5.9: Computation speed results for ACK3 implemented with different levels of optimization, for each image of our test set.

Figure 5 . 10 :

 510 Figure 5.10: Computation speed results for PKD6 implemented with different levels of optimization, for each image of our test set.

Figures 5 .

 5 Figures 5.11 and 5.12 show these ratios for ACK3 and PKD6, respectively. Results are summarized in Table5.2.

Figure 5 . 11 :Figure 5 . 12 :

 511512 Figure 5.11: Speed gain obtained by different optimizations on ACK3.

  Now we have optimized the speed of thinning algorithms, we can try to use them in our context of motion capture. In order to test their speed in a real context, we propose to use two data sets: the dancer data set (height views) from 4Drepository from INRIA and a home made data set with four views recording gesture of a hand.

Figure 5 . 13 :

 513 Figure 5.13: Thinning speed of fully optimized ACK3 and PKD6 on two data sets.

6. 1 . 2

 12 Palágyi and Kuba 6-directional Thinning In 1998, Palágyi and Kuba [130] proposed a 6-directional thinning algorithm producing curvilinear skeletons, which we denote by PKD6. Each iteration of PKD6 consists of 6 sub-iterations in which points are removed from the shape if and only if the configuration of their 26-neighborhood matches at least one mask of the mask set M d for the given orientation d. The mask set M U used for the UP direction is presented in Figure 6.1. Mask sets M D ,M N ,M S ,M E ,M W for the other directions (DOWN, NORTH, SOUTH, EAST and WEST ) are obtained by appropriate rotations and reflections.

Figure 6 . 1 :

 61 Figure 6.1: Base masks M 1 -M 6 and their rotations around the vertical axis form the set of masks M U for direction UP. Points marked by a black disk have to belong to the object, points marked by a white disk have to be out of the object, and at least one point marked by a cross in masks M 1 and M 5 has to belong to the object.

Figure 6 . 2 :

 62 Figure 6.2: Masks M α -M γ and their rotations around the vertical axis are added to the set of masks M U for direction UP in order to remove ending points.

  5 and M β . Masks M ′ 2 ,M ′ 3 ,M ′ 4 ,M ′ 6 and M ′ 7 are respectively the same as M 2 ,M 3 ,M 4 ,M 6 and M γ . Using the sets M ′ d with d ∈ {U, D, N, S, E, W }, we can now define τ d , as proposed in Algorithm 12.

Figure 6 . 4 :

 64 Figure 6.4: Thinning speed of fully optimized ACK3, PKD6 and D6I1D on two data sets.

Figure 6 . 5 :Figure 6 . 6 :

 6566 Figure 6.5: Skeletons resulting of different algorithms for several shapes. From the left to the right: algorithms ACK3, D6I1D, and PKD6. From the top to the bottom: shapes 4, 18, 23, 20 and 13.

Figure 6 . 7 :

 67 Figure 6.7: Examples of homotopic kernels obtained using D6U.

Figure 6 . 8 :

 68 Figure 6.8: Surface skeletons obtained using D6I2D.

Figure 7 . 1 ,Figure 7 . 1 :

 7171 Figure 7.1: Example of undirected graph representation. The circles represent the vertices, the letter in each circle is just used in a purpose of easily identify the associated vertex, and the lines represent the edges. The represented graph is (V = {a, b, c, d, e}, E = {{a, c}, {a, d}, {c, d}, {c, e}}).

Figure 7 . 4 Figure 7 . 2 :

 7472 Figure 7.2 illustrates the above definitions.

Figure 7 Figure 7 . 3 :

 773 Figure 7.3: Example of directed graph representation. The circles represent the vertices, and the arrows represent the arcs. The represented graph is (V = {a, b, c, d, e}, E = {(a, c), (d, a), (d, c), (c, e), (e, c)}).

Figure 7 . 4 :

 74 Figure 7.4: From the left to the right, a directed graph and its associated undirected graph.

Figure 7 .

 7 Figure 7.5 illustrates the above definitions.

4 Figure 7 . 5 :

 475 Figure 7.5: From the left to the right, four directed graphs: G 1 , G 2 , G 3 and G 4 . G 1 have two roots(b and e), G 2 and G 3 have only one root (e) and G 4 have no root. G 2 , G 3 and G 4 are antisymmetric, contrary to G 1 . G 3 is a rooted tree, contrary to the three others. G 3 and G 4 are rooted forests, composed respectively by one and two rooted trees.

Figure 7 . 6 :

 76 Figure 7.6: From the top to the bottom: an example of isomorphism between undirected graphs (dotted lines represent the bijection), an example of isomorphism between directed graphs, and an example of non isomorphism between directed trees (bold vertex cannot find any correspondence respecting the definition).

Figure 7 . 7 :

 77 Figure 7.7: Example of labeled graph representation.In order to avoid confusion with vertex identifier, the label associated to each vertex is write between brackets. In the case where the identifier must be represented, we use the representation shown on the left graph: the identifier with the label between brackets in exponent. In other case, we just write the label between brackets, like in the right graph.

Figure 7 .

 7 Figure 7.8 for an example).

Figure 7 . 8 :

 78 Figure 7.8: Example of weighted graph representation. Each number represent the weight of the closest edge.

1 .Figure 7 . 9 :

 179 Figure 7.9: Areas in the shape where we have to find usable spatial positions to represent in the data tree.

Figure 7 .Figure 7 . 10 :

 7710 Figure 7.10 shows the distribution of the tree classes in both symmetric and asymmetric skeletons. We can observe that for the asymmetric skeleton, the branches are composed by 2-points and extremities are always E-points. Furthermore, I-points are present only at intersections of branches.

Algorithm 13 :Figure 7 . 11 :

 13711 Figure 7.11: Left: example of connected components, each rectangle representing a connected component (CC). In this example, there is five E-points CC, seven B-points CC and two I-points CC. Right: data tree constructed from the skeleton at left. Dot lines show the CC associate to each vertex.

Figure 7 . 12 :

 712 Figure 7.12: Example of data tree construction for incomplete model. Left: hand skeleton, with in orange, a border point. Right, the data tree constructed from the skeleton.

Figure 7 . 13 :

 713 Figure 7.13: Examples of spurious branches and ghost limbs(in red) for a full human body shape(left) and a hand shape(right).

Figure 7

 7

Figure 7 . 14 :

 714 Figure 7.14: (a): initial data tree (spurious branches in red). (b): data tree after spurious branches removal (in blue, resulting useless 2-degree vertices). (c): data trees we want to take into consideration during the matching.

Figure 7 .

 7 Figure 7.15 for an example.

Figure 7 . 15 :

 715 Figure 7.15: (a): data tree with splitted vertices (in red). (b): model tree to match. Dotted lines represent the expected matching.

Figure 8 . 1 :Figure 8 . 2 :

 8182 Figure 8.1: Examples of edit operations for labeled rooted trees.

Proposition 14 .

 14 γ(a + b, c) ≤ γ(a, 0) + γ(b, c) Proof. Due to triangle inequality, we have γ(a + b, c) ≤ γ(a + b, b) + γ(b, c). Due to translation invariance, we have γ(a, 0) = γ(a + b, b). Thus, γ(a + b, c) ≤ γ(a, 0) + γ(b, c).

Figure 8 . 4 :

 84 Figure 8.4: An invalidate mapping: the ancestor condition is not respected (due to bold pair:c is an ancestor of d, contrary to g which is not an ancestor of h).

  if there is a mapping line between them in the alignment. See right part of Figure 8.5 for an example. Jiang et al. provided an algorithm for the alignment between labeled ordered rooted trees with O(S 1 S 2 (D 1 + D 2 ) 2 ) time complexity. In the unordered case, the problem is NP-complete, excepted if D 1 and D 2 are bounded. In this case, the complexity is in O(S 1 S 2 ).

Figure 8 . 5 :

 85 Figure 8.5: Left: an alignment mapping. Right: the associated alignment.

Figure 8 . 6 :

 86 Figure 8.6: Left: a mapping which is constrained. Right: a mapping which is not constrained.

  Mapping for labeled unrooted trees In [189], Zhang et al. proposed a definition of constrained mapping for labeled unrooted trees, using the definition of the center of three vertices. Let u, v, w be three vertices in a labeled unrooted tree T , the center of the three vertices u, v, w, denoted by center(u, v, w), is the only vertex belonging to the three paths π(u, v), π(u, w) and π(v, w) (see Figure 8.7 for an example).

Figure 8 . 7 :

 87 Figure 8.7: Illustrations of the center, represented by the bold circle.

Figure 8 . 8 :

 88 Figure 8.8: An invalidate constrained mapping: the center relationship preservation condition is not respected.

5 Figure 8 . 10 :

 5810 Figure 8.10: At left: a mapping which is top down. At right: a mapping which is not top down.

5 Figure 8 . 11 :

 5811 Figure 8.11: At left: a mapping which is bottom up. At right: a mapping which is not bottom up.

Figure 8 . 12 :

 812 Figure 8.12: Hierarchy of mappings.

  horizontal-merge Merges a consecutive subsequence of siblings v 1 , ..., v s into a single vertex v. The children of v 1 , ..., v s become the children of v. horizontal-split The complement of horizontal-merge. Splits a vertex v into a consecutive sequence of siblings v 1 , ..., v s . The children of v become children of v 1 , ..., v s . vertical-merge Merge a sequence of vertices v 1 , ..., v s , where par(v i+1 ) = v i , 1 ≤ i < s, into a single vertex v. The children of v 1 , ..., v s not in the sequence become the children of v. vertical-split The complement of vertical-merge. Split a vertex v into a sequence of vertices v 1 , ..., v s , where par(v i+1 ) = v i , 1 ≤ i < s. The children of v become the children of the sequence v 1 , ..., v s .

Figure 8 . 13 :

 813 Figure 8.13: Example of alignment for weighted trees. Top left: initial trees G 1 and G 2 . Top right: trees after edge insertions G ′ 1 and G ′ 2 , dotted lines representing the alignment. Bottom: associated alignment mapping of G 1 and G 2 .

Figures 9 . 2 Figure 9 . 1 :Figure 9 . 2 :

 929192 Figures 9.1 and 9.2, respectively.

e mergings on h and k 1 1Figure 9 . 3 :

 193 Figure 9.3: Examples of homeomorphism.

Figure 9 . 4 :Figure 9 . 5 :

 9495 Figure 9.4: Example of merging kernel of an undirected graph.

Figure 9 . 6 :

 96 Figure 9.6: Example of homeomorphic alignment.

Figure 9 . 7 :

 97 Figure 9.7: From the left to the right: a tree T , a subtree of T rooted in b, a subtree of T rooted in a and pruned in e, and a forest of T with origin b.

Algorithm 14 :

 14 Homeomorphic Alignment Distance for Rooted Trees Data: pattern rooted tree P , data rooted tree D Result: η(P, D) = η(F(P, r P ), F(D, r D )); // Prop.18 begin foreach p ∈ V P , in suffix order do foreach A ⊆ F(P, p) do Compute η(A, ∅); // Prop.19 foreach p a ∈ anc(p) \ {p} do Compute η(P (p, p a ), ∅); // Prop.19 foreach d ∈ V D , in suffix order do foreach B ⊆ F(D, d) do Compute η(∅, B); // Prop.19 foreach d a ∈ anc(d) \ {d} do Compute η(∅, D(d, d a )); // Prop.19 foreach p ∈ V P , in suffix order do foreach d ∈ V D , in suffix order do foreach A ⊆ F(P, p) do foreach B ⊆ F(D, d) do Compute η(A, B); // Prop.21 foreach p a ∈ anc(p) \ {p} do foreach d a ∈ anc(d) \ {d} do Compute η(P (p, p a ), D(d, d a )); // Prop.20 end

Figure 9 . 8 :

 98 Figure 9.8: A tree G and the rooted trees G a and G e .

1 .

 1 w) at the end of L end Let P = (V P , E P , ω P ) and D = (V D , E D , ω D ) be two weighted trees. The initialization of L P and L D is in O(|V P | + |V D |) time complexity. Observe that, for p ∈ V P , p ′ ∈N (p) |anc(P, p, p ′ ) \ {p}| = |V P | -As a result, (p,p ′ )∈L P |anc(P, p, p ′ ) \ {p}| = (|V P | -1) × |V P | .

Algorithm 17 : 23 end 9

 17239 Homeomorphic Alignment Distance between Rooted and Unrooted Trees Data: pattern rooted tree P , data unrooted tree D Result: η(P, D) begin list of couples of vertices L D ← computeOrder(D) foreach p ∈ V P , in suffix order do foreach A ⊆ F(P, p) do Compute η(A, ∅); // Prop.19 foreach p a ∈ anc(p) \ {p} do Compute η(P (p, p a ), ∅); // Prop.19 foreach (d, d ′ ) ∈ L D do foreach B ⊆ F(D d ′ , d) do Compute η(∅, B); // Prop.19 foreach d a ∈ anc(D, d, d ′ ) \ {d} do Compute η(∅, D da (d, d a )); // Prop.19 foreach p ∈ V P , in suffix order do foreach (d, d ′ ) ∈ L D do foreach A ⊆ F(P, p) do foreach B ⊆ F(D d ′ , d) do Compute η(A, B); // Prop.21 foreach p a ∈ anc(p) \ {p} do foreach d a ∈ anc(D, d, d ′ ) \ {d} do Compute η(P (p, p a ), D da (d, d a )); // Prop.20 Compute η(P, D); // Prop.18 and Prop.

Figure 9 . 9 :

 99 Figure 9.9: Examples of cut operations. Top row: on an unrooted tree. Bottom row: on a rooted forest.

Figure 9 .

 9 Figure 9.10: Example of bad homeomorphic alignment matching due to merging on model tree: the "thumb" edge is removed for a cost of 6, as the merging of "arm" edge and "palm" edge in both trees reduce the matching cost by 8.

10. 1 ≥ G ′ . 119 10. 1 . 2 Properties

 111912 Asymmetric Homeomorphism 10.1.1 DefinitionLet G = (V, E, ω) and G ′ = (V ′ , E ′ , ω ′ ) be two weighted graphs, G is homeomorphically greater than or equal to G ′ if and only if there exists an isomorphism between G ′ and a graph obtained by mergings on G. We denote this relation by G H In this section we denote by H = the homeomorphism relation, and by I = the isomorphism relation. 10.1.2.1 Usual definitions on binary relations A binary relation R on a set S is a collection of ordered pairs of elements of S. Let a and b two elements of S, the statement (a, b) ∈ R is denoted by aRb. The binary relation R is an equivalence relation if and only if for all a, b and c in S:

  R b and b R a then a = b (antisymmetry) • if a R b and b R c then a R b (transitivity) An order relation R is total if for all a and b in S, a R b or b R a. In the other case, R is a partial order relation. 10.1.2.2 Properties of asymmetric homeomorphism Contrary to isomorphism and homeomorphism which are equivalence relations on the set of graphs G, the asymmetric homeomorphism is an partial order relation on G. More precisely, as the asymmetric homeomorphism is a sub case of homeomorphism, there is no asymmetric homeomorphism relations between elements of different equivalence class by homeomorphism on G. If we consider a graph G, the equivalence class H = (G) is partially ordered by H ≥ (see Figure 10.1 for an example) and have a least element equal to MK(G): for all G ′ ∈ H = (G), G ′ H ≥ MK(G).

Figure 10 . 1 :

 101 Figure 10.1: Example of asymmetric homeomorphic relations on a set of homeomorphic trees.

  structural noise: we randomly add new vertices by three ways, corresponding to the different types of noise: splitting an existing vertex, and linking the two parts by a 0-weighted edge (splitted vertices), adding a new 2-degree vertex by the split of an edge (useless 2-degree vertices), and adding a new 1-degree vertex, linked by an randomly weighted edge (spurious branches and ghost limbs).

Figure 11

 11

Figure 11 . 1 :

 111 Figure 11.1: Illustration of structural noise generation.

Figure 11 . 2 :

 112 Figure 11.2: Illustration of precision measurement. Left: 0% of good matchings. Right: 100% of good matchings.

Figure 11 . 3 :

 113 Figure 11.3: Precision without weight variation. Top: pattern tree with 10 vertices. Bottom: pattern tree with 20 vertices. Left: without cut operation. Right: with cut operation.

Figure 11 . 4 :

 114 Figure 11.4: Precision with 10 percent of weight variation. Top: pattern tree with 10 vertices. Bottom: pattern tree with 20 vertices. Left: without cut operation. Right: with cut operation.

Figure 11 . 5 :

 115 Figure 11.5: Precision with 50 percent of weight variation. Top: pattern tree with 10 vertices. Bottom: pattern tree with 20 vertices. Left: without cut operation. Right: with cut operation.

Figure 11 . 6 :

 116 Figure 11.6: Example where alignment provide better result than homeomorphic alignment.

Figure 11 . 7 :

 117 Figure 11.7: Computation speed of different alignments for different sizes of trees. Top left: between two rooted trees. Top right: between a rooted tree and an unrooted tree. Bottom left:between two unrooted trees.

  P S D C) O(S P S D (C + H P H D )) O(S P S D (C + H D )) unrooted trees O(S P S D C) O(S P S D (C + S P S D )) O(S P S D (C + S D )) rooted tree with unrooted tree O(S P S D C) O(S P S D (C + H P S D )) O(S P S D (C + S D ))

. 2 summarizes

 2 the time complexities of the different alignments, taking into account the bounded maximal degree. If C >> H P + H D , the three alignments between rooted trees will have closed speed, as shown in Figure 11.7.

  .

Figure 11 . 8 :

 118 Figure 11.8: Distributions for a grid sided by 40. Left: dancer data set. Right: hand data set.

Figure 11 . 9 :

 119 Figure 11.9: Distributions for a grid sided by 80. Left: dancer data set. Right: hand data set.

Figure 11 . 10 : 6 Figure 11 . 11 :

 111061111 Figure 11.10: Distributions for a grid sided by 120. Left: dancer data set. Right: hand data set.

Figure 11 . 6 Figure 11 . 12 :

 1161112 Figure 11.13 for full human body application.

7 d=3 25 s=1 6 Figure 11 . 13 :Figure 11 . 14 :

 725611131114 Figure 11.13: Frequencies of asymmetric homeomorphic alignment of full human body rooted model tree with different unrooted data trees.

Figure 12 . 3 :

 123 Figure 12.3: A visual hull of hand where the skeleton finger branches are not proportional to the real lengths of the fingers.

  the 3d point associated to this vertex.The coordinate constraint on a m can be expressed as follows: let a d = (v d , w d ) be an arc of the data tree (possibly obtained by some merging), a d can match with a m if and only if P (v d ) and P (w d ) have relative positions satisfying the constraint. For example, for the head-torso constraint, a d can match with a m if and only if P y (w d ) > P y (v d ).We propose to integrate coordinate constraints in the alignment computation, by modifying the resizement cost formulation for constrained edges. We consider the case of asymmetric homeomorphic alignment for rooted trees (AHA between unrooted trees just consists in searching the best AHA for all rooted versions of the trees, please refer to Chapter 10 for more details). The cost of the resizement of a model arc a m = (v m , w m ), weighted by W m , to the weight W d of a data arc a d = (v d , w d ) is equal to +∞ if there is a coordinate constraint C associated to a m and the relative positions of P (v d ) and P (w d ) are not compatible with C. Otherwise, the cost is still γ(W m , W d ).

Figure 12 . 4 :

 124 Figure 12.4: All points in the grey area are between the two points represented by black crosses: at left, in the sense of R c , at right, in the sense of R p .

Figure 12 . 5 :

 125 Figure 12.5: All points in the grey area are between the two points represented by black crosses, in the sense of R b .

  Postulate B': a = b and a = c and b = c and a, b, c not equidistant ⇒ bac or cab or abc or cba or acb or bca. Postulate D': abc ⇒ a = b and a = c and b = c and a, b, c not equidistant.

  By definition, in any non-equilateral triangle, there is at least one edge with a length different from those of the other edges. We denote by d(a, b) the Euclidian distance between two points a and b. We denote by δ : R 3 × R 3 × R 3 → R the function: δ(a, b, c) = |d(a, b)d(b, c)| An interesting property is that if a,b and c are the three distinct vertices of a non-equilateral triangle, one and only one of the three values δ(a, b, c), δ(a, c, b) and δ(c, a, b) is strictly less than the two others. This can easily be proved by showing that the only case where two of the three values are equal is for isosceles triangle, and in that case, these values are strictly greater than the last one, which is equal to zero. Using this property, we can define a betweenness relation R u satisfying the four postulates. Let a, b and c be three distinct and non-equidistant points, b is between a and c (R u (a, b, c) is true) if and only if δ(a, b, c) ≤ δ(b, c, a) and δ(a, b, c) ≤ δ(b, a, c).

Figure 12 . 6 :

 126 Figure 12.6: All points in the grey area are between the two points represented by black crosses(in the sense of R u ). We can notice that the green cross is between the black ones, in the common conception of betweenness, but is not in the sense of R u . In the same way, the red cross is not between the black ones, in the common conception of betweenness, but is so in the sense of R u .

  First we have to introduce a new postulate providing concept of intuition for the common conception of betweenness. Definition 31. Let R a ternary relation on S. The couple (R 3 , R) is an intuitive betweenness system if: • it satisfies the postulates A,B',C and D'. • Postulate I: abc ⇒ orthogonal projection a ′ of a on the line (bc) is between b and c (in the sense, d(a ′ , b) + d(a ′ , c) = d(b, c)).

  Using these properties, we can define a betweenness relation R i satisfying the four postulates and which is intuitive. Let a, b and c be three distinct and non-equidistant points, b is between a and c (R i (a, b, c) is true) if and only if one of the following conditions holds: 1. d(a, c) > d(b, c) and d(a, c) > d(b, a).

2 .

 2 d(b, c) = d(b, a) and d(a, c) = d(b, c). This relation is illustrated in Figure 12.7.

Figure 12 . 7 :

 127 Figure 12.7: All points in the grey area are between the two points represented by black crosses(in the sense of R i ).

Figure 13 . 1 :

 131 Figure 13.1: Samples of the different used data sets: (a) dancer, (b) children, (c) dog and (d) hand.

Figure 13 .Figure 13 . 2 :

 13132 Figure 13.2 shows computation time and frame rate for different data sets and different voxel grid resolutions. It can be observed that for each data set, real-time (40 milliseconds or 25 iterations per second) is reached for voxel grid resolutions lower than 100 × 100 × 100. It is

Classes 3

 3 and 6 can be easily detected and rejected by our method, due to the lack of part positions. Classes 1 and 2 (good matchings and missing limb matchings) are considered as usable for initialization (we call them True Positive, shorted in TP), and classes 4 and 5 are considered as unusable for initialization (we call them False Positive, shorted in FP).

Figure 13 . 3 :

 133 Figure 13.3: Examples for the different classes of matchings: (a) good matching, (b) missing limb matching, (c) partial matching, (d) misplaced limb matching, (e) bad matching.

  Distance Historic for Dancer set (res. 50x50x50) Distance Historic for Dancer set (res. 150x150x150)

Figure 13 . 4 :

 134 Figure 13.4: Study of correlations between AHA distance and matching quality for dancer data set with different grid resolutions. Left: historic of AHA distances and matching qualities all along the sequence. Right: different ratios in regard of given AHA distance threshold. From top to bottom: for resolutions 50 3 , 100 3 and 150 3 .

Figure 13 . 5 :

 135 Figure 13.5: Study of correlations between AHA distance and matching quality for children data set with different grid resolutions. Left: historic of AHA distances and matching qualities all along the sequence. Right: different ratios in regard of given AHA distance threshold. From top to bottom: for resolutions 50 3 , 100 3 and 150 3 .

  betweenness relation to use for the betweenness constraints solving. For this purpose, we propose to compare the results for the hand data set, using on one hand the "ball-based" intuitive betweenness relation R b defined in Section 12.2.2.1, and on the other hand the "mathematically efficient" intuitive betweenness relation R i defined in Section 12.2.2.2. The measures are shown in Figure13.6 and Figure13.7, respectively. Distance Historic for Hand set using Rb(res. 50x50x50) Distance Historic for Hand set using Rb(res. 100x100x100) Distance Historic for Hand set using Rb(res. 150x150x150)

Figure 13 . 6 :

 136 Figure 13.6: Study of correlations between AHA distance and matching quality for hand data set with different grid resolutions, using betweenness relation R b . Left: historic of AHA distances and matching qualities all along the sequence. Right: different ratios in regard of given AHA distance threshold. From top to bottom: for resolutions 50 3 , 100 3 and 150 3 .

Figure 13 . 7 :

 137 Figure 13.7: Study of correlations between AHA distance and matching quality for hand data set with different grid resolutions, using betweenness relation R i . Left: historic of AHA distances and matching qualities all along the sequence. Right: different ratios in regard of given AHA distance threshold. From top to bottom: for resolutions 50 3 , 100 3 and 150 3 .

Figure 13 . 8 :

 138 Figure 13.8: Samples of output poses for different models. Model parts positions are represented in blue, with the associated label in red. Approximations of skeleton branches between the parts are drawn green.

Figure 13 . 9 :

 139 Figure 13.9: Sample of accuracy difference of pose estimation for different grid resolutions: (a) 50 × 50 × 50 resolution, (b) 150 × 150 × 150 resolution.

Figure 13 . 10 :

 1310 Figure 13.10: Examples of visual hull segmentations based on our method.

  1p 1 |, |q 2p 2 |) ≤ 1}.Figure 4.3: Left: 4-neighborhood of p represented by boxes labeled by 4 or p. Right: 8neighborhood of p represented by boxes labeled by 8 or p.

		4		8	8	8
	4	p	4	8	p	8
		4		8 8 8

Table 5 .

 5 1.

Table 5 . 1 :

 51 Look Up Tables for Crucial Point Detection

		related mask number of values	value length	Total size
				in bits in bytes
	LU T 2X	M 2X	2 26	6	1	64 Mb
	LU T 2S	M 2S	2 6	6	1	64 b
	LU T 1X	M 1X	2 26	12	2	128 Mb
	LU T 1S	M 1S	2 18	12	2	512 Kb
	LU T 0X	M 0X	2 26	8	1	64 Mb
	LU T 0S	M 0S	2 26	8	1	64 Mb

Notice that even if the masks proposed in the Figures 5.5, 5.6 and 5.7 do not require the consideration of more than 17 neighbors, the number of values of the LUT corresponds to full 26-neighborhood. It is due to the fact that LUT stores the results for each orientation of the masks, thus we have to consider not only the points considered for one orientation, but for all orientations at the same time.

Table 5 .

 5 2.

	speed up	0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 23 22 21 20	CCSU Gain LUT Gain Full Gain
		1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
		shape

Table 5 .

 5

			2: Summarize of Optimizations Speed Gain
		(results expressed as "average [min, max]")
		CCSU Speed Gain LUT Speed Gain CCSU+LUT Speed Gain
	ACK3	1.46 [1.04, 2.06]	4.62 [0.69, 12.25]	7.44 [0.73, 22.56]
	PKD6	1.71 [0.84, 2.93]	1.20 [0.68, 1.53]	2.04 [1.01, 3.86]

  The strategy used by Bertrand and Couprie consists in dynamically detecting the considered isthmuses (1D or 2D) and accumulating them in a constraint set. A point can be deleted if it is simple and not in the constraint set. By this way, a point detected as an isthmus in a given iteration cannot be removed later. Depending on the considered isthmuses, the method provides curvilinear skeletons (in case of 1D isthmuses) or surface skeletons (in case of 2D isthmuses). In both cases, this strategy provides skeletons with very few spurious elements.

proposed a sequential thinning based on isthmuses (see Section 4.2.2.2). Contrary to extremity points, isthmuses are not simple points. Due to this fact, isthmuses cannot be used in classical strategy in which a point is removable if and only if it is simple and does not match certain conditions. 57

  3 * N +H D * N )) time complexity (O(S 2 H D ) if the maximal degree is bounded).

Table 11

 11

	.1 summarizes the complexities of classic Jiang et al. alignment(α), homeomorphic
	alignment(η), and asymmetric homeomorphic alignment( η) for trees with bounded degree.

Table 11 . 1 :

 111 Time complexities of the different alignments.

		α	η	η
	rooted trees	O(S P S D ) O(S P S D H P H D ) O(S P S D H D )
	unrooted trees rooted tree with unrooted tree O(S P S D ) O(S P S D )	O(S 2 P S 2 D ) O(S P H P S 2 D )	O(S P S 2 D ) O(S P S 2 D )
		127

Table 11 . 2 :

 112 Time complexities of the different alignments, taking into account the bounded maximal degree.

  for an example) if all parts are positioned at appropriate locations; 2. an output pose is a missing limb matching (see Figure 13.3b for an example) if one part has no position, other parts being positioned at appropriate locations; 3. an output pose is a partial matching (see Figure 13.3c for an example) if several parts have no position, other parts being positioned at appropriate locations; 4. an output pose is a misplaced limb matching (see Figure 13.3d for an example) if at least one part is positioned at inappropriate location; 5. an output pose is a bad matching (see Figure 13.3e for an example) if almost all parts are positioned at inappropriate locations;6. an output pose is not a matching if at most one part is matched.

Acknowledgements

Final algorithm

We can compute Homeomorphic Alignment for unrooted trees (see Algorithm 16) with improved complexity, using this order of navigation. end Complexity Let T = (V, E, ω) be an unrooted tree. Each vertex of degree 1 (one) will be put in the queue, then, for each edge, each of its vertices will be put twice in the queue. As

3. ∃f ∈ A R , f being obtained by merging (j a , j) with other arcs. In this case, there is one and only one child j c of j, such (j, j c ) is merged with (j a , j) in f , and then all D(j ′ c , j), j ′ c ∈ C(j) \ {j c } are deleted.

4. (i a , i) is deleted, then only one P (i c , i), i c ∈ C(i) is aligned with D(j, j a ), the other being removed.

Cases 1,2,3,4 justify, respectively, the lines 1,2,3,4 of the expression of η(P (i, i a ), D(j, j a )) in the proposition.

The two last cases cannot lead to a better homeomorphic alignment: 5. The deletion of (i a , i) and (j a , j) cannot be preferred to the resizement (possible case 2), because γ(ω P (i a , i), 0) + γ(0, ω D (j a , j)) ≥ γ(ω P (i a , i), ω D (j a , j)).

6. If (j a , j) was deleted, and not (i a , i), then only one D(j c , j), j c ∈ C(j) is aligned with P (i, i a ), the other being removed. It is less expensive to merge (j a , j) with (j, j c ) (possible case 3), because γ(0, ω D (j a , j)) + γ(ω P (i a , i), ω D (j, j c )) ≥ γ(ω P (i a , i), ω D (j a , j c )).

Algorithm

Using the reformulations defined above, we can now design our algorithm following the bottomup approach (see Algorithm 18).

Complexity

We denote by N the maximum degree of a vertex in the pattern tree or in the data tree, H P and H D the maximal height of both trees and S the maximal size (number of vertices) of both trees.

We recall that a set of size n has 2 n subsets. The computation is in O(S 2 * (N * 2

Optional Constraints

The method proposed in the previous section already provides some good results. However there are still some special cases which can not be solved by our method, even with sufficient topological information. [67] A system is a near-betweenness system if it possess all but one of the five properties of a betweenness system.

In our case, S = R The new constraints can be integrated in our a priori models in order to solve the problematic cases, as shown in Figure 12.9.

Vertices: For the hand model, a betweenness constraint is set on the five fingers, in order to allow the pose estimation for a natural initial pose of the hand, for example when the hand is wide open.

For the full human body model, coordinate constraints are set on couples head-torso, and torsocrotch, allowing a high variety of initial poses, the only requirement being that the subject stands up and that the different parts of the body are not too close to each other.

Chapter 13

Results and Discussion

In this chapter, we present the results obtained with our pose initialization method, by considering the speed of our method in different situations (voxel grid size, subject), and by showing some samples of output poses.

Then, using empirical data, we discuss the choice of the good threshold to use for the acyclic homeomorphic alignment distance, and of the most adapted betweenness relation to use for betweenness constraints.

Finally, we propose some ways to use our method, in combination with other works found in the literature.

13.1 Implementation and Data Sets

Implementation

Our method has been tested on a computer with a processor Intel(R) Core(TM) 2 Quad Q8200 at 2.33 GHz, a GPU Nvidia(R) Geforce(TM) 9800 GT and 3 Go of RAM. Our program is implemented in C++, and the visual hull is computed on GPU, using GLSL.

The program has two modes:

• In online mode, input images are acquired from calibrated cameras (calibration can be done in a prior step, or can be loaded using previously saved data), and silhouettes are extracted using background extraction.

• In offline mode, input images are loaded from a data set directory, containing input images, silhouette images (or background images used to extract them), and all required parameter files (e.g. for camera calibration).

153

We can observe that the rate of true positives for a given AHA distance threshold is nearly invariant in regard of the chosen betweenness relation. However, the rate of false positives is sensibly higher when the R i betweenness is used. A quick observation of the measures allows us to see that the use of R b provides a larger amount of partial matchings than the use of R i , and a fewer amount of bad matchings, especially in case of low grid resolutions.

This fact can be easily explained by analyzing the differences between these two relations: for three distinct and non-equidistant points, there is always one point which is between the two others, in the sense of R i , even if the points are close to be equidistant. It is not the case for R b : when three points are close to be equidistant, none is between the others.

For low grid resolutions, positions of fingers are more approximative, leading to ambiguous cases (where three fingers are equidistant). Using R i , the method will always provide an order, even for these ambiguous cases, possibly leading to a bad matching. Using R b the method will provide an order only for non-ambiguous cases, the other leading to partial matchings, which can be detected and rejected.

The most adapted betweenness relation for our purpose is thus the "ball-based" betweenness relation R b , which is also the fastest to compute.

Concerning the AHA distance threshold, a good value seems to be 3.5, providing a good matching for between 10 and 20 percent of the frames, with a fail rate lower than five percent. Please remember that this data set is of very bad quality, with only four non-optimally placed cameras, and basic silhouette extraction. Considering this fact, the results are rather good.

Output Pose Samples

In order to show the quality of the output poses, the large amount of poses supported by our method and its adaptability, we recorded several outputs of our application for the different data sets, and present them in Figure 13.8.

Discussion

The choice of the grid resolution will depend on the application: a low grid resolution like 50 3 is suited for a fast computation, but less accurate pose estimations (see Figure 13.9) and higher rate of false positives. On the other hand, high resolution like 150 3 provides accurate good matchings and few false positives, but does not always allow a real-time use.

In our opinion, a 100 3 grid resolution is a good compromise (at least for the studied data sets), the accuracy being sufficient and the matching rates being close to those obtained with higher resolutions. Furthermore, this resolution allows a real-time use, even in the worst observed case.

Part V Conclusion

"The hardest thing is to go to sleep at night, when there are so many urgent things needing to be done. A huge gap exists between what we know is possible with today's machines and what we have so far been able to finish." Donald Knuth.