Rapid estimation of apple phenotypic parameters based on 3D reconstruction

In order to obtain the phenotypic parameters of apple quickly and accurately, which were commonly used as the basis of fruit sorting, a fast estimation method of apple phenotypic parameters based on three-dimensional (3D) reconstruction was proposed in this study. In this study, a three-dimensional model was constructed to estimate the phenotypic parameters of apple, such as volume, height, diameter, and fruit shape index. Firstly, an image acquisition system was built to capture sequence images of fruit with a binocular stereo vision system, and the images were extracted and matched using the Accelerated-KAZE algorithm to create the point cloud data. Secondly, the point cloud data were matched with the algorithm of Iterative Closest Point to establish a whole model of apple, and the surface reconstruction model of fruit was obtained by constructing irregular triangulation network. Finally, the apple phenotypic parameters were calculated by means of segmentation, surface complement and integral of the fruit model. Total of 200 apples were used as samples in the experiment. By this method, the phenotypic parameters of the apples were estimated based on their 3D reconstruction model, and the linear regression analysis was carried out between the estimated values and the real values. The results showed that R of the linear regression fitting of each parameter was higher than 0.90. Among them, the fitting of volume was the best with R of 0.97. In addition, the average errors of apple volume, height, fruit shape index, maximum diameter D and minimum diameter d were 8.73 cm, 1.43 mm, 1.28%, 0.90 mm, and 1.23 mm, respectively. According to the Chinese national standard of “fresh Apple”, the average error of the estimated result is within the range of allowable error. It indicates that the method of apple phenotypic parameter estimation based on 3D reconstruction has a high accuracy and practicability, and it can be used as the support for fruit sorting.


Introduction 
China is the largest producer of fruits and vegetables in the world, with an output of about 274 million t, which drives related industries to achieve an output value of more than 1 trillion CNY [1,2] . However, for a long time, China is in the status of emphasizing production and neglecting the processing. The postharvest processing rate is about 30% in China, which is less than 70% of that in developed countries, resulting in a 20% loss of postharvest yield due to decay [3,4] . There are two main reasons for this situation. On the one hand, there is insufficient understanding of intelligent quality classification of high-quality fruits and vegetables, low added value of commodities and poor market competitiveness. On the other hand, there is a lack of technical index system and standards for postharvest production, and serious lag in the research and development of related domestic equipment [5,6] . Therefore, there is an urgent need to strengthen the research on the technology of postharvest sorting and processing of fruits and vegetables.
The phenotypic parameters can objectively describe the appearance of fruit, which is one of the important bases in the fruit sorting process. In the past, the contact phenotypic parameter measurement had many problems such as low efficiency, low precision, easy damage and so on, and was gradually replaced by non-contact measurement.
The traditional researches of non-contact measurement were mostly carried out based on the method of two-dimensional image processing [7,8] . Among them, they mostly focused on apple [8,10] , kiwifruit [11,12] , navel orange [13,14] and so on. For example, Huang et al. [15] proposed an online apple classification method based on image feature fusion, which used the methods of discriminant tree and support vector machine to make fusion decision and achieve accurate classification. Fu et al. [16] studied the shape grade of kiwifruit and used the method of stepwise multiple linear regression to select the significant variables to predict the minimum diameter and volume of the equatorial section, and established the corresponding estimation model. Javadikia et al. [17] extracted 14 parameters from three varieties of oranges respectively, and estimated the quality of different varieties of oranges by using adaptive network fuzzy inference system (ANFIS) method. Wang et al. [18] adopted the iterative least median square ellipse fitting method to detect and classify the size of grapes on line. The above researches can detect and sort the fruits and vegetables by a single two-dimensional image.
Although they can realize industrialization and reduce labor intensity, there is a large classification error due to lack of information in the extraction of phenotypic parameters of fruits and vegetables in two-dimensional images, which is difficult to meet the requirements of practical use.
In recent years, the rise of 3D data acquisition devices such as stereoscopic camera [19][20][21] , depth camera [22,24] and LIDAR (Light Detection and Ranging) [25,27] has provided a new method for 3D reconstruction or 3D measurement of fruits and vegetables [28,30] . For example, Qiao et al. [31] detected the posture information of beans using binocular stereo cameras, and then determined the actual distance between beans. Liu et al. [32] established a laser speckle image acquisition system using laser speckle technology to distinguish the defects, pedicel and calyx of pear. Liu et al. [33] took strawberry plants as the research object and obtained the three-dimensional structure of strawberry canopy using the method of depth information segmentation and clustering, to further analyze the spatial structure of branches and leaves. Yao et al. [34] collected point cloud data of rape plants with a portable laser scanner to realize non-destructive measurement of rape leaf area by constructing a three-dimensional model by generating Delaunay triangulation. From the above researches, the extraction of phenotypic parameters of fruits and vegetables can achieve a good result based on the methods. However, the acquisition method or processing process of the target image, currently in common use, cannot well solve the problems in fruit and vegetable sorting industry. For instance: (1) only using a stereo camera to collect image, it cannot obtain all apparent information of the target, which will make the sorting result unreliable; (2) with depth camera, the image capture speed cannot meet the actual requirements of industrialization; (3) by a LIDAR, only the local contour of the target can be obtained, and the information such as color and texture is missing. Therefore, to meet the demand of fast, accurate and lossless acquisition of apple phenotypic parameters such as volume and shape in domestic fruit and vegetable sorting industry, a method for calculating multiple phenotypic parameters of apple based on complete three-dimensional model was proposed in this study. Firstly, a binocular stereo vision system was designed with an apple transport chain constituted by the dumbbell roller as the main body. Then, the sequence images of the apple were collected in the rotating state, and the three-dimensional model of the apple was constructed by using the method of three-dimensional reconstruction. Finally, the phenotypic parameters of apple such as volume, height, diameter, and fruit shape index were calculated by segmenting, filling surface and integrating the model. With the proposed method, it can achieve an accurate and reliable sorting result with an improved sorting speed.

