Nitrogen diagnosis based on dynamic characteristics of rice leaf image

Digital image processing is widely used in the non-destructive diagnosis of plant nutrition. Previous plant nitrogen diagnostic studies have mostly focused on characteristics of the rice canopy or leaves at some specific points in time, with the long sampling intervals unable to provide detailed and specific “dynamic features.” According to plant growth mechanisms, the dynamic changing rate in leaf shape and color differ between different nitrogen supplements. Therefore, the objective of this study was to diagnose nitrogen stress levels by analyzing the dynamic characteristics of rice leaves. Scanning technology was implemented to collect rice leaf images every 3 days, with the characteristics of the leaves from different leaf positions extracted utilizing MATLAB. Newly developed shape characteristics such as etiolation area (EA) and etiolation degree (ED), in addition to shape (area, perimeter) and color characteristics (green, normalized red index, etc.), were used to quantify the process of leaf change. These characteristics allowed sensitive indices to be established for further model validation. Our results indicate that the changing rates in dynamic characteristics, in particular the shape characteristics of the first incomplete leaf (FIL) and the characteristics of the 3rd leaf (leaf color and etiolation indices), expressed obvious distinctions among different nitrogen treatments. Consequently, we achieved acceptable diagnostic accuracy (training accuracy 77.3%, validation accuracy 64.4%) by using the FIL at six days after leaf emergence, and the new shape characteristics developed in this article (ED and EA) also showed good performance in nitrogen diagnosis. Based on the aforementioned results, dynamic analysis is valuable not only in further studies but also in practice.


Introduction
During recent years, spectral analysis and digital image processing technology have been widely used in plant nitrogen (N) diagnosis. This technology has greatly increased the PLOS  efficiency of monitoring the growth station of field crops in a timely and non-destructive manner, thus providing strong technical support for producing high yields in a favorable environment [1,2]. Hyperspectral imaging technology is a newly developed technique that combines spectral and spatial information simultaneously. Recent studies have demonstrated the potential of hyperspectral imaging technology for plant nutrition research [3,4]. However, the high cost of images and the strict requirements of the operating environment have limited its application in practice [5]. Alternatively, digital imagery is commonly used in N diagnosis, which not only records spectral and morphological information but also provides a low-cost alternative with suitable image resolution. Current studies have demonstrated the value of digital images in N diagnosis, in particular for characteristics extracted from digital images providing effective parameters in N diagnosis [6,7]. Moreover, digital imaging devices are portable and easy to operate for data acquisition, and subsequently allow digital image processing to provide a more practical context in terms of future development [8][9][10].
Rice plants with N deficiency have some specific symptoms such as small leaves, leaf etiolation from the tip [11][12][13]. These symptoms are closely related to N content as have been widely shown in previous studies [6,14]. Statistics show that the sampling interval is commonly around 10 days [7,10,15]; however, further studies have shown that some detailed changing processes have been missed within this time period [16]. This suggests that some specific and effective characteristics which appear during a short time period could be missed due to the long sampling interval, and for this reason a short sampling interval is considered important in identifying distinctive "dynamic features" in N diagnosis. Supported by previous studies, rice leaf images contain more specific information than that of canopy images, and they are deemed more suitable for N diagnosis compared to those of the canopy [15,17,18]. It is also important to note that N content strongly influences the changing rates in rice leaf characteristics, and it is responsible for their final condition. This indicates the possibility of "dynamic features" having a close relationship with leaf N content, and more importantly identifying a "hidden trait" within rice leaves allowing for a deeper understanding of growth and development [16]. For example, leaf etiolation is a specific symptom of N stress and the etiolation rate can effectively show leaf N stress levels. Therefore, analyzing the temporal and spatial changes in a rice leaf by monitoring leaf growth continuously is meaningful in research and practice.
Dynamic analysis is widely used in plant growth modeling and plant genotypic studies [16,19,20]. Duan [21] quantified canopy structure to characterize early plant vigor in wheat genotypes by analyzing a time series of images. Neilson [16] described the growth models of two sorghum hybrids under N and water deficiency by continuous monitoring. Similarly, Poiré [19] combined time series images with destructive harvest of shoots and roots to evaluate the N and phosphorus use efficiency of different genotypic plants. These studies demonstrate a plant's response to nutrition or water stress from the perspective of plant growth modeling or plant genotypic studies. Nevertheless, the perspective of plant N diagnosis has not yet been considered, and thus further research identifying specific characteristics and building diagnostic models is required for understanding the important dynamic characteristics of N diagnosis.
Based on the aforementioned, it is valuable and feasible to apply dynamic characteristics in N diagnosis. In this study, we focused on analyzing the dynamic changes in rice leaves for N diagnosis. To the best of our knowledge, this is the first attempt to use dynamic changing rates in rice leaf characteristics in N diagnosis. The objectives of this study were to (i) analyze the relationship between N supplement and the dynamic characteristics of rice leaves, (ii) explore an approach to quantify dynamic characteristics, and (iii) identify effective dynamic characteristics and establish N diagnostic model.

