Florent Grélard
email: florent.grelard@labri.fr

Fabien Baldacci

Anne Vialard

Jean-Philippe Domenger

New methods for the geometrical analysis of tubular organs

Keywords: Tubular organ analysis, Orthogonal planes, Centerlines, Voronoï covariance measure

This paper introduces new methods to study the shape of tubular organs. Determining precise cross-sections is of major importance to perform geometrical measurements, such as diameter, wall-thickness estimation or area measurement. Our first contribution is a robust method to estimate orthogonal planes based on the Voronoi Covariance Measure. Our method is not relying on a curve-skeleton computation beforehand. This means our orthogonal plane estimator can be used either on the skeleton or on the volume. Another important step towards tubular organ characterization is achieved through curve-skeletonization, as skeletons allow to compare two tubular organs, and to perform virtual endoscopy. Our second contribution is dedicated to correcting common defects of the skeleton by new pruning and recentering methods. Finally, we propose a new method for curve-skeleton extraction.

Introduction

Tubular organ analysis is essential for the physiological understanding of this kind of organs and for the characterization of related diseases. Many major diseases involve tubular organs. For instance, in the case of airway-trees, chronic obstructive pulmonary disease (COPD) has been reported to be one of the major causes of death [START_REF] Murray | Evidence-Based Health Policy-Lessons from the Global Burden of Disease Study[END_REF]. Another example is coronary heart disease, which is linked to obstructed vessels, and is the first cause of death in the US [START_REF] Members | Heart disease and stroke statistics : 2008 update: A report from the american heart association statistics committee and stroke statistics subcommittee[END_REF]. On another note, images of neurons acquired from laser-scanning microscopy help building the relationship between their morphology and their function [START_REF] Dima | Computer aided image segmentation and graph construction of nerve cells from 3d confocal microscopy scans[END_REF]. All the above-mentioned organs can be considered tubular, that is to say they are 3D objects with a circular or elliptical crosssection and elongation in one direction.

A lot of work and recent improvements have been made to segmentation techniques regarding airwaytrees [START_REF] Lo | Vessel-guided airway tree segmentation: A voxel classification approach[END_REF] and to neuron reconstruction [START_REF] Janoos | Robust 3d reconstruction and identification of dendritic spines from optical microscopy imaging[END_REF]. However, various charac-terization approaches still have defects with respect to tubular organs, as explained in Section 2.

The study of tubular organs relies on measurements such as lumen area, wall thickness, and diameter estimation, and has applications such as virtual endoscopy and shape matching. Measurements on tubular organs are precise only if they are taken from a cross-section orthogonal to the tube. The usual approach to obtain orthogonal planes consists in computing the curve-skeleton (or centerline). The cross-sectional plane is defined by the tangent at a given skeletal point on the skeleton. Precise measurements can then be performed in the computed cross-section. Thus, orthogonal plane estimation is central but is usually heavily reliant on the properties of the curve-skeleton.

Curve-skeletons, or centerlines, are simplifications of shapes consisting in a set of curves. They allow to make geometrical measurements such as tangent estimation [START_REF] Postolski | Tangent estimation along 3d digital curves[END_REF] or tortuosity [START_REF] Lang | Three-dimensional quantification of capillary networks in healthy and cancerous tissues of two mice[END_REF], and to perform virtual endoscopy [START_REF] Bauer | Extracting curve skeletons from gray value images for virtual endoscopy[END_REF]. However, orthogonal plane estimation and skeleton computation are impeded due to the organ's variability in diameter, to the presence of junctions (one tube divides in two or more) and to segmentation issues. Curve-skeletons can contain irregularities [START_REF] Kasthuri | Saturated reconstruction of a volume of neocortex[END_REF] (in gray) and skeleton computed by the method in Lee et al. (1994) (in red). Orthogonal planes are estimated using the λ-MST method of [START_REF] Postolski | Tangent estimation along 3d digital curves[END_REF] (in blue). The skeleton contains irregularities and spurious branches (see green points in the close-ups) which results in a deviated estimated tangent. As a result, orthogonal plane estimation is impacted.

(i.e. they are not smooth) which impact orthogonal plane computation (see Fig. 1). Moreover, spurious branches appear as a result of small irregularities on the object's surface which are kept through the skeletonization process. The orthogonal planes computed on faulty branches are not orthogonal to the tube. Thus, filtering the skeleton is often a necessary post-processing step. Pruning consists in the removal of irrelevant parts of the skeleton, yielding a skeleton representing the original shape more aptly. Once all the extra branches have been removed, there is a one-to-one correspondence between the parts of the skeleton and those of the volume, which allows to analyze the organ further. One of the main challenges of pruning approaches is having a complete and sensitive method, meaning a method which removes all spurious branches while keeping meaningful ones. This problem arises specifically for varying-diameter tubes, where the geometry of spurious branches may vary depending on the tube radius.

Once the spurious branches have been removed, it is possible to analyze the organ more easily. Shape decomposition is part of this process. It consists in breaking a complex shape into smaller and meaningful parts. In the case of the airway-tree, each bronchus is separated from another. Nevertheless, obtaining clear delineations in lower diameter tubes is not easily achievable.

Recentering is also key to filtering the skeleton. Indeed, not all skeletons are centered or smooth inside the shape, negatively impacting geometrical analyses and virtual endoscopy. Thus, for instance, a tangent vector computed at a non-centered skeleton point is potentially not aligned with the tube's axis.

The contribution of this paper is manifold. First, we introduce a robust orthogonal plane estimator. Our method is based on the Voronoi Covariance Measure (VCM), [START_REF] Mérigot | Voronoibased curvature and feature estimation from point clouds[END_REF]. This measure describes the shape of the Voronoi cells generated by the skeleton points or directly by the volume points. In the following section, we show how Voronoi cells have a simple and intuitive connection to orthogonal plane normals. The main advantage of our approach is its robustness to small irregularities in the curve. Moreover, unlike most orthogonal plane estimators, it does not necessarily rely on a curve, but it can be computed directly from a tubular volume.

Second, we propose original algorithms to improve existing curve-skeletons. We use the orthogonal planes to filter existing skeletons, which involves pruning and recentering steps.

Finally, we propose our own curveskeletonization approach. Our algorithm consists in a tracking procedure, which direction is given by the orthogonal planes. It consists in two steps: the first one is dedicated to tubular parts, the second one to junctions. We show our skeleton is complete and possesses few to no faulty branches.

The paper is organized as follows. In Section 2, we describe the state-of-the-art methods for curveskeleton extraction, and orthogonal plane estimation. Section 3 presents our method for robust orthogonal plane estimation. In Sections 4 and 5, existing skeletons are filtered and improved using orthogonal planes. In Section 6, we introduce our own centerline tracking algorithm.

Related work

The most common workflow to obtain robust measurements on tubular organs consists in: segmenting the organ; extracting the curve-skeleton from the volume; computing cross-sections from the curve-skeleton; and making measurements from the cross-sections. A workflow for airway-tree analysis following these four steps is described in [START_REF] Tschirren | Intrathoracic airway trees: segmentation and airway morphology analysis from low-dose ct scans[END_REF]. This section presents previous works aiming at extracting and filtering curve-skeletons and at estimating orthogonal planes.

There are various existing algorithms for curveskeleton extraction. [START_REF] Tagliasacchi | 3D Skeletons: A State-of-the-Art Report[END_REF] described a list of desirable properties for a curveskeleton to possess. In the context of tubular organs, not all properties are necessary. The curveskeleton must be complete, centered, thin and not impacted by noise.The curve-skeleton can be extracted directly from a gray-level image, but an implicit segmentation step is usually involved to capture the geometry of the organ [START_REF] Flasque | Acquisition, segmentation and tracking of the cerebral vascular tree on 3d magnetic resonance angiography images[END_REF].

Thus, in the context of varying-diameter tubes with junctions, we focus on curve-skeleton methods based on segmented sets of voxels. According to the survey by [START_REF] Cornea | Curve-skeleton properties, applications, and algorithms[END_REF] we can classify the curve-skeleton algorithms for discrete objects in the following classes: (a) thinning algorithms and (b) distance field or general field based algorithms.

Thinning methods consist in the removal of voxels from the boundary until a thin version of the shape is obtained. They are comprised, among many others, of the work of [START_REF] Lee | Building skeleton models via 3-d medial surface/axis thinning algorithms[END_REF], [START_REF] Palágyi | Quantitative analysis of pulmonary airway tree structures[END_REF] and [START_REF] Couprie | Discrete bisector function and Euclidean skeleton in 2D and 3D[END_REF] (see Fig. 2). The first two compute the skeleton on any given 3D object whereas the latter is applied to tubular shapes. The authors in [START_REF] Couprie | Discrete bisector function and Euclidean skeleton in 2D and 3D[END_REF] propose to compute the skeleton by filtering the medial axis (the centers of all the local maximal balls that can fit inside a shape) with an angle criterion. The resulting filtered medial axis is used as a constraint set for the computation of the skeleton by homotopic thinning. More importantly, and this is true for all thinning methods, it can create small faulty branches (see close-up in Fig. 2a), often removed automatically by a pruning step, which are, in our case, difficult to distinguish from small branches that are representative of the volume. The authors define two parameters in order to get rid of spurious branches. However, in practice, some small spurious branches remain, regardless of the parameter values.