Measurement of phenotypic parameters of apple
In this experiment, Red Fuji apple was used as the research object. In order to enrich the sample data and increase the reliability of the estimation results, 200 apple samples with a wide range of parameters such as size and shape were selected. A total of 5 phenotypic parameters of apple were measured in this study, including volume, height, transverse maximum diameter D, transverse minimum diameter d, and the fruit shape index β. Among them, the true values of height and diameter D and d were measured directly using Vernier calipers, with a resolution of 0.01mm. The fruit shape index β was calculated indirectly from the diameter. The true values of the volume were measured by the method of drainage, and two beakers of 1000 mL capacity and two standard cylinders of 250 mL capacity were used as measurement tools, with a resolution of 1 mL.

Imaging platform and data acquisition
In the experiment, the imaging platform was constituted by a notebook computer and an image acquisition system. The image acquisition system was mainly composed of two color industrial cameras, an apple conveyor chain, two direct current motors, a roller belt, a light source and a control system, as shown in Figure 1. The cameras are MV-CA00-11UC (Hikrobot, Hangzhou, Zhejiang, China). The light source is composed of 4 shadowless lamps with a size of 120×180 mm, and the apple conveyor chain is composed of several pairs of dumbbell rollers. The internal structure of the imaging system is shown in Figure 2.
When the image acquisition system is working, the apple conveyor chain driven by the DC motor moves forward at a uniform speed. It will contact with the below belt moving reverse when the dumbbell roller carrying the apple arrives under the camera, which will drive the apple rotating in a direction. Under the camera, the largest apple will rotate a circle within 1.5 s, while the camera continuously captures 32 images and transmits them to the laptop for real-time processing. When the largest apple rotates for one revolution within 1.5 s, the smaller apple may rotate for a revolution and a half. The captured apple sequence images of the smaller apple include the information of a revolution and a half, resulting in more information overlap and a waste of calculation. But for the processing of 3D reconstruction, the more information is acquired, the better the 3D point cloud model is, while the estimation accuracy of phenotypic parameters is not affected.
The processing platform was a Dell G7-7590-1865 laptop running Windows 10 (64-bit) on an Intel Core i7-9750 processor at 2.6 GHz, with 32 GB working memory, a 1 TB solid-state hard drive and an RTX2060 graphics processing unit (GPU). The implementation of the algorithm was completed with Matlab 2018b (MathWorks, Nadik, Massachusetts, USA). In this experiment, the industrial cameras were calibrated by the method proposed by Zhang [35] , with the advantage of achieving high accuracy and meeting the needs of 3D reconstruction. For camera calibration, a total of 28 images of the calibration board with a grid of 4 mm and accuracy of ±0.01 were taken by each camera, and then the internal and external parameters of the camera such as rotation matrix and translation vector were calculated using Camera Calibration Toolbox in Matlab.

