Non-Destructive Measurement of Three-Dimensional Plants Based on Point Cloud

In agriculture, information about the spatial distribution of plant growth is valuable for applications. Quantitative study of the characteristics of plants plays an important role in the plants’ growth and development research, and non-destructive measurement of the height of plants based on machine vision technology is one of the difficulties. We propose a methodology for three-dimensional reconstruction under growing plants by Kinect v2.0 and explored the measure growth parameters based on three-dimensional (3D) point cloud in this paper. The strategy includes three steps—firstly, preprocessing 3D point cloud data, completing the 3D plant registration through point cloud outlier filtering and surface smooth method; secondly, using the locally convex connected patches method to segment the leaves and stem from the plant model; extracting the feature boundary points from the leaf point cloud, and using the contour extraction algorithm to get the feature boundary lines; finally, calculating the length, width of the leaf by Euclidean distance, and the area of the leaf by surface integral method, measuring the height of plant using the vertical distance technology. The results show that the automatic extraction scheme of plant information is effective and the measurement accuracy meets the need of measurement standard. The established 3D plant model is the key to study the whole plant information, which reduces the inaccuracy of occlusion to the description of leaf shape and conducive to the study of the real plant growth status.


Introduction
With the development of protected vegetable production, vegetable seedlings are increasingly recognized by producers, which has become an important pillar of vegetable industry development.
In agriculture, information about the spatial distribution of crop growth is valuable for applications such as biomass and yield estimation or increasing field work efficiency in terms of fertilizing, applying pesticides and irrigation [1,2]. Many researchers have carried out experiments on plant growth detection based on computer vision technology, mainly covering three aspects: the detection of external growth parameters, fruit maturity, and the detection of nutritional components [3].
Plant height is an important phenotypic morphology parameters that can be used not only as an indicator of overall plant growth vigor, but a parameter to estimate traits [4,5]. The key technology to obtain plant height is to present the highest and the lowest points on the image (or point cloud), and obtain the height by distance formula [6][7][8]. Using the Euclidean distance [9,10] to present the plant height is commonly used on 2D images, extracting the plant skeleton on the region of interest (ROI) and detecting the lowest and highest points of the skeleton. The accuracy of pixel extraction in the interest

