Tiancheng Yang
email: tiancheng.yang@inrae.fr

Sylvain Jay

Yangmingrui Gao

Shouyang Liu

Frederic Baret

The balance between spectral and spatial information to estimate straw cereal plant density at early growth stages from optical sensors

Keywords: Plant density, spectral reflectance, spatial resolution, wheat, barley

Estimating straw cereal plant density at early stages is important for field crop management and phenotyping. Usual plant density estimation methods include manual counting and imagebased counting, both of which have limited throughput, due to the need for high spatial resolution images. In this study, we explored the potential of high-throughput estimations with spectral information. A large and diverse dataset was collected on micro plot field experiments, encompassing six sites, three leaf stages, and four species of straw cereals. Canopy spectral reflectance was acquired with a spectrometer, both in 0° or 45° view zenith angle, perpendicularly to the row direction. Two reflectance-based approaches were then tested. In the direct approach, density was directly estimated from reflectance using Gaussian process regression (GPR) and spectral bands selected based on Akaike's information criterion. In the indirect approach, the green fraction derived from high spatial resolution RGB images (GF_rgb) was first estimated from reflectance using GPR and selected bands, and then linearly related to density. These reflectance-based methods were compared to a classical image-based baseline method, which estimates density directly from GF_rgb.

An ablation study firstly showed the superiority of 45° observations, and the necessity to calibrate one model for each site, growth stage, and species. The band selection process recommended using no more than four bands as inputs to the GPR models. The resulting direct and indirect estimations had an overall relative error of 30%. The image-based baseline method had a lower error of 22% for submillimeter spatial resolutions, but it performed worse than reflectance-based methods when degrading the spatial resolution to more than 1 to 2 mm to mimic an increase in sensor altitude. These results showed that spectral information can compensate for spatial information and that spectral methods can potentially provide highthroughput and reasonably accurate estimates of straw cereal plant density.

Introduction

Plant density is a fundamental factor in the management and phenotyping of straw cereal crops since it may directly impact the final yield [START_REF] Valério | Seeding density in wheat: the more, the merrier?[END_REF]. However, [START_REF] Whaley | The physiological response of winter wheat to reductions in plant density[END_REF] showed that a larger number of tillers could compensate for a lower plant density in winter wheat. Furthermore, plant density generally reduces weed development by increasing competition for resources [START_REF] Carlson | Wild oat (Avena fatua) competition with spring wheat: plant density effects[END_REF][START_REF] Kristensen | Crop Density, Sowing Pattern, and Nitrogen Fertilization Effects on Weed Suppression and Yield In Spring Wheat[END_REF][START_REF] Lutman | A review of the effects of crop agronomy on the management of A lopecurus myosuroides[END_REF][START_REF] Olsen | How important are crop spatial pattern and density for weed suppression by spring wheat?[END_REF][START_REF] Tollenaar | Effect of crop density on weed interference in maize[END_REF][START_REF] Wilson | Predicting the competitive effects of weed and crop density on weed biomass, weed seed production and crop yield in wheat[END_REF]. In the context of plant phenotyping, plant density allows computation of the emergence rate, which is a valuable trait for breeders. Furthermore, in case of difficult emergence conditions, the knowledge of plant density helps the breeder to decide whether a microplot should be kept or not in an experiment. Finally, plant density is a key characteristic that can be used to assess other traits pertinent for breeders such as the growth stage depending on the number of leaves per plant, or the tillering coefficient.

Researchers have been looking for ways to replace laborious and time-consuming manual counting with high-throughput methods based on optical sensor data. These methods can be divided into two categories: (a) image-based methods and (b) reflectance-based methods.

(a) Image-based methods. On the one hand, many methods in this category begin by binarizing the image into vegetation pixels and background pixels based on RGB or multispectral features. Then, optional morphological analysis of the vegetation pixels is carried out. Finally, the results of classification and/or morphological analysis are used to estimate the number of plants.

The works of [START_REF] Gnädinger | Digital counts of maize plants by unmanned aerial vehicles (UAVs)[END_REF], [START_REF] Jin | Estimates of plant density of wheat crops at emergence from very low altitude UAV imagery[END_REF], [START_REF] Liu | Estimation of wheat plant density at early stages using high resolution imagery[END_REF], [START_REF] Liu | A method to calculate the number of wheat seedlings in the 1st to the 3rd leaf growth stages[END_REF], [START_REF] Roth | Repeated multiview imaging for estimating seedling tiller counts of wheat genotypes using drones[END_REF], [START_REF] Shrestha | Shape and size analysis of corn plant canopies for plant population and spacing sensing[END_REF], [START_REF] Tseng | Rice Seedling Detection in UAV Images Using Transfer Learning and Machine Learning[END_REF], and [START_REF] Wilke | Assessment of plant density for barley and wheat using UAV multispectral imagery for high-throughput field phenotyping[END_REF] used methods of this category. On the other hand, methods that do not rely on binarized images have been developed, and these methods are mainly based on deep learning. [START_REF] Shubhra | DeepWheat: Estimating Phenotypic Traits from Crop Images with Deep Learning[END_REF] employed a two-step, deep learning based method to estimate the number of wheat plants in an image: firstly, they segmented RGB images into plant patches with a deep learning module [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF], and then they estimated the amount of wheat plants within each patch using another deep learning module. Some researchers estimated the number of plants or plant organs on various species using deep learning regression, classification, or detection algorithms [START_REF] Liu | High-Throughput Rice Density Estimation from Transplantation to Tillering Stages Using Deep Networks[END_REF][START_REF] Lu | TasselNetV2+: A Fast Implementation for High-Throughput Plant Counting From High-Resolution RGB Imagery[END_REF][START_REF] Mukhtar | Wheat Plant Counting Using UAV Images Based on Semi-supervised Semantic Segmentation[END_REF][START_REF] Tseng | Rice Seedling Detection in UAV Images Using Transfer Learning and Machine Learning[END_REF][START_REF] Wu | Automatic Counting of in situ Rice Seedlings from UAV Images Based on a Deep Fully Convolutional Neural Network[END_REF]. Amongst them, deep learning detection is well-suited to the counting task, and all of these methods have the potential to be applied to plant counting for further density estimation.

(b) Reflectance-based methods. Compared to image-based research, there are fewer reflectance-based studies on plant density estimation. In general, the NDVI value is computed from reflectance measured from ground-based or satellite-borne spectral sensors, and linearly related to plant density. [START_REF] Aase | Determining Winter Wheat Stand Densities Using Spectral Reflectance Measurements[END_REF] showed that the NDVI value is a good proxy of plant density for winter wheat. Although they did not further explore the correlation between NDVI and density, their data showed the potential to create a good linear regression. [START_REF] Reyniers | Optical Measurement of Crop Cover for Yield Prediction of Wheat[END_REF] showed that the crop coverage of wheat obtained from spectral data is more related to sowing density in the early season. This result shows the possibility to estimate wheat seedling density from crop coverage . [START_REF] Habibi | Machine learning techniques to predict soybean plant density using UAV and satellite-based remote sensing[END_REF] combined the accurate deep learning method and the high-throughput reflectance-based method in a two-step soybean plant density estimation.

In the first step, deep-learning method was used to get the plant density with high accuracy, and the density value was used as input for the next step. In the second step, reflectance information and climate information were used to estimate plant density in high-throughput, with a moderate accuracy. [START_REF] Zhang | Stand density estimation based on fractional vegetation coverage from Sentinel-2 satellite imagery[END_REF] estimated the stand density of evergreen trees based on the linear relationship between fractional vegetation cover (FVC) and stand density. They first calibrated the FVC-density relationship on smaller scale areas of 1 hectare, and then applied this method to larger scale areas of about 100 hectares using NDVI calculated from Sentinel 2.

