3.5.3.3 General discussion

Even if the relief of the scene is not taken into account due to the used deformation model,
in our case, with flat ground, no difference arise. However, more complex deformation models
[Lombaert et al., 2012, Bookstein, 1989] could be used to improve the remaining error. But
could also, in some case, create large angular deformations caused by the proximity of key-
points, of course, it's possible to filter these key-points, which would also reduce the overall
accuracy.

Further research can be performed on each parameter of the feature extractors, for those
who need specific performance (time/precision), we invite anyone to download the dataset and
test various combinations. Otherwise feature matching can be optimized, at this stage, we use
brute-force matching with post filtering, but a different implementation that fulfill our spatial
properties should greatly enhance the number of matches by reducing false positives.

3.5.4 Conclusion

In this work, the application of different techniques for multi-spectral image registration was
explored using the Airphen camera. We have tested nine type of key-points extractor (ORB,
GFTT, AGAST, FAST, AKAZE, KAZE, BRISK, SURF, MSER) at different heights and the
number of control points obtained. As seen in the method, the most suitable method is the
GFTT (regardless of modalities 1, 2 or 3) with a significant number of matches 150 – 450 and
a reasonable calculation time 1.17 s to 3.55 s depending on the modality.

Furthermore, the best spectral reference was defined for each method, for example 570 nm
for GFTT. We have observed a residual error of less than 1 px, supposedly caused by the
difference of sensors nature (spectral range, lens).

Acknowledgments

The authors acknowledge support from European Union through the project H2020 IWM-
PRAISE<fnref n="3" /> (Integrated Weed Management : PRActical Implementation and Solutions for Eu-
rope) and from ANR Challenge ROSE through the project ROSEAU<fnref n="4" /> (RObotics SEnsorimotor
loops to weed AUtonomously).

We would like to thanks Combaluzier Quentin, Michon Nicolas, Savi Romain and Masson
Jean-Benoit for the realization of the metallic gantry that help us positioning the camera at
different heights.

<footnote n="3">3. https://iwmpraise.eu/</footnote>

<footnote n="4">4. http://challenge-rose.fr/en/projet/roseau-2/</footnote>
### 3.6 Conclusion de chapitre

La partie pré-traitement de notre chaîne fonctionnelle correspond à trois étapes. Premièrement, la calibration en hauteur de la caméra a permis d'obtenir la zone que nous souhaitons observer. Deuxièmement, les corrections géométriques des lentilles nous ont permis de limiter les erreurs entre les capteurs. Et troisièmement, le recalage des images multispectrales a permis de superposer les différentes images issues des six capteurs de la caméra Airphen.

Les deux premières corrections ont été présentées durant le chapitre précédent. Tandis que ce chapitre était dédié au recalage des images multispectrales. Cette étude a fait l'objet d'un article de conférence (VISAPP), présenté en Février 2020 à Malte (Rép. de Malte). La méthode proposée offre les performances nécessaires à l'exploitation des images multispectrales acquises en proxidétection, durant nos expérimentations en micro-parcelles et pendant les challenges RoSE. Toutes les images multispectrales acquises avec la caméra Airphen ont donc été transformées à l'aide de cette méthode avant leur annotation et utilisation.

Ce chapitre marque également la fin de la partie pré-traitement des données. Maintenant, que nous sommes capables d'exploiter les données images, nous allons nous intéresser à leur utilisation. La seconde étape de notre chaîne fonctionnelle, vise à détecter les surfaces et objets visibles d'une scène agronomique. L'approche traditionnelle s'appuie sur l'utilisation d'un indice de végétation sur lequel est appliqué une méthode de segmentation, afin d'obtenir les éléments à discriminer. Cette approche a été conservée, mais de nouveaux indices basés sur l'apprentissage profond sont présentés dans le chapitre suivant 4. Concernant la segmentation des objets, une nouvelle méthode a également été proposée. Elle concerne spécifiquement la détection des feuilles. Elle est présentée dans le chapitre 5. Pour rappel, la figure ci-dessous illustre notre chaîne fonctionnelle et où nous nous situons.

<figure><img src="image_102.png" /></figure>

Rappel de l'architecture de notre chaîne fonctionnelle ici la première étape en vert est terminée et en jaune l'étape suivante étudiée
<figure><img src="image_103.png" /></figure>
# 4 Indices de végétation

Dans le chapitre précédent, nous avons vu comment pré-traiter les données en vue de leur exploitation. La seconde étape de notre chaîne de traitement vise à répondre à la question suivante : comment discriminer la végétation du reste de la scène ? Pour répondre à cette question, l'approche classique repose sur l'utilisation d'indices de végétation, sur lesquels une méthode de segmentation simple est appliquée, telle qu'un seuillage. Comme nous l'avons vu dans le chapitre 2, la définition de tels indices (parmi lesquels le NDVI est le plus connu) repose essentiellement sur des connaissances empiriques [Bannari et al., 1995] (section 2.2.3). Ces indices utilisent souvent des données calibrées en réflectance, censées garantir leur stabilité.

Or de nombreux facteurs influencent de tels indices de télédétection. La calibration radiométrique et le capteur sont deux facteurs importants dans l'estimation de ces indices. [Lee et al., 2021] montre leurs effets au cours de la journée sur le NDVI. Cet indice peut varier d'un facteur de 20% après correction de ces deux paramètres. D'autres facteurs, comme les conditions atmosphériques, l'humidité de la végétation, l'humidité du sol, la couverture végétale globale, l'ombre et les différences de type de sol, jouent aussi des rôles importants dans l'estimation du couvert végétal. Par exemple, la lumière réfléchie par le sol et l'ombre obscurcissant partielle-ment une image, peuvent influencer les valeurs du NDVI dans une large mesure [Zhang et al., 2015], rendant difficile l'utilisation d'un tel indice.

Ces limitations ont été la raison du développement de nombreux indices de végétation, ajustés au sol (par exemple, Soil-Adjusted Vegetation Index, Modified Soil-Adjusted Vegetation Index) ou à l'atmosphère (Atmospherically Resistant Vegetation Index) [Bannari et al., 1995]. Ces modèles, plus robustes, restent néanmoins fortement limités et influencés par tous ces facteurs. L'utilisation d'un seuillage fixe est donc peu pertinente car il faut s'attendre à des valeurs différentes pour une même zone d'observation au fil du temps, en raison des changements environnementaux affectant l'observation (saisonniers, phénologiques, météorologique, ...).

Pour augmenter la performance de cette étape de traitement, certains utilisent une méthode de classification faisant appel à plusieurs indices usuels [Shishir and Tsuyuzaki, 2018, Gée et al., 2021]. D'autres en revanche, utilisent des algorithmes génétiques afin d'apprendre l'équation idéale permettant une segmentation optimale [Kabiri et al., 2011, Albarracín et al., 2020]. Dans les deux cas, cela revient à prendre en considération une plus large quantité d'information, ce qui augmente le potentiel de discrimination. En revanche, l'optimisation de ces indices utilise des méthodes, soit coûteuses en temps, soit ne garantissant pas un optimum, les rendant peu applicables à de larges bases de données.

Face à ce constat, une nouvelle méthodologie est proposée dans ce chapitre pour l'apprentissage de tels indices. Elle repose sur l'apprentissage d'approximateurs de fonctions. Ainsi, on cherche une fonction qui peut décrire toutes les autres. Cette approche peut être mise en œuvre en apprentissage profond. Ce qui permet d'exploiter le processeur graphique, permettant une évaluation rapide de larges bases de données et une optimisation d'équations (ou de réseau de neurones) par le mécanisme de rétropropagation. L'apprentissage profond peut-être vu comme une méthode de régression générique. Les indices qui ont été définis permettent de s'abstraire des conditions d'acquisition et de la calibration.
# DeepIndices : A new approach to remote sensing indices based on optimization and approximation of functions through deep-learning, with application to vegetation surfaces over uncalibrated data

<figure><img src="image_104.png" /></figure>

**Authors** Jehan-Antoine Vayssade, Jean-Noël Paoli, Christelle Gée and Gawain Jones

**Publication** 09/06/2021 – Remote Sensing Journal – Volume 13 Issue 12 Pages 2261

**Abstract** The form of a remote sensing index is generally empirically defined, whether by choosing specific reflectance bands, equation forms or its coefficients. These spectral indices are used as preprocessing stage before object detection/classification. But no study seems to search for the best form through function approximation in order to optimize the classification and/or segmentation. The objective of this study is to develop a method to find the optimal index, using a statistical approach by gradient descent on different forms of generic equations. From six wavebands images, five equations have been tested, namely : linear, linear ratio, polynomial, universal function approximator and dense morphological. Few techniques in signal processing and image analysis are also deployed within a deep-learning framework. Performances of standard indices and DeepIndices were evaluated using two metrics, the dice (similar to f1-score) and the mean intersection over union (mIoU) scores. The study focuses on a specific multispectral camera used in near-field acquisition of soil and vegetation surfaces. These DeepIndices are built and compared to 89 common vegetation indices using the same vegetation dataset and metrics. As an illustration the most used index for vegetation, NDVI (Normalized Difference Vegetation Indices) offers a mIoU score of 63.98 % whereas our best models gives an analytic solution to reconstruct an index with a mIoU of 82.19 %. This difference is significant enough to improve the segmentation and robustness of the index from various external factors, as well as the shape of detected elements.

**Keywords** image; precision agriculture; spectral indices; multi-spectral; deep-learning; vegetation segmentation
**4.1 Introduction**

An important advance in the field of earth observation is the discovery of spectral indices, they have proved their effectiveness in surface description. Several studies have been conduc- ted using remote sensing indices, often applied to a specific field of study like evaluations of vegetation cover, vigor, or growth dynamics [Jinru and Su, 2017, Mezera Jiří, 2019, Tanrıverdi et al., 2020, Elbeltagi et al., 2021] for precision agriculture using multi-spectral sensors. Some spectral indices have been developed using RGB or HSV color space to detect vegetation from ground cameras [Lee et al., 2020, Milioto et al., 2017, Hassanein et al., 2018]. Remote sensing indices can also be used for other surfaces analysis like water, road, snow [Dixit et al., 2019] cloud [Zhai et al., 2018] or shadow [Henrich et al., 2009].

There are two main problems with these indices. Firstly they are almost all empirically defined, although the selection of wavelengths comes from observation, like NDVI for vegetation indices. It is possible to obtain better spectral combinations or equations to characterize a surface with specific acquisitions parameters. It is important to optimize upstream the index, as the data transformation leads to a loss of essential information and features for classification [Zhang et al., 2015]. Most studies have tried to optimize some parameters of existing indices. For example, an optimization of NDVI (NIR - Red)/(NIR + Red) was proposed by [Gitelson, 2004] under the name of WDRVI (Wide Dynamic Range Vegetation Index) (αNIR - Red)/(αNIR + Red). The author tested different values for α between 0 and 1. The ROC curve was used to determine the best coefficient for a given ground truth. Another optimized NDVI was designed and named EVI (Enhanced Vegetation Index). It takes into account the blue band for atmospheric resistance by including various parameters G(NIR - Red)/(NIR + C₁Red - C₂Blue + L), where G, L are respectively the gain factor and the canopy background adjustment, in addition the coefficients C₁, C₂ are used to compensate for the influence of clouds and shadows. Many other indices can be found in an online database of indices (IDB: www.indexdatabase.de accessed August 10, 2019) [Henrich et al., 2009] including the choice of wavelengths and coefficients depending on the selected sensors or applications. But none of the presented indices are properly optimized. Thus, in the standard approach, the best index is determined by testing all available indices against the spectral bands of the selected sensor with a Pearson correlation between these indices and a ground truth [Liu et al., 2017, Kokhan and Vostokov, 2020]. Furthermore, correlation is not the best estimator because it neither considers the class ratio nor the shape of the obtained segmentation and may again result in a non-optimal solution for a specific segmentation task. Finally, these indices are generally not robust because they are still very sensitive to shadows [Zhang et al., 2015]. For vegetation, until recently, all of the referenced popular indices were man-made and used few spectral bands (usually NIR, red, sometimes blue or RedEdge).
The second problem with standard indices is that they work with reflectance-calibrated data. Three calibration methods can be used in proximal-sensing. (i) The first method use an image taken before acquisition containing a color patch as a reference [Yahui et al., 2019, Minařík et al., 2019], and is used for correction. The problem with this approach is that if the image is partially shaded, the calibration is only relevant on the non-shaded part. Moreover ideally the reference must be updated to reduce the interference of weather change on the spectrum measurement, which is not always possible since it's a human task. (ii) An other method is the use of an attached sunshine sensor [Gilliot et al., 2018], which also requires calibration but does not allow to correct a partially shaded image. (iii) The last method is the use of a controlled lighting environment [Chebrolu et al., 2017, Wu et al., 2019], e.g., natural light is suppressed by a curtain and replaced by artificial lighting. All of these approaches are sometimes difficult to implement for automatic, outdoor use, and moreover in real time like detecting vegetation while a tractor is driving through a crop field.

In recent years, machine learning algorithms have been increasingly used to improve the definition of presented indices in the first main problem. Some studies favor the use of multiple indices and advanced classification techniques (RandomForest, Boosting, DecisionTree, etc.) [Oldeland et al., 2010, Peña-Barragán et al., 2011, Nguyen-Robertson et al., 2012, Shishir and Tsuyuzaki, 2018, Elbeltagi et al., 2021, Lu et al., 2021]. Another study has proposed to optimize the weights in an NDVI equation form based on a genetic algorithm [Kabiri et al., 2011] but does not optimize the equation forms. An other approach has been proposed to automatically construct a vegetation index using a genetic algorithm [Albarracín et al., 2020]. They optimize the equation forms by building a set of arithmetic graphs with mutations, crossovers and replacements to change the shape of each equation during learning but it does not take into account the weights, since it's use calibrated data. Finally, with the emergence of deep learning, current studies try to adapt popular CNN architectures (UNet, AlexNet, etc.) to earth observation applications: [Lv et al., 2018, Gaetano et al., 2018, Fu et al., 2018, Kussul et al., 2017].

However there is no study that optimize both the equation forms and spectral bands weights. The present study explicitly optimize both of them by looking for a form of remote sensing indices by learning weights in functions approximators. These functions approximators will then reconstruct any equation forms of the desired remote sensing index for a given acquisition system. To solve the presented second problem, this study evaluates the functions approximators on an uncalibrated dataset containing various acquisition conditions. This is not a common approach but can be found in the literature [Bajwa and Tian, 2002, Bareth et al., 2016]. This will lead to creating indices that do not require data calibration. The deep learning framework has been used as a general regression toolkit. Thus, several CNN function approximators architectures are proposed. DeepIndices is presented as a regression problem, which is totally new, as is the use of signal and image processing.
## 4.2 Material and Data

### 4.2.1 Instrument Details

The images were acquired with the Airphen (Hyphen, Avignon, France) six-band multi-spectral camera (Figure 4.1). This is a multi-spectral scientific camera developed by agronomists for agricultural applications. It can be integrated into different types of platforms such as drones, phenotyping robots, etc.

<figure><img src="image_105.png" /><figcaption>**Figure 4.1** – AIRPHEN camera composed of 6 sensors</figcaption></figure>

The camera has been configured using the 450/570/675/710/730/850 nm bands with a 10 nm FWHM respectively denoted from $λ_0$ to $λ_5$. These spectral bands have been defined by a previous study [Louargant et al., 2017] for crop/weed discrimination. The focal length of each lens is 8 mm. The raw resolutions for each spectral band is 1280 × 960 px with 12 bit precision. Finally, the camera is equipped with an internal GPS antenna.

### 4.2.2 Image Dataset

The dataset were acquired on the site of INRAe in Montoldre (Allier, France) within the framework of the “RoSE challenge” founded by the French National Research Agency (ANR) and in Dijon (Burgundy, France) within the site of AgroSup Dijon. Images of bean and corn, containing various natural weeds (yarrows, amaranth, geranium, plantago, etc) and sowed ones (mustards, goosefoots, mayweed and ryegrass) with very distinct characteristics in terms of illumination (shadow, morning, evening, full sun, cloudy, rain, etc) were acquired in top-down view at 1.8 m from the ground. The Table 4.1 synthesis the dataset.

<table><caption>**Table 4.1** – Acquisition sources and global illumination.</caption><thead><tr><th>Source</th><th>Year</th><th>Corn</th><th>Bean</th><th>Illumination</th></tr></thead><tbody><tr><td>Dijon</td><td>2019</td><td>-</td><td>9</td><td>full sun, evening</td></tr><tr><td>Montoldre</td><td>2019</td><td>20</td><td>22</td><td>shadow, sunny, cloudy</td></tr><tr><td>Montoldre</td><td>2020</td><td>18</td><td>22</td><td>morning, cloudy, rainy</td></tr><tr><td>Total</td><td></td><td>38</td><td>53</td><td>=91</td></tr></tbody></table>

Manual annotation takes about 4 h per image to obtain the best quality of ground truth, which is necessary for use in regression algorithms. Thus the ground truth size is small and defined with very distinctive illumination condition. To simulate light variations effect on the ground truth images a random brightness (20%) and a random saturation (5%) are added to each spectral band during the training phase. As illustration the Figure 4.2 shows a false color reconstruction of corn crop in the field with various weeds and shadows on the corners of the image (not vignetting).
<figure><img src="image_106.png" /><figcaption>Figure 4.2 – False color in the left and corresponding manual ground truth on the right.</figcaption></figure>

### 4.2.3 Data Pre-Processing

#### 4.2.3.1 Images Registration

Due to the nature of the camera (Figure 4.1), a spectral band registration is required and performed with a registration method based on previous work [Vayssade et al., 2020] (with a sub-pixel registration accuracy). The alignment is refined in two steps, with (i) a rough estimation of the affine correction and (ii) a perspective correction for the refinement and accuracy through the detection and matching of key points. The result shows that *GoodFeatureToTrack* (GFTT) algorithm is the best key-point detector considering the $\lambda_{570}$ nm band as spectral reference for the registration. After the registration, all spectral images are cropped to 1200 × 800 px and concatenated to channel-wise denoted $\lambda$ where each dimension denoted $\lambda_d$ refers to each of the six spectral bands.

#### 4.2.3.2 Images Normalization

Spectral bands inherently have a high noise associated with the CCD sensor, which is a potential problem during normalization [Khanna et al., 2017]. To overcome this effect, 1% of the minimum and maximum signal is suppressed by calculating the quantiles, the signal is clipped on the given range and each band is rescaled in the interval [0, 1] using min-max normalization to obtain $\rho_d$:

$$ \rho_d = 0 \le \frac{\lambda_d - \min(\lambda_d)}{\max(\lambda_d) - \min(\lambda_d)} \le 1 \qquad (4.1) $$

The method also reduces the lighting variation. According to [Blackburn and Vignola, 2012] a little variation is observed in the spectral correction factors between clear and cloudy days. Thus, the correction has a limited impact on the scaling factor and should be managed by this equation. However, the displacement factor could not be estimated, thus the output images are not calibrated in reflectance.

#### 4.2.3.3 Enriching Information

In order to enrich the pool of information, some spectral band transformations are added, which allow to take into account spatial gradients and spectral mixing [Milioto et al., 2017] in the image. The choice is oriented towards seven important information in different respects:

The standard deviation between spectral band, noted $\rho_{std}$ can help to detect the spectral mixture. For example between two different surface like ground and leave which have opposite spectral radiance the spectral mixing make a pixel with linear combination, thus the standard deviation tend to zero [Louargant et al., 2017]. Three Gaussian derivatives on different orientations are computed $G_{xx}$, $G_{xy}$ and $G_{yy}$ over the standard deviation $\rho_{std}$ which give an important spatial information about the gradients breaks corresponding to the outer limits of
surfaces. These Gaussian derivatives are computed with a fixed Sigma = 1. The Laplacian computed over the standard deviation ρ<sub>std</sub>, the minimum and maximum eigenvalues of the Hessian matrix (obtained by Gaussian derivation Gxx, Gxy and Gyy), also called ridge are included. These transformations sould improve the detection of fine elements [Lin et al., 2014] such as monocotyledons for vegetation images.

All these transformations are concatenated to the channel-wise normalized spectral band input and build the final input image. In total seven transformations are added to the six spectral images for a final image of 13 channels, which will probably help the convergence.

4.2.4 Training and Validation Datasets

The input dataset is composed by spectral images I of size 1200 × 800 × 13 (or 6 if the “Enriching information” part is disabled) and a manual ground truth p of size 1200 × 800 × 1 where p ∈ {0, 1}. The desired output p̂ is a probability vegetation map of size 1200 × 800 × 1 where p̂ ∈ [0, 1]. This input dataset is randomly split into two sub-sets respectively training (80%) and validation (20%). All random seed is fixed at the start-up to keep the same training/validation dataset across all trained models which help to compare them. Keeping the same random seed also results in the same starting point between different new runs, making results reproducible on the same hardware.

4.3 Methodology

4.3.1 Existing Spectral Indices

From the indices database, 89 vegetation indices have been identified (Table 4.2) as com-
patible with the wavelengths used in this study (as near as possible), they will be tested and
compared to the designed DeepIndices. Five forms of simple equations have been extracted from
this database (a wide variety of indices are derived from these forms, generally a combination
of 2 or 3 bands) :

$$
\text{band reflectance} = \rho_i \tag{4.2}
$$

$$
\text{two bands difference} = \rho_i - \rho_j \tag{4.3}
$$

$$
\text{two bands ratio} = \rho_i \div \rho_j \qquad (4.4)
$$

$$
\text{normalized difference two bands} = (\rho_i - \rho_j) \div (\rho_i + \rho_j) \quad (4.5)
$$

$$
\text{normalized difference three bands} = (2\rho_i - \rho_j - \rho_k) \div (2\rho_i + \rho_j + \rho_k) \quad (4.6)
$$

By analyzing these five equations we can synthesize them into two generic equations (Linear
combination and Linear ratio) which take into account all spectral bands. Three other models
can generalize any function : the polynomial fitting, the continuous function approximations
by Taylor development, and the piecewise continuous function approximations trough morpho-
logical operators. These forms are interesting to optimize because they can approximate any
function. This optimization will lead to automatically defining new indices (DeepIndices). The
following subsections present these different models.
## 4.3.2 Deepindices : Baseline Models