Preprocessing of a Three-Dimensional Plant Model
In this experiment, a Kinect camera and an electric turntable were used to build a 3D plant model according to the following instructions. Fix the camera position and put the plant on the turntable, the height h of Kinect is set to 0.7 m, and the distance d between Kinect and plant is 1.0 m. Figure 2 shows the experiment devices which are placed at fixed positions, and marks the world coordinate system of the point cloud. Place the plant on the turntable and record the initial point cloud of the plant. Turn the turntable clockwise 45 • , then stop the turntable and record the data as a piece of point cloud. Stop running until the turntable rotates 360 • and 8 piece of point clouds are obtained. The original point cloud contains invalid data, so in this experiment the outlier points should be removed by statistical outlier filter [44], the distance d i from the query point to all the neighboring points is calculated, and the threshold standard deviation ε is set to judge whether the point is an outlier. Eight pieces of point clouds are registered by iterative closest point (ICP) [45] algorithm to obtain the initial three-dimensional point cloud model. In Figure 3, eight pieces of point clouds are obtained from different directions, related to the rotation angle of the turntable, angle 0 means the turntable rotates 0 • , angle 1 means rotating 45 • , angle 2 means rotating 90 • , angle 3 means rotating 135 • , angle 4 means rotating 180 • , angle 5 means rotating 225 • , angle 6 means rotating 270 • , angle 7 means rotating 315 • , angle n is recorded every 45 • (n ∈ [0, 7]) [46]. The point cloud of the turntable and the surrounding scene is marked as p b , the point cloud of ROI from each angle can be obtained by subtracting p b . Figure 3 shows the plant's point clouds obtained after rotating the turntable 8 times interval, and the three-dimensional points of ROI are color-rendered according to x-axis for displaying clearly. The reason for selecting eight directions to collect experimental data is that the two adjacent point clouds can fill the holes at the edge of the point cloud after registering, so the three-dimensional point cloud can reduce the influence of holes on accuracy.
Plant reconstruction is divided into two parts in this research. The first part is to register angle 0 with 7, 1, 2 and angle 4 with 3, 5, 6 separately. The former result is called Pos 0 and the latter result is called Pos 1 . The second part is to stitch the Pos 0 and Pos 1 point clouds to get a complete 360 • reconstruction model. Here, the commonly used stitching method is the ICP algorithm, and this method has many advantages. Not only the registration between point set is considered, but also the registration from point set to model and model to model is also considered in the ICP algorithm. The basic principle of the ICP algorithm is-in the target point cloud P (i) and the source point cloud Q (i) , find the nearest neighbor point (p i , q i ) according to certain constraints, and then calculate the optimal matching parameters R and t to make the error function minimum. The error function is E(R, t) in Equation (1) [45], where n is the number of nearest neighbor pairs, p i is the point in the target point cloud P, q i is the nearest point of the source point cloud Q, R is the rotation matrix, t is the translation vector. The traditional ICP algorithm can be summarized in two steps-(1) compute correspondences between the two point clouds.
(2) calculate the transformation that minimizes the distance between the corresponding points. The disadvantages [47] of this method are as follows-noise or abnormal data may cause the algorithm to fail to converge; the selection of the initial value has an important impact on the final registration results. The algorithm focuses on the similarity between the target and the source, and takes the minimum error as the evaluation standard, ignores whether there is a location relationship between the target point cloud and the source point cloud.
The improved ICP algorithm in this paper is more suitable for registering point clouds from different directions using a turntable. This stitching method first uses a rotation matrix and a translation vector to process the relationship between adjacent point clouds, then uses the ICP algorithm to iteratively obtain R, t. Obviously, the point clouds in eight directions obtained in this paper have positional relations. The point clouds in our experiment are formed by multiple rotations from the same plant, plant leaves from different directions are similar in shape, which can cause the ICP algorithm falling into local optimality, so this algorithm is improved. According to the position relationship of adjacent patch point clouds, firstly, rotate the source point cloud at angle θ 1 , then perform ICP registration of each point cloud, the following formula is in Equation (2), where (x t , y t , z t ) is the centroid of the target point cloud, (x s , y s , z s ) is the centroid of the source point cloud, ω 1 is the rotation matrix around the y-axis for angle θ 1 , and the value of ω 1 is based on the rotation relationship of target point cloud and source point cloud. After obtaining the centroid point three-dimensional vector of 8 pieces of point clouds, calculate the distance σ 1 between target point cloud and source point cloud, take target point cloud as the registration center, move source point cloud to the centriod coordinate point according to the distance σ 1 .
In this experiment, take angle 0 as the target point, rotate angle 7,1,2 to register the two point clouds respectively; take angle 4 as the target point, rotate angle 3,5,6 to register the two point clouds respectively. For example, take angle 1 as the source point cloud, and angle 0 as the target point cloud in Figure 4a, color green represents angle 0 , red represents angle 1 , the relationship between them is rotating angle 1 of 45 • anticlockwise along the angle 0 , so θ 1 is 45 • . The iterations number of matching target and source point cloud n is 1000 times. The traditional ICP algorithm cannot judge the relationship between source point cloud and target point cloud as in Figure 4b, it shows that the stem and leaves do not match correctly and the position relationships of point clouds are neglected in the matching process, so the stitching error is large. The location relationship between the source and target point cloud should be determined in this research, that is, applying the ICP algorithm after rotation θ 1 , and the results are shown in Figure 4c by our method, after 1000 iterations, the leaves and stem match well. According to the rotation direction and angle of the turntable, the problem of ignoring the relationship between target point cloud and original point cloud in ICP registration can be solved by calculating the rotation angle in advance when registering the adjacent point clouds. Similarly, after obtaining Pos 0 and Pos 1 point cloud, the ICP method cannot be used to match them directly, because Pos 0 and Pos 1 are in the different 3D coordinate systems, it should rotate Pos 1 by 180 • and keep it in the same coordinate system with Pos 0 .
For the Pos 0 and Pos 1 point clouds in Figure 5a, color red represents Pos 0 , green represents Pos 1 , half of the 360 • model has been registered separately, so the correct registering of them directly affected the accuracy of the whole model. Because of the randomness in the growth of leaves, the centroids of Pos 0 and Pos 1 could not overlap. If ICP registration is directly used for registering, it will fall into local optimal and mismatching, as shown in the Figure 5b, which cannot reflect the relationship between Pos 0 and Pos 1 point cloud correctly. For this situation, our approach is to extract the difference between Pos 0 and Pos 1 centroid, move the target point cloud to the source point cloud, that is the local movement between these two centroid vectors and then using the ICP method to obtain a complete plant point cloud model as shown in Figure 5c. The algorithm cannot fall into the local optimum, and the registration effect of Pos 0 and Pos 1 is better. Here, σ 2 is the difference between Pos 0 and Pos 1 's centroid vector, as in Equation (3): where (x pos0 , y pos0 , z pos0 ) is the centroid of Pos 0 , (x pos1 , y pos1 , z pos1 ) is the centroid of Pos 1 , ω 2 is the rotation matrix around the y-axis for angle θ 2 , and θ 2 is 180 • . After obtaining the centroid point three-dimensional vector of Pos 0 and Pos 1 , calculate the distance σ 2 of them, take Pos 0 as the registration center, move Pos 1 to the centriod coordinate point with Pos 0 according to the distance σ 2 . At this time, the registration and establishment of the 3D model is completed. To create complete models, glossy surfaces and occlusions in the point cloud must be accounted for. A solution is to use a resampling algorithm, which attempts to recreate the missing parts of the surface through high-order polynomial interpolation between surrounding data points. By performing the resampling process, small errors can be corrected, and the double-wall artifacts generated by registering 8 angles point cloud data together can be smoothed. Moving least squares (MLS) surface reconstruction [48] method can estimate normal vector based on polynomial reconstruction, and can also smooth and resample noisy data. In order to achieve the smoothness of the surface, the K nearest-neighbor radius [49] of the fitting is set as 0.5 mm, and the whole 360 • point cloud is smoothed. As shown in Figure 6 Through the above steps, a whole point cloud model with less data, accurate normal and curvature variance is obtained, which is beneficial to the following operations such as feature point extraction and feature boundary line collection.

