Abstract
Background
In recent years, imaging based, automated, noninvasive, and nondestructive highthroughput plant phenotyping platforms have become popular tools for plant biology, underpinning the field of plant phenomics. Such platforms acquire and record large amounts of raw data that must be accurately and robustly calibrated, reconstructed, and analysed, requiring the development of sophisticated image understanding and quantification algorithms. The raw data can be processed in different ways, and the past few years have seen the emergence of two main approaches: 2D image processing and 3D mesh processing algorithms. Direct image quantification methods (usually 2D) dominate the current literature due to comparative simplicity. However, 3D mesh analysis provides the tremendous potential to accurately estimate specific morphological features crosssectionally and monitor them overtime.
Result
In this paper, we present a novel 3D mesh based technique developed for temporal highthroughput plant phenomics and perform initial tests for the analysis of Gossypium hirsutum vegetative growth. Based on plant meshes previously reconstructed from multiview images, the methodology involves several stages, including morphological mesh segmentation, phenotypic parameters estimation, and plant organs tracking over time. The initial study focuses on presenting and validating the accuracy of the methodology on dicotyledons such as cotton but we believe the approach will be more broadly applicable. This study involved applying our technique to a set of six Gossypium hirsutum (cotton) plants studied over four timepoints. Manual measurements, performed for each plant at every timepoint, were used to assess the accuracy of our pipeline and quantify the error on the morphological parameters estimated.
Conclusion
By directly comparing our automated mesh based quantitative data with manual measurements of individual stem height, leaf width and leaf length, we obtained the mean absolute errors of 9.34%, 5.75%, 8.78%, and correlation coefficients 0.88, 0.96, and 0.95 respectively. The temporal matching of leaves was accurate in 95% of the cases and the average execution time required to analyse a plant over four timepoints was 4.9 minutes. The mesh processing based methodology is thus considered suitable for quantitative 4D monitoring of plant phenotypic features.
Background
In the coming decades, it is expected that mankind will need to double the quantity of food and biofuel produced in order to meet global demand [1]. To achieve this with existing resources, new plant characteristics need to be identified, quantified, and bred to obtain more productive plant varieties within existing environments. This will require a greater understanding of how the genetic makeup of plants determines their phenotype (visible traits) in high resolution and in high throughput. Performing plant phenomics involves screening large germplasm collections to facilitate the discovery of new interesting traits (forward phenomics), and analysing known phenotypic data in order to uncover the genes involved in their evolution and use these genes in plant breeding (reverse phenomics) [1]. Investigated plants are usually grown in thoroughly controlled conditions (growth chambers or glasshouses) and subjected to different environmental conditions and stresses (e.g. drought, salt, heat, etc.) with the primary aim of monitoring their phenotypic response using various measurements [2,3].
Common plant morphological traits of interest include parameters such as main stem height, size and inclination, petiole length and initiation angle, and leaf width, length, inclination, thickness, area, and biomass [14]. The usual procedure to collect these data consists of many laborious manual measurements, often requiring destructive harvests and thus multiple replicates of individual plant genotypes or varieties to allow successive harvests over time. A typical manual phenotypic analysis of 200 plants (daily objective) would require approximately 100 manhours of work (≃ 30 minutes per plant depending on the size and complexity), which is impractical. In light of the importance of gene discovery and agricultural crop improvement, the development of solutions to automate such a tedious task is imperative.
Highthroughput plant phenotyping aims to extend the standard approach by growing, measuring and analysing temporally thousands of plants [5]. In recent years, the plant phenotyping research has seen the emergence of highthroughput plant screening facilities [1,6]; however, few image and mesh processing solutions are available to analyse the large amount of data captured and extract yield determinants (i.e. plant, leaf, or root characteristics). Among existing solutions, PHENOPSIS [7] and GROWSCREEN [8,9], provide 2D imageprocessing based semiautomated solutions for leaf phenotyping (leaf width, length, area, and perimeter) and root data monitoring (number of roots, root area, and growth rate). LAMINA [10], another 2Dimage based tool for leaf shape and size probing proposes a leaf analysis for various plant species. Recent imageprocessing solutions, such as TraitMill [11] and HTPheno [12], provide a more general plant analysis and measure information such as plant height, width, centre of gravity, projected area and biovolume, and provide colorimetric analysis (e.g. greennessdifferences between plants). Due to the importance of rice as a primary food resource, imagebased solutions for rice phenotyping have been developed [6,13] and involve the measurements of parameters such as grain size (length, width, and thickness), panicle length, and number of tillers. In the past 2 years, fully automated imaging techniques for the highthroughput investigation of plant root characteristics (yield determinants) have been developed [1416] to analyse nondestructively phenotypic traits such as root average radius, area, maximum horizontal width, and length distribution.
The latest applications have introduced a third dimension to the plant analysis. Stereoimaging and mesh processing based systems, such as GROWSCREEN 3D [17], the 3D imaging and RootReader3D software platform[18], or the solution proposed in [19], have pioneered the explicit 3D analysis of leaves and roots, allowing more accurate measurements of leaf area, and extraction of additional volumetric data.
To date, the literature is distinctly dominated by 2D imageprocessing techniques for highthroughput phenotyping of plants [616]. The major limitation of these 2D solutions is the loss of crucial spatial and volumetric information (e.g. thickness, bending, rolling, orientation) when transposing available data from 3D to 2D. The recent introduction of new tools for plant analysis based on explicit 3D reconstructions [1719] (as opposed to inferred 3D based analysis [20,21], widely used since the 1960’s) promises to increase potential of highthroughput studies in terms of accuracy and exhaustiveness of the measured features, but available threedimensional solutions are currently focussed on a specific organ (e.g. leaves [17,19] or roots [18]), tailored to a particular image acquisition system [22], and tend to be qualitative (or applied) rather than providing quantitative information and estimates of accuracy. Hence, a clear need exists for a more generalised plant analysis based on increasingly explicit 3D models and in which the reliability of the measurements is questioned and quantitatively assessed.
In this paper, we present a novel meshbased technique developed for the highthroughput 3D analysis of plant aerialparts. A focus is made on the feasibility of accurately extracting plant phenotypic parameters from a 3D mesh acquired for the dicotyledonous crop cotton. In this initial study, meshes were reconstructed using a low cost commercial 3D reconstruction system [23]. The proposed methodology aims at a nonexhaustive, accurate, crosssectional (observation of a representative subset of a population at a fixed timepoint), and temporal investigation of the plant macroscopic phenotype. This requires advanced features such as plant mesh morphological segmentation [24,25], accurate plant data extraction [26], and plant organs tracking overtime. The mesh based methodology was tested on plant meshes reconstructed [23,27] for a set of six plants studied at four timepoints (i.e. 6×4 = 24 plant meshes).
Methods
To investigate the feasibility of a mesh based processing pipeline for the 3D analysis of plants, we initially developed the prototype crosssectional pipeline described on Figure 1. The following sections propose a brief presentation of the image acquisition and plant mesh reconstruction steps, and detailed descriptions of the mesh segmentation, temporal phenotypic analysis, and validation scheme.
Figure 1 . Prototype plant analysis pipeline. The horizontal axis (1) illustrates our prototype crosssectional (observation of a representative subset of a population at a specific timepoint) pipeline which is divided into 4 main stages: plant image acquisition, surface mesh reconstruction, morphological mesh segmentation (see Figure 2), and plant phenotypic parameters extraction. Monitoring the phenotypic parameters overtime (2) involves repeating the crosssectional pipeline throughout all the timepoints available and analysing the variations in the plant phenotype (the temporal pipeline is detailed on Figure 3).
Figure 2 . Automated segmentation pipeline. Illustration of the 4 successive steps of the segmentation pipeline: (1) coarse segmentation  takes as input the non segmented plant mesh and provides as output a plant mesh roughly segmented (1 region for the main stem (red) and N regions for the associated petioles and leaves (other colours)), (2) precise stem segmentation  takes as input the main stem previously segmented (region M) and partitions it into several internodes, (3) petiole segmentation  takes as input a given associated leaf and petiole (from step 1) and separates them into two distinct regions (process repeated for all the associated leaves/petioles), and (4) leaf segmentation  takes as input a given separated leaf and segments it into adaxial/abaxial surfaces and left/right parts (process repeated for all the leaves).
Figure 3 . Temporal analysis pipeline.The full temporal analysis involves repeating the same pairwise matching pipeline throughout all the time points. The pairwise matching (horizontal axis) is constituted of two steps: (1) a plant alignment, that takes as inputs the segmented plant meshes for the timepoints T_{x}and T_{x + 1}(see (a)) and rigidly aligns the mesh from timepoint T_{x}with the mesh from timepoint T_{x + 1}(see resulting mesh T_{x}in blue on (b)), and (2) a plant parts matching that matches the different plant organs using the aligned meshes (see example of matched leaves on (c)). The full temporal matching (vertical axis) consists of repeating this pairwise matching throughout all the timepoints in order to obtain the mapping of the different organs of the plant overtime.
Plant material
The prototype study involved acquiring and processing images for an initial set of six Gossypium hirsutum plants studied over four timepoints. Manual measurements, performed by X.S and S.B for each plant and each timepoint, were used to validate the accuracy and quantify the error on the meshbased phenotypic data estimation. The first three timepoints involved measuring invasively (but nondestructively) parameters such as main stem height, leaf width and leaf length using a measuring tape. For the last timepoint, measurements were collected after destructive harvest in order to optimise their precision. The petioles and leaves were cut from the main stem, laid flat on a table and carefully measured. Overall, a set of 384 measurements was manually collected (24 main stem height measurements, 180 leaf width measurements, and 180 leaf length measurements).
Plant data acquisition and mesh reconstruction
A manual data capture process similar to that described in [12] was used to collect multiple plant images from different viewing angles using a highresolution Pentax K10 SLR camera with a sigma 2040mm aspherical lens. Each cotton plant pot was placed at the centre of a rotating tray. The camera was fixed on a tripod during all the acquisition process. The rotating tray was manually turned and pictures were taken at each rotation angle (every degree). The acquisition process completed, 64 images were available (per plant and timepoint) for the multiview 3D reconstruction. An example of acquired plant image is shown on Figure 1, the image resolution was 3872x2592 pixels (≃ 10 Megapixels).
Plant 3D models (meshes) were created from the highresolution images using 3DSOM, a commercial 3D digitisation software [23]. The number of polygons constituting the reconstructed meshes fluctuated between 120000 and 270000.
The acquisition and mesh generation are not the primary focus of the current paper, however we acknowledge the “semiautomated” steps involved. An automated image acquisition platform [28] and a mesh reconstruction algorithm (based on [2931]) are under development and will allow full automation for future experiments.
Automated plant mesh segmentation
The identification of different plant organs is a critical stage in performing meshbased plant phenotyping and has proven problematic with 2D based image analysis solutions [1]. To complete this task, we developed an advanced mesh segmentation algorithm that partitions the plant mesh into morphological regions.
Mesh segmentation algorithms involve assigning a unique value (called a label) to all the points of the mesh (called vertices) that belong to the same region. A surface mesh is constituted of triangles that link the vertices together through their edges. Two vertices are said to be topologically connected (neighbours) if they share the edge of a triangle. Finally, a vertex comprises a normal vector equal to the average of the normal vectors of the neighbouring triangles.
Due to the complex and irregular morphology of plants, no generic mesh segmentation algorithm [24,25] is accurate and robust enough to identify the different plant parts (main stem, petioles, leaves). This paper introduces a “hybrid” segmentation pipeline that overcomes the morphological shape differences between cotton plants and various reconstruction inconsistencies due to occlusions (the most common being missing petioles, as they are occluded by the leaves in the images used by the reconstruction scheme). Our automated segmentation pipeline, illustrated in Figure 2, is constituted of 4 successive steps: a coarse segmentation, a stem segmentation, petioles segmentations, and leaves segmentations. All the operations described in the next paragraphs are fully automatic and do not require any manual input [28].
Step 1: Coarse segmentation
The purpose of this first step is to partition the plant into n+1 coarse regions (with n = number of leaves), one for the main stem (region M) and n for the pairs of petioles and leaves (regions N_{i}i=1,…,n). This is performed by a regiongrowing algorithm [24,32]. Regiongrowing algorithms start from a seed point (automatically selected based on prior criteria defined by the application) and gradually grow a region from neighbour to neighbour until a given criteria is met. Since the criteria to stop the growth of a region are userdefined, this generic approach is particularly convenient for coarse segmentations but often shows limitations when seeking accurate region delineation.
The scheme starts by defining a coarse region M as the main stem by fitting a curve c_{p} to the main stem from one extremity to the other and assigning to the region M all the vertices in a given planar radius of c_{p}. Remaining vertices are classified into the n regions N_{i}using a regiongrowing algorithm. The algorithm finds the first vertex that is not part of any region yet (at the start only one region is defined: M), uses it as seed point, and recursively grows a new region to all the eligible topological neighbours (creating a second region N_{1}). A neighbour is eligible if it does not belong to M or any of the regions N_{i}already created. The region stops to grow when there is no eligible neighbour remaining, i.e. all neighbours are labelled. The algorithm iterates through all the vertices of the mesh and grows a new region N_{i}each time it finds a vertex that does not belong to any of the regions M or N_{i} already created (new seed point). This scheme is robust to reconstruction issues such as holes in the mesh or detached mesh pieces, as a vertex does not need to be connected to the main mesh to become a seed. A typical result of this pass is shown in Figure 2.1.
Step 2: Main stem segmentation
The second segmentation step is based on a primitive fitting segmentation approach [33,34] that aims at refining the rough stem segmentation and partitioning the stem into different internodes (using the previously extracted region M). Primitive fitting algorithms involve finding a given shape (chosen based on the mesh structure) in a complex mesh and considering that all the vertices within the registered shape belong to the same region. In this work, the tubular shape fitting algorithm involves finding the tube parameters that minimise the point to surface distance to the region M.
The algorithm creates a tube around the curve c_{p}(see previous section) and optimises its radius to tighten the shape around the main stem (see Figure 4.a). Vertices inside the final tube fitted are considered part of the main stem region. The vertices bordering the tube (i.e. junctions between the petioles and the main stem) are used to partition M into different regions. A stem internode goes from one junction to the other (see Figure 4.a). Typical stem segmentations are shown in Figures 2.2 and 4.a.
Figure 4 . Schematic illustration of the segmentation.(a) shows the tubes used to segment the stem (blue) and one petiole (red). (b) illustrates the tube used to separate (segment) the petiole P_{i}from the leaf L_{i}. (c) illustrates the planar symmetry used to segment the leaf into two symmetric parts. In this particular case, the points p_{1}and p_{2}will belong to two different leaf regions as the angles α_{1}and α_{2}are signed differently .
Step 3: Petiole segmentation
The petioles are also segmented (and separated) from the leaves in each regions N_{i} (created in the first step of the segmentation) using tubular fitting. For each petiole and associated leaf in the regions N_{i}, we interpolate a curve along the petiole (using the local centre of mass of the vertices) and build a tube around it (see red tube on Figure 4.a/b). The tube follows the petiole and extends to the apex of the leaf. If we define B_{i} as the vertex outside the tube which is the closest to the main stem (leaf stalk), then all the vertices inside the tube which are closer to the main stem than B_{i} belong to the petiole region P_{i} (see Figure 4.b). Other vertices naturally belong to the leaf region L_{i}. In the case of a missing petiole (detected by a nontopological connectivity to the main stem), this step is skipped, and the region is processed using Step 4. Figure 2.3 illustrates a typical plant mesh after the petiole segmentations.
Step 4: Leaf segmentation
The leaf segmentation algorithm aims at obtaining a sagittal (left and right parts) and transversal (adaxial and abaxial surfaces: see Figure 5.d) segmentation. It has been designed to be robust to numerous natural leaf shape variations (bending, changes overtime) and erroneous leaf reconstructions due to occlusions during the 3D reconstruction (e.g. abnormal leaf thickness, leaves stuck together). A high accuracy of this stage is crucial for leaf width, length, area, and average thickness estimation. The proposed solution is based on two properties common to all the studied leaves: the symmetry and the vertices normal vector distribution (mainly pointing in two directions: away from the leaf adaxial or abaxial surface).
Figure 5 . Illustration of a leaf data extraction.(a) shows the leaf sagittal and coronal planes, and the points S_{i,1}, S_{i,2}, C_{i,1}, C_{i,2}which are used to compute the leaf width and length, (b) shows, in blue, the projection (onto the coronal plane) of the line fitted to the shape of the leaf from S_{i,1}to S_{i,2}. This projected line length is the estimate of the leaf width. (c) shows, in orange, the projection (onto the sagittal plane) of the line fitted to the shape of the leaf from C_{i,1}to C_{i,2}. This projected line length is the estimate of the leaf length. (d) represents a leaf transversally sliced into abaxial and adaxial surfaces (using the vertices normal vector clustering algorithm).
For the sagittal segmentation, a 2Dsymmetry based algorithm was found to be the most robust, accurate, and computationally efficient. It is obtained by projecting the vertices of the leaf onto the plane having the main stem axis as normal and comparing the sign of the angle between the vector going from the main stem (c_{p}) to the apex and the vector going from c_{p}to the considered vertex. The region to which a vertex belongs to depends on the sign of the angle (α_{1} and α_{2} in Figure 4.c). An illustration of this process is provided in Figure 4.c.
The transversal segmentation involves using a normal clustering algorithm that starts by computing as the average normal vector for the whole leaf L_{i}( will naturally point away from adaxial or abaxial surface). A first pass will sort the leaf vertices into two different regions (adaxial or abaxial) depending on the angle formed by their normal vector and (higher or lower than ). We then recompute using only the normal vectors of the vertices belonging to one of the two created regions and repeat the sorting process using the updated . This scheme is repeated until a natural convergence occurs (i.e. does not change between two iterations).
In the case where two leaves were merged together due to occlusion issues in the 3D reconstruction, the algorithm detects a largely greater volume, splits the region L_{i} into two regions, and performs the normal leaf segmentation on each region. To create the two regions, the algorithm detects the points f_{i,1} and f_{i,2} which are the furthest apart in the region L_{i}, computes the centroid f_{ci}of these two points, and uses as split plane the plane defined with f_{ci} as origin and the normalised vector from f_{ci}to f_{i,2} as normal. The mesh segmentation after this step is shown on Figure 2.4.
Phenotypic parameters of interest
For phenotypic analysis, important parameters are main stem height, size and inclination, petiole length and initiation angle, and leaf width, length, area, and inclination. This section presents the process used to extract these parameters from the segmented plant mesh, and focuses in particular on the leaf parameters, which are crucial indicators of the level of stress to which the plant is subjected to [2,3].
Main stem
The main stem height can be expressed as the height difference between the highest and lowest vertices of the region M. The normalised vector between these two vertices defines the main stem axis and the angle between this axis and the coordinates system upvector gives the inclination of the main stem. In this work the main stem length is defined as the length of the curve c_{p}fitted to the main stem.
Petiole
If c_{i} is a curve interpolated along the petiole P_{i}(using local vertices centre of mass), the length of the petiole can then be expressed as the length of c_{i}. In addition, if l_{i} and h_{i} denote the points of c_{i} that are the closest to c_{p} and the highest respectively, then the angle α between the main stem axis and the vector defines the petiole initiation angle (see αon Figure 4.a).
Leaf blade
For each segmented leaf L_{i}, we define L_{ci} as the centroid of the leaf, as the average of the vectors going from L_{ci} to the vertices belonging to the right part of the leaf, and as the vector going from L_{ci} to the tip of the leaf. Let and define the leaf sagittal and coronal planes (in which L_{ci}and are the origin and normal of the plane _{πi,x}) as displayed in Figure 5.a. Let S_{i,1}and S_{i,2} (resp.C_{i,1} and C_{i,2}) be the points on each side of π_{i,1} (resp.π_{i,2}) that maximise the distance to π_{i,1}(resp.π_{i,2}) as illustrated on Figure 5.b/c.
To estimate the leaf width (resp. length), we compute the length of the curve w_{i} (resp. l_{i}) interpolated to the leaf shape from S_{i,1}to S_{i,2} (resp.C_{i,1} to C_{i,2}) and projected onto π_{i,2}(resp.π_{i,1}) (to remove additional transversal length). Illustrations are provided in Figure 5.b/c. The projection of onto π_{i,1} is used as leaf axis. The angle between this axis and the main stem axis gives the leaf inclination. The leaf area can be estimated by averaging the areas of the adaxial and abaxial surfaces (see Figure 5.c), which are computed by summing the area of the triangles composing them. The leaf thickness can be estimated by averaging the distance between each vertex of the adaxial surface to the closest vertex on the abaxial surface.
Analysis overtime
This step of the pipeline involved monitoring the variations of the estimated plant parameters over time. Even though this is a straightforward process for stem height monitoring, the temporal petiole and leaf parameters analysis requires an efficient matching algorithm that tracks the different plant parts over time (orientation and size of the leaves change over time as a result of variations in growing conditions, making it difficult to find robust descriptors).
To perform this task, we developed the pipeline presented on Figure 3 which is based on the assumption that a plant organ position does not vary much between two close imaging dates. We apply the same “pairwise matching” algorithm (horizontal axis on Figure 3) throughout all the available timepoints (i.e. matching of T_{1} and T_{2}, of T_{2} and T_{3}, …, of T_{N−1} and T_{N}. See vertical axis on Figure 3) in order to obtain the sequences of leaves and petioles. The pairwise scheme is divided into two main steps: an alignment of the two plants and a parts matching algorithm.
Plants alignment
The plant at T_{x} is rigidly aligned with the plant at T_{x + 1}using a translation and a rotation around its main stem axis. The translation is performed using the vector going from centre of the plant at T_{x}(lowest point of the main stem region) to the centre of the plant at T_{x + 1}. The rotation is performed using the angle α that minimises the metric m_{α} defined in Eq. (1). In this equation , stands for the centroid of the leaf i, and for the number of leaves for the plant at timepoint T_{x}, and expresses the 3D Euclidean distance between and . An example of the alignment is shown in Figure 3.b.
Plant parts matching
Internodes matching
This step aims at matching the different internodes of a plant between two timepoints. Petioles may grow between two timepoints, and thus, split an internode of the plant at timepoint T_{x}into two parts, meaning that an internode from the plant at timepoint T_{x} can be matched with multiple internodes from the plant at T_{x + 1}. For the two plants, we rank the internodes by normalised height, retrieve their lowest and highest boundaries and match two together when there is an overlap between their boundaries. Figure 6 shows a schematic illustration of the process.
Figure 6 . Illustration of the process to match internodes. An internode from the plant at T_{x}can be matched with multiple internodes from the plant at T_{x + 1}in the case of a petiole appearance between the two timepoints.
Leaf blades and petioles matching
Using two aligned plants, we match the different leaves and petioles of the plants by solving an assignment problem. We build an adjacency matrix (comporting rows and columns) such that at a given position (ij) in the matrix, we store the distance between the centroids of the leaf i of the plant at T_{x} and of the leaf j of the plant at T_{x + 1}. Since the plants are now aligned, two leaves are eligible for pairwise correspondence only if they belong to a given angular range from each other, and we set to the distance between them in the adjacency matrix when this condition is not satisfied. The pairwise matching is performed using a simplified version of the Hungarian algorithm [35] that minimises the sum of the distances between the paired leaves. The petioles linked to the paired leaves are paired at the same time.
After this step, the morphological parts of the plant are matched over time.
Validation methodology
We were able to compute the relative error and the squared difference for a given automated measurement a_{i} with respect the manual measurement m_{i}by applying Eq. (2) and (3):
Let S_{a}={s_{a1},…,s_{aΩ}} and S_{m}={s_{m1},…,s_{mΩ}} denote the sets of automated and manual main stem height measurements (with Ω number of plants). Let also W_{a}={w_{a1},…,w_{aΦ}} and W_{m}={w_{m1},…,w_{mΦ}} denote the sets of automated and manual leaf width measurements (with Φ total number of leaves: 180). Using Eq. (2) and (3), we can then express the mean absolute percentage errors E_{s}, E_{w}, and the root mean square errors RMSE_{s}, RMSE_{w} on the main stem height and leaf width measurements by:
A similar analysis allows to compute E_{l}and RMSE_{l} for the leaf length measurements.
These errors were computed either using the whole datasets mentioned, or using the datasets trimmed from 10% of the outliers (5% of the best and worst relative errors). In addition, to be able to test the correlation between the automated and manual measurements, we calculated the squared Pearson productmoment correlation coefficient (^{R2}) [36] and the Intraclass Correlation Coefficient (ICC  Twoways random single measures) [3739]. The closer the ^{R2} and ICC coefficients are to 1, the stronger the correlation between two measurements.
Results
The results we obtained by applying our processing pipeline on the initial population of 6 Gossypium hirsutum plants studied over 4 timepoints are presented bellow.
Plant mesh segmentation
As illustrated by the segmentation results displayed in Figure 7, the segmentation pipeline performs well and meets its targeted expectation which is the identification of the morphological parts of the plant. The algorithm has proved to be robust to the mesh abnormalities caused by occlusions during the reconstruction step, including holes in the mesh and plant parts detached from the main mesh (typically, a missing petiole, see magnification on Figure 7 b/T2) or stuck together (see Figure 7 b/T3). The accuracy of the final mesh partitioning is limited by the quality of the 3D reconstruction and irregularities in the sagittal segmentation of a leaf could appear if the leaf is considerably rotated with respect to the main stem.
Figure 7 . Examples of plant segmentation resultsIllustration of the segmentation of two Gossypium hirsutum cotton plant (a) and (b) studied over 4 timepoints. As shown by the segmentation of the plant b at _{T2}and _{T3}the algorithm is robust to reconstruction issues such as floating mesh parts (typically leaves), and two leaves merges together.
Phenotypic parameters estimation
Our methodology allowed us to perform 384 phenotypic measurements on our initial population of 24 Gossypium hirsutum plant meshes (24 main stem height measurements, and 180 measurements for both leaf width and length measurements). Overall, the described methodology estimates plant phenotypic parameters with sufficient accuracy and reproducibility to be used as a surrogate for manual or destructive phenotypic analysis (Table 1). We noted mean absolute percentage errors of _{Es}≃9.34% (15.9 mm with an average stem height of 170.9 mm), _{Ew}≃5.75% (5.11 mm with an average leaf width of 88.9 mm), and _{El}≃8.78%(6.93 mm with an average leaf length of 78.9 mm) on the main stem height, and leaf width and leaf length measurements. The root mean square errors computed were RMS_{Es}≃19.043 mm, RMS_{Ew}≃7.287 mm, and RMS_{El}≃9.707 mm. These values provide the average difference between the mesh based measurements and the direct measurements in millimetres. The BlandAltman plot and the distribution of the relative error, presented in Figure 8.d/e, allow a more thorough analysis of the error and show that, even though most of the measurements were performed within an error range of [0%;10%] (see dotted red line on the Figure 8.e), many outliers remain in the analysis (vertical scattering on the BlandAltman plot). The presence of outliers is caused by imprecision in the mesh segmentation and/or erroneous plant reconstructions due to occlusions during the 3D reconstruction. The scatterplots and linear regressions displayed in Figure 8.a and 8.b allow to appreciate the strong correlations between the meshbased and direct leaf measurements. The plotted pointclouds fall into slightly scattered lines and the linear regressions are approaching the targeted reference (i.e. y=x). The squared Pearson correlation and intraclass correlation coefficients , IC_{Cw}≃0.974, , and IC_{Cl}≃0.967 calculated on the leaf width and length measurements concord with our previous statement of strong correlations as they approach 1. Finally, although the plot of Figure 8.c shows a more scattered pointcloud for the main stem height measurements, the correlation coefficients found were and IC_{Cs}≃0.941, which are acceptable precisions for our research.
Table 1. Main stem and leaf measurements analysis
Figure 8 . Statistical analysis: Comparison of meshbased measurements with manual groundtruthing(a), (b), and (c) present scatter plots of the different phenotypic parameters evaluated by our pipeline against manual measurements. The squared Pearson correlation and intraclass correlation coefficients computed for the main stem height, leaf width, and leaf length measurements were , , , ICC_{s}≃0.941, ICC_{w}≃0.974, ICC_{l}≃0.967, (d) is the BlandAltman plot of our datasets (i.e. the relative error against logarithm of the mean of two measurements), (e) illustrates the distribution of the error for each measurement type. The dotted red line represents the 10% relative error.
Temporal analysis
The automated temporal analysis of the plants was quite robust to the two major challenges: the growth of the plants over time and the changes in the topology and shape of the leaves over time. Correct matches of the different plant organs occurred in 95% of the cases (missing petioles were ignored). Illustrations of the results obtained by applying the pairwise matching pipeline for a plant studied over 4 timepoints are proposed in Figure 9.e. A dependency of the current matching scheme is that an organ needs to be accurately identified in order to be matched.
Figure 9 . Example of temporal analysis performed over 4 timepoints(a) presents the evolution of the number of leaves per plant, (b) illustrates the evolution of the main stem height for each plant in our initial set, (c) monitors the evolution of the average leaf area per plant, (d) monitor the evolution of the width of individual leaves for a given plant, and (e) illustrates the results of the pairwise matching process for a given plant (for each pairwise matching, the leaves matched are coloured similarly). These graphs illustrate an important contribution of our work, the temporal matching, as to date, no available solution allows to monitor the evolution of specific plant organs overtime while at the same providing a general analysis of the plant.
Computational cost
The automated 3D analysis was performed on a standard computer equipped with a processor Intel Core 2 Duo E8300 (2.83GHz). The mesh analysis involved an average execution time of 4.9 minutes and a total of 29.3 minutes for the complete analysis of six plants over four timepoints (with plant meshes decimated to 70000 triangles and without algorithm optimisations outside standard speeding techniques). Table 2 provides the full analysis of the computational cost of our algorithm. From Table 2 and Figure 9.a, we can denote that a higher the number of leaves will involve a slightly longer time to perform the analysis (due to repeated leaf segmentations and data computations, e.g. Plant 2). The time required to perform the full meshbased analysis (3D reconstruction excluded) is faster than any manual method can be, and the performed analysis provides additional information on the evolution of the data overtime. Future work will involve paralleling different stages of the algorithm (such as step 2 and 3 of the segmentation, the leaves segmentations, or the computation of the metrics from Eq. (1) in the plant alignment) and utilising cluster and cloud computing technologies to access more important computational resources.
Table 2. Analysis of the computational cost (in minutes)
Additional results
Additional results are provided in the website associated with this paper [see Additional file 1].
Addtional file 1. Website presenting the results. Website containing the results obtained by applying our method on the initial set of plant meshes. The different results are presented as tables containing links to the different webpages. The results of the segmentation and temporal matching between the different timepoints are available as images. Phenotypic parameters estimated by our method are available in the form of tables. In addition, a spreadsheet containing all the meshbased and manual measurements is available as a webpage and contains the statistical analysis presented in the paper.
Format: ZIP Size: 16.9MB Download file or display content in a new window
Discussion
As illustrated by Figure 9.a, 9.b, 9.c, and 9.d, that present a comparative study of the temporal evolution of phenotypic parameters for 6 Gossypium hirsutum plants, our methodology allows an accurate monitoring of the plants’ phenotypic traits overtime. By developing a hybrid mesh segmentation and analysis methodology for plant phenotyping, we have demonstrated that the automated temporal meshbased analysis of the plant aerial part is feasible (from the temporal broad plant analysis to the evolution of individual organs).
Nevertheless, our initial study has several limitations which should be acknowledge and will lead to further investigation and development.
As of today, the pilot study was limited in terms of the exhaustiveness of the phenotypic parameters estimated, but the explicit 3D reconstruction and robust identification of the morphological parts of the plant allow estimation of a large number of parameters of interest to plant biologists not easily extracted from 2D images with existing software platforms (accurate leaf area and biovolume rather than projected area, growth of individual leaves, organ quantification over time, leaf number / phyllochron, leaf angle). More phenotypic parameter extractions can be easily developed and incorporated to our pipeline as the biologists’ requirements evolve, allowing reuse of existing libraries of 3D models and the capacity to tailor the pipeline to new trait identification and quantification. Plant architecture is an important determinant of radiation use efficiency in crops and analysis of this trait in explicit 3D and over time has previously been an intractable problem with anything other than low throughput [1]. We should acknowledge, however, that tools for 3D analysis of roots based on inferred 3D “reconstructions” (i.e. 3D approximation using shapes such as tubes) exist and have been extensively used since the early 1960’s [20,21,40].
Although the methodology was solely tested on Gossypium hirsutum plants, it is expected that the method will be broadly and easily adaptable to other dicotyledonous crops such as canola, tomato, and low tillering monocotyledons with simple architectures such as corn. The pipeline can be easily adapted, and operators can be implemented and combined in order to increase the flexibility of the algorithm. Preliminary results (unpublished), obtained by reusing the two first steps of the segmentation pipeline (rough segmentation and stem segmentation) on corn, allowed to isolate the main stem, the leaves, and internodes, and allowed the direct computation of corn specific data. Due to the importance of rice and wheat as major food crops, the application of image based plant phenomics tools to grasses is of great interest. Significant development of our pipeline is needed to cope with occlusions due to the complex structure and the tillering observed in cereal crops. Their investigation will involve pushing the stateoftheart reconstruction and segmentation algorithms to their limit.
With respect to the accuracy on the phenotypic parameters, errors between 5 and 10% (involving ranges between 5mm to 7mm for the leaves) are acceptable for morphological scale phenotyping, reflecting the magnitude of errors already inherent in manual measurements and variations observed between individual plants of identical genetic makeup, and are low enough to distinguish changes in the relevant traits between two imaging dates during development (which is the aim of our research). Measurements for which the mean absolute error is above 10% (or over 10mm range, e.g. main stem measurements) will involve further work to improve the accuracy (for instance, the mean bias error  that characterises systematic over/under estimations  for the main stem height measurements was MBE_{s}≃9.8mm, against MBE_{w}≃−2.7mm and MBE_{l}≃3.1mm for the leaf width and length measurements, entailing that a systematic overestimation on the main stem height measurements is made). Finally, our current aim involves reducing the error on the measurements to less than 5%, which we believe is achievable by training our algorithms on phantom plant meshes (with phenotypic parameters exactly known) generated using existing plant modelling technologies [4,41].
While the focus of the current study has been the processing of meshes produced by a commercial 3D reconstruction product, major future work will involve improving the digitisation of plant structure and function by incorporating data other than visible light images into the 3D model. In addition to visible light cameras collecting multiple view geometries, PlantScan, a new screening platform recently developed in our laboratory [28] is equipped with LiDAR (Light Detection and Ranging sensors), infrared cameras, and multiwavelength cameras. The LiDAR cameras allow the reconstruction of accurate pointclouds (precision of 200 microns) that will be integrated in our probabilistic reconstruction scheme [2931] in order to improve the accuracy of the reconstructed plant meshes that currently limits the quality of the morphological segmentation and temporal analysis. These meshes will be overlaid with thermal infrared data and multispectral images data that provide colorimetric information (for chemical composition and photosynthetic functional analysis). Our laboratory expects to scan one plant every 7 minutes, making the current meshbased methodology (3D reconstruction excluded) suitable for highthroughput dicotyledonous plant analysis. As 3DSOM required an average processing time of 15 minutes to reconstruct suitable meshes, a special focus will be placed on the efficiency of the reconstruction scheme developed.
Conclusions
In this paper, we presented a hybrid meshbased methodology developed for highthroughput plant phenomics research. The proposed solution provides advanced meshprocessing features, including plant mesh morphological segmentation, accurate plant aerialpart phenotypic parameters estimation, and individual organ tracking and data monitoring overtime. Experiments involved testing our processing pipeline on an initial set of six Gossypium hirsutum plants analysed over four timepoints.
From the qualitative and quantitative results presented in the paper, we believe that the development of a mesh based methodology for highresolution and highthroughput plant phenomics platform is feasible and offers multiple advantages over current systems that use a small number of 2D images. The hybrid mesh segmentation presented allowed the identification of the different plant organs for all the test plants. The phenotypic parameter estimation algorithms allowed the retrieval of measurements such as main stem height and inclination, petiole length and initiation angle, and leaf width, length, area and inclination. By comparing 384 meshbased measurements with manual measurements, we observed errors ranging from 5.75% to 9.34% and correlations ranging from 0.887 to 0.974. The temporal organ tracking algorithm successfully matched plant organs between timepoints in 95% of the cases. Finally, the proposed analysis required only 4.9 minutes in average to analyse a plant over four timepoints. The meshbased analysis is thus considered a suitable mean to perform accurate and efficient 3D plant phenotypic analysis.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
XS, SB, and RF performed the manual data acquisition, 3D reconstruction and the manual measurements on the plants. AP and JF built the automated images and mesh based analysis pipeline (segmentation, data extraction, temporal analysis). AP, JF, and XS are the principal authors of this paper. All authors read and approved the final manuscript.
Availability and requirements
The different operators presented in this paper, as well as an initial set of plant meshes, are available for download from the PlantScan home page (Microsoft Windows 64bits installer): [28].
Acknowledgements
The analysis pipeline described in this manuscript was developed as part of an infrastructure grant received under the Commonwealth Government’s National Collaborative Research Infrastructure Strategy (NCRIS). AP was supported by a High Resolution Plant Phenomics Centre Studentship funded by the Australian Capital Territory Government, Canberra. Authors wish to thank Dr Francois Tardieu (INRA, Montpellier, FRANCE) for helpful comments and thoughtful review of the manuscript.
References

