Olivier Lézoray
email: olivier.lezoray@unicaen.fr

Manifold-based mathematical morphology for graph signal editing of colored images and meshes

This paper presents a framework for morphological processing of graph signals and investigates its usage for colored images and meshes editing tasks. The proposed method enables, with the help of the construction of a manifold-based ordering of color vectors, to define a new representation of graph signals in the form of an ordering of vectors and an index. The ordering relies on three steps: dictionary learning, manifold learning, and out of sample extension. This enables to formulate morphological operators for graphs signals and we demonstrate the performance of the proposed method on various colored image and mesh editing applications (simplification, abstraction, enhancement).

I. INTRODUCTION

Image filtering is recognized as one of the most important operation in image processing. In particular, structurepreserving filtering is an essential operation with a variety of applications in computational photography and image editing [START_REF] Farbman | Edge-preserving decompositions for multi-scale tone and detail manipulation[END_REF]. During the last decade, a lot of structure-preserving smoothing filters have been proposed [START_REF] Zhang | Rolling guidance filter[END_REF], [START_REF] Cho | Bilateral texture filtering[END_REF], [START_REF] Gastal | Domain transform for edge-aware image and video processing[END_REF], [START_REF] Xu | Image smoothing via L 0 gradient minimization[END_REF], and they aim at decomposing an image into prominent structures and fine-scale details, making it easier for subsequent image manipulation such as detail enhancement or visual abstraction. These filters were originally designed for images, but some have been extended to 3D meshes [START_REF] Fleishman | Bilateral mesh denoising[END_REF].

In this paper, we propose to investigate the use of Mathematical Morphology (MM) operators for such editing tasks. This has never been investigated before in literature. Morphological operators are non-linear vector-preserving filters (no new vectors are introduced in the processed image), and therefore they are not subject to the production of halos, which is a common problem in image editing tasks [START_REF] Farbman | Edge-preserving decompositions for multi-scale tone and detail manipulation[END_REF]. The construction of morphological operators relies on complete lattices [START_REF] Ronse | Why mathematical morphology needs complete lattices[END_REF] that impose the need of an ordering relationship between the elements to be processed. If MM is well defined for gray scale functions, there exists no general admitted extension that permits to perform morphological operations on vectors since there is no natural ordering of vectors. We have recently proposed in [START_REF] Lézoray | Complete lattice learning for multivariate mathematical morphology[END_REF] a framework for the construction of complete lattices for any kind of vector data. The latter learns a complete lattice from a modeling of the vectors' manifold. In this paper, we consider these works as a basis for morphological signal editing. In addition, to enable an application of the proposed framework both to colored images and 3D meshes, we consider a formalism based on graph signals [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF].

The paper is organized as follows. In Section II, we introduce a learned ordering of the vectors of a graph signal. In Section III we derive a graph signal representation from the learned ordering and define the associated morphological graph signal operators. Section IV shows how the proposed framework enables the editing of graph signals (with graphs in the form of 2D grid graphs and 3D meshes).

II. MANIFOLD-BASED COMPLETE LATTICE LEARNING

A. Notations

A graph G = (V, E) consists in a finite set V = {v 1 , . . . , v m } of vertices and a finite set E ⊂ V × V of edges. Let (v i , v j ) be the edge of E that connects two vertices v i and v j of V. The notation v i ∼ v j is used to denote two adjacent vertices. The considered graphs in the paper are undirected and unweighted. A graph signal is defined as a function that associates real-valued vectors to vertices of the graph. A graph signal is represented by the mapping f : G → T ⊂ R n where T is a non-empty set of vectors (we will consider only RGB color vectors, i.e., n = 3).

To each vertex v i ∈ G is associated a vector v i = f (v i ). The set T = {v 1 , • • • , v m }
denotes all the vectors associated to all vertices of the graph. We will use the notation T [i] = v i to denote the i-th element of a set.

B. Complete lattice

To perform morphological processing of graph signals, we build an ordering of the vertices' vectors T . Ordering all the values of the set T can be done with the use of an ordering relation within vectors. This amounts to dispose of a complete lattice (T , ≤), a key item for the definition of mathematical morphology operators [START_REF] Ronse | Why mathematical morphology needs complete lattices[END_REF]. Unfortunately there is no universal method for ordering vectorial data [START_REF] Barnett | The ordering of multivariate data[END_REF]. Many works have proposed specific orderings for color vectors, mainly with the use of lexicographic orderings [START_REF] Aptoula | A comparative study on multivariate mathematical morphology[END_REF], [START_REF] Aptoula | On lexicographical ordering in multivariate mathematical morphology[END_REF], but they use specific assumptions on the ordering of channels. Our proposal is to build on the ideas we presented in [START_REF] Lézoray | Complete lattice learning for multivariate mathematical morphology[END_REF]. One way to define an ordering relation between the vectors of a set T is to use the framework of h-orderings [START_REF] Goutsias | Morphological operators for image sequences[END_REF]. This corresponds to defining a surjective transform h from T to L where L is a complete lattice equipped with the conditional total ordering [START_REF] Goutsias | Morphological operators for image sequences[END_REF]. We refer to ≤ h as the h-ordering given by:

h : T → L and v → h(v), ∀(v i , v j ) ∈ T × T v i ≤ h v j ⇔ h(v i ) ≤ h(v j ) . (1)
Then, T is no longer required to be a complete lattice, since the ordering of T can be induced upon L by means of h [START_REF] Aptoula | Multivariate mathematical morphology applied to colour image analysis[END_REF]. When h is bijective, this corresponds to defining a space filling curve [START_REF] Chanussot | Bit mixing paradigm for multivalued morphological filters[END_REF] or equivalently a rank [START_REF] Lezoray | Rank transformation and manifold learning for multivariate mathematical morphology[END_REF]. It is obvious that the projection h cannot be linear [START_REF] Peyré | Manifold models for signals and images[END_REF] since a distortion of the space topology is inevitable. As a consequence, we choose to focus our developments on learning the vectors' manifold to construct h and deduce the complete lattice (T , ≤ h ). In addition, the same graph signal can be projected on many different graphs (see first line of Figure 1). This shows that the graph itself is not important for the definition of the complete lattice: the manifold of the vertices' vectors is much more important.

C. Manifold-based ordering

Our approach consists in learning the manifold of vectors from a graph signal with a nonlinear mapping and to define the ordering from this projection. To learn the manifold, we use Laplacian EigenMaps (LE), a technique for non-linear dimensionality reduction [START_REF] Lee | Nonlinear Dimensionality Reduction[END_REF]. Computationally, performing LE on the whole space of vectors of the graph signal is not tractable in reasonable time (especially for large sets), so we use a strategy that enables us to construct efficiently a hordering. We summarize its principle in the sequel.

Given a graph signal that provides a set

T of m vectors in R 3 , a dictionary D = {x 1 , • • • , x p } of p
m vectors is built by Vector Quantization [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF].

Manifold learning by Laplacian EigenMaps is performed on this dictionary. One starts by computing a similarity matrix K D that contains the pairwise similarities

K D (i, j) = exp - x i -x j 2 2 σ 2 (2)
between all the dictionary vectors x i . To have a parameterfree algorithm, we consider σ = max

(x i ,x j )∈D x i -x j 2 2 . The normalized Laplacian matrix L = I -D -1 2 D K D D -1 2 D
is then computed with D D the degree diagonal matrix of K D . Then, Laplacian Eigenmaps manifold learning consists in searching for a new representation Φ obtained by minimizing

1 2 ij Φ(x i ) -Φ(x j ) 2 K D (i, j) = T r(Φ T LΦ) (3)
under the constraint Φ T D D Φ = I. This cost function encourages nearby sample vectors to be mapped to nearby outputs. The solution is obtained [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF] by finding the eigenvectors Φ D of L. Therefore, the decomposition

L = Φ D Π D Φ T D is com- puted with corresponding eigenvectors Φ D = [Φ 1 D , • • • , Φ p D ] and eigenvalues Π D = diag[λ 1 , • • • , λ p ].
This obtained projection operator corresponds to constructing a h D -ordering from the data of the dictionary D and a new representation h D (x i ) is obtained for each element x i of the dictionary:

h D : x i → (φ 1 D (x i ), • • • , φ p D (x i )) T ∈ R p . ( 4
)
where φ k D (x i ) denotes the i-th coordinate of the k-th eigenvector. Such a strategy of modeling the manifold from a patch dictionary was also explored in [START_REF] Peyré | Manifold models for signals and images[END_REF]. This correspond to the construction of the complete lattice (D, ≤ h D ) with a h Dordering, and this ordering is only valid for the set of vectors of the dictionary. Since we need the complete lattice (T , ≤ h ), the reduced dictionary lattice is extended to all the vectors of the initial lattice T by Nyström extrapolation [START_REF] Talwalkar | Large-scale SVD and manifold learning[END_REF] of h D on T . To do so, we compute the similarity matrix K DT between sets D and T and the associated degree diagonal matrix D DT . The extrapolated eigenvectors are then obtained by

