Leaf shape simulation of castor bean and its application in nondestructive leaf area estimation

The leaf shape and leaf area measurement are crucial in plant growth modeling. The castor bean leaf is large, palm-shaped with multiple clefts. The leaf shape simulation and leaf area estimation were less studied. The circular model and nonrectangular hyperbolic model were developed to describe the standard leaf shape of castor bean in this study, providing a model for simulating the leaf shape and a nondestructive way for estimating the leaf area respectively. In addition, a formula was established to estimate the leaf area by the parameter of the standard leaf shape of castor bean. Based on validation results, the circular model fits the landmarks and nonrectangular hyperbolic model fits the lobe margins very well. The leaf area was accurately estimated by using the established formula. This study could provide a theoretical reference for leaf visualization, a nondestructive and easy way to estimate the leaf area for other complex leaves with multiple lobes.


Introduction
The leaf is an important organ of plants.Plant leaves participate in many physiological processes, such as transpiration, photosynthesis and so on.
The leaf area is an important determinant of light interception and consequently of transpiration, photosynthesis and plant productivity [1] .Thompson [2] stated that the area of an object was proportional to the square of its length (the area-length allometry).O'Shea et al. [3] studied the allometric relationship between the surface and the length of six species of fish.The estimates of the exponent of the first allometry ranged from 1.88 to 2.22.Montgomery [4] put forward a leaf area formula for corn: the leaf area (A) = k × leaf length (L) × leaf width (W).Where k is a constant to be fitted; The leaf length is defined as the distance from leaf apex to leaf base; The leaf width is defined as the maximum length of the segments perpendicular to the straight line passing through leaf apex and leaf base [5] .The Montgomery model has been demonstrated to be appropriate for calculating leaf area for other herbaceous and woody plants such as cotton, rice and sorghum [6][7][8] .These models would have many advantages, such as an expression in dimensional terms of length and width of a leaf, which would be easy to measure many plants with simple shape leaf, but not to easy measure other plants with complex shape leaf, such as caster bean leaf.
Recently, many studies have proposed parametric models to describe leaf shapes of some plants.These models can be used to calculate leaf area.For instance, Shi et al. [9] used the elliptical equation, the sectorial equation and the simplified Gielis equation to describe the leaves of Hydrocotyle vulgaris, Ginkgo biloba L., Polygonum perfoliatum L. and Oligostachyum sulcatum, and used these mathematical equations to derive the theoretical value of the Montgomery parameter.Dornbusch et al. [10] put forward a general parameter model to depict the leaf shapes of wheat, barley and maize.Gielis [11] proposed a superformula that could be used to describe the shapes of many abiotic and biotic objects, which has been simplified to fit to leaf dimensions.The simplified Gielis equation (SGE) has only two parameters.One is related to leaf length, and the other is related to the ratio of leaf width to length.Shi et al. [12,13] used the SGE to fit the leaf shapes of plants such as Indocalamus sp., Aucuba japonica var.variegata Dombrain, Chimonanthus praecox (L.) Link, Parrotia subaequalis (H.T. Chang) R. M. Hao & H. T. Wei, Phoebe sheareri (Hemsl.)Gamble and Pittosporum tobira (Thunberg) W. T. Aiton.However, these mathematical models would be not appropriate for complex caster bean leaf.
Castor bean (Ricinuscommunis L.) belongs to the family Euphorbiaceae and the genus Ricinus.Although, the castor bean was planted at the industrial scale in more than 30 countries, the leaf shape simulation and LA estimation were less studied.Jani and Misra [14] stated Montgomery model to estimate the area of Castor bean leaf, however, it would be not easy to know well the measure of the length and width of the leaf with different number of clefts.The castor bean leaf is large, palm-shaped with multiple clefts, which would not be described well by above mathematical model.In this study, a simulation model was constructed to describe geometric shape of castor bean leaf, and a formula of leaf area was derived by the model parameters.In addition, this model would be applied in the visualization of castor bean leaf and the identification of castor bean.