Extraction of Boundary Points on a Plant Leaf
In this section, we discuss how to segment leaves and stems on a complete plant model, extract the contour of each leaf, and calculate the phenotype of the leaf according to the three-dimensional points of the contour. In this paper, kdtree is used to organize the data of leaf point cloud and realize fast nearest neighbor retrieval based on fast library for approximate nearest neighbors (FLANN) [50]. The spatial topological relationship between data points is established to facilitate k-nearest neighbor search.
The pepper plants are used as the research objects to generate three-dimensional models by eight pieces of point cloud. On the basis of the 3D model, the stem and leaves should be divided into different parts, so the phenotypic characteristics of each leaf can be calculated. The leaves for measurement are extracted from pepper plants by the locally convex connected patches method (LCCP) [51]. The principle of LCCP method is, for the plant model, firstly to calculate the convexity-concavity relationship of adjacent patches. The relationship is judged by extended connectivity criterion and sanity criterion. Extended connectivity criterion uses the angle between the center line vector x 1 , x 2 and the normal vector n 1 , n 2 of the adjacent patches. The vector from n 1 to n 2 is set to t, the angle between t and n 1 is a 1 , and the angle between t and n 2 is a 2 . Obviously, if a 1 > a 2 , the relationship of the two patches p i , p j is concave, otherwise is convex. d is the difference between the vector x 1 , x 2 , s is the cross product between the vector n 1 , n 2 . If the relationship of the two patches is convex , it should be further judged by sanity criterion. When the angle between d and s is greater than the threshold δ, it can be sure that the relationship between them is convex, otherwise is concave.
After marking the concavity-convexity relationship of each adjacent small region, the region growth algorithm is used to cluster the small regions in Figure 7a into larger objects in Figure 7b, representing by small blocks of different colors. This algorithm restricted by the convexity of small regions can clearly distinguish the boundary of stem and leaves as shown in Figure 7b, that red line represents concavity and blue line represents convexity between adjacent patches. There is obvious concavity and convexity change at the junction of fruit stem and leaves, which meets the criteria of LCCP segmentation.
The convexity-concavity relationship of the plant point cloud is obtained by LCCP algorithm, which displays in different colors. From the effect of segmentation, the closer the leaves are to the bottom of the plant, as shown in Figure 7b, the relationship between leaves and stem are more obvious; while the closer the leaves are to the top of the plant, for example, the relationship between leaf 3 and leaf 4 is not obvious, so they could not be segmented correctly by the LCCP method. Because the leaves at the top of the plant are closer towards the stem, they are treated as a whole, so it is necessary to manually segment the leaves. For each divided leaf, the boundary condition [52] is used to determine whether the three-dimensional point of the leaf is an internal point or external contour point. For any point p, a tangent plane is established according to the point and its k neighborhood points, and project each point on this plane; then, judge whether the point p is a boundary feature point, repeat the above operations until all points are judged, and the boundary feature point set is obtained. When the target point is a boundary point, only a few or no points will appear in its upper region. According to this phenomenon, the discrimination mechanism of boundary feature points can be established. As shown in Figure 8, the upper part of point P is divided into two areas I and II by the vertical line and the horizontal line passing through point p. p is a boundary feature point can be distinguished according to whether there are data points in these two areas. Specifically, according to the boundary determination method, the edge point is the key to get the phenotype of leaf, taking k = 10 as an example, it can be divided into three situations: (1) there are no points in the either areas, p is the boundary point; (2) there are no points in one area, p is the boundary point; (3) there are points in both areas, find out the vector closest to the vertical line of these two areas, and the angle between them. If it is greater than the set threshold , p is the boundary point, threshold is π 2 in this paper. In all three cases, the points are boundary points and need to be preserved. According to this method, edge extraction is performed on the segmented 3D leaves to extract points. The boundary points of the divided leaf are the basis of calculating the length and width of each leaf. Taking the green pepper plant at seeding stage as the research object, the LCCP algorithm is used to segment the leaves and stem, the boundary point judgment method is used to get the internal point and the external point, and then the edge points of each leaf are obtained, the leaves and stem are displayed in different colors. There are five leaves, and the edge points are displayed in black in Figure 9, which are the key points to obtain the leaf's phenotypic parameters.