Φ = D -1 2 DT K T DT D -1 2 D Φ D (diag[1] -Π D ) -1 . (5)
Finally, the projection operator h :

T ⊂ R 3 → L ⊂ R p on the manifold is defined as h(x) = ( φ1 (x), • • • , φp (x))
T , and the complete lattice (T , ≤ h ) is obtained by using the conditional ordering on this new representation. Second line of Figure 1 shows for a given graph signal the learned dictionary, the learned manifold from the dictionary and the final manifold by Nyström extrapolation.

III. GRAPH SIGNAL MORPHOLOGICAL PROCESSING

A. Graph signal representation

Once the complete lattice (T , ≤ h ) is available, a new graph signal representation can be defined. Let P be a sorted permutation of the elements of T according to the manifoldbased ordering ≤ h , one has:

P = {v 1 , • • • , v m } with v i ≤ h v i+1 , ∀i ∈ [1, (m -1)]. (6)
This can also be written as P = PT with P a permutation matrix of size m × m. From this ordered set of vectors, an index graph signal can be defined. Let I : G → [1, m] denote this index graph signal. Its elements are defined as:

I(v i ) = {k | v k = f (v i ) = v i } . (7)
Therefore, at each vertex v i of the index graph signal I, one obtains the rank of the original vector f (v i ) in P, the set of sorted vectors. Given (I, P), a new representation of the original graph signal f is obtained. When a graph signal is encoded in this way, the information is not directly carried by I, but is stored in a separate piece of data called a palette: the set P of sorted vectors. The original graph signal f can be directly recovered since f 1 shows examples of obtained graph signal representations for the same graph signal on two different graphs.

(v i ) = P[I(v i )] = T [i] = v i . Third line of Figure

B. Morphological processing of graph signals

Now that a new representation has been proposed to represent graph signals, we present how morphological processing tasks can be performed with the latter. The erosion and dilation of a graph signal f at vertex v i ∈ G by a structuring element B k ⊂ G are defined as: A structuring element B k (v i ) of size k defined at a vertex v i corresponds to the set of vertices that can be reached from v i in k walks:

B k (f )(v i ) = {P[∧I(v j )], v j ∈ B k (v i )} (8)
δ B k (f )(v i ) = {P[∨I(v j )], v j ∈ B k (v i )} . ( 9
) 1 : G 1 → R n f 2 : G 2 → R n D h D : D ⊂ R n → R p h : T ⊂ R n → R p I 1 : G 1 → [1, m 1 ] P 1 I 2 : G 2 → [1, m 2 ] P 2
B k (v i ) = {v j ∼ v i } ∪ {v i } if k = 1 B k-1 (v i ) ∪ ∪ ∀v l ∈B k-1 (vi) B 1 (v l ) if k ≥ 2 (10
) As detailed in [START_REF] Heijmans | Graph morphology[END_REF], the number of vertices in a given k-hop neighborhood B k (v i ) is highly dependent on the vertex v i , but the associated erosion and dilation are symmetry preserving operators. The formulation of the proposed morphological operators shows that they operate on the index graph signal I, and the processed graph signal is reconstructed through the sorted vectors P that represent the learned complete lattice. It is easy to see that these operators inherit the standard algebraic properties of morphological operators [START_REF] Velasco-Forero | Mathematical morphology for vector images using statistical depth[END_REF]. From these basic operators, we can obtain all standard morphological filters for graph signals such a as openings

γ B k (f ) = δ B k ( B k (f )) and closings φ B k (f ) = B k (δ B k (f )).

IV. GRAPH SIGNAL EDITING APPLICATIONS

To illustrate the benefit of our approach, we provide several examples of its use for editing applications of graph signals. We consider two different types of graphs signals: color vectors assigned to 8-connected grid graphs or triangulated meshes [START_REF] Lézoray | Image Processing and Analysis with Graphs: Theory and Practice, Digital Imaging and Computer Vision[END_REF]. The first type of graph signal corresponds to the classical case of 2D color images, whereas the second case corresponds to colored 3D meshes. This last type of graph signal has recently emerged as a new imaging modality with the advent of low-cost 3D color scanners that enable to simultaneously acquire both the 3D coordinates and color of an object. Very few works have considered editing such a type of graph signal [START_REF] Hidane | Graph signal decomposition for multi-scale detail manipulation[END_REF].

We start by providing graph signal editing for the case of colored 8-connected grid graphs. Figure 2 presents an example of morphological image filtering with an opening by reconstruction γ rec B k with a k = 2-hop as structuring element. Such a processing performs an iterative geodesic dilation starting from an erosion and enables to reconstruct the contours of the objects that have not been totally removed by the erosion [START_REF] Wilkinson | Connected filtering by reconstruction: Basis and new advances[END_REF]. As it can be seen, our method produces excellent smoothing results with flat regions produced when similar colors occur while preserving the important object edges. In addition, no blurring is present on the object boundaries. This last property is essential for structure-preserving filtering and is a consequence of the vector-preserving property of morphological filter. This shows the interest of morphological operators for editing applications, that we investigate in the sequel.

Original image

Filtering result Such smoothed graph signals can be used for nonphotorealistic abstraction with simultaneous detail flattening and edge emphasizing. To do so, the graph signal is strongly smoothed with an Open Close Close Open (OCCO) filter defined as pixel wise average of open-close and close-open [START_REF] Aptoula | A comparative study on multivariate mathematical morphology[END_REF]:

OCCO B k (f ) γ B k (φ B k (f )) + φ B k (γ B k (f )) 2 (11)
From the smoothed image, edges are extracted, enhanced and composed back in the smoothed graph signal to augment the visual distinctiveness of different regions. Figure 3 present two examples of obtained results. Our proposed morphological filtering effectively removes texture from structure and the extracted structure can be used to obtain a non-photorealistic rendering of the scene. Figure 4 presents a first approach towards detail manipulation. A morphological contrast mapping is iteratively applied two times within a 2-hop. This enhances the local contrast of the graph signal f by sharpening its edges with the following transformation, similarly to a shock filter:

Original image Abstraction result

κ B k (f )(v i ) = δ B k (f )(v i ) if ∆ 1 B k (f )(v i ) ≤ ∆ 2 B k (f )(v i ) B k (f )(v i ) if ∆ 1 B k (f )(v i ) > ∆ 2 B k (f )(v i ) (12) with
∆ 1 B k (f )(v i ) = f (v i ) -δ B (f )(v i ) 2 ∆ 2 B k (f )(v i ) = f (v i ) -B (f )(v i ) 2 (13)
The obtained processed images can also be regarded as another way of getting image abstraction results by enhancing the local contrast with a sharpening of edges.

Usually, for detail manipulation, edge-aware transforms are considered [START_REF] Gastal | Domain transform for edge-aware image and video processing[END_REF]. We adopt the the strategy of [START_REF] Farbman | Edge-preserving decompositions for multi-scale tone and detail manipulation[END_REF] that consists in decomposing an image into a base layer and several detail layers. We propose the following multi-scale morphological decomposition of a graph signal:

     f -1 = f , f i = OCCO B k-i (f i-1 ), i ≥ 0, d i = d i-1 -f i , i ≥ 0, (14)
Original image Contrast mapping result where k gives the number of scales and B k-i is a sequence of structuring elements of decreasing sizes with i ∈ [0, k -1]. Indeed, it is clear that in order to extract the successive layers in a coherent manner, the sequence of scales should be decreasing and therefore the size of the structuring element decreases, expressed by B k-i . In terms of graph signal decomposition, this assumption has the following simple interpretation: as the process evolves, the successive steps extract more details from the original graph signal (similarly as [START_REF] Hidane | Graph signal decomposition for multi-scale detail manipulation[END_REF]). The layer f 0 should be interpreted as a first sketch of the f , while the residuals d i are to be understood as detail layers. The graph signal can then be represented by:

(∀k ≥ 0) f = k-1 i=0 f i + d k-1 . (15)
The f i 's thus represent different layers of f captured at different scales. In Figure 5, we apply three levels of decomposition (k = 3). The image is reconstructed from the obtained decompositions with specific coefficients for each layer. The base layer is kept unchanged but the two detail layers are boosted with factors 2 and 3. One can see that the image is recomposed without the presence of halos since the decomposition is a vector-preserving one. This shows that our framework can also be advantageously used for such edgeaware image detail manipulation, and this has never been explored before with morphological methods.

Up to now, we have presented results for graph signals that correspond to colors defined on grid-graphs, the classical case of images. Now we show results for much challenging data: 3D colored meshes. Figure 6 presents results for detail manipulation on two different meshes. The first (the horse) was obtained by using a 3D color scanner and the second (the fire hydrant) by photogrammetry. To manipulate the details of these mesh and modify the colors of the vertices, we consider this time a toggle contrast mapping filter:

