|
Artificial Medjool Date Fruit Bunch Image Synthesis: Towards Thinning Automation
Or Bar-Shira1, Yosef Cohen1, Tal Shoshan1, Avital Bechar1,2, Avraham Sadowsky3, Yuval Cohen4, Sigal Berman1,*
Published in Journal of the ASABE 66(2): 275-284 (doi: 10.13031/ja.15217). Copyright 2023 American Society of Agricultural and Biological Engineers.
1The Department of Industrial Engineering and Management, Ben-Gurion University of the Negev, Beer-Sheva, Israel.
2Institute of Agricultural Engineering, Agricultural Research Organization, Volcani Center, Rishon LeZion, Israel.
3Southern Arava Research and Development, Israel.
4Institute of Plant Sciences, Agricultural Research Organization, Volcani Center, Rishon LeZion, Israel
*Correspondence: sigalbe@bgu.ac.il
Submitted for review on 1 June 2022 as manuscript number ITSC 15217; approved for publication as a Research Article and as part of the Artificial Intelligence Applied to Agricultural and Food Systems Collection by Associate Editor Dr. Yuzhen Lu and Community Editor Dr. Yiannis Ampatzidis of the Information Technology, Sensors, & Control Systems Community of ASABE on 19 January 2023.
Highlights
Abstract. Medjool is a premium date cultivar, and the market demands high-quality fruits, for which specific horticultural practices, including timely and efficient fruitlet thinning, are required. Currently, thinning the fruitlets is one of the most labor-intensive tasks in the Medjool cultivation cycle, and there is a need to develop methods for automating the thinning process. An algorithm determining the required thinning is a prerequisite for advancing toward thinning automation. An annotated Medjool fruit bunch image dataset is necessary for developing such an algorithm using state-of-the-art machine learning methods. Acquiring such a dataset is difficult and costly. The difficulty can be alleviated by using synthetic images. However, current methods for generating synthetic plant images are unsuitable for Medjool dates due to their geometrical features. The current work suggests a method for generating artificial images of Medjool fruit bunches from a 3D model based on structural decomposition into plant parts and the use of Bezier curves. Nineteen model variables and their distributions were defined for fruit bunch model generation. The models and synthetic images generated based on the models were verified by two plant physiologists who are experts in Medjool date cultivation. Fruit-bunch features were extracted from the generated images and used for learning the required remaining length of the spikelets after thinning using kernel estimation. The estimation was tested for images of two whorl-period combinations (Top-Early and Middle-Middle). The average scaled absolute estimation errors for both periods were very low (less than 1%).
Dates are a health food product with an annual market of about 9.5 million tons worldwide (FAOSTAT, 2022). Medjool is a premium date cultivar with increasing annual demand, currently approximately 100,000 tons. Market expansion is impeded by the intensive labor required for cultivating the Medjool crop. One of the most labor-intensive tasks in the Medjool production cycle is thinning the fruit bunches, which requires approximately 3.5 person-hours per tree (25% of the total required annual labor) (Cohen and Glasner, 2015). Therefore, automating the thinning process is a promising avenue for reducing the labor needed in the Medjool production cycle. Developing an algorithm for determining the required thinning of the Medjool fruit bunch is a critical prerequisite for thinning automation (Shoshan et al., 2022).
Determining the required thinning, i.e., the required quantity of remaining fruitlets on the fruit bunch during the thinning period, for attaining a high-quality yield later on at harvest is challenging. Developing an algorithm for determining the required thinning based on images of fruit bunches within the canopy is complicated, partly because state-of-the-art algorithms for image processing require a large amount of annotated data suitable for training. For example, an annotated Medjool fruit bunch image dataset taken during the thinning period when the fruit bunches are located within the canopy has recently been published (Shoshan and Berman, 2021). The dataset includes 146 color images, which is a relatively small number for training deep neural networks.
The lack of datasets suitable for training advanced algorithms is typical in the agricultural domain. It is considered one of the main reasons for the lower performance of image-processing algorithms in agriculture and one of the main obstacles to introducing robotics into agriculture (Barth et al., 2020; Lu and Young, 2020). Researchers have suggested overcoming this lack of suitable image datasets by artificially generating agricultural image data, and several methods have been developed for such syntheses. For example, Barth et al. (2018) developed a synthesis method for large-scale semantic image segmentation datasets and implemented it to construct a greenhouse bell pepper image dataset. The synthesis is based on structurally decomposing the plant into parts and then generating randomized plant meshes by concatenating instances of plant organs. Taking a different approach, Benoit et al. (2014) simulated plant growth using an adapted root-growing simulation based on an L-systems process (Leitner et al., 2010). Where, L-systems, as iteratively created strings, are highly suitable for representing the growth of multiple link structures. However, these approaches are unsuitable for modeling Medjool fruit bunches due to their geometric features. The current research suggests a different modeling approach suitable for modeling date fruit bunches. To describe the proposed modeling method and demonstrate its suitability for Medjool fruit bunches, it is first necessary to present a short review of the palm structure and spatial organization and the Medjool cultivation practices.
Medjool fruit bunches grow in whorls within the date crown at the top of the palm tree (Chao and Krueger, 2007). They develop at the base of the date leaves (fronds) that have grown in the previous year. Farmers divide the fruit bunches into three typical whorls (top, middle, or bottom) according to their location in the crown. The higher bunches tend to be larger and develop earlier than the lower bunches. There are significant variations in bunch size within the same tree and between trees of different ages, with the variations also depending on growth conditions and horticultural management protocols. Each fruit bunch consists of a fruit stalk, i.e., the fruit bunch rachis, spikelets (unbranched rachillae), and fruitlets (fig. 1). The rachis grows upward, tilted with respect to the palm trunk, and as the fruitlets become larger and heavier, the rachis tilt angle increases. The spikelets develop in a helical arrangement from the top part of the fruit bunch rachis and are grouped by clusters (like fruit bunch whorls). The fruitlets develop from flowers along the spikelets. The lower sections of the spikelets, about 5 to 8 cm above the spikelet connection to the rachis, do not carry flowers or fruitlets (fig. 1).
The rachis and fruitlets can be modeled based on typical organ instances. However, due to their intricate, three-dimensional spatial arrangement, current agricultural image synthesis methods are unsuitable for modeling spikelets. The large variability in the 3D organization of spikelets is challenging to represent from samples of organ instances, and the spikelets do not resemble a multiple-link structure as modeled by L-systems. However, the spatial organization of each spikelet resembles a smooth 3D curve, and parametric equations can parsimoniously represent 3D curves. In computer graphics, Bezier curves are commonly used parametric curves. Bezier curves are specified by a set of discrete control points with a mathematical formulation defining a continuous curve through the control points (Singh, 1996). Bezier curves are used in many implementations, e.g., computer-aided design and computer-aided manufacturing (CAD/CAM) applications, graphic design programs, and graphic technologies (Pomax, 2022). The modeling method developed in the current work is based on structural decomposition into plant parts and 3D Bezier curve modeling and is suitable for capturing date fruit bunch characteristics. The usefulness of the synthesized models is demonstrated for analyzing Medjool date thinning.
 |