### 4.3.2.1 Linear Combination

To synthesize Equations (4.2) and (4.3), a simple linear equation such as $y = \sum_{d=0}^{N} \alpha_{d} \rho_{d}$ can be defined. This equation can be generalized to the 2D domain using a 2D convolution allowing consider the neighboring pixels. For a pixel at the position $[i, j]$ the convolution is defined by :

$$y[i, j] = \sum_{d=0}^{D} \sum_{h=0}^{N} \sum_{w=0}^{N} \rho_{d}[i - N/2 + h, j - N/2 + w] * H[h, w, d] \quad (4.7)$$

where $H$ defines neighborhood weights (corresponding to $\alpha_{i}$). $D$ is the number of dimensions (6 spectral bands + 5 transformations) and $N$ is the kernel size. The linear combination is given by $N = 1, D = 12$. The kernel weights are initialized by a truncated normal distribution centered on zero [Jang and Son, 2019], weights are updated during the training of the CNN through back-propagation and unnecessary bands should be set to zero. The interesting part is that increasing the kernel size $N$ allows to take into account the neighborhood of a pixel and should estimate more accurately the spectral mixing [Louargant et al., 2017]. Figure 4.3 shows the corresponding network.

<figure><img src="image_107.png" /><figcaption>Figure 4.3 – Linear combination model.</figcaption></figure>

### 4.3.2.2 Linear Ratio

To generalize Equations (4.4)–(4.6), a simple model based on the division of two linear combination is set. In the same way, this form is generalizable to the 2D domain and then corresponds to two 2D convolutions, one for the numerator, the other for the denominator. When the denominator is zero, the result is set to zero as well, to leverage the “not a number” output. The Figure 4.4 shows the corresponding network.

<figure><img src="image_108.png" /><figcaption>Figure 4.4 – Linear ratio model.</figcaption></figure>

### 4.3.2.3 Polynomial

According to the Stone-Weierstrass theorem any continuous function defined on a segment can be uniformly approximated by a polynomial function. Thus all forms of color indices can be approximated by a polynomial $y = \sum_{d=0}^{N} \alpha_{d} \rho_{d}^{\delta_{d}}$ of degree $N$. Setting the degree is a difficult task which may imply under-fitting or over-fitting. In addition un-stability can be caused by near-zero $\delta_{d}$. But since the segment is restricted to the domain $[0, 1]$ the Bernstein polynomials are a common demonstration and the equation can be wrote as a weighted sum of Bernstein basis polynomials $B_{N,i} = (1 - \rho)^{i} \rho^{N-i}$ which are more stable during the training. Moreover Bernstein Neural Network can solve partially differentiable equations [Sun et al., 2018]. For implementation reasons, two different layers are defined in the network (visible in the Figure
4.5). One for the Bernstein expansion limited to $B_{11,11}$ which takes the input image and produces different Bernstein basis polynomial, then each Bernstein basis is concatenated to the channel-wise and the linear combination is defined by a 2D convolution.

<figure><img src="image_109.png" /><figcaption>Figure 4.5 – Polynomial model with Bernstein expansions between $B_{4,1}$ and $B_{4,4}$.</figcaption></figure>

### 4.3.2.4 Universal Function Approximation

The Gaussian color space model proposed by [Geusebroek et al., 2000] shows that the spatio-spectral energy distribution of the incident light $E$ is the weighted integration of the spectrum $\rho_d$ denoted $E(\rho_d)$. Where $E$ can be described as a Taylor series and the energy function is convolved by different derivatives of a Gaussian kernel or structured receptive fields [Jacobsen et al., 2016]. This important point shows that Taylor expansions can decompose any function $f(x)$, especially for color decomposition and remapping, into:

$$f(x) = f(0) + f'(x)x + \frac{1}{2!}f''(x)x^2 + \frac{1}{3!}f'''(x)x^3 + o(x^3) \quad (4.8)$$

Here, the signature of the incident energies distribution of a remote sensing index associated to a surface can be reconstructed. An approach to learn this form of development is proposed by [Huang et al., 2016] which is commonly called DenseNet and then corresponds to the sum of the concatenation of the signal and these spatio-spectral derivatives

$$\mathbf{x} \rightarrow [\mathbf{x}, f_1(\mathbf{x}), f_2(\mathbf{x}), f_1(\mathbf{x})), \dots] \quad (4.9)$$

Various convolutions allow to learn receptive fields and derivatives in spectral domain when the kernel size $k$ is 1, and in spatio-spectral domain when $k$ is higher. Batch-Normalization are used to reduce the covariate shift across convolution output by re-scaling it and speed up the convergence. Finally the Sigmoid activation function is used and defined by

$$\text{Sigmoid}(x) = \frac{1}{1 + e^{-x}} \quad (4.10)$$

Sigmoid function allows to learn more complex structures and non-linearity of the reconstructed function. The number of derivative and receptive field are configurable with two parameters. The *depth* which corresponds to the number of layers in the network. And the *width* which refers to the number of outputs for each convolution. By default, the *depth* is fixed to 3 and the *width* is fixed to 5. The Figure 4.6 shows the corresponding universal function approximator network.
<figure><img src="image_110.png" /><figcaption>Figure 4.6 – Universal function approximation model (depth = 3, width = 5).</figcaption></figure>

### 4.3.2.5 Dense Morphological Function Approximation

As for the Taylor series, an approximation of any piecewise continuous function can be established by morphological operators such as dilatation and erosion [Mondal et al., 2019], respectively denoted $\rho \oplus s$ and $\rho \ominus s$ where s are the corresponding erosion or dilatation coefficients. Several erosion and dilation are defined for each spectral band i, then the expanded layer is defined as the channel concatenation of $z_i^+$ and in the same way for the erosion layer via $z_i^-$. Both are defined by

$$z_i^+ = \rho \oplus s_i = \max_k(\rho_k - s_{k,i}, 0) \qquad (4.11)$$

$$z_i^- = \rho \ominus s_i = \max_k(s_{k,i} - \rho_k, 0) \qquad (4.12)$$

To obtain the output $I = \sum_{i=0}^{N} z_i^+ w_i^+ + \sum_{i=0}^{N} z_i^- w_i^-$ of which the $w_i^+$ and the $w_i^-$ are the linear combination coefficients obtained by a 2D convolution. We chose to set the number of dilation and erosion neurons at 6. The Figure 4.7 shows the corresponding network.

<figure><img src="image_111.png" /><figcaption>Figure 4.7 – Dense-morphological model.</figcaption></figure>

### 4.3.3 Enhancing Baseline Models

#### 4.3.3.1 Input Band Filter (IBF)

To remove parts of the signal that may be dispensable, the addition of a low-pass, high-pass and band-pass filter upstream of the network are studied. A good example is provided by vegetation indices, only the high values in the green and near infra-red, and the low values in the red and blue characterize the vegetation.
This is the principle of the NDVI index. Due to the internal structure, the leaves reflect a lot of light in the near infrared, which is in sharp contrast to most non-vegetable surfaces. When the plant is dehydrated or stressed, the spongy layer collapse and the leaves reflect less light in the near-infrared, reaching red values in the visible range [Joshi et al., 2020]. Thus, the mathematical combination of these two signals can help to differentiate plants from non-plant objects and healthy plants from diseased plants. However, this index is then less interesting when detecting only vegetation and is strongly influenced by shade or heat.

We will therefore add a filter in the previous equations to remove undesirable spectral energies of each $ρ_d$ by using two thresholds a and b, which will also be learned. If it turns out that the whole signal is interesting, these two parameters will not change and their values will be a=0 and b=1. To apply the low-pass filter the equation $z = \max(\rho - a, 0) \div (1-a)$ is used and thus allows to suppress low values. For the high-pass filter the equation $w = \max(b - \rho, 0) \div b$ is applied to suppress high values. The band-pass filter it's the product of low and high-pass filters $y = z * x$. The output layer is the concatenation in the channel-wise of the input images, the low-pass, the high-pass and the band-pass filter which produce 4 × 13 = 52 channels. Finally to reduce the output data for the rest of the network, a bottleneck is inserted using a convolution layer, and generate a new image with 6 channels. This image is used by the rest of the network defined previously in Section 4.3.2. The Figure 4.8 shows the corresponding module inserted upstream of the network.

<figure><img src="image_112.png" /><figcaption>Figure 4.8 – Input Band Filter inserted at the beginning of the model.</figcaption></figure>

### 4.3.3.2 Spatial Pyramid Refinement Block (SPRB)

To take into account different scales in the image, the addition of a “Spatial Pyramid Refinement Block” at the downstream part of the network is studied. [Liu et al., 2015] showed that fusing the low to high-level features improved the segmentation task. It consists in the sum of different 2D convolutions whose core sizes have been set to 3, 5, 7 and 9. The results of all convolutions are concatenated and the final image output is given by a 2D convolution. The Figure 4.9 shows the corresponding module inserted downstream of the network.

<figure><img src="image_113.png" /><figcaption>Figure 4.9 – Spatial refinement block inserted at the end of a model.</figcaption></figure>
### 4.3.4 Last Activation Function

To obtain an index and facilitate convergence, we will only be interested in the values between 0 and 1 at the output of the last layer with the help of an activation function of type clipped ReLU defined by

$$ \text{ClippedReLU}(x) = \begin{cases} 1 & \text{if } x > 1 \\ x & \text{if } 0 < x < 1 \\ 0 & \text{if } x < 0 \end{cases} \qquad (4.13) $$

where x is a pixel of the output image. Each negative or null pixel will then be the unwanted class, greater or equal to 1 will be the searched class. The indecision border is the values between 0 and 1 which will be optimized. And then correspond to the probability that the pixel is the searched surface $P(Y = 1)$ or not $P(Y = 0)$. This is valid for the output prediction denoted $\hat{p} \in [0, 1]$ and the ground truth denoted $p \in \{0, 1\}$.

### 4.3.5 Loss Function

A wide variety of loss functions have been developed during the emergence of deep-learning (MSE, MAE, Hinge, Tversky, etc). A cross-entropy loss function is usually used when optimizing binary classification [Bokhovkin and Burnaev, 2019]. This loss function is not optimized for the shape. Recently, with deep neural network and for semantic segmentation [Rahman and Wang, 2016] has proposed a solution to optimize an approximation of the mean intersection over union (mIoU) and defined by

$$ \text{mIoU\_Loss} = 1 - \frac{p\hat{p}}{p + \hat{p} - p\hat{p}} \qquad (4.14) $$

The performance of this loss function seems more efficient than previous methods [Zhou et al., 2019, van Beers et al., 2019, Jadon, 2020]. We will then use it as a loss function.

### 4.3.6 Performance Evaluation

Commonly, accuracy and Pearson correlation are used to quantify the performance of remote sensing indices [Liu et al., 2017, Kokhan and Vostokov, 2020]. However this type of metrics does not take into account either the class ratio nor the shape of the segmentation. Correlation is also highly sensitive to non-linear relationship, noise, subgroups and outliers [Aggarwal and Ranganathan, 2016, Armstrong, 2019] making incorrect evaluation. According to [Shamir et al., 2019, Choi et al., 2019], the dice score and the mean intersection over union (mIoU) are more adapted to evaluate the segmentation mask. Defined by :

$$ \text{Dice} = \frac{2p\hat{p}}{p + \hat{p}} \qquad (4.15) \qquad \text{mIoU} = \frac{p\hat{p}}{p + \hat{p} - p\hat{p}} \qquad (4.16) $$

We will then use these two metrics for the performance evaluation. Prior to quantization, a threshold of 0.5 is applied to the output of the network to transform the probability into a segmentation mask. When $\hat{p}$ is lower than 0.5, it is considered as the background, otherwise it is considered as the object mask we are looking for. Other metrics are not considered because they are not always appropriate in case of segmentation or use in unbalanced data.
### 4.3.7 Comparison with Standard Indices

In order to make a fair comparison it is necessary to optimize each standard index. A minimal neural network is used to learn a linear regression. The network is thus composed of the spectral index, followed by a normalization $x = (x - min)/(min - max)$, then a 2D convolution with a kernel size of $k = 1$ is used for the linear regression. To perform the classification in the same way as our method, a ClippedReLU activation function is used. This tiny network is presented in the next Figure 4.10. Obviously the same metrics and loss function are used.

<figure><img src="image_114.png" /><figcaption>Figure 4.10 – Optimized model for standard indices</figcaption></figure>

### 4.3.8 Training Setup

The training is done through Keras module within Tensorflow 2.2.0 framework. All computation is done on an NVidia GTX 1080 which have 8111MiB of memory, this limits the number of simultaneous layers on the memory and so the size of the model. Each model is compiled with Adam optimizer. This optimization algorithm is primarily used with lookahead mechanism proposed by [Zhang et al., 2019a]. It iteratively updates two sets of weights : the search directions for the **fast weights** are chosen by inner optimizer, while the **slow weights** are updated every $k$ steps based on the direction of the **fast weights** and the two sets of weights are synchronized. This method improves the learning stability and lowers the covariance of its inner optimizer. The initial learning rate is fixed to $2^{-3}$. Batch size is fixed to 1 due to memory limitation. And the learning rate is decreased using *ReduceLROnPlateau* with factor = 0.2, patience = 5, min lr = 2e-6. The training is done through 300 iterations. Finally an EarlyStopping callback is used to stop the training when there is no improvement in the training loss after 50 consecutive epochs.
# 4.4 Results and Discussion

## 4.4.1 Fixed Models

All standard vegetation models have been optimized using the same training and validation datasets. Each of them has been optimized using a min-max normalization followed by a single 1 × 1 2D convolution layer and a last clipped ReLU activation function is used like the generic models implemented. The top nine standard indices are presented in Table 4.2. Their respective equations are available in Table 4.8 in Appendix 4.5.

<table><caption>**Table 4.2** – Synthesized standard indices performances : the nine best models are presented.</caption><thead><tr><th>Standard Index</th><th>Used ρ</th><th>mIoU</th><th>Dice</th></tr></thead><tbody><tr><td>Modified Triangular Vegetation Index 1</td><td>3</td><td>73.71</td><td>83.23</td></tr><tr><td>Modified Chlorophyll Absorption In Reflectance Index 1</td><td>3</td><td>73.68</td><td>83.22</td></tr><tr><td>Enhanced Vegetation Index 2</td><td>2</td><td>67.94</td><td>79.20</td></tr><tr><td>Soil Adjusted Vegetation Index</td><td>2</td><td>67.28</td><td>78.65</td></tr><tr><td>Soil And Atmospherically Resistant VI 3</td><td>2</td><td>65.86</td><td>77.61</td></tr><tr><td>Enhanced Vegetation Index 3</td><td>2</td><td>65.05</td><td>77.07</td></tr><tr><td>Global Environment Monitoring Index</td><td>2</td><td>65.04</td><td>77.01</td></tr><tr><td>Adjusted Transformed Soil Adjusted VI</td><td>3</td><td>64.96</td><td>77.00</td></tr><tr><td>NDVI</td><td>2</td><td>63.98</td><td>75.97</td></tr></tbody></table>

It is interesting to note that most of them are very similar to NDVI indices in their form. This shows that according to all previous studies, these forms based on a ratio of linear combination are the most stable against light variation. For example the following NDVI based indices are tested and show very different performances, highlighting the importance of weight optimization :

$$ 
\begin{aligned}
\text{NDVI} &= (\rho_5 - \rho_2) \div (\rho_5 + \rho_2) \\
\text{EnhancedVegetationIndex} &= 2.5 * (\rho_5 - \rho_2) \div (\rho_5 + 6 * \rho_2 - 7.5 * \rho_0 + 1) \\
\text{EnhancedVegetationIndex 2} &= 2.4 * (\rho_5 - \rho_2) \div (\rho_5 + \rho_2 + 1) \\
\text{EnhancedVegetationIndex 3} &= 2.5 * (\rho_5 - \rho_2) \div (\rho_5 + 2.4 * \rho_2 + 1) \\
\text{SoilAdjustedVegetationIndex} &= (\rho_5 - \rho_2) \div (\rho_5 + \rho_2 + 1) * 2 \\
\text{SoilAndAtmosphericallyResistant VI 3} &= 1.5 * (\rho_5 - \rho_2) \div (\rho_5 + \rho_2 + 0.5)
\end{aligned}
$$

The *Modified Triangular Vegetation Index 1* is given by $vi = 1.2*(1.2*(\rho_5-\rho_1)-2.5*(\rho_2-\rho_1))$ which shows that a simple linear combination can be as much efficient as NDVI like indices by taking one additional spectral band ($\rho_2 = \text{green}$) and more adapted coefficients. However, the other 80 spectral indices do not seem to be stable against of light variation and saturation. It is thus not relevant to present them.
4.4.2 Deepindices

Finally, each baseline model such as *linear*, *linear ratio*, *polynomial*, *universal function approximation* and *dense morphological function approximation* are evaluated with 4 different modalities of each kernel size *N* = 1, *N* = 3, *N* = 5 and *N* = 7. In addition *input band filter* (*ibf*) and *spatial pyramid refinement block* (*sprb*) are put respectively at the upstream and downstream of the network. Figure 4.11 shows that network synthesis. To deal with lighting variation and saturation a *BatchNormalization* is put in the upstream of the network in all cases. The *ibf* and *sprb* modules are optional and can be disabled.

<figure><img src="image_115.png" /><figcaption>**Figure 4.11** – Network synthesis with *ibf*, *evaluated index equation*, and *sprb*.</figcaption></figure>

When the input band filter (ibf) is enabled, the incoming tensor size of 1200 × 800 × 13 is transformed to a tensor of size 1200 × 800 × 6 and passed to the generic equation. When it is not, the generic equations get the raw input tensor of size 1200 × 800 × 13. In all cases the baseline model output a tensor of shape 1200 × 800 × 1. The spatial pyramid refinement block transforms the output tensor of the baseline model to a new tensor of the same size.

All models are evaluated with two metrics, respectively the dice and mIoU score. For each kernel size, the results are presented in Tables 4.3–4.6. All models are also evaluated with and without *ibf* and *sprb* for each kernel size.

<table><caption>**Table 4.3** – Score of DeepIndices with/without *ibf* and *sprb* for a kernel size of 1</caption><thead><tr><th rowspan="2">Model</th><th colspan="4">mIoU</th><th colspan="4">dice</th></tr><tr><th>baseline</th><th>ibf</th><th>sprb</th><th>ibf+sprb</th><th>baseline</th><th>ibf</th><th>sprb</th><th>ibf+sprb</th></tr></thead><tbody><tr><td>linear</td><td>78.58</td><td>79.63</td><td>78.88</td><td>78.12</td><td>87.56</td><td>88.34</td><td>87.57</td><td>86.93</td></tr><tr><td>linear-ratio</td><td>79.01</td><td>78.86</td><td>77.73</td><td>79.67</td><td>87.85</td><td>87.87</td><td>86.55</td><td>88.28</td></tr><tr><td>polynomial</td><td>70.08</td><td><b>80.03</b></td><td>74.47</td><td>79.32</td><td>80.53</td><td><b>88.61</b></td><td>84.07</td><td>88.03</td></tr><tr><td>universal-function</td><td>78.39</td><td>76.59</td><td>79.04</td><td><b>80.15</b></td><td>87.27</td><td>85.36</td><td>87.63</td><td><b>88.53</b></td></tr><tr><td>dense-morphological</td><td>76.15</td><td>78.86</td><td>75.96</td><td><b>80.00</b></td><td>85.26</td><td>87.80</td><td>85.15</td><td><b>88.54</b></td></tr><tr><td>Stats</td><td>—</td><td>2.35</td><td>0.78</td><td>3.01</td><td>—</td><td>1.90</td><td>0.50</td><td>2.37</td></tr></tbody></table>

<table><caption>**Table 4.4** – Score of DeepIndices with/without *ibf* and *sprb* for a kernel size of 3</caption><thead><tr><th rowspan="2">Model</th><th colspan="4">mIoU</th><th colspan="4">dice</th></tr><tr><th>baseline</th><th>ibf</th><th>sprb</th><th>ibf+sprb</th><th>baseline</th><th>ibf</th><th>sprb</th><th>ibf+sprb</th></tr></thead><tbody><tr><td>linear</td><td>78.89</td><td>78.21</td><td>78.53</td><td>79.76</td><td>87.66</td><td>87.16</td><td>87.35</td><td>88.36</td></tr><tr><td>linear-ratio</td><td>76.63</td><td>78.21</td><td>74.90</td><td>78.17</td><td>85.49</td><td>87.37</td><td>83.89</td><td>86.92</td></tr><tr><td>polynomial</td><td>72.83</td><td>79.31</td><td>73.20</td><td>79.13</td><td>83.06</td><td>88.13</td><td>82.78</td><td>87.82</td></tr><tr><td>universal-function</td><td>76.67</td><td>79.63</td><td>77.81</td><td><b>81.08</b></td><td>85.57</td><td>88.28</td><td>86.67</td><td><b>89.22</b></td></tr><tr><td>dense-morphological</td><td>76.54</td><td>79.39</td><td>75.65</td><td><b>80.29</b></td><td>85.43</td><td>88.17</td><td>84.40</td><td><b>88.66</b></td></tr><tr><td>Stats</td><td>—</td><td>2.64</td><td>-0.29</td><td>3.37</td><td>—</td><td>2.38</td><td>-0.42</td><td>2.75</td></tr></tbody></table>
<table><caption>**Table 4.5** – Score of DeepIndices with/without ibf and sprb for a kernel size of 5</caption><thead><tr><th rowspan="2">Model</th><th colspan="4">mIoU</th><th colspan="4">dice</th></tr><tr><th>baseline</th><th>ibf</th><th>sprb</th><th>ibf+sprb</th><th>baseline</th><th>ibf</th><th>sprb</th><th>ibf+sprb</th></tr></thead><tbody><tr><td>linear</td><td>77.80</td><td>78.83</td><td>78.92</td><td>79.92</td><td>86.91</td><td>87.67</td><td>87.61</td><td>88.24</td></tr><tr><td>linear-ratio</td><td>75.72</td><td>77.94</td><td>77.36</td><td><strong>80.08</strong></td><td>84.87</td><td>87.26</td><td>86.33</td><td><strong>88.43</strong></td></tr><tr><td>polynomial</td><td>73.11</td><td>79.92</td><td>73.69</td><td><strong>80.67</strong></td><td>83.29</td><td>88.58</td><td>83.31</td><td><strong>88.83</strong></td></tr><tr><td>universal-function</td><td>77.60</td><td><strong>80.63</strong></td><td><strong>80.31</strong></td><td><strong>80.63</strong></td><td>86.38</td><td><strong>89.02</strong></td><td><strong>88.53</strong></td><td><strong>88.71</strong></td></tr><tr><td>dense-morphological</td><td>74.89</td><td>79.74</td><td>76.04</td><td><strong>81.92</strong></td><td>83.84</td><td>88.42</td><td>85.09</td><td><strong>89.80</strong></td></tr><tr><td>Stats</td><td>–</td><td>3.59</td><td>1.44</td><td>4.82</td><td>–</td><td>3.13</td><td>1.12</td><td>3.74</td></tr></tbody></table>