Experimental design
Experiments were implemented during both 2014 and 2015, with the application of hydroponics to cultivate rice plants under different levels of N stress. Rice seeds (cultivar ZheYou-NO.1) were pre-germinated for 3 days, in moist sand at 26˚C. After 15 days' cultivation, seedlings were individually transplanted in 5 L polyvinyl chloride pots that contained precisely controlled nutrient solution. Experiments were conducted in a glasshouse at Zijingang campus, Zhejiang University (30˚17' N, 120˚05' E, Hangzhou, China). Plants were grown under natural light conditions. The temperature of the glasshouse was programmed to be 30˚C/25˚C (day/ night) and the relative humidity 50%. Nutrient solution formula from the International Rice Research Institute (IRRI) was used to cultivate the rice plants and replaced every 15 days. The experiment applied four different N level treatments, with five replications for each level. Four N levels (ammonium nitrate N1: 0 mg/L, N2: 57.20 mg/L, N3: 85.70 mg/L and N4: 114.30 mg/ L, respectively) via hydroponic solutions were applied to different pots. N1 represented extreme deficiency, N2 medium deficiency, N3 slight deficiency and N4 normal supplement. Every 5 days the pH of the nutrient solution in each pot was measured and adjusted to 5.5-6.5 using 1 mol/L NaOH.

Image acquisition and processing
For image acquisition, scanning was applied due to its advantages in data acquisition and processing, especially in closed environments that can effectively eliminate the influence of the environment and operator [18,22]. Rice leaves were scanned (EPSON GT20000, Epson Inc., Suwa, Nagano prefecture, Japan) every 3 days from 20 days after transplanting (DAT 20) to 44 days after transplanting (DAT 44) during both 2014 and 2015. The top four leaves (including the first incomplete leaf) of each rice plant were scanned at the sampling time and experiments were conducted nine times. Then, a total of 1920 rice leaf images were processed in MATLAB 2013b (MathWorks Inc., Natick, Massachusetts, USA). MATLAB provides programs that have been designed to extract leaf color and shape characteristics to describe the dynamic change of leaf extension and etiolation.

Characteristics extraction
2.3.1 Characteristics of leaf shape. Leaf area (LA) and leaf perimeter (LP) were chosen to describe the process of leaf shape change. These parameters were chosen because N deficiency results in stunted growth and therefore the process of leaf extension and etiolation would be distinct among N supplements. Hence, the changing rate in leaf area and perimeter are important characteristics for evaluation of leaf growth status.
Under N deficiency, the N element would be transferred from an old leaf to a young leaf, which would accelerate the process of leaf etiolation in the old leaf. Therefore, the etiolation area (EA) and the degree of leaf etiolation (ED) can be good indicators of N stress levels. Image segmentation was used to extract the etiolated part (Fig 1) and etiolation area was calculated. Based on the aforementioned, EA and ED were calculated to describe leaf etiolation status as follows (Eq (1)): 2.3.2 Characteristics of leaf color. Color characteristics were extracted from different parts of the leaves according to different leaf positions. Color characteristics of the first incomplete leaf (FIL) were extracted from the entire leaf, because the color of the FIL was generally uniform. However, for a fully expanded leaf, color characteristics were extracted from the leaf tip because this part was more sensitive to N stress than the entire leaf. According to previous research, the leaf tip characteristics performed well in N stress diagnosis [23]. Therefore, rice leaves were trisected along their long axes, and the color characteristics of the first one-third were chosen for further analysis. Color indices (Table 1) examined in previous studies that have good correlations with N content were applied to describe leaf color variations [6,15,24].

