Monitoring Plant Status and Fertilization Strategy through Multispectral Images

A crop monitoring system was developed for the supervision of organic fertilization status on tomato plants at early stages. An automatic and nondestructive approach was used to analyze tomato plants with different levels of water-soluble organic fertilizer (3 + 5 NK) and vermicompost. The evaluation system was composed by a multispectral camera with five lenses: green (550 nm), red (660 nm), red edge (735 nm), near infrared (790 nm), RGB, and a computational image processing system. The water-soluble fertilizer was applied weekly in four different treatments: (T0: 0 mL, T1: 6.25 mL, T2: 12.5 mL and T3: 25 mL) and the vermicomposting was added in Weeks 1 and 5. The trial was conducted in a greenhouse and 192 images were taken with each lens. A plant segmentation algorithm was developed and several vegetation indices were calculated. On top of calculating indices, multiple morphological features were obtained through image processing techniques. The morphological features were revealed to be more feasible to distinguish between the control and the organic fertilized plants than the vegetation indices. The system was developed in order to be assembled in a precision organic fertilization robotic platform.


Introduction
Environmental protection allied with health concerns represents an increasingly important trend in the consumer behaviour and has led to the development of green products and organic markets [1]. These markets have experienced exponential growth in recent years, especially in the organic markets that address multiple consumer concerns relating to health, food safety and environmental conservation [2]. Tomatoes are one of the vegetable crops with a higher demand for fertilization, with recommended doses close to 400 kg·ha −1 depending on soil type [3].
Deficient uses of organic fertilization can cause "abiotic diseases" on plants, decrease the suppressive effect of the soil, and the readiness of plants to defend themselves against attacks of plant pathogens [4]. Nutrient deficiency affects the plant growth, with consequences in the production and the profitability of the organic field. On the other hand, the indiscriminate use of Sensors 2020, 20, 435 3 of 21 use a bicameral system, binarization techniques, and morphological algorithms to differentiate between the crop and spontaneous plant species [24][25][26][27].
Both types of real-time precision spraying systems are commercially available or are in process of being commercialized, and have shown a positive reduction in the use of chemical compounds with a reduction of costs and an increase in yield levels (in the case of precision fertilization) compared to traditional sprayers [8].
Image processing techniques with multispectral cameras from visible to near-infrared spectrum are also being used to provide non-destructive plant phenotype image datasets. These approaches have allowed more precise and real-time, high throughput, and high-resolution data for the modelling and prediction of plant growth and morphological development in different conditions, with recent applications in plant health analysis [28][29][30].
More advanced prototypes were made in order to insert these types of sensors in autonomous vehicles. Several terrestrial robotic platforms were developed using multispectral cameras seeking for a more automatic, low-energy cost and accurate analysis of the crop parameters compared to conventional equipped tractors and unnamed aerial vehicles [11,31,32].
Some of these autonomous agricultural devices are capable of identifying weeds and control them using low-dose spraying, mechanical control or thermal control [33]. Other platforms can spray fertilizers in arable crops without the need for heavy tractors, reducing the compaction of the soil and the physical damage to the crops [29]. More recently, modular agricultural robots were developed for mapping different aspects of the cereal plants in parcels of breeding trials [31].
Besides the advances in the robotic agricultural platforms, most of them were developed for conventional orchards and arable crops, presenting a lack of technologies in the context of organic horticulture. For this reason, the SUREVEG team has developed the robotic prototype shown in Figure 1. Both types of real-time precision spraying systems are commercially available or are in process of being commercialized, and have shown a positive reduction in the use of chemical compounds with a reduction of costs and an increase in yield levels (in the case of precision fertilization) compared to traditional sprayers [8].
Image processing techniques with multispectral cameras from visible to near-infrared spectrum are also being used to provide non-destructive plant phenotype image datasets. These approaches have allowed more precise and real-time, high throughput, and high-resolution data for the modelling and prediction of plant growth and morphological development in different conditions, with recent applications in plant health analysis [28][29][30].
More advanced prototypes were made in order to insert these types of sensors in autonomous vehicles. Several terrestrial robotic platforms were developed using multispectral cameras seeking for a more automatic, low-energy cost and accurate analysis of the crop parameters compared to conventional equipped tractors and unnamed aerial vehicles [11,31,32].
Some of these autonomous agricultural devices are capable of identifying weeds and control them using low-dose spraying, mechanical control or thermal control [33]. Other platforms can spray fertilizers in arable crops without the need for heavy tractors, reducing the compaction of the soil and the physical damage to the crops [29]. More recently, modular agricultural robots were developed for mapping different aspects of the cereal plants in parcels of breeding trials [31].
Besides the advances in the robotic agricultural platforms, most of them were developed for conventional orchards and arable crops, presenting a lack of technologies in the context of organic horticulture. For this reason, the SUREVEG team has developed the robotic prototype shown in Figure 1. This prototype is a cart with a manipulator robot and multiple sensors and actuators. The proposed sensors are laser scanners to build three dimensional models of the plants, which will be used to study their growth, and a multispectral camera, which will be used to study the state and health of the crops. Moreover, a system to apply treatments to the soil will be installed in the prototype, consisting of a tank to store the liquid treatment and a sprayer to apply it. Finally, the manipulator robot will be available to accurately place sensors and actuators at the target points.
Having in mind this long-term goal of automatically monitoring crops and applying treatments, this work is focused on the collection of information about the plants by using multispectral imagery. This prototype is a cart with a manipulator robot and multiple sensors and actuators. The proposed sensors are laser scanners to build three dimensional models of the plants, which will be used to study their growth, and a multispectral camera, which will be used to study the state and health of the crops. Moreover, a system to apply treatments to the soil will be installed in the prototype, consisting of a Sensors 2020, 20, 435 4 of 21 tank to store the liquid treatment and a sprayer to apply it. Finally, the manipulator robot will be available to accurately place sensors and actuators at the target points.
Having in mind this long-term goal of automatically monitoring crops and applying treatments, this work is focused on the collection of information about the plants by using multispectral imagery. Specifically, the paper describes the acquisition of multispectral images of tomato plants with different levels of organic fertilization, as well as the estimation of their nutritional states in the early stages through multiple vegetation indices and morphological features by using computer vision techniques.

