Image-Based Dynamic Quantiﬁcation of Aboveground Structure of Sugar Beet in Field

: Sugar beet is one of the main crops for sugar production in the world. With the increasing demand for sugar, more desirable sugar beet genotypes need to be cultivated through plant breeding programs. Precise plant phenotyping in the ﬁeld still remains challenge. In this study, structure from motion (SFM) approach was used to reconstruct a three-dimensional (3D) model for sugar beets from 20 genotypes at three growth stages in the ﬁeld. An automatic data processing pipeline was developed to process point clouds of sugar beet including preprocessing, coordinates correction, ﬁltering and segmentation of point cloud of individual plant. Phenotypic traits were also automatically extracted regarding plant height, maximum canopy area, convex hull volume, total leaf area and individual leaf length. Total leaf area and convex hull volume were adopted to explore the relationship with biomass. The results showed that high correlations between measured and estimated values with R 2 > 0.8. Statistical analyses between biomass and extracted traits proved that both convex hull volume and total leaf area can predict biomass well. The proposed pipeline can estimate sugar beet traits precisely in the ﬁeld and provide a basis for sugar beet breeding.


Introduction
Sugar beet ranks as the second largest crop for sugar production in the world and the first in Europe. It is also the source for bioethanol and animal feed [1,2]. According to the prediction of per capita sugar consumption and population growth, the world sugar production is about 175 metric tons at present, and needs to be increased by one million tons every year to meet the demand of 230 million tons in 2050 [3,4]. To improve global sugar production, plant breeding programs need to find the desired sugar beet genotypes to meet the requirements [2,5]. However, plant phenotyping, as a key to plant breeding programs [6,7], is usually performed manually which is time-consuming and expensive [8,9]. Thus, plant phenotyping becomes a bottleneck in many breeding programs.
To speed up and improve plant phenotyping, two-dimensional (2D) solutions have emerged in recent years, such as some semi-automated solutions based on 2D image for phenotyping of leaf and root [10][11][12][13]. However, due to the loss of one dimension information (e.g., thickness, bending, rolling and orientation) when transposing the data from three dimensions (3D) to 2D [14][15][16], these solutions struggled to estimate plant phenotypic traits precisely. Nowadays, several sensor-based 3D technologies have been used in plant phenotyping and can be divided into two categories: active and passive sensors [17]. Active sensors emit its own light source independent of external light, such as hand-hold laser scanning (LS) [18,19], structured light (SL) [20,21], terrestrial laser scanning (TLS) [22,23] and time of flight (TOF) [24,25]. These active technologies have good resolutions and can be applied to different phenotyping scales. However, they are expensive, especially for LS and TLS [17,26].
Structure from motion (SFM), featured as lightweight and cheap [17,26], is a commonly used passive technology in phenotypic traits quantification [27,28], biomass estimation [29,30] and yield prediction [31]. A previous study achieved high precision for leaf area and volumetric shape of taproots of sugar beet in the laboratory [18]. To our knowledge, quantification of 3D architecture of sugar beet growing in the field is still missing. The conditions in the field are much more complicated than that in laboratory, such as serious occlusion between plants, wind and light condition etc. Many technologies (e.g., SL, TOF) have difficulty in obtaining complete point cloud of individual plant caused by occlusions between plants or overlapping between leaves in the field. Structure from motion technique can reduce the effect of occlusion by adding pictures from different viewpoints. The more data collected from different points of view, the more complete the point cloud of individual plant obtained [17].
The overall aim of this study is to dynamically quantify aboveground structure of individual sugar beet plant growing in the field based on SFM method. The specific goals are to: (a) reconstruct 3D point cloud of individual plant based on SFM method; (b) develop an automatic data processing pipeline to process point cloud of individual plant and extract phenotypic trait values, such as plant height, maximum canopy area, convex hull volume, total leaf area and individual leaf length; (c) explore the possibilities of using the extracted phenotypic trait values for sugar beet biomass estimation.