Calculation the Length, Width and Surface Area of Leaf
Calculate the length and width according to the contour of each leaf. On the three-dimensional points of the boundary, manually select p top at the top of the tip and p bottom at the bottom of the leaf to calculate the Euclidean distance between them as the length of the leaf; and select the maximum distance perpendicular to leaf length on the edge point as p side1 and p side2 respectively to calculate the Euclidean distance between them as the width of the leaf in Equation (4), where (x t , y t , z t ) is the coordinate of p top , (x b , y b , z b ) is the coordinate of p bottom , (x 1 , y 1 , z 1 ) is the coordinate of p side1 , and (x 2 , y 2 , z 2 ) is the coordinate of p side2 . l lea f is the length of leaf, w lea f is the width of leaf. In Figure 10, the blue three-dimensional points are the leaf contour obtained by boundary conditions, it is shown that the length and width of leaf can be obtained by using the Euclidean distance method, where length refers to the three-dimensional distance between the top and the bottom of the leaf, and width refers to the three-dimensional distance between the midpoint on both sides of the leaf. In the 3D plant model, (x, y) of the point cloud coordinates (x, y, z) can be regarded as a coordinate pair obtained by sampling the coordinates on the xOy plane. For grid sample points, it should generate two matrices of the same size with the vector x as the rows and the vector y as the columns, where the rows of x start from the minimum value x min to the maximum value x max of the model, and generate data every 0.1 mm, which integrate as Matrix XI; Similarly, the columns of y start from the minimum value y min to the maximum value y max , and generate data every 0.1 mm, which integrate as matrix YI. XI and YI form a uniform grid. XI is a row vector that determines a matrix with a fixed number of columns, YI is a column vector that determines a matrix with a fixed number of rows. Fit the surface (XI, YI, ZI) formed by z = f (x, y) to the scattered data in the vector (x, y, z). Perform cubic interpolation of data on the surface at the query point specified by (XI, YI) and return the interpolated value ZI to generate a smooth surface. The surface passes through the data points defined by (x, y) to form a complete mesh surface.
The leaf area is calculated according to the three-dimensional points. For the reconstructed leaf surface, the area is calculated by surface integral algorithm, and the area of the small rectangular block is added to get the surface area. When calculating the gradient, the step size needs to correspond to the small rectangular block's length d x and width d y . d x and d y in the Equation (5) are the slopes of rectangular blocks in x-axis and y-axis direction respectively. The area obtained from the integration of surface s on region D xy is: The fitted surface of leaf is as shown in Figure 11, the leaf is a closed surface, and the value of d x and d y are set at 0.0001, that means, the area of small square shall be obtained in every step of 0.1 mm. Figure 11. Fitting surface of leaf with double integral algorithm.

