Quality of forest plantations using aerial images and computer vision techniques 1

Geotechnology has provided several tools that allow the spatial and temporal variability of soils and plants to be investigated, leading to the consolidation of Precision Agriculture. The great challenge for studies using sensors mounted aboard Remotely Piloted Aircraft (RPA) lies in interpreting the high-dimensional data, since most sensors do not measure the biometric parameters of a plant directly. Therefore, the aim of the present study was to develop a methodology for using digital images (obtained by means of an airborne RGB sensor mounted aboard an RPA) in the quality control of forest plantations, specifically Eucalyptus (Eucalyptus ssp.), planted in a commercial area. A Phantom 4 Pro multirotor RPA was used, equipped with a 20 Megapixel RGB sensor, acquiring images with 80% and 60% longitudinal and lateral overlap, respectively. From the generated orthomosaic, a Test Area was outlined to be used in developing the processing routine based on computer vision techniques. In general, the proposed methodology maps the individual location of each plant in the orthomosaic, resulting in a mesh that allows the automatic generation of report maps of various silvicultural variables, such as plant count, planting failures, and spacing between rows and plants. In addition to high computer performance, with real-time processing, the methodology was highly accurate in correctly identifying more than 93% of plants in an area of more than 3,000 plants.


INTRODUCTION
Basic management principles that consider the variability of soils and plants in space and time can assist rural managers in decision-making in various agricultural production systems (ARTUR et al., 2014). This variability, when ignored, increases the inefficiency of using different inputs, reduces productivity and, as a result, raises the cost of production (BERNARDI et al., 2014).
The digital images used in agriculture can be obtained from several platforms, specifically, at three levels of acquisition: (a) orbital (Satellites), (b) aerial [Manned Aircraft (MA) and Remotely Piloted Aircraft (RPA)], and (c) terrestrial (portable sensors) (SHIRMASUCHI et al., 2014). These authors highlight various products derived from such images: estimating biomass, assessing water stress, the severity of disease, plant identification and count, and planting failures, among others.
During collection of the aero-photogrammetric data, some local variations should be considered, such as the relief (MANCONI et al., 2019) and lighting (SHI et al., 2016). It is recommended that the operation be planned and executed within the flight window (true solar noon, with an interval of plus or minus two hours), always taking into account the sensor settings (ISO, exposure time and aperture), as noted by Wang et al. (2019).
The great challenge for studies using airborne sensors (in RPA) lies in interpreting the high-dimensional data, since most sensors do not measure the physical, morphological or physiological parameters of plants directly (DUARTE;SILVA;TEODORO, 2018). There are many areas of research that need to be more detailed, in particular monitoring the models (until now, largely manually evaluated) and production forecasting (BALLESTEROS et al., 2014;BORGOGNO-MONDINO et al., 2018;KOH et al., 2019;POBLETE-ECHEVERRÍA et al., 2017). Ruza et al. (2017), state that, despite the advances in aero-photogrammetric activities employing RPA in the forestry sector, the production of applied and automated technical products has not progressed as quickly and still requires much further research. Özcan et al. (2017), and Park et al. (2017), also point out the importance of new studies that explore and present more feasible methods using images from RPA. For Yao, Qin and Chen (2019), with the development of the RPA, together with the onboard sensors, digital image processing, at its most simple, which consists of segmenting and classifying images in two and/or three dimensions, has also seen great advances, making it possible to detect and track objects on more precise scales (temporal and spatial). The above authors reaffirm that this has only been possible using techniques of computer vision and machine learning, affording significant gains in automated data analysis.
As such, the development and validation of methodologies that make it possible to extract and visualize the spatial variability of various plant attributes quickly and accurately (remote and automated) will contribute to the monitoring and quality control of agricultural and forest systems (FAN et al., 2018;GUERRA-HERNÁNDE et al., 2017).
Encouraged by these challenges, the aim of this study was to develop a methodology (using automated routines) for using aerial images (obtained by means of sensors mounted aboard RPA) in the quality control of forest plantations, specifically Eucalyptus (Eucalyptus ssp.), using computer vision techniques.

Study area
A commercial plantation (39.4812 ha) of Eucalyptus (Eucalyptus urophylla x E. grandis) was chosen, with a spacing of 3 m x 3 m, 180 days after planting, located in the district of Juramento, Minas Gerais. In order to develop the computer routine, a rectangle (Test Area, 3.2768 ha) was selected, centered on UTM coordinates E 663.461, N 8.131.368 and h = 1210 m (SIRGAS 2000, 23-k) ( Figure 1).