τ B k (f )(v i ) =      δ B k (f )(v i ) if ∆ B k (f )(v i ) > 0 B k (f )(v i ) if ∆ B k (f )(v i ) < 0 f (v i ) otherwise (16)
with the morphological Laplacian defined as

∆ B k (f )(v i ) = δ B k (f )(v i ) -2f (v i ) + B k (f )(v i ) . ( 17
)
We applied 5 iterations of such a contrast enhancing filter on a 2-hop. As it can be seen, the colored mesh has been simplified and the contrast around edges enhanced. This provides a way of manipulating such graph signals to enhance and obtain an abstracted result of a raw 3D scanning.

Original Toggle contrast mapping Finally, as for images, detail manipulation can be performed by iterative signal decomposition into a base and several detail layers. We consider, as for images, a morphological OCCO filter as a basis for the decomposition and decompose the color mesh graphs signal on 4 levels. The signal is then recomposed with specific coefficients for each layer. Figure 6 presents an exemple of such detail manipulation. For colored 3D meshes, this enables to manipulate the colors of the mesh vertices to sharpen it. Finally, it is important note that, it is the first time that a morphological method is proposed to sharpen 3D colored meshes and this shows the innovation of our proposal. The potential of this method is huge since many 3D scanners now also acquire a color per vertex, but this information is never used to manipulate the details of the point cloud. V. CONCLUSION This paper has detailed an approach for the morphological editing of graphs signals. To address this, a new graph signal representation has been proposed in the form of an ordering of vectors and a graph signal index. The ordering of vectors is interpreted as the construction of a complete lattice and a manifold-based approach is proposed towards this. A strategy combining dictionary learning, manifold learning and out of sample extension has been devised to learn the manifold used to construct the complete lattice of color vectors. From the new representation, morphological operators that operate on graphs are formulated and their use for graphs signal editing tasks has been studied. In particular, morphological multiscale decomposition of graph signals have been proposed for detail manipulation. In future works, we will investigate which morphological filter is best suited to produce the different layers of the proposed morphological multi-scale graph signal decomposition. Indeed, many potential morphological filters can be considered (opening by reconstruction, OCCO, leveling, alternate sequential filter, etc.) and this needs to be investigated (similarly to what has been done in [START_REF] Karantzalos | Comparing morphological levelings constrained by different markers[END_REF], [START_REF] Maragos | Leveling cartoons, texture energy markers, and image decomposition[END_REF] for images).

Fig. 1 .

 1 Fig. 1. First line: The same graph signal on two different graphs (planar grid graph and 3D mesh) with m 1 = 262144, m 2 = 10000. Second line: the dictionary D with p = 64 colors, the mapping h D (projection shown on the three first eigenvectors) learned from the dictionary D (dilated 5 times for visualization purposes), the mapping h (projection shown on the three first eigenvectors) interpolated to all the vectors of f 1 . Third line: the index graph signals I 1 and I 2 and the associated sorted vectors P 1 and P 2 for the two graph signals of the first line.

Fig. 2 .

 2 Fig. 2. Morphological image filtering (opening by reconstruction on a 2-hop), see text for details.

Fig. 3 .

 3 Fig. 3. Morphological image abstraction (OCCO on a 2-hop), see text for details.

Fig. 4 .

 4 Fig. 4. Morphological image contrast mapping (2 iterations of a contrast mapping on a 2-hop), see text for details.

5 .

 5 Morphological detail manipulation (3 levels of decomposition with a OCCO filter of decreasing size with an initial 3-hop size -reconstruction coefficients of each layer are 1, 2 and 3), see text for details.

Fig. 6 .

 6 Fig. 6. Morphological colored mesh toggle contrast mapping (5 iterations of a toggle contrast mapping on a 2-hop), see text for details.

Fig. 7 .

 7 Fig. 7.Morphological colored mesh detail manipulation (4 levels of decomposition with a OCCO filter of decreasing size with an initial 4-hop size -reconstruction coefficients of each layer are 1, 1.25, 1.5 and 1.75), see text for details.

ACKNOWLEDGMENT

This work received funding from the Agence Nationale de la Recherche, ANR-14-CE27-0001 GRAPHSIP.


Franck Hétroy-Wheeler

Eric Casella
email: eric.casella@forestry.gsi.gov.uk

Dobrina Boltcheva

Segmentation of tree seedling point clouds into elementary units

Keywords: terrestrial laser scanning, tLiDAR, 3D point cloud, segmentation, spectral clustering, tree seedling

This paper describes a new semi-automatic method to cluster TLS data into meaningful sets of points to extract plant components. The approach is designed for small plants with distinguishable branches and leaves, such as tree seedlings. It first creates a graph by connecting each point to its most relevant neighbours, then embeds the graph into a spectral space, and finally segments the embedding into clusters of points. The process can then be iterated on each cluster separately. The main idea underlying the approach is that the spectral embedding of the graph aligns the points along the shape's principal directions. A quantitative evaluation of the segmentation accuracy, as well as of leaf area estimates, is provided on a poplar seedling mock-up. It shows that the segmentation is robust with false positive and false negative rates around 1%. Qualitative results on four contrasting plant species with three different scan resolution levels each are also shown.

Introduction

Functional-structural plant models describe a plant as a collection of interconnected elementary units (internode, petiole, leaf-blade, see [START_REF] Godin | A method for describing plant architecture which integrates topology and geometry[END_REF]). Their goal is to help biologists understand the relationships between the plant structure and the biological and physical mechanisms underlying the plant growth [START_REF] Godin | Functionalstructural plant modelling[END_REF]). These models require an in situ validation on real plants, that can be done by measuring the three-dimensional (3D) characteristics of vegetation. Similarly the growing field of plant phenomics, which is concerned with the discovery and analysis of complex plant traits [START_REF] Furbank | Phenomics -technologies to relieve the phenotyping bottleneck[END_REF]; International Plant Phenotyping Network (2016)), requires the measurement of individual quantitative parameters such as leaf characteristics.

Destructive measurements have long been used but are time consuming and expensive. As a consequence, various kinds of sensors are being investigated for non-destructive and non-invasive plant metrology. For example, the use of different imaging techniques has been proposed for plant phenotyping, see [START_REF] Li | A review of imaging techniques for plant phenotyping[END_REF] for a review. The most popular imaging techniques are based on single-lens cameras [START_REF] Quan | Image-based plant modeling[END_REF]; [START_REF] Paproki | A novel mesh processing based technique for 3d plant analysis[END_REF]), time-of-flight cameras [START_REF] Chéné | On the use of depth camera for 3d phenotyping of entire plants[END_REF]; [START_REF] Alenyà | Robotized plant probing: Leaf segmentation utilizing time-of-flight data[END_REF]; [START_REF] Chaivivatrakul | Automatic morphological trait characterization for corn plants via 3d holographic reconstruction[END_REF]; [START_REF] Xia | In situ 3d segmentation of individual plant leaves using a rgb-d camera for agricultural automation[END_REF]) or multi-view stereo imaging systems [START_REF] Golbach | Validation of plant part measurements using a 3d reconstruction method suitable for highthroughput seedling phenotyping[END_REF]; [START_REF] Lou | Estimation of branch angle from 3d point cloud of plants[END_REF]; [START_REF] Rose | Accuracy analysis of a multi-view stereo approach for phenotyping of tomato plants at the organ level[END_REF]). All these methods allow one to reconstruct and measure single leaves, although some of them require manual interaction [START_REF] Quan | Image-based plant modeling[END_REF]; [START_REF] Golbach | Validation of plant part measurements using a 3d reconstruction method suitable for highthroughput seedling phenotyping[END_REF]; [START_REF] Rose | Accuracy analysis of a multi-view stereo approach for phenotyping of tomato plants at the organ level[END_REF]) or prior knowledge of the plant [START_REF] Chéné | On the use of depth camera for 3d phenotyping of entire plants[END_REF]; [START_REF] Chaivivatrakul | Automatic morphological trait characterization for corn plants via 3d holographic reconstruction[END_REF]).

Terrestrial laser scanning (TLS), a remote sensing technique, has become an increasingly popular technique to measure vegetation from grass to forest plant species, see e.g. [START_REF] Dassot | The use of terrestrial lidar technology in forest science: application fields, benefits and challenges[END_REF]; [START_REF] Lin | Lidar: An important tool for next-generation phenotyping technology of high potential for plant phenomics?[END_REF] for recent reviews. Compared to imaging techniques, TLS provides direct accurate 3D measurements. It has also proved to be more robust to diverse environments, in particular to changing lighting conditions [START_REF] Li | A review of imaging techniques for plant phenotyping[END_REF]; [START_REF] Lin | Lidar: An important tool for next-generation phenotyping technology of high potential for plant phenomics?[END_REF]). Thus, LiDAR seems more adapted to greenhouse and field conditions [START_REF] Tilly | Evaluation of terrestrial laser scanning for rice growth monitoring[END_REF]; [START_REF] Lin | Lidar: An important tool for next-generation phenotyping technology of high potential for plant phenomics?[END_REF]).