Plant materials and leaf description
The plant materials were collected from a garden at Wulong community of Beisheng in Liuyang, Hunan Province.They are located in a subtropical monsoon humid climatic zone (28.25°N, 113.41°W; the elevation is 65 m; the mean annual temperature is 17°C; the mean annual precipitation is 1400 mm).The soil in the garden is red-yellow mud paddy soil developed from Quaternary red clay.The contents of soil organic matter, total nitrogen, total phosphorus, total potassium, available nitrogen, available phosphorus and available potassium were 33.51 g/kg, 1.52 g/kg, 0.94 g/kg, 12.68 g/kg, 130.12 mg/kg, 30.78 mg/kg and 134.71 mg/kg respectively.Four varieties of castor beans were planted on the ridges of the garden with row spacing and plant spacing of 0.8 m which is the recommended planting density.
300 kg/hm 2 phosphorus diamine was used as a basal fertilizer in total during the plowing (i.e., before the planting).Based on the soil moisture content, the irrigation was used to prevent the plants from experiencing drought stress, and the weeding was performed regularly.
For convenience, the geometric parameters determined in the study are shown in Figure 1.In Figure 1, α i is the intersection angle between adjacent leaf veins; β i is the intersection angle between adjacent leaf veins and leaf valley lines.
The serial numbers refer to leaves counterclockwise from the main leaf vein.The leaf edge (green line) is the entire boundary of the leaf.The longest vein was denoted by blue line, and the joint point of all veins was denoted by the leaf center (blue dot).The leaf apexes (red dots) are determined from the joint point of the leaf vein and the leaf edge where N is the number of lobes.The leaf valley endpoint (yellow dots) is the dividing point of the lobes which is the nearest point from the lobe margin to the leaf center.The vein lines (pink lines) are drawn from the leaf apexes to the center.These lines are usually referred to as veins.The leaf valley lines (blue dashed lines) connect the leaf valley endpoints to the center.The lobe margin refers to the leaf boundary connecting a leaf apex and a leaf valley endpoint.Thus, each lobe has two edges.

Sample collection and data processing 2.2.1 Sample collection
Leaves were sampled randomly 45 d after branching in 2018.The samples include 333 leaves of four varieties (Zibi5, Yunbi2, Yunbi4, and Xiangbi1).Two hundred leaves were classified into 18 sample groups by variety and lobe number (Zibi5:6-11 lobes, Yunbi2:6-9 lobes, Yunbi4:7-10 lobes and Xiangbi1:6-9 lobes).As some groups (Zibi5:6 and 11 lobes, Yunbi2:6 lobes, Yunbi4:10 lobes, and Xiangbi1:6 lobes) had few leaves, 12 sample groups (four varieties: Zibi5, Yunbi2, Yunbi4, and Xiangbi1; three numbers of lobes: 7-9) were used in this study.In Table 1 the sample information was shown.Photos of samples can be found in the supplementary information (S1_Data).The sample group Zibi5 was used for the model development, discussion and application, and others were used for discussion and application partly.

Data processing
Samples were pressed and photographed against a background (1280×1024 pixel resolution) with the instrument described [15] .Leaf binary images were obtained after differential operation of the image and background.The leaf edge was extracted by using the boundary-tracking algorithm [16] .And the additional geometric parameters, such as leaf apexes and valley endpoints, edges of lobes and leaf area (LA), were extracted based on the leaf edge.The detailed steps are as follows: Step 1. Pixel coordinates of the leaf center were marked manually.
Step 2. Differential operation of the image and background [17] was carried out by using Equation (1), and the result was recorded as LIMd where LIM and BIM indicate the RGB data [18] of the leaf image and background respectively.
LIMd(i, j, k)=LIM(i, j, k)-BIM(i, j, k) (1) When k=1, 2,3, the values of LIM(i, j, k) and BIM(i, j, k) are the red, green and blue channels of the leaf image and background at pixel (i, j) respectively.
The image binar ization was performed by using Equation (2): where, LIMb is the binarized image and TH is the threshold adjusted to obtain a high-quality binarized image.
Step 3. Objects were extracted [19] by labeling the 8-connected objects, the leaf shape image was taken as the largest-area object, and the expansion operation and corrosion operation [20] were performed with the structure element of around plate with a radius of 3 pixels to smooth the leaf image.
Step 4. Leaf margin was obtained with boundary tracking [21] .Then, the pixel coordinates of leaf apexes and valley endpoints were determined by finding peaks with a minimum separation of 45 pixels [22] .The edges of the lobe were separated by leaf apexes and valley endpoints.After image holes were filled [23] , the leaf pixel area was calculated.
Step 5. Pixel coordinates of the leaf center, landmarks and point on the leaf margin were transformed to standard coordinates in a real coordinate frame (standard coordinate frame).The coordinate frame was shifted to ensure that its origin was overlapped with the leaf center.It was also rotated to ensure that its x-axis was overlapped with the longest vein and its unit was normalized to the length of the longest vein by transformation, translation and scaling.

