Growth Analysis of Wheat Using Machine Vision: Opportunities and Challenges

Crop growth analysis is used for the assessment of crop yield potential and stress tolerance. Capturing continuous plant growth has been a goal since the early 20th century; however, this requires a large number of replicates and multiple destructive measurements. The use of machine vision techniques holds promise as a fast, reliable, and non-destructive method to analyze crop growth based on surrogates for plant traits and growth parameters. We used machine vision to infer plant size along with destructive measurements at multiple time points to analyze growth parameters of spring wheat genotypes. We measured side-projected area by machine vision and RGB imaging. Three traits, i.e., biomass (BIO), leaf dry weight (LDW), and leaf area (LA), were measured using low-throughput techniques. However, RGB imaging was used to produce side projected area (SPA) as the high throughput trait. Significant effects of time point and genotype on BIO, LDW, LA, and SPA were observed. SPA was a robust predictor of leaf area, leaf dry weight, and biomass. Relative growth rate estimated using SPA was a robust predictor of the relative growth rate measured using biomass and leaf dry weight. Large numbers of entries can be assessed by this method for genetic mapping projects to produce a continuous growth curve with fewer replicates.


Introduction
Crop growth analyses have been integral to crop physiology over the last century and entailed calculation of net assimilation rates [1,2]. Blackman [3] introduced the concept of accumulation of dry weight, later interpreted as rate of change in weight or relative growth rate (RGR). Blackman [3] also introduced the efficiency index of dry weight production and drafted the compound interest law as a simple equation that best describes the growth of annual plants: (1) where W 1 = the final weight, W 0 = the initial weight, r = the efficiency index of dry weight production, t = time, and e is the base of natural logarithms. Since then, crop physiologists attempted to explain growth quantitatively using agricultural inputs such as fertilizer supply, environmental conditions, and periodic measurements of leaf area or dry weight as explanatory variables or responses at intervals throughout the plant growth from sowing to maturity. express growth traits was side projected area (SPA). One of the key steps in applying the RGB images for plant phenotyping is segmentation, which is to separate the targets of interest from the background [23]. Most of the current high-throughput plant phenotyping facilities leverage color index-based technologies, which establish segmentations by finding the difference in the color between the targets of interest and the background [24][25][26][27][28]. This method could be problematic, especially when the color in the target plants or background changes due to differences in species, varieties, or growth stage [29], or when shadows present or background scenes vary [30], making it challenging to adopt these techniques for high-throughput, automated image segmentation in the HTP facility data pipeline.
In this study, we used correlative measurements to demonstrate that low-throughput traits measured in consecutive time points can be predicted with high accuracyby incorporating high throughput surrogate SPA using this method of segmentation.