Furbank RT, Tester M: Phenomics  technologies to relieve the phenotyping bottleneck.
Trends Plant Sci 2011, 16(12):635. PubMed Abstract  Publisher Full Text

Granier C, Tardieu F: Multiscale phenotyping of leaf expansion in response to environmental changes: the whole is more than the sum of parts.
Plant Cell Environ 2009, 32(9):1175. PubMed Abstract  Publisher Full Text

Schurr U, Heckenberger U, Herdel K, Walter A, Feil R: Leaf development in Ricinus communis during drought stress: dynamics of growth processes, of cellular structure and of sinksource transition.
J Exp Bot 2000, 51(350):1515. PubMed Abstract  Publisher Full Text

Vos J, Evers JB, BuckSorlin GH, Andrieu B, Chelle M, de Visser PHB: Functionalstructural plant modelling: a new versatile tool in crop science.
J Exp Bot 2010, 61(8):2101. PubMed Abstract  Publisher Full Text

Eberius M, LimaGuerra J: HighThroughput plant phenotyping–data acquisition, transformation, and analysis.

Duan L, Yang W, Huang C, Liu Q: A novel machinevisionbased facility for the automatic evaluation of yieldrelated traits in rice.
Plant Methods 2011, 7:44. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Granier C, Aguirrezabal L, Chenu K, Cookson S, Dauzat M, Hamard P, Thioux J, Rolland G, BouchierCombaud S, Lebaudy A, et al.: PHENOPSIS, an automated platform for reproducible phenotyping of plant responses to soil water deficit in Arabidopsis thaliana permitted the identification of an accession with low sensitivity to soil water deficit.
New Phytologist 2006, 169(3):623. PubMed Abstract  Publisher Full Text

