Deep learning-based prediction of plant height and crown area of vegetable crops using LiDAR point cloud

Study area

In 2017 (January-May), a multi-crop drip-irrigated field experiment was established at the University of Agriculture Sciences (UAS), Bengaluru, India (Fig. 1). This experiment forms part of a larger project aiming to develop field-sensitive techniques using space technology and advanced pattern recognition for detecting crop type, growth condition, productivity, and ecosystem services. (Geographic coordinates:12°58′20.79”N, 77°34′50.31”E). Three different vegetable crops — tomato, eggplant, and cabbage — were cultivated in the field plots according to a factorial design. Thanks to its importance as a computer and electronics industry city, Bengaluru is one of the rapidly expanding cities in Asia and has been experiencing a host of transformations in land use and land cover conversions, and spatial distribution of grains, fruits, vegetables, and herbaceous crops. Geographically, Bengaluru lies off the Western Ghats of India, and the climate is mild, with a mean temperature of 29.2°C and annual precipitation of 873 mm.

Experimental design and data acquisition

Based on the important regional food production and consumption patterns, three different vegetables – tomato, eggplant, and cabbage- were selected for growing on the experimental plots. Aligned with the onset of the Indian monsoon, the crop growing experiments were conducted during the Kharif (June to October) crop growing season. The cultivars selected for the crops correspond to those typically used by farmers across south India. All the experimental datasets on the plants cultivated were acquired complying to the mandatory institutional guidelines applicable.

The layout of the experimental set-up (Fig. 1) consisted of 12 plots, each of size 12 m × 18 m. There were four replications for each crop plot. To assess the nutrient responses and possible discrimination of crops with specific reference to nitrogen level, each of these plots were further divided into subplots of size 6 m × 12 m. Within each plot, three different levels of mineral nitrogen (N) fertilizer were supplied to the subplots randomly. Corresponding to the standard dose of the region and indicated as ‘medium N’ level in the experiment, urea at the rate of 46, 50, and 60 kg N ha^-1 was applied to tomato, eggplant, and cabbage, respectively. The other two levels of nitrogen indicated as ‘low N’, and ‘high N’ correspond to the nitrogen application at the rate of 50% less and 50% more compared to the medium N level. Apart from the nitrogen, phosphorus (P) and potassium (K) were applied at sowing uniformly at the rate of 17.5 19.9, 16.6, and 41.5 kg K ha^-1 for tomato, eggplant, and cabbage, respectively. The nitrogen was applied in two instalments to optimize the uptake and minimization of nutrients leaching due to heavy rainfall. In some of the plots, a moderate level of surface N transport occurred due to the heavy rainfall events during the crop growing experiments. As the primary aim of the work was predictive modelling of crop structural attributes, we have not considered the differential N rates in the experimental plots.

Reference measurements of crop structural parameters: Two different types of datasets were acquired over the experimental farms: LiDAR point cloud and reference crop structural parameters. Concurrent with the LiDAR point cloud acquisitions, reference measurements of plant height and crown area were obtained by distributing sampling locations across the plots. The plant height was measured using a ruler to the nearest centimetre. For each reference height measurement, the heights of 30 plants height were averaged, ensuring distribution in each subplot. Linear distance measurements of of plant crowns in the N-S and E-W directions were taken at several locations in the plots. The plant-level reference crown area was computed from the digital reconstruction of plants from the digital photographs as per²⁴. The reference measurements were used as the training samples in the model development and to validate the predicted values.