Location and Growing Conditions
The experiment was carried out on a greenhouse in order to obtain the different spectral responses of plants. This greenhouse is in the facilities of the School of Agronomic, Food and Biosystems Engineering of the Technical University of Madrid, which are located in Madrid, Spain (40 • 26 19.9" N 3 • 44 15.7" W). The tomato (Mina F1 cv.) seedlings were initially cultivated in 3-l pots with height and diameter of 15 cm, filled with a mix of peat substrate (50%) and coconut fiber (50%). After 5 weeks of experiment, the plants were transplanted to 8-l pots with the same mixture.

Trial Design and Fertilization
The experiment was designed as randomized blocks, with four treatments to be applied to four groups of plants (T0, T1, T2 and T3). The trial design is described by Figure 2, which shows the specific moments when the treatments were applied, and the multispectral images were acquired. The fertilizations were carried out weekly using water-soluble organic fertilizer (3% of Nitrogen and 5% of potassium) obtained from beet vinasse and phosphorite, and registered for use in organic production (COMPO-Fertilizante Huerto y Frutales). The initial volume of irrigation and fertilization was 300 mL of solution. One week after transplant, plants were assigned to the different treatments according to the fertilizer label (T0: 0 mL, T1: 6.25 mL, T2: 12.5 mL and T3: 25 mL, which corresponds to T0: 0 g of N and 0 g of K, T1: 0.15 g of N and 0.30 g of K, T2: 0.30 g of N and 0.60 g of K, T3: 0.60 g of N, 1.20 g of K). The plants also received two supplements of vermicomposting: one in the first transplant (T0: 0 mL; T1: 75 mL; T2: 150 mL; T3: 300 mL) and another after the second transplant (T0: 0 mL; T1: 237.5 mL; T2: 475 mL; T3: 950 mL). Images of the plants were acquired with 7-days intervals and processed to compute vegetation indices and analyze them. Specifically, the paper describes the acquisition of multispectral images of tomato plants with different levels of organic fertilization, as well as the estimation of their nutritional states in the early stages through multiple vegetation indices and morphological features by using computer vision techniques.

Location and Growing Conditions
The experiment was carried out on a greenhouse in order to obtain the different spectral responses of plants. This greenhouse is in the facilities of the School of Agronomic, Food and Biosystems Engineering of the Technical University of Madrid, which are located in Madrid, Spain (40°26′19.9″ N 3°44′15.7″ W). The tomato (Mina F1 cv.) seedlings were initially cultivated in 3-l pots with height and diameter of 15 cm, filled with a mix of peat substrate (50%) and coconut fiber (50%). After 5 weeks of experiment, the plants were transplanted to 8-l pots with the same mixture.

