Pheno-Robot: An Auto-Digital Modelling System for In-Situ Phenotyping in the Field

In this study, we utilise a novel neural network, the 3D Object Detection Network (3D-ODN), as introduced in our previous study [19]. Specifically designed for processing point cloud maps in Bird’s-Eye View (BEV), the 3D-ODN consists of three integral network branches: point cloud subdivision, feature coding, and a detection head. First, the entire point cloud map is uniformly segmentated. The points within each segment are then projected onto the BEV perspective, where local semantic features are extracted. These features are then fed into the main branch, forming a single-stage network architecture dedicated to processing the point cloud features from the BEV angle. To improve the learning of multi-scale feature embeddings, the feature pyramid is used to facilitate the fusion of features of multi-scales. Finally, the recognition results are projected into 3D space based on the predicted height values. The identified instances are denoted as $\mathbf{T}=\{t_{1},t_{2},\cdots\rvert t_{i}\in\mathbb{R}^{2}\}$, representing the instances detected by the 3D-ODN (see Fig.3 (a)).

For each instance, eight nodes are designated to represent the four corners in two different directions of motion, as shown in Fig.3 (b). This orientation consideration aims to mitigate turning movements in narrow passages, which could otherwise lead to a higher failure rate of movement due to trapping into potential obstacles. Consequently, each instance $o$ is associated with eight nodes, denoted as $t_{i}=\{N^{i}_{1},N^{i}_{2},\cdots,N^{i}_{8}\rvert N^{i}_{j}\in\mathbb{R}^{2}\}$ (see Fig.3 (b)). To extract the line structure of farms, a line detection algorithm is used to group instances in the same row, denoted $\mathbf{g}$. Each row has two common access nodes at both ends, denoted $N_{l}$ and $N_{r}$ ($N_{l}$ and $N_{r}\in\mathbb{R}^{2}$). These nodes serve to establish connections between all instances within the group and to form links with other groups (see Fig.3 (a)). A group is thus characterised as $g_{i}=\{t_{1},t_{2},\cdots,N^{i}_{l},N^{i}_{r}\}$. The map is then represented as a set of groups denoted by $\mathbf{G}=\{g_{1},g_{2},\cdots\}$.

A Greedy-search-based path generation algorithm on the graph map is developed, as illustrated in Alg.1. The following are the relevant definitions:

• •

  $\mathbf{T}$ represents the goal set that requires phenotyping.

• •

  $\mathbf{V}$ represents the subgroups if their instances in $\mathbf{T}$.

• •

  $\mathbf{N_{start}}$ represents the robot position at the start.

• •

  $\mathbf{\Gamma}$ represents the global path that includes a sequence of nodes.

This step aims to compute the detailed path based on the global path, including three steps: path generation, optimisation, and interpolation.

Initial path Generation: Given the global path $\Gamma$, the first step is to generate an initial collision-free path between two nodes. Two scenarios are considered, including path generation for phenotyping data acquisition (where two adjacent nodes belong to the same instance) and paths between nodes of different instances. In the first case, a sequence of collision-free preset sample positions is established and the RRT algorithm is used to connect these positions to form the sampling trajectory between nodes. In the second case, the $\mathbb{A}^{*}$ algorithm is used to determine a collision-free path between nodes.

Trajectory Optimisation: The initial path $\Phi$ consists of a sequence of points $\{Q_{0},Q_{1},\cdots|Q\in\mathbb{R}^{2}\}$, and we conceptualise a trajectory $t\in[0,1]\rightarrow\Phi\subset\mathbb{R}^{2}$ as a continuous function mapping time to robot states. The objective function incorporates three key aspects of the robot’s motion. It penalises velocities to encourage smoothness, proximity to the environment to secure trajectories that maintain a certain distance from obstacles, and the distance between the state and the preset viewpoint. These terms are represented as $f_{s}$, $f_{c}$, and $f_{o}$, respectively. The objective function is:

We update the trajectory by the functional gradients $\bar{\nabla}f(\Phi_{i})$ using $\Phi_{i+1}=\Phi_{i}-l_{\mathbf{r}}\cdot\bar{\nabla}f(\Phi_{i})$, following the work [20], where $l_{\mathbf{r}}$ is the learning rate. The functional gradient of the objective in (2)-(4) is given by

where $\nabla c$ is the derivative of obstacles to the control points $Q_{j}$ by the $\partial c/\partial Q_{j}$, the definition of $\kappa$ is given as below.

B-spline Interpolation: The optimised trajectory is then parameterised by a piece-wise B-spline into a uniform curve $\Phi_{B}$. Given the determined degree $m$, and a knot vector $\{k_{0},k_{1},\cdot,k_{M}\}$, where $M=n_{q}+2m$. The parameterised uniform curve $\zeta(t)$ can be formulated as:

where $w_{i}$ is the weight of each control points $Q_{i}$, $R_{i,m}$ is the basis function for $Q_{i}$ at $m$ degree, $u$ is the control value of the curve. Then, the TEB-planner [21] is used to find the proper velocity to follow the trajectory.

The $\mathbf{FindNearestAndFeasible}$ function in Alg.2 evaluates the feasibility between two nodes from two aspects: orientations and traversability.

I. Orientation Checking: This term evaluates whether two nodes have similar directions, as the narrow tunnels between plants do not allow turn-round. The largest allowed orientation difference between two nodes is $1/3\pi$.

II. Traversability Analysis: The point cloud between two nodes is projected onto a 2D grid map. The traversability of each grid is assessed as the multiple risk factors, including:

(a) Collision risk: A risk factor quantified by the possibility of a grid belonging to obstacles, denoted $\theta\in[0,1]$. A terrain analysis [22] is utilised here.

(b) Slope risk: For each point on the terrain $T_{i}$, points on the map are selected by a cube box with sides of length $l_{s}$, which is represented as $\Omega_{i}=\{(p^{j}_{i})_{j=1:N_{i}}\rvert p_{i}^{j}\in\mathbb{R}^{3}\}$. The SVD is used to fit a plane $P_{i}$ from the $\Omega_{i}$ and get the normal vector $n_{z}\in\mathbb{R}^{3}$. The slop angle between terrain $T_{i}$ and vertical direction $n_{z}$ is obtained by $\mathbf{s}=\mathbf{arccos}\frac{||e_{z}\cdot n_{z}||}{||e_{z}||\cdot||n_{z}||}$.

(c) Step risk: This term evaluates the height gap between adjacent grid cells. Negative obstacles can also be detected by checking the lack of measurement points in a cell. The maximum height gap is denoted as $\lambda$.

The above three risks are combined by a weighted sum to obtain the weighted traversability value $\Upsilon$, as:

where the $\mathbf{s}_{crit}$ and $\lambda_{crit}$ are the maximum allowed slope angle and height gap, respectively.

III. Terrain analysis for trajectory optimisation: For robot operation, the odometry and KD-tree [23] are utilised to build the local terrain map (Fig.5 (a)). The detected obstacle points (Fig.5 (b)) are converted to polygons, denoted as $\{\mathcal{O}\subset\mathcal{R}^{2}\}$. Let the $\phi(Q_{i},\mathcal{O})$ describe the minimal Euclidean distance between obstacles and a control point. A minimum separation $\phi_{\mathrm{min}}$ between all obstacles and $Q_{i}$ is found by

The obstacles $\mathcal{O}$ is updated online for dynamic environments.

NeRF represent the 3D scene as a radiance field which describes volume density $\sigma$ and view-dependent color $c$ for every point $x$ and every viewing direction $d$ via a MLP [24]:

The ray tracing-based volume rendering is used to render the parameters of a 3D scene into colours $\hat{C}$ of a ray $r$, expressed as:

Where $T$ is the volume transmittance and $\delta$ is the step size of ray marching. For every pixel, the squared loss on photometric error is used for the optimization of MLP. When applied across the image, this loss $\mathcal{L}_{\mathrm{Rendering}}$ is represented as:

where $C(r)$ is the ground truth colours.

The artefact of ”white floaters” caused by rendering distortion is the common failure mode in NeRF learning. Autonomous data acquisition by robots can always lead to imperfect density and sparse views with fewer overlapping regions (Fig. 6 (a)), thus distorting the rendering of these regions. This is essential because, between these sparse views, NeRF’s training lacks sufficient information to estimate the correct geometric information of the scene, leading to significant and dense floats in the region that is close to the camera. To reduce those artefacts, an optimisation term $\mathcal{L}_{occ}$, which is designed to regulate the learning behaviour of NeRF [25], is used to penalize the dense fields near the camera via “occlusion” regularization, which can be expressed as:

where $\boldsymbol{\sigma_{K}}$ represents the density values of the $K$ points sampled along the ray, ranked in order of proximity to the origin, $\mathbf{m_{K}}$ is the binary mask vector that determines whether a ray sector will be penalised or not.

Given a predefined 3D region of interest, a set of spatial points $P=\{p_{1},p_{2},...,p_{n}\}$ is generated via dense volumetric sampling. For each point $p_{i}\in P$, it’s evaluated through the NeRF model to obtain The density values, $\sigma(p_{i})=\mathrm{NeRF}_{\sigma}(p_{i})$, form the basis for surface extraction. The Marching Cubes algorithm identifies the iso-surface by threshold of the density values:

Where $M$ is the resultant mesh and $\sigma_{threshold}$ is an optimal density value demarcating the object’s boundary.

For vertex $v_{j}$ in mesh $M$, a viewing ray $r_{j}$ is constructed and queried $\mathbf{c}(v_{j})=\mathrm{NeRF}_{\mathbf{c}}(v_{j},r_{j})$ in NeRF. Those values derived from the field are mapped onto mesh $M$, assigning colour to each vertex:

For a 2D texture representation, the vertex-coloured mesh undergoes UV unwrapping. To minimize distortion, Least Squares Conformal Mapping (LSCM) [26] is used, which is to minimize the conformal energy:

where $(u,v)$ are the 2D texture coordinates for each vertex in $M$, $\Omega$ represents the object’s surface, and $dA$ is a differential area element on the mesh’s surface, indicating that the energy is computed by integrating over the surface of the mesh.

The Instant-NGP [27], which makes use of a multi-resolution hashed encoding that can represent learned features of the scenes with tiny MLPs, is utilised. In detail, Instant-NGP operates on the premise that the object to be reconstructed is enclosed within multi-resolution grids. For any point $\mathbf{x}\in\mathbb{R}^{3}$ in various resolution grids, it obtains the hash encoding $h^{i}(\mathbf{x})\in\mathbb{R}^{d}$, where $d$ is the features’ dimension, $i$ is the level of tri-linear interpolation. The hash encodings of all levels are concatenated to form the multi-resolution feature $h(\mathbf{x})=\{h^{i}(\mathbf{x})\}_{i=1}^{L}\in\mathbb{R}^{L\times d}$.

This section evaluates the performance of the EPM in typical agricultural environments. Two different settings, shown in Fig. 7 (a) and (b), were selected for the system evaluation, and their corresponding semantic maps are shown in Fig. 7(c) and (d), respectively. The results show the precise instance-level understanding of the EPM of the point cloud map in agricultural scenarios, with an overall recall and precision for detection of 0.95 and 1.0, respectively. The average position error and bounding box error are measured to be 0.06m and 0.02m respectively. Particularly in agricultural landscapes with row structures, our method demonstrates the ability to detect such features and generate a graph map for robot autonomy. Furthermore, the method is versatile and supports both online and offline modes depending on specific requirements. For this study, we perform the semantic extraction process offline.

This section evaluates the performance of the MPM for automated phenotyping by robots. We conducted tests in two different environments, shown in Fig. 8 and 9, corresponding to the environments shown in Fig. 7 (a) and (b), respectively.

In the experiment, target plants requiring phenotyping were randomly selected. The results show that the presented MPM achieves a high success rate and meets the robot’s requirements in both scenarios, as shown in (a) of Fig. 8 and 9. Compared to conventional $A^{*}$ or Dijkstra planners, which may generate inappropriate global trajectories leading to turning movements in narrow channels, causing local motion planning failures for car-like robots, our graph-based global planner performs better. Beyond the global planner, the developed local trajectory planner also shows strong robustness in field environments, as illustrated in (b) and (c) of Figs. 8 and 9, respectively. A significant factor affecting the performance of the local trajectory planner is the traversability analysis. Given the complex terrain in agricultural environments, traversability may be prone to overestimating or underestimating risk areas during planning, leading to motion planning failures. In the experiment, the robot was optimised with an additional solid-state lidar at the front to improve its perception under complex conditions. In addition, a replanning mechanism is employed for local motion planning, with an update frequency of 5HZ and a forward planning distance of 10m. Overall, our system demonstrates robust performance in the field with appropriate parameter tuning. The maximum speed during sampling is 0.2m/s, while the maximum speed for other conditions is 1m/s.

This section assesses the performance of IPM. We compare the in-situ phenotyping models of both scenarios by using handheld data acquisition and robot-automated acquisition. The quality of the results is evaluated by Peak Signal to Noise Ratio (PSNR), as detailed in Table 1.

In the initial experiments, both vanilla Instant-NGP and our model showed fast convergence, taking less than 2 minutes when using hand-collected samples, resulting in a PSNR above 24 dB and indicating high-quality results. However, when using robot-collected samples, Instant-NGP took over 4 minutes to converge and achieved a PSNR below 23 dB. In contrast, our method achieved a PSNR above 23.5 dB and converged in 3 minutes, demonstrating the effectiveness of occlusion regularisation in NeRF training with sparse-view inputs, effectively mitigating geometric estimation errors in the renderings. Moving on to the greenhouse experiments, the quality of modelling using both Instant-NGP and our model showed a decrease, with PSNR values of around 23 dB in approximately 3 minutes when using handheld collected samples. This reduction in quality is mainly due to the complex geometry of the canopy leading to occlusion. When modelling with robotically collected samples, both models showed a decrease in model quality along with a longer training convergence time. In comparison, our model performed better under these field conditions. The hierarchy maps for both scenarios are shown in Figures 10 and LABEL:fig:enter-label-2, revealing global structure-aware maps of farms and detailed models for individual plants. These results highlight the ability of the robotic system to achieve high-quality in-situ phenotyping in complex agricultural environments, demonstrating resilience despite certain limitations.