LiDAR point cloud acquisition: LiDAR point cloud was acquired for five different sampling dates using a 3D terrestrial laser scanner. The crops were sown on 27 March 2017, and subsequent sampling dates were selected at intervals to capture key growth stages and changes in crop structure. Specifically, the sampling dates were scheduled for 18 April, 4 May, 22 May, 14 June of 2017 to encompass various stages of plant development, from early growth to maturity, allowing for comprehensive data collection and analysis across different growth phases. On the first sampling date, a Riegl TLS (Model: VZ-400; Make: Riegl Laser Measurement Systems GmbH, Austria) was used. For the remaining three sampling dates, the TLS used was that of FARO (FARO 350^S, FARO Technologies Inc., USA). Though the instruments were of different make, the scan density and the associated scan parameters were similar, and the resulting point cloud was seamless. The TLS scanner acquires laser returns in the electromagnetic spectrum’s specified near-infrared wavelength (at 1550nm). The range of sensing is 350m. To cover the entire agricultural field, the TLS was placed at 16 different scan positions spread across different view directions in the experimental plots. Between each pair of successive scans, about 20–25% overlap was maintained to enable co-registration of the multiple scans of the point cloud. The tripod’s height was adjusted to match the growth stage of the crops. To ensure acquiring a high-density point cloud, the scan duration was set to 10 min with a scan-point spacing of 6.1mm at 10m. To maintain point-level precision of geo-positioning of the point cloud, DGPS (Differential-GPS) measurements were acquired at each scan position, and the specified reference checkerboard targets were distributed in the plots. The same procedure was repeated for all the subsequent dates of point cloud acquisition. Figure 2 shows field photographs of a part of the experimental plots during datasets acquisition, and Fig. 3 depicts the acquisition of point clouds obtained from multiple angles.

Methodology

An overall view of the methodological process flow, broadly indicating the critical tasks and algorithms adapted for the realization of the goal of the work, is shown in Fig. 4.

Point cloud processing and crop canopy modelling

Due to the presence of plant support infrastructure and the possibility of multiple reflections from the nearby field boundary trees, some points could have out-of-the-context noise and outlier elevation values. We corrected the scans of point clouds by removing the noise and outlier removal by the neighbourhood thresholding method. Considering the DGPS-measured reference geo-coordinates of the sensor placement and checkerboard targets spread across the plots, the point clouds were co-registered using the multi-scan registration method suggested by²⁵. To help maintain point-level matching across the scans and tide over the lack of sturdy corners in the point cloud of agriculture landscapes, tie points were established based on the guided sphere fit algorithm²⁶. Each data file was geo-referenced to the World Geodetic System (WGS)-84 with UTM projection.

Generation of crop height model (CHM): As the structural parameters considered in this work pertain to the height and canopy of the plant, we classified the co-registered point cloud into the ground and above-ground point clouds. For this, we used the progressive triangular irregular network (TIN) densification (PTD) method²⁷. To reduce the computational challenges, we divided the point cloud into several tiles. The reference surface was constructed based on a triangular irregular network (TIN) built using the ground points identified in each tile of the point cloud. The unclassified points were added to the TIN triangles based on the criteria of minimum distance to the TIN facet and the angle between the TIN facet and the line joining the point with the closest vertex. Continuing the triangulation in an iterative manner, the point cloud was classified into two categories: ground or above-ground. To compare the elevation of the point cloud across different growth phases, the elevation of the point cloud was normalized to the field furnace level. The plant-level canopy height model (CHM) was generated by rasterization of the normalized above-ground point cloud.

Crown area estimation: Plant level crown area was estimated by identifying the crown top of each plant and the subsequent delineation of the crown area. The plant crown tops were identified using the variable window filter (VWF) method proposed by²⁸. For each pixel in the CHM, which represents plant height, the plant crown top was identified based on the local maxima defined by the radius of the window defined adaptively. Based on the centroid of intensity distribution indicated by the outcome of the VWF, the plant level crown area was modelled based on the semi-supervised segmentation of the raster CHM using the Watershed segmentation algorithm²⁹. The possibility of over-segmentation was minimized by removing the segments which do not have associated crown top values.

Following the retrieval of plant height using the developed canopy height model (CHM) and crown area estimation utilizing variable window filter (VWF) for the three sampling dates, these parameters were inputted into the TemporalCropNet architecture to predict crop height for the fourth sampling date. In the case of tomato, sampling of the first three dates was considered as tomato crop was harvested before 4 June 2017.