Dynamic characteristics analysis 2.4.1 Data processing of dynamic characteristics. 1) Growth pattern fitness
To identify the rice plant growth pattern response to N treatments, different types of mathematical models were applied to the leaf shape data such as the power, exponential and logistic functions. The suitability of growth models was evaluated using two statistical methods: coefficient of determination (R 2 ) and Akaike information criterion (AIC). R 2 was used to evaluate the goodness of model fitting and AIC was used to estimate the complexity of model. This process was carried out in IBM SPSS 22.0 (IBM, Armonk, New York, USA).
2) Quantification of dynamic characteristics Relative growth rate (RGR) is an essential parameter that is widely used in botanical science to evaluate plant biomass change [16,19]. It was chosen to describe the dynamic changing rate of rice leaf shape under different N stress. RGR (LA, LP, EA and ED) was calculated according to Hunt [25] as follows: where W1 and W2 are the initial and final characters (LA, LP, EA and ED) at the beginning (t1) and end (t2) of the measurement period, respectively (Eq (2)).  On the other hand, the average changing rate (ACR) derived from the absolute growth rate (AGR), another important parameter in botanical science, integrates plant color indices with time. It was used to describe the dynamic changing of leaf color. ACR was calculated as follows: where X1 and X2 are the initial and final characters (G, NRI, DGCI, ExG and ExR) at the beginning (t1) and end (t2) of the measurement period, respectively (Eq (3)) RGR (leaf shape) and ACR (leaf color) were calculated with the time interval between t1 and t2 being set at 3 days or 6 days. However, because of the senescence of old leaves, some data were missing during the late stage. Consequently, in this study, there were 7 data sets with a 3-day interval and 3 data sets with a 6-day interval used for further analysis (Fig 2).

Statistical analysis method.
The quantified dynamic characteristics of leaf shape and color were applied in parameter evaluation and model establishment.
(1) Selection of optimal characteristic parameters One-way analysis of variance (ANOVA) was conducted to estimate parameters and screen optimal parameters for further analysis. F-values and p-values were used to evaluate the intergroup and intra-group differences. A p-value less than 0.05 means the difference is significant and the indices were effective in classification. The optimal parameters to be used in model establishment were selected according to p-value and F-value.
(2) Model establishment and validation To calculate the optimal diagnostic time and leaf position, a total of 10 sets of data sets (7 using a 3-day interval and 3 using a 6-day interval) from different growth stages were used in model establishment and validation. However, to make a further exploration of dynamic characteristics, we combined the single data sets in time order to improve the diagnostic effect ( Table 2). LibSVM was adopted to build the N diagnostic model, and leave one out cross validation (LOOCV) was used to validate the model in MATLAB. Establishment of data sets. P1 to P8 represent the data set calculated using a 3-day interval. P1' to P4' represent a data set calculated using a 6-day interval.