Collection of the aero-photogrammetric data
The data were collected on 11/23/2018, by means of a Phantom 4 Pro multi-rotor RPA equipped with a 20 Megapixel RGB sensor, 13.2 mm x 8.8 mm, and a focal length of 9.1561 mm. The flight plan was prepared using the DroneDeploy ® software, to construct the orthomosaic with a longitudinal and lateral overlap of 80% and 60%, respectively. The altitude was programmed for 250 m, with a maximum speed of 15 m.s -1 and a Ground Sample Distance (GSD) of 7.0 cm (8 bits radiometric resolution). It should be noted that the flight was planned to cover an area of 396 m x 997 m, and that no control points were used. Finally, the mission was carried out at 11:00 am (with the aim of minimizing the effects of shading), and resulted in the collection of 62 images.
After the flight, the images were imported into the Agisoft Metashape ® software for geoprocessing and to generate the orthomosaic. Cropping of the Test Area and photointerpretation were carried out using the QGIS 3.14 software.
Quality of forest plantations through aerial images and computer vision techniques

Implementation of the automated routine
The routine was implemented using the Spyder development environment in the Python programing language, making use of the OpenCV, Numpy, Pandas, Math, Matplolib and Scipy libraries. Basically, the structure of the algorithm was divided into two main phases. The first phase consisted of opening the image of the Test Area and segmenting the Soil and Vegetation by means of the Vegetation Index (VI), as well as identifying the center of each plant canopy for counting. The second phase consisted in building a triangular mesh using the vertices obtained as the center of the canopy of each identified plant and constructing the corresponding occupation diagram. From the mesh and diagram, the following Report Maps (RM) were generated: a) area of influence of each plant; b) identification of planting failures, rows and the area between rows; c) boxplot of the spacing between plants and rows; d) spacing distribution between plants and between rows. Each RM was accompanied by descriptive statistics and a Euclidean area (ha) of the selected polygon. The workflow for implementing the proposed routine for the quality control of forest plantations is summarized in Figure 2.
In the first phase of the routine, once the image was imported (Test Area) and after photointerpretation, the number of classes, or contrast patterns, was defined based on variations in the predominant characteristics of the ground cover. The number of samples (template) per class was then defined. In this study, tests were carried out using two classes: 1) plants in soil with ground cover (SC), predominant in the study area; 2) plants in exposed soil (SE). Both were combined using up to three samples. It was sought to diversify the physical characteristics (color, texture, canopy size and shading) of each sample, with the aim of representing the variability.
Six tests were carried out, in the following combinations: 1) one sample in SC; 2) two samples in SC; 3) three samples in SC; 4) three samples in SC + one sample in SE; 5) three samples in SC + two samples in SE; 6) three samples in SC + three samples in SE ( Figure 3).
The capture and automatic display of each template was made possible through the mouse-event control function, the user just having to position the mouse cursor on the desired plant and, with a double-click followed by enter, it is selected and the template displayed.
The RGB bands, the Test Area and the template(s) are then split. This procedure is necessary to apply the Gaussian filter to remove noise from the image and, consequently, estimate the VI. In order to segment the regions of soil and plants, a threshold with a value of 0 on the VI response scale was defined where lower values were considered as soil and other values as plants (Eucalyptus), this result was taken as the input for the next step of the first phase, template matching. It should be noted that the VI was estimated using the VARIgreen metric, proposed by Gitelson et al. (2002), which uses only the visible spectrum (RGB) channels.
Template matching, a crucial technique in this study, was applied using the OpenCV library. This technique is widely used in computer vision to locate areas (or regions) of an image that correspond (are similar) to a model image (template) (BRADSKI, 2000). A mask (template (T)) runs over an image (I) comparing each position [x, y] by the process of convolution. The result of the comparison is stored in a matrix, R[x, y] ( Figure 4).
In this study, (I) and (T) correspond respectively to the VI from the Test Area and the sample(s). It should  be noted that each of the six samples (templates) was processed separately during template matching, and then, when combined, the final value for R was obtained from the arithmetic mean.
The technique of template matching, implemented in the OpenCV library, offers several metrics to calculate the response; in this study, 'TM-CCOEFF-NORMED' was adopted as the metric, given by the mathematical expression in Equation 1.  x, y]. The result of this procedure will be a matrix with values in [x, y] which range from -1 to +1, i.e. low to high correlation, respectively. Clusters of pixels that presented values greater than 0.3 were marked as plant canopies, where the maximum values were considered as the center of each canopy. The first phase ended after locating [x, y] each pixel representing the center of the canopy and outputting the map identifying the eucalyptus plants.
The second phase started with the coordinates [x, y] of each pixel from the previous phase, i.e. a set of vertices with the location of each plant. From this set of vertices, a triangular mesh and the occupation diagram were constructed, obtaining, respectively, an edge mesh of the triangles (which allowed the spacing between plants and rows to be estimated) and the area of influence (diagram) of each plant. To identify planting failures, any edges outside the normal spacing between plants were identified. Finally, the following RM were obtained: a) area of influence of each plant; b) identification of planting failures, rows, and the area between rows; c) Descriptive Statistics.
For computer-processing time, it is relevant to note that, using a laptop with the Linux operating system (Ubuntu 18.04), Intel ® Core (TM) i7-4500U CPU (4 core, 1.80GHz, 4096 KB L2 cache), 16 GB RAM, SSD storage, and GeForce GT 740M GPU, it took 96 seconds and 25% of the RAM to analyze the specified Test Area (using the above routine).