| Figure 1. Medjool date fruit bunch (early in thinning period). |
A fully grown Medjool palm has 20 to 30 fruit bunches. Each fruit bunch has 50 to 100 spikelets, and each spikelet carries 50 to 100 fruitlets. While some fruitlets shed naturally during fruit development, there is usually a significant excess of fruitlets per bunch, and thinning is required. The thinning aims to ensure a high-quality harvest where each fruit bunch yields 400 to 600 fruits. Ideally, thinning should be done by removing individual fruitlets along the spikelets with a shearing motion, allowing the remaining fruitlets to thrive. However, this method is typically not feasible due to time constraints since the entire thinning process must be accomplished within four weeks, starting four to five weeks after pollination. In practice, manual thinning removes the fruit bunch's inner (higher) spikelet cluster and shortens the remaining spikelets. In this process, the worker bundles the spikelets together with one hand and uses the other hand to shorten the spikelets (with pruning shears) near the bundling point. Repeating this operation 2 or 3 times may be necessary for large fruit bunches, bundling and cutting only some of the spikelets each time. The remaining length of the spikelets is manually determined based on multiple factors: the bunch whorl (top, middle, or bottom), the time of thinning, the distribution of the fruitlets on the spikelets, the expected natural fruitlet abscission, the girth of the rachis, the annual climate conditions, the tree's history, and the farmer’s long-term yield policy. Estimation algorithms suitable for manifold learning, e.g., deep neural networks (DNN) or kernel estimation (Cohen et al., 2021), can be applied to automate the estimation of the required remaining length of the spikelets. DNNs can scale well to high dimensions and handle correlations between output variables, but they require large training datasets to attain accurate representations. Gaussian kernels uniformly dispersed in the problem space can establish a smooth surface that depicts the mapping of the parameter space (Tamošiunaite et al., 2009). Kernel estimation facilitates a highly accurate representation that can be achieved with only a moderate training dataset size due to the ability to adapt kernel weights locally. On the downside, kernel estimation does not scale well to high dimensions and is not well suited to handling correlations between output variables. Hence, examining the parameter space before the estimation process is essential for kernel estimation. Examining the parameter space using dimensionally reduction methods like principal component analysis (PCA) can provide insights into the required structure of the kernels and the suitable network (Wold et al., 1987; Cohen et al., 2021). For example, when PCA indicates that features do not contribute to explaining the variability in the data, they can be excluded from the kernel equations.
The current paper presents a method for generating synthetic Medjool fruit bunch models and images. The presented work examines the usability of the synthesized images for training a kernel estimation model for automatically determining the required length of the spikelets (for leaving the required number of fruitlets on the bunch) after thinning. Preliminary work has been published in conference abstract form (Bar-Shira et al., 2019a, b).
Materials and Methods
The 3D Medjool Fruit Bunch Model
As Barth et al. (2018) suggested, the Medjool fruit bunch can be structurally decomposed into plant parts (fruit bunch organs), namely, the rachis, the spikelets, and the fruitlets. The fruitlets and rachis can be represented by shape primitives (basic geometric shapes). The fruitlets have an oval shape and are relatively small with respect to the fruit bunch size during the thinning period. Therefore, the fruitlets can be modeled as spheres. The rachis has an elongated elliptical shape that narrows towards its end, where the tip of the narrowed part is obscured from view by the spikelets. Rachis dimensions differ between fruit bunches. The rachis can be represented as an elliptic cylinder topped with a cone. Spikelets have an elongated, hose-like shape with a smooth, curved spatial formation. Thus, neither concatenating randomized basic shape meshes nor growing a series of links, as was done in previous works, is suitable for capturing the characteristics of the spikelets. An alternative modeling approach, suitable for spikelets, uses Bezier curves (Singh, 1996) to define their spatial shape as a deformable linear object. After the spikelet curve is determined, the elements of the model can be assembled to form realistic Medjool fruit bunches, i.e., the fruitlets can be attached to a spikelet, and spikelets can be attached to clusters along the rachis. In the current work, nine colors were defined for drawing the synthetic fruit bunch, three for each plant part (spikelets, rachis, fruitlets) at each thinning period (early, middle, late). The colors were empirically determined based on physical images of the respective plant parts and periods (Shoshan and Berman, 2021).
| Table 1. Medjool date fruit bunch model variable distributions. Distribution parameters are given here for the top whorl - early period and the middle whorl - middle period (the distribution parameters for the additional models are provided in the appendix). | ||||
| Index | Variable | Distribution function[a] | Distribution parameters | |
| Top whorl– early period | Middle whorl–middle period | |||
| p1 | Rachis width [cm] | U(a,b) | 2, 3 | 3, 4 |
| p2 | Rachis depth [cm] | U(a,b) · (p1) | 0.56, 0.76 | 0.56, 0.76 |
| p3 | Viewed length of the rachis [cm] | Constant | 20 | 20 |
| p4 | Hidden rachis length (between spikelets) [cm] | Triang(a, b, c) | 20, 25, 30 | 15, 23, 30 |
| p5 | Angle between rachis and trunk [°] | U(a, b) | 0, 10 | 15, 25 |
| p6 | Number of spikelets | Triang(a, b, c) | 80, 85, 100 | 70, 80, 100 |
| p7 | Number of spikelet clusters | Constant | 3 | 3 |
| p8 | Distance of clusters from top of rachis [cm] | Constant · (p4) | 0, 0.4, 0.6 | 0, 0.4, 0.6 |
| p9 | Dispersion of spikelet quantity between clusters | Constant · (p6) | 0.5, 0.3, 0.2 | 0.5, 0.3, 0.2 |
| p10 | Spikelet length [cm] | Triang(a, b, c) | 75, 80, 100 | 65, 79, 89 |
| p11 | Spikelet radius [cm] | U(a, b) | 2, 2.5 | 2, 2.5 |
| p12 | Max radial distance of spikelet end from rachis [cm] | Triang(a, b, c) | 10, 12, 15 | 40, 45, 50 |
| p13[b] | Number of fruitlets per spikelet | U(a, b) | 60, 70 | 30, 40 |
| p14 | Distance from spikelet base to first fruitlet [cm] | Triang(a, b, c) | 20, 27, 45 | 20, 27, 45 |
| p15 | Fruitlet radius [cm] | Triang(a, b, c) | 1, 1.2, 1.5 | 1, 1.2, 1.5 |
| p16[b] | Distance between two adjacent fruitlets [cm] | Exp(1/?) | 1.3 | 1.3 |
| p17 | Natural shedding probability [%] | Berr(p) | 16 | 44 |
| p18 | Bending of the fruit bunch [%] | Constant | 0 | 0.2 |
| p19 | 3rd control point radial distance (controls bending arc) | Constant | 1 | 0.85 |
[a]U (a,b ) – uniform distribution; Triang (a,b,c ) – triangular distribution; Berr (p ) – Bernoulli distribution; Exp (1/?) – exponential distribution. [b] Variables used for testing the generation process detailed below. | ||||
Random variables can be used to form a fruit bunch model suitable for capturing fruit bunch shape diversity within the general fruit bunch shape. Nineteen model variables were defined to capture the geometry of the fruit bunch (table 1). The model variables capture the characteristics of the rachis (width, thickness, length, bending, and orientation), the spikelets (quantity, dispersion, and length), and the fruitlets (quantity, dispersion, size, and natural shedding probability before the thinning). A distribution function was defined for each variable. The distribution parameters were adapted for three thinning periods (early, middle, and late relative to fruitlet development) and three whorls (top, middle, and bottom), leading to nine sets of values. Two examples of distribution parameters are given in table 1, and all the parameter sets are provided in the appendix. The distribution functions and their parameter values were determined based on interviews with plant experts and farmers and a visual comparison of the synthetic images with images of physical fruit bunches. Since the overall model is complex, the most straightforward distribution function that coincided with available knowledge for each feature was selected.
The spikelet clusters are modeled as deformable linear objects growing around the rachis (fig. 2), with the clusters beginning at three different distances from the top of the rachis depending on variable p8 (distance of clusters from top of rachis). For each cluster, spikelets are added around the rachis at intervals of 3°, where the probability of adding a spikelet at each interval is calculated according to the number of spikelets in the cluster, which is based on variables p6 (the number of spikelets) and p9(dispersion of spikelet quantity between clusters). For intervals of 3°, the maximal number of spikelets in each cluster is 120. According to the model, the expected number of spikelets in each cluster is defined by the vector S = p6*p9. Therefore, the probability of adding a spikelet at each interval of cluster i is S[i]/120. The length of the spikelets in the bottom (outer) cluster is determined based on variable p10 (spikelet length). The inner spikelets in the middle and top clusters are shorter than the spikelets in the bottom cluster. The model determines their length based on their higher initial growth location and the desired ending aligned with the spikelets from the bottom cluster, achieving the typical overall shape of the bunch that resembles a witch's broomstick.
 |