<table><caption>**Table 4.6** – Score of DeepIndices with/without ibf and sprb for a kernel size of 7</caption><thead><tr><th rowspan="2">Model</th><th colspan="4">mIoU</th><th colspan="4">dice</th></tr><tr><th>baseline</th><th>ibf</th><th>sprb</th><th>ibf+sprb</th><th>baseline</th><th>ibf</th><th>sprb</th><th>ibf+sprb</th></tr></thead><tbody><tr><td>linear</td><td>79.08</td><td><strong>80.29</strong></td><td>79.25</td><td><strong>81.49</strong></td><td>87.75</td><td><strong>88.57</strong></td><td>87.80</td><td><strong>89.42</strong></td></tr><tr><td>linear-ratio</td><td>78.43</td><td><strong>80.58</strong></td><td>77.85</td><td><strong>81.35</strong></td><td>87.04</td><td><strong>88.78</strong></td><td>86.68</td><td><strong>89.45</strong></td></tr><tr><td>polynomial</td><td>72.49</td><td><strong>80.79</strong></td><td>74.14</td><td><strong>81.21</strong></td><td>82.92</td><td><strong>88.99</strong></td><td>83.77</td><td><strong>89.27</strong></td></tr><tr><td>universal-function</td><td>78.49</td><td><strong>80.20</strong></td><td><strong>80.21</strong></td><td><strong>80.36</strong></td><td>87.38</td><td><strong>88.72</strong></td><td><strong>88.35</strong></td><td><strong>88.70</strong></td></tr><tr><td>dense-morphological</td><td>75.70</td><td><strong>80.35</strong></td><td>76.34</td><td><strong>82.19</strong></td><td>84.48</td><td><strong>88.70</strong></td><td>85.61</td><td><strong>89.94</strong></td></tr><tr><td>Stats</td><td>–</td><td>3.60</td><td>0.72</td><td>4.48</td><td>–</td><td>2.84</td><td>0.53</td><td>3.44</td></tr></tbody></table>

For all baseline models, the results (in term of mIoU) show that increasing the kernel size also increases performances. The gain performance between best models in kernel size 1 and 7 are approximately 2% and then correspond to the influence of spectral mixing. So searching for spectral mixing 3 pixels farther (kernel size 7) still increases performance. It could also be possible that function approximation allows to spatially reconstruct some missing information.

For all kernel sizes, the **ibf** module enhances the mIoU score up to 3.6%. So the **ibf** greatly prunes the unneeded part of the input signal which increases the separability and the performances of all models. The **sprb** module allows to smooth the output by taking into account neighborhood indices, but their performance are not always better or generally negligible when it is used alone with the baseline model.

The baseline **polynomial** model is probably over-fitted, because it's hard to find the good polynomial order. But enabling the **ibf** fixes this issue. However further study should be done to setup the order of Bernstein expansion.

The **dense morphological** with a kernel size of 5 and 7 using both **ibf** and **sprb** modules is the best model in terms of dice (≈90%) and mIoU score (≈82%). Followed by **universal function approximator** with a kernel size of 1 or 3 with both **ibf** and **sprb** modules (dice up to 89% and mIoU up to 81%). Further studies on the width of the universal function approximator could probably increase performance. According to [Mondal et al., 2019] it seems normal that the potential of **dense morphological** is higher although the hyper-parameters optimization of **universal function approximator** could increase their performance.
### 4.4.3 Initial Image Processing

To show the importance of the initial image processing, each model has been trained without the various input transformations, such as $ρ_{std}$, Gxx, Gxy, Gyy filters, Laplacian filter, minimum and maximum Eigen values. Table 4.7 shows the score of DeepIndices considering only kernel size of 1 in different model.

<table><caption>**Table 4.7** – Score of DeepIndices in different modalities for a kernel size of 1 without initial image processing</caption><thead><tr><th rowspan="2">Model</th><th colspan="4">mIoU</th><th colspan="4">dice</th></tr><tr><th>baseline</th><th>ibf</th><th>sprb</th><th>ibf+sprb</th><th>baseline</th><th>ibf</th><th>sprb</th><th>ibf+sprb</th></tr></thead><tbody><tr><td>linear</td><td>72.34</td><td>74.29</td><td>72.94</td><td>76.97</td><td>83.15</td><td>84.66</td><td>83.03</td><td>86.50</td></tr><tr><td>linear-ratio</td><td>73.72</td><td>70.51</td><td>73.30</td><td>71.55</td><td>84.10</td><td>82.36</td><td>83.19</td><td>81.57</td></tr><tr><td>polynomial</td><td>74.33</td><td>74.14</td><td>77.88</td><td>76.42</td><td>85.07</td><td>84.49</td><td>87.19</td><td>85.94</td></tr><tr><td>universal-function</td><td>74.24</td><td>74.42</td><td>75.46</td><td>76.25</td><td>84.36</td><td>84.49</td><td>85.16</td><td>85.86</td></tr><tr><td>dense-morphological</td><td>72.04</td><td>73.72</td><td>71.03</td><td>74.69</td><td>82.27</td><td>84.00</td><td>81.33</td><td>84.72</td></tr><tr><td>Stats</td><td>—</td><td>0.08</td><td>0.79</td><td>1.84</td><td>—</td><td>0.21</td><td>0.19</td><td>1.13</td></tr></tbody></table>

The results shows that none of optimized models outperforms the previous performance with the initial image processing (best mIoU at 80.15%). The maximum benefit is approximately 6% for mIoU score depending on the model and module, especially when using combination of **ibf**, **sprb** and small kernel size. Meaning that signal processing is much more important than spectral mixing and texture.

### 4.4.4 Discussion

Further improvements can be set on hyper-parameters of the previously defined equations, such as the degree of the polynomial (set to 11), the CNN depth and width for Taylor series (set to 3) and the number of operations in morphological network (set to 10). In particular the learning of 2D convolution kernel of Taylor series may be replaced by a structured receptive field [Jacobsen et al., 2016]. In addition it would be interesting to transpose our study with new data for other surfaces such as shadows, waters, clouds or snows.

The training dataset is randomly split with a fixed seed, which is used for every learned models. As previously noted, this is important to ensure reproducible results but could also favor specific models. Further work to evaluate the impact of varying training datasets could be conducted.
### 4.4.4.1 Model Convergence

Another way to estimate the robustness of a model against its initialization is to compare the model's convergence speed. Models with faster convergence should be less sensitive to the training dataset. As an example, the convergence speed of few different models is shown in Figure 4.12. The baseline model convergence is the same, as well as sprb module. However the speed of convergence also increases with the size of the kernel but does not alter subsequent observations. For greater readability only models with ibf are presented.

<figure><img src="image_116.png" /><figcaption>**Figure 4.12** – first 80 epochs of loss of generic models with ibf in kernel size of 1.</figcaption></figure>

An important difference in the speed of convergence between models is observed. An analysis of this figure allows the aggregation of model types and speed :

- Slow converging models : polynomials models converge slowly as well as the majority of linear or linear-ratio models.

- Fast converging models : universal-functions and dense-morphological are the fastest to converge (less than 30 iterations)

A subset of slow and fast converging models could be evaluated in term of sensitivity against initialization. It shows that the dense morphological followed by universal function approximator convergence faster than the other. Regardless of the used module nor kernel size.
4.4.4.2 Limits of Deepindices

Shadows can be a relatively hard problem to solve in image processing, the proposed models
are able to correctly separate vegetation from soil even with shadowy images, as shown in Figure
4.13. In addition, the Figure 4.16 in the Appendix 4.5 shows the impacts of various acquisition
factors, such as shadow, noise, specular or thin vegetation features.

<figure><img src="image_117.png" /><figcaption>**Figure 4.13** – Correct vegetation/soil discrimination despite shadows.</figcaption></figure>

Some problem occurs when there are abrupt transitions between shadowed and light areas
of an image as shown in Figure 4.14.

<figure><img src="image_118.png" /><figcaption>**Figure 4.14** – Vegetation/soil discrimination issue with abrupt transition between shadow and light.</figcaption></figure>

It appears that the discrimination error appears where the shadow is cast by a solid object,
resulting in edge diffraction that creates small fringes on the soil and vegetation. A lack of such
images in the training dataset could explain the model failure. Data augmentation could be
used to obtain a training model containing such images, from cloud shadows to solid objects
shadows. Further work is needed to estimate the benefit of such a data augmentation on the
developed models.

The smallest parts of the vegetation (less than 1 pixel, such as small monocotyledon leaves or
plant stems) cannot be detected because of a strong spectral mixture. This limitation is due to
the acquisition conditions (optics, CCD resolution and elevation) and should be considered as is.
As vegetation with a width over 1 pixel is correctly segmented by our approach, the acquisition
parameters should be chosen so that the smallest parts of vegetation that are required by an
application are larger than 1 pixel in the resulting image.
A few spots of specular light can also be observed on images, particularly on leaves. These spots are often unclassified (or classified as soil). This modifies the shape of the leaves by creating holes inside them. This problem can be seen on Figure 4.15. Leaves with holes are visible on the left and the middle of the top bean row. It would be interesting to train the network to detect and assign them to a dedicated class.

<figure><img src="image_119.png" /><figcaption>**Figure 4.15** – Vegetation/soil discrimination issue caused by specular lights on leaves.</figcaption></figure>

Next the location of the detected spots could be studied to re-assign them to two classes : specular-soil and specular-vegetation. To perform this step, a semantic segmentation could be set up to identify the surrounding objects of the holes specifically. It would be based on the UNet model, which performs a multi-scale approach by calculating, treating and re-convolving images of lower resolutions.

More generally, the quality of the segmentation between soil and vegetation strongly influences the discrimination between crop and weed, which remains a major application following this segmentation task. Three categories of troubles have been identified : the plants size, the ambient light variations (shades, specular light spots), and the morphological complexity of the studied objects.

The size of the plants mainly impacts their visibility on the acquired images. It is not obviously related to the ability of the algorithm to classify them. However, it leads to the absence of essential elements such as monocotyledon weeds at an early vegetation stage. A solution is proposed by setting the acquisition conditions to let the smallest vegetation part be over 1 pixel.

Conversely, the variations of ambient light should be treated by the classification algorithm. As previously mentioned, shadow management needs an improvement of the learning base, and specular light spots could be treated by a multi-scale approach. Their influence on the discrimination step should be major. Indeed, they influence the shape of the objects classified as plants, which is a useful criterion to discriminate crops from weeds. The morphological complexity of the plants can be illustrated by the presence of stems. In our case, bean stems are similar to weed leaves. This problem should be treated by the discrimination step. The creation of a stem class (in addition to the weed and crop classes) will be studied in particular.
**4.5 Conclusions**

In this work, different standard vegetation indices have been evaluated as well as different methods to estimate new DeepIndices through different types of equations that can reconstruct the others. Among the 89 standard vegetation indices tested, the MTVI (Modified Triangular Vegetation Index 1) gives the best vegetation segmentation. Standard indices remain sub-optimal even if they are downstream optimized with a linear regression because they are usually used on calibrated reflectance data. The results allow us to conclude that any simple linear combination is just more efficient (+4.87% mIoU) than any standard indices by taking into account all spectral bands and few transformations. The results also suggest that un-calibrated data can be used in proximal sensing applications for both standard indices and DeepIndices with good performances.

We therefore agree that it is important to optimize both the arithmetic structure of the equation and the coefficients of the spectral bands, that is why our automatically generated indices are much more accurate. The best model is much more efficient by +8.48% compared to the best standard indices and by +18.21% compared to NDVI. Also the two modules **ibf**, **sprb** and the initial image transformation show a significant improvement. The developed DeepIndices allow to take into account the lighting variation within the equation. It makes possible to abstract from a difficult problem which is the data calibration. Thus, partially shaded images are correctly evaluated, which is not possible with standard indices since they use spectrum measurement that change with shades. However, it would be interesting to evaluate the performance of standard indices and DeepIndices on calibrated reflectance data.

These results suggest that deep learning algorithms are a useful tool to discover the spectral band combinations that identify the vegetation in multi spectral camera. Another conclusion from this research is about the genericity of the methodology developed. This study presents a first experiment employed in field images with the objective of finding deep vegetation indices and demonstrates their effectiveness compared to standard vegetation indices. This paper's contribution improves the classical methods of vegetation index development and allows the generation of more precise indices (i.e., DeepIndices). The same kind of conclusion may arise from this methodology applied on remote sensing indices to discriminate other surfaces (roads, water, snow, shadows, etc).
**Author Contributions** Conceptualization, Vayssade, J.-A.; data curation, Vayssade, J.-A.; formal analysis, Vayssade, J.-A.; funding acquisition, Jones, G and Paoli, J.-N.; investigation, Vayssade, J.-A.; methodology, Vayssade, J.-A.; project administration, Paoli.J-N. and Jones, G; resources, Paoli.J-N. and Jones, G; software, Vayssade, J.-A.; supervision, Paoli, J.-N. and Jones, G; validation, Paoli, J.-N.; Gée, C. and Jones, G.; visualization, Vayssade, J.-A.; writing—original draft preparation, Vayssade, J.-A.; writing—review and editing, Vayssade, J.-A.; Paoli, J.-N.; Gée, C. and Jones, G.; All authors have read and agreed to the published version of the manuscript.

**Funding** This project is funded by ANR Challenge RoSE and the Horizon 2020 project IWM-PRAISE.

**Data Availability** Data in this study is publicly available at https://data.inrae.fr/dataset.xhtml?persistentId=doi:10.15454/DSQC8N, using Creative Common CC0 1.0 Public Domain Dedication licence

**Acknowledgments** We would like to thank Masson Jean-Benoit for the realization of the metal gantry which allowed us to position the camera at different heights, it was used in particular for the calibration of the camera and the band registration. We also thank Djemai Mehdi for the spelling correction of the English. And we thank Aubry Clément and Cozic Thibault of the company SITIA for their help in interfacing the camera with the used robot “Trecktor”.

**Conflicts Of Interest** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendices and supplementary material

$$b = \rho_0, g = \rho_1, r = \rho_2, e = \rho_3, u = \rho_4, n = \rho_5$$

<table><caption>**Table 4.8** – Top optimized fixed vegetation model equations</caption><thead><tr><th>Model</th><th>Equation</th></tr></thead><tbody><tr><td>Modified Triangular Vegetation Index 1</td><td>1.2 * (1.2 * (n - g) - 2.5 * (r - g))</td></tr><tr><td>Modified Chlorophyll Absorption In Reflectance Index 1</td><td>1.2 * (2.5 * (n - r) - 1.3 * (n - g))</td></tr><tr><td>Enhanced Vegetation Index 2</td><td>2.4 * (n - r)/(n + r + 1)</td></tr><tr><td>Soil Adjusted Vegetation Index</td><td>1.5 * (n - r, n + r + 0.5)</td></tr><tr><td>Soil And Atmospherically Resistant VI 3</td><td>1.5 * (n - r, n + r)</td></tr><tr><td>Enhanced Vegetation Index 3</td><td>2.5 * (n - r)/(n + 2.4 * r + 1)</td></tr><tr><td>Global Environment Monitoring Index</td><td>(2*(n²-r²)+1.5*n+0.5*r) * (1 - n/4) - (r-0.125)/(1+r)</td></tr><tr><td>Adjusted Transformed Soil Adjusted VI</td><td>a * ((n-r)/(n+r)) * ((n-r)/(n+r))</td></tr><tr><td>NDVI</td><td>(n-r)/(n+r)</td></tr></tbody></table>
<figure><img src="image_120.png" /><figcaption>**Figure 4.16** – Visual comparison between some relevant models. NDVI (63.98 mIoU), MTVI1 (73.71 mIoU), linear 1 baseline (78.58 mIoU), dense 7 ibf-sprb (82.19 mIoU). Blue indicates sure soil, red indicates sure vegetation, and the other colors indicate uncertainty.</figcaption></figure>
## 4.6 Conclusion de chapitre

Dans cet article, de nouveaux indices de végétation (Deep Indices) ont été proposés. En partant de la définition d'une large quantité d'équations d'indices de végétation (89), il a été possible de déterminer deux équations simples qui généralisent la majorité d'entre elles, tout en considérant l'ensemble des bandes spectrales. Pour aller plus loin, des méthodes reposants sur des approximateurs de fonctions (développements de Taylor, polynômes de Bernstein, et opérateurs morphologiques) ont également été testées, de même que l'ajout de transformations des données d'entrée (déviation standard, matrice hessienne, ...) ou de traitement du signal en amont des fonctions à optimiser (filtre passe-bas, passe-haut et passe-bande).

Ces indices ont été développés au travers d'un outil d'apprentissage profond permettant des temps de traitement minimes. De plus, une fonction de perte adaptée permet d'optimiser ces indices pour obtenir la segmentation la plus proche de celle souhaitée. Un ensemble conséquent de 89 indices standards a également été optimisé en suivant cette démarche. Cela permet de comparer indices standards et indices proposés avec des degrés d'optimisations similaires.

Ces indices ont été appris sur un ensemble de 91 images non calibrées et contenant des variations importantes des conditions d'acquisition. La taille du jeu d'apprentissage ne nécessite pas d'être plus conséquente, en revanche le choix des images et la représentation des différentes conditions d'acquisition dans des proportions spécifiques jouent un rôle important. Pour ce type de démarche, il est nécessaire de sous-représenter les images simples, acquises dans de bonnes conditions. L'utilisation de techniques d'augmentation de données (bruit de Perlin) renforce également la robustesse des indices appris.

Les résultats montrent que certains indices standards peuvent proposer des résultats intéressants en termes de segmentation sol/végétation après optimisation. Le MTVI (Modified Triangular Vegetation Index 1) donne ainsi les meilleurs résultats (73.71 mIoU). A contrario, le NDVI, l'indice le plus utilisé dans la littérature, montre des résultats moyens (63.98 mIoU) pour cette tâche. Ces indices standards, sont facilement dépassés par une simple régression linéaire (78.58 mIoU). Et ils sont largement moins performants que les indices proposés dans notre étude (82.19 mIoU). Certains artefacts ne sont toutefois pas traités par les Deep-Indices. C'est le cas de la réflexion spéculaire (figure 4.15), des aberrations chromatiques visibles en bordure d'ombre (figure 4.14) et des pixels présentant un mélange spectral (figure 4.16), souvent en bordure de feuille, ou sur les feuilles fines.

Grâce à ces nouveaux indices, ne nécessitant pas de calibration, nous avons répondu efficacement à la première question de l'étape de segmentation : comment discriminer la végétation du reste de la scène ? Ceci constitue un des apports scientifiques de cette thèse. Maintenant que la première étape de segmentation est achevée, la seconde question est : comment différencier les individus au sein d'un couvert végétal ? Pour cela, plusieurs échelles sont possibles : composantes connexes, plantes, pixels, .... Les méthodes utilisables dépendent de ce choix et le chapitre suivant, en jaune sur la figure ci-dessous, propose une réponse à cette problématique.

<figure><img src="image_121.png" /></figure>

Rappel de l'architecture de notre chaîne fonctionnelle
ici les premières étapes en vert sont terminées, et en jaune l'étape suivante étudiée
# 5 Détection des individus

Après avoir séparé la végétation du reste de l'image, l'étape suivante vise à isoler des individus au sein de la végétation. Les algorithmes pouvant être utilisés pour cette opération ont été introduits dans la section 2.3. Toutefois, le choix de l'algorithme adéquat dépend du type d'entité recherchée. Il est en effet possible de s'intéresser à une composante connexe, une plante, une feuille, un super-pixel ou encore à un simple pixel.

Dans le cadre du challenge RoSE, la densité des plantes est importante. Il existe souvent un recouvrement entre les plantes (figure 5.1). Si l'on définit l'entité comme étant une composante connexe, alors cette entité correspondra parfois à un amas de plantes (quand il y a recouvrement) et parfois à une plante isolée. Cela induira des erreurs d'étiquetage pour les amas de plantes, donc une classification grossière et l'impossibilité de détecter les adventices intra-rang, ce qui est l'objectif initial du challenge. A contrario, si nous définissons une segmentation au niveau du pixel, alors nous manquerons de propriétés pour une classification cultures/adventices performante.

Les problématiques de segmentation, de caractérisation et de classification sont fortement liées. En réalité, la définition d'une méthode permettant d'obtenir un niveau de segmentation adéquat pour la classification culture/adventices constitue l'un des verrous scientifiques majeurs de cette thèse.

<figure><img src="image_122.png" /><figcaption>**Figure 5.1** – Mélange de haricot (culture) et de moutarde (adventice) dans le challenge RoSE</figcaption></figure>

Ainsi, dans notre contexte d'utilisation, en proxidétection, l'échelle de la composante connexe et celle du pixel sont à rejeter. L'échelle du superpixel pourrait être intéressante, en revanche, les algorithmes permettant une telle segmentation, sont longs et peu fiables, car fortement influencés par un paramètre essentiel, celui de la taille des superpixels. En effet, trop petits, les superpixels sont trop nombreux et ne contiennent pas suffisamment d'information discriminante. Trop grands, ils sont alors constitués de mélanges des classes culture et adventice.