TLS generates unstructured sets of points where its laser beam is incident and reflected. Thus it gives a raw sketch of the spatial distribution of plant elements in 3D, but it lacks explicit and essential information on their shape and connectivity. The points need to be clustered into geometrically meaningful sets for further analysis and dendrometric measurements. For example, leaf-blade points need to be separated from petiole and internode points to assess leaf areas.

In this paper, we segment TLS data of small plants or tree seedling scans into their elementary units: internodes, petioles and leaf-blades. Our method only considers the 3D positions of the points (no intensity value or normal estimate is required). As a consequence, it can be applied to sets of 3D points generated by other techniques than TLS such as time-of-flight cameras. We focus on accurate segmentation so that individual elementary unit characteristics such as leaf area are estimated as accurately as possible.

Related work

Segmentation of 3D data is critical for many applications in science. Research has considered the segmentation of point clouds into basic geometric primitives (planes, cylinders, spheres, etc.) for various purposes, such as building or city modelling (see [START_REF] Haala | An update on automatic 3d building reconstruction[END_REF] for a survey), reverse engineering of mechanical objects (e.g. [START_REF] Bey | Reconstruction of consistent 3d cad models from point cloud data using a priori cad model[END_REF]; [START_REF] Li | Globfit: consistently fitting primitives by discovering global relations[END_REF]), or background subtraction, see [START_REF] Nguyen | 3d point cloud segmentation: a survey[END_REF] for a recent overview. These approaches are designed for man-made objects which can be almost completely decomposed into uniform geometric shapes. Yet since stems and leaves are not exactly cylindrical and planar shapes, their efficiency to robustly segment plants is questionable.

The problem we are interested in, that is the segmentation of plants into their elementary units, has been partially tackled in the literature. Recovering the branching structure of a plant or a leaf-on tree is a specific issue that has been addressed in order to estimate various wood parameters by e.g., [START_REF] Bayer | Structural crown properties of norway spruce (Picea abies [l.] karst.) and european beech (Fagus sylvatica [l.]) in mixed versus pure stands revealed by terrestrial laser scanning[END_REF]; [START_REF] Belton | Processing tree point clouds using gaussian mixture models[END_REF]. Some authors additionally propose to reconstruct the foliage, using some heuristics to position them, for example to create visually pleasing virtual 3D models of trees from TLS point clouds (e.g., [START_REF] Xu | Knowledge and heuristic based modeling of laser-scanned trees[END_REF]; [START_REF] Livny | Automatic reconstruction of tree skeletal structures from point clouds[END_REF]) or to derive global characteristics such as total leaf and wood areas [START_REF] Côté | The structural and radiative consistency of three-dimensional tree reconstructions from terrestrial lidar[END_REF]). The number of leaves as well as their individual location and shape are plausible but do not correspond to the actual tree.

Other works segment a plant into two clusters only, one for the stems and one for the leaves, mostly by classifying points according to local geometric features [START_REF] Belton | Processing tree point clouds using gaussian mixture models[END_REF]; [START_REF] Paulus | Surface feature based classification of plant organs from 3d laserscanned point clouds for plant phenotyping[END_REF]Paulus et al. ( , 2014)); [START_REF] Wahabzada | Automated interpretation of 3d laserscanned point clouds for plant organ segmentation[END_REF]). Other approaches use geometric distance information [START_REF] Tao | A geometric method for wood-leaf separation using terrestrial and simulated lidar data[END_REF]) or intensity information [START_REF] Douglas | Finding leaves in the forest: the dual-wavelength echidna lidar[END_REF]). Deriving a full segmentation of the plant from such a classification is possible in some cases, using prior knowledge about the plant or its organs (Paulus et al. (2014)). However this is not straightforward in general, especially when leaves are almost overlapping.

Geometric segmentation of a plant into elementary units has been proposed using either plant-specific prior knowledge (e.g., [START_REF] Kaminuma | Automatic quantification of morphological traits via three-dimensional measurement of arabidopsis[END_REF]) or a tedious interactive procedure [START_REF] Dornbusch | A method to extract morphological traits of plant organs from 3d point clouds as a database for an architectural plant model[END_REF]; [START_REF] Hosoi | 3-d modeling of tomato canopies using a high-resolution portable scanning lidar for extracting structural information[END_REF]Paulus et al. (2014)). In contrast, our approach only requires a minimum user interaction, and no prior knowledge. It is thus applicable to any species. [START_REF] Yin | Full 3d plant reconstruction via intrusive acquisition[END_REF] have recently proposed a destructive approach to accurately segment and reconstruct a pot plant. Their approach requires laser scanning the whole plant, manually cutting the plant leaves and laser scanning each leaf individually. This is time-consuming and obviously does not allow for tracking changes in the plant traits over time.

A approach analogous to that used here has recently been published by [START_REF] Lou | Estimation of branch angle from 3d point cloud of plants[END_REF], although they work on point cloud data generated using a multi-view stereo imaging system. The methodology is similar: a graph is first constructed to build neighbourhood relationships between the plant's points, then a spectral clustering approach is performed on this graph. However, we used more advanced graph construction and clustering techniques, as detailed in Sections 2.1.2 and 2.1.4.

Approach

The main contributions of this paper are:

(1) a semi-automatic approach to accurately segment a TLS 3D point cloud of a tree seedling into meaningful clusters of points (Section 2.1). More specifically, each cluster gathers points of an elementary unit of the tree seedling.The approach is global in the sense that leaves are not segmented first from branches, and robust to non uniform density in the point cloud;

(2) the assessment of the method robustness by quantitatively evaluating clustering results and individual leaf-blade areas on computer-generated scans from a plant mock-up (Sections 3.2.1 and 3.2.2).

Qualitative validation on real plant scans and parameter sensitivity analysis are also given (Sections 3.1 and 3.2.3, respectively).

Materials and methods

The segmentation algorithm underlying the approach is first explained in detail in Section 2.1. Reference data used in the experiments is then described in Section 2.2. The method used for statistical analysis of leaf area estimates is explained in Section 2.3.

Point cloud segmentation algorithm

In this section, the algorithm that analyses data collected from the TLS is described. The input point cloud data is merely a 3D location of points with no additional information. The algorithm is designed to cluster points into subsets corresponding to the plant elementary units. This algorithm is a three-stage process (Figure 1). The first stage, called Graph Construction, finds neighbouring points for each point of the raw TLS data (see Section 2.1.2 for details). The second stage, called Spectral Embedding, finds the major intrinsic plant directions, i.e. the main directions of each elementary unit (see Section 2.1.3). This allows us to define the distances between neighbouring points according to the intrinsic plant directions, rather than the usual Euclidean distance. For example, the distance between two points sampled on a leaf-blade with an ellipsoid shape corresponds to the distance between their projections on the leaf's midrib (ellipsoid's main axis). Thus, two points on both sides of the ellipsoid's main axis but with similar projections will appear close to each other. As a consequence, this stage of the algorithm transforms the raw TLS data into a cloud of points aligned along principal plant axes (see Figure 2). Finally, the third stage, called Spectral Clustering, uses the computed neighbouring relationships to decompose the shape into subsets of points according to the principal plant axes. All points in a subset are given the same label (in our experiments, a colour), and points in different subsets have different labels. During this stage, each elementary unit is thus split from the one it originates. For example, a leaf-blade is separated from its petiole (see Section 2.1.4). Since each point in the embedding space corresponds to a point in the Euclidean space, the segmentation of the input TLS data is automatically found by giving to each point the label of its associated point in the embedding space (reverse embedding).

Parameters

The algorithm uses three parameters, one for each stage:

(1) the minimum angle a between two neighbours of any point in the point cloud, for the graph construction; (2) an estimate d of the number of intrinsic directions in the plant, for the spectral embedding;

(3) the number c of desired subsets of points (elementary units), for the clustering stage.

Graph construction

The first stage of the method aims at recovering neighbouring information between points. This is a difficult task since the only information available is the 3-D location of the points.