Potential field methods generally produce very smooth skeletons. For instance, we can cite [START_REF] Cornea | Computing hierarchical curve-skeletons of 3D objects[END_REF], [START_REF] Hassouna | Variational curve skeletons using gradient vector flow[END_REF] where the skeleton is extracted from a potential field computed on the volume. In [START_REF] Cornea | Computing hierarchical curve-skeletons of 3D objects[END_REF], the authors compute a hierarchical skeleton based on a Newtonian potential field. Each boundary point is considered as an electric charge, repelling interior points of the volume. From the potential vector field, critical points (points where the magnitude of the force vector vanishes) are extracted and connected thanks to paths, yielding a first hierarchy in the skeleton. Then, two levels of hierarchy are added through high divergence value and highcurvature points. The main drawback of this kind of algorithm is that the resulting skeleton is not necessarily connected, can be incomplete (see Fig. 2b) and can contain spurious branches (see Fig. 2c).

In summary, the various classes of curve-skeleton algorithms have defects, including introducing irregularities, being incomplete, or having faulty branches. Pruning consists in removing faulty branches while keeping meaningful ones with respect to the object. The general idea is to define a significance measure for each point: this value corresponds to the relevance of a point in the skeleton, and is preferably low for spurious branches. Various significance measures have been defined and reviewed in [START_REF] Shaked | Pruning medial axes[END_REF]. These measures aim at determining if a branch stems from a small irregularity on the surface. For instance, maximal thickness of implied erosion [START_REF] Arcelli | Euclidean skeleton via centre-of-maximal-disc extraction[END_REF] estimates the part of a volume which is encoded by the branch. Another measure is dedicated to the ratio between the length of the branch and the number of closest boundary points to a branch point. A recent work intro- duced in Serino and Sanniti di Baja (2014) involves pruning 3D-curve skeletons. It considers the curveskeleton as a graph, where each branch is an edge. Pruning is done in a hierarchical manner. As a first step, all peripheral edges are examined. Then, a neighboring internal edge can be pruned only if all of its neighboring peripheral edges have been deleted. The authors use the combination of four significance measures to achieve pruning: a length criterion, a criterion estimating whether the branch is centered, and two criteria involving the amount of information the branch encodes in the shape. The thresholds for each criterion are set automatically by averaging measures over branches. This method has two main drawbacks with respect to our tubeshaped volumes. First, one criterion is relying on centeredness of the branch. However, certain algorithms cannot guarantee this property (see Fig. 3). Second, using the average of the significance measure might not be relevant in the context of varyingdiameter shapes.

Another post-processing step consists in recentering the skeleton. Points are not necessarily centered in the branch (see Fig. 3a), or branching points differ from one method to another (see [START_REF] Sobiecki | Comparison of curve and surface skeletonization methods for voxel shapes[END_REF] for example). Zigzags in the skeleton impact geometrical measurements leading, for example, to an overestimation of the length. [START_REF] Bradley | Post-processing techniques for making reliable measurements from curveskeletons[END_REF] smooth and recenter skeletons as a post-processing step. First, they compute the set of maximal balls centered at each point of the input skeleton. Then, they compute the set of minimally overlapping balls. The resulting centered and smooth skeleton consists in the centers of the balls linked by splines. Although the presented approach is interesting to make measurements in the context of noisy data, it can discard interesting small geometrical features as well. Bar-Figure 4: Tangent estimation (in red) at two curve points (in yellow) using the 3D λ-MST method [START_REF] Postolski | Tangent estimation along 3d digital curves[END_REF]. The tangent is given as a weighted average of the digital straight segments (in green). Due to irregularities on the curve, the tangent estimated at the point on the right does not correspond to the shape of the curve.

bieri et al. ( 2015) describe an approach to recenter a skeleton inside a shape. The new recentered point correspond to a centroid in the cross-section. The main drawback of this method is the orthogonal plane estimation is very sensitive to irregularities in the skeleton, and junctions are not properly processed.

Orthogonal plane estimation is usually based on the skeleton. The tangent at a point on the curve-skeleton is the normal of the orthogonal plane. Tangent estimation at a point p on a 3D discrete curve can be achieved in a naive way, by computing the vector defined by the difference between two points in the skeleton. The λ-MST estimator described in [START_REF] Postolski | Tangent estimation along 3d digital curves[END_REF] allows to estimate 3D tangent more robustly. It relies on maximal segments (longest digital straight segments) which capture the shape of the curve. The tangent at p corresponds to a weighted average of the direction vectors of all maximal segments passing through p (see Fig. 4). Both the naive and the λ-MST methods are insufficient as they are sensitive to irregularities in the discrete curve. Other methods such as [START_REF] Kang | Vessel crosssection determination based on nonrigid registration and electric field model[END_REF] use non-rigid registration of a deformable cylinder on the acquisition image, and compute orthogonal planes from the axis of this cylinder. However, the fitted cylinder has a constant diameter. Thus, this method is not suited for varying-diameter tubes.

This paper addresses the shortcomings of the afore-mentioned methods, and proposes a new orthogonal plane estimator which is more accurate than existing methods and does not necessarily rely on a curve-skeleton computation for certain types of tubular volumes. Moreover, we introduce automatic post-processing approaches for the curveskeleton, including pruning and recentering which tackle the issues mentioned above. Finally, we present a new algorithm for curve-skeleton extrac-tion on tubular volumes which yields complete and centered skeletons. 2011) on point clouds which was used to estimate normals. The VCM was also defined on sets of voxels in [START_REF] Cuel | Voronoi-based geometry estimator for 3D digital surfaces[END_REF] and was proven to be a reliable tool to estimate the surface normal. This measure was shown to be resilient to Hausdorff noise and to outliers [START_REF] Cuel | Robust geometry estimation using the generalized voronoi covariance measure[END_REF].

Let y be a point in a digital set. The VCM considers the covariance matrix of Voronoi cells in the neighborhood of y. The eigenvector with the largest eigenvalue in this matrix corresponds to the surface normal. Since the Voronoi diagram is global by nature and sensitive to noise and digitization effects, VCM computation is restricted to a small computation window taking into account a local neighborhood. This ensures robustness of the estimator and its multigrid convergence properties [START_REF] Cuel | Voronoi-based geometry estimator for 3D digital surfaces[END_REF].

We recall here the formal definition of the VCM. In the following, we refer to digital points, i.e. centers of pixels in 2D or voxels in 3D, as points. Let us denote by n the dimension, O ⊆ Z n the object, B(o, ρ) the points in the ball of radius ρ and centered at o. The R-offset of O is defined as O R = ∪ o∈O B(o, R). The R-offset bounds the Voronoi cells Vor(.), allowing to extract local information from the Voronoi diagram. Voronoi cells in a neighborhood around a point are taken into account: at a given point y, the computation window is defined as W = B(y, r) ∩ O, where r is called integration radius. The domain of integration for the VCM at y is given as:

DI O (y, r, R) = ∪ w∈W (Vor(w) ∩ O R )
The VCM is the covariance matrix of all vectors between the points in the domain of integration and their respective Voronoi sites (see Fig. 5a). The site for a point x is denoted by p O (x). Finally, the covariance matrix is expressed as: The VCM is defined on any compact set. Figure 5a shows its computation at a point on a 2D digital curve, and Fig. 5b shows its computation on a 2D digital object. Since the Voronoi cell of a point inside the object is reduced to the point itself, these points do not bear influence in the VCM computation. Only points on the border contribute to the VCM.

V O (y, r, R) = x∈DI O (y,r,R) (x -p O (x))(x -p O (x)) T (a) (b)
From this definition, it is clear the choice for the integration radius r, and R is of major importance. An empirical study for the VCM parameters r and R has been led and its results are summarized in Section 3.3.

Orthogonal plane

We propose to use the VCM to estimate orthogonal planes (as described in our work in [START_REF] Grélard | Precise cross-section estimation on tubular organs[END_REF]). The orthogonal plane at a given curve point is defined by its tangent to the curve. The two main elongation directions of a 3D Voronoi cell define the orientation of the orthogonal plane (see Fig. 6a). As explained in Section 3.1, the VCM estimates the shape of the expected Voronoi cell at a point by integrating several digital Voronoi cells. In 2D, the normal at a curve point is given as the largest eigenvector in the VCM. Estimating a precise orthogonal plane on a 3D discrete curve is the same problem as detecting a normal on a 2D curve point, only with a supplementary dimension. The two eigenvectors with the superior eigenvalues in the VCM give the basis for the orthogonal plane. The eigenvector with the lowest eigenvalue corresponds to the tangent to the curve.

Since the definition of the VCM holds for any compact set, the orthogonal plane estimation does not necessarily rely on a curve, and can be used directly from a set of voxels. As an example, the shape of the Voronoi cells on the surface of a cylinder are pictured in Fig. 6b. Let us denote by y the point where an orthogonal plane must be estimated. There is a correspondence between the shape of the Voronoi cells at surface points close to y and the expected orthogonal plane. Indeed, the union of all such Voronoi cells is aligned with the orthogonal plane. As a result, the conclusions drawn from a 3D discrete curve apply for a volume, and the basis of the orthogonal plane is given by the two eigenvectors with the superior eigenvalues. The main challenge in order to estimate orthogonal planes directly from the volume is to define which local points on the surface must be taken into account in the domain of integration.

Automatic setting of the integration radius

Let P(p) be an orthogonal plane at a point p and r the integration radius of the VCM used for the orthogonal plane computation, as described in Section 3.