En agriculture de précision, pour la détection des adventices, l'entité de la plante est fréquemment utilisée. En effet, d'un point de vue agronomique, la plante est l'entité la plus fine utilisée pour caractériser, entre autres, l'état sanitaire et la pression des adventices.

Les approches historiques de segmentation ne permettent pas de segmenter une entité plante de manière stable en conditions réelles car les images présentent de fortes variations. De plus, les tailles et formes des plantes varient fortement entre individus au sein d'une même espèce, ce qui rend tout processus d'optimisation préalable impossible sur des algorithmes tels que les superpixels. C'est pourquoi, aujourd'hui, l'essentiel de la communauté a recours à des algorithmes reposant sur de l'apprentissage profond.
Les algorithmes d'apprentissage profond pris en considération en agriculture de précision utilisent des réseaux de neurones convolutifs (CNN) associés à une proposition des régions d'intérêts. Ces algorithmes de proposition de régions examinent une image d'entrée à travers un ensemble de convolutions qui maximisent l'identification de l'emplacement potentiel d'un objet. Cette méthode repose sur une régression de "bounding box", telle que proposée par [Girshick et al., 2014] (RCNN) ou [Kim and Cho, 2020] (YOLO).

<figure><img src="image_123.png" /><figcaption>**Figure 5.2 –** Structure d'un réseau “You Only Look Once” (YOLO)</figcaption></figure>

Les premiers travaux de [Girshick et al., 2014] utilisaient ensuite ces régions d'intérêt pour extraire les caractéristiques de sortie du réseau de neurones, afin de les introduire dans une machine à vecteurs de support (SVM) pour la classification finale. Au fil des évolutions, la classification SVM a été remplacée par un réseau de neurones dédié à l'identification de la classe de l'objet [Girshick, 2015], permettant un apprentissage dit de "bout en bout". Cette régression des régions d'intérêt a toutefois deux inconvénients : (1) elle propose une grande quantité de faux positifs et (2) des sorties multiples, mais légèrement décalées, pour un même objet. Pour pallier ces problèmes, le "score detection" est utilisé pour supprimer la majorité des faux positifs, ainsi qu'un algorithme de suppression des non-maximum (NMS). Le NMS permet d'éviter la détection répétée de la même instance en combinant les régions (bounding box) qui se chevauchent en une seule boîte englobante pour chaque objet.

Hélas, ce type d'approche montre ses limites en environnement dense, du fait de l'algorithme NMS. Ainsi, détecter des plantes par cette méthode ne permet pas un résultat optimal dans le cadre du challenge RoSE. En réalité, de par les nombreuses occultations et recouvrements inter-individus, l'échelle de la plante est discutable. De ce constat, il a été décidé de travailler sur une entité plus petite : la feuille. De plus, nous proposons ici, le développement d'un réseau de neurones ne s'appuyant pas sur l'algorithme NMS pour contourner ces limitations.

En allant plus loin sur les réseaux de neurones définis pour les indices de végétation et la segmentation sol/végétation, il est possible de proposer un réseau de segmentation sémantique pour la segmentation des instances des feuilles. L'approche repose sur la détection des contours des feuilles [Morris, 2018]. Cette approche est peu utilisée mais elle peut être améliorée de différentes façons. La suite de ce chapitre porte sur les travaux qui ont été menés dans ce domaine pour développer une méthode de segmentation de feuille en couvert dense (chevauchement d'objets) et est l'objet d'une publication internationale dans la revue *Computers and Electronics in Agriculture* (travail en cours de révision).
## Pixelwise Instance Segmentation Of Leaves In Dense Foliage

<figure><img src="image_124.png" /></figure>

**Authors** Jehan-Antoine Vayssade, Gawain Jones, Christelle Gée and Jean-Noël Paoli

**Publication** Sent on 20/08/2021 – Computer and Enlectronics in Agriculture – under-review

**Abstract** Detecting and identifying plants using image analysis is a key step for many applications in precision agriculture (from phenotyping to site specific weed management). Instance segmentation is usually carried on to detect entire plants. However, the shape of the detected objects changes between individuals and growth stages. A relevant approach to reduce these variations is to narrow the detection on the leaf. Nevertheless, segmenting leaves is a difficult task, when images contain mixes of plant species, and when individuals overlap, particularly in an uncontrolled outdoor environment. To leverage this issue, this study based on recent Convolutional Neural Network mechanisms, proposes a pixelwise instance segmentation to detect leaves in dense foliage environment. It combines “deep contour aware” (to separate the inner of big leaves from its edges), “Leaf Segmentation trough classification of edges” (to separate instances with a specific inner edges) and “Pyramid CNN for Dense Leaves” (to consider edges at different scales). But the segmentation output is also refined using a Watershed and a method to compute optimized vegetation indices (DeepIndices). The method is compared to others running the leaf segmentation challenge (provided by the International Network on Plant Phenotyping) and applied on an external dataset of Komatsuna plants. In addition, a new multispectral dataset of 300 images of bean plants is introduced (with dense foliage, individuals overlapping, mixes of species and natural lighting conditions). The ground truth (e.g. the leaves boundaries) is defined by labelled polygons and can be used to train and assess the performance of various algorithms dedicated to leaf detection or crop/weed classification. On the usual datasets, the performances of the proposed method are similar to those of the usual methods involved in the leaf segmentation challenges. On the new dataset, their results are strongly better than those of the usual RCNN method. Remaining errors are bad fusion between neighboring areas and over segmentation of multi-foliate leaves. Structural analysis methods could be studied in order to overcome these deficiencies.

**Keywords** precision agriculture; remote sensing; leaf segmentation; dense foliage; boundary detection; semantic segmentation; CNN; multispectral
### 5.1 Introduction

In precision agriculture, one of the hardest tasks is the detection of crops and weeds by imagery. Imaging systems will play a significant role in the new generation of agriculture, from the genetic selection in phenotyping [Omari et al., 2020] to site-specific weed management [Louargant et al., 2017]. They are also used to record frequently and accurately the plant growth, crop yield, leaf area, etc [Gée et al., 2021]. These data are then used to quantify and evaluate the quality of production.

In proximal detection, this work is mainly done at the plant level and gives important agronomic information once reported at the plot level (number of weeds, crops and weeds locations, diseases, stress ...). The plant level is also needed for the new agricultural revolution, such as robotic weed management. Thus, studies try to detect crop and weeds plants by using a wide variety of techniques, which have been reviewed by [Wang et al., 2019a]. The instance segmentation is a key-step used before the task of classifying plants as crop and weeds. In CNN field, it is mostly based on a major class approaches [Hafiz and Bhat, 2020], such as Mask-RCNN. However, detecting the entire plants has 3 main limitations (i) when plants are too numerous, it is hard to distinguish individuals when they overlap, one of them is undetected or both are merged which is due to the underlying non-max suppression algorithm [Bonneau et al., 2020], (ii) small and thin elements are undetected and (iii) the number of variations for a plant is infinite : the number of leaves, their orientations, their sizes and other differences that radically change between individuals and growth stages, implies a large amount of data for the training process. These limitations logically impact the detection rate and may cause a significant quantity of miss-detection. To solve this problem, one approach that should be relevant is to base the detection not on the whole plant, but on the leaf. For this purpose, a pixel-wise instance segmentation is proposed. The idea is to separate the instances by detecting the edges of the leaves whose projected shadow gives an interesting gradient break, particularly in an uncontrolled outdoor environment.

#### 5.1.1 Related Works

Some recent works on leaf segmentation were found in controlled illumination environments with limited occlusion between individuals. Especially on an open dataset of Arabidopsis Thaliana [Scharr et al., 2016]. Other studies related to biomedical imagery and nuclei segmentation were also found with pixel-wise instance segmentation. These studies show the importance of defining one or more edges classes and a relevant loss function. These two factors are both related to the weight given to each sample of the training dataset in the estimated error for the optimization algorithm. But the way of how parts of the network are dedicated to the edge classification and thus on how the network focuses on instance separation. Edge detection therefore plays a significant role in pixelwise instance segmentation task.
The first studies were dedicated to the definition of specific class of edges to separate instances. Thus [Chen et al., 2016] proposed a novel approach named “Deep Contour-Aware” based on the semantic segmentation of two classes, one class is dedicated to the inside of glands, while the second is for the segmentation of glands edges. A bit later, [Morris, 2018] was the first to define a pixelwise instance segmentation for dense leaf detection. They proposed the “Pyramid CNN for Dense Leaves” architecture which is similar to U-Net (the most common CNN used for biomedical images segmentation). The network is dedicated to the detection of leaf boundaries. To facilitate the learning of edges at different scales, an auxiliary loss function is placed at each sub-scale of the pyramid. Finally, the instances are retrieved by using a superpixel algorithm.

[Cui et al., 2019] enhance the “Deep Contour-Aware” model proposed in 2016 by using a real U-Net architecture and data augmentation. Concerning edge classes based instances segmentation, [Bell and Dee, 2019] study shows the importance of separating edges into two classes. As the outer edges are dominant in the samples, the corresponding error on contiguous edges are less important. Thus, the edges of leaves are separated into two classes, one for the outer edges and one for inner edges (when leaves are touching or overlapping). The multi-scale loss function is still used and proposed by [Xie et al., 2020] for nucleus instance segmentation. They show that multi-scale loss helps to regularize the network and narrow down the perceptual distances and enlarge the semantic dissimilarity between individuals. In addition, a count ranking loss is used on the last feature layer of a fully-connected layer. This count ranking loss enforces the network to focus on the learning of samples containing crowded nuclei. This technique results in an implicitly regularized trained network, to be aware of individuals quantity.

### 5.1.2 Objectives

All these studies show that pixelwise instance segmentation technique is viable but still limited. They also underline the importance of choosing a loss function adapted to the defined semantic classes. Based on these related works, this paper proposes to merge most recent advances in the field of pixelwise instance segmentation.

First, the proposed method combines “deep contour aware” (to separate the inner of big leaves from its edges), “Leaf Segmentation through classification of edges” (to separate instances with a specific inner edges) and “Pyramid CNN for Dense Leaves” (to consider edges at different scales). Second, a new loss function is also introduced to limit under and over-segmentation. Third output is refined using a specific vegetation index based on previous work [Vayssade et al., 2021] and a watershed algorithm.

This method is applied to dense leaf segmentation, on images containing mixes of plant species and acquired in natural light. A specific multispectral dataset has been acquired, it is presented in the next section 5.2.1 and released publicly. The method is also evaluated on two common online leaf RGB image databases [Scharr et al., 2016], presented in section 5.2.2.

The proposed method is much more robust compared to previous pixel-wise instance segmentation methods and solves the issue of methods based on bounding-box regression : all small and thin elements are detected.
## 5.2 Material and data

### 5.2.1 Specific multispectral dataset

#### 5.2.1.1 Experimental plot

The data were acquired at the site of INRAE (Figure 5.3) in Montoldre (Allier, France, at 46°20′30.3″N 3°26′03.6″E) within the framework of the “RoSE challenge” founded by the French National Research Agency (ANR) in 2019. The aim of the Challenge is to objectively compare the solutions proposed by participants for intra-row weed control ([Avrin et al., 2020]). Within this context, the challenge provides to contestants an evaluation plan and a set of experimental plots of bean and maize plants. In addition various natural weeds (yarrows, amaranth, geranium, plantago, etc) and sown ones (mustards, goosefoots, mayweed and ryegrass) are managed to compare performances.

<figure><img src="image_125.png" /><figcaption>**Figure 5.3** – Aerial view of the experimentation plot located in Montoldre (now INRAE)</figcaption></figure>

#### 5.2.1.2 multispectral camera

The images were acquired with the Airphen (Hyphen, Avignon, France) six-band multispectral camera (visible on the upper-left of the Figure 5.4). This is a scientific camera developed by agronomists for agricultural applications. The camera embeds six sensors using six bandpass (450/570/675/710/730/850 nm) filters with a 10 nm FWHM each. The focal length of each lens is 8 mm. The raw resolutions for each spectral band are 1280 × 960 px with 12 bit precision. Due to the conception of the camera, spectral images are not aligned. Based on previous work [Vayssade et al., 2020], a method for registration has been developed with a registration accuracy down to sub-pixel. After the registration, all spectral images are cropped to 1200*800 px and concatenated to channel-wise where each dimension refers to a spectral band.

#### 5.2.1.3 Image acquisition and annotation

From the presented experimental plots, a set of images were acquired. The camera is attached in front of an hybrid autonomous tractor called “TREKTOR” launched by SITIA Company (Bouguenais, France) in 2019. The camera is setup to have a top-down view of crop rows, thus it is mounted on a pole in front of the platform allowing to remove visible part of the robot and at 1.8 m from the ground. The Figure 5.4 below shows the arrangement of the elements. Crops and weeds were between phenological state 3 and 4 which means they have between 2 and 6 leaves. The ground truth is defined on images by experts with polygons around each leaf boundary. In addition, polygons contain a label for their classification between crop and weed (not used in this study). The annotation was defined using the VIA annotation software [Dutta and Zisserman, 2019] and a total of 300 images of bean were annotated, 170 from acquisition made in June and 130 in October. These dataset is freely available at https://doi.org/10.15454/JMKP9S.
<figure><img src="image_126.png" /><figcaption>Figure 5.4 – The experimental set-up : the multispectral camera and the robotic plateform</figcaption></figure>

**5.2.2 Online image databases**

**5.2.2.1 Additional data from the Computer Vision Problems in Plant Phenotyping dataset (CVPPP)**

To compare the method proposed in this study with others, an additional dataset was
used. This dataset is proposed for the Leaf Segmentation Challenge (LSC [Scharr et al., 2017])
provided by the International Network on Plant Phenotyping (INPP), a very popular challenge
for data scientists. It is composed of RGB images of Arabidopsis Thaliana (783 images) and
Rosette (27 images) plants segmented into several leaves. The authors state that the images
were collected from multiple locations in a growth chamber experiment and divided into four
groups, named A1 through A4. In addition, the dataset is composed of various image sizes
(respectively 530 × 530, 565 × 565, 2048 × 2048, 441 × 441 for each sub-dataset), which have
been resize to 512 × 512 × 3.

<figure><img src="image_127.png" /><figcaption>**Figure 5.5** – Example of images from the CVPPP dataset (top) and their corresponding ground truth (bottom)</figcaption></figure>
### 5.2.2.2 Additional data from Komatsuna dataset

Similar to the LSC dataset, Komatsuna dataset consists in RGB images of plants taken in a top view. This dataset contains a large number of plant growth stages, it was designed to solve the problem of 3D phenotyping [Uchiyama et al., 2017]. It is used here primarily as a test set due to its similarity to the CVPPP dataset (training set), as shown in Figure 5.6. This dataset includes 900 images of size 480 × 480 × 3.

<figure><img src="image_128.png" /><figcaption>**Figure 5.6** – Example of RGB images from the Komatsuna dataset (top) and their corresponding ground truth (bottom)</figcaption></figure>

## 5.3 Methodology

We consider the leaf segmentation as a binary segmentation problem of boundaries as proposed by [Morris, 2018]. The main idea is to detect the sharp edges of leaves or to follow the projected shadow of a leaf on the one below it. This methodology section is split into three sub-tasks (5.3.1) the proposed CNN architecture to detect and separate leaves, (5.3.2) the specific loss function defined to limit under and over-segmentation, and (5.3.3) a simple watershed algorithm which takes the CNN output and a vegetation mask to refine the segmentation.

### 5.3.1 Proposed CNN architecture

Unlike recent CNN architectures, the proposed approach is slightly more decomposed like standard biomedical and agricultural computer vision pipelines [Perez-Sanz et al., 2017, Lottes and Stachniss, 2017]. Thus, the architecture (Figure 5.7) is composed of three **upstream modules** (IIT, IBF, UFA) that improves the input data and eliminates unnecessary information. This step composed of 3 **upstream** modules, was proposed in a previous work to construct a vegetation index ([Vayssade et al., 2021]). It is used to identify relevant spectral features on the input data to exploit the inter-channel relationships. After this stage, a **core network** is used to consider spatial information at different scales on the image, the core network returns four down-scaled feature maps. Finally, at the end of the network, three **downstream modules** (CoordConv, UFA, Classification) are proposed to fuse spectral and spatial information. Sigmoid activation function is used at the end of all convolution layers to learn more complex structures and allows non-linearity of the reconstructed function.
<figure><img src="image_129.png" /><figcaption>**Figure 5.7** – Diagram of the CNN network architecture and losses (dotted). Multiple arrow show concatenation as input layer.</figcaption></figure>

The network takes an input image of size 512 × 512 × 6, thus the learning and computation are done by a sliding window within the registered images of the Airphen camera of size 1200 × 800 × 6. The output layer is defined as a semantic segmentation of four classes. One class is dedicated to the inner of individuals, while three classes are dedicated to the detection of edges to keep aware of leaves instance. As mentioned by [Bell and Dee, 2019], one class is dedicated to outer of leaf (touching soil texture), while the second is dedicated to inner edges (touching/overlapping another leaf). Within our dataset, a third class appears with thin leaves which can be considered as a kind of edge. The ground truth is a set of polygons, drawn with the respect of these classes. Edge classes were empirically drawn with three pixel thickness. The following figure shows the input ground truth.

<figure><img src="image_130.png" /><figcaption>**Figure 5.8** – Example of input data : from left to right, the first three bands (RGB) of the Airphen multispectral camera, three edge classes and the inner of individual class</figcaption></figure>
**5.3.1.1 Upstream of the network**

It is composed of three **upstream** modules : Initial Image Transforms (IIT), Input Band Filter (IBF) and Universal Function Approximator (UFA).

**Initial Image Transforms (IIT)** In order to enrich the pool of information, spectral band transformations are added to take into account specific spatial gradients in the image and spectral mixing. Seven important transformations are considered. The standard deviation between spectral bands, noted ρ_std can help to detect the spectral mixture. For example, between two different surfaces like ground and leaf (which have opposite spectral radiance), the spectral mixing makes a pixel with linear combination, thus the standard deviation tends to zero [Louargant et al., 2017]. Three Gaussian derivatives on different orientations are computed. Gxx, Gxy and Gyy filters on ρ_std give an important spatial information about the breaks of gradient, and therefore about the outer limits of surfaces. The Laplacian, the minimum and maximum eigenvalues (of the Hessian matrix) also called ridge of the ρ_std seems to easily detect fine elements [Lin et al., 2014], such as monocotyledons for vegetation images. All these transformations are concatenated to the channel-wise normalized spectral band input and build the final input image. In the end we have six spectral images plus seven transformations for a final image composed of 13 channels.

**Input Band Filter (IBF)** To remove unneeded parts of the signal, low-pass, high-pass and band-pass filters are added. To apply the low-pass filter we use the equation $z = max(x-a, 0)/(1-a)$ and thus it allows to suppress low values. For the high-pass filter we apply the equation $w = max(b-x, 0)/b$ to suppress the high values. The band-pass filter is the product of low-pass and high-pass filters defined by $y = z * w$. The output layer is the concatenation in the channel-wise of the input images, the low-pass, the high-pass and the band-pass filter which produce 4\*13 = 52 channels. Finally, to reduce the output data for the rest of the network, a bottleneck is inserted using a 3 × 3 convolution layer, and it generates a new tensor with 16 channels.

**Universal Function Approximator (UFA)** To separate efficiently leaves from the background and to learn spectral features, a Universal Function Approximator is defined on the upstream of the network (Figure 5.9). The UFA is based on Taylor expansion theorem, an approach to learn this form of development in deep-learning is called DenseNet and then corresponds to the sum of the concatenation of the signal with these spatio-spectral derivatives. This was successfully used for vegetation segmentation [Vayssade et al., 2021]. Three parameters, such as the *depth* (number of convolutions), the *width* (number of filters denoted *W*) and *k* (kernel size) configure the network and were empirically fixed to *depth* = 3, *width* = 16 and *k* = 1. An auxiliary output is used here to maximize the class similarity on the upstream of the network and to extract pure spectral information.

<figure><img src="image_131.png" /><figcaption>**Figure 5.9** – Universal function approximator based on DenseNet. Multiple arrows shows concatenation as input layer</figcaption></figure>
### 5.3.1.2 Core network model

The proposed core CNN architecture is based on an advanced U-Net architecture named MFP-Unet (multi-feature pyramid U-net) proposed by [Moradi et al., 2019]. It is a neural network composed by several 2D convolution layers between different sub-scales. At each sub-scale, the spatial size is divided by two and the number of filters is multiplied by two. Sub-scales are obtained by Max-Pooling layers. Then, to retrieve the original size a 2D UpSampling layer is used. In this study the depth of the U-Net model is fixed to three down-scale (size 512, 256, 128, 64). The specificity of MFP-Unet is that all sub-scale feature maps are directly up-sampled to the initial size, concatenated to the channel-wise and then used for the classification (Figure 5.10). In addition according to [Morris, 2018] and [Xie et al., 2020] an auxiliary loss function is put at each sub-scale feature layer and it enforces the learning of edges at different scale, making the network more robust to spatial resolution. Losses at each prediction also shorten the back propagation path leading to faster convergence. All convolution layers use a kernel size of 3 × 3 and are followed by a Batch Normalization and a sigmoid activation function ([Moradi et al., 2019, Nwankpa et al., 2018]).

<figure><img src="image_132.png" /><figcaption>**Figure 5.10** – Synthesis of the core network based on MFP-Unet. Red arrows shows Max-Pooling. Blue arrows shows Conv+UpSampling. Black arrows shows Conv. Sub-scale ratio are labelled on corresponding layer. The final output is the concatenation of features UpSampling. Multiple arrows shows concatenation as input layer</figcaption></figure>

### 5.3.1.3 Downstream of the network

It is composed of three **downstream** modules : CoordConv and Universal Function Approximator (UFA) and Classification.