Walter A, Scharr H, Gilmer F, Zierer R, Nagel K, Ernst M, Wiese A, Virnich O, Christ M, Uhlig B, et al.: Dynamics of seedling growth acclimation towards altered light conditions can be quantified via GROWSCREEN: a setup and procedure designed for rapid optical phenotyping of different plant species.
New Phytologist 2007, 174(2):447. PubMed Abstract  Publisher Full Text

Jansen M, Gilmer F, Biskup B, Nagel KA, Rascher U, Fischbach A, Briem S, Dreissen G, Tittmann S, Braun S, De Jaeger I, Metzlaff M, Schurr U, Scharr H, Walter A: Simultaneous phenotyping of leaf growth and chlorophyll fluorescence via GROWSCREEN FLUORO allows detection of stress tolerance in Arabidopsis thaliana and other rosette plants.
Funct Plant Biol 2009, 36(11):902. Publisher Full Text

Bylesjö M, Segura V, Soolanayakanahally R, Rae A, Trygg J, Gustafsson P, Jansson S, Street N: LAMINA: a tool for rapid quantification of leaf size and shape parameters.
BMC Plant Biol 2008, 8:82. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Reuzeau C, Pen J, Frankard V, de Wolf J, Peerbolte R, Broekaert W: TraitMill: a discovery engine for identifying yieldenhancement genes in cereals.