Trial Design and Fertilization
The experiment was designed as randomized blocks, with four treatments to be applied to four groups of plants (T0, T1, T2 and T3). The trial design is described by Figure 2, which shows the specific moments when the treatments were applied, and the multispectral images were acquired. The fertilizations were carried out weekly using water-soluble organic fertilizer (3% of Nitrogen and 5% of potassium) obtained from beet vinasse and phosphorite, and registered for use in organic production (COMPO-Fertilizante Huerto y Frutales). The initial volume of irrigation and fertilization was 300 mL of solution. One week after transplant, plants were assigned to the different treatments according to the fertilizer label (T0: 0 mL, T1: 6.25 mL, T2: 12.5 mL and T3: 25 mL, which corresponds to T0: 0 g of N and 0 g of K, T1: 0.15 g of N and 0.30 g of K, T2: 0.30 g of N and 0.60 g of K, T3: 0.60 g of N, 1.20 g of K). The plants also received two supplements of vermicomposting: one in the first transplant (T0: 0 mL; T1: 75 mL; T2: 150 mL; T3: 300 mL) and another after the second transplant (T0: 0 mL; T1: 237.5 mL; T2: 475 mL; T3: 950 mL). Images of the plants were acquired with 7-days intervals and processed to compute vegetation indices and analyze them.

System Overview
The general system was composed of images acquired from the greenhouse experiment with different levels of fertilization. The evaluation system was composed of a bracket support stands clamp supporting the sensor and the power supply (battery). The platform has the possibility of height adjustment for the sensor. The distance between the sensor and the bottom of the pot was set at 0.7 m during the first 5 weeks and moved to 1.4 m after Week 5 (Figure 3), since the plants were too high to be observed from the initial distance. The higher distance between the camera and the plant makes the plant image smaller, so the leaf area represented by the pixel area became smaller after that modification. In order to model the data, regression analysis using the mean values was performed, so the results could be extrapolated for the next weeks.

System Overview
The general system was composed of images acquired from the greenhouse experiment with different levels of fertilization. The evaluation system was composed of a bracket support stands clamp supporting the sensor and the power supply (battery). The platform has the possibility of height adjustment for the sensor. The distance between the sensor and the bottom of the pot was set at 0.7 m during the first 5 weeks and moved to 1.4 m after Week 5 (Figure 3), since the plants were too Sensors 2020, 20, 435 5 of 21 high to be observed from the initial distance. The higher distance between the camera and the plant makes the plant image smaller, so the leaf area represented by the pixel area became smaller after that modification. In order to model the data, regression analysis using the mean values was performed, so the results could be extrapolated for the next weeks.
commercial system.
The Sequoia Camera has also a sunshine sensor that is used to calibrate the images depending on the sunlight. This makes it possible to compare photos over time, despite variations in light during photo shoots. The sunshine sensor is attached on the upper part to the system, facing the sky and correcting the signal and it also contains a GPS/GNSS module, a magnetometer and inertial measurement system ( Figure 3).
The images were taken with and without a black background placed to help the segmentation. The developed algorithms showed different efficiency in segmentation with or without the background. The pots were tagged with a coloured stamp in order to take the images in the same position during the experiment (Figure 3b). The sensor communicates with a smartphone or computer via Wifi protocol in order to store the acquired images.

Computational System for Image Analysis
The computational system used for the image analysis is composed of four major steps: preprocessing, calculation of vegetation indices, image segmentation, and morphological analysis.