Deep learning-based prediction of plant height and crown area

Deep learning (DL) techniques for image processing and pattern recognition in computer vision have reached maturity for operational applications, and a host of commercial products are developed using the standard optical RGB imagery^30,31,32. Several DL techniques have been experimentally applied on multispectral and hyperspectral imagery for several tasks, such as segmentation, classification, unmixing etc. Point clouds acquired from various platforms have been processed using DL techniques for object detection, semantic segmentation, object-based landscape modelling and visualization etc. Several DL-based standard point cloud processing models are available^12,33,34,35. Among the wide variety of model architectures and networks reported in the literature, convolutional neural networks (CNNs) and Recurrent Neural Networks (RNNs) are widely used for processing and analyzing remote sensing data. While CNNs are best suitable for instance-based processing tasks such as classification, RNNs are superior for processing time-series data.

The primary DL architecture adopted in this work is a variant of the RNNs. Unlike the feed-forward neural networks, RNNs use their internal state memory to process the data sequences. Recent studies confirm that DL networks such as Long Short-Term Memory (LSTM) and Gated Recurrent Unit (GRU) overcome the drawbacks of RNNs, such as vanishing gradient or exploding gradient problems. Based on the observations from the preliminary assessment of the performance of the LSTM and GRU for predicting the crop structural parameters independently, we used a stacked representation of LSTM and GRU. Combining multiple LSTM layers leads to greater model complexity and high-level time dependency. Stacked LSTM provides output for each time stamp and not the single output for all time steps. Similar is the case of GRU. Thus, a hybrid hierarchical model generated by stacking the LSTM and GRU helps improve the projection of information in latent dimension space, giving better prediction results. The stacked model was then fed to an individual fully connected layer (FCN) and then combined and given to the final FCN, where the feature concatenation occurs and the parameters are predicted. For ease of reference, we name the deep learning model developed for this specific purpose as ‘TemporalCropNet’.

The top-level architecture of the TemporalCropNet is shown in Fig. 5. The crop parameters– height and crown area pertaining to different data acquisition dates are given as input to the TemporalCropNet. In summary, the hybrid hierarchical model of stacked LSTM and GRU layers, followed by fully connected layers, forms the core of TemporalCropNet, demonstrating significant improvement in predicting crop structural parameters. For ready reference, a brief description of the recurrent neural network (RNN), long short-term memory (LSTM), and gated recurrent unit (GRU) is given in Fig. 5.

Recurrent neural networks (RNNs): Unlike the feed-forward neural network such as CNN which has only a finite receptive field, a recurrent neural network (RNN) can process sequential data using the internal state memory. This dynamic nature of RNN finds application in audio, speech analysis and several other temporal-based applications. The recurrent state is achieved by connecting the time steps’ output as input to the network and using the same weights to backpropagate through time. Therefore, RNN is a neural network that uses the shared weights for each time step. At the same time, RNN suffers from short-term memory when the sequence length becomes larger and suffers from vanishing gradient problems during backpropagation. Neural network weights are updated using gradients while backpropagating. The issue of vanishing gradient produced leads to minimal gradients and contributes nothing to the learning process, thus leading to short-term memory.