3D reconstruction and phenotypic parameters estimation
In this study, Accelerated-KAZE (A-KAZE) [36] algorithm was used to extract and match features of apple images, to ensure the accuracy of 3D reconstruction and avoid too many holes in 3D point clouds. The local feature algorithm which is real-time can be applied to the mobile devices and presents a good stability under the influence factors such as illumination change, scale transformation, image rotation, noise of view angle change and so on. The process of 3D reconstruction was shown in Figure 3. The method of Fast Explicit Diffusion (FED) [37] was used to construct nonlinear scale space in the algorithm of A-KAZE, which has a faster speed and higher accuracy than the method of Additive Operator Splitting (AOS) in KAZE. The constructed nonlinear scale space can effectively distinguish uniform region from edge region and retain more detail information in the image. Firstly, the features of the sample image were detected using the A-KAZE algorithm in the nonlinear scale space, that is, the determinant of all pixels in the sample image were calculated by the Hessian matrix to get the local maximum. The feature points were located accurately by the method of Taylor expansion subsequently. Secondly, in the process of feature description, an efficient Modified Local Difference Binary (M-LDB) [38] was used to obtain 486 bit binary descriptors, which enhanced the robustness of rotation and scale invariance compared with the original LDB, and increased the uniqueness after combining the scale space gradient information constructed by FED. After that, the method of Ratio was used to set a smaller threshold to match the feature points, and then the mismatching points were eliminated by random sampling consistency (RANSAC) [39] algorithm.
Finally, the essential matrix calculated according to the basic matrix and the camera internal parameter matrix was decomposed into rotation matrix and translation vector by Singular Value Decomposition (SVD) [40] . The triangulation method was used to calculate the three-dimensional coordinates of the actual point cloud, and the sparse point clouds were generated by adding successively the new sample images for iterative calculation. Then the sparse point clouds were diffused by the algorithm of Patched-based Multi-View Stereo (PMVS) [41] to create the dense point clouds of apple. In this study, the binocular vision system was used to collect apple sequence images which were reconstructed to obtain two groups of point cloud data. After rough registration based on RANSAC, the point cloud data was processed with the algorithm of Iterative Closest Point (ICP) [42] to achieve accurate point cloud registration and form a complete apple 3D point cloud model. The flowchart of point cloud registration was shown in Figure 4. Firstly, two groups of point cloud data were divided into source point cloud set and target point cloud set. The effective point set P was extracted from the source point cloud set. The point set Q, corresponding to point-to-point P, in the target point cloud set was calculated by the method of Euclidean distance threshold. The position coordinates of the center of point set P and Q were calculated and centralized to generate a new point set p i and q i . Principal component analysis was used to perform preliminary correspondence estimation for N pairs of points in the new point set, and then the rotation matrix and the translation vector were calculated to minimize the error function E. When the error function was less than the threshold, the iterative calculation was stopped; otherwise, the above calculation process would be repeated until the condition was met. The error function E was as follows: where, n is the number of nearest point pairs; p i is the point in the source point cloud P; q i is the nearest point corresponding to p i in the target point cloud Q; R is the rotation matrix, and t is the translation vector. Figure 4 Flowchart of single-iteration point cloud registration The complete apple 3D point cloud model, registered by ICP algorithm, was filled up with the Delaunay triangulation generated by Bowyer-Watson algorithm to realize surface reconstruction. And the surface was completed to approximate the apple surface with continuous triangular faces generated by the method of Triangulated Irregular Network (TIN) using irregularly distributed point cloud data points. The flowchart of surface reconstruction was shown in Figure 5. The first step was to generate an initial network by defining a rectangle that can cover the whole computing area. And the rectangle was divided into two initial triangles according to the diagonal. Based on the principle of Delaunay division, the boundary points were added to the initial triangles one by one to set up an expanded mesh, and then the triangles outside the computing area were removed from the extended mesh. The second step was to insert a new node into the initial mesh with point-by-point insertion method, and find all the triangles whose outer circles contained the newly added nodes and remove them. One cavity was generated, and keeping its nodes in joint with the newly added nodes created a new Delaunay triangular mesh. The third step was to update the data structure by using the newly generated triangles to replace the deleted triangles, and the rest were added in the end of the array. The last step was to repeat the above steps until all the nodes were joined into the new triangular mesh.
In the process of surface reconstruction, the Local Optimization Procedure (LOP) algorithm was used to optimize the new triangles, which can update the triangles generated before and effectively avoid the cross phenomenon of the cavity.