| Figure 2. Generating Medjool clusters based on the stochastic model. Left: Rachis and spikelets one from each cluster (Top, Middle, Bottom). Right: Rachis, spikelets, and fruitlets. The rachis is modeled as a cylinder and concentric cone. Spikelets are modeled as linear deformable objects with a spatial shape determined by Bezier curves. The spikelets are allocated to three clusters which start at three typical heights along the rachis cone (Bottom, Middle, Top) and are arranged about cone. Fruitlets are modeled as spheres connected around each spikelet. |
Variables p3 (viewed length of the rachis), p10(spikelet length), p12 (max radial distance of spikelet end from rachis), p18 (bending of the fruit bunch), and p19 (3rd control point radial distance) are used to define four Bezier points for a cubic Bezier curve, B(t), which define the spatial shape of each spikelet adapted change (eq. 1).
(1)
where V0, V1, V2, and V3 are Bezier points in 3-dimensional space, and Vi = [Vi,x, Vi,y, Vi,z]. The curve starts at V0, on the circumference of the rachis at the height appropriate for the cluster, and ends at V3. For all spikelets, V3,y, i.e., the height of V3, is defined by equation 2.
(2)
For the spikelets in the bottom (outer) cluster, the x and z coordinates of V3 depend on variable p12 (max radial distance of spikelet end from rachis) and on V0,j(j = x, z), such that the point V3 lies at the same angle as V0 and in a radius that equals p12. For spikelets in the higher (inner) clusters (middle, top), the x and z coordinates of V3 are similarly defined. However, the radial distance is adapted to adjust the overall form of the bunch. In the early period, fruit bunches in the top whorl are neatly arranged (the spikelets are not entangled, and almost all inner spikelets remain inner). For lower fruit bunches and as the season progresses, the order in the arrangement of the spikelets diminishes, and the spikelets become disheveled. For the top whorl - early period, the typical radius (of variable p12) was reduced by 34% for the middle cluster and by 60% for the top cluster. For the middle whorl - middle period, the typical radius was increased by 11% for the middle cluster but not changed for the top cluster. The Bezier control points, V1 and V2, determine the derivative at each end of the spikelet, respectively: V0 = 3(V1 -– V0), V3 = 3(V3 – V2). At their lower end, the spikelets grow along the direction of the rachis. As the season progresses, the spikelets start bending downward along their length due to the growing weight of the fruit. This typical shape is empirically attained based on the control points. To capture the disarrangement found in the field and to achieve a wide diversity of fruit bunch shapes, the coordinates of V1, V2, and V3 were increased randomly. In most cases, the added values were small (relative to the length of the spikelets), and in a few cases, the added values were large. For 70% of the spikelets, the coordinates were increased by up to 10%, with a 50% probability sampled separately for each coordinate (x, y, z). For the remaining 30% of the spikelets, the coordinates were increased by up to 60% with 50% probability sampled separately for each coordinate (x, y, z).
The fruitlets are drawn in a spiral around the Bezier curve, defining the shape of the spikelet. The fruitlets are attached to the points generated by the curve’s equation. The first fruitlet is attached at the first point whose distance from the starting point is larger than variable p14 (distance from the spikelet base to the first fruitlet). From that point until the end of the curve, fruitlets are drawn at each curve point with a probability according to variable p17(natural shedding probability). The number of Bezier curve points is selected such that the average number of fruitlets per spikelet (p13) and the distance between two adjacent fruitlets (p16) fit the model.
Generating Image Datasets for Learning the Required Length After Thinning
To demonstrate the usefulness of the constructed model, fruit bunch datasets for two whorl–period combinations were constructed based on the 3D model: the top whorl during the early period and the middle whorl during the middle period. A total of 1000 fruit bunches were sampled for each dataset, where each fruit bunch was sampled with a different random seed to create independent samples. The variable values for each fruit bunch were sampled from the model distribution and visualized using the open-source PyOpenGL library (PyOpenGL, n.d.). A screenshot of every generated fruit bunch was taken automatically, forming an image database. The image database for each of the two whorl–period combinations comprised 1000 images.
A thinning policy considering the available information at the thinning period, the future shedding of fruitlets, and targeting 400 to 500 fruits on each fruit bunch in the future at harvest was defined by two plant physiologists, both experts in Medjool cultivation. The thinning policy determined by the experts was to cut the spikelets along a thinning plane perpendicular to their growth direction such that an average of 12 fruitlets remain along the spikelets of the lower (outer) cluster. This thinning policy was determined for both whorl-period combinations. According to this policy, the remaining length of the spikelets after thinning was determined based on the simulated fruit bunch model and was visually marked with a disk placed orthogonally to the spikelets’ growth direction at the determined thinning length. This length was defined as the ground truth remaining spikelets length for training the estimation algorithms. The disk and the number of fruitlets on the bunch before and after thinning were indicated in images presented to the experts to facilitate verifying that the fruit bunches were correctly synthesized and that the thinning policy was suitable. The date experts visually verified the appearance of the fruit bunch and affirmed that the determined remaining length of the spikelet and the remaining number of fruitlets were appropriate. The synthetic fruit bunch images were manually placed within physical images to validate the model's similarity to physical images of Medjool fruit bunches and to enhance the visualization authenticity. This was done by rotating and scaling the synthetic bunch, applying a directed light filter based on the light in the physical image, and using a Gaussian blur filter on the rachis.
Examining the Usability of the Synthetic Images for Learning the Required Remaining Length of the Spikelet After Thinning
The required remaining spikelets length is influenced by time epoch, whorl, and the characteristics of the fruit bunch. These characteristics include the number of fruitlets on the fruit bunch and the length of the spikelets. However, since the spikelets are curved and the fruitlets do not start growing at the base of the spikelet, determining the required length based on a 2D image of the fruit bunch is challenging. During the season, the shape of the fruit bunch changes. The spikelets spread outward from the center of the bunch. The fruitlets become heavier, causing the spikelets to curve downward. These changes affect the area of the spikelet region and its angle. The size of the minimal oriented bounding box (OBB) about the spikelets is related to the size of the spikelet region, and the angle of the OBB is associated with the angle of the spikelets, both of which are critical components for determining the required length of the spikelets after thinning (Shoshan et al., 2022). Therefore, the size and angle of the minimal OBB were selected as parameters for estimating the required thinning length.
Synthetic 2D images of fruit bunches were made by taking screenshots of the generated fruit bunch models. Segmentation for extracting the fruit bunch from the screenshot was conducted based on color attributes, and a mask was created for each fruit bunch. The initial growth line on the rachis of the lowest spikelet cluster was found geometrically based on the segmented mask of the fruit bunch. Since the rachis is the lowest part of the fruit bunch, the growth line of the spikelets was found by examining the width of the segmented object mask starting from the lower part of the image until reaching a threshold value determined based on the maximal modeled width and depth of the rachis.
The minimal OBB was computed about the spikelets (excluding the rachis based on the initial growth line of the lowest spikelet cluster) using rotation matrices and the Feret diameter (Legland, 2014) at each angle between 0 and 180 degrees relative to the x-axis and choosing the OBB with the minimum area. The minimal OBB surrounding the spikelets was defined by its length (LS), width (WS), and angle (AS). The OBB features were extracted from the images and arranged in a database along with the fruitlet load (FL) and the required length for thinning (L) extracted directly from the fruit bunch model. The feature database for each of the two whorl–period combinations comprised 1000 feature tuples. Each database was subdivided into 800 samples for training and 200 samples for testing.
The two generated feature databases were used to examine problem complexity based on PCA and the estimation of L using kernel estimation. The estimation error of the required thinning length was computed as the scaled absolute error (SAE) between the output of the kernel estimation function and the ground truth value (eq. 3).
(3)
where yi is the true value, and y^i is the estimated value of sample i.
The results were compared using an ANOVA analysis with fruit bunch group (top whorl – early period, middle whorl – middle period) as the factor.
Results
Generating Image Datasets for Learning the Required Length After Thinning
For both datasets, generating Bezier curves with 250 points fit the model constraints. The generated fruit bunch images (fig. 3) and the required spikelet length after thinning (fig. 4) were validated by two plant physiologists, both in Medjool date cultivation. The experts reviewed and approved ten images of fruit bunches from each dataset.
 |
