Assessment of canopy vigor information from kiwifruit plants based on a digital surface model from unmanned aerial vehicle imagery

: Information about canopy vigor and growth are critical to assess the potential impacts of biotic or abiotic stresses on plant development. By implementing a Digital Surface Model (DSM) to imagery obtained using Unmanned Aerial Vehicles (UAV) , it is possible to filter canopy information effectively based on height, which provides an efficient method to discriminate canopy from soil and lower vegetation such as weeds or cover crops. This paper describes a method based on the DSM to assess canopy growth (CG) as well as missing plants from a kiwifruit orchard on a plant-by-plant scale. The DSM was initially extracted from the overlapping RGB aerial imagery acquired over the kiwifruit orchard using the Structure from Motion (SfM) algorithm. An adaptive threshold algorithm was implemented using the height difference between soil/lower plants and kiwifruit canopies to identify plants and extract canopy information on a non-regular surface. Furthermore, a customized algorithm was developed to discriminate single kiwifruit plants automatically, which allowed the estimation of individual canopy cover fractions ( f c ). By applying differential f c thresholding, four categories of the CG were determined automatically: (i) missing plants; (ii) low vigor; (iii) moderate vigor; and (iv) vigorous. Results were validated by a detailed visual inspection on the ground, which rendered an overall accuracy of 89.5% for the method proposed to assess CG at the plant-by-plant level. Specifically, the accuracies for CG category (i)- (iv) were 94.1%, 85.1%, 86.7%, and 88.0%, respectively. The proposed method showed also to be appropriate to filter out weeds and other smaller non-plant materials which are extremely difficult to be distinguished by common colour thresholding or edge identification methods. Assessment of canopy vigor ınformation from kiwifruit plants based on a digital surface model from unmanned vehicle ımagery. J & Biol Eng, 2019; 12(1): 165–171.


