In-Field Estimation of Orange Number and Size by 3D Laser Scanning

agronomy Article In-Field Estimation of Orange Number and Size by 3D Laser Scanning Valeriano Méndez 1, Antonio Pérez-Romero 2, Rubén Sola-Guirado 3, Antonio Miranda-Fuentes 3,4, Francisco Manzano-Agugliaro 5,* , Antonio Zapata-Sierra 5and Antonio Rodríguez-Lizana 4 1 Department of Applied Mathematics, E.T.S.I.A, Polytechnic University of Madrid, Ciudad Universitaria, Av. Puerta de Hierro, 2, Madrid 28040, Spain; [email protected] 2 Department of Graphic Engineering, E.T.S.I.A., University of Sevilla, Utrera Road, km 1, 41013 Sevilla, Spain; [email protected] 3 Department of Rural Engineering, Campus Rabanales, University of Cordoba, Ed. Leonardo da Vinci, Ctra. Nacional IV, km 396, 14014 Córdoba, Spain; [email protected] (R.S.-G.); [email protected] (A.M.-F.) 4Department of Aerospace Engineering and Fluids Mechanics, Area of Rural Engineering, E.T.S.I.A., University of Seville, Utrera Road, km. 1, 41013 Seville, Spain; ar[email protected] 5Department of Engineering, ceiA3, University of Almeria, Almeria 04120, Spain; [email protected] *Correspondence: [email protected]; Tel.: +34-950-215693 Received: 2 November 2019; Accepted: 11 December 2019; Published: 13 December 2019   Abstract: The estimation of fruit load of an orchard prior to harvest is useful for planning harvest logistics and trading decisions. The manual fruit counting and the determination of the harvesting capacity of the field results are expensive and time-consuming. The automatic counting of fruits and their geometry characterization with 3D LiDAR models can be an interesting alternative. Field research has been conducted in the province of Cordoba (Southern Spain) on 24 ‘Salustiana’ variety orange trees—Citrus sinensis (L.) Osbeck—(12 were pruned and 12 unpruned). Harvest size and the number of each fruit were registered. Likewise, the unitary weight of the fruits and their diameter were determined (N=160). The orange trees were also modelled with 3D LiDAR with colour capture for their subsequent segmentation and fruit detection by using a K-means algorithm. In the case of pruned trees, a significant regression was obtained between the real and modelled fruit number (R 2 = 0.63, p=0.01). The opposite case occurred in the unpruned ones (p=0.18) due to a leaf occlusion problem. The mean diameters proportioned by the algorithm (72.15 ± 22.62 mm) did not present significant differences (p=0.35) with the ones measured on fruits (72.68 ± 5.728 mm). Even though the use of 3D LiDAR scans is time-consuming, the harvest size estimation obtained in this research is very accurate. Keywords: orange tree; fruit recognition; K-means; LiDAR; HDS; GNSS; yield estimation; in-field 1. Introduction The estimation of a crop in its different growth stages is essential when making decisions about harvest, storage, transport, and marketing. In the case of fruit, this estimation is commonly based on manual counts, but they are time-consuming, expensive, and come with huge errors. Therefore, it is necessary to search for automatic counting alternatives, and it is an option to count from information obtained with a camera or a sensor such as LiDAR. Machine vision systems are crucial for automatic fruit detection, the challenges of which were suggested by Sarig [ 1 ], as it is a more intuitive approach. The fruit detection must be able to be carried out under different environmental conditions and with the restrictions of shade by leaves, branches, and other immature fruits. Agronomy 2019,9, 885; doi:10.3390/agronomy9120885 www.mdpi.com/journal/agronomy Agronomy 2019,9, 885 2 of 18 In the literature, diverse works of detection in different types of fruits or harvest can be found, such as almond [ 2 ], apple [ 3 – 7 ], cherryfruit [ 8 ], cucumber [ 9 ], mango [ 10 , 11 ], orange [ 12 , 13 ], pineapple [ 14 , 15 ], or tomato [16]. Fruit detection requires segmentation, shape selection, and identification phases [ 17 ]. Segmentation consists of filtering through a colour threshold of the components of the scene that can be considered fruit [ 18 ]. Different characteristics of perimeter, area, or compaction allow to select the shapes (blobs) and to identify the fruits one by one. In bidimensional models, circles can be detected by the Hough transform [ 19 ] or by adjusting circular contours [ 20 ]. The colour cameras allowed Harrell et al. [ 12 ] to implement a robotic orange harvesting system and Grasso and Recce [ 21 ] to perform RGB segmentation. Qureshi et al. [ 11 ] used clustering of K-nearest neighbour pixels. Qureshi et al. [ 22 ] presented a texture-based method for shape recognition using an over-segmentation of super-pixels from the local gradient calculation. Colour segmentation using charge coupled device cameras returns pixels with RGB graduations or colour composition by addition of the primary colours red, green, and blue, allowing rapid detection of ripe fruit [ 23 ]. However, its drawback is its false positives. In addition, shape detection has a high computational cost. Therefore, a suitable alternative is the use of a colour filter followed by shape detection to avoid these false positives. The colour and shape characteristics allow us to approach its count using 2D photos, filtering by colour or chromaticity, delimiting shapes by contour, and being able to adjust circles. All this allows us to count the total number of fruits and even estimate their size, a major aspect when assessing yield [24]. Regarding problems in fruit detection, they can be classified into illumination and shading by leaves and branches [ 25 ]. The fruit illumination is never uniform and depends on factors such as environmental conditions and the position of the fruit within the tree. The detection of shapes, on the contrary, is very affected by the shading, mainly of the leaves, branches, or other fruits. Citrus fruits are trees with high foliar density [ 26 ], so 70%–80% of fruits have half or more of their surface darkened. Plebe and Grasso [ 27 ] used a system of stereoscopic cameras assembled on telescopic robotic arms to detect a spherical shape. The local detection of the spherical shape from the telescopic arms reduces the occlusion of the sphere. A regression of the spherical object or simply separation of blobs can be achieved depending on the quality of the segmented point cloud of the fruits. Using a LiDAR sensor to perform a 3D scan of a tree could reduce the problem of shading by constructing a point cloud of the tree that was provided of depth. However, the richness of canopy may continue to restrict the complete observation of the tree fruits. Although the LiDAR sensor allows an optimal location, it is large, expensive, and time-consuming [28]. LiDAR technology from a cloud of points allows to reconstruct the detail of a tree and, therefore, measure its elements. However, the reconstruction of all the vegetative elements of the plant due to their complexity cannot be addressed with complete accuracy. In this sense, general measures have been addressed, such as the diameter at breast height (DBH) [29,30]; tree crown determination, particularly relevant for the adaptation of farming machinery [ 31 ] or adjustment of hydro-pneumatic spray dose [ 32 – 34 ]; or the woody structure recreation [ 35 – 38 ]. The methods that allow its estimation are essentially triangulations, adjustment by cylinders, and generation of the skeleton or use of voxels. The digital recreation of canopy is not approachable except the estimation of its surface as derived from the reconstruction of the own crown of the tree [39,40], but the detailed count of the number of leaves is still a huge problem nowadays. However, fruits and other vegetative elements, such as inflorescences, have characteristics that allow the scene to be segmented by means of a colour filter. This shadowing effect significantly reduces the density of points for trying to select identifiable blobs. For this reason, it is possible to consider using a combination of a LiDAR sensor and a colour photographic sensor. Therefore, the purpose of the present study was to evaluate: first, the capacity of a K-means algorithm for fruit detection based on a 3D model of pruned and unpruned orange trees, considering Agronomy 2019,9, 885 3 of 18 the shading problems that may appear in the trees, and second, to assess the harvest forecasting ability based on the algorithm compared to actual field measurements of harvested fruits 2. Materials and Methods 2.1. Orange Trees Plots The study was conducted on 26 and 28 December 2017 in a private farm of 14 ha of irrigated orange trees at harvest time. The farm is located in Puebla de la Parrilla (Hornachuelos, Cordoba, Spain), see Figure 1A. The studied crop is a fully developed orange tree plantation of the ‘Salustiana’ variety, with a plantation frame of 6 × 4 m. For the study, a total of 24 trees were selected, divided into two blocks of 6 × 2 (Figure 1B). The first block corresponds to trees pruned in March 2017, and the second block corresponds to trees not pruned in that year and, therefore, with larger leaf density and larger crown volume (Figure 2a,b). Agronomy 2019, 9, x FOR PEER REVIEW 3 of 19 The study was conducted on 26 and 28 December 2017 in a private farm of 14 ha of irrigated orange trees at harvest time. The farm is located in Puebla de la Parrilla (Hornachuelos, Cordoba, Spain), see Figure 1A. The studied crop is a fully developed orange tree plantation of the ‘Salustiana’ variety, with a plantation frame of 6 × 4 m. For the study, a total of 24 trees were selected, divided into two blocks of 6 × 2 (Figure 1B). The first block corresponds to trees pruned in March 2017, and the second block corresponds to trees not pruned in that year and, therefore, with larger leaf density and larger crown volume (Figure 2a,b). Figure 1. (A) Location of field trials in the south of Spain. (B) Detail of trial plots; Zone 1: 2 × 6 unpruned trees. Zone 2: 2 × 6 pruned trees. The Universal Transverse Mercator coordinates in the central points of the unpruned trial plots (Zone 1 in Figure 1B) are X = 306,597 m, Y = 4,180,494 m, zone = 30 north and datum ETRS89, with an elevation of 63 m above sea-level. The coordinates of the pruned plot (Zone 2 in Figure 1B) are X = 306,401 m and Y = 4,180,720 m; the rest of the parameters remain constant. The climate in this area is typically Mediterranean, with three to five hot and dry months (June to September). Figure 1. ( A ) Location of field trials in the south of Spain. ( B ) Detail of trial plots; Zone 1: 2 × 6 unpruned trees. Zone 2: 2 ×6 pruned trees. The Universal Transverse Mercator coordinates in the central points of the unpruned trial plots (Zone 1 in Figure 1B) are X=306,597 m, Y=4,180,494 m, zone =30 north and datum ETRS89, with an elevation of 63 m above sea-level. The coordinates of the pruned plot (Zone 2 in Figure 1B) are X= 306,401 m and Y=4,180,720 m; the rest of the parameters remain constant. The climate in this area is typically Mediterranean, with three to five hot and dry months (June to September). 2.2. Data 2.2.1. Yield and Orange Sampling The harvest of each tree was picked separately to obtain the exact number of oranges from each tree individually. The weight of the orange harvest was to be estimated according to the number of oranges detected by the K-means algorithm used (Section 2.3), so a sample of 10 oranges per test tree, selected at different heights and cardinal points, was taken. In each of them, the weight was determined using an electronic weighing scale. Its diameter was also measured with an electronic gauge in the field in order to verify its concordance with the orange size results provided by the algorithm. Finally, the harvest of each tree was individually weighted. Agronomy 2019,9, 885 4 of 18 Agronomy 2019, 9, x FOR PEER REVIEW 4 of 19 Figure 2. Views of the trial plots. (A) Zone 1, unpruned trees. (B) Zone 2, pruned trees. 2.2. Data 2.2.1. Yield and Orange Sampling The harvest of each tree was picked separately to obtain the exact number of oranges from each tree individually. The weight of the orange harvest was to be estimated according to the number of oranges detected by the K-means algorithm used (Section 2.3), so a sample of 10 oranges per test tree, selected at different heights and cardinal points, was taken. In each of them, the weight was determined using an electronic weighing scale. Its diameter was also measured with an electronic gauge in the field in order to verify its concordance with the orange size results provided by the algorithm. Finally, the harvest of each tree was individually weighted. 2.2.2. Three-Dimensional Modelling Using Laser Scanner (a) Data Acquisition A laser scanner, model ScanStation P20 of Leica Geosystems, was used to obtain the point clouds. It is a compact, high-speed pulse scanner with precision, range, and vision of surveying type a ) Figure 2. Views of the trial plots. (A) Zone 1, unpruned trees. (B) Zone 2, pruned trees. 2.2.2. Three-Dimensional Modelling Using Laser Scanner (a) Data Acquisition A laser scanner, model ScanStation P20 of Leica Geosystems, was used to obtain the point clouds. It is a compact, high-speed pulse scanner with precision, range, and vision of surveying type with an integrated camera that uses Waveform Digitising technology (WFD). These specifications allow 3D position accuracies of 3 mm at 50 m distance on observed objects and the measurement of up to one million points per second. Prior to the measurements with the scanner, it was necessary to place some checkpoints using 9 steel screws in both zones (see Figure 3). The presence of the control points is essential to make the adjustment of the point clouds taken from each independent position of the laser scanner. The coordinates of these points were made using GPS (Leica 1200) in RTK mode in order to obtain their locations integrated into the geodetic system ETRS89 (European Terrestrial Reference System 1989), the official reference system for Spain in the Iberian Peninsula. According to the manufacturer, of 10 mm + 1 ppm in the horizontal position and 20 mm +1 ppm in the vertical position, which means a joint error Agronomy 2019,9, 885 5 of 18 in three dimensions of approximately 22 mm +1 ppm. To improve accuracy, 10 measurements were made at each point. The Leica 1200 software allows a basic statistical analysis of the results obtained for each point, deleting those points in which the deviation from the mean had an error of more than 5 mm. 1 Figure 3. Distribution of checkpoints and scanner positions (Zone 2). Above the points materialized by the checkpoints, topographical Prism Poles with HDS (high definition surveying) targets were subsequently placed (see Figure 4). These points are intended to be measured from each position occupied by the 3D Laser Scanner. In this way, it is possible to integrate all the scans made and to obtain the correct georeferencing of the entire survey. Agronomy 2019, 9, x FOR PEER REVIEW 6 of 19 Figure 4. High definition surveying (HDS) targets on topographical Prism Poles. To ensure accuracy, it is necessary to measure at least two of these checkpoints from each position, but it is best to use three, to perform more robust adjustments. For this reason, the checkpoints were placed in groups of three, in the central area of the crop streets, as can be seen in Figure 3. In order to cover all the trees as optimally as possible, the scanner had to be strategically positioned (Figure 5). One of the main drawbacks of the laser scanner-based HDS measurement technique is that of the “shadow” areas caused by occlusions that prevent complete surveying of the entire surface from a single position, forcing repositioning and new scans to achieve optimal coverage [41]. In this case, in order to cover each tree, and given the environmental conditions, it was required to use four different positions of the scanner, as shown in Figure 5. Figure 5. Example of scanner positions to cover the top of a tree. In a 3D laser scan survey, the density of the obtained point cloud and, consequently, its resolution, is dependent on the scanning distance and angle [42]. In this study, considering the distribution of the scanner positions, the plantation framework and the size of the trees, the distances to the measurement objects were small (between 1 and 3 m). Figure 4. High definition surveying (HDS) targets on topographical Prism Poles. Agronomy 2019,9, 885 6 of 18 To ensure accuracy, it is necessary to measure at least two of these checkpoints from each position, but it is best to use three, to perform more robust adjustments. For this reason, the checkpoints were placed in groups of three, in the central area of the crop streets, as can be seen in Figure 3. In order to cover all the trees as optimally as possible, the scanner had to be strategically positioned (Figure 5). One of the main drawbacks of the laser scanner-based HDS measurement technique is that of the “shadow” areas caused by occlusions that prevent complete surveying of the entire surface from a single position, forcing repositioning and new scans to achieve optimal coverage [ 41 ]. In this case, in order to cover each tree, and given the environmental conditions, it was required to use four different positions of the scanner, as shown in Figure 5. Agronomy 2019, 9, x FOR PEER REVIEW 6 of 19 Figure 4. High definition surveying (HDS) targets on topographical Prism Poles. To ensure accuracy, it is necessary to measure at least two of these checkpoints from each position, but it is best to use three, to perform more robust adjustments. For this reason, the checkpoints were placed in groups of three, in the central area of the crop streets, as can be seen in Figure 3. In order to cover all the trees as optimally as possible, the scanner had to be strategically positioned (Figure 5). One of the main drawbacks of the laser scanner-based HDS measurement technique is that of the “shadow” areas caused by occlusions that prevent complete surveying of the entire surface from a single position, forcing repositioning and new scans to achieve optimal coverage [41]. In this case, in order to cover each tree, and given the environmental conditions, it was required to use four different positions of the scanner, as shown in Figure 5. Figure 5. Example of scanner positions to cover the top of a tree. In a 3D laser scan survey, the density of the obtained point cloud and, consequently, its resolution, is dependent on the scanning distance and angle [42]. In this study, considering the distribution of the scanner positions, the plantation framework and the size of the trees, the distances to the measurement objects were small (between 1 and 3 m). Figure 5. Example of scanner positions to cover the top of a tree. In a 3D laser scan survey, the density of the obtained point cloud and, consequently, its resolution, is dependent on the scanning distance and angle [ 42 ]. In this study, considering the distribution of the scanner positions, the plantation framework and the size of the trees, the distances to the measurement objects were small (between 1 and 3 m). Then, looking for a proper compromise between the resolution required for the work and the time available to complete it in the field, it was necessary to configure the area to be measured from each of them, especially in the parameters of the scan resolution and image taking. Covered angles scanning between 150 and 170 sexagesimal degrees were in the positions of the outer streets of the crop, but the central street was chosen to perform full circle scans to cover the trees on both sides of the axis of the street from one position. Note that images are not necessary to obtain point clouds, but by allowing them to be assigned the colours of reality captured by storing RGB values, they make it possible to identify differentiated elements [ 43 ], which is especially useful in this work, since ripe fruits (oranges) are clearly differentiated from the rest of the tree by colour. The scanning parameters are shown in Figure 6. The scans made from the outside streets had an average time of 5 min, while in the central street, working with full circle sweeps, the times were 13.72 min. Agronomy 2019,9, 885 7 of 18 Agronomy 2019, 9, x FOR PEER REVIEW 7 of 19 Then, looking for a proper compromise between the resolution required for the work and the time available to complete it in the field, it was necessary to configure the area to be measured from each of them, especially in the parameters of the scan resolution and image taking. Covered angles scanning between 150 and 170 sexagesimal degrees were in the positions of the outer streets of the crop, but the central street was chosen to perform full circle scans to cover the trees on both sides of the axis of the street from one position. Note that images are not necessary to obtain point clouds, but by allowing them to be assigned the colours of reality captured by storing RGB values, they make it possible to identify differentiated elements [43], which is especially useful in this work, since ripe fruits (oranges) are clearly differentiated from the rest of the tree by colour. The scanning parameters are shown in Figure 6. The scans made from the outside streets had an average time of 5 min, while in the central street, working with full circle sweeps, the times were 13.72 min. Figure 6. Scanning parameters. Once the scan was made in a certain position, the equipment moved to a new position, levelling the equipment and measuring the control points placed in the control points. This is necessary for the connection between each of the obtained point clouds, consuming an additional average time of 15 min. Thus, measuring each zone took 8 h. (b) Data Processing As the trees were scanned from various angles, from each position of the laser scanner, three of the previously established control points were measured, which allowed the software used (Cyclone) to be able to move and rotate the clouds of points obtained independently in the field so that they can be assembled in the common coordinate system, materialized by the control points. The connection between the point clouds was made using Cyclone software (Leica Geosystems, St. Gallen, Switzerland), achieving internal adjustments between these clouds with maximum errors of 1 mm and adjustments with the coordinates obtained for the control points (measured with GPS techniques) with errors ranging from 5 to 9 mm. The point clouds were exported in PTS format, which conforms to the following structure: the first column corresponds to the X coordinate, the second to the Y coordinate, the third to the Z coordinate, the fourth to the reflectivity and the last three correspond to the values of red, green, and blue (RGB). Figure 7 shows the point cloud generated for one of the scanned trees, and an enlarged detail showing some oranges. Figure 6. Scanning parameters. Once the scan was made in a certain position, the equipment moved to a new position, levelling the equipment and measuring the control points placed in the control points. This is necessary for the connection between each of the obtained point clouds, consuming an additional average time of 15 min. Thus, measuring each zone took 8 h. (b) Data Processing As the trees were scanned from various angles, from each position of the laser scanner, three of the previously established control points were measured, which allowed the software used (Cyclone) to be able to move and rotate the clouds of points obtained independently in the field so that they can be assembled in the common coordinate system, materialized by the control points. The connection between the point clouds was made using Cyclone software (Leica Geosystems, St. Gallen, Switzerland), achieving internal adjustments between these clouds with maximum errors of 1 mm and adjustments with the coordinates obtained for the control points (measured with GPS techniques) with errors ranging from 5 to 9 mm. The point clouds were exported in PTS format, which conforms to the following structure: the first column corresponds to the Xcoordinate, the second to the Ycoordinate, the third to the Zcoordinate, the fourth to the reflectivity and the last three correspond to the values of red, green, and blue (RGB). Figure 7shows the point cloud generated for one of the scanned trees, and an enlarged detail showing some oranges. Agronomy 2019, 9, x FOR PEER REVIEW 8 of 19 Figure 7. Partial view of one of the point clouds corresponding to an orange tree. 2.3. K-Means Algorithm Application 2.3.1. Data Segmentation The full reconstruction of the 3D point cloud obtained with the LiDAR sensor is voluminous and extremely slow [28]. The available methods allow the reconstruction of the woody structure or at the local level of other structures, such as leaves or fruits. However, a joint reconstruction of a fruit tree is not usually suitable because the combined mass of leaves, branches, fruits, and flowers does not allow to group and identify each of the elements. Therefore, it was decided to make a colour segmentation based on a direct filter by RGB (R ≥ 180, G ≥ 50, B ≤ 105) or according to the red chromaticity (r) (Equation (1))_ for those values with r ≥ 0.5 [25]. r = R / (R + G + B) (1) 2.3.2. Algorithm Segmentation determines a new cloud of points belonging to interceptions of tree fruits. Obturation problems prevent complete hemispheres from being detected, which can make it difficult to identify oranges. Instead, spherical shapes are obtained, which, in general, constitute a set of blobs that must be individually selected and subsequently identified. The selection is made according to an ISODATA process [44], which is a K-means where the number of groups or forms are not a priori defined but change during the process. Ridler and Calvard [45] use the ISODATA algorithm to segment images by grey level; Unay and Gosselin [46] or Shruthi [47] identified apples with that algorithm and using 2D colour images. The following steps apply to the process: Generate a KD-Tree to facilitate point searches. Generate a neighbourhood matrix with rn < 20 mm that is supported by the previous KD-Tree, where rn is the maximum neighbourhood distance. Initial generation of blobs grouping the neighbouring points to each other. Readjustment of shapes, merging or dividing the initial shapes by introducing the minimum or maximum dimensions of a fruit (rmin = 30 mm and rmax = 45 mm, respectively). Figure 7. Partial view of one of the point clouds corresponding to an orange tree. Agronomy 2019,9, 885 8 of 18 2.3. K-Means Algorithm Application 2.3.1. Data Segmentation The full reconstruction of the 3D point cloud obtained with the LiDAR sensor is voluminous and extremely slow [ 28 ]. The available methods allow the reconstruction of the woody structure or at the local level of other structures, such as leaves or fruits. However, a joint reconstruction of a fruit tree is not usually suitable because the combined mass of leaves, branches, fruits, and flowers does not allow to group and identify each of the elements. Therefore, it was decided to make a colour segmentation based on a direct filter by RGB (R ≥ 180, G ≥ 50, B ≤ 105) or according to the red chromaticity (r) (Equation (1))_ for those values with r≥0.5 [25]. r=R/(R +G+B) (1) 2.3.2. Algorithm Segmentation determines a new cloud of points belonging to interceptions of tree fruits. Obturation problems prevent complete hemispheres from being detected, which can make it difficult to identify oranges. Instead, spherical shapes are obtained, which, in general, constitute a set of blobs that must be individually selected and subsequently identified. The selection is made according to an ISODATA process [ 44 ], which is a K-means where the number of groups or forms are not a priori defined but change during the process. Ridler and Calvard [ 45 ] use the ISODATA algorithm to segment images by grey level; Unay and Gosselin [ 46 ] or Shruthi [47] identified apples with that algorithm and using 2D colour images. The following steps apply to the process: Generate a KD-Tree to facilitate point searches. Generate a neighbourhood matrix with r n <20 mm that is supported by the previous KD-Tree, where rnis the maximum neighbourhood distance. Initial generation of blobs grouping the neighbouring points to each other. Readjustment of shapes, merging or dividing the initial shapes by introducing the minimum or maximum dimensions of a fruit (rmin =30 mm and rmax =45 mm, respectively). The point distribution of the shapes is balanced by applying the K-means algorithm [ 48 ]. The objective is that the points are distributed among the different defined shapes minimizing the distances to the centre of gravity. Eliminate those shapes that have a number of points less than a given number (ns<9). The KD-Tree splits the point cloud with planes perpendicular to the coordinate axes so that half the points are on each side of the divisor plane. If “k” divisions are made in a d-dimensional point cloud, a d k+1 partition will be done, which reduces the computational costs of searches when generating the neighbourhood matrix in step 2. The computational cost of the process of generating a KD-Tree data structure is given by Equation (2). Cost =O (N ×log(N)) (2) where N is the total number of points in the initial cloud [ 49 ], and O indicates that the computational cost is less than N ×log (N). Given two points P i and P j of the point cloud, if the Euclidean distance between them is 0 <P i P j <r n , then the two points are neighbours. It is indicated by doing M ij =M ji =1, M being a matrix of dimension N × N, and taking the value 0 otherwise. In order to reduce memory consumption, an N-dimensional vector is adopted by storing the positions of neighbouring vertices in the position M i = “j, . . . ” and in the position Mj=“I, . . . ”. As an example, M 2 =“7, 8, 23” indicates that vertex 2 has vertexes 7, 8, and 23 as neighbours. Obtaining the neighbourhood matrix requires each point of the cloud to access the rest to check those Agronomy 2019,9, 885 9 of 18 that are at a distance less than r n . Since the cost of accessing a point on the KD-Tree according to Cormen et al. [ 49 ] is O (log(N)), the computational cost of obtaining the neighbourhood matrix is O (N ×log(N)). The generation of the initial blob groups all the neighbouring points to each other, obtaining a first distribution by compaction criterion. The process will go through the N elements of the neighbourhood matrix M i =“j,...” so that the points j,... will be neighbours to each other and form a blob, but it will be necessary to add to the form the neighbours, if any, of Mj. The readjustments in the process take place in the following stages, iterating as long as some form is affected, that is, as long as an ungrouping or grouping has taken place: Ungrouping of shapes that exceed the maximum size. Those shapes where the distance between its centre of gravity and the furthest point (outsider) exceeds the maximum distance (45 mm) are ungrouped. The ungrouping process consists of creating a new shape with the outsider point and then balancing points according to the K-means algorithm [ 48 ], which, in this case, would be reduced to 2-means. Grouping of two shapes in which the distance from the outsider point to the centre of gravity is less than the minimum distance (30 mm) and in which the grouping does not exceed the maximum distance (45 mm). Assuming there are Nb shapes, go through the list of shapes (with an O(N b ) cost) to identify those that have a dimension greater than the maximum. When this condition is met, the group is split with the creation of a new group with the outsider point. A 2-means operation is then performed to minimize the distances to the centres of gravity of the two resulting shapes. The 2-means clustering has the following cost: Cost =O(Nb×log(Nb)×log(N/Nb)). (3) As discussed below for step 5, but being a condition-based process, the total cost of the process will remain as: Cost =O(Nb×log(Nb×log(Nb)×log(N/Nb))), (4) which is under the level O (N ×log(N)). The grouping runs through the list of shapes, so it starts from an O(N b ) cost and checks those shapes that are nearby to see if they can join. The cost of this last operation is O(N/Nb), which results in an estimated O(N) cost. In step 5, the K-means algorithm is applied to the whole Nb set of existing shapes. Initially, the list of shapes is crossed and in each one of them, those in which the centres are closer, are selected. An O(Nb×log(Nb)) cost is estimated for this process. Subsequently, a point checking operation is executed between the two forms. This implies the route of the list of points of one of the forms with an average cost of O(N/N b ). And if a point had to be exchanged from one form to the other, the recalculation of the two centres of gravity with an O(N/N b ) cost is required. It is estimated as O(log(N/N b )) because it is a cost under condition. Therefore, it would be O(N b× log(N b ) × log(N/N b )), which, in general, will be below the O(N × log(N)) level. Finally, the debugging of forms that do not have significant points have an O(Nb) cost. Therefore, the expected cost of the algorithm is below the O(N ×log (N)) level. The process, therefore, consists of the following parameters: maximum neighbourhood distance (r n ), minimum distance of the shape (r min ), maximum distance of the shape (r max ), and minimum number of points of a shape (n s ). Nis the number of points obtained after filtering by colour, and Nb is the number of shapes in which the Npoints of the cloud are grouped. As the space dimension “d” is used (d=3) and “k” in the kd-tree structure is the number of levels in which the space is partitioned (k =2 is adopted), and “K” in K-means indicates the number of groups, so K=Nbin this context. Agronomy 2019,9, 885 16 of 18 5. Conclusions The use of 3D scanning and a K-means algorithm proved effective in estimating the crop structure and also the position and the size of the fruits. There is a significant linear correlation between the number of oranges in the pruned trees and those estimated with the K-means algorithm presented. The acquisition of this relationship between both variables is necessary because the direct use of the algorithm does not provide precise results. The algorithm provides accurate orange diameter estimates, which is of high interest in determining the size of the fruits of the harvest. An experimental linear relationship between the diameter of the fruits and their weight, with high determination coefficient, has been obtained. The knowledge of this experimental relationship, and of the relationship between the number of real oranges and the number modelled with the algorithm allows us to obtain crop estimates with a small margin of error. Although the use of a 3D laser scanning is a highly time-consuming method nowadays, it opens up new perspectives in the yield assessment of oranges and orange count. Likewise, due to its high precision, it could serve as a reference method to compare with other cheaper and faster methods in the future, such as photogrammetry-SFM (structure from motion) and automatic image detection techniques based on machine learning, or new faster 3D scanning equipment. In summary, the results of this study are interesting and needed as they show potential for uses of 3D laser scanning in the future for yield estimation and, thus, represent valuable basics for further study and development of 3D laser scanning. Author Contributions: Conceptualization, A.R.-L., V.M.; methodology, V.M. and A.P.-R.; resources: A.P.-R. and F.M.-A., data collection, A.P.-R., R.S.-G., A.M.-F., A.R.-L.; formal analysis: V.M., A.R.-L., A.P.-R., R.S.-G., A.M.-F.; software: V.M., A.P.-R.; writing—original draft, V.M., F.M.-A., A.R.-L., A.P.-R., A.Z.-S.; supervision, A.R.-L.; writing—review and editing, F.M.-A., R.S.-G., V.M., A.Z.-S., A.M.-F. Funding: This research received no external funding. Acknowledgments: The authors would like to thank the owners of the orange tree grove for his collaboration. Conflicts of Interest: The authors declare no conflict of interest. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the authors of this article. References 1. Sarig, Y. Robotics of Fruit Harvesting: A State-of-the-art Review. J. Agric. Eng. Res. 1993 ,54, 265–280. [CrossRef] 2. Underwood, J.P.; Hung, C.; Whelan, B.; Sukkarieh, S. Mapping almond orchard canopy volume, flowers, fruit and yield using LiDAR and vision sensors. Comput. Electron. Agric. 2016,130, 83–96. [CrossRef] 3. Grand d’Esnon, A. Robotic Harvesting of Apples. In Proceedings of Agri-Mation 1 Conference and Exposition; ASAE: St. Joseph, MI, USA, 1985; pp. 210–214. 4. Kassay, L.; Slaughter, D.C.; Molnar, S. Hungarian Robotic Apple Harvester; ASAE: St. Joseph, MI, USA, 1992. 5. Bulanon, D.M.; Burks, T.F.; Alchanatis, V. Study on Fruit Visibility for Robotic Harvesting; ASABE Paper No. 07-3124; ASABE: St. Joseph, MI, USA, 2007. 6. Tabb, A.; Peterson, D.; Park, J. Segmentation of Apple Fruit from Video via Background Modeling; ASAE: St. Joseph, MI, USA, 2006. 7. Bargoti, S.; Underwood, J.P. Image segmentation for fruit detection and yield estimation in apple orchards. J. Field Robot. 2017,34, 1039–1060. [CrossRef] 8. Tanigaki, K.; Fujiura, T.; Akase, A.; Imagawa, I. Cherry Harvesting Robot. In Proceedings of the 3rd IFAC International Workshop on Bio-robotics, Information Technology and Intelligent Control for Bioproduction Systems, Sapporo, Japan, 9–10 September 2006; pp. 254–260. 9. Van Henten, E.J.; Hemming, J.; Van Tuijl, B.A.J.; Kornet, J.G.; Bontsema, J. Collision-free Motion Planning for a Cucumber Picking Robot. Biosyst. Eng. 2003,86, 135–144. [CrossRef] 10. Payne, A.B.; Walsh, K.B.; Subedi, P.P.; Jarvis, D. Estimation of mango crop yield using image analysis–segmentation method. Comput. Electron. Agric. 2013,91, 57–64. [CrossRef] Agronomy 2019,9, 885 17 of 18 11. Qureshi, W.S.; Payne, A.; Walsh, K.B.; Linker, R.; Cohen, O.; Dailey, M.N. Machine vision for counting fruit on mango tree canopies. Precis. Agric. 2017,18, 224–244. [CrossRef] 12. Harrell, R.C.; Adsit, P.D.; Pool, T.A.; Hoffman, R. The Florida Robotic Grove Lab; ASAE Paper no. 88–1578; ASAE: St. Joseph, MI, USA, 1988. 13. Hannan, M.W.; Burks, T.F. Current Developments in Automated Citrus Harvesting; ASAE Paper no., 043087; ASAE: St. Joseph, MI, USA, 2004. 14. Chaivivatrakul, S.; Moonrinta, J.; Dailey, M.N. Towards Automated Crop Yield Estimation-Detection and 3D Reconstruction of Pineapples in Video Sequences. In Proceedings of the Fifth International Conference on Computer Vision Theory and Applications, Angers, France, 17–21 May 2010; pp. 180–183. 15. Chaivivatrakul, S.; Dailey, M.N. Texture-based fruit detection. Precis. Agric. 2014,15, 662–683. [CrossRef] 16. Kondo, N.; Nishitsuji, Y.; Ling, P.P.; Ting, K.C. Visual Feedback Guided Robotic Cherry Tomato Harvesting. Trans. ASAE 1996,39, 2331–2338. [CrossRef] 17. Clement, J.; Novas, N.; Gazquez, J.A.; Manzano-Agugliaro, F. An active contour computer algorithm for the classification of cucumbers. Comput. Electron. Agric. 2013,92, 75–81. [CrossRef] 18. Novas, N.; Á lvarez-Bermejo, J.A.; Valenzuela, J.L.; G á zquez, J.A.; Manzano-Agugliaro, F. Development of a smartphone application for assessment of chilling injuries in zucchini. Biosyst. Eng. 2019 ,181, 114–127. [CrossRef] 19. Whittaker, D.; Miles, G.E.; Mitchell, O.R.; Gaultney, L.D. Fruit location in a partially occluded image. Trans. ASAE 1987,30, 591–596. [CrossRef] 20. Pla, F. Recognition of partial circular shapes from segmented contours. Comput. Vis. Image Undst. 1996 ,63, 334–343. [CrossRef] 21. Grasso, G.M.; Recce, M. Scene analysis for an orange harvesting robot. AI Appl. 1997,11, 9–15. 22. Qureshi, W.S.; Satoh, S.I.; Dailey, M.N.; Ekpanyapong, M. Dense segmentation of textured fruits in video sequences. In Proceedings of the IEEE International Conference on Computer Vision Theory and Applications (VISAPP), Lisbon, Portugal, 5–8 January 2014; Volume 2. 23. Jim é nez, A.R.; Ceres, R.; Pons, J.L. A survey of computer vision methods for locating fruit on trees. Trans. ASAE 2000,43, 1911. [CrossRef] 24. Sharifi, M.; Rafiee, S.; Keyhani, A.; Jafari, A.; Mobli, H.; Rajabipour, A.; Akram, A. Some physical properties of orange (var. Tompson). Int. Agrophys. 2007,21, 391–397. 25. Hannan, M.W.; Burks, T.F.; Bulanon, D.M.A. A machine vision algorithm combining adaptive segmentation and shape analysis for orange fruit detection. Agric. Eng. Int. CIGR J. 2010,11, 1281. 26. Zapata-Sierra, A.J.; Manzano-Agugliaro, F. Controlled deficit irrigation for orange trees in Mediterranean countries. J. Clean. Prod. 2017,162, 130–140. [CrossRef] 27. Plebe, A.; Grasso, G. Localization of spherical fruits for robotic harvesting. Mach. Vis. Appl. 2001 ,13, 70–79. [CrossRef] 28. Gongal, A.; Amatya, S.; Karkee, M.; Zhang, Q.; Lewis, K. Sensors and systems for fruit detection and localization: A review. Comput. Electron. Agric. 2015,116, 8–19. [CrossRef] 29. Lefsky, M.A.; Harding, D.; Cohen, W.B.; Parker, G.; Shugart, H.H. Surface lidar remote sensing of basal area and biomass in deciduous forests of eastern Maryland, USA. Remote Sens. Environ. 1999 ,67, 83–98. [CrossRef] 30. Huang, H.; Li, Z.; Gong, P.; Cheng, X.; Clinton, N.; Cao, C.; Wenjian, N.; Wang, L. Automated methods for measuring DBH and tree heights with a commercial scanning lidar. Photogramm. Eng. Remote Sens. 2011 ,77, 219–227. [CrossRef] 31. Sola-Guirado, R.R.; Bayano-Tejero, S.; Rodr í guez-Lizana, A.; Gil-Ribes, J.A.; Miranda-Fuentes, A. Assessment of the Accuracy of a Multi-Beam LED Scanner Sensor for Measuring Olive Canopies. Sensors 2018 ,18, 4406. [CrossRef] [PubMed] 32. Miranda-Fuentes, A.; Rodr í guez-Lizana, A.; Gil, E.; Agüera-Vega, J.; Gil-Ribes, J.A. Influence of liquid-volume and airflow rates on spray application quality and homogeneity in super-intensive olive tree canopies. Sci. Total Environ. 2015,537, 250–259. [CrossRef] [PubMed] 33. Miranda-Fuentes, A.; Llorens, J.; Rodr í guez-Lizana, A.; Cuenca, A.; Gil, E.; Blanco- Rold á n, G.L.; Gil-Ribes, J.A. Assessing the optimal liquid volume to be sprayed on isolated olive trees according to their canopy volume. Sci. Total Environ. 2016,568, 296–305. [CrossRef] Agronomy 2019,9, 885 18 of 18 34. Miranda-Fuentes, A.; Rodr í guez-Lizana, A.; Cuenca, A.; Gonz á lez-S á nchez, E.J.; Blanco-Rold á n, G.L.; Gil-Ribes, J.A. Improving plant protection product applications in traditional and intensive olive orchards through the development of new prototype air-assisted sprayers. Crop Prot. 2017,94, 44–58. [CrossRef] 35. Verroust, A.; Lazarus, F. Extracting skeletal curves from 3D scattered data. Vis. Comput. 2000 ,16, 15–25. [CrossRef] 36. Mizoue, N.; Masutani, T. Image analysis measure of crown condition, foliage biomass and stem growth relationships of Chamaecyparis obtusa. For. Ecol. Manag. 2003,172, 79–88. [CrossRef] 37. Pfeifer, N.; Gorte, B.; Winterhalder, D. Automatic reconstruction of single trees from terrestrial laser scanner data. In Proceedings of the 20th ISPRS Congress: Geo-imagery Bridging Continents, Istanbul, Turkey, 12–23 July 2004; pp. 114–119. 38. Phattaralerphong, J.; Sinoquet, H. Tree Analyser: Software to Compute Tree Structure Parameters from Photographs. User Manual. PIAF-INRA. 2007. Available online: http://www2.clermont.inra.fr/piaf/eng/ download/download.php (accessed on 1 August 2019). 39. Phattaralerphong, J.; Sinoquet, H. A method for 3D reconstruction of tree crown volume from photographs: Assessment with 3D-digitized plants. Tree Physiol. 2005,25, 1229–1242. [CrossRef] 40. Bailey, B.N.; Ochoa, M.H. Semi-direct tree reconstruction using terrestrial LiDAR point cloud data. Remote Sens. Environ. 2018,208, 133–144. [CrossRef] 41. Lim, M.; Rosser, N.; Allison, R.; Petley, D. Erosional processes in the hard rock coastal cliffs at Staithes. Geomorphology 2010,114, 12–21. [CrossRef] 42. Garc í a-G ó mez, I.; de Gorostiza, M.; Moraza, A. L á ser esc á ner y nubes de puntos. Un horizonte aplicado al análisis arqueológico de edificios. Arqueol. Arquit. 2010,30, 25–44. [CrossRef] 43. Bened í , J.; Moreno, J.; Molina, C.; Arguedas, G.; Serrano, S.; Oliv á n, A. Utilizaci ó n de t é cnicas de l á ser esc á ner terrestre en la monitorizaci ó n de procesos geomorfol ó gicos din á micos: El manto de nieve y heleros en á reas de montaña. Cuadernos de investigación geográfica 2013,39, 335–357. 44. Ball, G.H.; Hall, D.J. ISODATA, A Novel Method of Data Analysis and Pattern Classification; Stanford Research Inst.: Menlo Park, CA, USA, 1965. 45. Ridler, T.W.; Calvard, S. Picture thresholding using an iterative selection method. IEEE Trans. Syst. Man Cybern. 1978,8, 630–632. 46. Unay, D.; Gosselin, B. Thresholding-based segmentation and apple grading by machine vision. In Proceedings of the 13th European Signal Processing Conference, Antalya, Turkey, 4–8 September 2005. 47. Shruthi, G. Automatic Detection and Classification of Apple-A Survey. In Proceedings of the Third International Conference on Current Trends in Engineering Science and Technology (ICCTEST 2017), Bangalore, India, 5–7 January 2017; Don Bosco Institute of Technology: Dbit, Bangalore, India. 48. Lloyd, S.P. Least squares quantization in PCM. IEEE Trans. Inform. Theory 1982,28, 129–137. [CrossRef] 49. Cormen, T.H.; Leiserson, C.E.; Rivest, R.L.; Stein, C. Introduction to Algorithms, 3rd ed.; The MIT Press: London, UK, 2009; pp. 248–300. 50. Goldberg, K.M.; Iglewicz, B. Bivariate extensions of the boxplot. Technometrics 1992 ,34, 307–320. [CrossRef] 51. R Core Team. R: A Language and Environment for Statistical Comput-Ing; R Foundation for Statistical Computing: Vienna, Austria, 2015. 52. Stein, M.; Bargoti, S.; Underwood, J. Image Based Mango Fruit Detection, Localisation andYield Estimation Using Multiple View Geometry. Sensors 2016,16, 1915. [CrossRef] [PubMed] 53. Liu, J.; Yuan, Y.; Zhou, Y.; Naz Syed, T. Experiments and Analysis of Close-Shot Identification of On-Branch Citrus Fruit with RealSense. Sensors 2018,18, 1510. [CrossRef] [PubMed] 54. Tabatabaeefar, A.; Vefagh-Nematolahee, A.; Rajabipour, A. Modeling of orange mass based on dimensions. Agric. Sci. Technol. 2000,2, 299–305. © 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).