Phenotypic parameters estimation
After establishing surface reconstruction model of apple, some phenotypic parameters of apple can be estimated and calculated, such as volume, diameter, height and fruit shape index commonly used in automatic grading equipment. Based on the Chinese national standard GB/T10651-2008 "fresh apple" [43] , apple height, diameter D and d, fruit shape index are important parameters of apple sorting, and apple volume can provide a more accurate basis for size sorting. To measure the volume of apple, the geometric information of 3D reconstruction model was used. Firstly, the apple 3D model was horizontally layered along the height of the fruit, and the hierarchical apple model was clustered by the k-means clustering algorithm [44] .
Secondly, the iterative progressive convex hull algorithm was used to extract the outer contour points to build an irregular polygon contour, which was used to establish the triangle mesh after approximating the edge contour of a real apple by the contour smoothing. And finally, the triangle meshes were used to make up a closed irregular platform, which can be added up to calculate the volume of the apple model. There were sunken areas at the upper and lower ends of apple model, thus the plane shape of the sliced layer can be regarded as ring-disk-ring from top to bottom. After patching the sliced layer, the upper and lower parts of the middle slice are quasi-circles, the area of which is calculated by Equation (2). The upper and lower part of the slice at both ends is shaped like a circular ring, the area of which is calculated by Equation (3). The relevant formulas are as follows: where, S 1 is the area of quasi-circles section in the middle of the model; S 2 is the area of circular ring section at both ends of the model; n is the total number of outer ring convex hull points of single apple model; x i and y i are the x and y coordinate of the i th point at the outer ring position; x i+1 and y i+1 is the x and y coordinate of the (i+1) th point, m is the total number of inner ring convex hull points of the single apple model; x k and y k are the x and y coordinate of the kth point at the inner ring position, and x k+1 and y k+1 is the x and y coordinate of the (k+1) th point at the inner ring position. Finally, the volume of apple 3D model was obtained by accumulating the volume of a single irregular platform by using Equation (4) [45] . The formula is as follows: where, V is the volume of apple model, z is the number of slices; s m is the slice area of layer m; s m+1 is the slice area of layer m+1, and h m is the height of m th layer. Fruit shape index is an index to estimate the quasi-round degree of apple cross section, and is defined as the ratio of minimum fruit diameter and maximum fruit diameter according to the Chinese national standard GB/T10651-2008 "fresh apple". The maximum diameter D and the minimum diameter d of the apple were calculated by finding the slice with the maximum transverse area, and carrying it out with edge detecting, contour finding, the minimum external moment drawing and size measuring. The height of the apple model was obtained by slicing longitudinally along the maximum diameter and repeating the above calculation process. The calculation formula of fruit shape index is as follows: where, β is the fruit shape index; d is the minimum apple diameter and D is the maximum apple diameter.

Visualization of model construction
In order to observe the real-time effect of data processing, the data processing results of each stage were visualized in the experiment. The image processing flowchart is shown in Figure 6. In Figure 6a, the camera was used to collect the calibration plate images to make image sets. The camera calibration toolbox in Matlab was used to detect the grid contour in the image, and then the camera calibration was completed by deleting the images with large errors, thus obtaining the internal parameters of the camera, such as rotation matrix and motion vector. The three-dimensional reconstruction process is shown in Figure 6b. The binocular vision system was used to collect the sequence images of directionally rotated apples, then the image features were extracted successively and the adjacent images were matched. After point cloud diffusion, the dense apple point cloud contour was obtained, and finally the smooth filtering processing was carried out. Figure 6c shows the process of point cloud registration and phenotypic parameter calculation. A complete apple point cloud model was established by capturing the key points to register the point cloud data, and then a continuous irregular triangular network was constructed to approximate the contour of the apple model through surface reconstruction. Then, the three-dimensional model was segmented and filled in both horizontal and vertical directions to estimate the phenotypic parameters such as volume, maximum diameter, minimum diameter, height and fruit shape index.