| Figure 3. Medjool fruit bunches Left: physical images; Middle: synthetic images: Right: synthetic images planted in physical images. Top: top whorl – early period; Bottom: middle whorl – middle period. |
Feature Extraction and Complexity Analysis
The results of the fruit bunch segmentation algorithm were manually verified for the synthetic images (fig. 5). The length (LS), width (WS), and angle (AS) of the OBB were successfully extracted for all the images in the two synthetic datasets (2000 images overall).
All five PCA components were required to explain more than 90% of the variance in the top whorl – early period (fig. 6, left). Four PCA components explained more than 90% of the variance in the middle whorl – middle period. In both whorl-period combinations, the cutting length, L, was related to several combinations of all input features (fig. 6, right). This finding suggests that the input features may be appropriate for predicting L and that dimensionality reduction is not straightforward using a linear model. For the top whorl – early period, the first component comprises large weights for WS, FL, and L, and the second component comprises large weights for AS, WS, and LS. Both components contribute similarly to explaining the overall variance. For the middle whorl – middle period, the first component comprises large weights for AS, WS, and LS, and the second component comprises large weights for FL and L and medium weights for WS and LS. The first component contributes nearly twice as much as the second component to explaining the overall variance. Based on this analysis and since FL is difficult to obtain in the field, AS, WS, and LS were retained for kernel estimation of L.
The Required Length of the Spikelets After Thinning
Three Gaussian kernels were uniformly distributed in each of the three feature axes (AS, WS, and LS), between the minimal and maximal feature values in the training database, leading to 27 Gaussian kernels in a three-dimensional space. The kernel weights and standard deviations were determined by minimizing the average square error over the training dataset (Shinamazaki and Shinomoto, 2010). The estimation was tested with the testing dataset.
The estimation errors of the required remaining spikelet length after thinning for the test set (200 images) are presented in figure 7. The average scaled absolute errors for each period were lower than 1% (Top-Early: 0.95, 0.69 SD, Middle-Middle: 0.88, 0.69 SD), with a very small number of outliers (8). The errors were similar in the two whorl-period groups.
Discussion
Modeling based on a structural decomposition of the Medjool fruit bunch into rachis, spikelets, and fruitlets and Bezier curves for modeling the spatial shape of the spikelets resulted in realistic models of Medjool fruit bunches requiring relatively low computation power. As found in the field, the fruit bunches in the top whorl – early period had a somewhat standardized, elongated appearance, while the fruit bunches in the middle whorl – middle period were more disheveled. Date experts verified that the model was suitable for generating fruit bunch images and establishing the required thinning length.
 |