Usual methods create neighbouring relationships, called edges, between any point p and either all points which lie within a sphere of radius ε centred at p, or the k nearest points (Figure 3 (a,b)). These methods are known as the ε-Neighbourhood and the k-Nearest-Neighbours methods, respectively [START_REF] Yang | Building connected neighborhood graphs for isometric data embedding[END_REF]; von Luxburg (2007)). ε and k are user-chosen parameters. ε-Neighbourhood is for example used by [START_REF] Belton | Processing tree point clouds using gaussian mixture models[END_REF], while the k-Nearest-Neighbours method is used by [START_REF] Côté | The structural and radiative consistency of three-dimensional tree reconstructions from terrestrial lidar[END_REF]; [START_REF] Lou | Estimation of branch angle from 3d point cloud of plants[END_REF]. These methods are convenient so long as the density of the point cloud is uniform, which is not the case for our TLS data. For non-uniform samples, many redundant edges may be created or relevant ones may be missed and the main problem is to find the right value for the parameters. This problem is shown on Figure 3 (a), where the ε-neighbourhood of a blue point is depicted for two different values of ε (in green and in red and green, respectively). Similarly, the k-nearest-neighbours, for k = 2 (in green) and k = 5 (in red and green), are shown on Figure 3 Note that it is critical for the next stage of our approach, both in terms of memory usage and computation time, to avoid redundant neighbouring relationships since during this stage we work with an adjacency matrix computed from the neighbouring graph. Therefore, the lower the number of neighbours for a given point the sparser the matrix thus the faster the computation. This is why we have developed a specific algorithm which is summarised in Algorithm 1 and Figure 3 (c). This algorithm starts by selecting a number k of candidate neighbours for every point p ∈ P . In practice, we choose k = 0.1% of the total number of points. Then, in order to select the neighbours within the set of candidates, it uses one parameter which is the minimum angle a at p allowed between two edges with endpoint p (see Figure 3 (c)). If many candidates lie in the same direction, only one (the closest) is selected as a neighbour of p. This prevents the creation of redundant, almost parallel, edges. On the contrary, this algorithm having found the closest point in a given direction will go on to look for other points farther away but in a distinctly different direction. Thus, this method captures all relevant edges and is robust to non-uniform density within the point cloud. Figure 3 (c) shows the result of the method for a = 45 • . This method discards the two red points D since the corresponding edges are within a small angle of existing edges (in green) but it does capture point 3 which is a neighbour in a clearly different direction.

The graph construction runs this algorithm for every point p in the point cloud. The resulting neighbouring graph is thus the union of the selected edges E(p) for every p ∈ P . Note that we do not look for mutual nearest neighbours: if q is computed as a valid neighbour of p but p is not considered as a valid neighbour of q, we still connect these two points. Our experimental results have shown that choosing the angle parameter a = 90 • is a good compromise in practice (see Section 3.2.3). This allows us to search for neighbours in the 3 cardinal directions around a point in 3D, which has been shown to be sufficient for building Figure 4 shows an example of a graph construction. In this example, the petioles were very sparsely scanned compared to the leaf-blades and the main branch. For the ε-neighbourhood and the k-nearest-neighbours methods, the minimum value of the parameter was chosen such that the resulting graph was connected. Notice how both methods, contrary to ours, create numerous redundant edges on the main branch. Other methods, which guarantee connectedness of the graph, are proposed by [START_REF] Yang | Building connected neighborhood graphs for isometric data embedding[END_REF]. However, their computational complexity (at least O(n 2 ), where n is the number of points) may become prohibitive in the context of this study. The approach proposed in Algorithm 1 reaches a O(n log n) complexity with appropriate data structures, i.e. a kd-tree for the k-nearest-neighbours searches and heaps for the priority queues.

Spectral embedding

In the second stage, the major intrinsic directions of the shape are recovered and the weights of the edges modified accordingly. This is done using a technique called dimension reduction or spectral embedding. Indeed, embedding a (discrete) shape into a low-dimensional spectral space is known to help recover its intrinsic features (see e.g. [START_REF] Reuter | Laplace-beltrami spectra as 'shape-dna' of surfaces and solids[END_REF]). In this work, we build on the Laplacian Eigenmaps framework of [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF], the main differences being March 16, 2016 International Journal of Remote Sensing plantscan the graph construction approach described above and the choice of the distance between neighbouring points. This framework is now described.

Let A be the adjacency matrix of the graph constructed in the previous stage. Points are numbered from 1 to n, A is a n × n matrix such that A(i, j) is equal to the weight of the edge connecting points i and j. A(i, j) = 0 if there is no edge between these points. The Euclidean distance between i and j is used as a weight. Let W be the diagonal valency matrix of the graph. W (i, i) is equal to the sum of the weights of edges with endpoint i. The matrix L = W -A is called the Laplacian matrix of the graph. The spectral embedding of the graph into a d-dimensional space is given by the d eigenvectors V 1 , . . . , V d of L associated with the first d non zero eigenvalues (in increasing order). Namely, the embedding coordinates of point number i are given by row i of the matrix whose columns are vectors V 1 , . . . , V d (von Luxburg ( 2007)).

It is known that the eigenvectors associated to the lowest non zero eigenvalues of L give the main "intrinsic" (curved) directions of the graph [START_REF] Lévy | Laplace-beltrami eigenfunctions: towards an algorithm that understands geometry[END_REF]). This property has previously been used for shape compression [START_REF] Karni | Spectral compression of mesh geometry[END_REF]), progressive reconstruction [START_REF] Lévy | Laplace-beltrami eigenfunctions: towards an algorithm that understands geometry[END_REF]) and deformation [START_REF] Dey | Eigen deformation of 3d models[END_REF]) purposes. The Laplacian spectral embedding is also known as the eigen skeleton of the input graph [START_REF] Dey | Eigen deformation of 3d models[END_REF]). Using this property makes sense in the context of this study, since a plant is a strongly anisotropic shape; the natural directions of the plant follow the directions of each stem, branch, petiole and the main directions of each leaf-blade. It is therefore expected that the spectral embedding of the graph aligns points into a curve, or at least a strongly anisotropic shape, that samples each elementary unit of the plant, as shown on Figure 2. It is easier to segment the spectral embedding of the graph into subsets of points than the TLS data, since it does not depend on the particular shape of the leaves. Moreover, geometrical noise accumulated during the acquisition process is implicitly altered by the spectral embedding.

Note that computing the eigen-decomposition depends on the number of edges in the graph. The lower number of neighbours a point has, the sparser the matrix is, thus the faster the computation is. This is why the algorithm described in Section 2.1.2 is used rather than the standard ε-neighbourhood or k-nearestneighbours methods. Figure 2 shows the 3-D embedding of a simple plant model with two leaves, thus having three main directions. Notice how the plant is nearly collapsed to a set of curves.

Spectral clustering

This stage clusters points into sets corresponding to the plant's elementary units (internodes, petioles and leaf-blades). To this aim, the point cloud is segmented according to its spectral embedding; the objective is thus to cluster together points of an elongated curve in the embedded shape. Since the embedded point cloud is almost a set of elongated curves (see Figure 2), segmentation techniques for stems such as e.g. Paulus et al. (2014) could be applied. However, they would not benefit from the point neighbourhood information retrieved from stage one of the approach (Section 2.1.2).

The usual clustering technique, applied in spectral space, is known as K-means clustering [START_REF] Von Luxburg | A tutorial on spectral clustering[END_REF]). It is used for example by [START_REF] Lou | Estimation of branch angle from 3d point cloud of plants[END_REF].

K-means clustering randomly selects K initial "means" among the points, with K being a user-defined parameter. Each point is assigned to the nearest mean. Then, for each cluster of points, the closest point to the centroid (centre of mass) of the cluster is computed and selected as the new mean. The process is iterated until convergence to stable mean positions is reached, which is generally fast. This technique is well adapted to isotropic data, i.e, a point cloud without any principal direction. This is obviously not the case in this study where the graph is embedded in spectral space almost as a set of elongated curves. More general approaches such as expectation maximisation could be used, but as K-means clustering they do not naturally benefit from the neighbourhood information (graph edges). Note that [START_REF] Lou | Estimation of branch angle from 3d point cloud of plants[END_REF] merge neighbouring clusters with similar normals, but this may lead to undersegmentation since different elementary units (e.g., two leaves) may have similar normals.

A new clustering method, more adapted to elongated shapes, is therefore proposed. This method is described in Algorithm 2 and Figure 5. The idea is to compute the main directions of the graph (in spectral space), as sets of edge-connected points which are called the segments. As many segments as the desired number c of clusters are computed. Finally, each point of the graph is labelled according to its closest segment. Note that c should be odd, by construction. In case the desired number of clusters is even, we recommend to segment in c + 1 clusters and merge two of them.

Data: Graph G = (V, E) (in spectral

Algorithmic details

Edges of the computed graph are weighted by a distance between their two endpoints called the commute-time distance, which represents the expected time for a random walk on the graph to travel from one point to the other and then return [START_REF] Qiu | Clustering and embedding using commute times[END_REF]). In our plant segmentation context, this is a more meaningful distance than the Euclidean distance. For example, points on two different leaf-blades connected by a few edges (see Figure 6) may have a short Euclidean distance and a large commute-time distance in the graph. Since we want such points to belong to different clusters, we want their distance to be large. Moreover, commute-time distance has been proved to be robust against noise for clustering purposes [START_REF] Qiu | Clustering and embedding using commute times[END_REF], Sec. 5.1).