Hartmann A, Czauderna T, Hoffmann R, Stein N, Schreiber F: HTPheno: An image analysis pipeline for highthroughput plant phenotyping.
BMC Bioinf 2011, 12:148. BioMed Central Full Text

Yang W, Xu X, Duan L, Luo Q, Chen S, Zeng S, Liu Q: Highthroughput measurement of rice tillers using a conveyor equipped with xray computed tomography.
Rev Sci Instrum 2011, 82(2):025102. PubMed Abstract  Publisher Full Text

Naeem A, French A, Wells D, Pridmore T: Highthroughput feature counting and measurement of roots.
Bioinformatics 2011, 27(9):1337. PubMed Abstract  Publisher Full Text

Yazdanbakhsh N, Fisahn J: High throughput phenotyping of root growth dynamics, lateral root formation, root architecture and root hair development enabled by PlaRoM.
Funct Plant Biol 2009, 36(11):938. Publisher Full Text

IyerPascuzzi A, Symonova O, Mileyko Y, Hao Y, Belcher H, Harer J, Weitz J, Benfey P: Imaging and analysis platform for automatic phenotyping and trait ranking of plant root systems.
Plant Physiol 2010, 152(3):1148. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Biskup B, Scharr H, Fischbach A, WieseKlinkenberg A, Schurr U, Walter A: Diel growth cycle of isolated leaf discs analyzed with a novel, highthroughput threedimensional imaging method is identical to that of intact leaves.
Plant Physiol 2009, 149(3):1452. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Clark R, Maccurdy R, Jung J, Shaff J, McCouch S, Aneshansley D, Kochian L: Threedimensional root phenotyping with a novel imaging and software platform.
Plant Physiol 2011, 156(10):455. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wang H, Zhang W, Zhou G, Yan G, Clinton N: Imagebased 3D corn reconstruction for retrieval of geometrical structural parameters.
Int J Remote Sensing 2009, 30(20):5505. Publisher Full Text