Phenotypic parameters measurement
A total of 200 red Fuji apples with diversified phenotypic parameters were selected as the test samples. The real values of apple volume were measured with the drainage method, and the values of the maximum diameter D, minimum diameter d and height of each apple were measured using the Vernier caliper.
The values of fruit shape index were calculated by the ratio of D and d. As shown in Figure 7, the values of phenotypic parameters of apple are presented directly by the way of dot-line statistical chart. It can be found from the chart that there is an obvious consistency among the volume, diameter, and height of the apples. As an important indicator of apple sorting, volume is the more intuitive representation of apple size than diameter and height. In addition, fruit shape index determined by the diameter D and d can accurately describe the longitudinal quasi-circle degree of the apple, and can be used to screen out the abnormal apples within a certain volume range as a supplement indicator of fruit sorting. The measured values of apple phenotypic parameters were analyzed using the method of mathematical statistics, and the results are shown in Table 1. From Table 1 and Figure 7, it can be seen that the volume values of apples used in the experiment ranged from 200 cm 3 to 600 cm 3 , the values of height range from 60 mm to 100 mm, the values of diameter D range from 70 mm to 105 mm and values of diameter d range from 65 mm to 100 mm, and the values of fruit shape index range from 0.9 to 0.99. The values of each index are distributed in a wide range, and their average values were in the middle of the range, indicating that the selected apples are evenly distributed. The rich types of sample apple, with different sizes and shapes, were used to avoid the singleness of samples and make the experimental results more persuasive and universal.

Phenotypic parameters estimation
The 3D surface reconstruction models of 200 apples were obtained by three-dimensional reconstruction method, and the apple phenotypic parameters (volume V, height h, maximum diameter D, minimum diameter d and fruit shape index β) were estimated after segmentation, layering and filling surface of the model. The estimated values and the measured values of apples were analyzed by linear regression analysis, and the results are shown in Figure 8. The values of determination coefficient R 2 and root mean square error (RMSE) were used to estimate the result of linear regression analysis. The R 2 indicates the fitting degree between the regression line and the observed values, ranging from 0 to 1. When the values of R 2 are closer to 1, it shows that the regression model has a strong fitting ability to the observed values. The RMSE is used to measure the deviation degree between the observed values and the real values. In Figure 8, the measured values, which were measured manually, were listed as the x-axis, and the estimated values calculated based on the 3D reconstruction method were listed as the y-axis. From Figure 8, it can be seen that the apple size distribution range is wide and uniform and there is a high correlation between the estimated values and the corresponding measured values of apple phenotypic parameters. The values of R 2 and RMSE of each phenotypic parameter were calculated and listed in Table 2. Among them, the correlation of volume was the highest, with R 2 of 0.9673 and RMSE of 0.03 mm. The correlation of fruit shape index was the lowest, with R 2 of 0.9035 and RMSE of 1.53 mm, due to the influence of diameter D and d. The values of correlations of other phenotypic parameters were between volume and fruit shape index. Also, the RMSE of the fruit shape index with the lowest correlation was still within the range of allowable error.  Figure 8 Regression results between the estimated values and measured values of apples" phenotypic parameters The average error values of apple phenotypic parameters were obtained by dividing the absolute value of the difference between estimated and measured values obtained from the above experimental process by the number of sample groups. The results are shown in Table 3. The average errors of apple volume, height, diameter D and d, and fruit shape index were 8.73 cm 3 , 1.43 mm, 0.90 mm and 1.23 mm, and 1.28%, respectively. According to the average error of the phenotypic parameters of apple in the Chinese national standard GB/T10651-2008 "fresh apple", the volume error is within ±10 cm 3 and the height and diameter error is within ±2 mm. In this experiment, the average error of apple shape index estimation is small, within the allowable error range.