The commute-time distance is similar to the Euclidean distance in spectral space, except each coordinate is divided by the corresponding eigenvalue. More precisely, the commute-time distance between points i and j is given by

k (i(k) -j(k)) 2 e(k)
, with i(k) and j(k) the k-th coordinates in spectral space of i and j, respectively (that is to say, the i-th and j-th coordinates of the k-th eigenvector V k of the Laplacian matrix of the graph, as explained above), and e(k) the k-th eigenvalue of the Laplacian matrix of the graph [START_REF] Qiu | Clustering and embedding using commute times[END_REF], Sec. 2.3).

Finding the main segment of a weighted graph is a typical issue in medial structure axis and skeleton-related problems. The main segment of a graph can be computed successfully by using a one-source shortest path algorithm from an endpoint of the graph, e.g. the Dijkstra's algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF]). This endpoint can be found as the farthest point to some random point [START_REF] Lazarus | Level set diagrams of polyhedral objects[END_REF]) and computed once again using a one-source shortest path algorithm. Other seg- ments are then computed the same way by taking all points of already computed segments as source points as in [START_REF] Hassan | Automatic localization and quantification of intracranial aneurysms. Computer Analysis of Images and Patterns[END_REF]. As a result, each point of the graph is linked to its closest point on the segments and its distance to this point is computed. See Figure 5 for an example.

Asymptotic computational complexity

As explained above, Algorithm 2 uses a one-source shortest path algorithm (c -1)/2+3 times. Then, a Disjoint Set data structure [START_REF] Cormen | Introduction to algorithms[END_REF]) is used to cluster and label the points according to their closest point on the segments. The computational complexity of Dijkstra's algorithm, using a heap data structure, is O(m + n log n), where n is the number of points in the graph and m is the number of edges. The complexity of cluster creation within a Disjoint Set framework and using relevant heuristics is O(n log n) [START_REF] Cormen | Introduction to algorithms[END_REF]). The computational complexity of Algorithm 2 is thus O(c(m + n log n)).

Reference data

Reference point clouds were obtained at various resolution levels, from two different ways. First, point clouds were generated from a virtual poplar seedling mock-up through a computer simulation of TLS (Section 2.2.1). Second, points clouds were acquired from four real plants using a Leica Geosystems HDS-6100 TLS device (Section 2.2.2).

Point cloud computations from a poplar mock-up

The 3D structure of a one-year-old single-stem seedling of poplar clone Trichobel (Populus trichocarpa Torr. & Gray x P. trichocarpa) was generated by the 3D Coppice Poplar Canopy Architecture model (3D CPCA) developed by [START_REF] Casella | A method for describing the canopy architecture of coppice poplar with allometric relationships[END_REF] (Figure 7, Table 1). The model is based on a multi-scale decomposition of a plant structure into components (axis and growth unit) described as a collection of metamers, themselves defined as a collection of elementary units (nodes plus internodes, petioles and leaf-blades) [START_REF] Godin | A method for describing plant architecture which integrates topology and geometry[END_REF]). For this study, axes (stem and branches) were divided as a sequence of conical frustums (a sequence of internode units), petioles were represented as cylinders and leaf-blades were regarded as planar objects. Each elementary unit was scaled to the appropriate geometric dimensions (e.g. height, base and top radius for a conical frustum) although a leaf-blade prototype was created and represented as a polygon with a set of 4 contiguous triangles to fit the leaf-blade shape and the allometric relationships between the leaf-blade area, the leaf-midrib length and the leaf-blade width. Each unit was them rotated and translated according to its orientation and location in the scene. Each unit was scaled so that no discontinuity between elements was possible, and there were no contact between laminae. Empirical functions and random deviation used in this study for the reconstruction of the plant architecture were as in [START_REF] Casella | A method for describing the canopy architecture of coppice poplar with allometric relationships[END_REF]. The resulting poplar mock-up consisted of 17 leaves, 17 petioles and 24 internodes. 2). The last image shows the point cloud generated without simulation of the occlusion: all hits from the laser source to objects were recorded.

(L) (M) (H) (U) (U)
Seedling Height (m) Nb. of leaves Total leaf area (10 This mock-up could then be scanned from any point of view, after having placed a virtual TLS in the scene. Point clouds were computed for three TLS positions around the mock-up and for four scanner resolution levels i.e. by simulating the characteristics and settings of a Leica Geosystems HDS-6100 TLS device (Table 2) used in this study for point cloud acquisitions from real plants (see next Section). The positions of the virtual TLS in the scene were computed for a distance of 3 meters from the base of the stem to the laser source, an elevation angle of 25 • and an azimuth angle of 0, 120 or 240 • . For each TLS resolution level, a point cloud was generated using a simple hit/not hit determination algorithm coded from a set of ray/objects (i.e. /cylinder, /cone and /triangle) intersection equations (see [START_REF] Haines | Essential ray tracing algorithms[END_REF]) by determining either the closest or all hits from the laser source along any simulated ray trajectory within the scene. Every recorded hit (x hit , y hit , z hit ) was then referred to both its related object in the scene (e.g. leaf-blade #) and the position of the laser source (x T LS , y T LS , z T LS ). Point clouds were generated both with and without simulation of the occlusion, i.e. recording either only the closest or all hits from the laser source to objects.

Point cloud acquisitions from real tree seedlings

TLS point clouds were acquired from tree seedlings of silver birch (Betula pendula Roth), horse chestnut (Aesculus hippocastanum L.), sweet chestnut (Castanea sativa Mill.) and red oak (Quercus rubra L.) (Table 1) using a Leica Geosystems HDS-6100 TLS device (Table 2). Each seedling was scanned inside a large area warehouse from three TLS positions around the plant (two for the Horse chestnut seedling) and for three scanner resolution levels. These species were chosen in order to get seedlings with varying structural and leaf geometrical complexities.

TLS point clouds of trees are usually noisy because of various interferences during the acquisition process, see e.g. [START_REF] Hebert | 3-d measurements from imaging laser radars: how good are they?[END_REF]. Each point cloud was thus filtered in order to remove most of the outliers, using the statistical outlier removal filter of the Point Cloud Library [START_REF] Rusu | 3d is here: Point cloud library (pcl)[END_REF]). For each point p, its k nearest neighbours were first retrieved, and the mean distance d of these points to p was computed. If d is outside an interval defined by the mean and the standard deviation of the mean distances to all points, then p is detected as an outlier and removed. We have set k so that to remove as many outliers as possible, without removing relevant points. k = 20 has been taken for the horse chestnut, k = 30 for the sweet chestnut, and k = 50 for the silver birch and the red oak, for all resolution levels. Between 0.2% (red oak) and 11% (silver birch) of the points have been removed by this filtering. Table 3 shows the number of point in each point cloud after filtering. Each resolution level contains about 4 times more points than the previous one.

Statistical analysis of leaf area estimates on the poplar mock-up

Leaf area (LA) has been estimated for each labelled leaf-blade, by projecting its points into the least-square fitting plane, computing the Delaunay triangulation of the projected points [START_REF] Edelsbrunner | Geometry and topology for mesh generation[END_REF]), projecting the points back to their original positions and summing the areas of the Delaunay triangles. This was done for the leaf clusters as labelled in the input data, both without and with occlusions, as well as for the clusters computed with the presented algorithm. The quality of the method was then determined by two parameters, the root mean square error (RM SE) and bias (b), defined as:

RM SE = n k=1 (ŷ pk -y ak ) 2 n b = n k=1 (y sk -y ak ) n
where n is the number of observations and ŷpk is the predicted average value from the regression line between the simulated y sk and the actual y ak values for the k th observation.

Results and discussion

The algorithm has been implemented in C++ and Matlab. Because of the generally complex structure of a plant, perfect clusters may not be created in a single run. In practice, the algorithm is first run with a low number of desired clusters (less than the actual number of metamers), then each cluster is segmented by running this algorithm again. A simple graphical user interface has also been implemented, which allows merging clusters by selecting a point in each cluster. The overall approach is thus semi-automatic.

Results of the segmentation process on the poplar mock-up and on the four tree seedlings are shown in Section 3.1. A quantitative validation is provided in Section 3.2. It includes an evaluation of the segmentation accuracy, a statistical analysis at the leaf area scale, and a parameter sensitivity analysis of the algorithm.

Qualitative results

The method has been tested on the five different tree seedlings (Tables 1 and3), for the low, medium and high resolution levels (Table 2), as well as the ultra high resolution level for the poplar mock-up.

First segmentation

Results of the first run of the algorithm are shown on Figure 8 for the poplar mock-up (ultra high resolution), the sweet chestnut (high resolution) and the red oak (medium resolution) seedlings. When a small number c of clusters is set, the algorithm usually segments the point cloud into connected subsets of elementary units, even when the point cloud is very noisy (e.g., the red oak). The higher value for c, the higher probability that a elementary unit (usually, a leaf) is segmented by the algorithm into several clusters (see Figure 9). We elaborate on the choice of c in Section 3. Once the initial point cloud has been segmented, the user can select any given cluster through the graphical interface and re-run the algorithm on this cluster. This is done interactively; no botanical knowledge is used in our approach and the user decides which subsets of points to segment and when to stop this process.