Measurement the Height of a Plant
The height of the plant is obtained by vertical distance along y direction, by comparing all points' three-dimensional coordinates, the maximum value y i max and minimum value y i min value along the y-directions are calculated in Equation (6). The height of the plant H is presented with flowerpot, H 1 is the vertical distance of the whole point cloud. Because the height of the flowerpot H 2 is included in the height H 1 , in order to get the true height of the plant, the actual height of the flowerpot should be subtracted here. In Equation (7), the height of the cuboid H 1 is the sum of the plant H and the flowerpot H 2 , and the height of the flowerpot has been subtracted in this experiment.
The height of the flowerpot H 2 can be obtained in Figure 7b after the result of LCCP segmentation. The concavity-convexity relationship between the flowerpot and the stem of the plant separate the two into different spaces. Therefore, the maximum value y maxpot and the minimum value y minpot of flowerpot in the y-axis direction can be used as the vertical distance of the flowerpot.
In Figure 12, the red arrow represents the x-direction, the blue arrow represents the y-direction, and the green arrow represents the z-direction, the external minimum cuboid is drawn by the main direction coordinate system. Then, for the y-axis of interest in this paper, the actual height difference between the cuboid's vertical distance and the flowerpot is obtained to get the height of the plant.

Experiment
The experimental results are analyzed by comparing the measured values with actual values. Plant height, leaf length, width and surface area are the main analysis indexes. According to the experiment in this paper, the absolute error (AE) of plant height, the mean absolute error (MAE), root mean square error (RMSE) of leaf length, leaf width and surface area of each plant are calculated. Ten seedling plants of pepper are selected as the test objects, denoted as plant n , (n ∈ [1,10]).
Based on Kinect V2, the 3D model of green pepper was constructed and the plant's phenotypic parameters including plant's height, leaf's length, width and area were measured. The AE of plant height is the difference between the actual value and the measured value, and the MAE of plant height is the average value of 10 plants' AE. The minimum actual value of plant height of 10 green peppers is 11.12 cm, and the maximum actual value is 21.09 cm. In Figure 13, the maximum AE of plant height is 0.66 cm, the minimum AE is 0.15 cm, and the MAE of plant height is 0.392 cm. The leaf phenotype is an important growth parameter for pepper seedling plants. In this paper, pepper plants' point clouds were segmented by the LCCP method, four leaves of each green pepper were selected as a whole for measurement, and recorded these as measured values. The average length, width, area errors of four leaves for each green pepper are denoted by MAE l , MAE w and MAE a . The maximum actual value of MAE l is 0.365 cm, and the minimum actual value of MAE l is 0.1825 cm; the maximum actual value of MAE w is 0.305 cm, and the minimum actual value of MAE w is 0.2175 cm; the maximum actual value of MAE a is 1.1125 cm 2 , and the minimum actual value of MAE a is 0.821 cm 2 . The MAE of the leaf's length, width and area of these 10 peppers are denoted as MAE lea f s length , MAE lea f s width and MAE lea f s area . In Figure 14, MAE lea f s length is 0.2537 cm, MAE lea f s width is 0.2676 cm and MAE lea f s area is 0.956 cm 2 . MAE l and MAE w are directly proportional to MAE a , because the leaf's area is related to the value of length and width, with the increase of MAE l and MAE w , the MAE a become larger. In Table 1 the measurement accuracy in this paper was compared with other references. In these methods, Reference [24,25] use binocular vision to obtain the point cloud, Reference [28] uses the 3D laser scanner to obtain the point cloud, and References [3,15] use Kinect to obtain point cloud of regions of interest (ROIs). In Reference [25], VisualSFM was used to register the point cloud of eggplant, pepper and cucumber, segmented these plants by region growing method, and measured the leave's parameters of these plants. In Reference [24] 3DSOM was used to register the point cloud of Gossypium hirsutum, segmented the plant by region growing and mesh method, and measured the leaves' parameters of the plant. In Reference [28], target ball extraction algorithm [44] was used to register the point cloud of apple tree branches, segmented the plant by LCCP and Kmeans clustering algorithm, and measure the leaves' parameters of the plant. In Reference [3], higher solution RGB images was used to register the point cloud of pepper, measured the leaves' parameters and plant's height. In Reference [15] segmented the cucumber plants using the Euclidean clustering method, and using the vertical instance method to obtain the height of the plants. Compared with other methods, MAE and RMSE of our method have lower comprehensive errors, and our method has advantages of leaf's surface area measurement. The measurement errors of plant height, leaf length and width are also at the same level with other measurement methods.