In the case of straw cereal crops, the image-based methods mentioned above often require high spatial resolution images to identify the small leaves observed at early growth stages, when the plants have no tiller and less overlap. For example, [START_REF] Jin | Estimates of plant density of wheat crops at emergence from very low altitude UAV imagery[END_REF] showed that plant density estimation performance decreases with coarser image spatial resolution and thus recommended using spatial resolutions finer than 0.4 mm. Similarly, [START_REF] Liu | Estimation of wheat plant density at early stages using high resolution imagery[END_REF], [START_REF] Liu | A method to calculate the number of wheat seedlings in the 1st to the 3rd leaf growth stages[END_REF], [START_REF] Shubhra | DeepWheat: Estimating Phenotypic Traits from Crop Images with Deep Learning[END_REF], and [START_REF] Mukhtar | Wheat Plant Counting Using UAV Images Based on Semi-supervised Semantic Segmentation[END_REF] used images of 0.2-0.5 mm spatial resolutions. Such spatial resolutions are usually obtained using a high-resolution camera and acquiring images at a low altitude, either from a UAV or from a ground-based system. However, imaging at low altitudes also reduces the throughput, which can be problematic for large fields that need to be sampled in a reasonable time. In this respect, reflectance-based methods, despite being less often used, present interesting advantages over image-based methods: not only does the canopy reflectance remain unchanged as the spatial resolution decreases according to the spectral linear mixing model [START_REF] Adams | Spectral mixture modeling: A new analysis of rock and soil types at the Viking Lander 1 site[END_REF][START_REF] Ritter | Chapter 4 -Lattice Algebra Approach to Endmember Determination in Hyperspectral Imagery[END_REF], but richer spectral information can also potentially compensate for the loss of spatial information. For example, [START_REF] Habibi | Machine learning techniques to predict soybean plant density using UAV and satellite-based remote sensing[END_REF] and [START_REF] Zhang | Stand density estimation based on fractional vegetation coverage from Sentinel-2 satellite imagery[END_REF] have shown high-throughput reflectance-based density estimation is feasible for larger plants (soybean and trees), but further investigation is needed for small crops such as straw cereals. A drawback of the reflectance-based method is that it can be affected by the detrimental influence of soil on canopy reflectance. Several solutions can be implemented to limit this influence. For example, sensing the canopy from a 45° view zenith angle increases the green fraction (GF, the proportion of green vegetation pixels in the sensor field of view) compared to the nadir, capturing more signal from the vegetation and thus increasing the sensitivity of the optical data to changes in plant density for such small plants [START_REF] Jay | Estimating leaf chlorophyll content in sugar beet canopies using millimeter-to centimeter-scale reflectance imagery[END_REF]. Also, using 45° view zenith angle will make the observed GF value less sensitive to the plant leaf inclination angle compared to using smaller angles than 45 ° [START_REF] Weiss | Review of methods for in situ leaf area index (LAI) determination[END_REF], and there will be less overlapping between rows compared to using angles larger than 45°. With 45° view zenith angle, using an azimuth direction perpendicular to the crop row further reduces the overlap between plants inside one row [START_REF] Baret | GAI estimates of row crops from downward looking digital photos taken perpendicular to rows at 57.5° zenith angle: Theoretical considerations based on 3D architecture models and application to wheat crops[END_REF]. This observation set was also used in the work of Liu et al. (2017) and [START_REF] Jin | Estimates of plant density of wheat crops at emergence from very low altitude UAV imagery[END_REF] for wheat seedling density. Besides changing the acquisition geometry, another solution to further limit the soil influence is to constrain the reflectance-based density estimation by first estimating GF, and second relating estimated GF to density. Indeed, canopy reflectance is strongly related to GF [START_REF] Baret | LAI, fAPAR and fCover CYCLOPES global products derived from VEGETATION: Part 1: Principles of the algorithm[END_REF][START_REF] Gitelson | Novel algorithms for remote estimation of vegetation fraction[END_REF], which is itself proportional to plant density when the plants are of similar size with little overlap such as in the case of early-stage plants [START_REF] Wilke | Assessment of plant density for barley and wheat using UAV multispectral imagery for high-throughput field phenotyping[END_REF].

The previous literature review on plant density methods shows that there is currently no comparison of the performances between image-based and reflectance-based methods achieved over the same dataset. Furthermore, the possible degradation of performances as a function of sensor spatial resolution for both types of methods is still lacking. Therefore, in this work, we developed two reflectance-based approaches to estimate plant density from nadir or 45° observations. In the first approach, plant density was estimated directly using a machine learning regression algorithm. In the second approach, GF was first estimated from spectral data, and then related to density. By introducing GF as a proxy, we wanted to make the estimation more interpretable. These two approaches were compared to a popular image-based method, trained on a large and diverse dataset comprising six sites, three leaf stages, and four straw cereal species. In summary, this research has the following objectives:

(a) Evaluate the performance of two reflectance-based methods and one image-based method for estimating cereal straw plant density in terms of accuracy and robustness to changes in sensor spatial resolution.

(b) Evaluate the added value of several strategies to improve the performance of reflectancebased methods, i.e., using 45° instead of 0° observations, and using GF as a proxy for density.

Materials and methods

The experiments

Microplot experiments were conducted in 2021 and 2022 at five sites in France (Avignon, Salin-de-Giraud, Gardanne, Greoux-les-Bains, and Mauguio) and one site in China (Nanjing) (Table 1). The size of each microplot was 1 m * 1.4 m in Avignon, Salin-de-Giraud, Gardanne, and Nanjing, whereas it was mainly 2 m * 12 m in Greoux-les-Bains and 1.4 m * 8 m in Mauguio.

These sites had different soil types and very different soil colors. Both dry and wet soils were included in the experiments. At the Salin-de-Giraud site, no herbicide was used, so there were more weeds. Three density treatments were applied the Avignon, Nanjing, Gardanne, and Salinde-Giraud plots. These treatments resulted in different density values ranging from 37 to 535 plants/m 2 , and most density values were between 100 and 450 plants/m 2 . At Avignon, the plants were sown earlier (September 23) than at the other sites, where more traditional sowing dates (October -January) were used (Table 1). Four cereal crop species (soft wheat, durum wheat, barley, and rye) were considered, as an extreme case of different crop varieties, to explore the possible effects of plant structure on density estimation.

For each site and sowing date, the actual plant density was measured manually on the first day of measurements (section 2.2.1). One to three measurements of spectral reflectance and RGB images were made from the one-to three-leaf stage (sections 2.2.2 and 2.2.3). The other 92 plots at the Greoux-les-Bains and Mauguio sites were larger (section 2.1).

Therefore, for each of these plots, two or three subplots of 1 m in length and two rows wide were selected to be representative of the plot. The densities of the subplots were calculated and averaged to represent the plant density of the plot.

Spectrometer measurements

In France, canopy reflectance data were collected with an SM-3500 spectrometer (Spectral Evolution, Massachusetts, US), with 737 bands ranging from 343 to 2517 nm, and with full-width at half maximum (FWHM) between 1.5 nm and 3.8 nm. In China, canopy reflectance data were collected with an ASD FieldSpec 4 spectrometer (Analytical Spectral Devices, Colorado, US), with bands ranging from 350 to 2500 nm, and with FWHM between 1.1 and 1.4 nm.

Spectral data were measured at 0° and 45° view zenith angles. In the 45° measurements, the azimuth was perpendicular to the direction of the rows. The acquisition geometry was designed such that (1) the area covered by the 25° field of view (FOV) of the spectrometer was similar in 0° and 45° situations, and (2) the area covered by the FOV was large enough to represent the plot while not exceeding the plot boundaries (Figure 1). The spatial resolution of spectrometer measurements, which is defined here as the length of the side of a square with the same area, is 736 mm. In practice, the spectrometers and cameras were held manually so there could be an accidental but slight error in height. The 45° angle was controlled by checking a device with bubble level.

To introduce some variation in the soil background in Avignon, Gardanne, and Nanjing, measurements were performed on both dry soil surface and wet soil surface when it was possible, i.e., when the soil was not already wet due to the rain. A first measurement was made on dry soil.

Then some water was poured onto the surface of the soil to change its color, and the second measurement was made. In the other sites (Greoux-les-Bains, Salin-de-Giraud and Mauguio),

water was not available so only one measurement per plot was performed.

Forty-nine canopy reflectance spectra were removed from the dataset, either because of inaccurate reflectance calibration, or because the shadows of nearby buildings and trees accidentally covered the plots. Among these 49 removed samples, 9 samples were acquired from 0° view zenith angle, and 40 samples were acquired from 45° view zenith angle. In total, there were 262 samples in the nadir view and 231 samples in the 45° view zenith angle for assessing the reflectance-based methods presented in section 2.3. Some preprocessing operations were applied to the spectral data. First, the bands of the ASD spectrometer were interpolated into the bands of the SM-3500 spectrometer. Then, a

