The experiment was carried out in a glass greenhouse at the School of Agriculture and Biology, Shanghai Jiao Tong University, China. The pakchoi (cv. Hua wang) were cultivated in a substrate containing available N (332 mg·k⁻¹), phosphorus (1.24 mg·kg⁻¹), potassium (118 mg·kg⁻¹), and organic matter (270.3 g·k⁻¹). The seeds were sowed in a plug tray on 18 November 2017. After a nursing of 18 days, the seedlings were transplanted to pots (volume of 1.47 L) with one plant per pot. The N treatments (CK, T1, T2, and T3 corresponding to N levels of 0, 0.134, 0.163, and 0.191 g N·pot⁻¹, respectively) were applied to each pot after one week of transplanting. The T3 treatment was based on the conventional N application rate. Three blocks were considered with 20 repetitions per block for 240 plants. Urea (Sinopharm Chemical Reagent Co., Ltd, Shanghai, China) was chosen as the N fertilizer. The urea was divided into six parts and applied every seven days. KH₂PO₄ (Sinopharm Chemical Reagent Co., Ltd, Shanghai, China) was selected as the phosphate and potassium fertilizer. The dosages of P and K in each treatment were the same, and the additive amount was 0.4613 g of KH₂PO₄ per pot. The substrate water content was controlled at a maximum water holding capacity of 75 ± 5% using the weighing method. The images and plant tissue were collected at different growth stages: G0, G1, G2, G3, G4, G5, and G6 corresponding to 0, 7, 14, 21, 28, 35, and 42 days after the first N application, respectively.

The greenhouse environmental data were monitored every five minutes using an automatic data logging system (PM-11 Phytomonitor, Bio Instruments S.R.L., Chisinau, Moldova). During the experimental period, the air temperature ranged from −0.70 to 37.95 °C with an average of 13.88 °C. The average relative humidity was 86.4%, reaching a maximum of 99.2% and a minimum of 29.5%. The mean solar radiation in the greenhouse was 1.69 MJ·m⁻²·day⁻¹.

A self-made platform, i.e., a sealed box (80 × 80 × 80 cm³), with a set of camera and supplement light was used for image acquisition (Figure 1). An industrial camera (MV-EM1400C, MicroVision, Shaanxi, China) was mounted at the top center of the box, and a lens with a 12 mm focal length (M1214-MP2, MicroVision, Shaanxi, China) was implemented. Ring-shaped light emitting diodes (LEDs) (illumination: 6341 lux, the correlated color temperature (CCT): 6138K) as the sole light source were mounted on top of the box and centered on the camera. The pakchoi were sent to the platform one by one, and their images were taken from the top view. The aperture and exposure time of the camera were set to F16 and 0.13 s, respectively, throughout the experiment. We collected 29 plant images before N treatment. After N application, the images were acquired from G1 to G6 growth periods with 13–16 replicates for each nitrogen level in a given period. As a result, the total of the plant images collected was 382. The images were stored in an uncompressed PNG format to avoid color artifacts due to compression algorithms.

The images analysis was performed using the OpenCV (3.4.4 version) and Python language (3.6.1). The details are given in the supplementary material section. The image processing included the computation of excess green index (ExG) and grayscale [48,49]. Figure 2 shows the process flow, and the details include: (1) ExG calculation: the original images were transformed by using 2G-R-B, and the two thresholds (where 40 identified the minimum value and 200 the maximum value) were set to identify the green part of the images (including plant and the green area of background). (2) Grayscale processing: the original images were pretreated by graying, and then the threshold value segmentation was used for removing the high brightness area. The grayscale value was set as 240. (3) Intersection calculation: we take the intersection of the mask obtained by the two operations above, and then the noise (small size) was eliminated based on the mathematical morphological open operation. (4) Denoising: the operations of denoising were done by using OpenCV. Removing objects smaller than the one tenth of the mask area was to remove the residual; filling contiguous holes smaller than one-twentieth of the mask area was to restore leaf veins deleted. Hence, the final images were acquired, and they are shown in Figure 3.