Genotypes and Growing Condition
Plant materials were two spring wheat genotypes, Yecora Rojo and Seri-82, developed by the International Maize and Wheat Improvement Center (CIMMYT). Yecora-Rojo (pedigree: Ciano 67//Sonora 64/Klein Rendidor/3/II-8156 = II23584) is a released cultivar (CItr 17414) in California. Seri-82 (PI 591774) is a derivative of Kavkaz with wide geographical adaptation (pedigree: Kavkaz/Buho sib//Kalyansona/Bluebird). The potting medium was a 1:2 volumetric ratio of fine clay (Profile) to organic bagged soil (Sungro). Each pot was filled with approximately 2.2 ± 0.06 kg freshly prepared medium, which filled the pot all the way, leaving 3 cm of the pot empty. Our saturation optimization experiment revealed that the addition of 2.8 kg water will bring each pot to saturation (field capacity). This is a 126% gravimetric water content at saturation. For watering, we irrigated the plants using 500 mL (500/2800 = 18% of water content in saturated state) of 80 ppm nutrient solution of Peters Excel Cal Mg 15-5-15, six days per week. Initially, three seeds were planted, and upon successful emergence of one seedling, the others were extracted from the soil. Plants were grown in a walk-in controlled environment at 60% relative humidity (RH) under 800 µmol of photosynthetically active radiation (PAR), provided by an array of ceramic metal halide lamps (Philips Elite Agro T12), with a 14 h photoperiod, where the day and night temperatures were 25 and 20 • C, respectively. Plants were grown in the controlled environment growth chamber for a total of 53 days.

Experimental and Imaging Procedure
Each genotype was planted in 40 pots, keeping one single plant in each pot. The 40 pots were randomly divided into eight groups of five pots each. In each sampling time, i.e., 21,25,30,35,39,44,49, and 53 days after planting (DAP), five pots were imaged to produce high throughput surrogates and immediately destructively analyzed to produce ground-truth phenotype data. The plants were imaged in a closed imaging booth. At each imaging time points, plants were carried into the imaging booth on an automated conveyer belt and then stood on a rotation table, where we captured side view images at 12 angles, each rotating 30 • . Plants were imaged using a custom-made RGB imager (Aris, Eindhoven, Netherlands), with a standard 5 megapixel CMOS sensor providing a 23 fps frame rate (Basler Ace, Ahrensburg, Germany). A monochrome camera (acA2440-20gm, Basler Ace, Ahrensburg, Germany) with a CMOS sensor of the same resolution and the same frame rate was equipped with a long-pass filter to capture chlorophyll fluorescence images. The imaging box was equipped with a blue light emitting diode (LED) array and a white LED array that provided the strobing light for image acquisition. Every time that a plant came into the imaging position on the rotation table, the blue LED array fired up to excite a chlorophyll fluorescence signal from the plant. The monochrome camera was synchronized with the blue LED light strobe for chlorophyll fluorescence image acquisition. The white LED then fired up, synchronized with the color camera to acquire an RGB image of the plant. The images that were taken by the two cameras were aligned by using a built-in image registration function.

Image Analysis
The image analysis processes started with creating the segmentation mask using the monochrome chlorophyll fluorescence image. Using the collected chlorophyll fluorescence image, Otsu's [31] algorithm was adopted to establish a threshold to separate foreground and background based on the histogram of the gray values in the image. A binary mask was generated using this threshold; any pixels with values smaller than the threshold were considered background and assigned a value of zero (0); otherwise, the pixel belonging to the plant was assigned a value of one (1). A sequence of binary image morphological operations, including opening and closing operations, followed the binarization to minimize undesirable noises before establishing the final mask of the plant. Using the established mask, traits such as the side-projected plant area, the overall plant height, and width of the plant, etc., were calculated. The side projected area (SPA) of the plant was calculated as the sum of the pixels in the mask. The plant height was calculated as the length of the rectangle that enclosed the plant mask. Initially, these measurements were represented in number of pixels. A set of pre-established calibration functions then converted the pixel numbers into physical units (e.g., cm and cm 2 ). All image analysis processes were implemented in MATLAB (Mathworks, Natick, MA, USA, Version 2019a). The average of SPA derived from 12 different angles was used to indicate each plant's SPA. According to this, for each sampling time point and each variety, there were five replicates of SPA data. The five plants of each genotype were then destroyed for low-throughput trait analysis.

Trait Measurement
Destructive sampling began nearly two weeks after plant emergence when plants reached a height of 20 cm. In this experiment, eight sampling time points were completed, including 21,25,30,35,39,44,49, and 53 days after planting (DAP), where in each time point five pots from each genotype were destructively characterized for total leaf area and total phytomass production. Traits measured destructively were leaf area, leaf dry matter, total phytomass (above-ground dry matter), heading date, tiller number, and number of spikes per plants. Leaf area (LA) was measured using an LI-3100 Area Meter (LI-COR, inc. Lincoln, NE, USA) and expressed as cm 2 plant −1 (or cm 2 ). Plant tissues were dried at 65 • C for 72 h. The dry weight of plant leaves, leaf dry weight (LDW), was measured after drying and expressed in mg plant −1 (or mg). Above-ground biomass (BIO), representing all dry matter above-ground, was also measured after drying and expressed in mg plant −1 (or mg). Relative growth rate (RGR) was measured based on biomass produced and defined as total plant dry weight increase in a time interval in relation to the initial weight and expressed as gg −1 day −1 [3]. In addition, we calculated RGR based on the increases in unit time of the side projected area and expressed it as mm 2 mm −2 day −1 .

Data Analysis
For demonstrating growth across the time points and differential growth between the two genotypes, we used two-way analysis of variance (ANOVA) with time point and genotype as main factors, with consideration of the interaction effect of genotype and time point. This analysis was performed on three low throughput traits, namely BIO, LDW, and LA, as well as the high-throughput trait SPA. For demonstrating how well SPA predicts other low throughput traits, we used simple linear regression by taking SPA as an independent variable and BIO, LDW, and LA as dependent variables in each model. For each model, the coefficient of determination (R 2 ) and the Pearson correlation coefficient (r) of SPA with each of the low throughput traits were reported. Two-way ANOVA was performed using "lm", "anova", and "summary" commands in R environment [32] with five replicates. All graphs were made using Microsoft Excel.

Low Throughput Phenotyping (LTP)
The progression of LTP traits are shown in Table 1. In Yecora-Rojo, LA increased from 160 to 810 cm 2 plant −1 . In Seri-82, LA increased from 118 to 1800 cm 2 plant −1 . A similar increasing trend was observed in LDW, with Seri-82 always being greater than Yecora-Rojo beyond 30 DAP. The BIO showed an increasing trend as well, but Yecora-Rojo BIO was significantly greater than Seri-82 BIO. The total biomass of Yecora-Rojo at 21 DAP was 920 mg, rising to 45,390 mg at the eighth time point (53 DAP). In Seri-82, BIO increased from 730 mg to 41,600 mg over the same period.
The other difference among the two genotypes was noted in days of heading and number of tillers. Heading occurred for Yecora-Rojo at 39 DAP, but not until between 49 DAP and 53 DAP time points for Seri-82 ( Figure 1). Yecora-Rojo produced a total of 25 tillers plant −1 , of which 24 were fertile. Seri-82 produced 21 tillers, with only 8 fertile tillers at 53 DAP. The tillering stage of the two genotypes was also different by about 7 d.
Sensors 2020, 20, x FOR PEER REVIEW 6 of 12 The other difference among the two genotypes was noted in days of heading and number of tillers. Heading occurred for Yecora-Rojo at 39 DAP, but not until between 49 DAP and 53 DAP time points for Seri-82 ( Figure 1). Yecora-Rojo produced a total of 25 tillers plant −1 , of which 24 were fertile. Seri-82 produced 21 tillers, with only 8 fertile tillers at 53 DAP. The tillering stage of the two genotypes was also different by about 7 d. The two-way ANOVA indicated that the individual effects of time point and genotype were highly significant on the four traits, indicating that the LA, LDW, BIO, and SPA increased in both genotypes continuously from the first time point to the last time point, and that the two genotypes showed different growth patterns ( Table 2). The interaction effect of time point and genotype was only significant for LDW and LA. Table 2. Significance levels of time points, genotype, and their interaction for the four traits evaluated via two-way analysis of variance.  The two-way ANOVA indicated that the individual effects of time point and genotype were highly significant on the four traits, indicating that the LA, LDW, BIO, and SPA increased in both genotypes continuously from the first time point to the last time point, and that the two genotypes showed different growth patterns ( Table 2). The interaction effect of time point and genotype was only significant for LDW and LA.

The Predictive Performance of High Throughput Phenotyping
The SPA showed the same trend as the LTP traits (Table 1). It increased from 12,542 to 81,479 mm 2 in Yecora-Rojo and from 10,882 to 105,419 mm 2 in Seri-82 during the same period. Our objective was to examine how SPA predicts LTP traits. Therefore, we fitted linear regression by taking SPA as predictor variable and the LTP traits, i.e., LA, LDW, and BIO, as dependent variables. The highest prediction accuracy was observed for LDW (R 2 = 0.98 for Yecora-Rojo and R 2 = 0.98 for Seri-82) followed by LA (R 2 = 0.89 for Yecora-Rojo and R 2 = 0.98 for Seri-82). BIO was linearly predicted with the least accuracy (R 2 = 0.86 for Yecora-Rojo and R 2 = 0.91 for Seri-82). While the relationship between SPA and LA (and LDW) was shown to be a linear fit, the relationship between SPA and BIO was not linear and resembled a quadratic or exponential trend. When SPA was used to predict LA, LDW, and BIO, the regression was highly significant. The correlation coefficient between SPA and LA, LDW, and BIO was found to be between 0.92 and 0.99 for Yecora-Rojo and between 0.91 and 0.99 for Seri-82 (Table 3). Table 3. Simple linear regression of leaf area (LA), leaf dry weight (LDW), and biomass (BIO) as predicted by side projected area (SPA). All regression models are significant (p < 0.01). R 2 is the coefficient of determination, and r is the Pearson coefficient of regression.

Relative Growth Rate Measurements on the Basis of Low-and High-Throughput Traits
Relative growth rate (RGR) is the increase of plant biomass relative to size. This parameter can be expressed for other traits such as LA, or any other quantitative traits that are subject to change over time. In our experiment, we measured relative growth rate for LDW, BIO, and SPA. The RGR values were relatively higher during 21-30 DAP, dropping to lower values over time (Table 4) Interestingly and with high accuracy, the RGR that was calculated based on SPA (RGR SPA ) values also followed the same trend for both genotypes. It decreased from the first interval through to the last interval (Table 4). One objective of our study was to determine whether RGR measured on the basis of HTP can replace the time and large sample consuming LTP based measurements. The results showed that the regression relationship was significant, and that RGR SPA could explain up to 93% of Sensors 2020, 20, 6501 8 of 12 the variations observed in RGR BIO (Figure 2). The coefficient of correlation of RGR BIO and RGR SPA was 0.97 for Yecora-Rojo and 0.91 for Seri-82, respectively.

Relative Growth Rate Measurements on the Basis of Low-and High-Throughput Traits
Relative growth rate (RGR) is the increase of plant biomass relative to size. This parameter can be expressed for other traits such as LA, or any other quantitative traits that are subject to change over time. In our experiment, we measured relative growth rate for LDW, BIO, and SPA. The RGR values were relatively higher during 21-30 DAP, dropping to lower values over time (Table 4) Interestingly and with high accuracy, the RGR that was calculated based on SPA (RGRSPA) values also followed the same trend for both genotypes. It decreased from the first interval through to the last interval (Table 4). One objective of our study was to determine whether RGR measured on the basis of HTP can replace the time and large sample consuming LTP based measurements. The results showed that the regression relationship was significant, and that RGRSPA could explain up to 93% of the variations observed in RGRBIO ( Figure 2). The coefficient of correlation of RGRBIO and RGRSPA was 0.97 for Yecora-Rojo and 0.91 for Seri-82, respectively.

Side-Projection Area Efficiency during Growth
Side-projected area is an efficient predictor of biomass up until heading (Figure 3, top panel). However, with the transitioning from vegetative growth to reproductive and emergence of heads, a linear function of SPA is a poor predictor of biomass, as indicated by a drop in the coefficient of determination (Figure 3, bottom panel), and underestimates the biomass. An exponentially fitted function provides a better prediction of biomass, as indicated by higher coefficients of determination. panel). However, with the transitioning from vegetative growth to reproductive and emergence of heads, a linear function of SPA is a poor predictor of biomass, as indicated by a drop in the coefficient of determination (Figure 3, bottom panel), and underestimates the biomass. An exponentially fitted function provides a better prediction of biomass, as indicated by higher coefficients of determination.