Savitzky-Golay filter with a window of 5 and a polynomial order of 2 was applied to reduce the influence of noise [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF]Virtanen et al., 2020). Afterward, only the bands from 343 to 1338 nm and from 1494 to 1798 nm were used in this experiment to avoid the atmospheric water absorption effect around 1400 nm and 1850 nm and to avoid the noise in bands with wavelengths longer than 2000 nm. Next, to accelerate the calculations after, the number of bands was further reduced. Starting from the first band of 343 nm, each time the shortest band that is greater than 10 nm away was chosen, thus making a band set of 343 nm, 354 nm, 365 nm, etc. This reducing process gave 140 bands as output.

RGB imagery

In France, a Sony Alpha 5100 camera (Sony, Inc. Minato, Tokyo, Japan.) with 24M pixels and a 45 mm focal length (in 35 mm equivalent focal length) was used to collect RGB images with a spatial resolution ranging from 0.1 and 0.3 mm at the ground level. In China, a Sony RX0 camera with 15M pixels and a 24 mm focal length (in 35 mm equivalent focal length) was used, and the images had a spatial resolution between 0.2 and 0.5 mm at the ground level. These ranges of spatial resolutions were caused by 1) the variability in pixel size in 45° images due to the variable distance between soil and camera within the imaged scene, and 2) the accidental error in camera height that was controlled manually. The cameras were held at the same place and with the same orientation as the spectrometer. Since the cameras had a larger FOV than the spectrometers, the image contained parts that did not belong to the target plot, so these parts were cropped during preprocessing. The images were collected on both dry soil and wet soil with the same method as in section 2.2.2.

Direct and indirect density estimation methods

from canopy reflectance with Gaussian process regression

Description of direct, indirect and baseline methods

Two density estimation methods based on canopy reflectance were compared in this work (Figure 2). The first one was a one-step direct estimation approach, in which canopy reflectance was related to plant density using a Gaussian process regression (GPR) model.

The second method was a two-step indirect estimation approach, in which GF was used as a proxy of plant density. The first step was to estimate GF from canopy reflectance with GPR. The ground truth values of GF (GF_rgb) were derived from high spatial resolution RGB images with the SegVeg deep learning segmentation method developed by our team [START_REF] Madec | VegAnn, Vegetation Annotation of multi-crop RGB images acquired under diverse conditions for segmentation[END_REF][START_REF] Serouart | SegVeg: Segmenting RGB images into green and senescent vegetation by combining deep and shallow methods[END_REF]. The second step was to estimate plant density from estimated GF, using a linear regression model with a zero intercept, since GF is zero when the density is zero.

In addition to the direct and indirect methods, a baseline method inspired by Wilke et al.

(2021) was applied, fitting the GF_rgb to the density with a linear regression model through the origin, to have a better understanding of the second step of the indirect method (Figure 2). This proportional relationship is based on the hypothesis that the individual plants are of similar sizes and the overlapping can be neglected, and this is the case for the plants in the early stages [START_REF] Wilke | Assessment of plant density for barley and wheat using UAV multispectral imagery for high-throughput field phenotyping[END_REF].

For either GF or plant density estimation, the GPR method here used a similar kernel as Verrelst et al. ( 2013) (Equation 1), except that we did not use Kronecker delta as the multiplier of noise:

(

) 1
where is a scaling factor, is the number of bands, and are the reflectance value of the and samples at the band , is the scale for the reflectance value of each band, and is the standard deviation of noise. Different sets of bands were used as inputs in the ablation study and the final performance evaluation. The ablation study employed five common bands (section 2.3.2), and the final performance evaluation used the optimal bands identified with forward band selection (section 2.3.3). The model should not use hundreds of bands as input, as more bands will contain more redundant information and noise, which are not really related to GF or plant density, and may cause the model to overfit more easily [START_REF] Verrelst | Spectral band selection for vegetation properties retrieval using Gaussian processes regression[END_REF]. used calibrated models for each growth stage and each cultivar, while [START_REF] Liu | Estimation of wheat plant density at early stages using high resolution imagery[END_REF] emphasized the importance of data acquisition in a certain time-window when the plants have 1.5 to 2 leaves.

For each view zenith angle (0° and 45°), we thus performed an ablation study to identify which factor(s), among species, site, and growth stage, should be differentiated when calibrating a regression model. We tested if the calibration of separate sub-models for one or several factor(s) would significantly improve the estimation performance obtained when not differentiating the above three factors. We tested all the possible factor combinations for each of the three estimation steps presented in section 2.3, i.e., spectra-density for direct estimation and spectra-GF and GF_rgb-density (baseline method) for indirect estimation. The estimation performance was quantified with the root mean squared error (RMSE), and the relative root-mean-square error (rRMSE). The rRMSE was obtained by dividing the RMSE by the mean value of the target variable (GF or density) of the full dataset. Due to the potentially low number of samples attained when differentiating several factors, the RMSE and rRMSE were computed using five-fold crossvalidation when the dataset had more than five samples, and leave-one-out cross-validation otherwise. The cross-validations were replicated ten times with different partitions of datasets at each time, to reduce the random error introduced by the random partition of small datasets.

Finally, the averages of RMSE and rRMSE were calculated for each factor combination, for each view zenith angle, and for each estimation step.

To simplify the ablation study, five common bands were used here as spectral reflectance input to the GPR regression model (section 2.3). These five bands were the blue (B, 475 nm), green (G, 560 nm), red (R, 668 nm), red edge (RE, 717 nm), and near-infrared (NIR, 842 nm) bands. This set of bands is widely used in commercial multispectral cameras, such as Rededge (Micasense, Washington, US), P4-multispectral (DJI, Shenzhen, China), Airphen (Hiphen, Avignon, France), MiniMCA (Tetracam, Bolton, UK), and Sentera-6X (Sentera, Minneapolis, US).

We took the bands of the Rededge camera as a reference.

Forward band selection based on Akaike's information criteria for final performance evaluation

Based on optimal factor differentiation determined using the ablation study described in section 2.3.2, we evaluated the plant density estimation performance obtained by exploiting the full spectrum instead of just five common bands. Since using 140 spectral bands as inputs to the GPR model would risk overfitting, a forward band selection method was used to determine the best input bands for each view zenith angle. These bands were optimized for Spectra-GF estimation (part 1 of the indirect method) and used for both spectrally-based estimations of GF and density. At each iteration of the forward band selection, the best band was chosen among all the candidate bands based on the corrected version of the Akaike Information Criterion (AICc)

and then added to the model. AICc takes into account both the error and the parsimony of the model, and a model with a low AICc value is preferred [START_REF] Burnham | Multimodel inference: understanding AIC and BIC in model selection[END_REF]. AICc was calculated as follows:

(2)

Where is the likelihood of the Gaussian Process estimation, is the number of parameters to be determined in the GPR model and is the number of samples used to calibrate this model. In GPR, the log-likelihood of the model was calculated using a Python sci-kit package [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF][START_REF] Williams | Gaussian processes for machine learning[END_REF], and the number of parameters was the number of bands plus two according to Equation (1). In this study, when multiple sub-models were calibrated on different sub-datasets due to factor differentiation (section 2.3.2), the sum of AICc values from multiple sub-models was calculated as a criterion. When one model was calibrated for the whole dataset, the AICc value of this model was used as a criterion. The optimal set of bands comes from the model with a minimum sum of AICc values. No cross-validation was applied to the models when calculating AICc since the choice of model with AIC and crossvalidation are asymptotically equivalent when maximum likelihood estimation is used [START_REF] Stone | An asymptotic equivalence of choice of model by cross-validation and Akaike's criterion[END_REF].

For the best view zenith angle, plant density estimation performances obtained with direct, indirect, and baseline methods were finally evaluated, using optimal factor differentiation and optimal band set. These performances were quantified using the coefficient of determination (R 2 ), the RMSE, and the rRMSE.

Impact of spatial resolution on density estimation

The impact of spatial resolution on image-based density estimation was studied using RGB images with degraded spatial resolutions as inputs to estimate GF and then density. The results were compared with those obtained with reflectance-based methods, which do not change with spatial resolution according to the linear mixing model of reflectance spectra [START_REF] Adams | Spectral mixture modeling: A new analysis of rock and soil types at the Viking Lander 1 site[END_REF].