Huang Q, Stockman GC: Generalized tube model: recognizing 3D elongated objects from 2D intensity images. In PROCEEDINGS CVPR IEEE. , ; 1993:104109.

Huang Q, Jain A, Stockman G, Smucker A: Automatic image analysis of plant root structures. In Pattern Recognition, 1992. Vol.II. Conference B: Pattern Recognition Methodology and Systems, Proceedings., 11th IAPR International Conference on. , ; 1992:569572.

Nathalie W, JeanChristophe P: Highcontrast threedimensional imaging of the Arabidopsis leaf enables the analysis of cell dimensions in the epidermis and mesophyll.
Plant Methods 2010, 6(17):1. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Baumberg A, Lyons A, Taylor R: 3D SOM–A commercial software solution to 3D scanning.
Graphical Models 2005, 67(6):476. Publisher Full Text

Shamir A: A survey on Mesh Segmentation Techniques.
Computer Graphics Forum 2008, 27(6):1539. Publisher Full Text

Golovinskiy A, Funkhouser T: Consistent Segmentation of 3D Models.

Cornelissen JHC, Lavorel S, Garnier E, Díaz S, Buchmann N, Gurvich DE, Reich PB, Steege HT, Morgan HD, Heijden MGaVD, Pausas JG, Poorter H: A handbook of protocols for standardised and easy measurement of plant functional traits worldwide.
Aust J Bot 2003, 51(4):335. Publisher Full Text