The integration radius r must be "well-adjusted" to capture the shape of the object locally, while being robust to shape irregularities. A parameter study correlating r to the radius of the tube ρ can be found in an online appendix (F. [START_REF] Grélard | [END_REF]. The results show the integration radius r should be chosen as ρ+2. The parameter study also highlights that r should be adjusted such that the domain of integration contains (a) surface points for which Voronoi cells are aligned with the expected orthogonal plane and (b) the least amount of irrelevant surface points. In the following paragraph, we use a ball as domain of integration in order to integrate all the relevant surface points, and we describe how to automatically obtain the minimal value for its radius r.

The size of the domain of integration is initialized with the error-free Euclidean distance transform value (DT ). Obviously, this value does not correspond to the expected value for r in the case of a tube with elliptic cross-sections, or when the orthogonal plane is computed at a point p near the boundary of the object. The process is the following: an initial orthogonal plane P 0 (p) is computed with r = DT (p). Let C be the cross-section, i.e. the connected component in O ∩ P(p) which contains p. If there is at least one point x in C, such that d(p, x) > r, then r is incremented. This process is repeated until no points in C are found at d > r. It allows to converge towards both the expected orthogonal plane and the minimal value for r to include all the necessary surface points in C.

Results

In this section, the VCM efficiency in computing orthogonal planes is compared to the λ-MST method. These methods have been implemented using the DGtal library (D. [START_REF] Coeurjolly | DGtal: Digital geometry tools and algorithms library[END_REF] and are available along with the data online1 . We compare the methods both on synthetic and real data.

Noisy synthetic data

The goal of this section is to compare our method to the λ-MST estimator on known volumes with altered surfaces. We have generated a slightly curved cylinder with constant diameter (see Fig. 7), and a straight elliptic cylinder with varying minor and major axis values (see Fig. 10a).

Irregularities on the surface of tubular organs are the source of skeleton distortions. In order to generate similar irregularities on the synthetic data we have added some noise on the surfaces of the synthetic object. Noisy versions of the volumes are produced using a simplified version of Kanungo's algorithm [START_REF] Kanungo | A statistical, nonparametric methodology for document degradation model validation[END_REF]. This method adds noise on binary images by switching the value of each voxel, according to its distance d to the object boundary, with a probability α d . Cavities and unconnected noisy voxels are then removed using morphological operators.

The first test to assess our method's robustness consists in extracting geometrical characteristics on the computed 2D orthogonal planes, and in comparing them to known values.

The intersection between each orthogonal plane and the unaltered digital curved cylinder should be a disk. In order to quantify how close to a disk the results are, two features are computed on the 2D shape resulting from the intersection: the area in number of pixels, and the roundness, given as 4A π * a 2 where a is the length of the major axis, and A the area. This value ranges from 0 (line) to 1 (perfect circle).

The corresponding expected values can be computed with the known radius of the curved cylinder (20 pixels). The results are obtained on 93 cross-sections. A sample of a few orthogonal planes are shown on the curved cylinder in Fig. 9. The mean values obtained with the two variants of the VCM method (see Figs. 8a and8b), are closer to the theoretical value than those obtained with λ-MST method on the computed skeleton. Furthermore the difference in standard-deviation between the two methods is significant. The coefficient of variation, defined as the ratio of the standard deviation to the mean, is 29% and 12% for the λ-MST, against 1.5% and 1.1% for the VCM, for the area and roundness respectively. This reflects the high-variability of the λ-MST method, as it finds a substantial number of incorrect orthogonal planes, whereas the results of the proposed method are consistent.

Another evaluation process consists in evaluating the difference between the normal direction of the estimated plane, and the known direction. This test has been performed on a straight elliptic cylinder for which the normal direction at each point of its initial centerline is constant (see Fig. 10a). For each point of the skeleton computed using the method in [START_REF] Lee | Building skeleton models via 3-d medial surface/axis thinning algorithms[END_REF], the angle defect between the computed normal and the expected normal is determined for all methods. Figure 10b shows the results Using VCM computed from both the centerline and the volume yields an angle defect close to zero, with a low standard deviation. The mean value of the angle defect are greater for the λ-MST, and again suffers from high-variability (standard deviation of 22 degrees against 1.1 and 0.46 degrees for the VCM on the curve and on the volume respectively). The λ-MST estimator does not perform well because DSS recognition is sensitive to slight pixel deviation in a curve: in some cases, short DSS are found which means the computed orientation is not representative of the actual tangent.

Real data

Estimation from the curve. We aim at evaluating our orthogonal plane estimator on real data, namely an airway-tree dataset from which a curve-skeleton is extracted. The bronchi have been segmented manually, and the skeleton was computed on the resulting volume with the method of [START_REF] Lee | Building skeleton models via 3-d medial surface/axis thinning algorithms[END_REF].

Parts of the skeleton are impacted by irregularities of the surface. This leads to some incorrect orthogonal planes estimation with the λ-MST method (see Fig. 11a) for these points. On the contrary, similarly as what was observed on synthetic data, the VCM method gives consistent results which are not affected by slight distortions in the skeleton (see Fig. 11b).

In junctions, there are various skeletal branches and this impacts the local computation. Thus, the orthogonal plane estimation is done branch by branch. Furthermore, the first planes estimated after a branching part are not of interest since not only do they contain a section from a tubular part, but they also contain a large part of the second beginning branch (see Fig. 11d).

Both results on synthetic and real data show our method is more precise than state-of-the-art algorithms to estimate orthogonal planes on a curve.

Estimation from the volume. Computation of orthogonal planes from the volume is of major interest since it does not rely on a curve-skeleton and its potential defects. Typical state-of-the-art methods are based on tangent estimation on a curve in order to estimate the orthogonal plane. Thus, only the results of our method are displayed here. The radius for the domain of integration was computed automatically according to the method described in Section 3.3.

Estimating orthogonal planes directly from the volume requires that the input volume is tubular. The orthogonal planes may indeed not be properly estimated in parts where the volume is not locally tubular, in the sense of the domain of integration of the VCM. This is the case of the colon segmentation, in areas where the colon folds (see Fig. 12a). In such areas, the orthogonal planes may not be properly positioned. However, on tubular volumes, such as neurite volumes, the orthogonal planes are properly positioned and very close to the planes estimated from the skeleton (see Fig. 12b). Our method is robust to surface irregularities or noise located on the surface of the volume (see close-up on Fig. 12b). Regarding the airway-tree (see Fig. 12c), the planes are properly oriented in tubular parts as well.

The orthogonal plane estimation from the volume is similar to the one from the skeleton for tubular shapes. This makes the volume estimator very interesting since it is not necessary to compute the curve-skeleton as a pre-processing step. To our knowledge, it is the first orthogonal plane estimator which works directly on the volume and takes into account its geometry.

Pruning

As explained in Section 2, a lot of curveskeletonization approaches on tubular objects produce irrelevant branches originating from surface irregularities. Skeletal branches are maximal sets of connected points in the skeleton which are not branching points (i.e. points having at least three neighbors). Pruning methods aim at keeping relevant branches while removing spurious branches. Spurious branches are not aligned with the tube's axis. In this section, we propose a new significance measure based on orthogonal plane computation to estimate the deviation of a skeleton branch from the tubular object direction. This measure can be used to perform automatic pruning of an existing curve-skeleton.

Method

The curve-skeleton is decomposed into a set of branches, in order to decide, for each branch, 2005), (b) on a neurite, data from the SNEMI3D challenge [START_REF] Kasthuri | Saturated reconstruction of a volume of neocortex[END_REF], with a close-up showing surface irregularities: our method is not sensitive to them; and (c) on an airway-tree.

whether to keep it or prune it. Since the curveskeleton is decomposed into branches, it must satisfy the two following properties : (a) it must be thin, in order to detect branching points based on their neighborhood and (b) it must be connected. These two properties are ensured by a lot of thinning approaches, and general-field algorithms which makes it applicable on a wide variety of skeletonization methods.

Figure 13: Difference between orthogonal planes computed from the curve (green plane) and the volume (pale red plane). On the left, the two orthogonal planes overlap, because the skeleton is aligned with the tube axis. On the right, the plane computed on a spurious branch is highly deviated compared to the one computed on the volume.

The significance measure consists in measuring the angle α formed by the two orthogonal planes defined, on one hand by the branch (a curve) and, on the other hand by the volume. The orthogonal plane normal of an undesired branch is directed at the object's surface, whereas its volumetric counterpart is along the shape. The angle formed by the two planes for spurious branches is large, whereas it is low for meaningful branches (see Fig. 13).

Orthogonal planes on the curve and the volume are computed with the VCM at each point in the branch (see Section 3). In order to obtain a robust estimation of the orthogonal plane on the volume, the integration radius r is automatically set using the method described in Section 3.3. Regarding the computation of the orthogonal plane from the curve, r depends on the branch length. Indeed, a low-resolution version of a given curve should have a lower r, in order to get into account only relevant local information. Experimentally, we choose r as half the branch length since it gives a good representation of the tangent to the curve.

The angle difference α is computed for each point in the branch. The average angle difference for a branch allows to discard it or keep it according to a threshold. Unlike methods presented in Section 2, setting a global threshold is not a problem in this case since the significance measure is independent of the local diameter of the shape. The choice for the threshold is discussed in Section 4.2.

All branches in the skeleton are analyzed during the pruning procedure. Only those having a significance measure above a given threshold and preserving the same number of connected components as in the initial skeleton when removed are deleted. The skeleton is modified and updated when a branch is deleted. Computation of the significance measure is done on each branch of the updated skeleton until no branch can be removed. This ensures that the pruned skeleton is not disconnected, and that all spurious branches across the skeleton are deleted.

Table 1: Comparison of our pruning method to the state-of-the art method in [START_REF] Serino | Pruning the 3d curve skeleton[END_REF]. Four measures are used to assess the completeness and sensitivity of our method compared to the state-of-the-art. The number of spurious branches deleted by the existing method and not deleted by ours (denoted by Err VCM1 ); the number of spurious branches deleted by our method but not deleted by the existing one (Err ST1 ); the number of relevant branches deleted by our method but preserved by the other (Err VCM2 ); and the number of relevant branches deleted by the existing one, but preserved by ours (Err ST2 ). These measures are computed on the two volumes of Fig. 14 and six supplementary volumes (results not shown): AW1, AW2 and AW3 (other airway-tree segmentations), PA (pelvic arteries segmentation) and N1 and N2 (neurite volumes from the SNEMI3D challenge [START_REF] Kasthuri | Saturated reconstruction of a volume of neocortex[END_REF]).

Volumes Fig. 14a

Results

Skeletons have been generated using the algorithm described by [START_REF] Cornea | Computing hierarchical curve-skeletons of 3D objects[END_REF] with parameters chosen specifically to yield spurious branches. Then, we apply our pruning scheme on the skeletons. In the following, interesting experimental values for the angle threshold have been found to be in the range of 20 degrees to 35 degrees. The optimal threshold depends on the minimum angle spurious branches make with the tube's axis. In practice, since the significance measure is computed for all the points in the volume, users may adjust the threshold on-the-fly and estimate which branches are removed in interactive time. Results on two varying-diameter tubular objects are presented in Fig. 14.

The quality of a pruning approach is determined by the number of spurious branches it deletes while preserving the maximum number of relevant branches. Out of the eight volumes in Table 1, only one spurious branch is not deleted. This particular spurious branch is aligned with the tube's axis, thus it cannot be discarded using our method. In addition, spurious branches intersecting with one another and forming a cycle can also be deleted.

Our pruning approach also preserves relevant branches. Only incomplete peripheral branches might get removed, because there is not enough information to estimate the significance measure robustly (see close-up in Fig. 14b). These branches are located in a junction and thus are not suitable to study tubular organs, so their removal does not impede further analysis.

Furthermore, our method is compared to the pruning approach introduced in Serino and Sanniti di Baja (2014). Results in Table 1 show our method is more complete as it deletes more spurious branches than the existing method (see Err ST1 and low values for Err VCM1 ). Moreover, our method is more sensitive, because it preserves more relevant branches (see Err VCM2 and Err ST2 ).

Recentering

Recentering is another post-processing step to improve curve-skeletons. Our recentering approach is also based on orthogonal plane estimation. We handle tubular parts and non-tubular parts, such as junctions, separately.

Method

In a tubular part, new recentered points are defined as the centers of mass of the cross-sections computed at each point of the input non-centered skeleton. Each skeletal branch is processed independently. The recentered skeleton points might be disconnected due to the fact that points are shifted from the input skeleton during the recentering procedure. However, the initial skeleton is connected, which implies that computed orthogonal planes are close from one another. As a result, recentered points remain close (generally at most one point apart from each other). For this reason, each pair of disconnected points is linked by a small digital segment, without altering the centeredness of the skeleton.

Junctions

Junctions are not tubular, thus they must be processed independently. The goal is to decompose a junction into tubular sub-volumes (see Fig. 15b). The junction can be regarded as the division of a parent tube into children tubes. Each child should be associated with its parent: the union of a child branch with its parent branch can be seen as a tubular object. Shape decomposition is required to achieve this. In the following, only junctions consisting of three tubes are described for the sake of simplicity.

(a) (b) From a given skeleton, we aim at finding the corresponding shape decomposition in order to create sub-volumes. Various state-of-the-art methods ex-ist for this purpose, such as the work presented in [START_REF] Palágyi | Quantitative analysis of pulmonary airway tree structures[END_REF], which consists in assigning a different label for each branch of the skeleton. Then, the voxels in the volume are given the label of the closest skeletal point. In the context of recentering, the produced cuts do not yield tubular parts (see Fig. 15a). As a result, orthogonal planes cannot be accurately computed with the VCM on the sub-volumes from this decomposition.

Thus, we propose a new shape decomposition scheme. The idea is to find orthogonal planes near a junction which delineate a part from another (see Fig. 15b). These planes are referred to as cutting planes in the following paragraphs.

The input skeleton must possess the same properties as those presented in Section 4.1 (i.e. thin and connected), and have a correspondence with the shape, which means it must be pruned. Shape decomposition is determined by junctions, because they correspond to zones where different parts meet. In a pruned skeleton, branching points are in correspondence with junction areas in the volume, and interesting cutting planes are in the neighborhood of branching points in the skeleton.

Orthogonal planes are used as cutting planes to decompose the shape. Only a few orthogonal planes are interesting for shape decomposition.

For each edge around a branching point b, the cutting plane should be positioned at a point where it delineates the tubular branch from the junction area. Orthogonal planes computed near branching points in a child branch intersect the other child branch (see Fig. 16a). In other words, the area of the cross-section is much larger in the junction than in the child branch (see Fig. 16b). The cutting point corresponds to the point where the area variation is the largest (see Fig. 16c). The area A(p) is computed as the number of points in the intersection between the orthogonal plane at p and the volume. To study the area variation, we study the factor σ(x) along the edge E.

σ(x) = max y∈N E (x) A(y) A(x)
where N E (x) corresponds to the direct backward and forward neighbors of x in the curve E.

Then the cutting point c in the edge E is defined as c = max x∈E σ(x). In a junction, only the cutting points of the child branches and their corresponding cutting planes are considered.

The cutting planes provided by the decomposition allow to separate each part but are not proper delineations to create tubular sub-volumes. Indeed, the union of two parts is not tubular in the junction area (see Fig. 17a, in green). The idea is to align the plane with the axis of the other child tube, so as to "extend" the child tube in the parent tube. The direction of the axis is given by the orthogo- Sub-volume (in green) obtained after rotating the plane P 1 around the center (in orange) so that the new cutting plane P 1 is orthogonal to the plane P 2 located in the other child branch.

nal plane normal (i.e. the tangent to the curve) computed in the other child tube. Thus, the first cutting plane P 1 , defined by the normal n1 , is rotated so that it is orthogonal to the other cutting plane P 2 , defined by the normal n2 . The direction of the rotation axis is given as the cross-product of the two normals. Then, the cutting plane normal is

n 1 = (n 1 ∧ n2 )∧ n2 .
The rotation axis is placed on a point in the cross-section defined by P 1 and which minimizes the distance to the cross-section defined by P 2 . After rotation of the planes, two tubular subvolumes S 1 and S 2 can be isolated for each branching point b (see Fig. 17b). The skeleton is recentered by computing the centers of mass in each subvolume.

Results

This evaluation aims at verifying whether our recentering scheme produces a recentered skeleton. Tubular volumes with junctions have been generated using a parametric curve. This parametric curve is the theoretical skeleton. Initial noncentered skeletons were generated using the algorithm described in [START_REF] Palágyi | Directional 3d thinning using 8 subiterations[END_REF]. Then, we use the Hausdorff distance d H (C, C ) defined as the maximal distance between a point in C and its nearest point in C. It allows to measure the dissimilarity between two curves. Here, two Hausdorff distances are of interest: d H (T, I) which is the Hausdorff distance between the theoretical skeleton T and the initial non-centered skeleton

(T, I R ) 1 1 1 √ 2 √ 2
I; and d H (T, I R ) which is the Hausdorff distance between the theoretical skeleton and the skeleton I R produced by our recentering method. Volumes are comprised of simple "Y"-junctions consisting in three tubes glued together with varying division angles (0.45π; 0.40π; 0.35π; 0.30π; 0.25π). The results are presented in Table 2 and show that the recentering procedure produces recentered skeleton from initial non-centered skeleton. Indeed, computed skeletons are at most at √ 2 from theoretical skeletons, which means the theoretical and recentered skeletons points all share a voxel edge. Moreover, our technique is independent on the junction angle.

The recentering procedure was tested on noncentered skeletons generated using the thinning algorithm described in [START_REF] Palágyi | Directional 3d thinning using 8 subiterations[END_REF]. Results on real data (airway-tree segmentation) show the skeleton is centered in junctions, and zigzags are properly smoothed (see Fig. 18).

Curve-skeleton computation

In this section, we show how to use orthogonal planes to compute a curve-skeleton on certain tubular volumes. Tubular parts and junctions are handled in two separate passes in our algorithm. The main idea of our algorithm is that the curveskeleton of the object O is the set of centers of mass of the cross-sections defined by the orthogonal planes. The curve-skeleton is obtained by tracking the centers of mass iteratively, as described in our work in [START_REF] Grélard | Centerlines of tubular volumes based on orthogonal plane estimation[END_REF].

Tracking algorithm

Tracking. It is computationally expensive and redundant to go through all points in the volume. A cross-section contains a unique curve-skeleton point so once this point is added to the skeleton, all the remaining points in the cross-section do not need to be processed further. These points are stored in a set M (marked points), and are not used in the remainder of the tracking procedure.

The starting point p 0 of the algorithm is the point in the volume which has the highest DT value. After the orthogonal plane computation at a point p i we obtain a center of mass g i in the connected component defined by the orthogonal plane and which contains p i (i.e. the cross-section). Then, we prop-agate to the next point p i+1 in the 26-neighborhood of g i with the normalized plane normal n P(pi) :

p i+1 = g i + n P(pi) until g i + n P(pi) ∈ M .
Since g i + n P(pi) ≈ g i+1 in regular tubular portions of the object, propagation is done either on the next center of mass or on its neighborhood (see Fig. 19). This ensures the algorithm prioritizes points which are well-centered and for which the orthogonal plane normal at p i is close to the plane normal at p i+1 . This process is done for the two normal directions (g i + n P(pi) and g i -n P(pi) ). Once points have been tracked in both directions and the next tracked point is in M , then a new point p i+1 is chosen as the one with the highest DT value.

In junctions, orthogonal planes are not properly defined. The following paragraphs describe how to specifically detect and handle junctions.

Junction detection

An approach towards junction detection is based on analyzing the local topology of the object, and has been described in [START_REF] Xiong | Tracking the motion trajectories of junction structures in 4d ct images of the lung[END_REF]. The authors define the spherical shell intersection (SSI) as the difference between two concentric balls intersected with the volume. The number of connected components of the SSI is equal to two in a regular tube and to the number of involved tubes in a junction (see Fig. 20a). Balls are defined using geodesic distance so as to consider the object locally. The authors show the choice for both the ball radii is crucial to find junctions accurately. Indeed, if the inner ball radius is too small, it might consider a junction area as a regular tube, and if it is too large, it might discard small but relevant protrusions in the volume. The ball radii depend on the local scale of the tube. In [START_REF] Xiong | Tracking the motion trajectories of junction structures in 4d ct images of the lung[END_REF], the authors set the radii empirically based on the distance transform value at the ball center. It makes sense to rely on the distance transform values at centered points for tubes with circular cross-sections because they are equal or close to the tube radius. However, this does not hold for tubes with deformations or elliptical cross-sections. Thus, we propose a way to set both the radii in a manner which makes sense regardless of the shape of the tube. The ball radius is set as the distance between the two farthest points in the cross-section.

All points p satisfying SSI(p) ≥ 3 are called 3shell points and are not added to the skeleton.

Junction processing

At this point, the skeleton is complete in tubes: for each tubular part in the volume, there is an associated part in the skeleton. Each skeletal part (i.e. connected component in the skeleton) must be connected to another in the corresponding junction area. This problem is very close to the one consisting in recentering an existing skeleton in junctions (see Section 5.2). However, skeletal part endpoints might not correspond to cutting points (they are located deeper inside the tube). As a result, if we use the recentering approach, the cutting plane would not be properly positioned, which means the resulting sub-volume would not be tubular and the resulting computed skeleton would not be centered.

Thus, skeleton points in the junction are computed by linking skeletal part endpoints with Bezier curves. First, endpoints must be grouped together if they belong to the same junction area. A junction area J is defined as a maximal connected set of points closer to a 3-shell point than to a skeleton point (see Fig. 20b). We have to connect the endpoints which belong to the same junction area and which are not terminal points, i.e. for which SSI > 1.

We define the reference branch as the skeletal branch to which other branches in the same junction must be linked (see Fig. 20c). The reference endpoint e ref is the point which normal direction is contrary to the most of the other endpoint normals, i.e.

for endpoints belonging to the same group G, e

ref = arg max e1∈G (card({e 2 ∈ G | n e1 • n e2 < 0})).
Finally, we link each point e belonging to a same group G to e ref using Bezier curves. The two control points in the Bezier curve are defined by the orthogonal plane normals at both endpoints. It is also possible to use our recentering procedure after having linked endpoints with Bezier curves to improve the skeleton in junctions if necessary.

All the Bezier curves meet at the reference endpoint, so they are post-processed so as to avoid creating thick parts. For a junction, all voxels of the different Bezier curves are added in parallel, from the endpoint e = e ref .

Once the Bezier curves intersect, a new common Bezier curve is computed, starting from the intersection point and ending at e ref . This ensures thinness in junctions.

Results

Our method's efficiency is evaluated by comparison to ground-truth curve-skeletons and to state-of- the skeleton is computed by our method. Our method is not sensitive to irregularities and major diameter variation. Our skeletons contain no faulty branches.

the-art methods namely euclidean skeleton thinning approach [START_REF] Couprie | Discrete bisector function and Euclidean skeleton in 2D and 3D[END_REF] and the potential field method [START_REF] Cornea | Computing hierarchical curve-skeletons of 3D objects[END_REF] presented in Section 2. Although the limits of both these methods have been discussed, they satisfy interesting properties and are commonly used in our field of application. Our method is tested on both synthetic and real data.

Synthetic data

Various simple tubular volumes have been generated using a parametric curve: a cylinder, a cylinder with varying-diameter, referred to as deformed cylinder, a junction, and a curved tube (see Table 3). Four noisy versions of these tubes were generated as well using a simplified version of Kanungo's method (see Section 3.4). For each volume, the expected curve-skeleton T is the parametric curve used to generate the volume. The expected skeleton is compared to the computed skeletons C using the Hausdorff distance d H . The skeleton is centered if d H (T, C) ≤ √ 3, that is to say if the points in the computed skeleton are neighbors of the points in the expected skeleton. Otherwise, the skeleton is thick or contains faulty branches. Results are shown in Table 3.

Various parameters for the state-of-the art methods are tested in order to obtain the lowest number of faulty branches and the most complete skeleton. Only results with the best parameters are shown here.

For all non-noisy volumes, all methods produce centered skeletons (d H (T, C) ≤ √ 3, results not et al., 2007), PF [START_REF] Cornea | Computing hierarchical curve-skeletons of 3D objects[END_REF] and ours) on different noisy volumes. Regarding noisy volumes, our method produces a skeleton which is the closest to the initial skeleton. The thinning method of [START_REF] Couprie | Discrete bisector function and Euclidean skeleton in 2D and 3D[END_REF] produces faulty branches which cannot be deleted regardless of the parameter values. As a result, the Hausdorff distance for these skeletons is large. Our method produces results which are close to those of the potential field method of [START_REF] Cornea | Computing hierarchical curve-skeletons of 3D objects[END_REF]. Indeed, skeletons produced by both methods are centered. A visual inspection on the junction volume shows our skeleton contains less irregularities than the potential field skeleton and that it is more complete (i.e. there are skeletal points at the extremities of the volume). Moreover, it is not sensitive to irregularities and to major diameter variation. Results for the various methods on synthetic data are available in the online appendix (F. [START_REF] Grélard | [END_REF].

Data

Real data

Our method was designed to work on tubular organs. We applied it applied on neurites, airwaytrees and arteries.

Figure 21 shows the results on various volumes of neurites. Compared to the state-of-theart method, our skeleton is thin and produces no faulty branches. Figures 22a, and22b show the resulting skeletons on two different airway-trees. Figure 22c shows resulting skeletons on pelvic arteries, which also contain junctions. A general visual inspection shows the skeleton does not have the defects of the state-of-the-art methods illustrated in Section 2. Close-ups show the resulting skeleton can capture local information in junctions, as well as in small branches and be exempt of faulty branches throughout the majority of the volume. Moreover, our skeleton is well centered and does not suffer from major irregularities on the surface (see Fig. 22b). In addition to the fact the skeleton is obtained automatically, these properties are interesting in light of the various applications mentioned in Section 1. More results are available online showing our method performs better than state-of-theart methods (F. [START_REF] Grélard | [END_REF]. One limit of our approach is that it does not work for organs such as the colon which can be locally non-tubular, as already explained in Section 3.4.2.

Conclusion and prospects

In this article, we have presented several new and original approaches for tubular organ analysis. We presented a method to estimate orthogonal planes is more accurate than state-of-the-art methods and is not sensitive to irregularities on the surface. Orthogonal planes can not only be estimated from a curve, but also directly from a set of voxels. This has helped defining new tools for the study of tubular organs. All these tools, along with the source code, are available online 1 . We presented an automatic pruning method removes close to all the spurious branches, and the associated significance measure is not dependent on the local scale of the shape, which makes it suitable for varying-diameter tubular organs. Pruned skeletons are recentered, which makes virtual endoscopy and geometrical measurements more accurate. Finally, we use the centers of mass of the cross-sections to compute our own curve-skeleton. The resulting skeleton is centered, thin and complete, which is critical to capture the shape of varying-diameter tubes aptly. These methods could be used in an automated workflow in order to make geometrical measurements on a large tubular organ database. The results could eventually help establish the relationship between the structure of the organ and pathologies, and better understand its physiology.

Acknowledgements

This work received the financial support from the Laboratory of Excellence ANR-10-LABX-57 (TRAIL). for comparison with state-of-the-art methods. The skeleton is complete, that is to say there is a branch in the skeleton for each bronchus in the volume. This is true even in cases where bronchi have low diameter (cases where only the skeleton in red is visible). Moreover, junctions are properly processed through the second pass of our algorithm. Our skeleton is not sensitive to irregularities on the surface of the volume (see close-up on the volume in Fig. 22b).

Figure 1 :

 1 Figure 1: Neurite volume from the SNEMI3D dataset[START_REF] Kasthuri | Saturated reconstruction of a volume of neocortex[END_REF] (in gray) and skeleton computed by the method in[START_REF] Lee | Building skeleton models via 3-d medial surface/axis thinning algorithms[END_REF] (in red). Orthogonal planes are estimated using the λ-MST method of[START_REF] Postolski | Tangent estimation along 3d digital curves[END_REF] (in blue). The skeleton contains irregularities and spurious branches (see green points in the close-ups) which results in a deviated estimated tangent. As a result, orthogonal plane estimation is impacted.

Figure 2 :

 2 Figure2: (a) Skeleton (in red) computed on a airway-tree segmentation using the thinning method described in[START_REF] Couprie | Discrete bisector function and Euclidean skeleton in 2D and 3D[END_REF]: parameter values were chosen to obtain the most complete skeleton with the least amount of faulty branches, but there are some regardless, see close-up view. (b) (c) Skeleton computed using potential field[START_REF] Cornea | Computing hierarchical curve-skeletons of 3D objects[END_REF]: two levels of hierarchy on a skeleton obtained with two different parameter values. Obtaining a complete skeleton containing all the branches in the volume as well as no faulty branches is not possible.

Figure 3 :

 3 Figure 3: Pruned branches (in yellow) of a skeleton computed using Cornea et al. (2005) on tubular volumes (in gray). The pruning method is described in Serino and Sanniti di Baja (2014). Their centeredness criterion leads to errors (see circled areas): (a) some relevant branches are pruned and (b) some spurious branches are kept.

3.

  Orthogonal plane estimation 3.1. Voronoi Covariance Measure All of our contributions are based on the Voronoi Covariance Measure (VCM). It is a measure, first introduced in Mérigot et al. (

Figure 5 :

 5 Figure 5: Digital VCM computation at point y on (a) a discrete 2D curve O and (b) a 2D object O. Solid lines (in gray and orange) show the link between points and their associated Voronoi sites. All the lines with the same Voronoi site define a Voronoi cell. (a) The VCM is the covariance measure of the Voronoi cells in the domain of integration DI O (y) (associated lines in orange). (b) Voronoi cells are null for points inside the object O, such that only border points have a Voronoi cell. Orange lines are integrated for the VCM at y.

Figure 6 :

 6 Figure 6: Correspondence between Voronoi cells and orthogonal planes. (a) Voronoi cells (in gray) computed at curve points are aligned with orthogonal planes (in blue). (b) VCM computation on a volume (a cylinder) at a point y. Similarly as in the 2D case (see Fig. 5b), only surface points x have a non-null Voronoi cell Vor(x). The Voronoi cells Vor(x) are aligned with the orthogonal plane at y.

Figure 7 :

 7 Figure 7: (a) Noisy tubular volume (in gray) generated fitting balls of constant radius at each point of the digitized initial centerline (in red). (b) Computed skeleton using the method in Lee et al. (1994) containing irregularities.

  Figure 8: (a) Area in pixels and (b) roundness mean values for all the orthogonal planes found for the noisy tube volume of Fig. 7. The λ-MST on the computed skeleton (white) yields results with high-variability, whereas the λ-MST method applied on the initial centerline (L-MST on IC) and the two VCM methods applied on the computed skeleton and the volume respectively are consistent with the theoretical values.

Figure 9 :

 9 Figure 9: Examples of orthogonal planes (in blue) obtained on a noisy tube-like volume with the volumetric VCM method.

Figure 10 :

 10 Figure 10: (a) Noisy elliptic cynlinder with varying minor and major axes, generated along a centerline with constant normal direction. The computed skeleton (Lee et al., 1994) is shown in red. (b) The angle defect (in degrees) between the computed normal and the expected normal shows our method outperforms the λ-MST tangent estimator.

Figure 11 :

 11 Figure 11: Orthogonal planes computed on the skeleton of an airway-tree with (a)-(c) λ-MST and (b)-(d) VCM. The orthogonal planes computed by the λ-MST method are impacted by irregularities in the curve, those computed by our method are not.

Figure 12 :

 12 Figure 12: Orthogonal planes (in blue) computed on various tubular volumes with our method. (a) On a colon, data courtesy of Cornea et al. (2005), (b) on a neurite, data from the SNEMI3D challenge[START_REF] Kasthuri | Saturated reconstruction of a volume of neocortex[END_REF], with a close-up showing surface irregularities: our method is not sensitive to them; and (c) on an airway-tree.

Figure 14

 14 Figure 14: (a) and (c) Initial skeletons generated using Cornea et al. (2005) on various varying-diameter tubular objects. (b) and (d) Resulting pruned skeletons containing little to no faulty branches. Close-up in Fig.14bshows examples of deleted relevant branches (in cyan), however, these branches are not complete and not suitable to study the tubular organ.

Figure 15 :

 15 Figure 15: Methods to create sub-volumes. (a) Shape decomposition obtained using Palágyi et al. (2006). The completeness of the decomposition is dependent on the skeleton. The close-up (black circle) highlights that the decomposition does not always provide tubular parts (see dark green part). (b) Computation of a sub-volume. The cutting plane (in green) isolates the sub-volume, consisting in a parent and a child branch (dark grey), from another child (light grey).

Figure 16

 16 Figure 16: (a) Computation of orthogonal planes at two points on the same edge. The number of volume points (i.e. the area) in the cross-section varies a lot in junctions. (b) Typical area variation (y-axis) observed along an edge. Points are traversed in an ordered manner from a branching point to the end-point in the skeleton daughter branch. (c) Corresponding σ values for the area values of Fig. 16b. The cutting point is the point where the σ value is the highest (shown by an arrow).

Figure 17 :

 17 Figure17: (a) Sub-volume (in green) computed using orthogonal planes on an input skeleton. Using the decomposition scheme does not allow to recenter skeleton points. (b) Sub-volume (in green) obtained after rotating the plane P 1 around the center (in orange) so that the new cutting plane P 1 is orthogonal to the plane P 2 located in the other child branch.

Figure 18 :

 18 Figure 18: (a) and (c) Initial non-centered skeletons generated using Palágyi and Kuba (1999) on airway-tree segmentations. (b) and (d) Corresponding resulting recentered skeletons.

Figure 19 :

 19 Figure19: Tracking performed by our algorithm. From a determined center of mass g 0 , we propagate to the next tracked point p 1 thanks to the plane normal n P(p 0 ) .

Figure 20 :

 20 Figure 20: Junction detection and junction processing. (a) Spherical shell intersection (SSI, in dark green) allow to detect junctions. The number of connected components of the shells is equal to two in a tubular part (p T ) and to three in a junction (p J ). (b) Method to group endpoints from skeletal parts (in red) in junctions in order to link them. Junction areas J 1 , J 2 are groups of points which are closer to a 3-shell point (squares) than to a skeletal point (red lines). We group skeletal part endpoints together based on their distance to a junction area. (c) Determining the reference endpoint e ref to which other endpoints must be linked is done by comparison of the direction of the normal with the other endpoints' orthogonal plane normals. The normals also allow to determine the position of the control points (not shown) and trace the Bezier curve (in green). Here, only one endpoint is linked to the reference endpoint for sake of clarity.

Figure 21 :

 21 Figure 21: Skeletons (in red) on neurite volumes (in gray) from the dataset of the SNEMI3D challenge (Kasthuri et al., 2015). On the left, (Figs. (a) and (c)) the skeleton is computed by the thinning method of Couprie et al. (2007). On the right, (Figs. (b) and (d)) the skeleton is computed by our method. Our method is not sensitive to irregularities and major diameter variation. Our skeletons contain no faulty branches.

Figure 22

 22 Figure 22: (a) (b) Centerline (in red) extraction on two different airway tree volumes (in gray) and (c) pelvic arteries acquired in MRA. The volume for the pelvis arteries is cropped due to the large image size. Fig. 22a corresponds to the Fig. 2 in Section 2for comparison with state-of-the-art methods. The skeleton is complete, that is to say there is a branch in the skeleton for each bronchus in the volume. This is true even in cases where bronchi have low diameter (cases where only the skeleton in red is visible). Moreover, junctions are properly processed through the second pass of our algorithm. Our skeleton is not sensitive to irregularities on the surface of the volume (see close-up on the volume in Fig.22b).

Table 2 :

 2 Hausdorff distances between (middle row) a theoretical centered skeleton T and an initial non-centered skeleton I and (bottom row) between T and a skeleton I R computed by our recentering procedure on various "Y"-shaped volumes.

		0.45π 0.40π 0.35π 0.30π 0.25π
	d H (T, I)	3.74	2.82	3.74	4.47	5.91
	d H

Table 3 :

 3 Hausdorff distance for each method (Thin (Couprie

https://github.com/fgrelard/ OrthogonalPlaneBasedTools


Throughout Europe, urban authorities have gained more responsibilities and resources. They seek to maximise their effectiveness by fostering new forms of cooperation among a large variety of actors and interests, in order to address urban problems and set their own policy priorities.

Why have some cities in Europe been able to achieve a shift from the car-oriented city towards more sustainable modes of transport? How can other parts of Europe and the EuroMed draw lessons from CREATE?

In order to answer these fundamental questions, the CREATE (Congestion Reduction in Europe, Advancing Transport Efficiency) project's main objective is to reduce road congestion in European cities by encouraging a switch from cars to more sustainable modes of transport. Among other valuable lessons, the CREATE project shows how strong urban governance contributes to the shift towards the liveable city.

Find

  out more about the project: www.create-mobility.eu This part of the study was coordinated from Sciences Po, Centre d'études européennes in Paris. For more information, contact Dr. Charlotte Halpern: charlotte.halpern@sciencespo.fr CREATE has received funding from the European Union's Horizon2020 research and innovation programme under grant agreement N. 636573.

  CHALLENGES » shift from the car-oriented city towards planning the liveable city » transport and mobility policies to ensure social and spatial inclusion » foster transport policy change outside cities THE ROLE OF THE EU » regulations » long-term policy goals » the promotion of knowledge and policy solutions Local governments use different methods to achieve their goals. Each city has a different history and local context. They develop differently and find the best solutions for their local circumstances. According to the local context, policy innovations in transport come from either /or: small-scale experiments, policy innovations, trial and error » promotion of economic incentives » offer new mobility services » develop alternative street uses A COMBINATION OF HIGH-AND LOW-TECH SOLUTIONS by mobilising transport organisations » developing new systems while modernizing old infrastructures GAIN MORE AUTONOMY through capacity building » finances » information and data management » authority » human resources and expertise » communication

citizens and other stakeholders private sector urban governance GOVERNING MOBILITY STRONG URBAN GOVERNANCE FOR LIVEABLE CITIES LESSONS FROM THE CREATE PROJECT

	1	URBAN AUTHORITIES FACE COMMON TRANSPORT CHALLENGES	2	URBAN GOVERNANCE IS CHANGING AND BECOMING STRONGER
		» address contradictory demands
		» seek trade-offs betwen environmental, economic and social
		objectives
		» overcome resistance and protest
		national level	local/metropolitan
			governments
		binding laws and
		funding, policy
	priorities, and large
		projects
		pressure
		and financing of private
		projects


Brigitte Bidegaray-Fesquet

Laurent Fesquet
email: laurent.fesquet@imag.fr

Levels, peaks, slopes. . . which sampling for which purpose?

There are many ways to nonuniformly sample a signal, the most widely used, but not only, way to perform it being level crossing. The targeted application, and also the characteristics of the signal, can lead to various choices of samplings, to ensure a good representation of the input signal or robustness with respect to noise.

I. INTRODUCTION

Nonuniform samples can stem from an external source, and in this situation it is usually considered as a drawback for the subsequent processing. It is so when then sampler is not always available, or when it is perturbed (jitter, missing samples. . . ).

Here we consider the case when nonuniformity is wanted, first as a way to compress the signal. This reduction is an effective way to reduce power consumption in electronic systems since there exists a direct relation between activity and dynamic power. Considering today's system power requirements, mitigating the system activity is becoming an important issue for autonomous communicating objects or mobile systems. Indeed, less data to process, store and transmit is the way to decrease power consumption. It can be processed on systems that are event-based and naturally produce and treat nonuniform samples. Another characteristics of such a setting for embedded systems can also be a lower electromagnetic emission and also a stronger protection to spying attacks since the signature of asynchronous operations is less clear than synchronous ones.

Once uniformity is discarded, they are many ways to be nonuniform, more than ways to be uniform. Apart from the less regular mathematical structure, this variety is one reason why it is more difficult to give general results and propose generic algorithms for nonuniform samples, since the way samples have been captured also carries information.

We explore here a few classical techniques to do this, assuming we know the application and general features of the signal (eg. bandwidth, range, or statistical regularity), but no more. We do not address here the case when the system could fix levels or bounds from a deeper knowledge of the signal, or adapt them during the processing through a learning procedure. An already well studied sampling from the point of view of applications is level crossing sampling [START_REF] Marvasti | Nonuniform Sampling. Theory and Practice, ser. Information Technology: Transmission, Processing and Storage[END_REF]- [START_REF] Fontaine | Filtrage de signaux à échantillonnage irrégulier[END_REF] but a lot remains to do from a theoretical point of view [START_REF] Bidégaray-Fesquet | Data driven sampling of oscillating signals[END_REF]. We can use instead peak sampling, or combine both [START_REF] Greitans | Combined peak and level-crossing sampling scheme[END_REF]. Other systems do not focus on values of the system but on their variations. This is the case of derivative level-crossing sampling which is described in [START_REF] Martínez-Nuevo | Derivative level-crossing sampling[END_REF]. Here we explore what we call slope sampling, for which the sampling criterion is some Lebesgue norm distance to an oversampled signal and for which samples are taken when the slope varies significantly.

II. VARIOUS TYPES OF SAMPLING

In this section we discuss various types of samplings and fix the notations. To show the samples, we first choose a short fragment of an ECG signal, which is a sporadic and noisy signal.

A. Level crossing

Since we do not want to adapt here the sampling to the signal, we only suppose that we know the range of the signal and consider a predefined set of levels that encompasses this range. Each time the continuous signal crosses one of these levels, a sample is taken together with the time elapsed since the previous sample was taken. It forms the couple (x n , δtx n ), where x n and δtx n are respectively the amplitude of the signal at the nth sampling time and the time delay (see Figure 1). In fact, this approach is dual to the classical sampling that captures uniformly in time the samples. Indeed, in the classical case, the instants are perfectly known and the amplitudes are quantized, while, in nonuniform sampling, the instants are quantized and the amplitudes are perfectly known. This has been studied for several years from the theoretical point of view [START_REF] Mark | A nonuniform sampling approach to data compression[END_REF]- [START_REF] Sayiner | A level-crossing sampling scheme for A/D conversion[END_REF] and practical implementations have been proposed several years later [10]- [START_REF]Asynchronous level crossing analog to digital converters[END_REF]. Note that all the discussion in this paper about level crossing sampling will also apply to send-on-delta sampling which only differs in a global shift of the levels (when they are uniformly distributed) and a decimation of successive samples of equal amplitude.

The main advantage of such a sampling scheme is its ability to adapt the number of samples to the signal variations. In the active signal parts, the signal often crosses the levels inducing a lot of samples, whereas, in the silence parts, the signal stays for a while between two levels and no sample is taken. Therefore, for smoothly varying or sporadic signals, the number of samples can drastically be reduced compared to the classical sampling technique [START_REF] Mark | A nonuniform sampling approach to data compression[END_REF], [START_REF] Aeschlimann | Asynchronous FIR filters: Towards a new digital processing chain[END_REF], [START_REF] Qaisar | Adaptive rate filtering for a signal driven sampling scheme[END_REF]. This is demonstrated on our ECG fragment in Figure 2. The next step after sampling is often an interpolation to provide (at least theoretically, for algorithm construction) a continuous-time signal. It has been noticed that for many applications zeroth-or first-order hold reconstruction is sufficient and that no practical gain is obtained by using a higher-order time interpolation of samples [START_REF] Grybos | Event-driven sampling of signal with quadratic prediction[END_REF]. We see in Figure 2 that there are many consecutive samples with the same amplitude. In each series of such samples only one or two samples are needed to lead to exactly the same zeroth-or first-order reconstruction and therefore subsequent processing. The samples of Figure 2 can thus be further decimated to remove unneeded same amplitude consecutive samples.

x n-1 x n δtx n
A possible way to capture less samples, and thus decrease as much as possible the activity, is to reduce the number of levels. Nevertheless, applications need to fulfill conditions such as accuracy and relevance of the data, which limits the reduction of the level number. But it is also corollated to the level position and an appropriate positioning of the thresholds can drastically reduce their number. It has been demonstrated by [START_REF] Pelleter | Non-uniform sampling pattern recognition based on atomic decomposition[END_REF], [START_REF] Roa | Designing ultra-low power systems with non-uniform sampling and event-driven logic[END_REF]. In this work, the signal and the application are first studied. Then a method is proposed where a classifier is used for extracting the optimal threshold number and placement. The paper demonstrated that, with only four thresholds optimally placed, it is possible to reduce by two orders of magnitude the number of samples. This shows the potential of such an approach where the level positioning is not necessarily uniform and accurately chosen with an a priori knowledge of the signal and the application.

B. Peak detection

Peak detection consists in taking a sample each time a local minimum or maximum is detected. On the same model as for level crossing, it leads to samples of the form of couples, (y n , δty n ), y n denoting the amplitude of the signal at the nth sampling time and δty n the time delay. This is demonstrated on our ECG fragment in Figure 3.

The analysis of an ECG diagramme is typically the detection of characteristic points and among them peaks, but clearly peak sampling is not well adapted to noisy signals. Although it is not designed as a filter, level crossing sampling naturally filters high frequency components of the signal which have a small amplitude (it is typically the case of noise). Peak sampling does exactly the contrary and it exacerbates noise. Moreover a monotonous signal, even with a large variation, would lead to no sample, since there would be no local extremum.

One advantage of this sampling is however that it only depends on the signal and not on a set of predefined levels like in level crossing sampling or quantum like in send-ondelta sampling.

C. Levels and peaks

Apart from successive equal amplitude samples, we see another feature on the samples of Figure 2: in the very rapid large amplitude peak (named as QRS complex) there is a number of unnecessary samples in the rapidly increasing and rapidly decreasing parts, when keeping the lower and higher samples would lead to a very similar first-order reconstruction. Following [START_REF] Greitans | Combined peak and level-crossing sampling scheme[END_REF] we can combine the advantages of level crossing and peak samplings. We first sample with level crossing. It already yields a first compact and denoised version of the signal. Then we detect local maxima and minima in the level crossing samples and keep only these samples. This once more yields (z n , δtz n ) amplitude-delay couples, that are depicted in Figure 4.

We keep the main drawback of peak detection: if there is no peak in the original signal (monotonous signal, or small peaks that are not detected with the chosen levels) there still is no sample with this technique. But there are some gains, in particular in the relevance of the captured data. Indeed, level crossing is limited in the number of levels in order to be competitive with Nyquist sampling. If it is coupled with peak sampling, we can relax this constraint. For a regular signal, in parts that are predominantly linear, adding more levels would not lead to more samples since linear parts will be decimated.

To that, we have more chance to catch the real value of the peaks, since more levels (in the same range) means less peak clipping.

In the three previous types of sampling, for example in the case of a constant signal, no sample can be taken for a long time, inducing a large cumulative error between the real signal and the sampled one. In particular in level crossing sampling, the error which is considered is mainly a L ∞ -error. In the sequel, we propose a sampling type that takes into account other Lebesgue-norm errors.

In the previous samplings the samples are taken when some event occurs (level-crossing and/or change of slope sign). It could be performed on a continuous-time signal, or on a finely sampled signal that mimics a continuous signal for the purpose of numerical experiments.

D. Slopes

Slope sampling supposes that there exists a pre-existing sampling (regular or not), but original samples are not necessarily stored, they can just be captured and kept or not by the nonuniform slope sampler. As for level-crossing where levels should be predefined, there is also an external parameter to the input signal, which is a tolerance ε, and the choice of a norm. The principle is shown in Figure 5. The point labelled 1 is the first sample of the input signal (or the last sample that has been kept during the process). Then the two following points (labelled 2 and 3 are explored). If the distance between the signal described by points (1,2,3) is ε-close to the signal described only by poins (1,3), then point 2 is discarded. For the L 1 -norm this corresponds to measuring the surface of the dark gray triangle. Then point 4 is explored and if the signal described by points (1,2,3,4) is ε-close to the signal described only by points (1,4), point 3 is also discarded, and so on. When the distance is too big, the last but one point is kept and plays the role of point 1 in the subsequent processing.

This type of sampling, as level crossing, has a tendency to filter small amplitude noise. Depending on the type of input signal it can lead to a drastic reduction of the number of samples. Contrarily to level-crossing, it ensures that even if the signal has only small variations, samples will eventually be taken since small errors accumulate. We denote by (w n , δtw n ) the induced amplitude-delay couples. The samples obtained for the ECG fragment are depicted in Figure 6 for a 3 10 -2 tolerance on the L 1 -norm.

III. NUMBER OF SAMPLES

In regular sampling, the number of samples depends only on the sampling frequency and the total duration of the signal. One of the goals of nonuniform sampling is the reduction of the number of samples, but this number is much less predictable from simple parameters and features of the signal. The example chosen above is a noisy biological signal, which is obviously not the target of peak sampling. The initial signal is regularly sampled with a t s = 5 10 -4 s period, leading to 2801 samples on the fragment of interest. The number of samples and compression percentages corresponding to Figures 2, 3, 4, and 6, are given in Table I.

In [START_REF] Bidégaray-Fesquet | Data driven sampling of oscillating signals[END_REF], theoretical results are obtained on the number of samples that can be obtained through a level crossing sampling procedure with an assumption on the Hölder regularity of the input signal, leading to an upper bound of the number of samples in the worst case. Such an analysis would be impossible for peak sampling, for which there is no compression at all in the worst case. The level and peak sampling is always better in terms of compression than simple level crossing sampling. The worst case in level crossing sampling is obtained when all levels are crossed successively and it is typically the case when the subsequent peak sampling will cancel samples.

For samplings which add extra parameters to the signal itself -i.e. all but peak sampling-the impact of the variation of the parameters can be studied. For example level crossing leads to peak clipping and a small variation in the position of the levels (offset) can make samples appear or disappear. This is one of the reasons why the problem of the right placement of levels is an important issue.

The impact of noise is another issue. We would want the sampling, and especially the number of samples not to depend too much on noise. The case of additive noise for level crossing sampling has been studied in [START_REF] Bidégaray-Fesquet | Noise and offset in the iir filtering of eventbased sampled data[END_REF]. As it can also be seen in Figure 2 a noisy system can cross many times levels with no special information carried by all these nonuniform samples. They can of course be decimated a posteriori. Level and peak sampling does not share this drawback since the decimation is part of the sampling process. We can therefore expect the number of samples to be very stable with respect to noise, if of course it is zero mean and its standard deviation remains below the quantum threshold used for the levels.

We explore further the impact of an additive noise in the next section.

IV. IMPACT OF NOISE ON THE NUMBER OF SAMPLES

To test the impact of noise, we perform the following test on the input signal

f (t) = 1 2 + 1 2 sin(2πt).
This type of signal is not the one targeted by nonuniform sampling, since it is not sporadic and uniform sampling at a Nyquist rate would suit more here. Our aim is to compare the nonuniform methods with respect to noise and to mimic a continuous signal by starting with an oversampled uniform signal. We therefore consider 100 periods of this signal, uniformly sampled at 100 Hz, thus yielding 10000 original samples. We then compare the four sampling methods described in this paper and also add a noise with zero mean and a standard deviation σ. We take the quantum q used for level crossing (with 8 levels) as a unit for σ and explore values of σ = αq for α ∈ [0, 1]. The first test case consists in comparing the number of samples yield by the various types of nonuniform sampling. The results are given in Figure 7. The fact that we use 100 periods of the signal provides some averaged results which are quite similar from one run to another. Thus Figure 7 is quite generic and not specific to a given run.

We see that the number of samples is very constant with slope sampling. At least for σ > .1q, the number of samples is linear for level crossing and level and peak sampling. Of course the number of samples for level and peak sampling is lower than for level crossing sampling since it is a decimation of the later. Peak sampling seems to saturate for high values of σ. For σ = 0 it yields exactly the same number of samples as level and peak sampling, namely 198 samples. It increases rapidly for peak sampling since noise, whatever its standard deviation, induces peaks. Let us notice that no genuine compression of the signal is obtained for level crossing and peak sampling for large values of σ since the number of samples have to be compared with the original value of 10000.

V. IMPACT OF NOISE ON THE RECONSTRUCTION

The same test-case as in the previous section is carried out to study the quality of reconstruction. We now plot the mean over one period L 1 -error between the noisy input signal and the nonuniform signal. To compute this the nonuniform signal is linearly interpolated at the times of the input signal. The results are given in Figure 8, which is once more generic due to the number of studied periods.

For large values of σ all the methods are more or less equivalent in terms of error. Level crossing, peak and slope sampling have a linear behavior away from σ = 0. We expect slope sampling to saturate for even larger values of σ since its criterion is indeed the L 1 error. A little noise is quite favorable for peak sampling, since it gives the opportunity to capture a few samples away from the minima and maxima of the sinusoid. Level and peak sampling does not give very good results for small values of σ. If we now take into account both criteria, quality of reconstruction and number of samples, we see that slope sampling outperforms the others (of course it also depends of the number of levels for level crossing and the tolerance for slope sampling). Indeed we obtain very comparable errors but with much fewer samples for slope sampling.

VI. APPLICATIONS

The above criteria are natural from a mathematical point of view but not always for applications. The number of samples is very important for applications since it triggers the energy that will be necessary to use the samples for subsequent processing. The quality of the reconstruction may not be so important because accuracy is not always the goal of the applications.

In [START_REF] Pelleter | Non-uniform sampling pattern recognition based on atomic decomposition[END_REF], [START_REF] Roa | Designing ultra-low power systems with non-uniform sampling and event-driven logic[END_REF], the nonuniform samples are used for pattern recognition. Then accuracy in the samples is not necessary if the pattern is correctly analyzed. These works show that the main issue is the good placement of the levels for the specific application.

In [START_REF] Qaisar | A signal driven adaptive resolution short-time Fourier transform[END_REF] nonuniform samples are only used to detect the activity of the system. The input signal is sampled nonuniformly. It allows to collect data on time delays. From this information, a classical sampler yielding uniform samples is switched on during the time-windows when delays in the nonuniform samples are small. These uniform samples are then processed classically. Clearly this application does not need to be very precise in terms of signal reconstruction from the non-uniform signal.

Slope sampling has been introduced to test for the main frequency of a signal. After the capture of the samples an evaluation of this frequency was computed. It is a very costly procedure for nonuniform samples since there is no FFT-like algorithm at hand. The very low number of samples given by slope sampling even in the presence of noise was necessary to make this tractable. This last type of application is not designed for embedded systems but for off-line computations. In the case of embedded systems that really need a low number of samples for power saving reasons, another criterion is the capability of implementing the sampling in the systems.

VII. CONCLUSION

We have explored four different ways to nonuniformly sample signals. We have compared them from the point of view of the number of samples and the robustness of this criterion when noise is added, as well as the reconstruction error. We demonstrated the potentiality for reducing the amount of data with level-crossing sampling but also using slopes, peaks and a combination of peak and level-crossing techniques. We also indicated their limitations to be practically employed.

Another point of view can be taken when considering the integration of these sampling schemes in an electronic system. Some of these sampling techniques have already been implemented in hardware [START_REF] Greitans | Combined peak and level-crossing sampling scheme[END_REF], [START_REF] Roa | Designing ultra-low power systems with non-uniform sampling and event-driven logic[END_REF], [START_REF] Beyrouthy | An event-driven fir filter: design and implementation[END_REF] and dedicated analogto-digital converters have also been designed, fabricated and tested by Allier [START_REF] Allier | A new class of asynchronous A/D converters based on time quantization[END_REF], Shell [START_REF] Shell | A continuous-time adc/dsp/dac system with no clock and with activity-dependent power dissipation[END_REF] and Akopyan [10]. Levelcrossing sampling can benefit from the architecture of the existing classical analog-to-digital converters. Even if all the others sampling schemes have not been completely demonstrated in hardware, the implementation of peak detection or slope-crossing is not a real issue and we can have no doubt, that if these sampling techniques are adequate, a circuit will be designed accordingly.

Today's challenge is clearly to develop a framework able to help designers and practitioners to choose the appropriate sampling scheme for their applications. Indeed, it exists plenty of ways of sampling and that can enhance the integrated system quality by reducing their activity and power consumption but also the amount of data to store and to exchange through the internet. Sampling differently offers several supplementary degrees of freedom to the designers and is probably the key for smart autonomous systems requiring extremely low-power as Le Pelleter [START_REF] Roa | Designing ultra-low power systems with non-uniform sampling and event-driven logic[END_REF] demonstrated in his proof of concept applied to a sub-block of a medical implant. Moreover, if we go further, sampling techniques can also be devised for other purposes such as mitigating the circuit electromagnetic emission or security issues.

Fig. 1 .

 1 Fig. 1. Principle and notations for level crossing sampling.

Fig. 2 .

 2 Fig. 2. Demonstration of level crossing sampling on a small ECG fragment.

Fig. 3 .

 3 Fig. 3. Demonstration of peak sampling on a small ECG fragment.

Fig. 4 .

 4 Fig. 4. Demonstration of level and peak sampling on a small ECG fragment.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Principle and notations for "slope sampling".

Fig. 7 .

 7 Fig. 7. Number of samples in terms of noise standard deviation (in q units: σ = αq).

Fig. 8 .

 8 Fig. 8. L 1 error in terms of noise standard deviation (in q units: σ = αq).

TABLE I .

 I NUMBER OF SAMPLES AND COMPRESSION FOR AN ECG

FRAGMENT.

ACKNOWLEDGMENT

This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01).