| Figure 4. Medjool fruit bunches with required remaining length of spikelets marked for expert review. Top: top whorl – early period. Bottom: middle whorl – middle period. Numbers of fruitlets before and after thinning are given below each image. |
 |
| Figure 5. Segmentation and minimal oriented bounding box. Left: top whorl – early period. Right: middle whorl – middle period. |
The features of the OBB extracted from the images facilitated training a kernel estimator for determining the required remaining spikelet length after thinning. The average absolute scaled errors were lower than 1% for both whorl–period groups. For physical Medjool fruit bunches, Shoshan et al. (2021) attained absolute scaled errors below 20% for spikelet remaining length estimation. The estimation was accomplished using an adapted form of the VGG-16 convolutional neural network (Simonyan and Zisserman, 2014), pre-trained based on the ‘ImageNet’ dataset (Deng et al., 2009). The estimation errors with the physical images are much higher than in the current work. This difference can be due to obstructions (by leaves and other fruit bunches) and lighting issues in the images of physical fruit bunches that are not sufficiently represented in the current model. Visualization authenticity of the synthetic images can be enhanced by embedding the generated models within physical images. Such an embedding can be automated, facilitating the use of synthetic images for enhancing annotated databases of physical images.
The fruitlet load is readily available through the synthetic fruit bunch model. However, it is challenging to estimate the fruitlet load in the field. Due to their small size and color similarity to other plant parts, the fruitlets are very hard to visually identify. While the fruitlet load is related to the area of the fruit bunch through spikelets length and number, this relation is not linear due to spikelet bending, fruitlet fallout, and since the lower sections of the spikelets do not carry fruitlets. The fruitlet load may be estimated more accurately in realistic settings using additional sensors, e.g., sonar sensors (Finkelstein et al., 2017) or hyperspectral vision systems (Yang et al., 2008).
Conclusions
The stochastic Medjool fruit bunch 3D model developed in the current work facilitated the generation of multiple Medjool fruit bunch models and images. Date cultivars vary greatly in terms of structural properties and cultivation process requirements. The current work specifically targets Medjool dates due to their premium value in the date market. However, the developed synthesis and analysis methods could be adapted to fit other date cultivars.
Integration of model variables with smooth Bezier curves facilitated a parsimonious representation of major plant features. An alternative representation based on the concatenation of small spikelet sections, as suggested by Barth et al. (2018) or Benoit et al. (2014), would have resulted in a less accurate and more complex model.
  |