Data analysis
The performance of the implemented routine was evaluated by comparing the number of identified plants and those counted manually (photointerpretation of the Test Area). Three performance measurement metrics (Equation 2, 3 and 4) were then applied, as per Fan et al. (2018). (2) (3) The performance measurements use the following counts: (TP) true positive, or the number of correctly identified eucalyptus plants; (TN) true negative, or the number of correctly identified planting failures; (FP) false positive, or the number of points incorrectly identified as eucalyptus, but being weeds and/or shade; and (FN) false negative, or the number of unidentified eucalyptus plants.
The sensitivity (Sb) reflects the ability of the algorithm to detect eucalyptus plants, while specificity (Sp) is a measure of the efficiency of the algorithm in identifying 'non-eucalyptus' plants. Overall accuracy (Ac) is a global measure of the performance of the proposed method.
In order to increase the rigor of evaluating the algorithm, two more evaluation metrics were used in this study. The first, producer accuracy (Pacc), is recommended by Lavrač, Flach and Zupan (1999) and Moranduzzo and Melgani (2013), and represents the percentage of correctly identified eucalyptus plants, Equation 5.
In measuring Pacc, N indicates the actual number of eucalyptus plants in the Test Area. The second measurement, recommended by Armstrong and Collopy (1992) and Stine et al. (2004), is the relative error (Er), by which the performance of the proposed method was determined. For this measurement, Np = TP + FP is used, which is the number of eucalyptus plants detected by the automated routine, Equation 6.