Experimental Design and Data Acquisition
The experimental site is located at the Liangcheng experimental station, Inner Mongolia, China (40.502261 N, 112.145304 E) with an altitude of 1459.24 m. There were 20 planting plots in the experimental area of 88 m 2 (8 m × 11 m). The area of each plots was 1.2 × 2.3 m. Twenty genotypes of sugar beets were planted. The sowing time was 20 May 2019. The plant spacing was 0.25 m, and row spacing was 0.4 m. The spacing between the plots was 0.5 m. The field management of all plots was the same.
The sugar beet plants were sampled at 69 (T1), 96 (T2) and 124 (T3) days after emergence. In each period, 20 plants were randomly sampled. We used Canon 800D digital camera (Canon, Inc., Tokyo, Japan) to capture multi-view images of individual plant with 6000 × 4000 resolution in situ. Photographing was taken in two circles around the plant and the overlapping area between the two adjacent images was about 70-80%, in order to capture the feature points across-the-board from side view and top view ( Figure 1). After taking images, the manual measurements of individual plant were performed, including height and fresh biomass of aboveground. The plant height was measured as the vertical distance between the base of the stem to the highest point of the plant. Fresh biomass was obtained using an electronic balance.

Data Process
The process flow of multi-view images was shown as was taken in two circles around the plant and the overlapping area between the two adjacent images was about 70-80%, in order to capture the feature points across-the-board from side view and top view ( Figure 1). After taking images, the manual measurements of individual plant were performed, including height and fresh biomass of aboveground. The plant height was measured as the vertical distance between the base of the stem to the highest point of the plant. Fresh biomass was obtained using an electronic balance.