Dynamic change of leaf expansion.
Model fitting results of the power, exponential and sigmoidal logistic functions are shown in Table 3. All three models produced a curve with an R 2 greater than 0.80. In particular, the sigmoidal logistic function produced a curve an R 2 near 0.99 which was higher than that of the others. In contrast to the R 2 , a lower AIC value represents a higher quality model. The AIC of the sigmoidal logistic function ranged from -17.3333 to 1.2889 in the leaf area growth model, and from 4.3837 to 14.5591 in the leaf perimeter model, both of which were obviously less than those of the other models. Utilizing the aforementioned information, the sigmoidal logistic function was chosen for subsequent analysis.
The extension process of leaves from emergence to full expansion containing different N supplements are shown in Fig 3. Our results indicate that the N supplement influences the color and size of rice leaves, with rice leaves becoming greener and larger with a higher N supplement (Fig 3A), and the growth rate of leaf area and leaf perimeter increased with the increasing N supply (Fig 3B). Overall, it took approximately 10days to reach the peak of leaf area (perimeter). Furthermore, different variation among the four treatments could also be observed in the RGR changes. As shown in Fig 3C, RGR decreased with leaf extension and the leaves with higher N supply initially showed greater RGR (leaf area and perimeter) values. Consequently, bigger decline in RGR value could be observed in higher N treatment.

Dynamic change in leaf color during etiolation.
The etiolation process of the leaf tip showed that the leaves became "yellower" and "narrower" with decreasing N supply, and the etiolation area increased more rapidly under low N treatments compared to that of high N treatments (Fig 4). The dynamic changing of leaf color was observed utilizing the "slice function" in MATLAB which plotted the relationship between RGB value and time (DAT). This showed that leaves with higher N content produced a lower RGB value. During leaf etiolation, the RGB value increased at different rates and the increasing RGB value meant leaves were losing their "green" and becoming "yellower"

Optimal characteristic parameters in discrimination
Leaf samples were divided into four groups according to leaf position: the first incomplete leaf (FIL), and the 1 st , 2 nd , and 3 rd fully expanded leaves from the top of the rice plant. Table 4  shows the F value ranged from 0.04 to 44.67. Characteristic parameters whose F values were greater than 2.92 (the p values were less than 0.05) had an obvious difference among the different N treatments. From the distribution of parameters whose p-values were less than 0.05, most were extracted from the 3 rd leaf, followed by the FIL, 2 nd leaf and 1 st leaf. It was observed that leaves from different leaf positions had different growth statuses; in fully expanded leaves (which includes the 1 st , 2 nd and 3 rd leaves), dynamic changes in color characteristics were more sensitive to N stress than shape characteristics during the leaf etiolation stage. For the FIL, it was the shape characteristics that showed more distinction among the N treatments during the leaf extension stage.
The P1, P2, P3 and P1' data sets for the FIL were established during the leaf extension stage, and their shape indices (LA and LP) had higher F values than those of the color indices. Meanwhile, the P3P4P5 and P2' data sets of the 3 rd leaf were built during the leaf etiolation stage, and the color indices (DGCI, G, ExG, ExR, and NRI) and etiolated indices (ED, EA) performed better than LA and LP. Similarly, the color characteristics of the 1 st leaf and 2 nd leaf in the P6, P7 and P3' data sets had a higher F value. Table 2, 18 data sets were established in every leaf position. Every data set was used to establish a diagnostic model. Table 5 showed the leaf position that achieved the best validation accuracy for each data set. Overall, the FIL and the 3 rd leaf performed better than the other leaves, and the combined data sets achieved higher accuracy than single data sets.

Selection of optimal leaf position. As indicated in
In single data sets, "dynamic features" of P1, P2 and P1' were more distinctive in FIL than in other leaves as the FIL was captured during the extension stage, P1' from the FIL got the best accuracy among them (training accuracy 77.3%, validation accuracy 64.4%). For the P3, P5, P2' and P3' data sets, the 3 rd leaf showed obvious etiolation, dynamic characteristics of the 3 rd leaf were more sensitive than those of the others. During the late stage (P6 and P7), the upper leaves (FIL and 1 st leaf) performed better than the lower leaves (the 2 nd leaf and 3 rd leaf). In combined data sets, the FIL did best in the combined P1P2 and P1-P7 data sets, and P1-P7 achieved the highest accuracy (training accuracy 80.5%, validation accuracy 72.7%). For the combinations that contained P3-P6 (or P2'-P3') data sets, the older leaves performed better than the others.