RESULT AND DISCUSSION
The study presented a highly relevant characteristic, applicability, considering that the routine was developed and tested under practical conditions for the crop under evaluation (commercial plantation). In terms of obtaining aero-photogrammetric data, the RGB orthomosaic, with a GSD of 7 cm, allowed the objects under analysis, plants and soil, to be discriminated. In images with high spatial resolution, planting and soil patterns are detectable, affording great potential for their discrimination and characterization (DELENNE et al., 2010).
With the adopted flight characteristics, it was possible to cover 39.4812 ha in approximately nine minutes. For the model of RPA used, this represents a minimum operational efficiency of 80 ha battery -1 . The time taken to collect such information is also of fundamental relevance to the forest survival inventory, since, when necessary, replanting is recommended 15 to 30 days after the original planting. As this interval is extended, the smaller the individual volume of each tree, and, consequently, the smaller its contribution to the final total volume of the forest (OLVIEIRA et al., 2014). For extensive forest plantations, such information is of great relevance.
For the developed routine, it is important to note that the user only has one moment of manual interference, which is the choice of sample(s) (template) that will be used as input for template matching. In this step, when necessary, the user should take into account the diversity of the characteristics (color, texture, shape and size) of the template, as different templates may contribute to the final result of matching.
Among the proposed (tested) combinations, there was, in general, no significant difference in performance metrics, especially for combinations of two classes. This confirms that the threshold value of zero, defined to segment the Plants and Soil using the VI, was efficient ( Figure 5). Calculation of the VI, followed by a process of noise removal (application of a Gaussian filter), made it possible to segment the plants and soil, increase gain in the correlation process of template matching and, as such, contributed to the identification of greater quantities of correctly identified eucalyptus plants (TP), and the number of correctly identified planting failures (TN), as well as reducing the number of points incorrectly identified as eucalyptus (FP), and the number of unidentified eucalyptus plants (FN).
The calculated VI (from the RGB channels) does not require the use of the infrared channel. In practice, this represents a saving for the user, as a second sensor is not needed (modified RGB, multispectral or even hyperspectral). Although the RGB sensor does not have the ability to decipher narrow spectral characteristics, as it is only three spectral bands in width, it gives good results in obtaining references for plant growth parameters (shape, color and canopy size) and in identifying weeds (BALLESTEROS et al., 2014;BURKART et al., 2018).
The image (R), resulting from the Test Area correlation (I), in this case the VI-RGB Map and the template(s) (T), was submitted to a process to identify the maximum values in each group of pixels (canopy), with the highest value (center of the canopy) being taken, (Figure 6).
Each group (representing the canopy of each plant) comprises a minimum value for the magnitude of the correlation, known as the cut-off threshold. Definition of this threshold represents the minimum acceptable magnitude, defined by the user, to correlate (I) and (T). As shown by Koh et al. (2019), low cut-off values can result in overestimation when identifying and counting the plants, whereas underestimation may result from high values.
The process of defining the ideal value for the cut-off threshold is iterative and often exhaustive. It is part of the calibration process of the computer routine to define the value that allows the greatest level of correspondence with the template to be obtained, i.e. the greatest overall accuracy of the algorithm. Koh et al. (2019), used template matching to identify and count saffron (Carthamus tinctorius) in the early stages of development, using images obtained with RPA, defining 0.5 as the value of the cut-off threshold. In the present study, a threshold of 0.3 was used.
Once the cut-off threshold was defined, the position [x, y] of each pixel with the highest value was identified in each group of pixels with values greater than 0.3. The map identifying the eucalyptus plants was then generated, Figure 7.
In terms of overall accuracy in identifying and counting the plants, the developed routine presented a mean value of 0.9353, with a Pacc of 0.9489. The relative error (mean) was 0.11%. Using deep artificial neural networks to assess performance in detecting tobacco plants by means of images from RPA and automated routines, Fan et al. (2018), presented an algorithm with an overall accuracy of 0.9370 and a relative error of 3.94%. The routine showed a high capacity for detecting eucalyptus plants, with an Sb (mean) of 0.9856. However, effectiveness in identifying 'non-eucalyptus' plants, i.e. the (mean) Sp value, was compromised, due to the amount of FP presented. The high FP values are due to the occurrence of weeds or even to shadows of the eucalyptus plants that showed spectral behavior similar to that of the above-mentioned crop (Figure 8).
It was seen that when only one sample was used the total number of plants was underestimated, unlike in the other combinations, albeit with greater overall accuracy. Abstracting from the metrics Pacc and Er, the best estimate for TP was reached using the second combination (two templates in SC). However, FN values increased (Table 1). TN estimates remained constant for all combinations.   The second phase of the algorithm consisted in generating the triangular mesh, which allowed the spacing between plants and between rows to be estimated. On the other hand, the occupation diagram allowed the area of influence of each plant to be identified and estimated ( Figure 9A). Using the first combination (one sample in SC), with 3316 detected plants (Np), a (mean) occupied area of 9.46 m 2 plant -1 was estimated with a coefficient of variation of 12.88%, i.e. a mean occupation 5.11% greater than planned. Although the planting project was planned at 3 m x 3 m (theoretical final stand of 3454 plants), the field survey, followed by the analysis, showed underutilization of the land, with a mean spacing between plants and between rows of 2.79 m and 3.6 m, respectively.
The information generated by the automated routine is useful for monitoring the plantation in a specific area of interest. The variations that occur in the spacing between plants (CV of 6.73%), and especially between rows (CV of 13.73%), can be seen from the RM. The estimated number of plants for replanting was 23, i.e. a survival rate of 99.31%.
The information contained in the histogram ( Figure 9D) also made it possible to define the intervals of the respective spacings, as well as the cut-off threshold (based on the mean and standard deviation) of the occurrence of planting failures. Once this threshold was defined, the edges between the rows and even the failures, were drawn by controlling the angle of each edge on the triangular mesh. It is important to point out that the above angle is a function of the planting arrangement.
Several advances and innovations stand out in this study: (i) an automated routine was presented, scientifically validated, which provides various parameters related to the quality control of forest plantations; (ii) the routine requires minimal user intervention, simply selecting the sample(s) from the area of interest; (iii) the routine was developed on open-source platforms and in Python (a free programming language); (iv) the routine shows the excellent performance of an orthomosaic from RGB channels, i.e. a low-cost sensor.