Realistic low-resolution images were generated by successively degrading the spatial resolution by a factor of two according to the method proposed by [START_REF] Velumani | Estimates of maize plant density from UAV RGB images using Faster-RCNN detection model: impact of the spatial resolution[END_REF]: a

Gaussian filter with a sigma of 0.63 and a window size of 9 was first applied to the image, followed by motion blur with a kernel size of 3 and an angle of 45, and resizing to half of its height and half of its width. Further degradation of the spatial resolution was achieved by repeating the operations multiple times, resulting in ground sampling distances (GSDs) that were 2, 4, 8, and 16 times as large as the original size. Therefore, the average spatial resolutions of the original and generated image sets were 0.2, 0.4, 0.8, 1.6, and 3.2 mm.

For each spatial resolution, original or generated RGB images were used as inputs for the baseline method (section 2.

Relationship between the RGB-derived GF and Density

The relationship between RGB-derived GF (GF_rgb) and density was explored. For example, the data for barley in Avignon were plotted in Figure 4. For different view zenith angles and different growth stages, the linear relationships between GF and density were strong, with rRMSE values of density estimation not exceeding 20%. However, the slope strongly differed across growth stages and view zenith angles, i.e., it decreased from one-leaf to three-leaf growth stages, and from 0° to 45° view zenith angles.

The calibration relative RMSE values of density estimation from RGB-derived GF for all species, sites, and stages are shown in Table 2. Only the Avignon, Gardanne, and Nanjing sites had data for the three growth stages. For these three sites, the data showed strong relationships between GF and density, yet some differences could be observed: the rRMSE averaged over these three sites, the four species and the two view angles was smaller for Stage 2 (8%) and

Stage 3 (9%) compared to Stage 1 (13%), while 0° and 45° zenith angles had more similar rRMSE values (10% vs 9%, respectively). The data from Salin-de-Giraud, Greoux-les-Bains, and

Mauguio often had larger rRMSE values, showing weaker relationships between GF and density for these sites.

Ablation study on the effects of different factors

For spectra-GF and GF_rgb-density estimations for both view zenith angles, the performance strongly varied with factor combinations, with the rRMSE values ranging from 28% to 49%, and from 21% to 56% (Table 3). For spectra-GF, the best performance was obtained by calibrating one sub-model per site that include the four species and the three growth stages together (average rRMSE of 33% for 0° and 28% for 45°). For GF_rgb-density, the best performance was obtained by calibrating one sub-model for each site, each species, and each growth stage (average rRMSE of 23% for 0° and 21% for 45°). Not differentiating the species only slightly degraded the performances (average rRMSE of 28% for 0° and 24% for 45°), while the other factor combinations led to significantly worse GF_rgb-density estimation results. On the other hand, the spectra-density estimation performance was less variable concerning factor combination: the rRMSE value ranged from 30% to 42% for 0° and 45° view zenith angles. At 0°, the best performances were obtained by calibrating one sub-model per site (average rRMSE of 35%), while at 45°, it was better to calibrate a general model including six sites, four species and three growth stages (average rRMSE of 30%).

The effect of each factor can be shown by comparing the rRMSE before and after differentiating this factor (Table 3). For spectra-GF estimation, differentiating sites led to an obvious improvement in accuracy (row 7 versus row 8), but differentiating species (row 4) or stage (row 6) were not as effective. With finer differentiation of sub-models, the overall estimation accuracy degraded (row 1). For GF_rgb-density estimation, differentiating growth stage and site made a significant improvement to the estimation (row 5 versus row 8), and further differentiating species yielded the best accuracy (row 1). For spectra-density estimation, the rRMSE values were less variable, and the best factor combinations differed for 0° and 45°. For 0° observations, differentiating sites slightly improved the general model (row 7 versus row 8), while for 45° observations, calibrating a general model had the best accuracy (row 8). For spectra-density estimation, further differentiation led to a decrease in accuracy.

Overall, the estimation was more accurate at the 45° view zenith angle than at 0° in 22 of 24 cases in Table 3. 3.4 Band selection with AICc for Spectra-GF estimation Forward band selection was applied to Spectra-GF estimation using the sum of AICc values as a criterion, with each AICc value corresponding to one sub-model per site as recommended by the above ablation study (section 3.3). The AICc slightly decreased when adding one to four bands to the GPR model, and then increased more and more rapidly when adding more bands (Figure 5). The minimum sum of AICc values was reached at four bands for 0° view zenith angle, and three bands for 45° (Figure 5). These bands were 684, 759, 1128, and 1780 nm for 0° view zenith angle, and 419, 759, and 1548 nm for 45°. The site-specific GF estimations with selected bands got low RMSE values of 0.018 and 0.025 for 0° and 45°view zenith angles, respectively (Figure 6). However, due to the low GF values considered, these RMSE corresponded to moderate relative RMSE (rRMSE) values of 32% and 26% for 0° and 45° zenith angles, respectively. The comparison between selected bands and the five common bands showed subtle differences in rRMSE, i.e., 32% vs 33% for 0°, and 26% vs 28% for 45°.

Accuracy of specifically calibrated estimation for different sites, stages, and species

The density estimations were made based on the factor combinations chosen in section 3.3, and the band combination chosen in section 3.4. The 45° view zenith angle was chosen because it had a lower rRMSE than 0° with the chosen factor combination and band combination (Table 3).

Table 4 shows the average results obtained over ten replicated cross-validations, while Figure 7 shows scatter plots and residual plots obtained for one of these ten cross-validations. Note that data from all species have been included in the results (Table 4, Figure 7), but the different species were not marked to keep the results clear.

Overall, the direct reflectance-based estimation method (Spectra-Density) got similar results

to the indirect method (Spectra-GF-Density) with RMSE values close to 77 plants/m 2 (Table 4, Figure 7A,C). The baseline image-based estimation method (GF_rgb-Density) performed better, with an average RMSE value of 54 plants/m 2 (Table 4, Figure 7E). Note that estimated density values obtained with the direct method had a significantly smaller standard deviation (68 plants/m 2 ) than those obtained with the indirect (105 plants/m 2 ) and baseline methods (103 plants/m 2 ), both of which were comparable to the standard deviation of true density values ( 101plants/m 2 ). All three methods (direct, indirect, and baseline) tended to overestimate the density value when the true density value was low and to underestimate when the true density value was high. This trend was most evident in the direct method, less evident in the indirect method, and least evident in the baseline method (Figure 7B,D,F).

The estimation performance significantly differed across the different sites (Table 4). For example, the direct and indirect estimations led to smaller RMSE values between 38 and 61 plants/m 2 for Greoux-les-Bains and Mauguio, while they led to higher RMSE values between 92 and 102 plants/m ² for Gardanne, Salin-de-Giraud, and Nanjing. Generally, direct and indirect estimations were similar across sites. On the other hand, the RMSE values obtained with the baseline method were significantly lower than those obtained with reflectance-based methods for Avignon, Gardanne, and Nanjing, and similar for Salin-de-Giraud, Greoux-les-Bains, and Mauguio.

The differences among stages were checked through the comparison between pairs of RMSE values. For example, when considering the six sites and three methods, there were 7 out of 18 cases where RMSE could be computed for both Stage 1 and Stage 2 (Table 4). In 6 out of these 7 cases, the density estimation was more accurate at Stage 2 than at Stage 1. A similar paired comparison also showed that Stage 3 was better than Stage 1 in 8 out of 10 cases, while

Stage 2 was better than Stage 3 in 2 out of 3 cases (Table 4). The slope of the relationship between GF and density for barley in Avignon changed with view zenith angles and growth stages (Figure 4). Furthermore, the ablation study on the GF-Density relationship showed that, for a given view zenith angle, it was important to differentiate not only the stage but also the site, to significantly improve the density estimation accuracy based on GF (Table 3).

Especially at the 45° view zenith angle, the site factor had a large and complex effect on the GF-Density relationship. At least three sources of variation related to the site factor could be identified. First, for the same species and the same growth stage, the plant vigor could change according to soil and climate conditions, e.g., the air temperature or the soil type, thus affecting the plant architecture and GF values. For example, plots in Avignon had GF values approximately twice as large as those in Gardanne, even though both plots had the same species, the same growth stage, and similar density (Figure 10A,C). The cause could be the higher temperature in Avignon or differences in water availability during the experimental period (section 2.1). The second source of variation related to the site factor was the presence of weeds, which could artificially increase the estimated GF value. This effect was particularly important in sites with a large number of weeds such as Salin-de-Giraud (section 2.1, Figure 10B, F) because the GF values were very small at such an early stage (Figure 6). Finally, the third source of variation was the variability in soil roughness that could change the size and number of visible leaves at early stages, especially at the 45° view zenith angle. For example, plants were entirely visible in Avignon and Salin-de-Giraud where the soil surface was flat (Figure 10A, B), while only parts of the plants were visible in Nanjing where the soil surface was rougher (Figure 10D).

The growth stage factor also significantly affected the GF-Density relationship for a given view zenith angle, because this factor changed the size of plants. Later growth stages meant larger plants, and thus the GF values were larger while the density remained the same.

Compared to site and growth stage factors, differentiating the species factor had a more marginal yet positive effect on the GF-Density relationship (Table 3). As for growth stages, the different species could have different plant architectures, i.e., not only different leaf sizes but also different leaf orientations. For example, barley leaves were wider than wheat leaves in our experiments and that could lead to differences in GF-Density relationship.

The large diversity in growth stages, sites, and species in our dataset (Table 1) and the above results allow us to further discuss the results obtained by [START_REF] Gnädinger | Digital counts of maize plants by unmanned aerial vehicles (UAVs)[END_REF], [START_REF] Wilke | Assessment of plant density for barley and wheat using UAV multispectral imagery for high-throughput field phenotyping[END_REF] and [START_REF] Liu | Estimation of wheat plant density at early stages using high resolution imagery[END_REF]. First, the poor relationship between GF and density for maize observed by [START_REF] Gnädinger | Digital counts of maize plants by unmanned aerial vehicles (UAVs)[END_REF] was probably due to the nondifferentiation of three growth stages, four cultivars, and six cultural practices, all of which led to strongly different GF values for the same density. Second, our results confirm those of Wilke et al.

(2021), i.e., accurate plant density estimates can be obtained thanks to GF estimates when differentiating species and growth stages. However, our results further demonstrate the critical influence of the factor of site, which could not be observed by [START_REF] Wilke | Assessment of plant density for barley and wheat using UAV multispectral imagery for high-throughput field phenotyping[END_REF] since they only had one site and one year. Finally, our results are in agreement with those of [START_REF] Liu | Estimation of wheat plant density at early stages using high resolution imagery[END_REF], in which GF was one of the most important inputs to the plant density estimation model. Liu et al.

(2017) also emphasized the need for site-specific calibration models but did not separate wheat cultivars, probably because of the fewer differences observed among wheat cultivars as compared to among straw cereal species in our study (soft wheat, durum, barley, rye). Our study further demonstrates the importance of growth stage differentiation, since [START_REF] Liu | Estimation of wheat plant density at early stages using high resolution imagery[END_REF] only considered one stage. 4.2 GF is estimated more accurately using sitespecific models based on a few spectral bands acquired from the 45° view zenith angle

The spectra-GF relationship was analyzed in the ablation study (Table 3, "Spectra-GF" column). Differentiating only the site factor gave the best result, indicating that the site factor was the most important. Indeed, the site factor could change the canopy reflectance through different soil colors (Table 1, Figure 10) and reflectances, resulting in different canopy reflectances, even for the same GF. What made this change even stronger was the fact that the soil fraction was much larger than the vegetation fraction in the early stages. This explains the need for sitespecific models to estimate GF. On the other hand, the species and growth stage factors affect the spectral reflectance through a change in vegetation structure and more specifically, mainly through a change in GF. This fact could keep the spectra-GF relationship generally unchanged.

Furthermore, differentiating the species and growth stages in addition to the sites has a negative effect on the estimation (Table 3, Spectra-GF), probably because this led to too small training datasets and unstable GPR performance.

Using the 45° spectral observations for GF estimation generally performed better than using 0° observations (Table 3, Figure 6). One possible reason is the larger projection area of the plants in the 45° view zenith angle could help the small plants in the early stages to be more easily detected by spectrometer and camera [START_REF] Jin | Estimates of plant density of wheat crops at emergence from very low altitude UAV imagery[END_REF][START_REF] Liu | Estimation of wheat plant density at early stages using high resolution imagery[END_REF] (Figure 11).

The forward band selection method with AICc (Figure 5) showed the importance of choosing the right number of bands as inputs to the estimation model. When the number of bands was too low, the model was not able to properly separate the influences of the different factors causing variations in canopy reflectance, and thus not able to properly estimate GF. When the number of bands was too high, the high AICc values indicated that the model was less likely to reflect the true relationship [START_REF] Burnham | Multimodel inference: understanding AIC and BIC in model selection[END_REF]. Models with too many input variables could be more sensitive to the noise, and more prone to overfitting. The best number of chosen bands was three for the 45° view zenith angle and four for the 0° view zenith angle. In both cases, the model got a balance between low error and parsimony.

However, the interpretation of selected bands for 0° (684, 759, 1128, and 1780 nm) and 45° view zenith angles (419, 759, and 1548 nm) was difficult, even if bands in the red (684 nm) and near-infrared (759 and 1128 nm) domains are often used for vegetation remote sensing due to the strongly different responses of soil and vegetation in these spectral ranges. Despite the already large dataset collected, more data would thus be needed to confirm these band selections. Alternatively, using a common band set with five bands: B-G-R-RE-NIR led to slightly poorer performance than using selected optimal band sets Figure 6). On one hand, this result shows that forward band selection was effective in selecting an optimal band set for a specific dataset because this method yielded the optimal performance. On the other hand, it shows that the common band set was sufficient for practical use, as it performed similarly to the selected band sets.