Long Short-Term Memory (LSTM): The cell state and gates are core structural elements of the sequential network, LSTM. A high-level representation of LSTM is similar to that of RNN but differs in terms of the internal representation of the memory cell. The input to the sequential network is a sequence of vectors (x₁, x₂,… x_n), where x_t is the input at time t. The three basic units of the LSTM are forget gate (f_t), input (y_t), and output gate (p_t). The forget gate determines to keep or discard the information from the previous time step. From the input gate, LSTM learns new information that enters the cell. The output gate passes the updated information in the cell from the present time step to the next. Similar to RNN, LSTM also has hidden state d_t, d_t-1, representing the current and next step’s hidden unit. Also, the cell state is represented by c_t and c_t-1 for the current and the next time-step, respectively. The cell state and hidden states are the two states of LSTM. All the gates undergo the sigmoid activation function, which gives the output value for each state and gate ranging between 0 and 1. In addition to that, LSTM also has a temporary cell state C’_t, where the hyperbolic tangent activation function is applied that outputs values between -1 and + 1. In The final cell state C_t, the forget gate f_t and input gate y_t determine how much information should be kept in the current time step. The new hidden state h_t for the next time step is determined by C_t and p_t, which determines the amount of information the memory cell has at the output step. For all cell states and gates, in addition to the weight parameter, bias value is also added to each equation from (1)—(5). As represented in Fig. 4, the hidden state dₜ and the cell state C_t for an LSTM cell can be calculated as follows.

Here, y, f, p is called the input, forget, and output gates, respectively; W is the recurrent connection at the previously hidden layer and the current hidden layer, T is the weight matrix connecting the inputs to the current hidden layer. The GRU is an RNN architecture and is similar to LSTM units. The GRU comprises of the reset gate and the update gate instead of the input, output and forget gate of the LSTM. The reset gate determines how to combine the new input with the previous memory and analyses how much information is to be embedded with the current information. The update gate defines how much of the previous memory to keep around.

Gated Recurrent Unit (GRU): LSTM has a number of parameters which result in a number of operations to be performed. Thus, the variant of LSTM, GRU, is a simplified, compact representation of LSTM in terms of the parameters and operations performed. Unlike LSTM, GRU has only two gates: reset (s_t) and update (u_t), instead of the input, output, and forget gate of LSTM and a hidden state h_t. The reset gate determines how to combine the new input with the previous memory, and the update gate defines how much of the previous memory to keep around. Also, the two gated combine the input x_t and d_t-1 (information from previous state t-1). Since GRU has lesser tensor operations, it is faster to train the network than LSTM. The computation of z_t, d_t is given in Eqs. (6) – (9) expressed below.

The Stacking Mechanism:

The stacked LSTM and GRU layers can be expressed as:

where \({H}_{t}^{L}\) and \({H}_{t}^{G}\) in Eqs. 10 and 11 denote the hidden states of the LST and GRU layers at time t respectively, and L and G representes the number of stacked layers.

The outputs from the stacked LSTM and GRU layers are first fed into individual fully connected layers (FCNs). These individual FCN outputs are then concatenated and passed through a final FCN layer. Each stacked model’s output H_t is processed through its respective FCN:

where ϕ is the activation function, \({W}_{f}^{L}\) and \({W}_{f}^{G}\) are the weight matrices, and \({b}_{f}^{L}\)​ and \({b}_{f}^{G}\) are the bias vectors of the individual FCNs for LSTM and GRU, respectively in Eqs. 12 and 13.

The outputs \({y}^{L}\) and \({y}^{G}\) are concatenated and passed through the final FCN for prediction:

where \(\widehat{y}\) is the predicted output, \({W}_{final}\) is the weight vector of the final FCN in Eqs. 14 and 15.

The training, test and validation datasets were taken as 60, 20, 20% of the total dataset. The predictions from the DL model were validated by comparing with the actual plant-level measurements. So as to ensure the robustness of the model to fit into unseen data, we have devised a fivefold cross validation strategy. The quality of the matching is represented by comparing the symmetric mean absolute percentage error (SMAPE) expressed in Eq. (16) ³⁶.

Symmetric mean absolute percentage error,

Visualization of possible underestimation or overestimation relative to the one-to-one comparison of the measured and predicted values of the crop structural parameters is expressed using another graphical error metric, Logarithmic deviation³⁶, given by

where y_i is the estimated parameter, x_i is the measured parameter for the measurement pair.

Ethical Approval

This study has not used data or samples pertaining to any humans or animals. No ethical committee approval is required for this study as such. This study has used data acquired on cultivated plants. We would like to declare that the study has complied with all the relevant guidelines applicable.