Models and validation 2.3.1 Circular curve for landmark distribution
The distribution of landmarks (leaf apexes and valley endpoints) determining the profile of leaves is an important factor affecting the area of a palm-shaped leaf.As the geometric characteristics of castor bean leaves were nearly circular (Figure 1b), the distribution of leaf apex and valley positions was simulated based on a circular curve (Equation ( 3)) as follows: (x -x 0 ) 2 -(y -y 0 ) 2 = r 2 (3) where, (x 0 , y 0 ) is the circular center; r is the radius.

Lobe margin model
The lobe margin model describing the curves of half of the lobe margin in the coordinate systems was demonstrated in Figure 1b.In this model, the distance of a point on the lobe margin to a leaf vein was defined as a function of the distance from the leaf apexes to the projection of the lobe margin point on the leaf vein (Figure 1b).
In Equations ( 4) and ( 5), the l and H represent the length of the vector projected from the point to leaf apexes on the leaf vein and the distance of the point on the lobe margin from the leaf vein respectively.The light response curve of the photosynthetic rate was obtained by Equation ( 4) where , H m , and θ are the initial slope of the corresponding curvature, asymptotically final value and curvature respectively.Equation ( 5) was modified from a logistic model in which H 0 , H m and r represent the initial value of the "S-type" curve, final value and parameter depicting the rate of the curve change respectively.The term -H 0 is introduced to ensure that the leaf apex is located above the veins, that is, H(0)=0.

Simulation of standard leaf shape
The parameters of each leaf in Equation (3) were obtained by fitting the coordinates of leaf apexes and valley endpoints, and the parameters of each leaf in Equations ( 4) and ( 5) were fitted based on the coordinates of the edges of each lobe.The standard leaf shape was simulated as follows: The landmark distribution was simulated by the circular model with the mean parameters and the mean angle between the adjacent leaf vein line and valley line (Table 2); Then, the lobe margin model was run by Equations ( 4) and ( 5) with the mean parameters.The result of standardized leaf simulation was shown in Figure 2.

Leaf area calculation
The area (S) of each standard leaf was calculated by integrating the polar coordinates.Thus, the LA was calculated by using Equation (6). 6) where, L is the vein length; r is the radius of the leaf; s is the area of standard leaf and k is the regression coefficient of s×(r×L) 2 and actual leaf area.

Validation
The calculations in this study were carried out in MATLAB R2016a, and the performance of the model simulation was determined by the normalized root mean square error (RMSEn): When the RMSEn value is less than 10%, two values are desirably consistent.When the RMSEn values are 10%-20%, 20%-30% and greater than 30%, a good simulation, an acceptable simulation and a poor simulation were indicated respectively [26] .

Parameters of landmark distribution
The parameters (circular center and radius) in Equation ( 3) were determined by fitting the standard coordinates of landmarks for each leaf with 7-9 lobes and using nonlinear regression [27] .Based on the results (Table 2), the y-coordinates of the circular center of leaf apexes are small, i.e., the circular center is roughly located on the main leaf vein.In addition, for the circle of leaf apexes, the x-coordinates of the circular center in the coordinate frame are 0.21, 0.21, and 0.22, respectively, and the radii are 0.79, 0.78, and 0.77 for 7-9 lobes respectively.In other words, the center of the circle moves slightly away from the leaf center and the ratio of the radius to the leaf length slightly increases with the increase of lobe number.For the circle of leaf valley points, the x-coordinates of the circular center are 0.10, 0.10, 0.12, respectively, and the radius were 0.35, 0.39, and 0.41, respectively.In other words, the center of the circle moves slightly away from the leaf center and the ratio of the radius to the leaf length slightly increases with the increase of lobe number.Moreover, the result is the same for the leaf valley points.The RMSE of circumferential curve fitting is less than 10% for both leaf apexes and valley endpoints.