Final segmentations

Qualitative final segmentation results on all scans are shown on Figures 10,11, 12 and 13. The accompanying videos also show the segmentation results for the high resolution level point cloud of each of these five seedlings.

These results show that overall, despite large occlusions in real scans (see e.g. Figures 12 (d) and 13 (b)), the method correctly segments the point cloud into sets of individual leaf-blades, petioles and stem sections. Internodes can be detected when both ends are delimited by petioles and/or incident stems, otherwise they are merged. The method is insensitive to the leaf anatomy. It behaves correctly for both simple, small (e.g. sweet chestnut) and complex, large (e.g. red oak) leaves, as well as for both planar and curved leaves. However, a compound leaf is segmented into its leaflets, as shown for the horse chestnut, as each leaflet corresponds to a different intrinsic direction. Figures 10 and 11 show that the resolution level does not have a strong influence on the segmentation, as will be demonstrated in Section 3.2. The algorithm is also robust to non uniform density within a point cloud, as shown for example on Figure 12 (a). Finally, the method is insensitive to the noise level. Even when points are spread over the boundaries of a unit (leaf or stem), they are included into the correct cluster (see Figures 12 (a) and 13 (c)). This is also shown by the following experiment.

Robustness to acquisition noise

In order to test the robustness of the approach, a raw scan of the sweet chestnut (high resolution level) from a single viewpoint has also been segmented. 72754 points belong to this unfiltered point cloud. Results are shown on Figure 14, to be compared with Figures 11 (H) and 12 (c). Points are correctly assigned to their corresponding cluster, except on ambiguous areas (for example between two neighbouring leaves).

(L) (M) (H) (U)
Figure 10. Segmentation results for the poplar mock-up with occlusions. The letter indicates the resolution level. On each point cloud, all points with the same colour belong to the same cluster.

Evaluation

Segmentation accuracy

The number of points correctly assigned to each elementary unit has been retrieved on the poplar mock-up. We call false positive (FP) for a given cluster a point that is labelled as part of this cluster by the segmentation algorithm, while it does not belong to this cluster in the input poplar mock-up point cloud. A point is a cluster false negative (FN) if it is not labelled as part of this cluster, while it actually belongs to it. A cluster's false positive rate (FPR) is the ratio of false positives over the actual number of points in the cluster. We define false negative rates (FNR) the same way.

It is worth mentioning that all points were labelled by the algorithm. This is because the constructed graph contains all points of the point cloud, and the spectral segmentation algorithm browses the whole graph. Figure 15 (a) shows that a large majority of leaf points are assigned to the correct leaf cluster, the worst case being cluster 5 (the bottom red leaf on Figure 10) in the ultra high resolution level point cloud with 5.77% of points assigned to another cluster. As shown on Figure 15 (b), leaf false positive rates are similar to false negative rates. The maximum is reached for cluster 2 in the ultra high resolution point cloud, which correspond to the bottom dark blue leaf on Figure 10, with 5.84% of false positives. False positives and negatives usually occur near the junction of a leaf to its petiole. The segmentation is not always accurate for the petioles and the internodes. This is explained by the fact that several sets of internodes and/or petioles are not fully segmented into elementary units, thus points of different internodes are assigned to the same cluster. This is for example the case of clusters 35 to 39, which correspond to internodes of the poplar's basis stem. Since no geometrical feature enables to split the stem into its internodes, and since the algorithm does not use any botanical knowledge, points are not segmented into internode clusters and remain in one global cluster, in green on Figure 10. Results are summarised in Table 4, in which we have computed means and standard deviations of the number of points over leaf, petiole and internode clusters, respectively. It shows that for leaves, false positive and negative rates remain below 3.4%. However, since it is difficult to unravel some petioles or internodes to their adjacent units from a pure geometrical point of view, points of neighbouring petiole or internode clusters are often pooled together. As a result, many petiole or internode clusters have no point assigned, leading to huge false positive and negative rates. Table 4 also shows that the resolution level has little impact on the segmentation accuracy, although results are slightly better for low resolution point clouds than for high resolution ones.

Statistical analysis of the leaf area estimates

Results of the leaf area estimates for the poplar mock-up are shown in Table 5. They show that the resolution level has a stronger influence on leaf area estimates than our segmentation method. Our estimates are always close to the estimates computed for the correct clusters. Sampling underestimates leaf areas because our area estimation method creates a piecewise linear surface which boundary is defined by points labelled as belonging of the leaf. Since these points are actually fully inside the leaf and not on its boundary, and since leaves of the poplar mock-up are approximated by convex flat surfaces, the computed surface is smaller than the actual one. The higher resolution, the smaller underestimate, since the boundary points for the Delaunay triangulation are closer to the actual leaf boundary. In case of occlusions some points may be missing in a leaf cluster, leading to a smaller surface estimate, thus again an underestimate of the leaf area.

It can also be noticed on Table 5 that our approach tends to slightly overestimate leaf areas with respect to the estimate of the actual segmentation. This is mainly due to the fact that a false positive point may easily add a large area to the estimate, since the Delaunay triangulation will create big triangles between this point and other points in the cluster. This is a counterbalancing effect to the underestimates of the sampling and the resolution.

Sensitivity analysis