**CoordConv** The core network model produces a concatenation of 4 layers of 16 features (4 × 512 × 512 × 16) which results of a layer of size 1 × 512 × 512 × 64. A coordinate layer [Liu et al., 2018a] is also concatenated and allows to consider the mapping between the coordinates in (x,y) Cartesian space to one-hot pixel space. Three variables are appended, the normalized x and y coordinate and the radial coordinate $\sqrt{(x - 0.5)^2 + (y - 0.5)^2}$. This module improves the results removing noise, ground moisture and it fixes few small holes in the segmentation mask.

**Universal Function Approximator (UFA)** This UFA – also presented on the upstream – is used to accurately mix features coming from various scales as well as the Cartesian coordinates. This module reconstructs the mapping function from the Cartesian space to a spatio-spectral feature of size 1 × 512 × 512 × 16. However, in contrast with the upstream UFA, the kernel size *k* was fixed to 3 to take into account neighboring.
**Classification and Auxiliary output** The classification is done through a small network composed of two 1 × 1 convolution layers. Followed by a “Pyramid Pooling Module” to consider different scales before the outputs and smooth the boundary prediction. [Zhao et al., 2017] showed that fusing the low to high-level features improved the segmentation task. It consists in the sum of different 2D convolutions whose kernel sizes have been set to 3, 5, 7 and 9. The number of filters is the same as the number of classes : 4 (defined in 5.3.1). The result of each convolution is concatenated and the final image output is given by a 2D convolution. In addition all convolutions are followed by a Batch Normalization and a sigmoid activation function. The figure 5.11 shows that sub-network.

<figure><img src="image_133.png" /><figcaption>**Figure 5.11** – Auxiliary output and classification module. Multiple arrows shows concatenation as input layer</figcaption></figure>

### 5.3.2 Loss function

A wide variety of loss functions have been developed during the emergence of deep-learning. Recently, [Rahman and Wang, 2016] proposed a solution to optimize an approximation of the mean Intersection over Union (mIoU) which seems to be optimal for binary segmentation [Zhou et al., 2019]. The loss function using ground truth (p) and the prediction (p̂) is defined by :

$$ \mathbf{mIoU}(p, \hat{p}) = 1 - \frac{p\hat{p}}{p + \hat{p} - p\hat{p}} \qquad (5.1) $$

This loss function was used on each auxiliary. In addition the loss is computed separately on each class, weighted (denoted W<sub>C</sub>) and summed. The result of this function is :

$$ \mathbf{Aux}(p, \hat{p}) = \sum_{C=0}^{4} W_{C} \times \mathbf{mIoU}(p_{C}, \hat{p}_{C}) \qquad (5.2) $$

In the above equation, the weight W were empirically set to [0.175, 0.526, 0.211, 0.09]. Meaning that the second class (inner edges) is prioritized (to separate inner instance). Then it is outer edges that allow to separate the “big” leaves from small or thin leaves. Finally the thin leaves (mostly spectral mixing) and inner big leaves (essentially vegetation minus boundaries) have the smallest weight values because these classes should be easier to learn.

In recent CNN architectures for instance segmentation, the loss function does not take into account the number of detected instances or the shape of the segmentation. This aspect is only evaluated after learning, e.g., using a symmetric best dice metric. This implies that we can not guaranty that the network is well learned on crowded scene, where instance is generally merged. One problem is that until recently, instances could not be retrieved directly during the learning phase, this is due to the “non-maximal suppression” algorithm required after the CNN or the time required for the association between the detected instances and the ground truth.
In this paper, we introduce a new loss function considered at the downstream of the network. The main idea is to take into account an approximation of the segmentation quality of each leaf. To estimate the segmentation quality, the undetected, under/over segmented and fused objects can be evaluated, trough a sorted histogram of the number of pixels associated to each connected component for both prediction and ground truth, as showed in the next figure 5.12 for both prediction (orange) and ground truth (blue).

<figure><img src="image_134.png" /><figcaption>**Figure 5.12** – Sorted histogram of the number of pixels associated to each component, the blue (resp. red) line represents the true (resp. predicted) number of pixels for each component.</figcaption></figure>

This figure shows that when the prediction curve is lower (orange) than the ground truth curve (blue), it means that there is over-segmentation. This arise when a bigger element is split in two, when the borders of the shape are trimmed or when zero pixels of the shape are detected (undetected). On the other hand, when the prediction curve is higher than the ground truth, it means that the contours of the shape are roughly detected or if it greatly exceeds the ground truth, then we are in the presence of fused shapes. Based on these curves, a loss function can be construct the deal with over and under segmentation on shape criterion, thus a custom absolute percentage error was defined :

$$ \text{Leaf}(p, \hat{p}) = \frac{\sum_{i=0}^{N} |\text{leaky\_relu}(h_i(p) - h_i(\hat{p}))|}{h_{max}(p) + 1} \qquad (5.3) $$

$$ \text{leaky\_relu}(x) = \begin{cases} x & \text{if } x \ge 0 \\ x * 0.2 & \text{if } x < 0 \end{cases} \qquad (5.4) $$

On the above equation N is the number of components, $h_i(p)$, and $h_i(\hat{p})$ is respectively the number of pixels of a connected component in the ground truth and on the predicted segmentation within the sorted histogram. While $h_{max}(p)$ is the maximum number of pixels of a component in the ground truth.
The leaky_relu, is used to explicitly prioritise the learning on under-segmentation rather than over-segmentation which allows to prioritize the merged objects. This was defined because conventional losses did not give good results in dense vegetation cover, causing a large segmented area that is detected as a single entity. Note also that over-segmentation is less problematic, since it occurs mainly around the borders of the leaves, which are recovered later, through a watershed algorithm (section 5.3.3). It is the first study to suggest this type of loss.

The downstream loss is defined by DownAux($p, \hat{p}$) = Aux($p, \hat{p}$) + Leaf($p, \hat{p}$). Finally the global loss considers the upstream auxiliary loss, each of the 4 feature scale auxiliary loss and the downstream loss. Thus the global loss is defined as the weighted sum of all auxiliaries losses (in the same order) where the weights $W$ were empirically defined with $W = [0.01, 1.0, 0.1, 0.1, 0.01, 4.0]$.

### 5.3.3 Refinement with vegetation mask and watershed

The used U-Net architecture is good at detecting “big” elements on the scene but lacks precision on very small and thin elements. A method to produce an optimized vegetation mask was proposed in a previous work [Vayssade et al., 2021]. Using this mask provides better performance than adding a specific class to the network. Thus our previous work is used here to get the best foreground/background segmentation mask as input of the watershed. It is also learned on a specific dataset with more illumination conditions and should be more robust, especially for thin elements.

The proposed network generates 4 classes. The first two are associated to “big” leaf boundaries (denoted Edges = Outer + Inner). The third one is small and thin leaves (denoted Thin), and the fourth is the inner of big leaves denoted Big. The watershed seed is defined with the following equation 5.5 to generate the seed mask and can be seen in the figure 5.13 :

$$Seed = Thin + Big - Edges \qquad (5.5)$$

<figure><img src="image_135.png" /><figcaption>**Figure 5.13** – The refinement of the CNN output through watershed and Deep-Indices. The Seed of the equation 5.5 can be seen on the “watershed seed” image.</figcaption></figure>
### 5.3.4 Training setup

The image dataset is randomly split into a training set (80%) and a validation set (20%). However the initial seed is kept for reproducibility. In addition a random data augmentation is used during the training to increase the dataset variability. A random vertical and horizontal flip is considered as well as a perlin simplex noise [Bae et al., 2018] (of size 512 × 512), set with 2 modal in range [0.7 – 1.3] which multiplies the number of input images. Low values add shadows, while high values add specular effects. The training is done through Keras module within Tensorflow 2.6.0 framework. All computations are done on an NVidia RTX 3060 which have 12GiB of memory and due to the size of the network only one image is computed at once.

### 5.3.5 Evaluation metrics

There is a large number of possible evaluation metrics for instance segmentation, which have been reviewed in [Scharr et al., 2016]. However, we decided to keep the evaluation metrics used in the Leaf Segmentation Challenge (LSC) as a reference for comparison [Scharr et al., 2017, Kulikov et al., 2018, Bell and Dee, 2019, Ward and Moghadam, 2020]. Therefore, we use the foreground/background DICE metric to evaluate the separation of large leaves from the ground (denoted *FgBgDice*). We estimate the average accuracy of leaf segmentation by the symmetric best DICE score among all objects (leaves). The best DICE score among all objects (leaves) to estimate the average accuracy of leaf segmentation is denoted *BestDice*. The *AbsDiffFG* estimates how good the algorithm is at segmenting the leaves. And finally, *DiffFG* estimates the efficiency of the algorithm for counting leaves. *SBD* for Symmetric best DICE is extract to estimate the average leaf segmentation accuracy. All these metrics are common and presented by [Scharr et al., 2016]. In addition, to compare datasets results, a new metric is introduced named *NAbsDiffFG* defined by the division of *AbsDiffFG* and the mean of leaves count in the dataset.

## 5.4 Results

As previously defined, all datasets were split into training and validation sets with a ratio of 80-20%. Once training is done on the defined setup and using the loss function, the evaluation metrics are extracted : the *FGBGDice* metric is used to evaluate the soil versus vegetation segmentation, while the *BestDice* metric is used to evaluate the instance segmentation. Each connected component is associated to its best corresponding ground truth based on a dice score. Then the metric is defined by the mean dice score of the best match. Other metrics show relevant information about over and under segmentation. The following subsections are dedicated to each dataset, from the simplest to the most difficult to segment.

### 5.4.1 Komatsuna dataset

The first one is the Komatsuna dataset. It is composed of RGB data of growing Japanese Mustard. This dataset is interesting for its controlled illumination conditions. In addition, the number of leaves it quite the same for all the growing stages, ranging from 3 to 6 leaves. This is important regarding our loss function which take into account the quantity and quality of each segmented leaf. However, this dataset does not contain small and thin leaves, thus the third class was replaced by a foreground class. The following Table 5.1 shows the evaluation on this dataset.
<table><caption>**Table 5.1** – Evaluation metrics for the Komatsuna dataset</caption><thead><tr><td>Metric</td><td>Training</td><td>Validation</td></tr></thead><tbody><tr><td>FgBgDice</td><td>0.9732</td><td>0.9715</td></tr><tr><td>BestDice</td><td>0.8796</td><td>0.8565</td></tr><tr><td>SBD</td><td>0.8713</td><td>0.8517</td></tr><tr><td>DiffFG</td><td>-0.2639</td><td>-0.4444</td></tr><tr><td>AbsDiffFG</td><td>0.3944</td><td>0.5222</td></tr><tr><td>Number of images</td><td>720</td><td>180</td></tr><tr><td>Mean leaf count</td><td>4.9194</td><td>4.9722</td></tr><tr><td>NAbsDiffFG</td><td>0.0802</td><td>0.1050</td></tr></tbody></table>

The result on this dataset shows a relatively perfect score of the soil and vegetation segmentation with a *FGBGDice* of 0.97. Few errors remain as shown in figure 5.14, a green bottle is visible on the left side of some images, which seems to indicate that the green component has a major impact on vegetation detection. Since RGB data are used, this error should be visible in most studies. This interesting element also demonstrate that the used *CoordConv* layer does not play its role in the spatial management of this element.

<figure><img src="image_136.png" /><figcaption>**Figure 5.14** – Visual comparison of Komatsuna dataset after the training, the red arrow in the input image shows a green bottle</figcaption></figure>
As shown by the *DiffFG* score, our method mostly detects the right number of leaves. This score is small and negative and shows that few leaves are split, which occurs when a big leaf mostly overlaps a smaller one. A visible bottle also has an impact on the other metrics, shown in figure 5.14. However, most of the errors for the *BestDice* comes from under and over-segmentation. The under-segmentation occurs on the leaves stems. The stem that connects the leaf to the plant is usually not well detected, and this can play a significant role in lowering the *BestDice* and *SymmetricBestDice* scores. The stem can be undetected (under-segmentation) or may be identified as another leaf (over-segmentation). It can also cover some small leaves and divide them in two, which is expected in pixelwise segmentation. To benchmark our study, the following table 5.2 reports the results of few previous studies using the same metrics.

<table><caption>**Table 5.2** – Comparison of the ResNet 101 based on Detectron II framework of facebook and the UPGen based on Mask-RCNN. Our solution is the most effective in both metrics</caption><thead><tr><th>Study</th><th>SBD</th><th>AbsDiffFG</th></tr></thead><tbody><tr><td>Our method's results (2021)</td><td><strong>0.8565</strong></td><td><strong>0.5222</strong></td></tr><tr><td>[Ward and Moghadam, 2020] ward 2020</td><td>0.7776</td><td>-</td></tr><tr><td>[Gomes and Zheng, 2020] gomes 2020</td><td>0.7700</td><td>0.8800</td></tr></tbody></table>

As comparison our results tackle one of the states of the art ResNet 101 based on Detectron II framework of facebook [Gomes and Zheng, 2020]. And the proposed UPGen methods proposed by [Ward and Moghadam, 2020].

### 5.4.2 Leaf Segmentation Challenge dataset

The second dataset is interesting for its controlled illumination conditions. Unlike the Komatsuna dataset, this one includes a greater diversity of leaf quantity, between 2 and 16 leaves. This is important to show the robustness of our loss function. In the same regard as the Komatsuna section, this dataset does not contain small and thin leaves. Thus third class was replaced by a foreground class, which is also used by the competition evaluation. In addition, the testing dataset is not publicly available, thus a separated evaluation was performed on the online competition website<fnref n="1" />. The training and validation are defined by splitting the publicly available dataset. This step is used to reduce the over-fitting and retrieve stable parameters. The next table summarizes the number of images for each training, validation and test datasets.

<table><caption>**Table 5.3** – Number of images for each sub LSC dataset</caption><thead><tr><th>Dataset</th><th>A1</th><th>A2</th><th>A3</th><th>A4</th><th>Total</th></tr></thead><tbody><tr><td>Training</td><td>102</td><td>25</td><td>22</td><td>499</td><td>648</td></tr><tr><td>Validation</td><td>25</td><td>6</td><td>5</td><td>125</td><td>161</td></tr><tr><td>Test</td><td>33</td><td>9</td><td>56</td><td>168</td><td>266</td></tr></tbody></table>

<footnote n="1">1. competitions.codalab.org/competitions/18405</footnote>
The following Table 5.4 shows the overall evaluation on this dataset (merged A1, A2, A3, A4). The test results can be found in the online leader-board. The results are fairly the same between training, validation and test datasets in term of *FGBGDice* and *BestDice*. However a notable divergence with the testing dataset for the others metrics is visible. The *DiffFG*, and *AbsDiffFG* for training and validation indicate an over-segmentation, while the test dataset highlights an under-segmentation. This over-segmentation probably enhances the *SBD* score while the *BestDice* remains stable.

<table><caption>**Table 5.4** – Evaluation metrics of the overall LSC dataset</caption><thead><tr><th>Metric</th><th>Training</th><th>Validation</th><th>Test</th></tr></thead><tbody><tr><td>FgBgDice</td><td>0.9489</td><td>0.9522</td><td>0.9489</td></tr><tr><td>BestDice</td><td>0.7659</td><td>0.7707</td><td>0.7608</td></tr><tr><td>SBD</td><td>0.7608</td><td>0.7655</td><td>0.8047</td></tr><tr><td>DiffFG</td><td>-1.9634</td><td>-2.1223</td><td>3.5628</td></tr><tr><td>AbsDiffFG</td><td>2.1707</td><td>2.3050</td><td>6.1257</td></tr><tr><td>Number of images</td><td>648</td><td>161</td><td>266</td></tr><tr><td>Mean leaf count</td><td>13.8390</td><td>13.1463</td><td>-</td></tr><tr><td>NAbsDiffFG</td><td>0.1568</td><td>0.1753</td><td>-</td></tr></tbody></table>

As discussed in the materials and data section, this database is composed of four sub-databases with different cameras and plants. The learning was done independently on each of them, but they were merged for the presentation in the previous table. The test on the online evaluation platform returns the results for each of them, summarized in the table 5.5. First of all, the sub-datasets A2 and A3 contain a very small amount of images, respectively 25 and 22 for the training sets. Moreover the A3 dataset is composed of images of size 2048 × 2048 re-scaled to 512 × 512. This imply a huge loss of information, especially on the leaf boundaries. In addition, this A3 dataset contains hard shadows. These three factors explain most of the errors for this dataset. The quality, acquisition conditions, and plants of the other three data sets (A1, A2, A4) are similar. Only the background and the amount of images mainly changes, which is also reflected in the evaluation score. The figure 5.18 also shows unlabeled leaves in the center of the plant, moreover they are even defined as background.

<table><caption>**Table 5.5** – Evaluation metrics of the independent sub LSC dataset</caption><thead><tr><th>Metric</th><th>A1</th><th>A2</th><th>A3</th><th>A4</th></tr></thead><tbody><tr><td>FGBGDice</td><td>0.9641</td><td>0.9294</td><td>0.9692</td><td>0.9416</td></tr><tr><td>BestDice</td><td>0.6905</td><td>0.6944</td><td>0.6783</td><td>0.7964</td></tr><tr><td>SBD</td><td>0.8047</td><td>0.7895</td><td>0.7317</td><td>0.8294</td></tr><tr><td>DiffFG</td><td>2.0909</td><td>1.7778</td><td>21.2857</td><td>0.0892</td></tr><tr><td>AbsDiffFG</td><td>3.0000</td><td>3.1111</td><td>14.9285</td><td>1.6250</td></tr></tbody></table>

For all sub-dataset the *FGBGDice* score is slightly less than the Komatsuna dataset. Most pictures of the A1 present wide area of green moisture on the ground, visible on the Input inside the figure 5.15. A2 contains few very small weeds and few small spots of moisture, but the performances of *FGBGDice* is mainly due to the quantity of images for the learning, as showed by the figure 5.16 boundaries are miss-classified. For the same reason, it is also visible for the A3 dataset shown by the figure 5.17. It's also important to notice that leaves from an outside plant are visible in few images and detected by the algorithm, however these leaves are unlabelled, resulting an over-segmentation.
<figure><img src="image_137.png" /><figcaption>**Figure 5.15** – Visual example of LSC results for A1</figcaption></figure>

<figure><img src="image_138.png" /><figcaption>**Figure 5.16** – Visual example of LSC results for A2</figcaption></figure>
<figure><img src="image_139.png" /><figcaption>Figure 5.17 – Visual example of LSC results for A3</figcaption></figure>

<figure><img src="image_140.png" /><figcaption>Figure 5.18 – Visual example of LSC results for A4</figcaption></figure>
The following table 5.6 reports method's results of different previous studies, to benchmark our method's results on this LSC dataset.

<table><caption>**Table 5.6** – Comparison of our solution with state-of-the-art challengers in this dataset.</caption><thead><tr><th>Study</th><th>SBD</th><th>AbsDiffFG</th></tr></thead><tbody><tr><td>[Ward and Moghadam, 2020] ward 2020</td><td><strong>0.8800</strong></td><td>-</td></tr><tr><td>[Kulikov et al., 2018] kulikov 2018</td><td>0.8040</td><td>2.00</td></tr><tr><td>[Gomes and Zheng, 2020] gomes 2020</td><td>0.7700</td><td><strong>0.88</strong></td></tr><tr><td>Ours method's results (2021)</td><td>0.7608</td><td>6.12</td></tr><tr><td>[Pape and Klukas, 2015] pape 2015</td><td>0.7440</td><td>2.60</td></tr><tr><td>[Scharr et al., 2016] scharr 2016</td><td>0.6830</td><td>3.80</td></tr></tbody></table>

These results show that the studied method is less efficient than [Ward and Moghadam, 2020, Kulikov et al., 2018] and [Gomes and Zheng, 2020]. Due to the small size of the training sample (102, 25, 22, 499 images for data subsets A1, A2, A3, A4, respectively), it can be assumed that this is due to the increase in data they use in order to expand their data set. It seems the data augmentation based on Perlin noise is insufficient. Nevertheless, our method is better than the following two : [Pape and Klukas, 2015, Scharr et al., 2016], since they don't use any data augmentation. However it can be noted that the *AbsDiffFG* value for the studied method is much higher than for the others. This implies an important over-segmentation, resulting in particular from the quality of the A3 dataset as it can be seen in the table 5.5.

### 5.4.3 Airphen dataset

The last one presented is our multispectral dataset, which contains variable acquisition conditions (sunny, cloudy, rainy, etc.), a variable number of leaves (from a few to hundreds), and contains very small leaves that touch or overlap the others. The next table 5.7 show the results of our method applied our dataset.

<table><caption>**Table 5.7** – Evaluation metrics for the Airphen dataset</caption><thead><tr><th>Metric</th><th>Training</th><th>Validation</th></tr></thead><tbody><tr><td>FgBgDice</td><td>0.9791</td><td>0.9785</td></tr><tr><td>BestDice</td><td>0.6744</td><td>0.6678</td></tr><tr><td>SBD</td><td>0.6670</td><td>0.6638</td></tr><tr><td>DiffFG</td><td>-41.3583</td><td>-49.4667</td></tr><tr><td>AbsDiffFG</td><td>46.7500</td><td>53.2333</td></tr><tr><td>Number of images</td><td>240</td><td>60</td></tr><tr><td>Mean leaf count</td><td>265.83</td><td>295.80</td></tr><tr><td>NAbsDiffFG</td><td>0.1758</td><td>0.1799</td></tr></tbody></table>

It can be seen that the *FGBGDice* score shows adequate soil/vegetation segmentation, as demonstrated in our previous study [Vayssade et al., 2021]. The scores *DiffFG* and *AbsDiffFG* show the presence of over-segmentation. This is due to two issues : the presence of small leaves and the advanced stage of mustard development. Indeed, the method sometimes confuses small leaves and large leaves, which implies in some cases under-segmentation for monocotyledons or over-segmentation for dycotyledons, probably generated by an imprecise annotation of classes. In the second case, this is due to the proposed loss function that forces over-segmentation, resulting in the detection of leaflets instead of leaves in the case of advanced mustard develop-
ment, as showed in the figure 5.19. These errors imply that the *BestDice* and *SBD* scores are not as good as on the Komatsuna dataset and LSC Challenge, studied previously.

<figure><img src="image_141.png" /><figcaption>**Figure 5.19** – Example of the visual results for the Airphen dataset</figcaption></figure>

To show the difference and complexity of this dataset, the method developed by [Ward and Moghadam, 2020] was also learned. This method is based on the Mask-RCNN. However, the method uses data augmentation based on a library that only supports 8-bit unsigned integers. Thus, data augmentation was disabled because the dataset uses a 32-bit float format. The following table 5.8 shows these results.

<table><caption>**Table 5.8** – Evaluation metrics for the Airphen dataset using Mask-RCNN</caption><thead><tr><td>Metric</td><td>Training</td><td>Validation</td></tr></thead><tbody><tr><td>FgBgDice</td><td>0.6266</td><td>0.6111</td></tr><tr><td>BestDice</td><td>0.2271</td><td>0.2157</td></tr><tr><td>SBD</td><td>0.2266</td><td>0.2149</td></tr><tr><td>DiffFG</td><td>-186.0292</td><td>-217.0667</td></tr><tr><td>AbsDiffFG</td><td>186.5542</td><td>217.1333</td></tr><tr><td>NAbsDiffFG</td><td>0.7018</td><td>0.7341</td></tr></tbody></table>

These results show that the method proposed by [Ward and Moghadam, 2020] does not correctly detect leaves in a dense foliage environment. The *FgBgDice* score shows that the soil/vegetation discrimination is weak. As well as the *BestDice* and *SBD* scores which shows poor detection of individuals segmentation mask. The *DiffFG* and *AbsDiffFG* scores confirm these results and show a large amount of undetected elements. As announced in the introduction, the methods based on object detection uses bounding box regressions and non-max suppression which strongly limit the detection in dense environment. Moreover the part of the network allowing to obtain the segmentation mask uses a fixed low resolution resulting in coarse segmentation masks. As shows in the figure 5.20.
<figure><img src="image_142.png" /><figcaption>Figure 5.20 – Example of the Watershed output versus Mask-RCNN for Airphen dataset</figcaption></figure>

## 5.5 Discussion

In this study, we used different datasets. The Komatsuna dataset is composed only of mustard leaves at different stages of development acquired indoors under controlled acquisition conditions with a single camera. There are between 3 and 6 leaves per image for an average of 4.9 leaves. The LSC dataset is decomposed into three sub-datasets of Arabidopsis-Thaliana and one of Rosette plants all at different stages of development. They are acquired indoors under controlled acquisition conditions with their own camera. There are between 2 and 16 leaves per image for an average of 13.5 leaves. Our dataset is composed of bean plants and weeds. Beans are at a single stage of development while weeds are at different stages of development. The acquisitions were made outdoors under uncontrolled conditions. There are between 4 and 777 leaves per image with an average of 271.83 leaves. This dataset is now available online<fnref n="2" /> and may help other studies working on leaves segmentation in natural and complex situations.

We are interested in the metric *NAbsDiffFG* to compare the performances of our method on the different datasets. We see that for our dataset it is 0.1758, 0.1799 for the LSC dataset and 0.1050 for the Komatsuna dataset. We thus obtain good results for the Komatsuna dataset which is the simplest, we obtain an average result for our dataset which is the most complex and finally weaker results on the LSC dataset. This allows us to deduce that the lighting conditions as well as the number of leaves do not influence the performances. On the other hand, the use of several cameras, as on the LSC dataset, could be at the origin of the decrease in performance due to a variation in spatial resolution from one camera to another. We deduce that our method, designed to be used on complex data such as our dataset, gives results comparable to classical methods on other types of datasets.

The CNN used in our method deliberately gives a coarse segmentation in order to maximize the detection and separation of large leaves. The smaller ones are therefore poorly detected (over-segmentation or undetected). These results are reflected in the fact that there is little over-segmentation and under-segmentation of large leaves due to the loss function introduced which takes into account the number of elements detected and their surface. However two problems are raised with the use of this loss function. First of all, it is subject to error jumps in order to avoid merging objects, which makes its optimization difficult and thus requires several learning sessions for an optimal result.

<footnote n="2">2. https://data.inrae.fr/dataset.xhtml?persistentId=doi:10.15454/JMKP9S</footnote>
On the other hand, avoiding merging causes over-segmentation in some cases, such as the mustard in our dataset that detects leaflets. At this time we cannot say whether this is problematic. It is possible that the leaflet scale is more relevant in terms of segmentation than the leaf for advanced stages of plant development.

With the watershed used next, we improve the coarse detection provided by the CNN by extending the detected areas to the ideal soil/vegetation segmentation mask defined by the method developed in a previous study [Vayssade et al., 2021]. This step allows the identification of the smallest leaves. Although the smallest elements, with a diameter of 0 – 2 pixels, are still difficult to detect due to spectral mixing [Louargant et al., 2017]. This is an important result as it allows to configure an acquisition system for a specific minimum size of leaves to detect. In this work, the remaining small leaves are already segmented as vegetation and can be considered as stable due to their sizes. Some features are still relevant to be extracted (such as the size of the leaf and its distance from the crop row) and should help in discriminating crop from weeds, other may suffer from the spectral mixing and become irrelevant (spectral signature, texture...).

There is still an important type of error. Indeed, a single pixel can be at the origin of errors leading to a bad fusion between two areas. It would therefore be interesting to study structural analysis methods in order to overcome these deficiencies. A possible method would be to vectorize the contours and look for an algorithm to reconnect or split the contours for instance according to convexity singularities.

The main interest of the proposed method is its efficiency on mixes of plant species, acquired with natural light. It will be integrated in a processing chain dedicated to the discrimination of crops and weeds in agronomic scenes. Indeed, the detected leaves can be classified according to a large set of criteria (spectral signature, morphological characteristics, texture, distance from the crop row). The underlying hypothesis is that these criteria are more stable at the scale of the leaf than at the scale of the plant. However, this approach has certain limitations when leaves overlap others, as the detected shapes would be heterogeneous. Multifoliate leaves could also be difficult to characterize. In that case, detecting leaflets instead of leaves may be more relevant.

## 5.6 Conclusion

The presented work shows that the CNN network enhances the quality of the segmentation based on multispectral images. Indeed, the background is well removed due to the upstream network with IBF and UFA modules with an accuracy of 95 – 98% of mIoU. Our method is effective in the majority of cases, such as the segmentation of unifoliate leaves like Arabidopsis-Thaliana and early developmental stages of plants. However, our method is not effective in the advanced stages of plant development, especially on mustard which has highly segmented leaves. In this case our method detects leaflets instead of leaves. Their identification is nevertheless relevant for the phenotyping or classification of weeds.

We have seen that the developed method presents better SBD performances +[1.66 – 7.89]% compared to studies that do not use data augmentation. However, on small datasets, it presents lower SBD performances –[4.32 – 11.92]% compared to recent studies that use it. To make the method more consistent we will therefore focus on data augmentation. To improve the detection of classes it would be interesting to improve the upstream modules with attention mechanisms.
This method would allow to correct the illumination of the images. As well as the use of
a *MIRNET* [Zamir et al., 2020] would allow to eliminate the noise. Finally, to improve the
downstream modules, the use of *Selective Kernel Convolution* [Zamir et al., 2020] would allow
a better fusion of the multi-scale information instead of *Universal Function Approximator*. We
will then try to minimize the under-segmentation detected in our method, by using the *Deep
Watershed Transform for Instance Segmentation* method [Bai and Urtasun, 2017].

A small performance loss for the developed method is seen on our dataset compared to Komatsuna and LSC dataset. It would be interesting to evaluate the reason(s) leading to this loss as it may come from the uncontrolled acquisition conditions, the multispectral nature of the images, the size differences between leaves or the important number of leaves in each images.

In all cases, this approach should lead to an enhancement of features extraction which
may improve crop/weed classification. These increased performances would lead to a better
tracking of the weed flora. These algorithms show promising and robust results in natural
acquisition conditions. The segmentation results are obtained fast enough to be used in a real-
time crop/weed discrimination setting and could be embedded on a Unmanned Ground Vehicle
(UGV) for quick and localised intervention. Applied on images acquired from an Unmanned
Aerial Vehicle (UAV) this could be used for tool assisted plot management to help farmer in
their decision making.

**5.7 Further research**

We defined the CNN architecture from the state of the art, adding components in an increas-
sing development cycle. Thus, we have not established the contribution of each block on this
paper. The first 3 modules (iit, ibf, ufa) and the last sprb modules have been widely developed
in our previous study about vegetation segmentation, each module have showed a contribution.
But further research is undergoing to show the impacts of different modules on the upstream
and downstream of the network, and the proposed loss function, and loss weights. This task
must be automated trough an hyper-parameter optimization process which we did'nt explored
wet. And the watershed algorithm used to refine the output must be added to the neural net-
work, and part of this optimization process. Finally more advanced data augmentation could
be explored.

**Data availability**

The Airphen dataset used in this study have been released and it's publicly available at
https://data.inrae.fr/dataset.xhtml?persistentId=doi:10.15454/JMKP9S, using Creative
Common CC0 1.0 Public Domain Dedication licence

**Funds & Conflicts of Interest**

This project has received funding from the European Union's Horizon 2020 research and innovation program [grant agreement ID: 727321] (project acronym: IWMPRAISE) and from the French National Research Agency Challenge RoSE [grant agreement ID: ANR-17-ROSE-0002] (project acronym: ROSEAU). The authors declare no conflict of interest.
## 5.8 Conclusion de chapitre

Les résultats de cette étude montrent que l'approche proposée pour la segmentation des feuilles est intéressante sur plusieurs jeux de données accessibles en ligne et sur nos données. Elle génère d'excellents résultats sur les données de Komatsuna. En revanche, la méthode n'est pas parmi les meilleures lorsqu'elle est appliquée au jeu de données du LSC (Leaf Segmentation Challenge), mais elle reste quand même compétitive. Pour nos données multispectrales, dans le cadre du challenge RoSE, les résultats montrent qu'une approche traditionnelle, telle que les mask R-CNN qui reposent sur une utilisation de l'algorithme NMS ne permettent pas de détecter efficacement les individus en environnement dense. Contrairement à la méthode qui a été développée dans cet article. Cependant, bien que ces résultats représentent une avancée dans le domaine du "pixelwise instance segmentation", des recherches complémentaires doivent être entreprises.

Concernant l'échelle de segmentation (i.e. la feuille), ces résultats suggèrent qu'une entité légèrement plus fine serait préférable d'un point de vue segmentation. En effet la visualisation des résultats sur les données du challenge RoSE, montre que les feuilles unifoliées/unilobées sont bien détectées, ce qui n'est pas le cas des feuilles multifoliées ou multilobées. Ainsi, l'échelle foliole serait probablement plus opportune.

Le présent chapitre 5 et le chapitre précédent 4 constituent l'étape de segmentation de notre pipeline. Maintenant, que l'on est capable de détecter les individus, nous nous intéresserons à l'extraction de propriétés et la classification, comme nous pouvons l'observer sur la figure ci-dessous.

<figure><img src="image_143.png" /></figure>

Rappel de l'architecture de notre chaîne fonctionnelle
ici les premières étapes en vert sont terminées, et en jaune les étapes suivantes étudiées
# 6 Extraction de propriétés

Ce chapitre aborde le cœur de cette thèse “Approche multicritère pour la caractérisation des adventices”. C’est-à-dire l’extraction et l’exploitation de propriétés permettant la discrimination culture/adventices. Ces dernières, reposent sur l’exploitation de divers descripteurs comme les textures, les formes, ou les propriétés spatiales (présentées en section 2.4). Ces descripteurs sont couramment utilisés [Griepentrog et al., 2005, Gao et al., 2013, Xue et al., 2017, Louargant et al., 2018, Su, 2020, Mekhalfa and Yacef, 2021, Wu et al., 2021a] pour caractériser et classer les plantes détectées dans une image. En revanche, les études utilisant ces descripteurs sont basées sur un petit ensemble de propriétés – souvent moins de 20 – ainsi que sur des cultures et des matériels d’acquisition non concordants entre les études. Ainsi, ces études ne permettent pas une estimation précise des potentiels des différents descripteurs, appliqués à un sujet donné. Il est donc difficile de trouver la meilleure combinaison de critères et leur impact sur la qualité de la classification à partir d’une bibliographie.

Parallèlement, dans le chapitre précédent 5, concernant la segmentation, nous avons défini une approche alternative, reposant sur la détection des feuilles et non des plantes. Cela permet la caractérisation du couvert végétal d’une façon plus fine que ce qui est actuellement établi. L’objectif de l’étude qui suit est donc d’analyser, d’intégrer et d’évaluer les spécificités de chaque critère extractible à l’échelle de la feuille, à partir d’images prises en conditions terrain.

Un ensemble de 25 types de propriétés a été défini, dont 3 spatiaux, 8 de forme et 14 de texture. Chaque type de propriété regroupe un sous-ensemble de propriétés partiellement corrélées entre elles. Par exemple le type “shape particle” contient un sous-ensemble de 18 propriétés de forme, tandis que le type “texture haralick” est composé de 416 propriétés de texture. Le total de propriétés extraites dans cette étude monte à 3545. Étant donnée la quantité de propriétés, deux problèmes se posent :

1) Comment optimiser les propriétés comportant des paramètres ? Par exemple calculer des propriétés comme les DoG (Difference Of Gaussian) nécessite de fixer 4 paramètres préalablement (sig, lm, gm, ps). Ceci est un problème d'optimisation des hyperparamètres et dans l'exemple précédent, les 4 paramètres sont dans un espace continu, ce qui signifie qu'une recherche exhaustive ne peut être établie. Pour réaliser ces optimisations, des algorithmes plus complexes doivent être utilisés pour échantillonner l'espace des paramètres efficacement et automatiquement. C'est le premier point qui est abordé dans l'étude qui suit.

2) Comment sélectionner efficacement le meilleur sous-ensemble ? Pour répondre à cette question, les chercheurs dans le domaine de l'agriculture de précision utilisent principalement une métrique pour estimer dans quelle mesure une caractéristique est corrélée au problème, puis les meilleures caractéristiques sont utilisées. Cependant, il peut exister une meilleure combinaison, par exemple, si l'on combine la distance aux rangs et une caractéristique de texture, l'efficacité pourrait être supérieure à celle de deux critères de texture dont chacun est plus discriminants que le rang, car la distance au rang est probablement complémentaire à une propriété de texture. C'est en tout cas ce que suggèrent les études utilisant une distance au rang [Louargant et al., 2018]. C'est le second point qui est abordé dans l'étude qui suit.
# Weed discrimination at leaf scale : towards the characterization of the vegetation cover : A large comparison of shape, spatial and textural features

**Authors** Jehan-Antoine Vayssade, Gawain Jones and Jean-Noël Paoli

**Publication** 08/01/2022 – Computer and Electronics in Agriculture – under-review

**abstract** In the context of crop and weeds discrimination, different methods are used to detect and classify plants from an acquisition system. Various estimators and descriptors like textures, shapes, spatial properties, or vegetation indices are commonly used to characterize plants within an image.

However, the available studies are based on disparate criteria, plants, and acquisition materials. The criteria are often defined by parameters set empirically according to the situation, and the processing scale (cluster, plant, leaf, pixel) has a significant impact on classification performance. Thus, these studies do not allow an accurate estimation of the potential of criteria combinations, if applied to a new study. In addition, the detection of plants remains a difficult problem in image processing. Indeed, plant segmentation approaches are often based on unreliable methods when the coverage rate increases. Currently, an emerging solution is leaf segmentation, allowing to characterize dense vegetation cover.

Thus, the objective of this study is to: (1) experimentally evaluate the discrimination potential of each criterion at the leaf scale, using images taken in field condition; (2) optimize the parameters of these criteria; and (3) determine the best combination of criteria to use.

A literature review is conducted to determine the set of criteria that could be used. A huge set of 3545 criteria is studied. For criteria defined by parameters, an hyperparameter optimization algorithm is involved. Then, an algorithm is defined to select the best subsets of features. This algorithm evaluates the performance of each subset of criteria on a ground truth dataset, using an appropriate performance metric. Finally, a classification of the vegetation cover is proposed, using the best performing subset. Results show the importance of selecting a smaller set of properties (at most 20 features among the 3545 available) and associating different feature types (for instance spatial with textural and morphological features).

**Keywords** leaf classification; crop and weed discrimination; remote sensing; multi-spectral
### 6.1 Introduction

Crop and weed discrimination is a long-standing concern in modern agriculture, as weeds are generally considered a nuisance since they have a possible negative impact on crop yield or harvest quality. To have a reliable knowledge of the risks of yield loss over time, locally acquired experiences on 110 weed control trials under comparable conditions were conducted in France between 1993 and 2015 on three major annual crops (soft wheat, rapeseed and sunflower). A reduction of 25 q/ha, 4 q/ha and 3.5 q/ha was respectively observed [Cordeau et al., 2016]. Weeds are thus competing with their neighbours through competition for nutrients (nitrogen, potassium, phosphorus, carbon), light and water.

To face this problem the main weed control approach is based on the application of herbicide with sprayers. The methods have evolved from uniform chemical weeding of the plot to localized weeding in precision agriculture and the use of new technologies and digital tools. In this impetus, today's precision agriculture is seeking to discriminate even finer elements like the plant in real time. It allows to weed out undesirable individuals, to treat diseased crops or to supply nutrients for crops to reduce the use of synthetic products and the ecological impact of these products.

Current research on crop/weed discrimination are based on criteria extraction and exploitation and provide good discrimination results, up to 77-97% [Ahmed et al., 2016] depending on the criteria, acquisition conditions and observed plants. The extraction and exploitation of such criteria should be included in a computer vision pipeline that follows the key steps of image analysis including preprocessing, segmentation [Vayssade et al., 2022], feature extraction and classification [Perez-Sanz et al., 2017]. Figure 6.1 illustrates this general workflow.

<figure><img src="image_144.png" /><figcaption>**Figure 6.1** – key steps of image analysis for crop/weed discrimination with (a) acquisition enhancement (b) vegetation indices (c) instance segmentation (d) feature extraction (e) classification</figcaption></figure>

Major reviews of weed mapping techniques have been carried out by [Monteiro et al., 2019], [Gao et al., 2018], [Mostajer Kheirkhah and Asghari, 2019] and [Wang et al., 2019a] and show multiple approaches to weed mapping and feature extraction. For example, during the extraction stage, the criteria can be classified according to the nature of the information. We can distinguish (i) shape criteria which are based on different representations of the contours, [Ahmed et al., 2012, Pérez et al., 2000], (2) spatial criteria based on the distance of the individual to another element of the image, such as the distance to the nearest crop row [Louargant et al., 2018], (3) texture features based on the analysis of the spatial distribution of the pixels (may include the soil) [Rojas et al., 2017] and, finally (4) spectral indices at specific points on the surface can also be used [Liu et al., 2019], which could also be considered as a subset of texture features. However, each type of criterion has different effectiveness and limitations. Spectral
and texture features are more affected by illumination variations, while shape features are affected segmentation quality. Moreover, because of these limitations, studies usually combine these types of features to take advantage of all the benefits.

However, although these criteria seem to be efficient, they may loose effectiveness in real conditions because they are limited by the previous stages of the pipeline like in Figure 6.1. In particular by the plant instance segmentation step, which is still a scientific challenge despite the advances in deep-learning. In this context, we have proposed in a previous study to detect not the whole plant but the leaves which gives promising results in terms of instance segmentation and allows to work in dense foliage [Vayssade et al., 2022]. Moreover, the leaf offers more stable criteria : there is less intra-individual variability than with a plant. For example, depending on the growth stage, the plant will contain a variable number of leaves with different orientations and sizes, these leaves may also have a partial overlap between them leading to an infinite number of cases. Thus the extracted criteria are viable for a given stage of development, for given crops and specific acquisition conditions.

We will see that some of these issues do not occur at the leaf level. The extracted criteria are also more understandable at leaf scale : shape criterion, such as length or perimeter is less subject to interpretation on the leaf level than on the plant level. Texture criterion, such as gradients, gives the veins of the leaf, but there's no obvious correspondence at the plant scale. Thus the study of these criteria at the leaf scale could be promising to discriminate leaves between crops and weeds. In addition, such leaf-scale criteria have been used successfully for tree species discrimination, [Cerutti et al., 2013, Cerutti et al., 2014] showing an accuracy up to 90% with only shape criteria at leaf scale for species classification.

Therefore, this study focuses on two points, evaluating criteria that could be used at the leaf level and determining the best criteria combination to discriminate leaves between crops and weeds in field condition. Through this study, we also propose a Python framework allowing the extraction of these criteria available at github.com.

## 6.2 Material and data

### 6.2.1 Experimental plot

Data is acquired at the site of INRAE in Montoldre (Allier, France, at 46°20′30.3″N 3°26′03.6″E) within the framework of the ANR Challenge RoSE in 2019, visible on Figure 6.2. The objective of the Challenge is to objectively compare the solutions proposed by participants [Avrin et al., 2020, Regnier et al., 2021]. Within this context, the challenge provides an evaluation plan to contestants and a set of experimental plots of bean and corn plants. In addition various natural weeds (yarrows, amaranth, geranium, plantago, ...) and sowed ones (mustards, goosefoots, mayweed and ryegrass) is managed to compare performances.
<figure><img src="image_145.png" /><figcaption>**Figure 6.2** – Aerial view of the experimentation plot located in Montoldre (now INRAE)</figcaption></figure>

### 6.2.2 Multispectral camera

The images were acquired with the Airphen (Hyphen, Avignon, France) six-band multispectral camera (Figure 6.3). This is a multispectral scientific camera developed by agronomists for agricultural applications. The camera has been configured using the 450/570/675/710/730/850 nm bands with a 10 nm FWHM. The focal length of each lens is 8 mm. The raw resolutions for each spectral band are 1280x960 px with 12 bit precision. Due to the conception of the camera, spectral images are not aligned in space. A registration method based on previous work for this camera, with a registration accuracy down to sub-pixel is used [Vayssade et al., 2020]. After the registration, all spectral bands are rescaled to 1200 × 800 px and concatenated to channel-wise where each dimension refers to a spectral band.

<figure><img src="image_146.png" /><figcaption>**Figure 6.3** – The Airphen six-band multispectral camera</figcaption></figure>

### 6.2.3 Image acquisition and annotation

From the presented experimental plots a set of images were acquired. The camera is attached in front of an hybrid autonomous tractor called “TREKTOR” launched by SITIA (Bouguenais, France) in 2019. The camera is setup to have a top-down view of crop rows, thus it is placed at the end of a 2 m pole in front of the platform allowing to remove visible part of the robot and at 1.8 m from the ground to see all three crop rows. Crops and weeds were between phenological state 3 and 4 which means they had between 2 and 6 leaves. The ground truth is manually defined on images by experts with polygons around each leaf boundaries. In addition, polygons contain a crop/weed classification label. These annotation are performed using the VIA annotation software [Dutta and Zisserman, 2019] and a total of 300 images of bean were annotated, 170 in June and 130 in October. This dataset is now available online<fnref n="1" /> and may help other studies working on leaves classification. In addition the crop rows for each image are also manually set. The total number of leaves is 90 538, with 31 524 for crops and 59 019 for weeds.

<footnote n="1">1. https://data.inrae.fr/dataset.xhtml?persistentId=doi:10.15454/JMKP9S&version=1.0</footnote>
### 6.3 Methodology

The objective of this paper is to extract a wide variety of features and evaluate their performances for classifying leaves into crops and weeds, in addition to defining the best subset of feature space. The methodology followed in this study is as follows :

1. In order to evaluate criteria that could be used at the leaf level to discriminate between crop and weed leaves in field condition, a literature review was conducted. These criteria are organized into shape, spatial and textural features and presented in the section 6.3.1.

2. Most of the features have parameters that impact the extraction and classification performances of these features. To solve this problem, we propose to optimize these parameters by using an hyperparameter optimisation algorithm to efficiently sample the parameter space, it is described in section 6.3.2.

3. The results of the literature review indicate that there is no consensus on which of these criteria are the best for discriminating between crop and weed (at plant or leaf scale). In order to determine the best combination of criteria to use in discriminating between crop and weed leaves in field condition, an algorithm to select the best subsets of features is used and described in section 6.3.3.

4. The performance of each method is evaluated on the ground truth dataset, using Area Under the Curve (AUC) as the performance metric. The best performing feature subset is presented to classify leaves into crops and weed. These subsets are extracted for each type of feature (locally) and among all feature (globally), this is available in the results section 6.4.

5. Finally a classification of the vegetation cover is proposed using the best performing features subsets and the results are discussed in section 6.4.3.

#### 6.3.1 Literature review

In computer vision, shape, spatial and textural features are commonly used for describing and classifying objects. In this section we will present a review of these features. To facilitate the extraction of features and their merging, they have been programmed using the OOP<fnref n="2" /> paradigm in Python. Each feature is thus named according to its type (shape / spatial / spectral) and internal properties. Thus the following subsections present each of them according to its type.

##### 6.3.1.1 Shape Criteria

Shape properties are based on a structural analysis of the contours of related regions, it mostly consists of detecting morphological and anatomical traits. The shape can be coded as a set of positions (pixels) that go around the boundary. From this description, various information can be used. This type of analysis is found particularly in granulometry [Hentschel and Page, 2003], geomatics [McGarigal, 1995], leaf classification [Cerutti et al., 2014] or even animal posture detection [Vayssade et al., 2019, Bonneau et al., 2021]. These criterion are the most intuitive, easy to implement, and unaffected by lighting. Each of the following paragraphs explain an extraction method as well as the resulting properties used for discrimination.

<footnote n="2">2. Object-Oriented Programming</footnote>
**Shape Ellipse** [Watcharabutsarakham et al., 2012] have mentioned that the major and minor axis can be used to classify leaves. These properties can be obtained through an elliptical Hough algorithm as described by the article. From the retrieved ellipse, few more information can be computed such as the distance between the ellipse center and the detected shape center. Three factors of eccentricity, the directrix, the angle and the focal distance. Some of these properties were also used by [Lottes et al., 2016] for the classification of plant covers.

**Shape Particle** Based on the work of [Hentschel and Page, 2003], these shape metrics have been used to characterise and classify particles. [Vayssade et al., 2019] used it to classify goats activities within pasture, [Bonneau et al., 2021] to estimate the sow posture with an accuracy near deep-learning methods (94% vs 96%). This set of 18 characteristics are based on a variety of ratios of the Feret diameter, area and perimeter. Some of them are known as : Equivalent Diameter, Eccentricity, Aspect Ratio, Form Factor, Roundness or Convexity. Others are unna-med features. The Feret diameter is defined as (1) the longest distance between two points on the contour, and (2) the shortest distance across the shape. Some of these properties are also used by [Lottes et al., 2016] and [Saha et al., 2017].

**Shape Solidity** The solidity factor is widely used for weed plant classification [Ahmed et al., 2012]. The solidity factor is defined as a ratio between the segmentation mask area and convex hull area. Thus, it is bound from 0 to 1. This factor indirectly computes the holes area of the regions. When applied to weed discrimination, this factor allows to discriminate some specific species like carrot from others. When carrots are segmented, the segmentation mask is sparse, filled by a large amount of holes, due to leaves structures. This does not occur for dicotyledon plants, and it is less impacted by other monocotyledons weeds. At leaf scale, such criteria provide information about the contour (leaf teeth) as the inner leaf mask should be filled. This approach has been used for species classification [Thyagharajan and Raji, 2019] and weed discrimination [Lottes et al., 2016].

**Shape Metrics** This set of 8 properties comes from the documentation of the shape metrics tools of the University of Connecticut <fnref n="3" />, these are used for shape classification of lands [Angel et al., 2010] and building [Pászto et al., 2015]. Few properties, such as Cohesion, Traversal, Range and Viable Interior was not implemented due to time complexity. The others 8 implemented metrics are named : Proximity, Exchange, Spin, Perimeter Index, Depth, Girth, Dispersion and Detour.

**Shape Fragstat** The book by [McGarigal, 1995] shows some of the features used at the landscape, class, or parcel scale for classification or shape similarity measures. The fragstats documentation shows a small set of criteria <fnref n="4" /> that can be exploited at the leaf scale such as : Perimeter-Area Ratio, Shape Index, Fractal Dimension Index, Linearity Index, Related Circumscribing Circle, Contiguity Index and the Perimeter Area Fractal Dimension. These features have been selected as shape criterion in our study because they reflects shape complexity for a given spatial scale and correspond to : Perimeter-area-ratio, Shape Index, Fractal Dimension Index, Linearity Index, Related Circumscribing Circle, Contiguity Index and Perimeter Area Fractal Dimension.

<footnote n="3">3. http://clear.uconn.edu/tools/Shape_Metrics/method.htm</footnote>

<footnote n="4">4. http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Shape%20Metrics/SHAPE%20METRICS.htm</footnote>
**Shape Skeletonize** The skeletonization is commonly used in shape recognition techniques. It transforms the input mask (i.e the shape) to a new structural representation of the shape [Erdem and Tari, 2010]. This representation captures part of internal hierarchy. As instance, the [Lottes et al., 2016] study use the skeleton length for weed discrimination. In this study more features are extracted from that structural representation. These properties are defined as : the number of *end_point*, the number of internal *branches*, the number of *intersections* and the total length of the skeleton (*pixel_count*).

**Shape Angle** The distribution of angles between points around the contours can be exploited to measure the similarities between shapes [Belongie et al., 2002]. This method was used to classify symbols (digits and letters). The procedure is defined as follow : for a point of the contour, we compute the relative angles of all the other points that are in a given radius. The radii and angles are stored in a “shape matrix” of dimensions : *radial bins × angle bins*. Each element (i, j) of the matrix contains a counter which corresponds, for a given point, to the number of points which are in this bin of radius *i* and in the bin of angle *j*. Finally, the procedure is repeated for each point of the contour to feed the “shape matrix”. This one is used as a discriminating variable.

**Shape Hu Moments** The Hu Moments have been proposed long time ago by [Hu, 1962] and are part of structural analysis. The contours of the related regions are transformed into “moments”, which are then transformed into Hu Moments. The importance of this extraction comes from its translation, rotation, and scaling invariance. A similar element, e.g. two different bean leaves at various positions, rotation or scale will then be represented identically. This assertion remains true if intra-species and intra-individual variability are not taken into account. Hu Moment have been used for leaf classification challenge [Rhouma et al., 2017].

### 6.3.1.2 Spatial Criteria

The spatial properties are based on the extraction of global information related to the agricultural plot. These extraction methods are based on a distance which is calculated between the extracted global information and the centroid of the related region. The distribution of the values is strongly correlated to the vegetation density, allowing to discriminate the weeds, mostly between rows.

**Spatial Row** The distance between the plant and the crop line is one of the most used spatial descriptors when such prior knowledge is observable [Louargant et al., 2018, Koot, 2014, Bah et al., 2018b]. The number of techniques to extract such rows within an agronomic image is large, most of them are based on Hough Transform applied to line detection, as [Louargant et al., 2018]. Recent studies in deep-learning shows that CNN can be used to detect crop row in agricultural acquisition [Bah et al., 2019]. As the objective of this section is to evaluate the impact of the distance between plants and crop row (not the row detection performance), this study uses the crop rows from the ground truth.
**Spatial Blob** In most cases, plants can be classified into two types : monocotyledons and dicotyledons. At an early stage, when the plant emerges from the ground, its appearance is totally different according to this initial type. Monocotyledons usually looks like a filament, while dicotyledons look like an oval shape when seen from above. With this information, a specific shape detection algorithm can be implemented. Such as Blob detection, which aims to detect groups of connected pixels in an image that share a common property (shape, color, ...). The objective of Blob detection is to identify and mark these regions. These marks are used to calculate the Euclidean distance between the centroid of leaf and the nearest blob. Thus, the distribution of values is strongly correlated to the density of vegetation. This is an hypothesis we want to test and we have not seen any article using this type of approach.

**Spatial Corner** The detection of points of interest can be used to discriminate weeds. This is what is proposed in the article of [Xu et al., 2020], using an algorithm of detection of corners in the image. From these points of interest, the distance and the angle between the center of the plant and the closest points of interest are used as discriminating properties. In reality, this detection finds the leaves apexes. Here only the distance between the centroid and the nearest points of interest is computed since the angle is already present in Shape Ellipse properties.

### 6.3.1.3 Textural, Color and Spectral Criteria

Texture properties are vast and include many things, like image transformations (Fourier, Wavelets, color spaces, ...), color, spectral properties or histograms to characterise observed surfaces. The article [Mekhalfa and Yacef, 2021], have tested texture and color properties on RGB images of soybeans. The proposed color characteristics, are the means and standard deviations of the image in two color spaces (RGB and HSV). In terms of texture properties, i.e. an analysis of the spatial distributions of the colors, the paper relies on GLCM (Gray-Level Cooccurrence Matrix), Haralick and LBP (Local Binary Pattern) features. According to the study, the results are based on these properties and gave classification rates between 89.43% and 96.17%. Here, we have extended this work by adding their location to the extracted properties. Thus, all the extractions defined hereafter are automatically done on each individual bands, on NDVI and on the standard deviation (std) between spectral bands (which is also an image).

**Spectral Signature** The spectral signature is defined as the value (at the center of the shape and at the center of its bounding box) of all spectral bands, the NDVI value and the standard deviation between spectral bands.

**Spectral Stats** As shown by [Lottes et al., 2016], statistical properties of the underlying bounding box (that contains the shape) can be computed and used for weed discrimination. According to this article, a set of nine properties for each spectral bands is extracted separately as well as a vegetation index (NDVI) and the standard deviation (std). These properties are based on color moments (Table 6.1).

<table><caption>**Table 6.1** – Synthesis of colors moments</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>min</td><td>max</td><td>range</td></tr><tr><td>mean</td><td>std</td><td>median</td></tr><tr><td>skewness</td><td>kurtosis</td><td>entropy</td></tr></tbody></table>
**Spectral HoG** The HoG descriptors for Histogram of Oriented Gradient was used by [Saha et al., 2017] for weed discrimination. It is often used in computer vision, both for object detection and texture discrimination because of its properties of geometric invariance. The HoG descriptors are based on gradient intensity and edge directions. From this information, a histogram is computed and used as descriptor. Applied to a leaf, these properties allow to characterise leaf veins [Tsolakidis et al., 2014]. This features are generally robust to illumination changes and color differences due to the leaf maturity.

**Spectral Hu Moment** Moments and Hu Moments are generally applied to segmentation masks or contours. It retrieves shape properties, such as perimeter, centroid, etc. Computing these features on spectral bands and NDVI could include more discriminating information. This is an assumption that we wish to verify and there is no anterior of such use.

**Spectral Gabor** Gabor filters have been widely used for texture analysis in mono-bands images. By extension, [Jain and Healey, 1998] have proposed to use properties based on the results of Gabor filters applied to multispectral texture classification. They use a set of symmetric and circular Gabor filters, with three octave scales and four orientations to compute properties, such as averages and energy. These filters have been used in the field of weed discrimination by [Ishak et al., 2009]. In the current study, a symmetric and circular Gabor filters are applied with a single (optimized) octave scale and four orientations ([0, 40, 90, 158]) as defined by [Jain and Healey, 1998]. However, extracted properties from such filters are based on color moments (Table 6.1) on each Gabor orientation filter.

**Spectral LBP** The LBP descriptor, for Local Binary Pattern is a feature that encodes texture information. The general principle of LBP is to compare the level of a pixel luminance with the levels of its neighbours (0 if inferior, and 1 else). Depending on the neighbour position, a weight is applied to the positive value ([1, 2, 4, 8, 16, 32, 64, 128]), the LBP value is then the sum of weighted positives values. The last feature, is the histogram of the leaf texture for each spectral bands and NDVI. This gives an account of information relating to regular patterns in the image, in other words texture.

**Spectral CSLBP** The CS-LBP operator for Center-Symmetric Local Binary Pattern is an extension of the previously defined LBP operator [Banerji et al., 2012]. For each pixel, the absolute and symmetric difference of the neighborhood is computed, if the value is higher than a fixed threshold $T = 0.1$ the element takes the value 1 (0 otherwise). The set of results is coded as LBP through a weighted sum depending on the neighbour location. These features are considered to be more robust against rotation than standard LBP. This is an assumption that we wish to verify.

**Spectral OCLBP** The OC-LBP (Opposite Color LBP) operator is an extension of the LBP operation. [Banerji et al., 2012] propose an LBP extension to the colorimetric domain by considering the inter channel relationship instead of local and spatial relationship of a pixel. For each color pair $(u, v)$, $u > v$ the inter-channel properties are considered by $C = I(x, y, u) - I(x, y, v))$. With $I$ the multispectral image, $x, y$ the pixel position and $u, v$ the color channels to compare. The obtained image is defined in the same way as LBP : a weighted sum of positives values. These properties have been used by [Waghmare et al., 2016] for grape plant diseases detection.
**Spectral GLCM** This is a square matrix of dimension $N_g$, where $N_g$ is the number of gray levels in the image. The [i, j] element of the matrix is generated by counting the number of times a pixel of value i is adjacent to a pixel of value j, and then dividing the entire matrix by the total number of comparisons made. Each entry is thus considered as the probability that a pixel of value i is adjacent to a pixel of value j denoted $p(i, j)$, defining for us the co-occurrence matrix $G$. This have been used by [Bakhshipour et al., 2017] for weed discrimination.

**Spectral Haralick** The basis of these features is the GLCM (grayscale co-occurrence matrix). From that, Haralick have defined 14 texture descriptors that can be computed. This estimator was tested for canopy discrimination by [Wong et al., 2013].

**Spectral Polyfit** We introduce a new type of feature based on polynomial fitting. Multispectral images may lack information due to the selection of the spectral domain (The Airphen camera uses 6 spectral bands). To recover the original spectrum distribution across the surface, a simple solution that can be used is curve fitting [Brown, 2006], by using the value of observed pixel and their distances to the center of the shape. The result of the fitting is the polynomial parameter, which can be used as a discriminant variable.

**Spectral Shi Tomasi** The Shi-Tomasi feather is a kind of key point detector using the eigenvalue of the second moment of the color matrix [Harris et al., 1988]. On these values, a threshold is applied to determine if a pixel can be a corner or not. This algorithm has an interesting step, it transforms the input color image space into a gray scale color space, through this eigenvalue decomposition. Here, instead of finding the local maxima of this transformation, which are the key points. We propose to use a small histogram to describe the observed surface, so the values of corresponding to a corner accumulate in a specific bins of the histogram, while flat surfaces accumulate on other bins, etc. This transformation is applied on the first 3 bands of our data (450 nm, 570 nm, 650 nm) which are nearest to a RGB color space.

### 6.3.2 Optimization

Some features have parameters, like the number of bins in case of histogram analysis, the *inner* and *outer* radius of “shape angles” and so on. Some parameters are discrete numbers, but others are in a continuous space, like four parameter of the “texture gabor” which requires to set them before using it (sigma, lm, gm, ps). Setting these parameters has a critical impact on classification performances, they can be set empirically. But to optimize the performances of extracted features, these parameters must be optimized properly and automatically. This is known as an hyperparameter optimization problem.

There are many ways to sample the parameter space. One approach is grid search, where the parameter space is divided into a number of discrete regions (or grids) and a search is performed in each point of the grid [Lawler and Wood, 1966]. Another approach is random search, where a random point is chosen in the parameter space and the performance if evaluated [Bergstra and Bengio, 2012]. However, they can also lead to suboptimal solutions, because both methods are brute force approaches that scale poorly with the hyperparameters number. Therefore, it is important to choose an appropriate sampling method for the problem at hand. To achieve this optimization, more complex algorithms must be used to sample the parameter space efficiently.
We propose to use an asynchronous early stop method based on the continuous halving algorithm to solve this challenge, this is the ASHA algorithm (Asynchronous Successive Halving). The ASHA algorithm is a randomized, asynchronous and incremental algorithm for sampling the parameter space. It is based on the continuous halving algorithm, which is a greedy algorithm that divides the parameter space into two halves and selects the best half at each step. The ASHA algorithm uses this approach to iteratively divide the parameter space into smaller subsets. At each step, it selects the best subset based on a performance metric and stops when the subset is small enough that the performance is no longer improving. This approach is more efficient than grid search or random search, because it can find good solutions without exploring the entire parameter space. The asynchronous part of this algorithm allows it to run on multiple processors, which can speed up the search [Hertel et al., 2020].

In addition, to sample the parameter space correctly, each parameter has been manually bounded to limit the search space. The ASHA algorithm has been implemented in the Sherpa software package [Hertel et al., 2020] in Python. Finally, the optimization is done through 10000 epochs. At each epoch a 2-fold cross-validation is used for training and validating a decision tree classifier. This should enforce the learning of more stable properties and avoid overfitting, because the first fold contains the majority of data obtained in June, while the second fold contain the data obtained in October. The decision tree classifier is chosen to use simple decision criteria, rather than complex model, to extract most relevant features. The performance of the extracted features is therefore the performance of the underlying learned classifier (average of the cross-validation procedure). The Figure 6.4 synthesizes the optimization procedure :

<figure><img src="image_147.png" /><figcaption>**Figure 6.4** – Hyperparameter optimization procedure</figcaption></figure>

As an example, Table 6.2, show 5 epochs for the optimization of spectral Gabor features (parameters : sigma, lm, gm and ps), the objective value is the value to be maximized, it corresponds to the performances of classification model evaluated by the *Area Under the Curve* metric [Fawcett, 2006]. After the optimization process, all features are extracted.

<table><caption>**Table 6.2** – Five epochs of the optimization procedure through ASHA algorithm. The objective column show the model performance</caption><thead><tr><th>Epoch</th><th>sigma</th><th>lm</th><th>gm</th><th>ps</th><th>Objective</th></tr></thead><tbody><tr><td>1616</td><td>1.865</td><td>5.720</td><td>0.409</td><td>0.403</td><td>0.621</td></tr><tr><td>1617</td><td>1.838</td><td>5.732</td><td>0.363</td><td>0.421</td><td>0.604</td></tr><tr><td>1618</td><td>1.894</td><td>4.960</td><td>0.431</td><td>0.335</td><td>0.556</td></tr><tr><td>1619</td><td>3.363</td><td>8.253</td><td>0.083</td><td>0.690</td><td>0.465</td></tr><tr><td>1620</td><td>7.705</td><td>4.655</td><td>0.986</td><td>0.090</td><td>0.517</td></tr></tbody></table>
### 6.3.3 Feature Selection

Now that all the features are optimized and extracted, the main question of this article can be explored. What are the best features to classify leaves into crops and weeds. Since the number of features is huge (3545), as shown in Table 6.8, a manual estimation of each feature performances cannot be established properly. To answer this question, researchers in this field typically use a metric to estimate how well a feature fits the problem (correlation coefficient or root mean squared error), and then the best features are used. However, a better combination might exist, e.g., merging two high rank criteria might be less accurate than merging a high rank with a low rank. As an example, row distance and a texture characteristic, the efficiency may be greater than that of two texture criteria, because row distance is likely complementary to texture property. This is an assumption we wish to verify.

Instead of evaluating each feature through a metric, a data mining technique is proposed, such as "Sequential Forward Selection" (SFS) [Jain and Zongker, 1997]. It is an algorithm that analyzes each feature to check whether it is going to be useful for the classification task. To put it in another way the motivation behind this algorithm is to automatically select the subset of features most relevant to the problem. Thus, it reduces the initial *d*-dimensional feature space to a *k*-dimensional feature subspace, where *k* < *d*. The figure 6.5 shows how the algorithm works:

<figure><img src="image_148.png" /><figcaption>Figure 6.5 – Sequential Feature Selection Algorithm</figcaption></figure>

## 6.4 Results and Discussion

### 6.4.1 Feature optimisation