Note that, despite the small RMSE of around 0.02 obtained for GF estimation, the rRMSE was moderate (between 26% and 33%, Figure 6) because the overall GF values were also small for the early stages. This explains the moderate plant density performances obtained with the indirect method (Table 4, Figure 7), which used estimated GF values as inputs to the GF-Density linear model. The direct (Spectra-Density) and indirect (Spectra-GF-Density) methods had similar error values, but they showed different features. The direct method was slightly more accurate (Figure 7, Table 4). Another advantage of the direct method was that the best results were obtained using a general model calibrated with all the data, which could be convenient for practical use. However, this result is counter-intuitive, since the reflectance-based methods should be based on GF as a proxy, and they should not be able to handle plants of different growth stages with only one model. Therefore, this result should be confirmed using more data. A notable feature of the direct method was that its standard error of estimated values was much less than that of the imagebased method. The direct method indeed overestimated plots with low-density values and underestimated plots with high-density values more severely than the other two methods (Figure 7). As a result, the direct method is more likely to fail for plots with extremely low or high densities.

For the indirect method, the estimation was more interpretable. On one hand, ablation studies could be made separately for the first step (Spectra-GF estimation) and the second step (GF-density estimation). The ablation studies suggested different ways to make local calibration for these two steps (Table 3), thus allowing us to make further analyses in sections 4.1 and 4.2.

On the other hand, the error of density estimation could be tracked in each step. For example, the Salin-de-Giraud site and the Nanjing site showed great differences in the source of error (Table 4).

For the Salin-de-Giraud site, a large proportion of error came from the second step (GF-density estimation), probably because of the detrimental influence of weeds (section 2.1). For the Nanjing site, the error from the second step (GF-Density) was small, indicating that the large error of the indirect method (Spectra-GF-Density) mainly came from the first step (Spectra-GF estimation). In the Nanjing site, the rugged soil and small seedlings could make GF values smaller, such that an error of 0.02 in GF estimation corresponded to a high relative error, leading to a high error in density estimation.

The accuracy of reflectance-based methods was lower than that of the baseline method based on submillimeter spatial resolution images (Table 4, Figure 7). In our study, reflectancebased methods got the best density estimation results at Stage 3, either with direct or indirect estimation, with an rRMSE value of around 28%. As a comparison, the baseline image-based method achieved rRMSE values from 9% to 24%, which were consistent with other studies based on submillimeter spatial resolution wheat images that reached relative errors between 9% and 17% [START_REF] Jin | Estimates of plant density of wheat crops at emergence from very low altitude UAV imagery[END_REF][START_REF] Liu | Estimation of wheat plant density at early stages using high resolution imagery[END_REF][START_REF] Wilke | Assessment of plant density for barley and wheat using UAV multispectral imagery for high-throughput field phenotyping[END_REF]. The superiority of image-based methods could be explained by the possibility to remove the detrimental influence of soil from the vegetation signal thanks to the high spatial resolution. Therefore, when submillimeter spatial resolution images are available, image-based estimation methods are recommended because they usually yield higher accuracy than reflectance-based methods.