We now detail some experiments on the sensitivity of the method to the three parameters. The algorithm has been run on the poplar mock-up with different values for all three parameters, see Tables 6, 7 and 8. We have computed the false positive and false negative rates for each set of parameters, as well as the variation of the estimated total leaf area (-1% means that the estimated total leaf area is 1% lower than the actual leaf area, which is 0.02617m 2 ). According to the experiments made (Table 6), the total leaf false positive (FPR) and negative (FNR) rates and signed leaf area errors only vary by 0.05%, 0.12% and 2.18%, respectively, with respect to the angle a. As explained in Section 2.1.6, computation time should be affected by the number of edges in the graph, which in turn depends on the value chosen for the angle parameter a. However, as shown in Table 6, although the number of edges exponentially decreases with the angle (see also Figure the total computation time does not vary much with a. In particular, the computation time for stage 3 is always 91s. This the theoretical computational complexity analysis (Section 2.1.6). We explain this counter-intuitive result by the fact that in practice, since our data is a set of elongated shapes in the embedding space, Dijkstra's algorithm does not update the shortest paths much and many edges of the graph are not used. Its complexity in practice is thus close to O(n log n) rather than O(m + n log n).

As a conclusion, the method is rather insensitive to parameter a. However, in the case the graph is to be stored in a file, we advise to choose a value of a = 90 • to reduce its size (see Figure 16). According to our experiments (not shown here), a value of a greater than 90 • may lead to a disconnected graph.

The total computation time linearly increases with respect to the number d of intrinsic directions (Table 7). This parameter does not affect much the total false positive and negative rates, which only varies by 0.07% and 0.17%, respectively, and the total leaf area error, which only varies by 1.03%. Therefore, it is not necessary to set a high number of intrinsic directions. Our experiments indicate that d = 5 or d = 10 are good guesses in most of the cases.

Our experiments (Table 8, Figure 9) show that choosing a large number c of clusters may lead to over-segmentations of leaves. On the poplar mock-up, the bottom leaf is segmented in two different clusters from c = 9 (not visible on Figure 9 since this leaf is side-view), and this is also the case for a second leaf from c = 25. To overcome this problem, we suggest to first set a small value for c. According to our experiments, c ∼ 25% of the total final number of clusters is generally a good guess. If some elementary units are nonetheless over-segmented, we provide a graphical interface to easily select and merge the corresponding clusters.

For the examples shown on Figure 8, c was set to 23%, 27% and 26% of the final number of clusters, respectively (c = 11, 9 and 5 for 48, 33 and 19 final clusters).

Computation time linearly increases with respect to the number c of desired clusters, as shown in Table 8. This is consistent with the previously explained computational complexity analysis (Section 2.1.6). Note that indicated computation times are for the first iteration only.

Discussion and conclusion

We presented here a semi-automatic method to segment a TLS point cloud of a small plant into its elementary units (internodes, petioles and leaf-blades). Qualitative results on four real tree seedlings show that such small scale plants are successfully split into leaf, petiole and stem components. The only two exceptions are compound leaves which are segmented in their leaflets, and adjacent internodes on a stem which may not be separated if no geometrical feature (bud, incident stem) is available. The method does not use any prior botanical knowledge, therefore can be applied in a wide variety of cases. Quantitative results on leaves show that the method is robust (around 1% labelling error) and leads to accurate leaf area estimates.

Only three parameters are used by the method. Only one of them, namely the desired number c of clusters, has an actual influence on the results. According to our tests on four different seedlings, results are insensitive to the branching structure and the leaf anatomy. Acquisition noise during the scanning process is also robustly handled, as shown on the red oak scans. The scan resolution also has little effect on the segmentation result, but has a strong influence on the leaf area computation.

Since no botanical knowledge is used by the algorithm, computed elementary units are not explicitly labelled as leaves, petioles or leaf-blades. This could be done in an additional step with a Principal Component Analysis as in [START_REF] Belton | Processing tree point clouds using gaussian mixture models[END_REF] or feature based histograms as in [START_REF] Paulus | Surface feature based classification of plant organs from 3d laserscanned point clouds for plant phenotyping[END_REF].

The proposed method being semi-automatic, it is suited for small plants such as tree seedlings but may be time consuming for more complex structures. In order to enhance the quality of the segmentation with a large number c of clusters, thus to reduce the interaction time for large scale trees, two improvements are planned in the future. First, we are currently working on the correction of the acquisition noise during the scanning process, in order to reduce the number of points sparsely sampled between elementary units. Filtering the input scans in a pre-processing step, which has been done in this paper, is not a perfect solution since it removes points and thus leads to underestimates of the leaf areas. Second, we plan to enhance the graph construction process (first step of the algorithm), in order to decrease the number of edges between two non adjacent elementary units from a botanical point of view (e.g., two leaves, as in Figure 6). Then, the algorithm will be tested on more complex structures such as full-scale trees.

Figure 1 .

 1 Figure 1. Pipeline of the segmentation method.

Figure 2 .

 2 Figure 2. Example of the spectral embedding process. (a) Input scan. (b) Computed graph. (c,d) Point cloud embedded into a 3-D spectral space (two different views; see also the accompanying video). Colours are set to roughly indicate the elementary units.

Figure 3 .

 3 Figure 3. Examples of graph construction for the methods: (a) ε-neighbourhood of a point in blue, for two different values of ε. (b) k-nearest-neighbours, for k = 2 and k = 5. (c) Proposed method, for a = 45 • .

Data:

  Point cloud P , a point p ∈ P , a user-chosen angle parameter a (in radian) Result: Set E(p) of the edges of the graph with endpoint p E(p) := ∅; Compute the k nearest neighbours of p in P , and put them in a priority queue Q ordered by increasing distance to p; for p ∈ Q do if ∃e ∈ E(p) such that angle(pp , e) < a then Discard p ; end else Put the edge pp in E(p); end end Algorithm 1: Building the neighbouring edges of a single point p in the cloud a connected graph with as few as possible redundant edges.

Figure 4 .

 4 Figure 4. Example of graph reconstruction by (a) the ε-neighbourhood method (with ε = 0.006m); (b) the k-nearest-neighbours method (with k = 8) and (c) the method used in this study (a = 90 • ).

  space), desired number c of clusters Result: Segmentation of V into disjoint sets {Cluster[1], . . . , Cluster[c]} Source := farthest point to a random point of G; i := 1; Segment[i] := ComputeShortestPaths(Source,G); while i < c do Segment[i + 1] := ComputeShortestPaths(Segment[1..i],G); p := point of Segment[1..i] connected to Segment[i + 1]; j := number of the segment to which belongs p; Remove successive points of Segment[j] from p to one of its end and add them to Segment[i + 2]; i+ = 2; end ComputeShortestPaths(Segment[1..c],G); for p ∈ V do p := closest point of Segment[1..c] from p; j := number of the segment to which belongs p ; Add p to Cluster[j]; end Algorithm 2: Proposed graph segmentation method (applied in spectral space).

Figure 5 .

 5 Figure 5. The segmentation process. (a) Input graph (in spectral space) and selected source point (in red, together with the path from the initial random point). (b,c,d) Computation of successive segments. (e) Shortest paths from each remaining point to the segments. (f ) Computed clusters.

  Figure 6. Points A and B are on two different leaf-blades of the poplar mock-up with high TLS resolution level (see Section 2.2.1). Their commute-time distance in the graph is large while their Euclidean distance is small.

Figure 7

 7 Figure 7. 3-D point cloud images of the poplar seedling mock-up used in this study for the low (L), medium (M), high (H) and ultra high (U) TLS resolution levels (Table2). The last image shows the point cloud generated without simulation of the occlusion: all hits from the laser source to objects were recorded.

Figure 8 .

 8 Figure 8. Segmentation results after the first run for (a) the poplar mock-up with occlusions, (b) the sweet chestnut and (c) the red oak, with c = 11 (a), c = 9 (b) and c = 5 (c), respectively.

Figure 9 .

 9 Figure 9. Result of the first iteration of the algorithm for various numbers c of clusters.

Figure 11 .

 11 Figure 11. Segmentation results for point clouds of four contrasting plant seedlings with three different TLS resolution levels each. From top to bottom: birch, horse chestnut, sweet chestnut and red oak. On each point cloud, all points with the same colour belong to the same cluster.

Figure 12 .Figure 13 .Figure 14 .Figure 15 .

 12131415 Figure 12. Segmentation results (top view). (a) Birch. (b) Horse chestnut. (c) Sweet chestnut. (d) Red oak. All are high TLS resolution level point clouds.

  Figure 16. Number of edges in the graph and file size with respect to the chosen angle a.

Table 1 .

 1 -2 m 2 ) Structure parameters of the tree seedlings.

	Poplar mock-up	0.462	17	2.617
	Birch	0.650	7	1.296
	Horse chestnut	0.607	9	16.567
	Sweet chestnut	0.465	19	4.530
	Red oak	0.547	10	5.553

Table 2 .

 2 Characteristics and settings of the Leica HDS-6100 terrestrial laser scanner used in this study for point cloud acquisitions.

	Single shot phase-shift technology with single return signal
	Wavelength (nm)			650 -690
	Range (m)		0.3 -50 at 18% albedo
	Spot size at exit (m)			0.003
	Beam divergence ( • )			0.0126
	Pre-set scanner resolution levels	Low (L) Medium (M) High (H) Ultra high (U)
	Angular sampling resolution ( • )	0.072	0.036	0.018	0.009
	Maximum point spacing at 3m (m) 0.0038	0.0019	0.0009	0.0005

Table 3 .

 3 Number of points for each point cloud (after filtering).

Table 4 .

 4 Number of points and false positives (FP) and negatives (FN) for the poplar mock-up. SD stands for standard deviation.

	Resolution	Occlusions Segmentation Slope Intercept	RMSE	Bias	Mean LA (± SD)
					(10 -4 m 2 ) (10 -4 m 2 ) (10 -4 m 2 )	(10 -4 m 2 )
		No	Actual	0.99	-1.17	0.28	-1.34	14.06 (7.7)
	Low (L)	Yes	Actual	1	-1.63	0.38	-1.64	13.75 (7.7)
		Yes	Computed	1	-1.19	0.90	-1.25	14.15 (7.8)
		No	Actual	1	-0.33	0.10	-0.39	15.01 (7.7)
	Medium (M)	Yes	Actual	1	-0.53	0.14	-0.52	14.87 (7.7)
		Yes	Computed	1.01	-0.64	0.37	-0.43	14.96 (7.9)
		No	Actual	1	-0.09	0.02	-0.13	15.27 (7.7)
	High (H)	Yes	Actual	1	-0.21	0.08	-0.18	15.21 (7.8)
		Yes	Computed	1	-0.20	0.15	-0.17	15.22 (7.8)
		No	Actual	1	-0.04	0.01	-0.04	15.35 (7.7)
	Ultra high (U)	Yes	Actual	1	-0.09	0.04	-0.06	15.33 (7.8)
		Yes	Computed	1.03	-0.37	0.55	+0.12	15.51 (8.0)

Table 5 .

 5 Statistical analysis of leaf area estimates on the poplar mock-up, for all resolution levels. Mean LA (± SD) control = 15.39 (7.7) 10 -4 m 2 . SD stands for standard deviation.

Table 6 .

 6 Influence of parameter a on the poplar mock-up (H), with d = 10 and c = 11. TLA stands for total leaf area.

Table 7 .

 7 Influence of parameter d for the poplar mock-up (H), with a = 90 • and c = 11.

	c	5	7	9	11 15 19 25 31
	Nb. of overseg. leaves	0	0	1	1	1	1	2	2
	Computation time (s) 329 340 351 361 384 409 441 478

Table 8 .

 8 Influence of parameter c for the poplar mock-up (H), with a = 90 • and d = 10.

Acknowledgements

The authors would like to express their sincere gratitude to the Forestry Commission, the University of Grenoble Alpes and Inria for funding this work, as well as to Rémy Cumont for his participation in coding the segmentation algorithm and Dr Elisa Hétroy-Wheeler for proof-reading the paper.

Funding

The Forestry Commission, the University of Grenoble Alpes (through an AGIR project) and Inria (through the Action de Recherche Collaborative PlantScan3D).