Niem W, Wingbermuhle J: Automatic reconstruction of 3D objects using a mobile monoscopic camera. In 3D Digital Imaging and Modeling, 1997. Proceedings., International Conference on Recent Advances in IEEE. , ; 1999:173180.

Franco J, Boyer E: Fusion of multiview silhouette cues using a space occupancy grid. In Computer Vision, 2005. ICCV 2005. Tenth IEEE International Conference on, Volume 2 IEEE. , ; 2005:17471753.

Kolev K, Klodt M, Brox T, Cremers D: Continuous global optimization in multiview 3d reconstruction.
Int J Comput Vision 2009, 84:80. Publisher Full Text

Hosoi F, Omasa K: VoxelBased 3D modeling of individual trees for estimating leaf area density using highresolution portable scanning lidar.
Geoscience Remote Sensing, IEEE Transact on 2006, 44(12):3610.

Vieira M, Shimada K: Surface mesh segmentation and smooth surface extraction through region growing.
Comput Aided Geometric Des 2005, 22(8):771. Publisher Full Text

Attene M, Falcidieno B, Spagnuolo M: Hierarchical mesh segmentation based on fitting primitives.
Visual Comput 2006, 22(3):181. Publisher Full Text

Mortara M, Patané G, Spagnuolo M, Falcidieno B, Rossignac J: Plumber: a method for a multiscale decomposition of 3D shapes into tubular primitives and bodies. In Proceedings of the ninth ACM symposium on Solid modeling and applications. , ; 2004:339344.

Rodgers J, Nicewander W: Thirteen ways to look at the correlation coefficient.

Shrout P, Fleiss J: Intraclass correlations: uses in assessing rater reliability.

Bartko J: The intraclass correlation coefficient as a measure of reliability.

Danjon F, Reubens B: Assessing and analyzing 3D architecture of woody root systems, a review of methods and applications in tree and soil stability, resource acquisition and allocation.
Plant Soil 2008, 303:1. Publisher Full Text

Pradal C, DufourKowalski S, Boudon F, Fournier C, Godin C: OpenAlea: a visual programming and componentbased software platform for plant modelling.
Funct Plant Biol 2008, 35(10):751. Publisher Full Text