Selection of best diagnostic time.
Since the FIL and the 3 rd leaf performed better than the others in N diagnosis, we further combined these two leaves in the model establishment to improve diagnostic accuracy (Table 6).
In single data sets, the diagnostic accuracy (both 3-days interval and 6-days interval) decreased with time and performed poorly. The best accuracy appeared in the P2' data set, because the optimal characteristics of the FIL and 3 rd leaf were in different data sets, the "a" represents p-value <0.05 "b" represents p-value < 0.01 "c" represents p-value < 0.001. accuracy did not improve significantly. Subsequently, we combined optimal data sets for further analysis. Data sets P1P2P3 achieved the optimal accuracy (training accuracy 84.1%, validation accuracy 72.7%), but the performance of the other combinations was poor.

Discussion
Our results show that although data acquisition in dynamic analysis is more time-consuming compared to that of other methods, it still has several distinct advantages that allow effective characteristics to be identified in nitrogen diagnosis.

Importance of "dynamic features" in N diagnosis
(1). In plant N deficiency studies, the sampling interval was too long to detect detailed "dynamic features", and some specific characteristics were missed [26,27]. According to (2). Time series analysis of plant images using non-destructive technology provides both quantitative and qualitative information to understand physiological phenomena such as plant growth pattern responses to environmental change [28,29]. The results here are in agreement with previous studies and this consequently could be applied in plant growth studies [16].
(3). Dynamic analysis is important for plant studies to identify "hidden traits," particularly in terms of nutrition diagnosis [16]. In this study the dynamic change in leaf etiolation was quantified using EA and ED, with the results (Table 4) indicating that ED and EA achieved good performance in dynamic analysis.
(4). As is known, "dynamic features" exist over the entire life of a rice leaf; this indicates that N diagnosis can be implemented earlier. According to Tables 5 and 6, using the FIL can achieve an acceptable validation accuracy of 64.4% at DAT 26, and another acceptable validation accuracy of 72.7% calculated at DAT29 was achieved integrating FIL and the 3 rd leaf. This is earlier when compared to commonly used methods, implying that diagnosis could be implemented at the emergence of other leaves which would be more meaningful in practices. (5). To some extent, dynamic analysis of specific characteristics is simple and feasible. Our results show that the dynamic characteristics of the 3 rd leaf tip performed well in both ANOVA and modeling. Indicating that using leaf tip information from an old leaf can significantly improve diagnostic operability and efficiency. (6). Aside from the utilization of N stress diagnosis, dynamic analysis also provides a novel method to identify different types of nutrition deficiency, as indicated in Figs 3 and 4 where leaf extension rate and etiolation rate changed with N supplement. According to plant growth mechanism, nitrogen, phosphorus and potassium deficiency always result in stunted growth, but it differs in terms of new leaf growth rate and old leaf senescence rate. Based on the aforementioned, dynamic characteristics would be distinctive under nitrogen, phosphorus and potassium stress and could be ideal indices for diagnosis.

Application of "dynamic features" in N diagnosis
Leaf shape and spectral information acquired at some particular points in time have been widely used in previous studies, where results have shown that the 3 rd leaf or 4 th fully expanded leaf could be an ideal indicator in N diagnosis [12,30,31]. Our findings (Figs 3 and 4) showed the FIL and 3 rd leaf were ideal indicators in dynamic analysis as dynamic characteristics were mainly reflected during the leaf extension and etiolation stages. Dynamic characteristics of the FIL (area, perimeter) and the 3 rd leaf (ED, G, ExR, and ExG) showed obvious differences among the N treatments during different growth stages (Table 4). Employing RGR and ACR to quantify the dynamic characteristics, provided an in-depth understanding of the dynamic process of leaf extension and etiolation in regard to N supplement.
It is known to all that knowing plant nutrition status timely and precisely is beneficial to high crop yield and quality. In this study, the N status has been effectively identified by dynamic leaf characteristics which provides valuable information for field management. With the popularization of digital camera, portable scanner and smartphone, farmers can obtain plant dynamic information easily by using these devices. Based on this, the plant nutrition status could be identified by combining with the diagnostic model. Moreover, to determine the optimal fertilization, soil fertility is also needed in quantification of N supplement. By considering the plant nutrition status and soil fertility status, the farmer could adjust the fertilization properly and timely which would be a further application of leaf dynamics.