An extensive list of phenotypic features, including color, texture, and geometry, were extracted. The extraction of the color features involved the transformation of RGB (Red, Green, Blue), LAB (L is luminosity, A is the range from magenta to green, and B is the range from yellow to blue, and HSV ( Hue, Saturation, Value) color spaces. The texture attributes were extracted based on the grey co-occurrence matrix algorithm. The equations used were obtained from the study by Haralick and Shanmugam [50] and are described in detail in the supplementary material section. The geometrical features, such as the convex hull and boundary box, were extracted based on the plant contour. Table 1 lists the descriptions of the phenotypic features.

The fresh weight and quality, i.e., chlorophyll, flavonoid, nitrate, and soluble protein, were measured at the last harvesting, with data representing the mean of three replicates (three plants per replication) per treatment.

No. 2 leaves counted from the tip of the shoots were chosen to measure the quality traits. Dualex SCIENTIFIC+TM (FORCE-A, Orsay, France) was used for the chlorophyll and epidermal flavonoid measurements. The chlorophyll was used as given by the device because it was calibrated in μg·cm⁻² units. The absorbance in the ultraviolet (UV) region at 375 nm due to the presence of flavonol in the epidermis was used as given by the device without transformation. The soluble protein content was assayed using the Coomassie brilliant blue method, whereas the nitrate content was measured based on the method given by Xiong et al. [51].

The fresh shoot samples were heated at 105 °C for 20 min and dried at 75 °C for 72 h in an air-dry oven (DHG-9030A, Hangzhou, China). The dried samples were ground and sieved (60 mesh), and the nitrogen concentration was measured using Elementar (ELIII, Hanau, Germany).

The data were analyzed to determine Nc based on the method proposed by Justes et al. [52]. First, the analysis of variance (ANOVA, SPSS-16 software package) was used to identify the data points for which N does not limit growth (non-N-limiting) or when it is not in excess (N-limiting) from the experimental data. The N-limiting growth treatment is defined as a treatment in which a supplement of N application leads to a significant increase in the shoot biomass. The non-N-limiting growth treatment is defined as a treatment in which a supplement of N application does not lead to an increase in the shoot biomass. The regression between the shoot biomass (g·plant⁻¹) and the N concentration (mg·g⁻¹ DM) was conducted with the data obtained from the N-limiting growth treatments. The Nc value at the sampling point is calculated using the corresponding shoot biomass obtained from the non-N-limiting growth treatment. Furthermore, the NNI of the pakchoi at each sampling date was determined by dividing the total N concentration of the shoot (Nt) by Nc (NNI = Nt/Nc) based on a previous report on corn [53].

To investigate the robustness and applicability of the three models under different conditions, the test dataset was reorganized into different scenarios, different growth stages (seedling stage and harvest period), and N nutrition status (N excessive or N deficiency). The model performances in these scenarios were evaluated. The simulation values were obtained using the predictive models, and the coincidence degree was analyzed using the linear regression method of SPSS 16.0. The relationship between the measured NNI and the simulated NNI was assessed based on the criteria: (1) adjusted determination coefficient (R²), i.e., the percentage of variance of NNI of the models; (2) relative standard error (RE), i.e., the ratio of the standard error of regression estimate to the mean of the actual value; (3) accuracy, i.e., the linear correlation coefficient between the simulated and actual values [54].

After 42 days of cultivation under different N levels, the aerial parts of the pakchoi exhibited obvious growth differences (Table 2). The fresh weight significantly increased with N application. However, the pakchoi supplied with higher N had lower fresh weight under tested N levels. The phytochemicals in the leaf were also affected significantly (p < 0.05) by N application. N fertilization significantly increased the chlorophyll content and decreased the flavonoid content; however, no differences were found among the tested N treatments. The availability of N had a prominent effect on the nutrition quality. At harvesting, the nitrate content increased with the N application, and the degree of elevation is positively correlated to the nitrogen level. Similarly, the soluble protein content increased with the N application; however, the content after T1 treatment was significantly higher than those after the other two treatments.

The accumulation of biomass and N content in the shoots of the pakchoi was measured. Table 3 lists the results. The biomass increased continuously with the extension of the growth cycle; however, it was regulated by the nitrogen application rate. During the 7–28 days after treatment, the T3-treated pakchoi accumulated the highest biomass; until the 35th day, the T2-treated pakchoi had the highest biomass; on the 42nd day, the T1-treated pakchoi had the highest biomass. Simultaneously, the response of plant N concentration to N application was observed on the 7th day. However, the N concentration did not exhibit any significant difference between the N treatments until the 28th day, and it was always higher than that in the blank control. On the 28th day, the N concentration positively correlated to the nitrogen application rate. On the 35th and 42nd days, the N concentration had no significant difference between T2 and T3; however, it was higher than that in T1.

Under the tested N application rates, the shoot biomass in the T1-treated pakchoi was lower than that in the T3-treated pakchoi during the previous 28 days. Therefore, T1 belonged to the N-limiting growth treatment, and the relationship between plant N concentration and the shoot biomass was established with the data in T1. The Nc value at the sampling point was calculated using the corresponding shoot biomass value of the T1 or T2-treated pakchoi, and the Nc values were 56.08, 63.36, 69.16, 72.14, 76.44, and 79.46 g/kg. Table 4 lists the calculated NNI of the pakchoi at each sampling date.

In order to increase the efficiency of model development, phenotypic features selection is necessary to identify remarkable features which are relevant to plant N concentration. The significance of 65 extracted phenotypic features on classifying N application rates were investigated by using ANOVA (Figure 4). For better illustration, the p value was transformed into −log10 (p). There were 33, 29, 44, 59, 57, and 60 traits with p < 0.05 (i.e., −log10 (p) > 1.301) at corresponding sampling date, respectively. Finally, 23 phenotypic traits with P < 0.05 were selected at the six growth stages (G1–G6) to develop a model. The selected features contained 16 color-related (numbered as 1, 3, 6, 8, 11, 13, 16, 18, 26, 28, 31, 33, 36, 38, 41 and 43 in Table 1 respectively), three texture-related (numbered as 48, 48 and 51 in Table 1 respectively), and four geometry-related features (numbered as 52, 56, 61 and 62 in Table 1 respectively), and they are marked with asterisks in Table 1.

The models were constructed to quantify the ability of image-based features in statistically predicting the NNI. We combined the NNI measurements with the phenotypic features and then divided them into a training dataset and a test dataset. A model was trained on the training dataset and then applied to the test dataset to predict the NNI. The models were developed using three widely used machine learning methods, and the RMSE was used to select the optimal model with the lowest values. The RF model was developed with ntree as 300 and mtry as 2. The SVR model was trained with a kernel parameter cost of 16 and a regularization parameter gamma of 0.03125. In the NN model, the decay and size were tuned to 0.01 and 5, respectively. Cross-validation results of three modeling algorithms applied to the training set, and the values of RMSE ± SD were 0.083 ± 0.008, 0.076 ± 0.004 and 0.085 ± 0.010 for RF, SVR and NN, respectively. The results show that the SVR model exhibits the best predictive ability for NNI.

The relationship between the simulated NNI (predicted using the regression models) and the measured NNI was assessed (Figure 5). The goodness of fit was quantified by R², whereas the predictive ability was evaluated in terms of the RMSE and MAE. The R² values of the three models reached approximately 0.900, showing an excellent goodness of fit. The RMSEs of the RF, SVR, and NN models were 0.082, 0.073, and 0.086, respectively, whereas the MAE values were 0.056, 0.045, and 0.061, respectively, indicating the effectiveness of predicting the plant NNI based on the phenotypic features. Among the three modeling algorithms, the SVR exhibited the lowest RMSE and MAE values, indicating that the SVR model predicted the NNI most accurately.

As mentioned previously, the image-based features can be broadly classified into three categories: color-, texture-, and morphology-related traits. For each type of trait, we constructed a degenerate model using the corresponding traits as the predictors and compared the capability of each trait in predicting the NNI. The developments of three degenerate models were just based on the selected features (marked with asterisks in Table 1). Figure 6 shows the results. The color-related traits exhibited the best predictive performance among the three categories; however, they showed slightly lower performance than the full model in which all the traits were considered. Interestingly, the predictability of the other types of traits (such as the texture- and morphology-related traits) was substantial, indicating that these traits may act as unforeseen factors in NNI prediction.

We investigated the contribution of each feature in predicting the NNI using the full model (Figure 7). In the RF model, the contribution of a feature was determined as an increase in the prediction error (%IncMSE) when phenotypic data for this feature are permuted, thus indicating the contribution of the feature after considering its intercorrelation in the model. We found that the top five most important traits in the RF model in predicting NNI included both color and texture-related features. In the SVR and NN models, the relative importance of each trait was provided in the algorithms themselves. We found that the top five most important traits in predicting NNI only included color-related traits in the SVR model and color- and morphological-related traits in the NN model.

Crop management practices such as the application of N fertilizer and alternative irrigation systems are often practiced to improve N absorption in crops. Management decisions, such as the appropriate selection of fertilization method and nitrogen application rate, are vital for improving the yield and quality in vegetable production systems. The results of this study showed that the yield and soluble protein content increased, whereas the nitrate content decreased under T1 treatment, indicating that the optimal fertilization strategy of low amount and high frequency could help reduce the nitrogen fertilizer input while ensuring the yield and quality of pakchoi. In a conventional fertilization method, a basal fertilizer was employed; however, the topdressing was ignored, leading to a significant amount of residual N in the soil [4]. The conventional strategies for N application increase the input and the loss of N, and they can cause nitrate pollution in groundwater and farmland [2,55,56], which is not favorable to a sustainable agricultural production.

The plant N concentration in field crops has been related to the aerial biomass as follows: Nc=aDM−b, the parameter a represents the N concentration in the shoot biomass (t·ha⁻¹), and the parameter b represents the coefficient of dilution describing the relationship between N concentration and shoot biomass [5,10]. Various studies have indicated that plant N concentration decreases with the growth cycle. However, the plant N concentration in pakchoi remained high and even increased with the accumulation of shoot biomass, probably because of the fertilization strategy of low amount and high frequency in the pot experiment. This fertilization strategy might promote N absorption and utilization, leading to a faster accumulation of N than shoot biomass. Therefore, the relationship between Nc and the aerial biomass is related to the N strategy. Different critical nitrogen concentration curves should be established for different nitrogen fertilization strategies applied to the crop. These results reveal that topdressing could promote the absorption of N fertilizer by plants and contribute to the reduction in N fertilization.

Three types of phenotypic features (color, texture, and geometry-related features) based on visible light imaging were extracted, and the capability of each feature in predicting the NNI was compared. The color-related traits showed the best predictive performance among the three categories for NNI prediction, which might be attributed that the number of color features is 4 times bigger than the number of textures and geometry (16 features vs. 3 features). However, our result was in agreement with Rorie et al. [57] who found a close relationship between a dark green color index and leaf N concentration. The association between the color features and leaf N is mainly due to phytochemicals such as chlorophyll. Thus, leaf N prediction can be performed by non-destructive chlorophyll measurement. However, chlorophyll meters only provide a single point measurement; they cannot work at high resolutions and cannot provide proper information about the spatial structures of the product [48]. In addition, our results showed that the texture-related and morphology-related traits exhibited a good NNI predictability, and a significant difference was observed (P < 0.05) in phenotypic features including color, texture, and morphology. This indicates that there exist limitations in the conventional evaluation of plant nitrogen nutrition status based on only the leaf color features, and that the ability of the texture- and geometry-related traits to predict plant N cannot be ignored.

In this study, the color, texture, and morphological-related features were used as input parameters. RF, SVR, and NN models were developed to evaluate the NNI of pakchoi. The R² values of the three models reached approximately 0.900, and the RMSE value was lower than 0.1. The results revealed that phenotypic imaging combined with a machine learning algorithm can be used to effectively evaluate the nitrogen nutrition status of plants.

Many algorithms exist in the field of machine learning, each having different principles. Thus, the contribution of the same phenotypic feature in the different algorithms can be different. In the top five most contributing features, the RF model included H_median, ASM, H_mean, energy, and b_mean; the SVR model included H_mean, H_median, b_median, b_mean, and G_mean; the NN model included r, equivalent_diameter, H_mean, conture_area, and H_media. Different algorithms have different dependences on the same phenotypic features, resulting in the local adaptability of each algorithm. Therefore, to select an appropriate NNI model under different scenarios, such as growing periods and soil fertility, it is necessary to determine the robustness and applicability of each algorithm.

In our research, the evaluation effects of the three models were compared under different conditions. The results showed that none of the models could effectively predict the nutritional status when the plant N is excessive, indicating that camera vision in the visible range has a poor resolution when the plant N concentration is close to the normal range. Our results are in agreement with those of Tei et al. [58] who concluded that plants with surplus N reflect less visible light than N-deficient plants. Plateau responses occurred at relatively high N contents, suggesting saturation. However, it is unclear under what exact conditions saturation would occur [59,60]; this requires further elucidation. It is necessary to choose other optical sensors or non-destructive testing techniques to evaluate the N surplus status of a plant. For example, Padilla et al. [22] found that chlorophyll readings and canopy reflectance indices could sensitively respond to the absorption of N by muskmelon. Under N-deficit conditions, the RF model performed best with higher R² and lower RE values. We assessed the NNI prediction models throughout the entire growth cycle of pakchoi. The precise prediction of NNI at the seedling stage indicated that the prediction models can be used to guide N topdressing. The precise prediction of NNI at the harvest stage is beneficial for assessing the N uptake and calculating the N consumption by the current crops. The RF model exhibited outstanding robustness and applicability under different scenarios and thus can be used in practice.