| Figure 6. Principal component analysis. Top: top whorl – early period. Bottom: middle whorl – middle period. Right: cumulative variance explained by principal components. Left: weights of original features in each of principal components. L - cutting length, AS - spikelet bounding box angle, WS - spikelet bounding box width, LS - spikelet bounding box length, FL - fruit load. |
Synthesized fruit bunch images may be integrated within physical date palm images. Alternatively, future work can synthesize entire Medjool palm trees, by modeling additional plant organs, i.e., fonds and trunk. Realistic palm images can then be rendered by integrating the developed PyOpenGL modules with 3D graphics software, such as Blender (Blender Foundation, 2002). The dissimilarity between synthetic and physical images can be further reduced using advanced image processing methods, e.g., by applying adversarial networks (Barth et al., 2020). Such augmented synthetic images of date palms can be used to enhance annotated image databases of physical fruit bunches and date palms. The improved databases can then be used to train advanced motion planning algorithms required for developing date cultivation robots.
 |
| Figure 7. Box plots of average scaled absolute error of estimated required remaining length of spikelets after thinning. Square represents interquartile range (IQR), middle line of box represents median, whiskers represent variability outside upper and lower quartiles, and dots represent outliers (values above Q3 plus 1.5 times IQR or below Q1 minus 1.5 times IQR). |
Appendix A – Distribution Parameters
The distribution parameters for the top whorl – early period and the middle whorl – middle period are given in table 1. The distribution parameters of all nine whole-period combinations are given in table A1.
| Table A1. Medjool date fruit bunch distribution parameters. | ||||||
| Index | Parameter | Probability Function[a] | Whorl | Epoch | ||
| 1st | 2nd | 3rd | ||||
| p1 | Rachis width [cm] | U(a, b) | 1st | (2, 3) | (2, 3) | (3, 4) |
| 2nd | (3, 4) | (3, 4) | (4, 5) | |||
| 3rd | (4, 5) | (4, 5) | (5, 6) | |||
| p2 | Rachis depth [cm] | U(a, b) · (p1) | 1st | (0.56, 0.76) × (p1) | ||
| 2nd | ||||||
| 3rd | ||||||
| p3 | Viewed section of rachis [cm] | Constant | 1st | 20 | ||
| 2nd | ||||||
| 3rd | ||||||
| p4 | Hidden Rachis length [cm] | Triang(a, b, c) | 1st | (20, 25, 30) | (22, 25, 32) | (29, 35, 40) |
| 2nd | (17, 19, 25) | (15, 23, 30) | (25, 32, 37) | |||
| 3rd | (10, 13, 15) | (14, 20, 24) | (22, 25, 30) | |||
| p5 | Angle between rachis to tree trunk [°] | U(a, b) | 1st | (1, 10) | (10, 20) | (15, 30) |
| 2nd | (10, 17) | (15, 25) | (45, 60) | |||
| 3rd | (10, 20) | (30, 40) | (60, 90) | |||
| p6 | Number of spikelets | Triang(a, b, c) | 1st | (80, 85, 100) | ||
| 2nd | (70, 80, 100) | |||||
| 3rd | (50, 60, 65) | |||||
| p7 | Number of spikelets clusters | Constant | 1st | 3 | ||
| 2nd | ||||||
| 3rd | ||||||
| p8 | Distance of spikelets clusters from top of presented rachis | Constant · (p4) | 1st | Bottom cluster: 0 × (p4) Middle cluster: 0.4 × (p4) Top cluster: 0.6 × (p4) | ||
| 2nd | ||||||
| 3rd | ||||||
| p9 | Dispersion of spikelet quantity between the clusters | Constant · (p6) | 1st | Bottom cluster: 0.5 × (p6) Middle cluster: 0.3 × (p6) Top cluster: 0.2 × (p6) | ||
| 2nd | ||||||
| 3rd | ||||||
| p10 | Spikelet length [cm] | Triang(a, b, c) | 1st | (75, 80, 100) | ||
| 2nd | (65, 79, 89) | |||||
| 3rd | (50, 60, 65) | |||||
| p11 | Spikelet radius [cm] | U(a, b) | 1st | (2, 2.5) | ||
| 2nd | ||||||
| 3rd | ||||||
| p12 | Max radial distance of spikelet end from rachis [cm] | Triang(a, b, c) | 1st | (10, 12, 15) | (35, 40, 45) | (70, 80, 100) |
| 2nd | (20, 24, 30) | (40, 45, 50) | (50, 70, 90) | |||
| 3rd | (20, 21, 24) | (24, 30, 35) | (50, 70, 85) | |||
| p13 | Number of fruitlets per spikelet | U(a, b) | 1st | (60, 70) | (40, 50) | (25, 35) |
| 2nd | (50, 60) | (30, 40) | (17, 25) | |||
| 3rd | (35, 45) | (25, 35) | (12, 22) | |||
| p14 | Distance from the spikelet base to the first fruitlet [cm] | Triang(a, b, c) | 1st | (20, 27, 45) | ||
| 2nd | ||||||
| 3rd | ||||||
| p15 | Fruitlet radius [cm] | Triang(a, b, c) | 1st | (1, 1.2, 1.5) | (1.5, 1.8, 2) | (1.8, 1.8, 2.5) |
| 2nd | (0.9, 1.1, 1.2) | (1, 1.2, 1.5) | (1.5, 1.8, 2) | |||
| 3rd | (0.7, 0.8, 0.9) | (0.9, 1.1, 1.2) | (1, 1.2, 1.5) | |||
| p16 | Distance between two adjacent fruitlets [cm] | Exp(1/?) | 1st | (1.3) | ||
| 2nd | ||||||
| 3rd | ||||||
| p17 | Natural fallout probability [%] | Berr(p) | 1st | 16 | 44 | 60 |
| 2nd | ||||||
| 3rd | ||||||
| p18 | Bending of the fruit bunch [%] | Constant | 1st | 0 | 0.2 | 0.2 |
| 2nd | ||||||
| 3rd | ||||||
| p19 | 3rd control point radial distance (controls bending arc) | Constant | 1st | 1 | 0.85 | 0.70 |
| 2nd | ||||||
| 3rd | ||||||
[a]U(a, b) – uniform distribution; Triang(a, b, c) – triangular distribution; Berr(p) – Bernoulli distribution; Exp(1/?) exponential distribution. | ||||||
References
Bar-Shira, O., Cohen, Y., Shaubi, T., Sadowky, A., Cohen, Y., Bechar, A., & Berman, S. (2019a). Learning motion parameters for a Mejdool-date thinning robot. Proc. XXXVIII CIOSTA & CIGR Section V Conf.
Bar-Shira, O., Shoshan, T., Cohen, Y., Sadowky, A., Cohen, Y., Bechar, A., & Berman, S. (2019b). Recognition and motion planning for a Mejdool-date thinning robot. Proc. The Israeli Robotic Conf. (ICR).
Barth, R., Hemming, J., & Van Henten, E. J. (2020). Optimising realism of synthetic images using cycle generative adversarial networks for improved part segmentation. Comput. Electron. Agric., 173, 105378. https://doi.org/10.1016/j.compag.2020.105378
Barth, R., Ijsselmuiden, J., Hemming, J., & Van Henten, E. J. (2018). Data synthesis methods for semantic segmentation in agriculture: A Capsicum annuum dataset. Comput. Electron. Agric., 144, 284-296. https://doi.org/10.1016/j.compag.2017.12.001
Benoit, L., Rousseau, D., Belin, É., Demilly, D., & Chapeau-Blondeau, F. (2014). Simulation of image acquisition in machine vision dedicated to seedling elongation to validate image processing root segmentation algorithms. Comput. Electron. Agric., 104, 84-92. https://doi.org/10.1016/j.compag.2014.04.001
Blender Foundation. (2002) “blender.org - Home of the Blender Project - Free and Open 3D Creation Software.” blender.org, www.blender.org. Accessed 8 Jan. 2023.
Chao, C. T., & Krueger, R. R. (2007). The date palm (Phoenix dactylifera L.): Overview of biology, uses, and cultivation. HortScience, 42(5), 1077-1082. https://doi.org/10.21273/hortsci.42.5.1077
Cohen, Y., & Glasner, B. (2015). Date palm status and perspective in Israel. In J. M. Al-Khayri, S. M. Jain, & D. V. Johnson (Eds.), Date palm genetic resources and utilization: Volume 2: Asia and Europe (pp. 265-298). Dordrecht: Springer. https://doi.org/10.1007/978-94-017-9707-8_8
Cohen, Y., Bar-Shira, O., & Berman, S. (2021). Motion adaptation based on learning the manifold of task and dynamic movement primitive parameters. Robotica, 39(7), 1299-1315. https://doi.org/10.1017/S0263574720001186
Deng, J., Dong, W., Socher, R., Li, L. J., Li, K., & Fei-Fei, L. (2009). ImageNet: A large-scale hierarchical image database. Proc. 2009 IEEE Conference on Computer Vision and Pattern Recognition, (pp. 248-255). https://doi.org/10.1109/CVPR.2009.5206848
FAOSTAT. (2022). Crops and livestock products. Rome, Italy: United Nations FAO. Retrieved from http://www.fao.org/faostat/en/#data/QCL/visualize
Finkelshtain, R., Bechar, A., Yovel, Y., & Kósa, G. (2017). Investigation and analysis of an ultrasonic sensor for specific yield assessment and greenhouse features identification. Precis. Agric., 18(6), 916-931. https://doi.org/10.1007/s11119-016-9479-0
Legland, D. (2014). Feret diameter and oriented box v1.1. Retrieved from http://www.mathworks.com/matlabcentral/fileexchange/?30402-feret-diameter-andoriented-box
Leitner, D., Klepsch, S., Bodner, G., & Schnepf, A. (2010). A dynamic root system growth model based on L-Systems. Plant Soil, 332(1), 177-192. https://doi.org/10.1007/s11104-010-0284-7
Lu, Y., & Young, S. (2020). A survey of public datasets for computer vision tasks in precision agriculture. Comput. Electron. Agric., 178, 105760. https://doi.org/10.1016/j.compag.2020.105760
Pomax. (2022, September 3). A primer on Bézier curves. Pomax.github.io. Retrieved January 8, 2023, from https://pomax.github.io/bezierinfo/#introduction
PyOpenGL. PyPI. (n.d.). Retrieved January 8, 2023, from https://pypi.org/project/PyOpenGL/
Shimazaki, H., & Shinomoto, S. (2010). Kernel bandwidth optimization in spike rate estimation. J. Comput. Neurosci., 29(1), 171-182. https://doi.org/10.1007/s10827-009-0180-4
Shoshan, T., & Berman, S. (2021). Medjoul-date images - annotated, Mendeley Data, V1. https://doi.org/10.17632/k7xk2nwgrh.1
Shoshan, T., Bechar, A., Cohen, Y., Sadowsky, A., & Berman, S. (2022). Segmentation and motion parameter estimation for robotic Medjoul-date thinning. Precis. Agric., 23(2), 514-537. https://doi.org/10.1007/s11119-021-09847-2
Simonyan, K., & Zisserman, A. (2014). Very deep convolutional networks for large-scale image recognition. arXiv preprint arXiv:1409.1556. https://doi.org/10.48550/arXiv.1409.1556
Singh, N. (1996). Systems approach to computer-integrated design and manufacturing. NY: John Wiley & Sons.
Tamosiunaite, M., Asfour, T., & Wörgötter, F. (2009). Learning to reach by reinforcement learning using a receptive field based function approximation approach with continuous actions. Biol. Cybern., 100(3), 249-260. https://doi.org/10.1007/s00422-009-0295-8
Wold, S., Esbensen, K., & Geladi, P. (1987). Principal component analysis. Chemom. Intell. Lab. Syst., 2(1), 37-52. https://doi.org/10.1016/0169-7439(87)80084-9
Yang, C., Everitt, J. H., & Bradford, J. M. (2008). Yield estimation from hyperspectral imagery using spectral angle mapper (SAM). Trans. ASABE, 51(2), 729-737. https://doi.org/10.13031/2013.24370