Introduction 
Kiwifruit is an important cash value crop and their cultivation has become a major agricultural industry that promotes the economic development of the burgeoning kiwifruit producing regions in the northwest of China, the northeast of New Zealand, as well as Italy. However, traditional manual methods to manage a kiwifruit orchard are time-demanding and usually based on the monitoring of sentinel plants with extrapolation of management decisions without considerations of the natural spatial and temporal variability of orchards. Plant biotic and abiotic stresses can result in a partial or complete canopy and plant losses within an orchard. Furthermore, these effects can lead to significant economic losses depending on the severity of the specific detrimental effects which have a direct relationship with the production and quality of kiwifruits in the short and long term. Therefore, it is urgent to develop fast, cost-effective and near real-time methods to obtain more accurate information from kiwifruit orchards that reflect accurately the spatial and temporal variability at the plant level for efficient decision making and crop management to maximize fruit quality and yield, especially in a changing climate.
Remote sensing technologies using different platforms, such as satellite, airborne and more recently Unmanned Aerial Systems (UAS), have been increasingly applied in the agricultural industry worldwide. Within the optical/visible domain (400-700 nm; RGB imagery), it offers a potential to extract spatial information from the canopy growth (CG) of kiwifruit plants in a non-destructive manner, with the objective of the development of decision support systems [1,2] . The UAS have been transformed nowadays in affordable and flexible solutions that can provide images at high spatial, temporal, and spectral resolution. UAS allows planning flights at a date that is very close to targeted phenological stages, allowing the detection of plant biotic, abiotic stresses or for pre-harvest characterization [3][4][5] . Furthermore, UAS has been one of the technologies with the highest rate of development and implementation for precision agriculture applications in recent years [6,7] . In addition, most of the visible/near-infrared/thermal infrared indices that can be calculated from UAS remote sensing do not require atmospheric corrections, compared to satellite remote sensing [8] due to low altitude surveys required (50-100 MAGL).
For kiwifruit orchards, higher spatial/temporal resolution based UAS remote sensing can provide a potential and efficient way to identify and measure CG. One of the main issues in row-planted crops is the individual plant characterization, which requires the isolation of pixels belonging to canopies from those belonging to soil or any other vegetation material such as weeds or cover crops. One common approach is to differentiate spectral properties between rows and inter-rows, which require multispectral cameras, it was first explored by applying a threshold on vegetation indices computed from multispectral images [9][10][11] and, recently, machine learning algorithms [12,13] . The analysis based on spectral bands can result in some problems associated to the same plant material presenting different spectra (different growth level), separation of plant material and shade, different plant material with similar spectra (such as the case of weeds or cover crops in the inter-row) and mixed information per pixel such as soil/plant associated to the specific spatial resolution from the imagery [14] . Another approach is by recognizing the edges of the canopy using different edging algorithms applied to multispectral or infrared thermal imagery [15] . However, for multispectral, this method will encounter similar problems described before and for infrared thermal it could include in the analysis heat effect from soils and other objects (training wires and poles) specially in the edges of canopies. An alternative approach is a system based on plant shape and architecture, such as by the implementation of the DSM, which would discriminate more efficiently the regions of interest or plants through its canopy height, which are significantly different to ground and other objects height, such as grass, weeds, and cover crops.
However, sub-meter evapotranspiration estimation systems from remote sensing data also require specialized personnel for aircraft operation, data acquisition and analysis.
Some of these applications require higher instrumentation costs compared to high definition RGB cameras. Other studies, based on both the spatial and the spectral properties of crop rows, applied a Fourier transform to the red band of RGB images as a differentiation criteria [21] . However, these remote sensing methods for CG assessment are also based on multispectral remote sensing, which also includes the drawbacks mentioned before.
The kiwifruit plant is a vine-like plant that needs to be trained and planted in rows with fruits positioned at breast level for easy picking at harvest. A Digital Surface Model (DSM) is based on identifiable ground objects including buildings, forest trees and crop plants, among others, which can provide elevation information of various ground objects to show surface undulation conditions and crop growth conditions [22] . The DSM approach has been successfully applied to major crops [23] , for example, Adam et al. [24] , used the Structure from Motion (SfM) algorithm to model grapevine canopy structure in a vineyard site in the Texas Hill Country, and analyzed the data using a stepwise regression model to attempt to predict LAI, results showed a moderate R 2 value of 0.6. Grigorijs et al. [25] used point clouds derived from SfM matching techniques obtained from UAS to detect individual trees, measured tree heights, and provided RGB estimates in Australian tropical savannas. Marie and Fré dé ric [11] developed an algorithm for vineyard structural characteristic estimation based on dense point clouds derived from the RGB colour model images acquired with a UAV. Finally, Su et al. [22] used DSM acquired from UAV RGB images to assess affected and missing grapevine canopies affected by frost.
This work describes the development of a DSM for a kiwifruit orchard derived from the RGB colour model based images acquired using a UAV through three-dimensional (3D) reconstruction with SfM algorithm techniques.
Furthermore, photogrammetry algorithms were applied to assess the CG at the plant-by-plant level to be finally classified into four classes of CG: (i) missing plants; (ii) low vigor; (iii) moderate vigor; and (iv) vigorous.

Experimental site
The experimental field is located at the kiwifruit orchard belonging to Zhouzhi prefecture in the north Qinling Mountain Area of Shaanxi Province (34.14°N, 108.09°E) (Figure 1). The climate of the region is classified as a temperate arid with a mean annual precipitation of about 180 mm. The frost-free season corresponds to 200 d per year from April to October. The soil is mainly clay in texture with low water infiltration velocity. After many years of development, this area has become a renowned region for kiwifruit in China. The experimental kiwifruit orchard consisted in 5-year-old kiwifruit plants, cultivar Green kiwifruit, planted in a S-N row direction with spacing between plants and between rows of 1.5 m×4.2 m respectively with a total of 1587 plants per hectare.