Conclusions
The research objects of this paper are the plants in the early growing period. Taking green pepper plants as an example, each leaf plays an important role in the growth of plants. The information of the sum of single leaf on the whole plant is useful for assessing the growth state. The established three-dimensional model is the key to study the whole plant information. Using depth camera and turntable to complete point cloud acquisition and registration is the basic operation for plant segmentation, height measurement and a leaf's phenotypic features extraction. This paper processes the point clouds according to the characteristics of rotation by equal intervals, which can avoid the process of registration falling into the local optimal solution and ignoring the overall information. The plant model obtained by rotation and registration can reduce the inaccuracy of leaf shape description due to different shooting angles.
The use of the LCCP algorithm to segment the plant point cloud model is a popular three-dimensional point segmentation method. This method avoids the process of selecting initial points of region growing algorithm. The LCCP algorithm uses concave-convex features to segment points, while there is a natural convexity-concavity relationship between leaves and stems. This is another reason why we choose this segmentation algorithm to segment the plant's point cloud but the convexity-concavity relationship of the adjacent leaves' point clouds is not obvious, sometimes the young leaves on the top of the plant need to be selected manually. For plants in the growing period, the Leaf Area Index has a greater influence on the growth state of the plant. At this time, the leaves of the plant are dense and the shape of each leaf is larger. Obtaining each leaf individually is impossible, because there is often occlusion between dense leaves, which cannot be reduced by changing the shooting angle. Therefore, it is necessary to find new phenotypic features that can represent plant growth instead of leaves' parameters.
Our experiment was done indoors, and the shooting camera and the object were placed horizontally, so the three-dimensional points of y-axis coordinate system was used to determine the plant height. At present we have not reproduced the experiment outdoors. Because the outdoor scene is complicated as the plants are generally in the ground, it is difficult to complete the 3D modeling of plants according to the proposed rotation method. However, this paper proposed a general method, which can be extended to the related fields of 3D model building indoors.

Conflicts of Interest:
The authors also thank the editor and anonymous reviewers for providing helpful suggestions for improving the quality of this manuscript.

Abbreviations
The following abbreviations are used in this manuscript: