Reflectance Estimation from Multispectral Linescan Acquisitions under Varying Illumination—Application to Outdoor Weed Identification ^†

Abstract

To reduce the amount of herbicides used to eradicate weeds and ensure crop yields, precision spraying can effectively detect and locate weeds in the field thanks to imaging systems. Because weeds are visually similar to crops, color information is not sufficient for effectively detecting them. Multispectral cameras provide radiance images with a high spectral resolution, thus the ability to investigate vegetated surfaces in several narrow spectral bands. Spectral reflectance has to be estimated in order to make weed detection robust against illumination variation. However, this is a challenge when the image is assembled from successive frames that are acquired under varying illumination conditions. In this study, we present an original image formation model that considers illumination variation during radiance image acquisition with a linescan camera. From this model, we deduce a new reflectance estimation method that takes illumination at the frame level into account. We experimentally show that our method is more robust against illumination variation than state-of-the-art methods. We also show that the reflectance features based on our method are more discriminant for outdoor weed detection and identification.

1. Introduction

Nowadays, one of the biggest topics in precision farming is to increase yield production while reducing the quantity of chemicals. In order to optimize the application of herbicides in crop fields, less toxic and expensive weed control alternatives can be considered due to recent advances in imaging devices. During the last decade, sophisticated multispectral sensors have been manufactured and deployed in crop fields, leading to weed detection [1,2,3].

Multispectral cameras collect data over a wide spectral range and they provide the ability to investigate the spectral responses of soils and vegetated surfaces in narrow spectral bands. Two main categories of devices can be distinguished in multispectral image acquisition. “Snapshot” (multi-sensor or filter array-based) devices build the image from a single shot [4]. Although this technology provides multispectral images at a video frame rate, the few acquired channels and low spatial resolution may not be sufficient for fully exploring the vegetation spectral signatures. “Multishot” (tunable filter or illumination-based, push-broom, and spatio-spectral linescan) devices build the image from several and successive frame acquisitions [5,6,7]. Despite being restricted to still scenes, they provide images with a high spectral and spatial resolution. We use a multishot camera, called the “Snapscan” to acquire outdoor multispectral radiance images of plant parcels in a greenhouse under skylight [8]. From this radiance information, reflectance is estimated as an illumination-invariant spectral signature of each species. Several methods have been proposed for computing reflectance thanks to prior knowledge regarding cameras or illumination conditions [9,10,11]. In field conditions, the typical methods first estimate the illumination by including a reference device (a white diffuser or a color-checker chart) in the scene [2,12,13,14]. Subsequently, reflectance is estimated at each pixel p by the channel-wise division of the value of the radiance image at p by the pixel values that characterize the white diffuser or the color-checker white patch. In [15], an extension to the multispectral domain of four algorithms traditionally applied to RGB images is proposed for estimating the illumination. In [16], a Bragg-grating-based multispectral camera acquires outdoor radiance images and reflectance is estimated from two white diffusers, one along the image bottom border and another fully visible one. A coarse reflectance is computed first, and then rescaled using illumination-based scaling factors. Finally, the resulting reflectance image is normalized channel-wise by the average reflectance value that is computed over pixels of the second white diffuser that is present in the scene. In [17], a multispectral camera is used in conjunction with a skyward pointing spectrometer to estimate the reflectance from the acquired scene radiance.

These methods require additional devices and knowledge regarding the spectral sensitivity functions (SSFs) of the sensor filters. They also often assume constant incident illumination throughout a few seconds. However, in outdoor conditions, illumination may vary significantly during the successive frame acquisitions (scans) that last for several seconds.

In this paper, we propose a reflectance estimation method that is robust to illumination variations during the multispectral image acquisition that was performed by the Snapscan camera. In Section 2, we provide details regarding multispectral radiance image acquisition with this device and propose an original model for such image formation. In Section 3, this model is first used to study reflectance estimation under constant illumination. We also show that, when outdoor illumination varies during the frame acquisitions, the assumption about spatio-spectral correlation does not hold. Based on this model, we then propose a new method for estimating the scene reflectance from a multispectral radiance image acquired in uncontrolled and varying illumination conditions (see Section 4). Section 5 presents an experimental evaluation of the proposed reflectance estimation, and Section 6 shows the results of weed/crop segmentation using the estimated reflectance.

2. Multispectral Radiance Image Acquisition by Snapscan Camera

In this section, we first detail how the Snapscan camera achieves radiance measurement and frame acquisition. Subsequently, we explain how a multispectral image is obtained from the successively acquired frames. We propose an original multispectral image formation model that handles how illumination is associated to both the considered band and pixel because the radiance that is associated to a given spectral band at a pixel is measured in a frame that is acquired at a specific time. From this new model, we show that the spatio-spectral correlation assumptions do not hold when illumination varies during the frame acquisitions.

2.1. Radiance Measurement

The Snapscan is a multispectral camera manufactured by IMEC that embeds a single matrix sensor that is covered by a series of narrow stripes of Fabry-Perot integrated filters. It contains $B=192$ optical filters whose central wavelengths range from ${\lambda }^0=$ 475.1 nm to ${\lambda }^{B-1}=$ 901.7 nm with a variable center step (from 0.5 nm to 5 nm). Specifically, each filter of index $b\in [0,B-1]$ is associated with five adjacent rows of 2048 pixels that form a filter stripe, and it samples a band from the visible or near infra-red spectral domain according to its SSF $T^b\left(\lambda \right)$ with a full width at half maximum between 2 nm and 10 nm.

The Snapscan camera acquires a sequence of frames to provide a multispectral image. During frame acquisitions, the object and camera both remain static while the sensor moves and illumination may change. Therefore, the measurement of the radiance that is reflected by a given lambertian surface element s of the scene varies according to the frame acquisition time t, although s is projected at a fixed point q of the image plane. Let us denote, as $E_t\left(\lambda \right)$$\in [0,1]$, the relative spectral power distribution (RSPD) of the illumination at t and assume that it homogeneously illuminates all of the surface elements of the scene. The radiance that is reflected by s and refracted by the camera lens projects onto the image plane at q as a stimulus ${\mathcal{L}}_{t,q}\left(\lambda \right)$:

$${\mathcal{L}}_{t,q}\left(\lambda \right)=E_t\left(\lambda \right)·R_q\left(\lambda \right)·A_q\left(\lambda \right),$$

(1)

where $R_q\left(\lambda \right)$$\in [0,1]$ is the spectral reflectance of the surface element s that is observed by q, and $A_q\left(\lambda \right)$$\in [0,1]$ is the optical attenuation of the camera lens at q. All of these functions depend on the wavelength $\lambda$. The sensor moves forward on the image plane according to the direction perpendicular to the filter stripes (see Figure 1a). Between two successive frame acquisitions, it moves by a constant step $v=5$ (in pixels) that is equal to the number of rows in each stripe. Therefore, the radiance that is measured at q is filtered by a different Fabry–Perot filter of index $b_{t,q}$ ($0\le b_{t,q}<B$) at each acquisition time t. The radiance at q is fully sampled over N frame acquisitions, provided that each of them measures the radiance there, i.e., $N\ge B$. Let the coordinates of point q be ${(x_q,y_q)}_C$ in the camera 2D coordinate system $(O,x,y)$ whose origin O corresponds to the intersection between the optical axis and image plane. The unit vectors of x and y are given by the photo-sensitive element size (i.e., axis units match with pixels), and y is oriented opposite to the sensor movement. At a given point q, the filter index $b_{t,q}$ can then be expressed as:

$$b_{t,q}=t+⌊\frac{y_q-y_0}{v}⌋,$$

(2)

where $y_0$ is the coordinate along y of the first filter row at first acquisition time $t=0$. Note that the light stimulus ${\mathcal{L}}_{t,q}$ is only associated to a filter at a given point q when $t_q^0\le t<t_q^B$. The lower bound $t_q^0=⌊(y_0-y_q)/v⌋$ is the acquisition time at which the first optical filter of the sensor observes ${\mathcal{L}}_{t,q}$. The upper bound $t_q^B=⌊(y_0-y_q)/v⌋+B$ is the time at which all of the sensor filters have observed ${\mathcal{L}}_{t,q}$.

Besides, at a given time t, the coordinate $y_q$ of point q that is associated to a photo-sensitive element of the sensor satisfies:

$$y_0-t·v\le y_q<y_0+(B-t)·v$$

(3)

since $0\le b_{t,q}<B$. Given these restrictions, the radiance ${\mathcal{S}}_{t,q}$ that is then measured at q by the sensor at acquisition time t is expressed as:

$${\mathcal{S}}_{t,q}=Q\left(\tau {\int }_{\Omega }{\mathcal{L}}_{t,q}\left(\lambda \right)·T^{b_{t,q}}\left(\lambda \right)d\lambda \right),$$

(4)

where Q is the quantization function according to the camera bit depth, $\tau$ is the integration time of the frames, and $\Omega$ is the working spectral domain. Note that $\tau$ is set to the highest possible value that provides no saturated pixel.

2.2. Frame Acquisition

The radiance that is measured at q is stored by the camera as a pixel value $f_{t,q}={\mathcal{S}}_{t,q}$ in frame $f_t$ (see Figure 1b). We define the coordinate system $(O^{\prime },x^{\prime },y^{\prime })$ attached to the sensor, such that origin $O^{\prime }$ is the first (top-left) photo-sensitive element location, axis $y^{\prime }$ corresponds to y, and $x^{\prime }$ is parallel to x, in order to compute the coordinates of q relative to the frame. In this frame system, the coordinates of q are ${(x_q^{\prime },y_q^{\prime })}_F={(x_q,y_q-y_0+t·v)}_F$. Note that Equation (3) allows us to check that $0\le y_q^{\prime }<B·v$.

Conversely, any given pixel $p{(x_p^{\prime },y_p^{\prime })}_F$ of a frame $f_t$ is mapped to the coordinates in the camera coordinate system as:

$${(x_p,y_p)}_C={(x_p^{\prime },y_p^{\prime }+y_0-t·v)}_C,$$

(5)

at which the stimulus ${\mathcal{L}}_{t,p}\left(\lambda \right)$ of a surface element radiance is filtered by the filter of index $b_p=\lfloor y_p^{\prime }/v\rfloor$. From this point of view, each frame pixel value is, therefore, also expressed as:

$$f_{t,p}={\mathcal{S}}_{t,p}=Q\left(\tau {\int }_{\Omega }{\mathcal{L}}_{t,p}\left(\lambda \right)·T^{b_p}\left(\lambda \right)d\lambda \right).$$

(6)

Before the frame acquisitions, the Snapscan uses its internal shutter to acquire a dark frame $f_{dark}$ whose values are subtracted pixel-wise from the acquired frames. Therefore, we assume that the pixel value that is expressed by Equation (6) is free from thermal noise.

Let us also point out that, at two (e.g., successive) acquisition times $t_1$ and $t_2$, the sensor is at different locations. Therefore it acquires the values $f_{t_1,p}$ and $f_{t_2,p}$ from the stimuli ${\mathcal{L}}_{t_1,p}$ and ${\mathcal{L}}_{t_2,p}$ of two different surface elements at a given pixel p whose coordinate $y_p$ in the camera system is time-dependent (see Equation (5)). Besides, the stimuli ${\mathcal{L}}_{t_1,p}$ and ${\mathcal{L}}_{t_2,p}$ are filtered by the same filter whose index only depends on the pixel coordinate $y_p^{\prime }$ in the frame system. Equations (4) and (6) model the radiance that is measured at a given point in the image plane and stored at a given pixel of a frame, respectively. Both of the equations take account of illumination variation during the frame sequence acquisition, but differently take the sensor movement into account. Indeed, the filter index changes at a given point of the image plane during the frame acquisition (see Equation (4)), whereas the observed surface element changes at a given pixel in the successive frames (see Equation (6)).

2.3. Stripe Assembly

We now determine the first and last acquisition times of the frame sequence that is required to capture an object of interest whose projection points on the image plane are bounded along the y axis by $q_a{(x_{q_a},y_{q_a})}_C$ and $q_m{(x_{q_m},y_{q_m})}_C$, with $y_{q_a}>y_{q_m}$. Given the initial coordinate $y_0$ of the sensor along y, we can compute the first and last frame acquisition times $t_{q_a}^0$ and $t_{q_m}^{B-1}$, so that the measured radiances at the points between $q_a$ and $q_m$ are consecutively filtered by the B sensor filters (see the top part of Figure 2). The acquisition of the multispectral image from the frame sequence ${\left\{f_t\right\}}_{t=t_{q_a}^0}^{t_{q_m}^{B-1}}$ takes account of the spatial and spectral organizations of each frame. A frame $f_t$ is spatially organized as juxtaposed stripes of v adjacent pixel rows. A stripe $f_t^b$, $b=0,\cdots ,B-1$, of v adjacent pixel rows contains the spectral information of the scene radiance that is filtered according to the SSF $T^b\left(\lambda \right)$ of filter b centered at wavelength ${\lambda }^b$. All of the stripes that are associated with filter b in the acquired frames are stacked by the assembly function ⨁ to provide a stripe assembly defined as:

$${\left\{f_t^b\right\}}_{t=t_{q_a}^0}^{t_{q_m}^{B-1}}\stackrel{\mathrm{def}}{=}⨁\left({\left\{f_t\right\}}_{t=t_{q_a}^0}^{t_{q_m}^{B-1}},b\right)={[f_{t_{q_m}^{B-1}}^b,\cdots ,f_{t+1}^b,f_t^b,f_{t-1}^b,\cdots ,f_{t_{q_a^0}}^b]}^{⊺}.$$

(7)

The size of each stripe assembly is 2048 pixels in width and $N·v$ pixels in height, where $N=⌊\left(t_{q_m}^{B-1}-t_{q_a}^0\right)/\Delta ⌋+1$ is the number of acquired frames and $\Delta$ is the frame acquisition period.

To form the multispectral image ${\mathbf{I}}^{\left(B\right)}={\left\{I^b\right\}}_{b=0}^{B-1}$ of the object of interest, only the scene part that is common to all stripe assemblies is considered by the camera (see the bottom part of Figure 2). Specifically, the retained stripes in the b-th assembly are acquired between $t_{q_a}^b$ and $t_{q_m}^b$ to form each channel $I^b$:

$$I^b={\left\{f_t^b\right\}}_{t=t_{q_a}^b}^{t_{q_m}^b}.$$

(8)

The multispectral image ${\mathbf{I}}^{\left(B\right)}$ has its own coordinate system. For convenience, in the sequel, we denote a pixel as $p(x,y)$ in this system, since the camera and frame coordinate systems are not used any longer.

2.4. Formation Model of a Multispectral Image Acquired by Snapscan Camera

We can now infer an image formation model for multishot linescan cameras, such as the Snapscan. At any pixel p, the radiance value $I_p^b$ that is associated to a channel index $b\in [0,B-1]$ is acquired at $t=t_p^b$, with $t_{q_a}^b\le t_p^b\le t_{q_m}^b$ (see Equation (8)). It results from the light stimulus ${\mathcal{L}}_{t_p^b,p}$ that was filtered according to $T^b$ (whose index dependence upon p is dropped by stripe assembly step), and is therefore defined from Equations (1) and (6), as:

$$I_p^b=Q\left(\tau {\int }_{\Omega }{E}_{t_p^b}\left(\lambda \right)·R_p\left(\lambda \right)·A_p\left(\lambda \right)·T^b\left(\lambda \right)d\lambda \right).$$

(9)

The term $E_{t_p^b}\left(\lambda \right)$ shown in Equation (9) points out that illumination is associated to both a channel index and a pixel. These dependencies may weaken the spatio-spectral correlation assumptions of the measured scene radiance.

Spectral correlation relies on the assumption that the SSFs that are associated to adjacent spectral channels strongly overlap. Thus, radiance measures at a given pixel in these channels should be very similar (or correlated). Let us consider the radiance values in two channels $b_1$ and $b_2$ at a given pixel p. Even if the SSFs $T^{b_1}\left(\lambda \right)$ and $T^{b_2}\left(\lambda \right)$ strongly overlap (and are equal in the extreme case), the illumination conditions at $t_p^{b_1}$ and $t_p^{b_2}$ are different, hence $I_p^{b_1}\ne I_p^{b_2}$.

Spatial correlation relies on the assumption that the reflectance across locally close surface elements of a scene does (almost) not change. Thus, under the same illumination, the radiance measures at their associated pixels within a channel are correlated. Let us consider two pixels, $p_1(x_{p_1},y_{p_1})$ and $p_2(x_{p_2},y_{p_2})$, which observe surface elements of a scene with the same reflectance $R_{p_1}\left(\lambda \right)=R_{p_2}\left(\lambda \right)$ for all $\lambda \in \Omega$. If $|y_{p_1}-y_{p_2}|\ge v$, then the radiances at $p_1$ and $p_2$ are acquired at different times $t_{p_1}^b$ and $t_{p_2}^b$ associated to different illumination conditions $E_{t_{p1}^b}$ and $E_{t_{p2}^b}$, hence $I_{p_1}^b\ne I_{p_2}^b$.

Therefore, the spatio-spectral correlation assumption does not hold in the image formation model of the Snapscan camera when illumination varies.

3. Reflectance Estimation with a White Diffuser under Constant Illumination

This short section introduces how to estimate reflectance by a classical (white diffuser-based) method and how the result should be post-processed to ensure its consistency.

3.1. Reflectance Estimation

In order to estimate spectral reflectance from radiance images that were acquired under an illumination that is almost constant over time, one classically uses the image ${\mathbf{I}}^{\left(B\right)}\left[\mathrm{WD}\right]$ of a white diffuser acquired in full field beforehand and assumes that:

(i)

The illumination is spatially uniform and it does not vary during the frame acquisitions, thus $E_{t_p^b}\left(\lambda \right)=E\left(\lambda \right)$ for all $b\in [0,B-1]$ and $p\in {\mathbf{I}}^{\left(B\right)}$, and Equation (9) becomes:

$$I_p^b=\tau ·\sum _{l=0}^{B-1}E\left({\lambda }^l\right)·R_p\left({\lambda }^l\right)·A_p\left({\lambda }^l\right)·T^b\left({\lambda }^l\right).$$

(10)

Note that the quantization function Q is omitted here, since the different terms are considered as being already quantized.

(ii)

Each of the Fabry-Perot filters has an ideal SSF $T^b\left(\lambda \right)=\delta (\lambda -{\lambda }^b)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}\lambda ={\lambda }^b,\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$ such that Equation (10) becomes:

$$I_p^b=\tau ·E\left({\lambda }^b\right)·R_p\left({\lambda }^b\right)·A_p\left({\lambda }^b\right).$$

(11)

Reflectance is then derived for any pixel p that is associated to a spectral band centered at ${\lambda }^b$ as:

$$R_p\left({\lambda }^b\right)={\displaystyle \frac{I_p^b}{\tau ·E\left({\lambda }^b\right)·A_p\left({\lambda }^b\right)}}.$$

(12)

The white diffuser is supposed to be perfectly diffuse and reflect the incident light with a constant diffuse reflection factor ${\rho }_{\mathrm{wd}}$. Hence, for ${\mathbf{I}}^{\left(B\right)}\left[\mathrm{WD}\right]$, we can write:

$${\rho }_{\mathrm{wd}}={\displaystyle \frac{I_p^b\left[\mathrm{WD}\right]}{{\tau }_{\mathrm{wd}}·E\left({\lambda }^b\right)·A_p\left({\lambda }^b\right)}},$$

(13)

where ${\tau }_{\mathrm{wd}}$ is the frame integration time of ${\mathbf{I}}^{\left(B\right)}\left[\mathrm{WD}\right]$. Plugging Equation (13) into (12) yields the reflectance image that is estimated from a B-channel radiance image ${\mathbf{I}}^{\left(B\right)}$:

$${\widehat{R}}_p^b={\rho }_{\mathrm{wd}}·{\displaystyle \frac{I_p^b}{I_p^b\left[\mathrm{WD}\right]}}·{\displaystyle \frac{{\tau }_{\mathrm{wd}}}{\tau }}.$$

(14)

This reflectance estimation model implicitly compensates the vignetting effect, since the white diffuser and object (scene of interest) occupy the same (full) field of view. Accordingly, $I_p^b$ and $I_p^b\left[\mathrm{WD}\right]$ are affected by the same optical attenuation whose effect vanishes after division.

The estimated B-channel reflectance image ${\widehat{\mathbf{R}}}^{\left(B\right)}$ should then undergo two post-processing steps: spectral correction and negative value removal.

3.2. Spectral Correction

Each of the Snapscan Fabry–Perot filters is designed to sample a specific spectral band from the spectrum according to its SSF $T^b\left(\lambda \right)$. However, because of the SSFs and optical properties of some filters (angular dependence [18], high-energy harmonics), several spectral bands are redundant, which limits the accuracy of the spectral imaging system. This leads to redundancy in spectral bands and introduces spectral information bias. Therefore, the reflectance image with $B=192$ spectral channels is spectrally corrected and only $K=141$ channels are kept in practice.

The spectral correction of ${\widehat{\mathbf{R}}}^{\left(B\right)}$ provides a spectrally corrected K-channel reflectance image ${\widehat{\mathbf{R}}}^{\left(K\right)}$ that is expressed at each pixel p as:

$${\widehat{\mathbf{R}}}_p^{\left(K\right)}=\mathbf{M}·{\widehat{\mathbf{R}}}_p^{\left(B\right)},$$

(15)

where $\mathbf{M}$ is the sparse $K\times B$ correction matrix that is provided by the calibration file of our Snapscan camera. The linear combinations of the channel values of ${\widehat{\mathbf{R}}}^{\left(B\right)}$ according to Equation (15) are designed by the manufacturer to remove the redundant channels and attenuate second-order harmonics. This spectral correction provides new centers ${\left\{{\lambda }^k\right\}}_{k=0}^{K-1}$ for the bands (referred to as “virtual” bands by IMEC) that are associated to the image channels, but the spectral working domain $\Omega =[475.1\mathrm{nm},901.7\mathrm{nm}]$ is unchanged.

3.3. Negative Value Removal

The acquired radiance image contains negative values due to dark frame subtraction, when the value of a dark frame pixel is higher than the measured radiance at this pixel. This generally occurs in low-dynamics channels, where the central wavelengths are in the range [$475.1$ nm, $560.4$ nm] (before spectral correction). These negative values may lay on vegetation pixels and corrupt reflectance estimation at these pixels. Because we intend to classify vegetation pixels, this could lead to unexpected prediction errors. Negative values also occur—for even more pixels—in the spectrally-corrected reflectance image ${\widehat{\mathbf{R}}}^{\left(K\right)}$ (see Equation (15)), because the correction matrix M contains negative coefficients.

Negative values have no physical meaning and they must be discarded. Because our images mostly contain smooth textures (vegetation, reference panels, soil), we consider that, unlike radiance, reflectance values are highly correlated over close surface elements. Thus, we propose correcting negative values in image ${\widehat{\mathbf{R}}}^{\left(K\right)}$ by conditionally using a $3\times 3$ median filter, as:

$${\widehat{R}}_{\mathrm{ref},p}^k=\left\{\begin{array}{cc}\underset{3\times 3}{\mathrm{median}}\left\{{\widehat{R}}_p^k\right\}\hfill & \mathrm{if}{\widehat{R}}_p^k<0,\hfill \\ {\widehat{R}}_p^k\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(16)

where ${\widehat{R}}_{\mathrm{ref},p}^k$ is the final reflectance value at pixel p for channel k. Because we consider the reflectance that is estimated by this model (Equations (14)–(16)) as a reference, it is denoted as ${\widehat{\mathbf{R}}}_{\mathrm{ref}}^{\left(K\right)}$.

4. Outdoor Reflectance Estimation with Reference Devices in the Scene

Because illumination varies during the acquisitions of outdoor scene images, the reflectance estimation method that is described by Equation (14) is not adapted to linescan cameras, such as the Snapscan. In such a case, one solution is to use several reference devices [16].

As a first reference device, we use a white diffuser tile mounted on the acquisition system, so that the sensor vertically observes a portion of it (see Figure 3a). Therefore, the pixel subset $\mathcal{WD}$ contains (about 10%) right border pixels that represent the white diffuser, as shown in Figure 3b. Because $\mathcal{WD}$ spans all the image rows, we further extract a small white square $\mathcal{WS}$ that represents a sample of this reference device. Each acquired image also contains a GretagMacbeth${}^{™}$ ColorChecker that is principally used to assess the performances that are reached by reflectance estimation methods. The pixel subset $\mathcal{WP}$ representing the ColorChecker white patch is used as a second reference device by the double white diffuser ($\mathrm{dwd}$) method [16] that we have adapted to our Snapscan acquisitions, as described in Appendix A.

Although the vignetting effect only depends on the intrinsic camera properties, this method corrects it in each acquired image. In Section 4.1, we propose performing this correction by the analysis of the white diffuser image ${\mathbf{I}}^{\left(B\right)}\left[\mathrm{WD}\right]$. Subsequently, we present state-of-the art estimation reflectance methods that only require one reference device in the scene, but assume that the illumination is constant during the frame acquisition. We finally propose a single-reference method to estimate reflectance in the case of varying illumination.

4.1. Vignetting Correction

The $\mathrm{dwd}$ method acquires a full-field white diffuser image before each scene image acquisition in order to correct the vignetting effect [16]. However, this procedure can be cumbersome, since it requires an external intervention in order to place/remove the full-field white diffuser. Other methods that are presented in the following only require correcting it only once. Because we consider that vignetting only depends on the intrinsic geometric properties of the camera, we propose correcting it thanks to the analysis of a single full-field white diffuser image ${\mathbf{I}}^{\left(B\right)}\left[\mathrm{WD}\right]$ acquired in a laboratory under controlled illumination conditions.

The vignetting effect refers to a loss in the intensity values from the image center to its borders due to the geometry of the sensor optics. To highlight how the vignetting effect would affect radiance measurements, let us rewrite Equation (9) under the Dirac SSF assumption as:

$$I_p^b=\tau ·E_{t_p^b}\left({\lambda }^b\right)·R_p\left({\lambda }^b\right)·A_p\left({\lambda }^b\right).$$

(17)

In order to compensate for the spatial variation of $A_p\left({\lambda }^b\right)$, we compute a correction factor at p, because it requires no knowledge regarding the optical device behavior [19]. Being deduced from the full-field white diffuser image ${\mathbf{I}}^{\left(B\right)}\left[\mathrm{WD}\right]$, the correction factor is channel-wise and pixel-wise computed as:

$$C_p^b={\displaystyle \frac{{\underline{I}}^b\left[\mathrm{WD}\right]}{I_p^b\left[\mathrm{WD}\right]}},$$

(18)

where ${\underline{I}}^b\left[\mathrm{WD}\right]$ is the median value of the m pixels ($m=11$ in our experiments) with the highest values over $I_p^b\left[\mathrm{WD}\right]$, which discards saturated or defective pixel values. The correction factors are stored in a B-channel multispectral image, denoted as $\mathbf{C}$.

Because $\mathbf{C}$ is deduced from a single white diffuser image, it would be corrupted by noise (even after thermal noise removal during the frame acquisitions). Thus, we propose to directly denoise $\mathbf{C}$ by convolving each of its channels $C^b$ with an $11\times 11$ averaging filter $\mathcal{H}$:

$${\tilde{C}}^b=\mathcal{H}\ast C^b.$$

(19)

The vignetting effect in the B-channel radiance image ${\mathbf{I}}^{\left(B\right)}$ is corrected channel-wise and pixel-wise using the smoothed correction factors:

$${\tilde{I}}_p^b={\tilde{C}}_p^b·I_p^b,$$

(20)

where $I_p^b$ and ${\tilde{I}}_p^b$ are the intensity values before and after vignetting correction. This procedure should reduce noise while preserving image textures. We assume that the attenuation is spatially uniform after vignetting correction (i.e., $A_p\left({\lambda }^b\right)·{\tilde{C}}_p^b={\alpha }^b\in \mathbb{R}$ for any given channel index b and pixel p), such that each value of the vignetting-free radiance image is expressed from Equation (17) as:

$${\tilde{I}}_p^b=\tau ·R_p\left({\lambda }^b\right)·E_{t_p^b}\left({\lambda }^b\right)·{\alpha }^b.$$

(21)

4.2. Reflectance Estimation with One Reference Device under Constant Illumination

The illumination $E_{t_p^b}\left({\lambda }^b\right)$ that is associated to p can be determined using the radiance measured at a white diffuser pixel $p_{WD}\in \mathcal{WD}$. To determine illumination thanks to a single white diffuser as reference device, the methods in the literature often assume that illumination is constant, i.e., $E_{t_p^b}\left({\lambda }^b\right)=E_{t_{p_{WD}}^b}\left({\lambda }^b\right)=E\left({\lambda }^b\right)$. Equation (21) then becomes ${\tilde{I}}_p^b=\tau ·R_p\left({\lambda }^b\right)·E\left({\lambda }^b\right)·{\alpha }^b$ and, specifically, ${\tilde{I}}_{p_{WD}}^b=\tau ·{\rho }_{\mathrm{wd}}·E\left({\lambda }^b\right)·{\alpha }^b$, since the white diffuser has a homogeneous diffuse reflection $R_{p_{WD}}\left({\lambda }^b\right)={\rho }_{\mathrm{wd}}$ (95% in our case). The reflectance at p is then deduced from ${\tilde{I}}_p^b$ and ${\tilde{I}}_{p_{WD}}^b$ as:

$$R_p\left({\lambda }^b\right)={\rho }_{\mathrm{wd}}·{\displaystyle \frac{{\tilde{I}}_p^b}{{\tilde{I}}_{p_{WD}}^b}}.$$

(22)

To be robust against spatial noise, the white-average ($\mathrm{wa}$) method [2,20] averages all of the values over the white diffuser pixel subset $\mathcal{WS}$ (see Figure 3b) and estimates the reflectance at each image pixel as:

$${\widehat{R}}_{\mathrm{wa},p}^b={\rho }_{\mathrm{wd}}·{\displaystyle \frac{{\tilde{I}}_p^b}{{\displaystyle \frac{1}{\left|\mathcal{WS}\right|}}{\sum }_{s\in \mathcal{WS}}{\tilde{I}}_s^b}},$$

(23)

where $\left|·\right|$ is the set cardinal.

Similarly, the max-spectral ($\mathrm{ms}$) method [15] assumes that the pixel with maximum value within each channel can be considered to be a white diffuser pixel for estimating the illumination. While ignoring the diffuse reflection factor, reflectance is estimated at each pixel in each channel by the $\mathrm{ms}$ method, as:

$${\widehat{R}}_{\mathrm{ms},p}^b={\displaystyle \frac{{\tilde{I}}_p^b}{{max}_{s\in X}{\tilde{I}}_s^b}},$$

(24)

where X contains all of the image pixels, except $\mathcal{WD}$, and those of the ColorChecker.

The $\mathrm{wa}$ and $\mathrm{ms}$-based B-channel reflectance images undergo spectral correction and negative value removal (see Equations (15) and (16)) to provide the final K-channel reflectance images ${\widehat{\mathbf{R}}}_{\mathrm{wa}}^{\left(K\right)}$ and ${\widehat{\mathbf{R}}}_{\mathrm{ms}}^{\left(K\right)}$.

4.3. Reflectance Estimation with One Reference Device under Varying Illumination

In varying illumination conditions, the Snapscan acquires each row at a given time, hence under a specific illumination (see Section 2.4). Hence, reflectance can no longer be estimated, as in Equation (14). Instead, we propose determining the illumination that is associated to each row of the vignetting-free image ${\tilde{\mathbf{I}}}^{\left(B\right)}$ from the white diffuser pixel set $\mathcal{WD}$ [21]. The underlying assumption is that illumination is spatially uniform over each row at both the white diffuser and scene pixels (that may be not verified in the case of shadows).

Based on this row uniformity assumption for illumination, we estimate reflectance from ${\tilde{\mathbf{I}}}^{\left(B\right)}$ in a row-wise manner, as follows. At pixel p with spatial coordinates $x_p$ and $y_p$, Equation (21) can be rewritten as:

$$R_p\left({\lambda }^b\right)={\displaystyle \frac{{\tilde{I}}_p^b}{\tau ·E_{t_{y_p}^b}\left({\lambda }^b\right)·{\alpha }^b}}.$$

(25)

To determine the illumination $E_{t_{y_p}^b}\left({\lambda }^b\right)$ that is associated to the row of p for channel index b, we use a white diffuser pixel $r_{WD}\in \mathcal{WD}$ located on the same row as p. At $r_{WD}$, the reflectance is equal to the white diffuser reflection factor ${\rho }_{\mathrm{wd}}$, and Equation (21) provides the vignetting-free radiance as:

$${\tilde{I}}_{r_{WD}}^b=\tau ·{\rho }_{\mathrm{wd}}·E_{t_{r_{WD}}^b}\left({\lambda }^b\right)·{\alpha }^b.$$

(26)

Because p and $r_{WD}$ are located on the same row, $t_{y_p}^b=t_{r_{WD}}^b$ and $E_{t_{y_p}^b}\left({\lambda }^b\right)=E_{t_{r_{WD}}^b}\left({\lambda }^b\right)$ according to the assumption regarding the spatial uniformity over each row. Therefore, Equation (26) can be rewritten as:

$$E_{t_{y_p}^b}\left({\lambda }^b\right)={\displaystyle \frac{{\tilde{I}}_{r_{WD}}^b}{\tau ·{\rho }_{\mathrm{wd}}·{\alpha }^b}},$$

(27)

which can be considered to be an estimation of the illumination that is associated to pixel p. For robustness sake, we propose computing it from the median value ${\tilde{\underline{I}}}_{\mathcal{WD},y_p}^b$ of the m highest pixel values that represent the white diffuser subset $\mathcal{WD}$ in $y_p$, rather than from a single value ${\tilde{I}}_{r_{WD}}^b$. Plugging Equation (27) in (25) yields our row-wise ($\mathrm{rw}$) reflectance estimation at pixel p for channel index b:

$${\widehat{R}}_{\mathrm{rw},p}^b={\rho }_{\mathrm{wd}}·{\displaystyle \frac{{\tilde{I}}_p^b}{{\tilde{\underline{I}}}_{\mathcal{WD},y_p}^b}}.$$

(28)

In practice, setting $m=11$ pixels is a good compromise for accurately estimating the illumination for each row and each channel.

5. Experiments about Outdoor Reflectance Estimation

We now present the experimental setup and metrics that were used to objectively evaluate the estimated reflectance. The accuracy results are obtained and discussed for the previously described estimation methods as well as for the extra training-based method described in the present section.

5.1. Experimental Setup

An acquisition campaign that was conducted in a greenhouse under skylight (see Figure 4a) provided 109 radiance images of $2048\times 2048$ pixels × 192 channels of 10-bit depth. Among the targeted plants are crops (e.g., beet) and weeds (e.g., thistle and goose-foot). The images were acquired at different dates of May and June 2019, and different day times (see Figure 4b). Figures 6a and 7a show a RGB rendering of two of them with the D65 illuminant.

All of the images contain a GretagMacbeth${}^{™}$ ColorChecker that is composed of 24 patches. From ${\mathbf{I}}^{\left(B\right)}\left[\mathrm{WD}\right]$ and a radiance image ${\mathbf{I}}^{\left(B\right)}\left[\mathrm{CC}\right]$ of our ColorChecker acquired in a laboratory under controlled illumination, we estimate the K-channel reference reflectance image ${\widehat{\mathbf{R}}}_{\mathrm{ref}}^{\left(K\right)}\left[\mathrm{CC}\right]$ of the ColorChecker according to Equations (14)–(16). From ${\widehat{\mathbf{R}}}_{\mathrm{ref}}^{\left(K\right)}\left[\mathrm{CC}\right]$, we compute the K-dimensional reflectance vector (see Figure 5b,c) of each patch ${\mathcal{P}}_j^{}$ as:

$${\widehat{R}}_{\mathrm{ref},{\mathcal{P}}_j}^k\left[\mathrm{CC}\right]={\displaystyle \frac{1}{|{\mathcal{P}}_j|}}\sum _{p\in {\mathcal{P}}_j}{\widehat{R}}_{\mathrm{ref},p}^k\left[\mathrm{CC}\right],$$

(29)

where $|{\mathcal{P}}_j|$ is the number of pixels that characterize the considered patch.

Among the 24 color patches of the ColorChecker chart, we use a learning subset ${\mathcal{P}}^l$ of 12 patches for the learning procedure and the remaining 12 test patches ${\mathcal{P}}^t$ for testing the quality of reflectance estimation (see Figure 5). The learning patches of ${\mathcal{P}}^l$ are selected using an exhaustive search.

Among the 2,704,156 tested combinations, we retain the one that provides the lowest (mean absolute) reflectance estimation error (see Section 5.4).

The test subset ${\mathcal{P}}^t$ is used to assess the performances that are reached by reflectance estimation methods, and the learning one ${\mathcal{P}}^l$ is fed into a training-based reflectance estimation method, as described in the following.

5.2. Training-Based Reflectance Estimation

The linear Wiener ($\mathrm{wn}$) estimation technique can be applied to estimate reflectance thanks to a learning procedure [23]. It is based on a matrix $\mathbf{G}$ that transforms radiance spectra into reflectance. From any radiance image ${\mathbf{I}}^{\left(B\right)}$ in the database, we compute the spectrally-corrected vignetting-free radiance image ${\tilde{\mathbf{I}}}^{\left(K\right)}$ while using Equations (15) and (20), and then estimate the K-channel reflectance image as:

$${\widehat{\mathbf{R}}}_{\mathrm{wn},p}^{\left(K\right)}=\mathbf{G}·{\tilde{\mathbf{I}}}_p^{\left(K\right)}.$$

(30)

To compute $\mathbf{G}$, we use the spectra of the ColorChecker learning patches (${\mathcal{P}}^l$ subset) that are represented in each of our images. The estimation matrix G that is associated to each input radiance image is determined as:

$$\mathbf{G}={\mathbf{T}}_{ref}·{\mathbf{T}}_{rad}^{⊺}{\left({\mathbf{T}}_{rad}·{\mathbf{T}}_{rad}^{⊺}\right)}^{-1},$$

(31)

where ${\mathbf{T}}_{ref}$ and ${\mathbf{T}}_{rad}$ are the $K\times 12$ matrices that are formed by horizontally stacking the centered and transposed reference reflectance vectors (from ${\widehat{\mathbf{R}}}_{\mathrm{ref}}^{\left(K\right)}\left[\mathrm{CC}\right]$) and radiance vectors (from the current image ${\tilde{\mathbf{I}}}^{\left(K\right)}$) of the learning patches, and ${}^{⊺}$ denotes the transpose.

5.3. Evaluation Metrics

To evaluate the accuracy of reflectance estimation, we use the patches of the ColorChecker test subset ${\mathcal{P}}_{}^t$ (see Figure 5). Let ${\widehat{\mathbf{R}}}_{\ast ,{\mathcal{P}}_j^t}^{\left(K\right)}$, $\ast \in \left\{\mathrm{rw},\mathrm{wa},\mathrm{ms},\mathrm{dwd},\mathrm{wn}\right\}$ denote the reflectance image that is estimated for patch ${\mathcal{P}}_j^t\in {\mathcal{P}}_{}^t$ by either the proposed $\mathrm{rw}$ method (see Equation (28)) or the four implemented state-of-art methods (see Equations (23)–(30)).

This vector is compared to the reference reflectance ${\widehat{\mathbf{R}}}_{\mathrm{ref},{\mathcal{P}}_j^t}^{\left(K\right)}\left[\mathrm{CC}\right]$ of the same patch computed according to Equation (29). The spectra of the ColorChecker patches should be similar (and ideally superposed) to their laboratory counterparts when outdoor reflectance is well estimated.

We objectively assess each estimated reflectance image thanks to the mean absolute error (MAE) and angular error $\Delta \theta$ of each test patch ${\mathcal{P}}_j^t\in {\mathcal{P}}^t$ given by:

$$\mathrm{MAE}({\widehat{\mathbf{R}}}_{\mathrm{ref},{\mathcal{P}}_j^t}^{\left(K\right)}[\mathrm{CC}],{\widehat{\mathbf{R}}}_{\ast ,{\mathcal{P}}_j^t}^{\left(K\right)})={\displaystyle \frac{1}{K}}\sum _{k=1}^K\left|{\widehat{R}}_{\mathrm{ref},{\mathcal{P}}_j^t}^k\left[\mathrm{CC}\right]-{\widehat{R}}_{\ast ,{\mathcal{P}}_j^t}^k\right|,$$

(32)

and:

$$\Delta \theta ({\widehat{\mathbf{R}}}_{\mathrm{ref},{\mathcal{P}}_j^t}^{\left(K\right)}[\mathrm{CC}],{\widehat{\mathbf{R}}}_{\ast ,{\mathcal{P}}_j^t}^{\left(K\right)})=arccos\left({\displaystyle \frac{〈{\widehat{R}}_{\mathrm{ref},{\mathcal{P}}_j^t}^k[\mathrm{CC}],{\widehat{R}}_{\ast ,{\mathcal{P}}_j^t}^k〉}{{∥{\widehat{\mathbf{R}}}_{\mathrm{ref},{\mathcal{P}}_j^t}^{\left(K\right)}\left[\mathrm{CC}\right]∥}_2·{∥{\widehat{\mathbf{R}}}_{\ast ,{\mathcal{P}}_j^t}^{\left(K\right)}∥}_2}}\right),$$

(33)

where ${∥·∥}_2$ is the Euclidean norm. When $\Delta \theta$ between two vectors (spectra in our case) is equal to zero, it means that these two vectors are collinear.

5.4. Results

We compute the mean absolute error ${\mathrm{MAE}}_{\ast }$ and angular error ${\Delta \theta }_{\ast }$ averaged over all of the test patches of all reflectance images estimated from the whole database to obtain aggregated metrics.

Table 1. The reflectance estimation errors. Bold shows the best result and italics the second best one.

Method	Illumination-Based				Training-Based
	$\mathbf{rw}$	$\mathbf{wa}$	$\mathbf{ms}$	$\mathbf{dwd}$	$\mathbf{wn}$
${\mathrm{MAE}}_{\ast }$(%)	4.315	5.883	14.670	3.236	3.628
${\Delta \theta }_{\ast }$ (rad)	0.046	0.046	0.309	0.063	0.063

presents the results for the five tested methods.

The $\mathrm{MAE}$ and $\Delta \theta$ are complementary metrics and they, respectively, highlight two important properties: the scale and shape of the estimated spectra. Indeed, while $\mathrm{MAE}$ is mainly sensitive to the scale of the estimated spectra, $\Delta \theta$ especially focuses on the shape of the spectra, because it is a scale-insensitive measure. Consequently, there might be no correlation between the results that were obtained by the $\mathrm{MAE}$ measure and those obtained by $\Delta \theta$.

No method provides the best results according to the two metrics, as we can see from Table 1. Indeed, the $\mathrm{wn}$ and $\mathrm{dwd}$ methods provide better results than $\mathrm{rw}$ and $\mathrm{wa}$ according to the $\mathrm{MAE}$, but the $\mathrm{rw}$ and $\mathrm{wa}$ methods provide better results in terms of $\Delta \theta$.

The $\mathrm{ms}$ method provides the worst results, because it only analyzes pixels of background and vegetation that strongly absorb the incident light in the visible domain. Hence, the biased illumination estimation in this domain affects the performance of $\mathrm{ms}$ method. It is worthwhile to mention that the $\mathrm{wn}$ method performance might also be biased, since it uses some of the ColorChecker patches as training references (to build estimation matrix $\mathbf{G}$), while the other patches of the same chart are used to evaluate the reflectance estimation quality.

Among illumination-based methods that analyze a single reference device, $\mathrm{rw}$ provides similar results to $\mathrm{wa}$ in terms of $\Delta \theta$, as well as better $\mathrm{MAE}$ results. This shows that taking account of the illumination variation during the frame acquisitions improves the reflectance estimation quality.

6. Multispectral Image Segmentation

Now, we evaluate the contribution of our proposed $\mathrm{rw}$-based reflectance estimation method for supervised crop/weed detection and identification. For this experiment, we focus on the beet (crop) that must be distinguished from thistle and goose-foot (weeds). First, vegetation pixels are detected and ground truth (labels) regarding vegetation pixels is provided by an expert in agronomy (Section 6.1). In order to evaluate the robustness of each considered feature against illumination conditions, we use a data set composed of 37 radiance (13 single-species and 24 mixed) images that we split into a learning and test set, denoted as ${\mathcal{S}}^{learn}$ (23 images) and ${\mathcal{S}}^{test}$ (14 images) (Section 6.2). The illumination conditions are various in the two sets and ${\mathcal{S}}^{test}$ mostly includes images that are acquired on different days from those of ${\mathcal{S}}^{learn}$ (see Figure 4). Note that, as a consequence, vegetation in the learning and test image sets may not be exactly at the same growth stages. We first compare the discrimination power of reflectance features provided by our $\mathrm{rw}$ method against radiance features to assess each reflectance estimation method for crop/weed identification and detection. Subsequently, we compare it with reflectance features that are estimated using each of the four considered state-of-the-art methods ($\mathrm{wa},\mathrm{ms},\mathrm{wn}$, and $\mathrm{dwd}$) (Section 6.3, Section 6.4 and Section 6.5).

6.1. Vegetation Pixel Extraction and Labelling

Only vegetation pixels are analyzed because we aim to detect/identify crops and weeds. They are distinguished from the background (white diffuser, ColorChecker, and soil pixels) using the normalized difference vegetation index (NDVI) [24]. We compute the NDVI values from the $\mathrm{rw}$-based reflectance image ${\widehat{\mathbf{R}}}_{\mathrm{rw}}^{\left(K\right)}$, since the $\mathrm{rw}$ method considers illumination variation, but the images provided by any other reflectance estimation method should yield similar vegetation pixel detection results. We consider p to be a vegetation pixel if its NVDI value is greater than a threshold $\gamma$:

$${\displaystyle \frac{{\widehat{R}}_{\mathrm{rw},p}^{139}-{\widehat{R}}_{\mathrm{rw},p}^{67}}{{\widehat{R}}_{\mathrm{rw},\mathrm{p}}^{139}+{\widehat{R}}_{\mathrm{rw},p}^{67}}}\ge \gamma ,$$

(34)

with the Snapscan “virtual” band centers ${\lambda }^{67}=678.2\mathrm{nm}$ and ${\lambda }^{139}=899.2\mathrm{nm}$. Setting $\gamma =0.45$ experimentally provides a good compromise between under- and over-segmentation of vegetation pixels. Noisy vegetation pixels are filtered out as much as possible by morphological opening. The vegetation pixels are then manually labelled by an expert in agronomy to build the segmentation ground truth for each multispectral image.

6.2. Learning and Test Vegetation Pixels

From the learning set ${\mathcal{S}}^{learn}$, we randomly extract $\mathcal{N}$ learning pixels per class. For a given class ${\mathcal{C}}_i$, $i=[0,\cdots ,N_{\mathcal{C}}-1]$, the number of extracted learning pixels per image depends on the number of images where class ${\mathcal{C}}_i$ is represented in ${\mathcal{S}}^{learn}$ (occurrences). Among the 23 learning images, the beet (crop) class appears in 17 images, thistle in nine images, and goosefoot in 12 images.

In the test set ${\mathcal{S}}^{test}$, beet, thistle, and goosefoot are represented, respectively, in 12, 10, and four images. For the weed detection task, we extract $2\mathcal{N}$ learning pixels, half for crop and half for weed class. Because we merge thistle and goosefoot prototype pixels to build a single weed class, we extract $\mathcal{N}/2$ learning pixels for thistle and $\mathcal{N}/2$ for goosefoot.

Each pixel is characterized by a K-dimensional ($K=141$) feature vector of reflectance (or radiance) values. The reflectance/radiance images are averaged channel-wise over a $5\times 5$ pixel window to reduce noise and within-class variability.

Table 2. The number of learning and test pixels for crop/weed detection (left sub-column) and beet/thistle/goosefoot identification (right sub-column) ($\mathcal{N}=400,000$ pixels in this experiment).

		${\mathcal{S}}^{\mathbf{learn}}$ (23 Images)			${\mathcal{S}}^{\mathbf{test}}$ (14 Images)
Class${\mathcal{C}}_{\mathbf{i}}$		#Occurrences	# Learning Pixels per Occurrence		#Occurrences	#Test Pixels per Class
Crop	Beet	17	$\mathcal{N}/17$		12	5,714,326
Weed	Thistle Goosefoot	9 12	$\mathcal{N}/9$ $\mathcal{N}/12$	$(\mathcal{N}/2)/9$ $(\mathcal{N}/2)/12$	10 4	6,744,633	5,461,013 1,283,620

shows the number of learning and test pixels per class for weed detection and beet/thistle/goosefoot identification. All of the available pixels in ${\mathcal{S}}^{test}$ are used to assess the generalization power of a supervised classifier.

6.3. Evaluation Metrics

The classical accuracy score can be a misleading measure to evaluate a classifier performance when the number of test pixels that are associated to each class is highly skewed (like in Table 2) [25] (p. 114). A classification model that predicts the majority class for all test pixels reaches a high classification accuracy. However, this model can also be considered as weak when misclassifying pixels of the minority classes is worse than missing pixels from the majority classes. In order to overcome this so-called “accuracy paradox”, the performance of a classification model for imbalanced datasets should be summarized with appropriate metrics, such as precision/recall curve [25] (pp. 53–56, 114). Although some metrics may be more meaningful and easy to interpret, there is no consensus in the literature for choosing a single optimal metric. In our case, we want to correctly detect weed pixels without over-detection, because this would imply spraying crops with herbicides. Therefore, the performance of our classification model on both crop and weed detection/identification should be comparable. For this purpose, we use the per-class accuracy score and the weighted overall accuracy score. We also compute the F1-score that combines the precision and recall measures. These three measures should summarize the classification performance of imbalanced sets of test pixels well.

Let us denote the true test pixel labels as y and the set of predicted labels as $\widehat{y}$. The per-class accuracy score for class ${\mathcal{C}}_i$, $i=[0,\cdots ,N_{\mathcal{C}}-1]$, is:

$${\mathrm{Accuracy}}_{{\mathcal{C}}_i}={\displaystyle \frac{1}{|y_{{\mathcal{C}}_i}|}}\sum _{j=0}^{|y_{{\mathcal{C}}_i}|-1}\delta ({\widehat{y}}_{{\mathcal{C}}_{ij}}=y_{{\mathcal{C}}_{ij}}),$$

(35)

where $y_{{\mathcal{C}}_{ij}}$ and ${\widehat{y}}_{{\mathcal{C}}_{ij}}$ are the true and predicted labels for the j-th test pixel of class ${\mathcal{C}}_i$, respectively, and $|y_{{\mathcal{C}}_i}|$ is the number of test pixels of class ${\mathcal{C}}_i$.

The weighted overall accuracy for binary and multiclass classifications is defined as:

$$\overline{\mathrm{Accuracy}}={\displaystyle \frac{{\sum }_{i=0}^{N_{\mathcal{C}}-1}{\omega }_{{\mathcal{C}}_i}·{\mathrm{Accuracy}}_{{\mathcal{C}}_i}}{{\sum }_{i=0}^{N_{\mathcal{C}}-1}{\omega }_{{\mathcal{C}}_i}}},$$

(36)

where ${\omega }_{{\mathcal{C}}_i}=1/\left|y_{{\mathcal{C}}_i}\right|$ is the weight that is associated to class ${\mathcal{C}}_i$ and computed as the inverse of its size, so as to handle imbalanced classes.

Because the F1-score privileges the classification of true positives pixels (weed pixels in our case), we compute the overall $\overline{\mathrm{F}1}$-score as the population-weighted F1-score, so that the performances over all classes are considered:

$$\overline{\mathrm{F}1}={\displaystyle \frac{{\sum }_{i=0}^{N_{\mathcal{C}}-1}{\omega }_{{\mathcal{C}}_i}·\mathrm{F}{1}_{{\mathcal{C}}_i}}{{\sum }_{i=0}^{N_{\mathcal{C}}-1}{\omega }_{{\mathcal{C}}_i}}}.$$

(37)

The F1-score of class ${\mathcal{C}}_i$ is computed as:

$$\mathrm{F}{1}_{{\mathcal{C}}_i}=2·{\displaystyle \frac{{\mathrm{Precision}}_{{\mathcal{C}}_i}·{\mathrm{Recall}}_{{\mathcal{C}}_i}}{{\mathrm{Precision}}_{{\mathcal{C}}_i}+{\mathrm{Recall}}_{{\mathcal{C}}_i}}},$$

(38)

where

$${\mathrm{Precision}}_{{\mathcal{C}}_i}={\displaystyle \frac{{\sum }_{j=0}^{|y_{{\mathcal{C}}_i}|-1}\delta ({\widehat{y}}_{{\mathcal{C}}_{ij}}=y_{{\mathcal{C}}_{ij}})}{|{\widehat{y}}_{{\mathcal{C}}_i}|}},$$

(39)

and

$${\mathrm{Recall}}_{{\mathcal{C}}_i}={\displaystyle \frac{{\sum }_{j=0}^{|y_{{\mathcal{C}}_i}|-1}\delta ({\widehat{y}}_{{\mathcal{C}}_{ij}}=y_{{\mathcal{C}}_{ij}})}{|y_{{\mathcal{C}}_i}|}}.$$

(40)

6.4. Classification Results

The parametric LightGBM (LGBM) and non-parametric Quadratic Discriminant Analysis (QDA) classifiers are applied for supervised weed detection and identification problems. The choice of these two non-linear classifiers is motivated by their processing time during the learning and prediction procedures and their fundamentally different decision rules. Indeed, LGBM is a parametric tree-based classifier that requires a learning procedure to model a complex classification rule, whereas QDA is a simple non-parametric classifier that is based on Bayes’ theorem to perform predictions. For LGBM, we retain the default parameter values (learning rate of 0.05, 150 leaves) and use the log loss function as the learning evaluation metric. LGBM uses a histogram-based algorithm to bucket the features into discrete bins, which drastically reduces the memory and time consumption. The number of bins is set to 255 and the number of boosting operations to 100. Additionally, the feature fraction and bagging fraction parameters are set to 0.8 to increase LGBM speed and avoid over-fitting.

Table 3. Crop/weed detection and beet/thistle/goosefoot identification results with the QDA and LGBM classifiers. Bold shows best result and italics second best one.

			Radiance	Reflectance Features
		Classifier	Feature	Illumination-Based				Training-Based
				rw	$\mathbf{wa}$	$\mathbf{ms}$	$\mathbf{dwd}$	$\mathbf{wn}$
Crop/weed detection	${\overline{\mathrm{Accuaracy}}}_{}$ (%)	QDA	71.6	76.0	73.3	79.2	76.9	48.1
		LGBM	77.0	86.1	84.7	81.8	83.4	71.0
	${\overline{\mathrm{F}1}}_{}$ (%)	QDA	67.9	76.0	73.0	77.8	76.5	48.0
		LGBM	76.2	85.4	84.4	81.6	83.1	71.3
Beet/thistle/goosefoot identification	${\overline{\mathrm{Accuracy}}}_{}$ (%)	QDA	34.4	45.1	41.8	26.0	37.3	47.8
		LGBM	39.5	49.7	50.3	46.5	46.0	34.0
	${\overline{\mathrm{F}1}}_{}$ (%)	QDA	34.3	41.7	40.1	25.4	33.7	27.0
		LGBM	39.7	47.1	44.4	41.4	42.5	32.9

shows the classification results that were obtained with LGBM and QDA classifiers for each considered feature. Figure 6 and Figure 7 show the color-coded vegetation pixel classification of two test images using the LGBM classifier in weed detection and identification tasks, respectively.

Let us first compare the classification performance of reflectance against radiance features for the weed detection task. From the results that are given in Table 3, we can see that reflectance features estimated by illumination-based methods ($\mathrm{rw}$, $\mathrm{wa}$, $\mathrm{ms}$, and $\mathrm{dwd}$) provide better classification results than radiance features in terms of the average F1 and accuracy scores, whatever the classifier. The worst classification results are obtained with reflectance features that are estimated using the $\mathrm{wn}$ method. Training-based methods, such as $\mathrm{wn}$, can provide an accurate reflectance estimation of scene objects whose optical properties are close to those of the training samples. In our case, the optical properties of vegetation are very different from that of the training ColorChecker patches. Thus, $\mathrm{wn}$ provides inaccurate reflectance estimations at vegetation pixels, which affects its classification performance. Figure 6 illustrates the satisfying $\overline{\mathrm{Accuracy}}$ and $\overline{\mathrm{F}1}$ scores that are obtained thanks to the analysis of illumination-based reflectance features by LGBM. Indeed, this figure shows that weed is globally well detected by these methods.

For weed identification, the classification performances of all the features are degraded, because they provide weak performances on the goosefoot class (see Figure 7). This lack of generalization might be due to the high within-class dispersion (since we consider vegetation at various growth stages) and/or the physiological vegetation changes.

Let us now compare the classification performances of the reflectance features. The best overall classification results are obtained by our proposed $\mathrm{rw}$ method that performs well with both classifiers and reaches the highest average F1 and accuracy scores with LGBM for weed detection (85.4% and 86.1%, respectively). The $\mathrm{wa}$ method provides good classification results, better than those that were obtained by $\mathrm{dwd}$, although the latter accounts for illumination variation during the frame acquisitions. The computation of illumination scaling factors to compensate for illumination variation may explain this poor performance, as well as the loss of spectral information (saturated reflectance values) in the near infra-red domain that is caused by illumination normalization (see Equation (A4)).

6.5. Experimental Conclusions

The experiments with this outdoor image database allow us to compare the performances of different reflectance features according to the estimation quality and pixel classification. The evaluation results are summarized by separately studying weed detection and identification.

Table 4. The ranking of reflectance estimation methods for crop/weed detection. Bold shows best result and italics second best one.

Evaluation Criterion			Method
			$\mathbf{rw}$	$\mathbf{wa}$	$\mathbf{ms}$	$\mathbf{dwd}$	$\mathbf{wn}$
${\overline{\mathrm{MAE}}}_{\ast }$ (%)			3	4	5	1	2
${\overline{\Delta \theta }}_{\ast }$ (rad)			1	1	5	3	3
Crop/weed detection	${\overline{\mathrm{Accuracy}}}_{\ast }$	QDA	3	4	1	2	5
		LGBM	1	2	4	3	5
	${\overline{\mathrm{F}1}}_{\ast }$	QDA	3	4	1	2	5
		LGBM	1	2	4	3	5
${\sum }_{\diamond =1}^6{\mathrm{Rank}}_{\diamond ,\ast }$			12	17	20	14	25

and

Table 5. The ranking of reflectance estimation methods for crop/weed identification. Bold shows best result and italics second best one.

Evaluation Criterion			Method
			$\mathbf{rw}$	$\mathbf{wa}$	$\mathbf{ms}$	$\mathbf{dwd}$	$\mathbf{wn}$
${\overline{\mathrm{MAE}}}_{\ast }$ (%)			3	4	5	1	2
${\overline{\Delta \theta }}_{\ast }$ (rad)			1	1	5	3	3
Beet/thistle/goosefoot identification	${\overline{\mathrm{Accuracy}}}_{\ast }$	QDA	2	3	5	4	1
		LGBM	2	1	3	4	5
	${\overline{\mathrm{F}1}}_{\ast }$	QDA	1	2	5	3	4
		LGBM	1	2	4	3	5
${\sum }_{\diamond =1}^6{\mathrm{Rank}}_{\diamond ,\ast }$			10	13	27	18	20

show the rank ${\mathrm{Rank}}_{\diamond ,\ast }$ obtained by each reflectance estimation method ∗ according to each evaluation criterion ⋄ used in Table 1 and Table 3. The method with the lowest total rank is considered to be the best one, since it satisfies several criteria.

The total ranks of $\mathrm{ms}$ and $\mathrm{wn}$ methods are the highest ones, because they provide the worst results for either estimation quality ($\mathrm{ms}$) or classification performance ($\mathrm{wn}$). On the one hand, the $\mathrm{dwd}$ method that uses two reference devices to cope with illumination variation provides the second best total rank for weed detection (see Table 4). On the other hand, the $\mathrm{wa}$ method that uses one reference device, but assumes that illumination is constant, gives the second best total rank for weed identification (see Table 5). Our $\mathrm{rw}$ method, which row-wise analyzes one single reference device in order to take account of illumination variation, reaches the best total ranks for both weed detection and identification problems. These experiments suggest that $\mathrm{rw}$-based reflectance features are relevant for weed identification under variable illumination conditions. Their performance should also be confirmed with other crop species, such as bean and wheat.

7. Conclusions

This paper first proposes an original image formation model of linescan multispectral cameras, like the Snapscan. It shows how illumination variation during the multispectral image acquisition by this device impacts the measured radiance that is provided by a Lambertian surface element. Our model is versatile and it can be adapted to model the outdoor image acquisition of several multispectral cameras, such as the HySpex VNIR-1800 [26] or the V-EOS Bragg-grating camera used in [16]. From this model, we propose a reflectance estimation method that copes with illumination variation. Because such varying conditions may affect the reflectance estimation quality, we estimate illumination at the frame level using a row-wise ($\mathrm{rw}$) approach. We experimentally show that the $\mathrm{rw}$ method is more robust against illumination variation than the state-of-the-art methods. We also show that $\mathrm{rw}$-based features are more discriminant to target outdoor supervised weed detection and identification, and they provide the best classification results. The accuracy of weed recognition systems and their robustness against illumination can be improved using reflectance features. This allows for precision spraying techniques to be considered in order to get rid of weeds using fewer quantities of chemicals. This study enables to make a step towards sustainable agriculture. As future work, segmentation will be extended to other plant species (such as wheat and bean) and growth stages. Spectral feature selection and texture features extraction will also be studied to improve the crop/weed identification performance.