Discussion
The overall goal of this study was to validate the manual and destructive growth analysis of wheat using non-destructive machine vision surrogates. Explaining crop growth by using high throughput phenotyping facilitates crop physiology measurements and accelerates crop genetics

Discussion
The overall goal of this study was to validate the manual and destructive growth analysis of wheat using non-destructive machine vision surrogates. Explaining crop growth by using high throughput phenotyping facilitates crop physiology measurements and accelerates crop genetics research for growth-related traits. In this study, machine vision supported by automated plant RGB imaging was used to produce surrogates for wheat growth with high accuracy. The accuracy of digital growth analysis was demonstrated by correlation between HTP and LTP ground-truth measurements. The correlation results reported in this study were also reported by previous researchers. Rajendran et al. [33] reported a positive relationship between image-based surrogates and shoot biomass (R 2 = 0.94) and leaf area width (R 2 = 0.92) in Triticum monococcum species, and Ahmad et al. [34] reported high correlations between manual and digital methods on leaf area index estimation for wheat (R 2 = 0.53-0.93), barley (R 2 = 0.58-0.96), and oat (R 2 = 0.72-0.97). Similar results were obtained from other studies for wheat and barley [35]; soybean and corn [36]; Arabidopsis rosettes [22]; and forest tree seedlings [37].
Additionally, our results indicated that HTP can reduce the number of replicates and the magnitude of experimental error generated by plant-to-plant variation. For example, destructive sampling for eight time points and five replicates per time point requires 40 plants and yet may not be able to depict the changes in the rate of growth because of low resolution. In contrast, non-destructive image-based growth analysis can rely on fewer pots (e.g., five) and can provide a daily analysis of growth. Therefore, the overall benefits will be the time and cost that are needed for completing the growth rate measurements. This is significant for enabling assessment of large germplasm for crop improvement [38,39].
High throughput phenotyping can offer non-destructive methods that enable measurements on the same plant over time without the need to measure the plant destructively and thereby making it possible to track the growth and development of the plant over time or effects on growth responses in different environmental conditions. By establishing a high-dimensional data matrix that characterizes the plant's growth and development, HTP enables in-depth analysis of the interactions among distinct traits, which provides opportunities to identify key parameters that correlate with desirable traits [40], such as drought tolerance. More importantly, the machine vision can eliminate confounding effects of plant-to-plant variation in repeated measurements of growth analysis. The ability to control the measurements of individual plants with higher accuracy measurements can open the opportunity for measurements of short lived mortal genetic mapping population (such as F2 or backcross mapping populations), where each individual of the population has a unique genetic identity. Another benefit is that HTP enables us to fill the existing gap in plant physiology by decreasing the interval between two consecutive measurements, where a trait measurement from a given plant can be repeated over a very short time, mimicking continuous measurement for continuous growth. This will enable the capture of small changes and their timing in plant growth.
Despite these advantages, we identified limitations in using HTP for digital growth analysis. One area of improvement for surrogate-based modeling of biomass is the time of transitioning from the vegetative to reproductive stage. After heading, SPA was not predictive of organ biomass and underestimated the biomass if modelled linearly, because spikes have higher specific weight than leaves and stems as they progress through the grain-fill process. Another area that needs to be approached cautiously is interpretation of results obtained from pot and single-plant measurements to field-scale analysis, because the true grain yield expression is always measured in the plant community in field and in direct interaction with the environment. Using single plant analysis clearly misses the effect of plant-plant interactions and competition.

Conclusions
In this study, we presented a method for accurate estimation of plant LA, LDW, and BIO using an HTP surrogate. We converted images taken from plants into numeric values that represented growth in each time point. This approach provided an accurate and practical model for the estimation of wheat leaf traits and shoot biomass as a substitute for the conventional destructive method with high accuracy. The results indicated that the method can reduce the number of replicates and the error generated by plant-to-plant variation with very short interval times compared to the destructive method. In this study, we calculated RGR using a low throughput method based on LA (RGR LA ), LDW (RGR LDW ), and BIO (RGR BIO ) and showed that RGR based on the HTP surrogate (RGR SPA ) is well correlated with those measured by LTP. HTP has shown to be promising for monitoring the growth and development of the plant over time and for studying plant responses to environmental stresses. However, methodological improvements are needed to better model the plant growth parameter to capture biomass changes after transitioning from the vegetative to reproductive stage.