In practice, density estimation could be done for one growth stage, instead of all three growth 2021), where the 3-leaf stage was better than earlier stages for wheat or barley density estimation using GF as a proxy. The indirect reflectance-based method and image-based method in our work were similar to the method of [START_REF] Wilke | Assessment of plant density for barley and wheat using UAV multispectral imagery for high-throughput field phenotyping[END_REF]. Conversely, contrasting perspectives can be found in the work of [START_REF] Jin | Estimates of plant density of wheat crops at emergence from very low altitude UAV imagery[END_REF][START_REF] Liu | A method to calculate the number of wheat seedlings in the 1st to the 3rd leaf growth stages[END_REF], where morphological analyses on binary vegetation images were applied, and where 1-leaf and 2-leaf stages were preferred. At an earlier stage (e.g., 1-leaf), plants have less overlap, which could facilitate morphological analysis, while at a later stage (e.g., 3-leaf), plants are larger so they are easier to detect. This may explain the difference in the best growth stages with different methods.

Our study emphasizes an important limitation of methods that exploit the relationship between density and GF: the estimation is affected by site, stage, and possibly species. A possible solution for this problem is to calibrate one model for each site, stage, and species. To avoid laborious manual counting that is necessary for calibration, image-based counting methods, e.g., based on deep learning detection algorithms [START_REF] Liu | High-Throughput Rice Density Estimation from Transplantation to Tillering Stages Using Deep Networks[END_REF][START_REF] Shubhra | DeepWheat: Estimating Phenotypic Traits from Crop Images with Deep Learning[END_REF][START_REF] Tseng | Rice Seedling Detection in UAV Images Using Transfer Learning and Machine Learning[END_REF], could be used because they have higher estimation accuracy, and their lower throughput would be enough to build a small dataset for calibration. However, density estimation approaches based on GF would still require substantial effort in data collection, not only to train each model with a sufficient number of samples, but also to capture new data at the right growth stage (or time window) so that these data can be used as input to the trained model.

4.4 Higher spectral resolution can somehow compensate for a lower spatial resolution to estimate plant density over large fields

Estimations based on high spatial-resolution images have better accuracy, but it is hard to collect high-resolution images with high throughput. Many of the studies on wheat and rice seedling density estimation used ground-based or UAV images that were acquired at no more than 20 m height [START_REF] Jin | Estimates of plant density of wheat crops at emergence from very low altitude UAV imagery[END_REF][START_REF] Liu | Estimation of wheat plant density at early stages using high resolution imagery[END_REF][START_REF] Liu | A method to calculate the number of wheat seedlings in the 1st to the 3rd leaf growth stages[END_REF][START_REF] Liu | High-Throughput Rice Density Estimation from Transplantation to Tillering Stages Using Deep Networks[END_REF][START_REF] Shubhra | DeepWheat: Estimating Phenotypic Traits from Crop Images with Deep Learning[END_REF][START_REF] Wilke | Assessment of plant density for barley and wheat using UAV multispectral imagery for high-throughput field phenotyping[END_REF]. Calculating with the experimental plan of Jin et al. ( 2017), a UAV taking photos at 10 meters height could cover 0.17 hectare in 1 hour. That could be an insufficient throughput when sampling large fields of several hectares. In this case, one solution would be to increase the altitude of the UAV, thus decreasing the image spatial resolution.

However, when the spatial resolution gets coarser, the accuracy of image-based methods decreases (Figure 9), while the accuracy of reflectance-based methods should not change. The decrease in performance observed for the image-based method is due to the increasing number of mixed soil/vegetation pixels (Figure 8) and is consistent with the results of [START_REF] Jin | Estimates of plant density of wheat crops at emergence from very low altitude UAV imagery[END_REF] who recommended using spatial resolutions lower than 0.4 mm. On the other hand, canopy reflectance should remain unchanged when degrading the spatial resolution if the canopy is spatially homogeneous, according to the linear mixing model of reflectance spectrum [START_REF] Adams | Spectral mixture modeling: A new analysis of rock and soil types at the Viking Lander 1 site[END_REF][START_REF] Ritter | Chapter 4 -Lattice Algebra Approach to Endmember Determination in Hyperspectral Imagery[END_REF], thus the performance of reflectance-based methods should be stable for different spatial resolutions. In our experiment, this performance was obtained with a spatial resolution of about 736 mm, corresponding to the length of the side of a square with the same area as the spectrometer's field of view (Figure 1). According to Figure 9, the performances of reflectance-based methods would exceed those of the baseline image-based method when the GSD is larger than 1 to 2 mm for one-leaf to three-leaf stages, respectively. According to the experimental settings of [START_REF] Roth | Repeated multiview imaging for estimating seedling tiller counts of wheat genotypes using drones[END_REF], [START_REF] Wilke | Assessment of plant density for barley and wheat using UAV multispectral imagery for high-throughput field phenotyping[END_REF][START_REF] Jin | Estimates of plant density of wheat crops at emergence from very low altitude UAV imagery[END_REF], a spatial resolution of 2 mm can be obtained by flying the UAV at about 10 m to 30 m height for focal lens length ranging from 20 mm to 60 mm, respectively. This result would thus support the use of reflectance-based estimation for UAV altitudes of more than 10 m to 30 m above the ground depending on the optics, gaining higher throughputs in density estimation. Further studies are needed to test this method on larger fields and UAV reflectance data.

Note that the SegVeg model [START_REF] Madec | VegAnn, Vegetation Annotation of multi-crop RGB images acquired under diverse conditions for segmentation[END_REF][START_REF] Serouart | SegVeg: Segmenting RGB images into green and senescent vegetation by combining deep and shallow methods[END_REF] used to identify vegetation pixels was trained on images with spatial resolution ranging from 0.3 to 2 mm (section 2.4), so the simulated GSD of 3.2 mm in Figure 9 was slightly outside of the applicable range of SegVeg. While this may cause some uncertainties in the determination of the spatial resolution where both reflectance-based and image-based methods perform the same (only for Stage 2 and

Stage 3), this will not change the general trend already observed from 0.2 to 1.6 mm, i.e., that the performance of the image-based method decreases with increasing GSD.

By exploiting the spectral information, and especially a combination of NIR and visible bands where the responses of soil and vegetation are strongly different, it becomes possible to somehow compensate for lower spatial information [START_REF] Jacquemoud | PROSPECT+SAIL models: A review of use for vegetation characterization[END_REF]. This is consistent with the results of [START_REF] Wilke | Assessment of plant density for barley and wheat using UAV multispectral imagery for high-throughput field phenotyping[END_REF], who demonstrated that a multispectral camera could provide better plant density estimates based on image thresholding than an RGB camera, despite the lower spatial resolution of multispectral images (7 mm vs. 2 mm). In addition, our work shows that, for even coarser spatial resolutions for which soil and vegetation cannot be discriminated, canopy reflectance can be a reasonable proxy of GF and plant density.

Conclusion

In this study, the straw cereal plant density at early stages was estimated from spectral reflectance measured at the nadir or from 45° view zenith angle, with (indirect method) or without (direct method) using GF as a proxy. The results were compared to those obtained with a popular image-based method, using a large and diverse dataset including different sites, species, and growth stages. According to the ablation study performed with five common spectral bands, the spectra-GF estimation (first step of indirect estimation) was site-specific and stage-specific; the GF-density estimation (second step of indirect estimation) was site-specific, stage-specific, and species-specific; the spectra-density estimation (direct estimation) was not specific at all. Using a 45° view zenith angle showed slightly better performance on average so 45° was chosen. Using only three spectral bands selected by minimizing the AICc criteria, the direct and indirect estimations had similar relative errors of around 30% (RMSE = 76 plants/m 2 ), while better performance was obtained with the image-based method when using submillimeter image spatial resolutions (RMSE = 54 plants/m 2 ). However, a study on downsampled images showed that reflectance-based estimation outperformed image-based estimation when the GSD of images was larger than a threshold between 1 to 2 mm depending on growth stages, thus reflectancebased estimation has a better potential for high-throughput estimation of straw cereal plant density.