Effectiveness of N diagnosis using "dynamic features"
In this article, dynamic nature of rice leaf is quantified by different time intervals. As we can see, dynamic characteristics quantified using a 3-days interval and a 6-days interval have their respective advantages in diagnosis. The shorter time interval contributes to early diagnosis and the identification of "hidden traits" which exist in short time, while the longer time interval could help us to further improve the diagnosis accuracy to some extent.
Estimates showed that previous studies have typically collected data near DAT40 when N stress symptoms became notable [15,32,33]. Therefore, utilizing "dynamic features" to discriminate N stress levels, our methodology greatly increased the possibility of identifying N stress levels during an earlier stage when stress symptoms were not significant. Exploiting "dynamic features" allowed achievement of a validation accuracy of 64.4% (training accuracy 77.3%) in the FIL at DAT 26, while the 3 rd leaf produced higher accuracy at later stage. Because the FIL and the 3 rd leaf were the optimal leaves, integrating these two leaves effectively improved the diagnostic accuracy, especially during the early stage data sets P1P2 and P1P2P3 (Table 6). Although the diagnostic accuracy is not as high as that of previous researches, it shows the feasibility of using dynamic characteristics in early diagnosis. In our opinion, sacrificing a certain extent of accuracy for earlier diagnosis is more valuable.

Future work
According to plant growth mechanisms, plant responses to N deficiency embody not only in leaf characteristics, but also in plant physiology and biochemistry. In this paper, we mainly focused on the responses of leaf characteristics to N deficiency. From the perspective of plant physiology and biochemistry, the adaption to N deficiency in metabolism usually results in etiolation and stunted growth. During this process, the expression of N metabolism related genes in the shoots (nitrate reductase (NR1), nitrate reductase (NR2), glutamine synthetase (GS2), ferredoxin-dependent glutamate synthase (Fd-GOGAT), glutamate dehydrogenase (GDH2), glutamate dehydrogenase (GDH3), etc.) are upregulated under short-term N deficiency. In contrast, their expression would decrease under long-term N deficiency [34][35][36]. In this sense, the responses of the N metabolism related genes to N deficiency differ with N content. Therefore, exploring the relationship between gene expression, leaf characteristics and N content would contribute to the exploration of effective traits for N diagnosis, thus improving the diagnostic effect.

Conclusion
This study provided a detailed and feasible research method of using dynamic leaf characteristics in rice N diagnosis. We demonstrated the effectiveness and feasibility of dynamic analysis, in particular for quantifying specific "dynamic features" and allowing for earlier diagnosis of N stress.
The results showed that dynamic characteristics have distinctive differences among N treatments, particularly the dynamic changing of the FIL (area and perimeter) and the 3 rd leaf (etiolated indices and color indices) which could be ideal indicators in N diagnosis. In particular, leaf etiolation was quantified using two characteristics developed in this study (EA and ED) and achieved good performance in feature selection.
When applying "dynamic features" in modeling, N stress levels can be discriminated earlier than commonly used methods. Utilizing a single leaf, the FIL performed better than the other leaves during the early stage, achieving an acceptable validation accuracy of 64.4% at DAT26. Furthermore, diagnostic accuracy can be obviously improved during the early stage using combined parameter sets of the FIL and the 3 rd leaf.
Finally, this study, from dynamic analysis to diagnostic model establishment, demonstrated the value of "dynamic features" in N diagnosis, thus helping future research both in terms of innovative applications and field management.