UAV Imagery Acquisition
In this study, a DJI Phantom 3 Vision plus quadcopter (SZ DJI Technology Co., Ltd., Shenzhen, China) was used as the survey platform. This aircraft is an affordable small sized four-axis quadrotor aerial vehicle with a flight control system, a three-axis stabilizing gimbal and WiFi communication capabilities. The quadcopter specifications are: 1280 g of weight, 4480 mA· h of battery capacity, a maximum speed of 15 m/s, and maximum flight time of around 25 min. To reduce the complexity of the UAS operation, the Pix4Dcapture software (Pix4D, Lausanne, Switzerland) was used to control the flight and capture the RGB images. The Pix4Dcapture allows using the UAS as a mapping and measuring tool by defining autonomous mapping flights through pre-defined waypoints for data acquisition.
A parrot sequoia camera (Parrot, France) was mounted on the UAV to be used as the remote sensor ( Figure 2). The camera specifications are: 16 MP RGB sensor, four narrowband and synchronized 1.2 MP monochrome sensors. The camera can be used to take pictures of agricultural fields in several spectral bands which measure the state of the vegetation: Green (550 nm wavelength, 40 nm bandwidth), Red (660 nm wavelength, 40 nm bandwidth), Red Edge (735 nm wavelength, 10 nm bandwidth) and Near Infrared (790 nm wavelength, 40 nm bandwidth). It is an automatically calibrated camera because of the sunshine module and has an integrated GPS/GNSS to locate the camera when photos are being taken. In addition, users can connect a smart mobile device wirelessly to the UAV and the remote camera via WiFi to control the drone and to obtain real-time transmission of images and videos. Remote images were acquired through UAV fights that were performed on May 16th, 2017 during the kiwifruit growing period, the images were recorded on a flash card, allowing 60% to 80% overlapping along and between the tracks, according to the nominal 25.2 km/h UAV flight speed at 40 m flight altitude. For this study, only the RGB bands were considered, which can be accessed from any ubiquitous digital camera.

Geospatial analysis
The SfM technique facilitates the establishment of DSM and consists in a digital reconstruction technique, which is overlapped are extracted from a series of images captured from different angles of view to achieve a 3D reconstruction according to the image feature matching algorithm [17] . With the camera parameters and 3D data, the SfM technique can be used to perform motion estimation according to the geometric relationships between 2D images of multiple view angles. Then the SfM optimizes the motion calculation and determines 3D points with bundle adjustment to finally obtain a dense 3D point cloud through point cloud extension [13] . Finally, the DSM was built on the basis of the estimated camera position and the images themselves. In this research, the Pix4Dmapper Professional software was used for this task, Pix4Dmapper is a 3D reconstruction software based on the principle of SfM and the mosaicking process was fully automated involved the following phases: (1) Initial Processing, (2) Point Cloud and Mesh, (3) DSM, Orthomosaic and Indexing. After feature extraction and matching of UAS aerial images, a 3D point cloud was generated and the DSMs were obtained in Tag Image File Format (TIFF).
To extract the planting pattern and characteristics of kiwifruit orchard, the GIS software ArcGIS (Esri, Inc., Redlands, UAS) was used for georeferencing and crop the DSM for the selected study site. Figure 3 presents an overview of the algorithms and flow chart used to extract CG information from the kiwifruit orchard. The whole analysis process start with the image binarization from the overlapping RGB image acquired with the UAS (Figure 3, Step A). This is achieved by obtaining the DSM with the Pix4Dmapper software and then applying an adaptive threshold algorithm to extract plant positions and to filter weed and soil background. Then, RGB image is rotated and cropped to adjust planting row to a horizontal position to establish anchor points and masks as marks for further assessment (Figure 3, Step B). Finally, the canopy vigor of the whole kiwi orchard is assessed (Figure 3, Step C), firstly, the binary image is rotated, cropped and generating reference points and masks with the same position as anchor points and masks in the RGB image to calculate canopy cover fraction (f c ). Then, every single mask is judged according to the specific threshold value.  Figure 3 shows the DSM generated from the RGB images obtained with the UAS. The highest point for discrimination within the research region was 2 m above the ground. There was a distinctive difference in height between the kiwifruit trees (white stripped area) compared to objects in the inter-row.