Data Process
The process flow of multi-view images was shown as Figure 2, including following steps on point clouds of individual plant: (1) Reconstruction and preprocessing (Figure 2a  of sugar beet after filtering; (g) Plant structure extraction of sugar beet, green part means maximum canopy area; (h) Point cloud of sugar beet after segmentation; (i) Leaf structure extraction, the red line represents the extracted leaf length.

Reconstruction and Preprocessing of Point Clouds
The 3DF Zephyr Aerial software was used to reconstruct point cloud of individual sugar beet from multi-view image set based on multi-view stereo and structure from motion (MVS-SFM) algorithm (Figure 2a,b). The SFM algorithm can generate sparse point clouds from multi-view image set. During the process, a scale-invariant feature transform (SIFT) keypoints detector [32], which is featured as invariant to image translation, scaling and rotation, was used to detect keypoints in multi-view images.

Reconstruction and Preprocessing of Point Clouds
The 3DF Zephyr Aerial software was used to reconstruct point cloud of individual sugar beet from multi-view image set based on multi-view stereo and structure from motion (MVS-SFM) algorithm (Figure 2a,b). The SFM algorithm can generate sparse point clouds from multi-view image set. During the process, a scale-invariant feature transform (SIFT) keypoints detector [32], which is featured as invariant to image translation, scaling and rotation, was used to detect keypoints in multi-view images. An approximate nearest neighbor (ANN) algorithm [33] was adopted to match keypoints from multi-view images. Then, projection matrix and 3D point coordinates were correspondingly calculated. The MVS algorithm can turn sparse point clouds into dense point clouds.
For reconstructed point cloud, the coordinate is unitless. A magic cube (Figure 2c) with 5.5 cm side length was adopted to obtain actual size of sugar beet point cloud. In addition, point cloud of individual sugar beet (Figure 2c) was extracted from sugar beet population using CloudCompare (http://www.danielgm.net/cc/; 3D point cloud and mesh processing software) manually.

Coordinates Correction of Point Cloud
In order to correct the coordinates of point cloud as z-axis perpendicular to the ground, we need to find the rotation matrix from the ground normal vector to the z-axis vector at current coordinate system. Then, the rotation matrix is applied to point cloud of individual plant. First of all, the ground plane equation (Equation (1)) of point cloud was calculated using the Random Sample Consensus (RANSAC) algorithm [34], listed as follows: Then, according to the ground plane equation, we can easily obtain the vector where m and n are the length of where E is a third order unit matrix, θ is a rotation angle and the unit vector of rotation axis where (x, y, z) is origin coordinate and (x1, y1, z1) is corrected coordinate.

Filtering and Segmentation of Point Cloud
As the influence of light, wind and other factors in the field environment the point cloud will inevitably contain some noise. The extraction of plant structure will be affected. Therefore, the filtering needs to be performed to remove noise points and outliers in order to obtain a smoother plant dense points cloud. There are three classes of points in reconstructed point cloud: plant, soil and shadows ( Figure 2e). We selected G minus R (G-R) as color filter and calculated the G-R mean values of point cloud in three periods (T1, T1, T3). We found that the color filter divided point cloud into three categories. The G-R values represent plants, shadows and soil from larger to small value ( Figure 3). The value 7 of (G-R) between plants and shadows points was chosen to eliminate soil and shadows points from plant. Outliers were removed based on statistical filter ( Figure 2f).  The region-growing algorithm [35] was utilized to segment plant into individual leaves ( Figure 2h). The algorithm grows from the point with the smallest curvature (seed point) and classifies neighbor points with small normal deviation into the same category. If the normal deviation of a neighbor point exceeds the threshold but the curvature deviation remains within the threshold, the point will be the next seed point. If the curvature deviation and normal deviation of a neighbor point both exceed threshold, the point will be discarded.
For plants with flat leaves, the region-growing algorithm performs well. While for plants with curved leaves, certain parts of curved leaves will be discarded. Therefore, clustering based on Euclidean distance was used to cluster discarded points on the leaves which they originally belonged to.

Structure Extraction of Plant Aboveground
Five phenotypic traits were extracted from the processed point cloud of individual plant, including plant height, maximum canopy area, convex hull volume, total leaves area and individual leaf blade length. The base of the stem to the highest point of the plant is defined as the plant height ( Figure 2g). Maximum canopy area is calculated from the product of the maximum length and maximum width of the canopy (the green part in Figure 2g). Total leaf area is the sum of individual leaf area of the whole plant. These traits can be classified into two categories: (a) global traits and (b) organ-level traits. Global traits quantify the overall plant morphology, such as plant height, maximum canopy area, convex hull volume and total leaves area. Organ-level traits focus on the individual organs of plant, such as leaf blade length.
As the coordinates of point cloud of individual plant have been corrected, we can estimate plant height using: The region-growing algorithm [35] was utilized to segment plant into individual leaves ( Figure 2h). The algorithm grows from the point with the smallest curvature (seed point) and classifies neighbor points with small normal deviation into the same category. If the normal deviation of a neighbor point exceeds the threshold but the curvature deviation remains within the threshold, the point will be the next seed point. If the curvature deviation and normal deviation of a neighbor point both exceed threshold, the point will be discarded.
For plants with flat leaves, the region-growing algorithm performs well. While for plants with curved leaves, certain parts of curved leaves will be discarded. Therefore, clustering based on Euclidean distance was used to cluster discarded points on the leaves which they originally belonged to.

Structure Extraction of Plant Aboveground
Five phenotypic traits were extracted from the processed point cloud of individual plant, including plant height, maximum canopy area, convex hull volume, total leaves area and individual leaf blade length. The base of the stem to the highest point of the plant is defined as the plant height ( Figure 2g). Maximum canopy area is calculated from the product of the maximum length and maximum width of the canopy (the green part in Figure 2g). Total leaf area is the sum of individual leaf area of the whole plant. These traits can be classified into two categories: (a) global traits and (b) organ-level traits. Global traits quantify the overall plant morphology, such as plant height, maximum canopy area, convex hull volume and total leaves area. Organ-level traits focus on the individual organs of plant, such as leaf blade length.
As the coordinates of point cloud of individual plant have been corrected, we can estimate plant height using: where z max and z min are the z-coordinate of the tallest and the shortest points among the point clouds of individual plant, respectively. The quick hull algorithm [36] was adopted to build a convex hull for point cloud of individual plant as shown with the dark green hull in Figure 2g). The calculation of each leaf area is based on the leaf mesh. The following steps were performed: (a) Use the Moving Least Squares (MLS) [37] algorithm for point cloud smoothing and improve normal estimation; (b) Build a mesh for leaf point cloud using the marching cubes algorithm [38]. (c) Calculate the sum of each triangle area using Helen's formula (Equation (8)) as the leaves' area: where S is the area of a triangle, p is the half of perimeter of triangle and a, b and c are the lengths of the three sides of the triangle. Accurate estimation of leaf blade length largely relies on the endpoints of leaf midrib. In this study, the leaf base point is defined as the point closest to the main stem. However, there is no obvious main stem for sugar beet. Thus, the line parallel to the Z axis and passed through the midpoint of the bounding box is defined as the main stem. After obtaining the base point, the tip point of leaf blade is defined as the point farthest from a line, which is parallel to the Z axis and passing through the base point. Once finding the base and tip points of the leaf, we needed to rotate the leaf base and tip points to the same axis to reduce leaf midrib dimension from 3D to 2D. First of all, the leaf blade was translated so that the base point of leaf blade was at the coordinate origin. The point cloud of leaf blade needs to be rotated α radians around the z-axis (anti-clockwise) so that the tip point is on the y-axis. The α can be calculated by following formula: where x and y are the x-axis and y-axis of tip point, respectively. Finally, polynomial regression-fitting algorithm [39] was used to smooth midrib of leaf blade and shown with the red line in Figure 2i.

Performance Evaluation
Sixty sugar beets with one or two leaves randomly selected from each plant were used for performance evaluation. Linear regression analyses were performed between the manual measurements and the estimation values for plant height and leaf blade length. The following statistics were calculated for performance evaluation: correlation coefficient (r) and root mean square error (RMSE, Equation (10)): where n is the number of samples; y i and x i are the estimation value and manual measurement, respectively.
Univariate linear regression was used to find the relationship between extracted plant traits (total leaf area and convex hull volume) and fresh biomass. The r value was calculated to evaluate the relationships between the fresh biomass and those estimation traits.  Figure 4 showed the histograms of G-R mean value of sugar beet at different growth periods. The color filter G-R performed well for separating sugar beet from soils and shadows at three periods. There were small differences of the separating thresholds among three periods (8.5 at T1 period, 6.5 at T2 period and 7.5 at T3 period). Therefore, we selected a G-R difference of 7 as threshold to filter all point cloud of individual plant at three periods.

Segmentation of Point Cloud
Point clouds of individual plant at three growth stages were presented in Figure 6. Both leaf number and total leaf areas increased rapidly from stages T1 to T2, as T1 is a fast-growing stage for sugar beets. When it came to stage T3, most of the leaves became curved and withered. In the current study, we classified 20 genotypes of sugar beets into two plant types according to their appearances ( Figure 6). If leaves were more scattered, they can be easily distinguished, especially in stage T3 (Figure 6a). If leaves gathered together, it was very difficult to separate them from each other,

Segmentation of Point Cloud
Point clouds of individual plant at three growth stages were presented in Figure 6. Both leaf number and total leaf areas increased rapidly from stages T1 to T2, as T1 is a fast-growing stage for sugar beets. When it came to stage T3, most of the leaves became curved and withered. In the current study, we classified 20 genotypes of sugar beets into two plant types according to their appearances ( Figure 6). If leaves were more scattered, they can be easily distinguished, especially in stage T3 (Figure 6a). If leaves gathered together, it was very difficult to separate them from each other, especially in stage T2 (Figure 6b).

Segmentation of Point Cloud
Point clouds of individual plant at three growth stages were presented in Figure 6. Both leaf number and total leaf areas increased rapidly from stages T1 to T2, as T1 is a fast-growing stage for Remote Sens. 2020, 12, 269 8 of 17 sugar beets. When it came to stage T3, most of the leaves became curved and withered. In the current study, we classified 20 genotypes of sugar beets into two plant types according to their appearances ( Figure 6). If leaves were more scattered, they can be easily distinguished, especially in stage T3 (Figure 6a). If leaves gathered together, it was very difficult to separate them from each other, especially in stage T2 (Figure 6b).   Using the regional-growing algorithm, point clouds of individual leaves were initially separated, shown in Figure 7. Within three stages, the segmentation method performed well at T1 and T3, as the leaves were scattered without much overlaps. However, the leaves were not segmented well at T2 stage. As for point cloud at T2 stage, the middle part of the leaves overlapped significantly (leaves in the blue circle in Figure 7b). Thus, it was hard to separate them into individual leaves. The points of the curved leaves were abandoned after initial segmentation, shown in Figure 7a with a red circle. Considering the close distance between points within the same leaf, the clustering algorithm was further utilized to efficiently find those lost points and keep the integrity of segmented leaf (Figure 8).
(c) T3 growth stages using regional-growing algorithm. Different colors represent individual leaves at different position on the plant. The red points on each leaf were discarded points at the first segmentation. Leaves in red circle were the curved leaves which lost many points after segmentation. Leaves in blue circle were overlapped with each other at stage T2 and cannot be segmented properly into individual leaves.

Accuracy Assessment
In order to estimate the accuracy of our method, the relationship between estimated plant height and measured height was compared and shown in Figure 9a. Good correlations were found between estimated and measured plant height (R 2 = 0.88, RMSE = 5.2 cm). However, most estimated heights were higher than measured values. The estimation accuracy of leaf blade length was also performed and shown in Figure 9b. The high R 2 (0.83) and low RMSE (1.77 cm) were found. The fitted line was close to the 1:1 line.

Accuracy Assessment
In order to estimate the accuracy of our method, the relationship between estimated plant height and measured height was compared and shown in Figure 9a. Good correlations were found between estimated and measured plant height (R 2 = 0.88, RMSE = 5.2 cm). However, most estimated heights were higher than measured values. The estimation accuracy of leaf blade length was also performed and shown in Figure 9b. The high R 2 (0.83) and low RMSE (1.77 cm) were found. The fitted line was close to the 1:1 line. Several points were far away from 1:1 line in Figure 9a,b, meaning that certain mismatches existed between estimated and measured values. To figure out these mismatches, two cases are shown in figure  10. In the first case, estimated height was much higher than measured height (Figure 10a,b). In the second case, estimated height is different from measured height (Figure 10c,d). Several points were far away from 1:1 line in Figure 9a,b, meaning that certain mismatches existed between estimated and measured values. To figure out these mismatches, two cases are shown in figure 10. In the first case, estimated height was much higher than measured height (Figure 10a,b). In the second case, estimated height is different from measured height (Figure 10c,d). Figure 11 showed the extraction results of leaf blade length for two kinds of leaves. Common leaves with upward facing posture (Figure 11a-d) can be fitted well on leaf blade midrib (Figure 11d) with our method. Therefore, the leaf blade length can be calculated well. The twisted leaves (Figure 11e-h) were generally located at the lower part of the plant, not exposing to the sun. They always twisted their petiole and changed their facing directions. Our method can obtain smooth curves for both kinds of leaves (Figure 11c,g). However, our method was not performed well for twisted leaves (Figure 11h).
The proportion of twisted leaves to total leaves was very small (three out of 70 leaves in our accuracy assessment). Therefore, our current method did not consider this situation. Several points were far away from 1:1 line in Figure 9a,b, meaning that certain mismatches existed between estimated and measured values. To figure out these mismatches, two cases are shown in figure  10. In the first case, estimated height was much higher than measured height (Figure 10a,b). In the second case, estimated height is different from measured height (Figure 10c,d).  Figure 11 showed the extraction results of leaf blade length for two kinds of leaves. Common leaves with upward facing posture (Figure 11a-d) can be fitted well on leaf blade midrib (Figure 11d) with our method. Therefore, the leaf blade length can be calculated well. The twisted leaves ( Figure  11e-h) were generally located at the lower part of the plant, not exposing to the sun. They always twisted their petiole and changed their facing directions. Our method can obtain smooth curves for both kinds of leaves (Figure 11c,g). However, our method was not performed well for twisted leaves (Figure 11h). The proportion of twisted leaves to total leaves was very small (three out of 70 leaves in our accuracy assessment). Therefore, our current method did not consider this situation.

Dynamic Estimation of Plant Phenotypic Traits and Biomass
The dynamic changes of plant height, maximum canopy area and convex hull volume were monitored from reconstructed point cloud ( Figure 12). The mean value of extracted plant height increased from stages T1 to T2, and peaked at stage T2. Then, plant height decreased slightly from stages T2 to T3. Mean value of maximum canopy area and convex hull volume increased gradually from T1 to T3. Phenotypic traits here varied greatly at stage T3 compared with those at stages T1 and T2. Five points deviated more compared with others regarding plant height, canopy area and convex hull. These were only a very small part compared with the number of experimental plants.

Dynamic Estimation of Plant Phenotypic Traits and Biomass
The dynamic changes of plant height, maximum canopy area and convex hull volume were monitored from reconstructed point cloud ( Figure 12). The mean value of extracted plant height increased from stages T1 to T2, and peaked at stage T2. Then, plant height decreased slightly from stages T2 to T3. Mean value of maximum canopy area and convex hull volume increased gradually from T1 to T3. Phenotypic traits here varied greatly at stage T3 compared with those at stages T1 and T2. Five points deviated more compared with others regarding plant height, canopy area and convex hull. These were only a very small part compared with the number of experimental plants.

Dynamic Estimation of Plant Phenotypic Traits and Biomass
The dynamic changes of plant height, maximum canopy area and convex hull volume were monitored from reconstructed point cloud (Figure 12). The mean value of extracted plant height increased from stages T1 to T2, and peaked at stage T2. Then, plant height decreased slightly from stages T2 to T3. Mean value of maximum canopy area and convex hull volume increased gradually from T1 to T3. Phenotypic traits here varied greatly at stage T3 compared with those at stages T1 and T2. Five points deviated more compared with others regarding plant height, canopy area and convex hull. These were only a very small part compared with the number of experimental plants.  The fresh weight prediction based on the extracted traits was investigated and shown in Figure 13. At stage T1 and T2, both leaf area and convex hull volume had close relationship with fresh weight with R 2 more than 0.80 (Figure 13a,c). The values of leaf area, convex hull volume and fresh weight at T2 were obviously higher than that at stage T1 as the growth of the plants. When incorporating calculated leaf area and convex hull at stage T3, the fitting results with biomass declined (Figure 13b,d). The values of R 2 fell to 0.82 and 0.78 for the relationship of fresh biomass with leaf area or with convex hull volume. The fresh weight prediction based on the extracted traits was investigated and shown in Figure 13. At stage T1 and T2, both leaf area and convex hull volume had close relationship with fresh weight with R 2 more than 0.80 (Figure 13a,c). The values of leaf area, convex hull volume and fresh weight at T2 were obviously higher than that at stage T1 as the growth of the plants. When incorporating calculated leaf area and convex hull at stage T3, the fitting results with biomass declined (Figure 13b,d). The values of R 2 fell to 0.82 and 0.78 for the relationship of fresh biomass with leaf area or with convex hull volume.

Processing of Point Cloud of Individual Plant
The segmentation plays a very important role in extraction of plant phenotyping but still remains a big challenge [40]. In this study, regional growing algorithm was used to segment plant into individual leaves initially. The algorithm, based on the surface smoothness constraint, was often adopted to segment plants with few flat leaves, such as cucumber, eggplant [41] and maize [42]. However, due to the complex environments in the field, some leaves of sugar beet were curved, which brought a big

Processing of Point Cloud of Individual Plant
The segmentation plays a very important role in extraction of plant phenotyping but still remains a big challenge [40]. In this study, regional growing algorithm was used to segment plant into individual leaves initially. The algorithm, based on the surface smoothness constraint, was often adopted to segment plants with few flat leaves, such as cucumber, eggplant [41] and maize [42]. However, due to the complex environments in the field, some leaves of sugar beet were curved, which brought a big challenge to the algorithm (Figure 8a). Clustering based on Euclidean distance solved this problem well, and made the segmented leaf to retain complete point cloud (Figure 8b). In addition, the segmentation of plant into individual leaves becomes difficult for the middle and later growth stages because of the overlapping between leaves [39]. In our study, less overlapping leaves were found at stages T1 and T3 compared with stage T2. At stage T2, the overlapping between leaves was so serious that even the manual separation cannot perform well ( Figure 7b). As far as we know, fully automatic segmentation has only been realized on broad-leaf plants with few leaves [14] and narrow-leaf crops with few leaves and without tillers [43].
The correction of point cloud of individual plant is essential for estimating plant phenotypic trait, such as plant height and canopy ground cover, etc. The success of correction largely relies on the detection of the ground. The RANSAC algorithm can detect ground and soil plane of pot well [30] in an indoor pot experiment. In our study, the RANSAC algorithm was utilized to detect ground plane ( Figure 2d). As the RANSAC algorithm always found the largest plane with the most point cloud as the ground, another smooth plane was detected as base ground due to the uneven ground for field environment as shown in Figure 10c. With the wrong detection of the ground, plant phenotypic traits cannot be estimated well. To avoid this deficiency, the ground point cloud should be preserved as much as possible when the plant population is divided into individual plants (Figure 2b,c).

Extraction of Plant Phenotypic Traits
Plant height, as a basic phenotypic trait, can be used not only as an intuitive indicator of plant growth state, but also for the estimation of plant biomass and yield [44,45]. In this study, the ground truth of plant height was measured as the vertical distance between the bases of the stem to the highest point of the plant. Estimated plant height was calculated from bounding box of point cloud of individual plant after coordinates correction and filtering (Figure 2g). The estimated plant height was highly correlated with the measured height (R 2 = 0.88), but was a little bit higher than the measured value ( Figure 9a). It was found that the field ground was uneven, many sugar beet leaves fell on the ground which was lower than base of the stem (Figure 10a). Therefore, the height of the bounding box was often higher than the measured plant height (Figure 10b). Another obvious differences between the measured and calculated plant height were the inaccurate coordinate correction of the plant (Figure 10c,d).
Leaf blade length has been estimated in the literature based on meshes of leave blades [14,26,41]. It was very convenient and accurate to estimate the leaf blade length based on the point cloud of leaf blade. The most important step was to find the base and tip of the leaf blade accurately to calculate leaf blade length. Principal component analysis (PCA) was often adopted to find endpoints of individual leaves [43], but its accuracy was greatly affected by the shape of the leaves. Due to the complex shape of sugar beet grown in the field, we considered the relative position of leaves in the plant as a criterion to find the endpoints of leaves. Researchers have used the shortest path algorithm to find the midrib of leaf blade based on the endpoints on the leaf blade [43], while we rotated the Z-axis so that the endpoints of leaf were on the same axis to find the midrib of leaf blade. Our method can perform well for most of the leaves (Figure 11d). But certain deviation occurred between extracted curve and midrib of leaf blade for twisted leaves, (Figure 11h). Though the twisted leaves occupied less and was always located at the lower part of the plant, a more comprehensive method was still needed to estimate the length of all types of leaves.

Estimation of Biomass Using Plant Phenotypic Traits
Biomass estimation is important for biological research and agricultural management [46]. Many studies have shown that the total leaf area was highly consistent with biomass [47][48][49]. The correlation between total leaf area and fresh biomass was figured separately in this study with stages T1 to T3 (Figure 12b) and stages T1 to T2 (Figure 12a), and a better performance was found for the latter (R 2 = 0.88) compared with the former (R 2 = 0.82). Sugar beet grows very fast from stage T1 to T2 and is gradually aged from stage T2 to T3. There were other phenotypic traits involving plant area, such as leaf project area and surface area, that have been used to estimate biomass [29]. The correlations between these phenotypic traits and biomass may be lower than between total leaf area and biomass, but it is easier to measure. In addition, researchers have found that stem volume has more advantages than total leaf area in estimating biomass at the later stage of plant growth [30].
The convex hull and the concave hull are used to describe the volume [36,50]. As researches have shown that the concave hull has poor performance for biomass estimation [29], we only studied the correlation between the convex hull and biomass. Similar to the total leaf area, the estimation of stages T1-T2 is better than that of T1-T3, but the difference is not significant. The biomass estimation with convex hull is obviously inferior to that of total leaf area ( Figure 12). Lati [51] used the convex hull to estimate biomass of Black Nightshade. The R 2 is much higher than ours (0.88 versus 0.80). Our sugar beets were scattered in stages T1 and T3, but concentrated in stage T2. Only total leaf area and convex hull were used to explore the relationship with fresh weight in the study. More phenotypic traits will be used to predict the fresh and dry weight in the future.
Phenotyping studies of individual plants are usually carried out indoors, but rarely in the field because of the complex environment in the field. The SFM method can get the accurate point cloud of individual plant in the field, which is convenient to extract the organ-level traits accurately. Terrestrial light detection and ranging (lidar) can perform well for the phenotyping of individual plant in the field [52], but very expensive, with more than 100 times the price than common camera. Unmanned Aerial Vehicle (UAV) phenotyping has developed rapidly in recent years. UAV mounted with an rgb camera, lidar or multispectral camera, were always used to obtain point cloud or spectral information more quickly. However, UAV based phenotyping focuses on canopy scale with average height of plants, leaf area index and lodging area. For example, Hu [53] used UAV with an rgb camera to obtain DSM (Digital Surface Model) of sorghum in the field, and estimated the average height of sorghum. Su [54] used UAV with an rgb camera and multispectral camera to estimate average height, canopy leaf area index and lodging area of corn plants. Lei [55] used UAV with lidar to inverse leaf area index of maize and studied the effect of leaf occlusion on inversing leaf area index of maize. As far as we know, UAV with sensors cannot estimate individual organs phenotypic traits very well.
Our research also has limitations. Due to the influence of wind, the shake of plants will bring errors to the reconstruction of 3D model. Zhu [56] used a black and white chessboard pattern of cloth and iron frame as windshield equipment to improve the quality of 3D reconstruction. The influence of wind is very common in the field research. Minimizing the time interval of each data acquisition may be a feasible solution. In addition, the SFM method is convenient for its data collection but time-consuming for its post-process. With the development of GPU computing, the post-processing time has been greatly reduced, but still cannot meet the needs of real-time reconstruction. As the simultaneous localization and mapping (SLAM) is developing rapidly, a more efficient reconstruction algorithm is to be expected in the future. In the next study, more genotypes and phenotypic traits of sugar beet will be extracted for linking genetic with phenotypic traits and provide a basis for sugar beet breeding.

Conclusions
In this study, the 3D models of sugar beets for 20 genotypes were reconstructed by SFM method at three growth stages in the field. We developed an automatic pipeline for data processing, including processing of point cloud of individual plant and extraction of plant traits. Segmenting individual plant from plant population manually was the first step of processing of point cloud of individual plant. Then, the soil plane was detected using the RANSAC algorithm, as a criterion to correct the coordinates of point cloud of individual plant. We selected G minus R (G−R) as a color filter to separate plant from soil and shadows and used statistical filter to remove outliers. Finally, regional-growing algorithm was utilized to segment plant into individual leaves and Euclidean distance clustering was used to improve quality of segmented leaves. Extraction of plant traits is other important part of our pipeline. Global traits, such as plant height, maximum canopy area and convex hull volume, were extracted based on the overall structure of plant. Organ-level traits, such as leaf blade length and total leave areas, were extracted based on individual leaves. Plant traits were also used to explore their relationship with sugar beet biomass. Automated measurements from the pipeline showed agreement with the corresponding manual measurements, validating the accuracy and utility of our pipeline. Both the total leaf area and convex hull volume showed potential for biomass prediction (R 2 = 0.78-0.88). The pipeline proposed in this study can perform well for estimating plant phenotyping in the field and provide basis for breeding programs.