Discussion
The contrast experiment was carried out using a system of twodimensional color imaging, which was composed of two cross-placing cameras shown in Figure 9. Among them, one was in the horizontal direction and the other one was in the vertical direction. The calibration board was used to calibrate the images, with an accuracy of ±0.01 mm and a grid width of 2 mm. Then, the distance between the horizontal and vertical directions of the image is measured according to Zhang"s calibration method, and the actual distance represented by each pixel is 0.0417 mm and 0.0426 mm, respectively. Each apple image was processed by graying and binarization to count the horizontal and vertical pixels of the apple area, and to calculate the values of height h, diameter D and d. The images captured from the top of apples presented a hole in the apple area after binarization, thus the image filling method was used, as shown in Figure 10. Finally, the average radius of the apple was calculated by the height and diameter, and the volume of the apple was estimated by the sphere formula. Based on the method, the phenotypic parameters of 200 sample apples were estimated and analyzed to calculate average errors of each index, and the results were compared with the results of 3D reconstruction. As shown in Table 4, the average errors of the phenotypic parameters directly estimated by the two-dimensional color image were significantly larger than that using the 3D reconstruction method.
The average error of volume was maximum with a value of 23.28 cm 3 , which is nearly three times larger than the method proposed in this study.  Figure 9 Apple two-dimensional color image acquisition platform Vol. 14 No. 5 187 a. Height estimation procedure b. Maximum and minimum diameter estimation procedure Figure 10 Processing flow of comparative experiment The reasons of a larger measurement error using the two-dimensional color image are as follows: (1) there are hollows at the top and bottom of the apple and the fruit shape is generally an inverted cone, which leads to a high estimated value if the sphere volume formula is used directly. (2) It needs to strictly control the vertical and horizontal relationship between the apple and the camera in the measurement of diameter and height. The position offset of camera changes the angle of the apple image and deform the diameter and height of the apple in the image. Finally, there will be a high deviation between the estimated values and the measured values. Besides, with this method, the accuracy of estimation is difficult to guarantee and the operation is more difficult and tedious when the number of samples is large. Thus, it is difficult to use and popularize this method in the apple sorting industry. In this study based on 3D reconstruction model, it is only needed to collect the sequence images of the apple to rebuild the surface contour, and estimate the values of apple phenotypic parameters by segmenting the model. The process is simple, the estimation results are accurate, and the mass operation can be realized. In addition, in terms of time, the method proposed in this paper needs to process 64 images each time, which leads to a large amount of computation, thus the estimated time is slightly longer than that of the comparison method. In the future work, the image quality and calculation ability of the algorithm will be improved to meet the actual production requirements.

Conclusions
(1) In this study, a stereoscopic vision system was established with a pair of dumbbell rollers used for carrying the apple and driving the apple to rotate in a direction. The sequence images of the apple are quickly collected by two color industrial cameras as sample data. Then the 3D model of apple was conducted based on the three-dimensional reconstruction method, the image features were extracted and matched using the A-KAZE algorithm to realize the spatial position transformation of point cloud, and the ICP algorithm was used for the point clouds registration to generate a complete 3D point cloud model of the apple.
(2) Several phenotypic parameters of apple were estimated based on the three-dimensional model after surface reconstruction, such as the volume, height, diameter, and fruit shape index by segmenting, layering and filling the three-dimensional model. Then linear regression was used to estimate the correlations between the measured and estimated values of apple phenotypic parameters. Among them, the R 2 of volume, height, diameter D and d and fruit shape index were 0.9673, 0.9332, 0.9484, 0.9391 and 0.9035, respectively. The average errors of volume, height, diameter D and d and fruit shape index were 8.73 cm 3 , 1.43 mm, 0.90 mm, 1.23 mm, and 1.28%, respectively. However, the error estimated by the processing method based on two-dimensional color image is relatively high, and the errors of volume, height, diameter D and d, and fruit shape index are 23.28, 2.16, 1.54, 2.68 and 2.73%, respectively.
(3) In this study, the apple post-harvest sorting was taken as the research object and the estimation of the phenotypic parameters of apple were carried out to provide an effective implementation scheme for the fruit sorting based on external quality. The experimental results indicated that the proposed method had a high accuracy in calculating the phenotypic parameters of apple and it was easy to operate for batch processing, which provided a reference basis for the post-harvest sorting of apple.