Adaptive threshold algorithm to obtain plant position (Step A)
The thresholding method was then exported as a binary image. At this stage, the DSM allows only qualitative analysis of canopies height in relation to ground, therefore, the area corresponding to missing plant cannot be assessed directly via the DSM of study region because of the complicated surface structures of the undulating ground in the kiwifruit orchard. Figure 4 also shows the kiwifruit trees and their distribution in the non-regular or sloped fields, such as those encountered in the trial site. This presented a problem since a common discriminating threshold will over or underestimate the real height of the whole canopy (i.e. Line1, Line2, and Line3 corresponding to three different thresholds).
As a result, a local threshold method was implemented, the adaptive threshold algorithm, which was denoted as T(x, y) was implemented by setting a subset of pixels in a square area Ω x,y of size I w to calculate the local threshold by subtracting a constant C obtained from the average elevation of the neighborhood pixels from the pixel point (x, y). The local thresholds are then used to filter the kiwifruit tree from the background, the constant C is set to extract canopies which are lower than the average height in the sliding window. The size of I w and the value of C was confirmed according to the actual situation of the experiment. Using src(x, y) to represent the pixel value from each point in the original image, dst(x, y) represent the value of the pixel in the resulting binary image, with the binarization process represented mathematically by Equation (1).  Figure 4 (a) DSM extracted from the of kiwifruit orchard; (b) Diagram representing terrain slope changes within the orchard and the effect on a fixed height thresholding Figure 5 Resulting binary image obtained by using the adaptive threshold algorithm for height selection for a kiwifruit orchard studied corresponding to the polygon in blue As shown in Figure 5, the binary image obtained after the application of an adaptive thresholding algorithm for height selection contains all the boundaries from kiwifruit plants. This method can also remove the gradient background effectively and extract kiwifruit plants in non-uniform or slope ground conditions.

Establishment of reference points and assessment masks to separate canopy information from the inter-row material (Step B)
The homogeneous plantation pattern presented in this study can be used to separate the information from the processed image form the kiwifruit orchard in a plant-by-plant scale. Therefore, individual kiwifruit canopies can be detected by creating a set of reference points by linear interpolation, each point corresponding to a detection area of interest corresponding to a whole kiwifruit canopy.
Based on this principle, the initial row positions and end points of 12 sampled rows (red crosses; Figure 6a) were assigned in the RGB image and were used to calculate the row orientation, the assigned rows were selected at intervals of two rows(with a total of 36 rows). Using θ to represents the orientation between the kiwifruit planting row and the image X coordinate axis, x 1 and x 2 are the horizontal components and y 1 and y 2 are the vertical component of the verges of rows, then the orientation from each selected row can be computed by Equation (2) [26] .
The orientation from the west-east of linear features (i.e., rows in this case) was given by θ. The average angle of the 12 rows obtained is shown in Table 1. Once the averaged angle of the rows was calculated, the image was rotated by the averaged θ (Figure 6b).  The rotated image requires the creation of reference points of interest to identify individual plants. The reference points refer to stem position from plants or the centroid of the designated planting area. In order to achieve this, it was necessary to combine the known data of planting density for rows and inter-row distances of 1.5 m×4.2 m respectively. With this information, it is possible to automatically obtain the number of rows and the kiwifruit trees cultivated in each row, which corresponded to 36 and 100, respectively. As a result, the pixel position corresponding to the centroids can be computed by the pixel interval of 36 rows and 100 points per row through linear interpolation to generate points and masks automatically. The rotated image (Figure 6b) was then used to calculate automatically the position of areas corresponding to individual plants that allows the calculation and storage of any information related to pixel count or channel (RGB) data using programming loops (Figures 7a and 7b).