Distribution of the angles between leaf veins and valley lines
As shown in Figure 3, the angles between leaf veins and leaf valley lines were bilaterally symmetrical with the longest vein.The angles increased and the vein deviated from the longest vein for a leaf with 7 lobes; The angles increased for a leaf with 8 lobes although the last angle decreased for a leaf with 8 lobes; The angles increased with the increase of lobe number to 8 lobes and then decreased for a leaf with 9 lobes.

Parameters of the lobe margin function
The l i and h i represent length of the vector projected from the point to leaf apexes on the leaf vein and the distance of the point on the lobe margin from the leaf vein respectively.The lobe margin model was established by fitting half of the lobe margin based on the normalized data: l i = l i /max(l i ), h i = h i /max(h i ).The fitting effect for a leaf with 8 lobes was shown in Figure 4, and the parameters and validation were shown in Table 3.
Based on the RMSEn in Table 3, the fittings of the improved logistics model and the nonrectangular hyperbolic model are the same (both less than 10%).Thus, the models performed equally well.
Furthermore, the variation coefficients of the nonrectangular hyperbolic model parameters are lower than those of the improved logistic model parameters.It is indicated that the fitting of the nonrectangular hyperbolic model is better than that of the improved logistic model.The parameters of the nonrectangular hyperbolic model, including , H m , and θ, i.e., the initial slope, final value and curvature, were estimated to be 0.79, 0.51 and 0.96 with 95% confidence intervals (CIs) of 0.74-0.83,0.45-0.57and 0.89-1.03for Zibi5 with 7 lobes, 0.79, 0.49 and 0.97 with 95% CIs of 0.75-0.82,0.46-0.52 and 0.95-1.00for Zibi5 with 8 lobes, and 0.76, 0.50 and 0.97 with 95% CIs of 0.71-0.81,0.47-0.53and 0.95-0.99 for Zibi5 with 9 lobes respectively.

Leaf area calculation
As shown in Table 4, the formula for calculating the leaf area of different lobes was established.The area estimation based on main leaf vein length and shape parameters has a lower RMSEn and more stable coefficients than that based only on main leaf vein length.Thus, the formula of leaf area calculation based on main leaf vein length and shape parameters is better than that based only on main leaf vein length.

Application and significance
The castor bean leaf has a palm-like shape and a piecewise function, which may be a good option to depict the shape.An improved logistic model and a nonrectangular hyperbolic model were used to fit the leaf boundary, which showed great accuracy in fitting the lobe margin.Furthermore, the circular radius is almost in the interval [0.75, 0.80] for the leaf apexes and the interval [0.35, 0.41] for the leaf valley points (Table 5).The radius slightly increased with the increase of lobe number for all sample groups, except the leaf apex of Yunbi2 and the leaf valley point of Xiangbi1.In summary, the results of this study may have theoretical significance for the visualization of castor bean leaves.Nondestructive and accurate methods of LA estimation are significant, especially under the field conditions [28] .In this study, the empirical formula with leaf area, standard leaf shape, parameters of leaf radius and the number of lobes was established for LA calculation (where r is the radius of the standard leaf shape and L is the length of the leaf).Based on the general results, the LA of castor bean can be estimated more accurately from the leaf radius.

Conclusions
To describe the leaf shape of castor bean, a standard simulation model of castor bean leaf shape was developed by using statistical and mathematical modeling methods.The model includes sub models, i.e., a circular model describing the landmark distribution, and functions used to depict lobe margins.The simulated standardized leaf is similar in morphology to areal castor bean leaf.This study would provide a new method for morphological simulation and leaf area of complex shape leaf with lobes.It would also have important reference significance for other crops digital research and phenotype model construction.

Figure 1
Figure 1 Description of Ricinuscommunis L. leaf used in this study

Figure 2
Figure 2 Standardized shape of a blade with 7-9 lobes a. 14 angles for leaves with 7 lobes b. 16 angles for leaves with 8 lobes c. 18 angles for leaves with 9 lobes Note: The angles are bilaterally symmetric.

Figure 3
Figure 3 Angles between adjacent leaf veins and leaf valley lines

Figure 4
Figure 4 Lobe margin fitting for 8 lobes

Table 2 Parameters of circles fitted for leaf apexes and leaf valley points
Note: LN: the lobe number of a cleft leaf; RMSE: the root mean square error; Var: the variance.