The proposed indirect plant density estimation method could be applied to UAV multispectral images to get high-throughput density estimates. This potential was supported by two reasons.

First, the commonly used band set of commercial multispectral cameras (B, G, R, RE, NIR) performed almost as well as chosen bands from the spectrometer in spectra-GF estimation,

showing there is enough information in this band set for the density estimation task. Second, this reflectance-based method is robust to a degradation in spatial resolution, which means it allows higher flying altitudes and higher throughput, while keeping the same accuracy.

A general model of direct density estimation may be possible, but this somewhat unexpected result did not explain its capability in handling different growth stages. This will need to be confirmed with a larger dataset with different varieties, different sites, and different growth stages as factors.

Figure 1 :

 1 Figure 1: Acquisition for 0° and 45° spectral measurements. (A) The real scene at the Avignon site. The 45° measurement view direction is perpendicular to the row. (B) The geometry design of the

Figure 2 :

 2 Figure 2: Workflow of direct, indirect, and baseline methods. GPR means Gaussian process regression. LR means linear regression.

  3.1), and the density estimation RMSE was calculated. These RMSE values were then compared with the RMSE values obtained with methods based on spectral reflectance. Note that the SegVeg segmentation model of[START_REF] Serouart | SegVeg: Segmenting RGB images into green and senescent vegetation by combining deep and shallow methods[END_REF] was trained on images of 0.3 to 2 mm GSD, thus potentially causing some uncertainties in the segmentation results obtained at 3.2 mm spatial resolution.

Figure 3

 3 Figure 3 Stacked histograms of density (A, B, C, G, H, I) and green fraction values observed at 45°

Figure 4 :

 4 Figure 4: Relationships between RGB-derived GF and barley plant density at the Avignon site for

Figure 5 :

 5 Figure 5: Results of forward band selection from 140 bands, showing AIC c as a function of the number of input bands included in the GPR model for (A) 0° and (B) 45° view zenith angles. The models were calibrated for each site. For 0° view zenith angle, bands selected at minimum AIC c were 684, 759, 1128, and 1780 nm (C). For the 45° view zenith angle, bands selected at minimum AIC c were 419, 759, and 1548 nm (D). Numbers of bands greater than 25 were not shown because the AIC c value could be invalid in that case.

Figure 6 :

 6 Figure 6: Spectra-GF estimation with five common bands or with optimal band sets for (A) 0° and (B) 45° view zenith angles. One estimation was made for each site and the overall results are shown.

Figure 7 :

 7 Figure 7: Scatter plots for (A) direct and (C) indirect reflectance-based estimations, and (E) baseline image-based estimation; and the residual plots of (B) direct, (D) indirect, and (F) baseline estimations.The estimations were obtained at the 45° view zenith angle and based on the results of the ablation study (section 3.3) and band selection (section 3.4). In scatter plots (A, C, E), growth stages and sites are respectively shown using markers and colors, while species are not differentiated.The sample size for each category (n) is also provided. In residual plots (B, D, F), each blue point shows the residual for a data point, and each orange point shows the mean residual in its neighboring area of 50 plants/m 2 width. The estimations correspond to one of the ten replicated cross-validations.

Figure 8

 8 Figure 8 Patches of images with different spatial resolutions as input (RGB) and output (binary) of the SegVeg model in the first step of the baseline method (section 2.3.1). The GSD values are marked below each column of images.

Figure 9 : 1

 91 Figure 9: Impact of RGB image spatial resolution on the density estimation RMSE obtained with the baseline image-based method (dotted blue line), for the three growth stages (Stage 1: (A), (D); Stage 2: (B), (E); Stage 3: (C), (F)) and the two view zenith angles (0° : (A), (B), (C); 45°: (D), (E), (F)).This baseline method was calibrated for each site, each stage, and each species. For comparison, the RMSE obtained with the direct (solid blue line) and indirect (solid orange line) reflectance-based methods applied to 736 mm spatial resolution spectrometer measurements from section 3.5 were also shown. RMSE values were averages obtained over ten replicated cross-validations. Stage 1: from 1 to 1.6 leaves; Stage 2: from 1.7 to 2.3 leaves; Stage 3: from 2.4 to 3 leaves.

Figure 10 :

 10 Figure 10: Examples of RGB images acquired at 45° view zenith angle for different sites: (A) Avignon, (B) Salin-de-Giraud, (C) Gardanne, (D) Nanjing, and the corresponding binary images: (E), (F), (G), (H). Plants in these plots were soft wheat, being in similar growth stages (1.8~2 leaves), with similar plant density (113~133 plants/m 2 ).