Among the 25 types of properties, only 12 of them have hyperparameters that need to be optimized. For each of them the parameters and the underlying classifier performances have been extracted in the below Table 6.3 (for spatial and shape feature) and Table 6.4 (for spectral features). To show the importance of the feature optimization process, we added the column **Worst**, which is the least discriminant parameter.
<table><caption>**Table 6.3** – Hyperparameter optimization through ASHA algorithm, the objective column shows the performance of the underlying classifier</caption><thead><tr><th>Type</th><th>Parameters</th><th>Values</th><th>Objective</th><th>Worst</th></tr></thead><tbody><tr><td rowspan="4">Shape Angles</td><td>binsRadius</td><td>2</td><td rowspan="4">62.81</td><td rowspan="4">50.00</td></tr><tr><td>bins_theta</td><td>22</td></tr><tr><td>radius_inner</td><td>0.9827</td></tr><tr><td>radius_outer</td><td>9.6870</td></tr><tr><td rowspan="4">Spatial Corner</td><td>blockSize</td><td>2</td><td rowspan="4">59.52</td><td rowspan="4">51.49</td></tr><tr><td>k</td><td>0.0995</td></tr><tr><td>ksize</td><td>1</td></tr><tr><td>threshold</td><td>0.09926</td></tr><tr><td rowspan="4">Spatial Blob</td><td>area</td><td>0</td><td rowspan="4">52.00</td><td rowspan="4">51.23</td></tr><tr><td>circularity</td><td>0</td></tr><tr><td>inertia</td><td>0</td></tr><tr><td>convexity</td><td>1</td></tr></tbody></table>

<table><caption>**Table 6.4** – Feature optimization through ASHA algorithm. To show the importance of the feature optimisation process, we added the column **Worst**, which is the least discriminant parameters. The column **Objective** shows the best performance of the model.</caption><thead><tr><th>Type</th><th>Parameters</th><th>Values</th><th>Objective</th><th>Worst</th></tr></thead><tbody><tr><td rowspan="2">Spectral DoG</td><td>sigma_a</td><td>17.7569</td><td rowspan="2">72.39</td><td rowspan="2">46.75</td></tr><tr><td>sigma_b</td><td>0.0475</td></tr><tr><td rowspan="4">Spectral Gabor</td><td>gm</td><td>37.51</td><td rowspan="4">71.68</td><td rowspan="4">49.95</td></tr><tr><td>lm</td><td>3.8094</td></tr><tr><td>ps</td><td>0.6573</td></tr><tr><td>sigma</td><td>8.7758</td></tr><tr><td rowspan="2">Spectral GLCM</td><td>normed</td><td>false</td><td rowspan="2">68.96</td><td rowspan="2">61.55</td></tr><tr><td>symmetric</td><td>true</td></tr><tr><td rowspan="2">Spectral Zernike</td><td>m</td><td>1</td><td rowspan="2">60.14</td><td rowspan="2">55.12</td></tr><tr><td>n</td><td>1</td></tr><tr><td>Spectral Haralick</td><td>size</td><td>128</td><td>69.66</td><td>52.99</td></tr><tr><td>Spectral HOG</td><td>size</td><td>256</td><td>54.96</td><td>51.08</td></tr><tr><td>Spectral LBP</td><td>nbins</td><td>256</td><td>69.17</td><td>53.16</td></tr><tr><td>Spectral Shi-Tomasi</td><td>nbins</td><td>14</td><td>67.23</td><td>51.84</td></tr><tr><td>Spectral CSLBP</td><td>nbins</td><td>22</td><td>69.53</td><td>55.33</td></tr></tbody></table>

Table 6.3 and Table 6.4, show a noticeable improvement of the feature performances through this optimization procedure. Few features are strongly enhanced, from 50% to 62.81% for “Shape Angles”, or from 46.75% to 72.39% for “Spectral DoG”, showing the necessity of hyperparameters optimisation. Irrelevant features start to appear here, for features that have a hyperparameter, this is the case for “Spatial Corner”, “Spatial Blob” or “Spectral HOG” that does not exceed 60% of AUC with a decision tree classifier.

Once the hyperparameters optimization is performed, all features (with or without hyper-parameters) can be extracted. Thus, the next step is to use the feature selection algorithm to select the best criteria from all extracted features. In the next section, we will show the results of the feature selection algorithm.
### 6.4.2 Feature selection by type

For this study, the SFS algorithm is applied independently on each feature type (locally), and on all extracted features (globally). This allows us to check which features are the most relevant for their categories as well as globally. For the SFS algorithm, the number of features included in the subspace is considered a stopping criterion. Here, the stopping criterion have been set to *k* = 10 for each feature type and *k* = 20 for all features. The result of the procedure is available on the next Table 6.5 for spatial feature, in Table 6.6 for shape feature and in Table 6.7 for spectral features. Finally figure 6.6 shows best global criteria. On all of these Tables, the best subspace size is highlighted in bold for each feature type, and the green cells indicate the best feature type to use for a specific feature space size.

<table><thead><tr><th>Feature Count</th><th>1</th></tr></thead><tbody><tr><td>Spatial Blob</td><td>51.80</td></tr><tr><td>Spatial Corner</td><td>56.59</td></tr><tr><td>Spatial Row</td><td>66.90</td></tr></tbody></table>

**Table 6.5** – Best features for spatial features. Best results for a specific feature count are highlighted

Since these 3 feature types extract only one feature each, Table 6.5 shows only one performance for each feature type. But it also shows that among the 3 spatial feature “Spatial Row”, which is the distance of the leaf to the crop row, has the most relevant performance. This is expected, as it allows to discriminate inter-row weeds. However, the weed coverage in October were high and overlapped the crop row, reducing this feature performance. The worst is “Spatial Blob”, which may be due to a lack in the optimization procedure for this feature, or more likely to the lack of discriminating information in this analysis, as for “Spatial Corner” that offers moderate performances. The next Table 6.6 is for shape criteria.

<table><thead><tr><th>Feature Count</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>Full</th></tr></thead><tbody><tr><td>Shape Angles</td><td><strong>65.73</strong></td><td>65.72</td><td>65.73</td><td>65.72</td><td>65.72</td><td>65.73</td><td>65.72</td><td>65.73</td><td>65.73</td><td>65.72</td><td>68.38</td></tr><tr><td>Shape Ellipse</td><td>60.84</td><td>67.47</td><td>69.46</td><td><strong>69.76</strong></td><td>69.56</td><td>69.48</td><td>69.28</td><td>69.10</td><td>68.80</td><td>68.54</td><td>67.58</td></tr><tr><td>Shape Fragstats</td><td>67.85</td><td><strong>70.10</strong></td><td><strong>74.83</strong></td><td>73.62</td><td>73.10</td><td>72.55</td><td>—</td><td>—</td><td>—</td><td>—</td><td>72.55</td></tr><tr><td>Shape Hu-moment</td><td>56.87</td><td>59.79</td><td>63.30</td><td>64.87</td><td>65.52</td><td><strong>65.78</strong></td><td>65.74</td><td>—</td><td>—</td><td>—</td><td>65.74</td></tr><tr><td>Shape Metrics</td><td>57.01</td><td>62.74</td><td>66.88</td><td>67.81</td><td>68.52</td><td>68.76</td><td><strong>68.98</strong></td><td>68.89</td><td>—</td><td>—</td><td>68.89</td></tr><tr><td>Shape Particle</td><td>67.88</td><td>70.15</td><td>74.50</td><td><strong>76.08</strong></td><td>76.04</td><td>75.99</td><td>76.03</td><td>75.95</td><td>75.83</td><td>75.73</td><td>73.13</td></tr><tr><td>Shape Skeletonize</td><td><strong>68.81</strong></td><td>68.81</td><td>68.65</td><td>68.52</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>68.52</td></tr><tr><td>Shape Solidity</td><td><strong>59.04</strong></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>59.04</td></tr></tbody></table>

**Table 6.6** – Best classification features for each feature type through ASHA algorithm. Each line presents the results for a specific feature type. The feature count is the number of feature(s) combined for a specific feature type. The presented results are the optimal combinations for each feature counts. Column “full” is the result obtained with all the features for a specific feature type. Best results for a type of features are in bold, best results for a specific feature count are highlighted.

As shown in Table 6.6, “Shape Solidity” that have been found relevant at plant scale [Ahmed et al., 2012] is less relevant at leaf scale. The first feature selected for “Shape Skeletonize” is the pixel-count. It refers to the length of the skeleton, which is the most discriminating property among the others of this feature type. Among all shape feature, the “Shape Particle” is the most
important one, the four best features are respectively “inv-area”, “m1”, “m2”, and “form-factor” with a score of 76.08%. While the second best shape property is “Shape Fragstats” that is used to classify land plot, its 3 best features are respectively the ‘perimeter-area-ratio’, the ‘related-circumscribing-circle’, and the ‘perimeter-area-fractal-dimension’ with a score of 74.83%. The four first features selected for the “Shape Ellipse” seem quite good, such as “center-distance”, “minor-axis”, “angle”, and “eccentricity 2” with an accuracy score of 69.76%. The other ones seem to been less relevant for crop and weed discrimination at leaf scale. Finally, Table 6.9 in the supplementary material, shows the name of the ten top best shapes features.

<table><thead><tr><th>Feature</th><th>Count</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>Full</th></tr></thead><tbody><tr><td>Spectral LBP</td><td></td><td>69.81</td><td>70.71</td><td><strong>70.72</strong></td><td>70.71</td><td>70.71</td><td>70.71</td><td>70.71</td><td>70.71</td><td>70.71</td><td>70.71</td><td>72.61</td></tr><tr><td>Spectral OCLBP</td><td></td><td>61.69</td><td>65.09</td><td>68.26</td><td>69.34</td><td>69.94</td><td>70.34</td><td>70.47</td><td>70.66</td><td>70.56</td><td><strong>70.60</strong></td><td>73.03</td></tr><tr><td>Spectral CSLBP</td><td></td><td><strong>70.36</strong></td><td>70.36</td><td>70.36</td><td>70.36</td><td>70.36</td><td>70.36</td><td>70.36</td><td>70.36</td><td>70.36</td><td>70.36</td><td>68.12</td></tr><tr><td>Spectral Zernike</td><td></td><td>57.26</td><td>68.67</td><td>72.58</td><td>73.42</td><td>74.44</td><td>75.42</td><td>75.80</td><td>76.19</td><td>76.37</td><td><strong>76.60</strong></td><td>76.91</td></tr><tr><td>Spectral Shi-Tomasi</td><td></td><td>71.11</td><td>67.74</td><td>70.02</td><td>70.29</td><td>71.25</td><td>72.01</td><td>72.59</td><td>72.84</td><td>72.94</td><td><strong>73.03</strong></td><td>72.00</td></tr><tr><td>Spectral Signature</td><td></td><td>66.76</td><td>67.37</td><td>74.40</td><td>75.79</td><td>77.25</td><td>78.18</td><td>79.23</td><td>79.35</td><td><strong>79.58</strong></td><td>79.53</td><td>63.65</td></tr><tr><td>Spectral Hu Moment</td><td></td><td>68.29</td><td>68.71</td><td>70.48</td><td>74.21</td><td>81.70</td><td>83.09</td><td>83.15</td><td>83.05</td><td><strong>83.22</strong></td><td>83.19</td><td>71.29</td></tr><tr><td>Spectral Stats</td><td></td><td>68.97</td><td>72.58</td><td>79.09</td><td>81.61</td><td>83.36</td><td>84.37</td><td>84.98</td><td>85.65</td><td>86.29</td><td><strong>86.54</strong></td><td>74.58</td></tr><tr><td>Spectral Polyfit</td><td></td><td>59.57</td><td>62.80</td><td>66.01</td><td>66.94</td><td>67.70</td><td>68.40</td><td>68.97</td><td>69.54</td><td>70.48</td><td><strong>70.63</strong></td><td>70.20</td></tr><tr><td>Spectral HOG</td><td></td><td>55.99</td><td>64.73</td><td>69.83</td><td>72.01</td><td>73.93</td><td>74.83</td><td>75.44</td><td>75.95</td><td>76.37</td><td><strong>76.67</strong></td><td>74.40</td></tr><tr><td>Spectral DoG</td><td></td><td>70.54</td><td>74.04</td><td>78.36</td><td>80.76</td><td>82.54</td><td>83.53</td><td>84.64</td><td>85.42</td><td>86.03</td><td><strong>86.33</strong></td><td>75.76</td></tr><tr><td>Spectral Gabor</td><td></td><td>69.50</td><td>74.48</td><td>80.07</td><td>82.23</td><td>83.53</td><td>84.52</td><td>85.14</td><td>85.60</td><td>86.15</td><td><strong>86.40</strong></td><td>74.92</td></tr><tr><td>Spectral GLCM</td><td></td><td>59.57</td><td>66.22</td><td>73.84</td><td>76.87</td><td>78.93</td><td>80.12</td><td>80.97</td><td>81.61</td><td>82.13</td><td><strong>82.56</strong></td><td>79.52</td></tr><tr><td>Spectral Haralick</td><td></td><td>63.80</td><td>75.85</td><td>78.89</td><td>81.16</td><td>82.56</td><td>83.69</td><td>84.62</td><td>85.01</td><td>85.33</td><td><strong>85.51</strong></td><td>81.30</td></tr></tbody></table>

**Table 6.7** – Best features individually for each feature type for classification through ASHA algorithm. Each line presents the results for a specific feature type. The feature count is the number of feature(s) combined for a specific feature type. The presented results are the optimal combinations for each feature counts. Column “full” is the result obtained with all the features for a specific feature type. Best results for a type of features are in bold, best results for a specific feature count are highlighted.

As shown in Table 6.7, the best features are the ones related to the texture and color of leaves. If a single feature should be used then, the best feature is the “Spectral Shi-Tomasi”, which is related to the roughness of the texture, since it is the first step of corner detection, however the remaining features are less discriminant than the others feature types. Among these spectral properties, the LBP, OCLBP and CSLBP were expected to produce better performances (they max out at around 70%), this means that the histogram construction is not producing relevant features. The same remarks can be established for the Zernike features even if their performances are slightly better (around 76.60%).

“Spectral Signature” is commonly used for vegetation classification, it uses the spectral value of the center of the shape (or window) as a criteria. The performances of this feature are better than the features previously discussed but max out at around 80%. The next feature, which is the “Spectral Hu Moment”, is a transformation of the spectral value of the center of the shape through Hu Moment, this transformation shows a significant increase in the performances of such features, from 79.53% for “Spectral Signature” to 83.19% for “Spectral Hu-Moment”, in addition this transformation seems more interesting on spectral properties than shape properties as seen in Table 6.6.
Features “Spectral Polyfit” and “Spectral Stats” are relatively close by definition but show a significant difference on their performances. “Spectral Stats” extracts statistical properties (Table 6.1) of the underlying sub-texture which contains the leaf and the soil around it, while “Spectral Polyfit” uses the spectral information of the sub-texture to reconstruct the spectral distribution of the leaf surface, then the polynomial values are used as discriminating information. The initial assumption was that “Spectral Polyfit” (which is an attempt to reconstruct the initial spectrum distribution from the discrete values) should have been more relevant than “Spectral Stats”.

“Spectral HOG” for Histogram of Oriented Gradient, is generally used to retrieve gradient information. In case of leaf analysis, this transformation detects leaf veins [Larese et al., 2014], which is a criterion to discriminate leaf species. Since this analysis doesn’t produce good enough performances, it is possible that the spatial resolution is too low for this kind of methodology.

The remaining “DoG”, “Gabor”, “GLCM” and “Haralick” features offer good performances as expected. We can notice that “Spectral Gabor” is the best of the four when ten features are used, while “DoG” and “GLCM” show slightly less accurate results. Finally, Table 6.10 in the supplementary material, shows the name of those ten top best spectral features.

### 6.4.3 Best of all features

<figure><img src="image_149.png" /><figcaption>**Figure 6.6** – Best features among all extracted feature for classification through ASHA algorithm. Best overall result is highlighted.</figcaption></figure>

As shown previously (Table 6.7), the best starting feature is the second bin of “Spectral Shi-Tomasi” histogram, thus it is selected firstly and its performance level is found in the first point of Figure 6.6 (lower-left). The second selected feature is “Spatial Row”, it is expected that this feature should be on the top as it is a complementary with the color/texture features. At this
stage, we can notice that this combination is less accurate 74.98 % than two texture Haralick feature 75.85 %. The next three features are “Spectral Stats 2_min”, “Spectral Stats 5_max” and “Spectral Gabor 3_0 mean”. Which is the minimum spectral value of the underlying sub-texture of the 650 nm spectral band for the first feature. In the same way, it is the maximum for the 850 nm spectral band. The Gabor filter is applied on all spectral bands, the discriminating criterion is one from the 450 nm. This leads to a performance of 85.95 %. The next feature is “Shape Metrics Detour”, with a performance of 87.41 % which beats any previous classification model.

To conclude on this results, the best classification model includes 1 spatial feature, 6 shape features and 13 spectral features. In the previous Figure 6.6, only 20 features are shown, if more features are extracted, the performances begin to decrease (21 features = 91.16 %), this means that the performance reaches a vertical asymptote. To show this a model have been learnt with all features (3545), the cross-validation score is 79.97 %.

### 6.4.4 Visual results

Figure 6.7 shows two rows with field beans plants. Weeds are located between these two rows and on each of them. Most of plant and weed leaves are correctly classified (plant leaves are colored green, and weed leaves are colored blue). However, this image suggests a problem remains for the smallest leaves that are in the center of the plant. These plants leaves, colored in red, are classified as weed leaves. This classification is induced by shape features, which are less relevant for this category of leaves than for the others.

<figure><img src="image_150.png" /><figcaption>**Figure 6.7** – This figure shows the classification of leaves, from an image taken in June 2019. The blue color corresponds to well classified weeds. The green color corresponds to well classified crops. While the purple and red color correspond respectively to weeds and crops poorly classified.</figcaption></figure>
Figure 6.8 focuses on specific parts of plant rows. Weeds on figure 6.8.a are mostly monocots, whereas those visible on figure 6.8.b are mostly dicots. These two figures illustrates the quality of classification, which is equivalent for the two weed groups. Figure 6.8.c presents a dense foliage, with highly developed wild mustard. On this figure, the major part of mustard leaves is correctly classified. The classification of field bean leaves depends on their visibility. Several partially occluded leaves are classified as weed leaves (in red). Conversely, figure 6.8.d presents sparse foliage, where many bean plants have not emerged. In that case, the small bean leaves are also classified as weed leaves.

<figure><img src="image_151.png" /></figure>

**Figure 6.8** – This figure shows the classification of leaves, from parts of four images taken in June 2019. The blue color corresponds to well classified weeds. The green color corresponds to well classified crops. While the purple and red color correspond respectively to weeds and crops poorly classified. Here we see that leaves between crop rows are correctly classified, as are monocots. On the other hand, the emerging leaves of crops and some weed type leaves of dicotyledons are incorrect.

All of these images illustrate the difficulty in classifying leaves which, being in the row of crops, have particular shapes, because they are partially obscured or because they are not sufficiently developed.
### 6.5 Conclusion

The properties we evaluated were identified from a thorough literature review. We distinguished spatial properties, based on notions of distances at the image scale, shape properties allowing to characterize the morphological features of the individuals. And, finally, properties of textures, colors and spectra were proposed to extract information related to the composition and internal structures of leaves.

From this survey, it came out that some methods have parameters that influence the discrimination potential of these properties. Therefore, an algorithm was proposed to optimize these hyperparameters before the extraction and evaluation of these features. It has been shown that the optimization of these hyperparameters plays a crucial role in the performance of these criteria. For example, “Shape Angles” properties went from 50% to 62.81% of discrimination potential after optimization, as well as the “Spectral DoG” properties from 46.75% to 72.39%, etc.

After the optimization of these properties, the set of properties was extracted from the dataset for evaluation. In this evaluation, it was shown the importance of selecting a smaller set of properties, the overall score is better with 20 features than all the features (3545). It is also important to notice the contribution of each feature decreases as the number of features increases, resulting in a less interesting “performance over computing time” ratio. For each property type, the top 10 properties were defined. Then the top 20 among all properties combined were defined. These results show a crop/weed discrimination performance up to 91 percent with 20 properties, while a performance of 79.97% is observed using the 3545 properties. Indeed, the usual classification techniques do not seem to be able to handle large dimensions of properties. But here, the performances of classifications proposed by smaller sets are superior to the use of all the criteria. Thus, the definition of this minimal set also allows to decrease the computation time.
## 6.6 Further research

Today, deep learning may be the only competitive alternative that could improve the results of each feature type. Following the approach that has been undertaken in a previous article [Vayssade et al., 2021], an alternative based on function approximators is possible. We could then define the best shape properties from a fixed representation of the contour. The same is true for texture and spatial properties. In reality, a UNet model already synthesizes the set of texture, color, spectrum and a kind of shape features if used in this crop/weed discrimination framework. Moreover, the performance of the extraction of these features does not allow a real time use, which reinforces the need to study deep learning.

In this study, feature selection is only based on classification performances of features. Computation speed is neither estimated nor taken into account. However, depending on the application, a compromise between computation speed and classification performance could be used. The selection of the best features for real-time weed detection and elimination is of particular interest. Figure 6.6 shows that the major classification improvement comes from the 10 first features (from 71 to 90%), the next 10 features improve the classification result of a bit more than 1%. This observation shows that a study focused on performance over computing time should be relevant. This is also a key step in developing on-board weed control systems in autonomous vehicles where computing power can be restricted and time matters the most.

Considering precision farming applications, this approach could be replicated to select the best features for a specific set of crop(s) and weeds. The objective could change from an application to another and should be linked to the classification accuracy requirements and the time available. A real-time application may focus on computing time with less consideration for classification performance. A UAV approach may allow more computing time as the result is already delayed from the acquisition. In both case, the image resolution may be an important setting : too high and the computing time may be too important, too low and the classification performance may drop too much. This resolution may also impact the feature contribution to the overall classification with spectral or shape feature requiring highest resolution than spatial features.