Image Pre-Processing
The Sequoia sensor produces tagged image file format (TIFF) images with a size of 1280 × 960 pixels. The camera was designed to take images from a minimum distance of 30 m to the target [26]. When the images were taken with a shorter distance a displacement occurs between the 4-channel images due to an unexpected geometry between the sensors. In order to correct the displacement of the images, an algorithm that shifts an image by a specified number of pixels in either the x-or ydirection (or both) was used.
The parameters used to define the number of pixels and the direction that each image should be shifted were obtained using a sample image with a referential point and then applied subsequently to all collected data. The red image was arbitrarily selected as the reference, and then the green, near-b) The camera used for image acquisition was a multispectral camera (model Sequoia, Parrot Drones, Paris, France, 2017) which was originally designed for use in agricultural Unmanned Aerial Vehicles (UAV's). Its internal sensor is composed of four spectral bands which register the reflected light coming from the vegetation and can be used to distinguish plant vigor based on reflectance levels [34]. These bands are: green (550 nm wavelength, 40 nm bandwidth), red (660 nm wavelength, 40 nm bandwidth), red edge (735 nm wavelength, 10 nm bandwidth) and near infrared (790 nm wavelength, 40 nm bandwidth). Additionally, a 16-megapixel RGB camera is also fitted into the commercial system.
The Sequoia Camera has also a sunshine sensor that is used to calibrate the images depending on the sunlight. This makes it possible to compare photos over time, despite variations in light during photo shoots. The sunshine sensor is attached on the upper part to the system, facing the sky and correcting the signal and it also contains a GPS/GNSS module, a magnetometer and inertial measurement system ( Figure 3).
The images were taken with and without a black background placed to help the segmentation. The developed algorithms showed different efficiency in segmentation with or without the background. The pots were tagged with a coloured stamp in order to take the images in the same position during the experiment (Figure 3b). The sensor communicates with a smartphone or computer via Wifi protocol in order to store the acquired images.

Computational System for Image Analysis
The computational system used for the image analysis is composed of four major steps: pre-processing, calculation of vegetation indices, image segmentation, and morphological analysis.

Image Pre-Processing
The Sequoia sensor produces tagged image file format (TIFF) images with a size of 1280 × 960 pixels. The camera was designed to take images from a minimum distance of 30 m to the target [26]. When the images were taken with a shorter distance a displacement occurs between the 4-channel images due to an unexpected geometry between the sensors. In order to correct the displacement of the images, an algorithm that shifts an image by a specified number of pixels in either the xor y-direction (or both) was used.
The parameters used to define the number of pixels and the direction that each image should be shifted were obtained using a sample image with a referential point and then applied subsequently to all collected data. The red image was arbitrarily selected as the reference, and then the green, near-infrared (NIR) and red-edge images were shifted in order to be in the same position as the red image. The parameters used for the distances (0.7-1.4 m) are shown in Table 1. The central area of the images was defined as area of interest (AOI) and the consequent processes were used in this area in order to reduce the computational process and optimize the plant recognition and binarization. The images were clipped with AOI-x: 200:900 and AOI-y 200:1000, reducing the size of the image from 960 × 1280 to 700 × 800 squared pixels ( Figure 4).
Sensors 2020, 20, x FOR PEER REVIEW 6 of 21 infrared (NIR) and red-edge images were shifted in order to be in the same position as the red image. The parameters used for the distances (0.7-1.4 m) are shown in Table 1. The central area of the images was defined as area of interest (AOI) and the consequent processes were used in this area in order to reduce the computational process and optimize the plant recognition and binarization. The images were clipped with AOI-x: 200:900 and AOI-y 200:1000, reducing the size of the image from 960 × 1280 to 700 × 800 squared pixels ( Figure 4).

Vegetation Indices
Vegetation indices (VIs) computed from multispectral images are quite simple and effective parameters for the qualitative and quantitative evaluation of vegetation vigour, cover and growth dynamics among other applications. There is a vast implementation of these indices using different airborne and satellite platforms with recent advances using unmanned aerial vehicles (UAV) and tractors. Due to the complexity of the instrumentation platforms, light spectra combinations and resolutions used, there is no unified mathematical expression that attends all applications of VIs [26].
According to the same authors [35], customized algorithms have been studied for several applications combining visible light radiation, mainly the spectral region correspondent with the green region from vegetation, and nonvisible spectra in order to obtain proxy quantifications of the vegetation surface. Therefore, for precise measurement applications, the VIs are optimized and

Vegetation Indices
Vegetation indices (VIs) computed from multispectral images are quite simple and effective parameters for the qualitative and quantitative evaluation of vegetation vigour, cover and growth dynamics among other applications. There is a vast implementation of these indices using different airborne and satellite platforms with recent advances using unmanned aerial vehicles (UAV) and tractors. Due to the complexity of the instrumentation platforms, light spectra combinations and resolutions used, there is no unified mathematical expression that attends all applications of VIs [26].
According to the same authors [35], customized algorithms have been studied for several applications combining visible light radiation, mainly the spectral region correspondent with the green region from vegetation, and nonvisible spectra in order to obtain proxy quantifications of the vegetation surface. Therefore, for precise measurement applications, the VIs are optimized and usually constructed according to the specific application requirements, and a validation procedure is needed, along with customized methodologies.
Eleven VIs retrieved from the literature were adopted to find correlations between the fertilization status of the tomato plants and the spectral response. One of the vegetation indices were also used as a filter parameter for binarization of the plant leaves despite the background (image segmentation).
The VIs were chosen based on the literature and in the available bands present at the sensor. The list of indices, abbreviations, formula and traditional application of the VIs can be observed in Table 2. For use in this study, some closest Sequoia reflectance bands were substituted for the traditional narrowband wavelengths. Index to assess chlorophyll content in areas of complete leaf coverage.
Note: For use in this study, some closest Sequoia reflectance bands were substituted by the traditional narrowband wavelengths as described by [36] and showed in the Table 2.

Plant Extraction (Segmentation)
Digital image processing and computer vision approaches are powerful tools for plant analysis because they allow plant physiological and physical features to be measured non-destructively with high temporal and quantitative resolution [37]. After the image acquisition, the image analysis process starts with the extraction of the numerical data that describes the object in the image. First the background pixels must be separated from the object of interest through a process called object segmentation. The accuracy of this process decreases as the image quality decreases [38]. Multiple Sensors 2020, 20, 435 8 of 21 plant segmentation methods have been proposed to handle the variation among different sensors including the use of adaptative thresholding and machine learning [39].
Other systems were developed using a fixed threshold and image standardization method [38]. Normalized difference vegetation index (NDVI) has been widely used for distinguishing between plant and soil and plant species through satellite imagery. This occurs because the near infrared (NIR) combined with the red image provides a significantly higher reflectance than soil in natural light conditions [40]. However, NDVI is known for presenting problems related to saturation at extreme values. Aiming to design a robust segmentation algorithm that could be used in different soil conditions, the fixed threshold was avoided. Instead, we used an automatic histogram segmentation that extracted the 92% higher reflectance values from the NDVI images when the NDVI values were superior to 0.1. Another filter based in the NIR image was tested but without reliable results ( Figure 5).
Sensors 2020, 20, x FOR PEER REVIEW 8 of 21 Other systems were developed using a fixed threshold and image standardization method [38]. Normalized difference vegetation index (NDVI) has been widely used for distinguishing between plant and soil and plant species through satellite imagery. This occurs because the near infrared (NIR) combined with the red image provides a significantly higher reflectance than soil in natural light conditions [40]. However, NDVI is known for presenting problems related to saturation at extreme values. Aiming to design a robust segmentation algorithm that could be used in different soil conditions, the fixed threshold was avoided. Instead, we used an automatic histogram segmentation that extracted the 92% higher reflectance values from the NDVI images when the NDVI values were superior to 0.1. Another filter based in the NIR image was tested but without reliable results ( Figure  5). The images were converted from unit8 to double type in order to allow the pixels to admit decimal values and calculations. All the image processing analysis was performed in MATLAB 2013b (The MathWorks Inc., Natick, MA, USA). After this step, the bwarefilt function of MATLAB was applied. This function allows extracting all connected components from the binary image, where the area of the objects is in the specified range producing a new binary image containing only the objects that meet the criterion. In this case, the criterion was to extract the larger object of the image (with more connected components). This allowed the system to eliminate noises and unconnected pixels. The resulting binary image (mask) ( Figure 5) was multiplied by the vegetation indices images, in order to extract just the values of VIs present in the plant tissue ( Figure 6).

Morphological Analysis
Nutrient deficiency in tomato plants can cause several symptoms besides the change in reflectance of leaves. It is common to observe changes in the plant growing behavior, specially related with morphological aspects of the leaves and shoots.
Nitrogen deficiency can cause restricted shoot growth and small erect leaflets. Phosphorous deficiency can cause restricted shoot growth and small stiffed curved leaves. Deficiency of potassium can cause scorched and curled symptoms in old leaves. Zinc and Iron deficiencies can cause stunting. Boron and Calcium deficiencies can cause changes in the leaflets making them curved and deformed. Copper deficiency symptoms can be observed in the margins of leaflets and younger leaves that curl into a tube shape, the terminal leaves can become very small, stiff and contorted and the stem growth become somewhat stunted [41]. The images were converted from unit8 to double type in order to allow the pixels to admit decimal values and calculations. All the image processing analysis was performed in MATLAB 2013b (The MathWorks Inc., Natick, MA, USA). After this step, the bwarefilt function of MATLAB was applied. This function allows extracting all connected components from the binary image, where the area of the objects is in the specified range producing a new binary image containing only the objects that meet the criterion. In this case, the criterion was to extract the larger object of the image (with more connected components). This allowed the system to eliminate noises and unconnected pixels. The resulting binary image (mask) ( Figure 5) was multiplied by the vegetation indices images, in order to extract just the values of VIs present in the plant tissue ( Figure 6).
Sensors 2020, 20, 435 9 of 21 from the background. The vegetation indices results were stored by date and treatment in image format seeking a visual confirmation of the segmentation process and a visual representation of the spectral responses (Figures 6-8).
The images were taken in the same position during the weeks in order to observe the growth behavior of the plants in the different treatments. Both aspects could be observed during this time: morphological traits and vegetation indices (Figure 8).

Morphological Analysis
Nutrient deficiency in tomato plants can cause several symptoms besides the change in reflectance of leaves. It is common to observe changes in the plant growing behavior, specially related with morphological aspects of the leaves and shoots.
Nitrogen deficiency can cause restricted shoot growth and small erect leaflets. Phosphorous deficiency can cause restricted shoot growth and small stiffed curved leaves. Deficiency of potassium can cause scorched and curled symptoms in old leaves. Zinc and Iron deficiencies can cause stunting. Boron and Calcium deficiencies can cause changes in the leaflets making them curved and deformed. Copper deficiency symptoms can be observed in the margins of leaflets and younger leaves that curl into a tube shape, the terminal leaves can become very small, stiff and contorted and the stem growth become somewhat stunted [41].
From the binary plant datasets, it is possible to measure plant size, shape, area and other features in an automated way and correlate plant phenotypes with experimental treatments [42,43].
Several other morphological properties can be extracted from binary images using the regionprops function of MATLAB. This function uses the distribution of the pixels to calculate shape parameters and has been recently used for image analysis in medical studies [44], industry [45], plant pathology [46] and plant growth analysis.
These parameters can provide reliable information about the changes in the plants' area, perimeter, shape and growth behaviour.
Not all parameters resulting from the regionprops function were used, and the description of the morphological properties used for the tomatoes plants analysis can be observed in Table 3. Table 3. Morphological properties calculated for the tomato plants at early stages images with different levels of organic fertilization.

Morphological Property Description
Area Actual number of pixels in the region, returned as a scalar.

Convex Area
Number of pixels in the image that specifies the convex hull, with all pixels within the hull filled in (set to on), returned as a binary image (logical). The image is the size of the bounding box * of the region.

Eccentricity
Eccentricity of the ellipse that has the same second-moments as the region, returned as a scalar. The eccentricity is the ratio of the distance between the foci of the ellipse and its major axis length. The value is between 0 and 1. (0 and 1 are degenerate cases. An ellipse whose eccentricity is 0 is a circle, while an ellipse whose eccentricity is 1 is a line segment.)

Diameter Equivalent
Diameter of a circle with the same area as the region, returned as a scalar.
Euler Number Number of objects in the region minus the number of holes in those objects, returned as a scalar.

Extent
Ratio of pixels in the region to pixels in the total bounding box *, returned as a scalar. Computed as the area divided by the area of the bounding box *.

Filled Area
Number of on pixels in filled image, returned as a scalar.

Orientation
Angle between the x-axis and the major axis of the ellipse that has the same second-moments as the region, returned as a scalar. The value is in degrees, ranging from −90 degrees to 90 degrees.
Major Axis Length Length (in pixels) of the major axis of the ellipse that has the same normalized second central moments as the region, returned as a scalar.
Minor Axis Length Length (in pixels) of the minor axis of the ellipse that has the same normalized second central moments as the region, returned as a scalar Perimeter Distance around the boundary of the region returned as a scalar. This function computes the perimeter by calculating the distance between each adjoining pair of pixels around the border of the region Solidity Proportion of the pixels in the convex hull that are also in the region, returned as a scalar. Computed as area/convex area * Bounding box: Smallest rectangle containing the region.

Statistical Analysis
The experiment used randomized blocks, with four treatments and six repetitions. Each plot was composed of six pots. The statistical analysis carried out through the ANOVA and Tukey's HSD test with 99% of confidence for average comparisons using the statistical and machine learning toolbox from MATLAB software. ANOVA and Tukey's HSD test were also performed with 95% confidence, but since there it was no difference in terms of significant factors, the 99% was chosen for this study. For the significant predictor factors, a functional boxplot and a regression analysis (using the week means) was performed in the same software.

Results and Discussion
Images were collected weekly (one image per plant) during the course of eight weeks. The segmentation of the plant tissue from the background was needed to quantify the spectral response from each plant and to calculate the mean values of the evaluated parameters.
This process also enabled the morphological analysis of the plant growth behavior. Multiple segmentation codes were tested to remove the background, and the code based on the normalized difference vegetation index (NDVI) was the most efficient.
All the other bands (red, red edge and green) alone were not able to differentiate the plant tissue from the background. The vegetation indices results were stored by date and treatment in image format seeking a visual confirmation of the segmentation process and a visual representation of the spectral responses (Figures 6-8).   The images were taken in the same position during the weeks in order to observe the growth behavior of the plants in the different treatments. Both aspects could be observed during this time: morphological traits and vegetation indices (Figure 8).
The average values of the different vegetation indices per plant were calculated and compared during each week of the experiment. Parameters that showed significant difference and correspondent p-values are displayed below (Table 4). Aiming to identify the most reliable parameters to predict the fertilization level of organic tomato, boxplots were created with the whole amount of data, without week distinction. After that, a comparison test (Tukey's HSD) was performed to find which parameters are more robust to use in an automatic and nondestructive fertilization level analysis model.
The graphical representation of the boxplots and Tukey's comparison test can be seen in Figure 9. In the case of the morphological aspects, some of the parameters presented significant results with distinction between the control (T0) and the treatments (T1, T2, and T3) (Figure 9). On the other hand, none of the vegetation indices presented a significative difference (with a 99% or 95% level of confidence) to predict the nutritional level of organic tomato plants independent of the week after transplant ( Figure 10).
The results presented in Table 4, as well as Figures 9 and 10 indicate that the morphological parameters are more related to the fertilization level of tomato plants at early stages than the spectral responses. These parameters were selected to create regression models in order to predict the nutritional status of the tomato crop.
The regression models were created using the MATLAB function curvefit, the mean values of the control (T0), and the recommended dose (T3), for the first five weeks and excluding outliers. The selected parameters were: area, filled area, convex area, perimeter, equivalent diameter, major axis length and minor axis length (Table 5). Sensors 2020, 20, x FOR PEER REVIEW 13 of 21 Figure 9. Boxplots and Tukey Test with 99% of probability using all data (50 days) representing morphological parameters of tomato plants with different levels of organic fertilization. T0, T1, T2, and T3 were the four groups of plants treated with different amounts of fertilizer and vermicompost (see Section 2).   Seeking for a more robust, simple and reliable prediction model, functional boxplots were created using the mean values of each plant T0 and T3 for the first four weeks and plotted together. The functional boxplots represent the mean data of each week without the outliers with an interval of confidence of 95% ( Figure 11).

Discussion
The system used for image acquisition showed to be effective for tomato plants at early stages. The same system could be adjusted for tomato plants at late stages, as well as for other plant species and seedlings in different stages (depending on the height adjustments). The pre-processing technique showed to be adequate for the overlay of the Parrot Sequoia images in the different heights tested 0.70 m and 1.40 m.
Plant segmentation is the cornerstone of the image analysis at plant level. Diverse authors used multispectral images to access the nutritional level of the plants. These works are usually focused in small zones or areas, or in multiple plants in a parcel as the minimum management unit. In order to perform precise fertilization at the plant level, the segmentation and extraction of the plant of interest must be also precise. The algorithm developed using the normalized difference index was able to distinguish plant tissue from soil, pot, floor and metal platform extracting the whole plant apart of the background (Figures 4-6). The use of NDVI for vegetation extraction and monitoring is widely used for satellite imagery [47][48][49] and this study proves that can also be used with high-resolution plant images. The code was applied in different conditions with good efficacy in extract the plant. The visual inspection was needed to guarantee that the plants were well-segmented and also provides the distribution of the vegetation indices intensity through the leaf tissue. The segmentation also allows morphological analysis during the time, and this system can be useful for digital phenotyping projects. However, the code can present limitations for morphological analysis of plants with presence of senescent leaves since the reflectance of yellow and/or brown leaves are different from the active vegetative tissue.
The vegetation indices used in this study did not show significant differences between the treatments during the initial stages of the plant development, but there were significant differences in the 6th week after the transplant between T0 and T2 and T0 and T1. Among all the studied vegetation indices, the NDVI, GRVI, OSAVI, SR, and MSR were sensible to these differences. Different results were found by [50], that using the crop circle ACS 470 and the SPAD 502 sensors concluded that the vegetation indices based on the red band (NDVI and red vegetation Index) were the best predictors for nitrogen levels in tomatoes maintaining the relationships during the crop.
According to the authors [50], the SPAD readings, the GNDVI, and the GVI were good predictors of the nitrogen status of the crop in the beginning of the cycle, with low accuracy in the later part of the crop cycle. Using the same sensor (Circle ACS 470) and guided sampling method, [51] developed a yield prediction map of tomato crop. Their study used the NDVI index as a predictor parameter and the correlation between the vegetation index and the yield were of 0.67 and 0.71 in two different regions of Spain. Another study conducted by [52] showed the possibility to use vegetation indices to improve the use of nitrogen fertilization in tomatoes cultivated in protected environment. Applying the spectral sensors ASD Field Spec HandHeld 2, and the Minolta SPAD 502, the authors compared the spectral response of tomato plants in the reference plot and with lower levels of nitrogen fertilization.
The use of NDVI and the SPAD measurements as a nitrogen level parameter allowed to reduce significantly the use of nitrogen fertilization maintaining the same yield levels and fruit quality compared to the reference plot [52]. On the other hand, during three years of field trials with different levels of nitrogen fertilization [53] evaluated vegetation indices of tomato plots through the use of a multispectral radiometer (MSR-87, Cropscan Inc., Rochesters, MA, USA). In their study [53], sixteen vegetation indices were analyzed and the NDVI did not show a strong correlation with nitrogen concentration on leaves. It was concluded that the best vegetation index for nitrogen evaluation in tomato was the relation between NIR and the reflectance present in the range of 560 nm (NIR/R560). Some factors can be responsible for the divergence observed between literature results.
The first factor is the type of fertilizer, in all the previous studies, mineral nitrogen fertilization was used. This type of fertilizer is more soluble and less stable and can be absorbed and translocated faster than the organic amendments. The crop stage during the measurements can also be a factor. In the previous studies, the measurements started after the 29th day after transplant where the plants were developed more, and the deficiencies symptoms were more present. The type of sensor is also an important factor because this was the first study to use the Parrot Sequoia sensor at plant level to monitor nutrient deficiency in tomatoes. [54], showed the differences in measurements in different commercial spectral sensors for nitrogen management, and among the studied factors was concluded that the measuring distance, the device temperature and the light intensity can influence the performance of the sensors.
The most interesting results were obtained through the analysis of the morphological traits of the tomato plants. Seven of the twelve aspects analyzed showed significant difference between the control and the fertilized treatments. Only three weeks after the transplant, the plants started to exhibit a significant difference in the plant area, convex area, filled area, minor axis length and equivalent diameter discriminating T3 and T2 of T1 and the control. In the same period, the perimeter was significantly different between the recommended dose T3 and the control, without difference between T2 and T1. The length of the major axis could discriminate between the control against the treatments but not between the treatments already in the third week. In the following weeks, from the fourth week to the end of the study, all the seven morphological traits (area, convex area, filled area, equivalent diameter, perimeter, minor axis length and minor axis length) showed significant differences between the control and the treatments, except for the major axis length in the seventh week after transplant.
These results are in agreement with the results described in the literature. Mooy et al. [55] showed the positive influence of soluble organic fertilizer in tomato plant height and the number of leaves. Other studies also demonstrated that the nitrogen and phosphorous fertilization influences the fresh shoot weight, plant height, stem diameter, leaf number, and leaf area of tomato seedlings [56][57][58][59]. In field trials, the nitrogen fertilization produced an increase in the leaf area index and in the above-ground weight of tomato plants [60] and the use of NPK induce an increase in plant height and in the number of leaves [61]. Besides the slight difference in the morphological parameters between the fertilized treatments T1, T2, and T3, there was no significant difference with 1% of probability, because, for the functional boxplot analysis, the treatment with the recommended dose (T3) was chosen for comparison with the control (T0).
The use of functional boxplots as a predicting model is interesting because they can exclude outliers and can be used as a threshold parameter [62]. As can be observed, the functional boxplots displays the range of the observed values until the 4th week after transplant and not until the Week 7. At Week 5 the plants were too high to be observed with 0.70 m height as was planned for the robotic prototype, so the camera bracket was set for a 1.4 m height. This adjustment in the camera distance cause differences in the relationship between morphological parameters and the number of pixels.

Conclusions
The computer vision methodology developed in this article was suitable for monitoring the development of tomato var. Mina F1 at early stages with different levels of organic fertilization in a protected environment. The methodology involves consecutive steps, and using an automatic and non-destructive approach provides several morphological and spectral aspects of the plants that can be used in a robotic platform.
After the pre-processing, the NDVI threshold showed to be crucial and effective for plant extraction. Among the several spectral (vegetation indices) responses observed, the NDVI, OSAVI, simple ratio, GRVI, and modified simple ratio were the only ones that showed a small difference between treatments in the sixth and seventh weeks after transplant, but without practical applications.
Studying the evolution of a crop by using multiples images allowed to extract several morphological features. The automatic morphological analysis was able to distinguish between the control and the fertilized treatments from the third week after transplant until the end of the experiment with 99% of confidence according to the Tukey's honest significant differences test.
Among the studied morphological parameters, area, convex area, filled area, perimeter, major axis length, minor axis length and equivalent diameter showed to be useful to distinguish between the control and fertilized treatments. Convex area was the parameter that presented the highest difference between control and the fertilized plants according to the functional boxplot analysis.
This was the first study to relate the use of morphological features extracted automatically with the organic fertilization level in vegetable plants. The image analysis system was able to distinguish plant tissue from soil, floor, and metal structures. Therefore, most of the image processing steps can be extrapolated to other vegetable cultures, although future analysis is necessary to evaluate its robustness. Moreover, the Parrot Sequoia camera fulfilled the requirements to be used in the robotic platform and besides be originally developed for aerial images this study proved the capacity of its use at short distances.
Additionally, the same methodology can be used for analyzing other aspects related with morphology features in tomato plants, such as disease severity, hydric and salinity stress or any environmental modification that can cause aboveground changes in the morphology or spectral responses of vegetable plants. This technique has potential to be used in many branches of plant health studies and provides an innovative and automatic way to measure the severity of abiotic disorders.
In a subsequent investigation, we are developing an extended version of the methodology to be applied using the robotic platform at the field level. Trials are being conducted in an experimental line with different vegetable crops to study the best image acquisition interval in a scenario with movement. Other trials using the robotic platform are being conducted in an experimental field area with different crops such as onion, pumpkins, tomatoes, and zucchini with different levels of fertilization. In these trials, the multispectral camera is being integrated with laser sensors (Lidar) for a multiscale analysis of the crops. Funding: Transnational funding bodies, being partners of the H2020 ERA-net project, CORE Organic Cofund, and the cofund from the European Commission.