Figure 11 (

 11 Figure 11 (A) Reflectance spectra acquired from 0°and 45° view zenith angles, and

stages (Stage 1 : 1 .

 11 0~1.6 leaves; Stage 2: 1.7~2.3 leaves; Stage 3: 2.4~3.0 leaves). In our study, the comparison between pairs of RMSE in section 3.5 shows that Stage 2 and Stage 3 were better than Stage 1 in density estimation with the 45° zenith angle, with either the direct or indirect reflectance-based method. A supportive perspective can be found in the work of Wilke et al. (

Table 1

 1 Site, date, and species of the experiment. "2.5 leaves" means that the third leaf was not fully expanded, and the third leaf length was about 50%

	151
	152			of the second leaf length.
	Site and sowing date Soil type	Soil color	Number	Species	Density values in	Number of leaves
			(dry/wet soil	of plots		mean±std	(date of
			surface)			(#plants/m 2 )	measurements)
	Avignon, France.	Clayey,	White (dry),	18	soft wheat,	307 ± 130	1.0 (2021/Oct/2),
	2021/9/23	calcareous,	brown (wet)		durum wheat,		2.0 (2021/Oct/8),
		fluvisol			barley		3.0 (2021/Oct/15)
	Nanjing, China.	Sandy	Yellow (dry),	9	soft wheat,	171 ± 91	1.0 (2021/Nov/22),
	2021/11/5		brown (wet)		rye,		2.0 (2021/Dec/3),
					barley		3.0 (2021/Dec/17)
	Greoux-les-Bains,	Clayey,	Light Taupe	32	soft wheat,	239 ± 27	1.5 (2021/Nov/18),
	France.	alluvium	(wet)		durum wheat,		2.5 (2021/Nov/29)
	2021/10/28				barley
	Gardanne, France.	Silty clayey,	Red	9	soft wheat,	300 ± 120	1.5 (2022/Jan/3),
	2021/11/19	calcareous,	(dry and wet)		durum wheat,		1.8 (2022/Jan/11),
		alluvium			barley		3.0 (2022/Feb/9)
	Salin-de-Giraud,	Sandy,	Dark grey	28	soft wheat,	212 ± 103	1.0 (2021/Dec/20)
	France.	calcareous	(wet)		durum wheat,		1.8 (2022/Jan/4)
	2021/11/22				barley
	Mauguio, France.	Calcareous,	Light Taupe	13	soft wheat	296 ± 66	3.0 (2022/Jan/7)
	2021/11/19	fluvisol	(dry)

Table 2

 2 Calibration relative RMSE values (in %) of density estimation with RGB-derived GF for each site, stage, and species. "/" means there were no 377 data collected for the case. Stage 1: from 1 to 1.6 leaves; Stage 2: from 1.7 to 2.3 leaves; Stage 3: from 2.4 to 3 leaves.

	378
			Number of	rRMSE for 0° view zenith angle (%)	Number of	rRMSE for 45° view zenith angle (%)
	Specie s	Site	samples for 0° and all	Stage 1	Stage 2	Stage 3	samples for 45° and all	Stage 1	Stage 2	Stage 3
			stages				stages
		Avignon	19	19	10	12	10	20	10	13
		Salin-de-Giraud	18	43	16	/	18	46	30	/
	Soft	Gardanne	12	8	14	14	12	14	7	7
	wheat	Greoux-les-Bains	22	22	/	23	22	23	/	26
		Mauguio	35	27	/	17	34	19	/	21
		Nanjing	15	18	12	14	12	6	7	13
		Avignon	19	14	12	8	12	15	10	10
	Durum	Salin-de-Giraud	7	/	38	/	7	/	41	/
	wheat	Gardanne	12	3	5	7	12	7	2	5
		Greoux-les-Bains	19	30	/	15	19	17	/	12
		Avignon	19	20	11	12	12	17	7	9
		Salin-de-Giraud	3	/	20	/	3	/	22	/
	Barley	Gardanne	12	3	11	5	12	17	6	4
		Greoux-les-Bains	20	15	/	13	19	17	/	13
		Nanjing	15	8	6	4	14	9	5	7
	Rye	Nanjing	15	17	5	6	13	10	8	6

Table 3

 3

	Factor Combinations	Spectra→GF	GF_rgb→Density	Spectra→Density
	Diff-Spc	Diff-Stg	Diff-Site	0°	45°	0°	45°	0°	45°
	√	√	√	47*	42*	23*	21*	42*	41*
	√	√	×	47*	46*	35*	38*	42*	38*
	√	×	√	43*	34*	47*	38*	40*	38*

rRMSE values (in %) obtained for the different factor combinations, the three estimation steps (Spectra-GF, GF_rgb-density, and Spectra-Density, see section 2.3.2), and the 0° and 45° view zenith angles. "Diff-Spc" means "differentiate species", "Diff-Stg" means "differentiate growth stages", and "Diff-Site" means "differentiate sites". For each column, the best average rRMSE is in bold. For the sake of simplicity, only five common bands were used here as inputs to the GPR models. The asterisks (*) denote the values obtained using leave-one-out cross-validation (section 2.3.2). The difference in average values of GF and density used to compute relative RMSE within each column was ignorable.

Table 4

 4 RMSE values (in plants/m2 ) obtained for plant density estimation at 45° view zenith angle with the direct and indirect reflectance-based methods, and the baseline image-based method. The factor combinations and band combinations were chosen as described in section 3.3 and 3.4. RMSE values were computed per growth stage and per site, and by grouping all growth stages and/or all sites. All species available were used to compute each RMSE value. The symbol "/" in the cell means the number of samples is not sufficient for the estimation.

	The RMSE values were the mean values calculated over ten replicated cross-validations. Stage 1: from 1 to 1.6 leaves; Stage 2: from 1.7 to 2.3 leaves; Stage
		3: from 2.4 to 3 leaves. Direct method (Spectra-Density)		Indirect method (Spectra-GF-Density)	Baseline method (GF_rgb-Density)
		Stage 1	Stage 2	Stage 3	All	Stage 1	Stage 2	Stage 3	All	Stage 1	Stage 2	Stage 3	All
				stages				stages				stages
	Avignon	/	52	88	71	/	70	88	79	/	33	38	36
	Salin-de-Giraud	137	66	/	100	134	74	/	102	137	62	/	99
	Gardanne	119	76	94	97	125	122	67	100	70	25	22	40
	Greoux-les-Bains	58	/	26	45	38	/	38	38	51	/	50	51
	Mauguio	52	/	48	51	54	/	71	61	61	/	60	60
	Nanjing	101	103	87	96	/	89	94	92	/	16	25	21
	All sites	84	77	71	77	79	86	71	78	75	39	41	54

Acknowledgements

This work was supported by the Agence Nationale de la Recherche, programme Investissements d'avenir, ANR-11-INBS-0012 (Phenome), the Young Scientists Fund of the National Natural Science Foundation of China (No. 42201437), the Young Scientists Fund of the Natural Science Foundation of Jiangsu Province, China (No. BK20210411), the National Key R&D Program of China (No. 2021YFD2000105, No. 2022YFE0116200), and the Project of Seed Industry Revitalization in Jiangsu Province, China (JBGS[2021]007). One author of this article, Tiancheng YANG, was funded by the Chinese Scholarship Council (Number 202006850003). We

Declaration of Generative AI and AI-assisted technologies in the writing process

During the preparation of this work, the authors used ChatGPT, DeepL, and Grammarly in order to improve the grammar. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.


Maxime Derex

Robert Boyd

Technical reasoning alone does not take humans this far

While we see much utility in Osiurak and Reynaud's in-depth discussion on the role of what they term technical reasoning in cumulative culture, we argue that they neglect the time and energy costs that individuals would have to face to acquire skills in absence of specific sociocognitive abilities.

MAIN TEXT

) never claimed that reasoning plays no role in CC. Rather they argued that the improvement of culturally evolving technology is not necessarily tied to individuals' level of understanding. Indeed, Derex et al.'s experiment shows that, over successive overlapping generations, participants produce increasingly efficient solutions despite exhibiting no improvement in causal understanding. This does not mean that causal reasoning cannot play a role in that process, but it shows that increases in efficiency are not necessarily accompanied (or even powered) by changes in individuals' understanding. These results illustrate the effects of the selective retention of beneficial modifications across generations and stress the roles of social learning and population-level processes in the emergence of adaptive cultural traits. That is not to say that the ability to reason has nothing to do with CC. Asking whether CC could occur in absence of specific reasoning abilities is an entirely different question. Osiurak and Reynaud's proposal that CC necessarily depends on species-specific ability to technically reason about phenomena is both timely and welcome. Unfortunately, the evidence and arguments that the authors bring to bear in support of their hypothesis are weak.

First, Osiurak and Reynaud overestimate the ability of individuals to extract relevant information by observing artefacts alone. They base their claims on a few experimental studies that showed that even minimal social learning mechanisms (such as reverse engineering) can result in cumulative improvements. According to them, "if signs of [CC] in reverse-engineering conditions are observed in humans, then this suggests that human nonsocial cognitive skills are sufficient for the emergence of [CC]". There are several problems with this argument. First, experimental investigations of CC rely on relatively simple tasks that participants can solve in a short period of time (Miton & Charbonneau, 2018). As a result, the amount of information that individuals can typically extract from the observation of such simple experimental artefacts is unrealistically high and whether technical reasoning skills allow individuals to infer substantial amounts of missing information about more ecologically valid artefacts remains to be demonstrated. Moreover, Osiurak and Reynaud neglect the fact that experimental settings typically allow sustained and undisturbed observations maximizing the effectiveness of reverse engineering. Even if we assume that technical reasoning skills allow individuals to infer missing information, the usefulness of such capabilities in natural settings might be highly limited in absence of specific sociocognitive abilities that give individuals appropriate access to cultural artefacts. Finally, we would like to point out that investigating information acquisition requires properly controlling for knowledge previously acquired by social learning. Indeed, Osiurak and Reynaud's argument that information extraction is mediated by individuals' level of expertise suggests that they conflate information that was acquired during a specific learning event with information that was acquired prior to this learning event. To take Osiurak and Reynaud's own example, showing that physics graduate students listening to Einstein retain more about the theory of relativity than individuals with no knowledge of physics might say less about what individuals actually learnt than what they previously knew.

Another problem with Osiurak and Reynaud's argument is their claim that much can be learned by reverse engineering if learners can alternate between periods of social and individual learning and that socio-cognitive skills only boost CC. The fact that people can learn much through trial-and-error learning does not imply that socio-cognitive skills are unnecessary to CC. Even if we assume that individual learners could build a traditional Inuit kayak from a pile of driftwood and seal skins without learning from others (which is unlikely to say the least) that does not mean that they will do so in absence of appropriate social support. More effective social learning strategies do more than just change the rate of CC. Individuals constantly face intense trade-offs and so must allocate their time and energy strategically. When learning costs are too high, individuals might not be able to afford to acquire complex skills by themselves. Experiments with children, for instance, show that they have difficulty innovating even simple tools even though they manufacture them easily after being exposed to social demonstrations (Beck et al. 2011). Moreover, being able to acquire a few skills through a combination of observational learning and individual leaning is not close to good enough. In the Arctic, kayaks are only useful if individuals can also develop warm clothes, harpoons and all other tools that their survival depends upon. Without appropriate socio-cognitive capabilities, the acquisition of these skills would require an investment in terms of time and energy that is way beyond what individuals can afford.

Finally, we would like to point out that experimental research has shown that toddlers are more likely to infer causal connections when sequences of events are the result of human actions rather than when they occur "naturally" without involving human interventions (Meltzoff et al. 2012). This means that reasoning skills are likely to at least partially result from socio-cognitive abilities that increase individuals' opportunities to witness valuable events. Furthermore, it suggests that hypotheses based on the role of reasoning skills should not be considered as alternative to those based on socio-cognitive

skills. A more fertile approach may be to study how both types of abilities reinforce each other.