Binary image generation to assess of kiwifruit CG (Step C)
Reference points and masks with the same position as RGB image were then set in the binary image. From each sub-image corresponding to individual kiwifruit plants a f c can be calculated automatically, extracted and stored as pixels of canopy over total pixels of the sub-image. Figure 7a shows the resulting binary image of the kiwifruit orchard study area, which can be used as a mask to be overlapped with the processed image to obtain specific plant-by-plant areas or RGB information. Figure 8b shows the  The canopy vigor classification was quantified as i) missing plants using the criteria based on inexistent o minimal canopy pixels encountered in the specific analyzed region (f c ≤0.08) , ii) low vigor (0.08< f c ≤0.3), iii) moderate vigor (0.3< f c ≤0.6) and iv)vigorous (f c > 0.6)). The missing plant criteria was based on an extremely low level of vegetation, the low vigor trees related to plants that are partly affected by biotic/abiotic factors, the moderate vigor level was related to well-growing condition and less affected canopies, while the vigorous level corresponded to high vegetative growth that can be related to overgrowth or a dense biomass, which can be a result from over-irrigation or fertilization.
Using these criteria and thresholds, the sub-images were classified by colour automatically (Figure 10) as: (i) missing plants (white filling corresponding to 1517 plants or 42.14% of total plants); (ii) low vigor (yellow filling corresponding to 713 plants or 19.81% of total plants); (iii) moderate vigor (bright blue filling corresponding to 992 plants or 27.56% of total plants); and (iv) vigorous trees (green filling corresponding to 378 plants or 10.50% of total plants). Most of the missing kiwifruit plants in white filled boxes were surrounded by affected kiwifruit trees (low vigor), which makes sense if they were affected by biotic/abiotic stresses. On the contrary, vigorous canopies in green filled boxes are surrounded by the vigorous canopy and moderate vigor canopy gathered in the middle part of studied kiwifruit orchard. The accuracy test for the identification and classification method proposed was evaluated by comparing results using the algorithms proposed against the real growing condition obtained by visual inspection of expertise from the study area. Results showed high accuracy in the estimation of canopy cover by the method proposed with 88.0% correct identification of vigorous, 86.7% of moderate vigor, 85.1% of low vigor, 94.1% of missing plants and an overall accuracy of the method of 89.5% (Table 2). These results can be also applied to improve the accuracy of the DSM method implemented to separate canopies from the background, since over and underestimations of canopy cover were minimized. It is important to note that the method proposed produces an automatic classification of vigor of kiwifruit plants in a map form, but also the classification is numerical. Therefore, growers can have the specific row and plant with the respective classification on a list, which could help significantly to perform targeted managements for plant replacement, canopy management, fertilization and irrigation at a plant-by-plant scale. This method can also be used to obtain other canopy architecture parameters, such as leaf area index (LAI) and canopy porosity (ϕ) using the same algorithms for cover photography and computer application (VitiCanopy) developed by De et al. [27] , Arachchige et al. [28] also used 'Canny edge detection' and 'Otsu's' methods to derive canopy porosity to quantifying the severity of phytophthora root rot disease using RGB digital photographs.

Conclusions
A method for automated assessment of canopy vigor for a kiwifruit orchard using DSM from Unmanned Aerial Imagery was presented.
The proposed method effectively reduced the complexity of background information in kiwifruit orchard by filtering topography using an adaptive threshold algorithm and accurately extracts kiwifruit rows. Conversely to methods based on classification of multispectral or RGB images, the method proposed in this paper is not sensitive to the presence of weed between the rows.
Since the method proposed produces an automatic classification of vigor of kiwifruit plants in a map form, but also the classification is numerical, growers can have the specific row and plant with the respective classification on a list to target management practices. Other canopy architecture parameters, such as leaf area index (LAI) and canopy porosity (ϕ) can also be obtained using the same algorithms for cover photography and computer application.
However, the accuracy of the results can be affected by the quality of Unmanned Aerial Imagery, since the quality of DSM rest with the flight condition such as UAV speed, camera resolution, field of view and wind. Besides, the threshold of f c need to adjust when implementing the method in other orchard because of the difference between rows and inter-row, as well as the canopy size.