Integrating Real Tree Skeleton Reconstruction Based on Partial Computational Virtual Measurement (CVM) with Actual Forest Scenario Rendering: A Solid Step Forward for the Realization of the Digital Twins of Trees and Forests

Abstract

Digital twins of forests (trees) are computational virtual recreations of forests (trees) in which the entity distributions and physical processes in real-world forests (trees) are duplicated. It is expected that conventional forest science and management can be undertaken in a digital twin of forests (trees) if the recreation of a real-world forest (tree) has accurate and comprehensive enough information. However, due to the various differences between the current tree model and the real tree, these envisioned digital twins of the forests (trees) stay a theoretical concept. In this study, we developed a processing strategy that partially integrated computational virtual measurement (CVM) process into the tree modeling workflow. Owing to the feature of CVM, partial tree skeleton reconstruction procedures were considered to have higher mechanical objectivity compared to conventional mathematical modeling methods. The reason was that we developed a novel method called virtual diameter tape (VDT), which could provide a certain percentage of modeling elements using CVM. Technically, VDT was able to virtually measure diameters and spatial distribution of cross-sectional area of trees, including the basal area, from point clouds. VDT simulated the physical scenario of diameter tapes, observing point clouds of trees. Diameter and the cross-sectional area of stem and branches were obtained by two consecutive physical measurement processes, one in the forest sample site and another in the virtual space. At the same time, VDT obtained better or a similar accuracy compared to the mathematical methods, i.e., Hough transform-based methods, using the same data sets. The root-mean-square deviation (RMSE) of retrieval of diameter at breast height (DBH) using VDT was 1.02 cm, while DBH obtained from three conventional methods varied from 1.29 cm to 1.73 cm. Based on VDT measurement results, tree skeleton reconstruction and actual forest scenario rendering of our sample plots were further implemented. Beyond the visual consistency, we believe that our work might be a small and solid step in the technological evolution from tree models to the digital twin of forests (trees).

1. Introduction

1.1. Form Tree Models to Digital Twin of Forests

The construction of realistic three-dimensional (3D) tree models is becoming an attractive research focus [1,2,3]. Making a detailed tree model depends on two technological bases. The first and fundamental technology is LiDAR (light detection and ranging) scanning, which records the spatial information of trees at a very high sampling rate [4]. Compared to measures in sample plots using conventional instruments [5], the volume of raw measurement data increases exponentially, which supports tree modeling with high level of detail (LoD) [6]. The second dependency for making a detailed tree model is the computational method. In contrast to the computational virtual measurements (CVM) we employed in this work [7], most of the current methods are mathematically based. For example, quantitative structure models (QSMs) were developed by Raumonen et al. in 2013 [8]. The QSM’s approach considers that a tree can be decomposed into a huge number of little cylinders. A tree model is then generated as a collection of little cylinders. The QSM’s approach utilizes the classical technical route of calculus with great success. Since then, improved QSM methods such as TreeArchitecture [9], Simpleforest [10], and AdQSM [11], have been developed. The implementations of QSMs are booming [12,13,14]. In addition to QSM methods, conventional voxel tree modeling methods are continuously evolving [15]. However, it is more likely to be used as a tool for the estimation of tree and forest parameters [16]. One of the reasons, we think, is that voxel modeling methods are slightly weaker than QSMs in the visual presentation of the model.

Nowadays, increasing studies of the digital twin are reported in a variety of different fields [17,18]. However, the concept and clear definition of the digital twin of trees and the digital twin of forests are rarely discussed [19]. From our point of view, the deviation between the 3D tree model and the digital twin of the tree is one of the difficulties. What can be expected is that accurate 3D tree models could act as a precursor of digital twin of trees as well as virtual trees. In our previous study, we envisioned that a virtual tree should be an exact digital replica of the corresponding tree in reality [7,20]. However, due to various theoretical and morphological complications, e.g., the topological vagueness of tree structures reflected by remote sensing data, complex internal relevance [21,22,23,24] and the internal logic used in tree modeling [25,26], it seems technically difficult to construct an ideal virtual tree that perfectly maps a real tree, especially in the aspects of bole diameter estimation and tree skeleton reconstruction [27]. Meanwhile, a clear definition of virtual trees or digital twins of trees has not yet reached a unified consensus in society today. In light of recent trends [28,29,30], the term “realistic” is more likely to refer to visual consistency between real trees and virtual trees [31]. Any further analysis based on visually realistic tree models may raise concerns about objectivity. This is why we cannot regard the existing tree models as equivalent to virtual trees or digital twins of trees.

Equation (1) describes our view of the relationships between digital twins of trees and tree models. A tree model can be considered as a digital twin of tree if every predicational component made by the tree modeling algorithm can be verified. In such a situation, no false competent of trees can be added to the tree model. On this basis, a tree model with higher-level objectivity is constructed, which is an ideal state, i.e., the digital twin of a tree in reality.

$$digital twin of tree=\underset{n\to 0}{\lim } tree model\left(n\right)$$

(1)

where $n$ refers to the amount of unverifiable predictions made by algorithms during tree modeling.

Based on Equation (1), we further assume the technical description of the digital twin of forests. As shown in Equation (2), a digital twin of forest firstly consists of each digital twin of tree and other digital twins, e.g., the digital twin of undergrowth and soil surface. In addition to static expressions, it should also contain the realistic simulation of basic physical and physiological processes of trees and other objects in the forest. In this way, a digital twin of forest, as a toolbox, can support virtual experimentation and produce meaningful data that helps people to understand what could happen in reality.

$$digital twin of forest=\sum f\left(digital twin of tree\right)+\sum f\left(other digital twin\right)$$

(2)

where $f\left(x\right)$ refers to the function of the simulation of basic physical and physiological processes as toolboxes.

1.2. A Step Forward for the Realization of the Digital Twin of Trees

In our previous study, we developed the theory of computational virtual measurement (CVM) and an implementation of CVM, i.e., the virtual water displacement (VWD) method [7]. Compared to conventional mathematical procedures used in processing LiDAR point clouds, CVM simulates physical principles of measuring instruments and measuring procedures to measure point clouds of tree. In terms of algorithmic mechanisms, CVM is superior to mathematical procedures, e.g., tree modeling. This is because the CVM procedure is an observational process and does not make any predictions.

How to improve the mechanical objectivity of the tree modeling process? CVM provides a novel perspective of the utilization of algorithmic logic instead of the need for external reference. In this study, we developed virtual diameter tape (VDT) as a new implementation of CVM. VDT measured the point cloud of trees for the diameter and area of cross-sections of trees. Accordingly, the measurement results provided confident raw references for tree skeleton reconstruction. In addition to the tree skeleton reconstruction, actual forest scenario rendering was further implemented. Thus, compared to tree modeling using pure mathematical processes, this integrated workflow earned a small step toward better mechanical objectivity.

The mechanical objectivity of VDT depends on the precise duplication of the physical scenario of the diameter tape. The diameter tape was developed in earlier times and is still in use nowadays [32]. Explicit and simple physical processes are employed by these measuring instruments, while predictive mathematical processes are excluded. For example, a diameter tape makes physical contact with a tree stem in the position of the diameter at breast height (DBH) [33]. The measurement process is subdivided by us into several steps, and all of them are physical processes. In the beginning, the diameter tape encloses an arbitrary-shaped area, including the stem inside. Then, the diameter tape starts to shrink in space. In this process, the shape of the enclosed area gradually approaches the measurement object. In other words, areas that do not belong to the cross-sectional area are gradually excluded by the movement of the diameter tape. Theoretically, when all the redundant areas are excluded, we can obtain an exactly true DBH. In the final step, the termination signal of the measurement completion is also physical. First, diameter tape has no chance to cross the inner part of the stem. Second, the diameter tape is in contact with each point on the surface of the stem and cannot be further constricted by external forces from human operators. Despite the fact that human influence on physical processes has led to the need for non-predictive mathematical processes, such as the calculation of the mean value [33,34], we believe that the high objectivity of measurement using diameter tape is granted by its own processes. The key is that measurements using diameter tapes are purely physical processes, with no need for any mathematical preset of the measured objects.

As a widely used data source for tree modeling, LiDAR scanners use completely different physical mechanisms from that used by diameter tape to collect tree information. The sensor in LiDAR scanners records the distance between the sensor and the laser-irradiated spot [35]. Compared to diameter tapes, the beam emitted by the sensor travels in a straight line and cannot encircle the stem. Accordingly, the spatial information of the outline of the stem is unable to be determined in the raw measurement data of a single beam. We believe that this is the first decline in objectivity. As compensation for the loss of physical mechanism, the LiDAR sensor projects dense beams on the measuring target [36]. It uses rotating mirrors to deflect the beam to cover the necessary stem area [37]. The motion of the rotating mirror is controlled by a gear drive mechanism step by step [38]. Each shift of the rotating mirror leads to an independent measuring process. Point clouds are stored as a discrete data set containing the morphologies of each tree component, which are always affected by the instability of apparatus performance and disturbances raised by external environments. Restoring the continuous contours on the bole surface from the discrete data set requires the assistance of computer algorithms [39] and additional mathematical derivation [40]. Conversely, employing a diameter tape [41] causes a cumbersome experimental measurement process, especially for branches at the higher parts of forest canopy [42]. Moreover, manual measurement always reduces overall mechanical objectivity and results in unavoidable artificial errors.

The primary problem of the mathematical procedure is that there is a deviation between the shape of the natural stem and the human-defined stem using mathematical approaches. For example, commonly used approaches based on the concepts of circle fitting the approximate contour of the basal area of the stem [43,44], such as the Monte Carlo method [45], Hough transformation [46], random sample consensus (RANSAC) [47], and 3D modeling [48], are always restricted by the pre-definition of the regular shape of stem with detrimental factors of original shape restoration. Meanwhile, due to the lack of objectivity at the mechanical level, there is an obligation to provide justifications, which leads to additional workloads [49].

In the evolution from tree models to digital twins of trees, while waiting for the overall evolution of current LiDAR scanning mechanisms and tree modeling methods, we could try different technological routes. With this work, we employed CVM as an efficient tool for measuring the diameter of the bole and branches. Meanwhile, the phenotypic traits of the tree skeleton can be graphically portrayed in pursuit of forest plot scenario reconstruction with high fidelity and enhanced immersion in a parallel spatiotemporal continuum. Meanwhile, trees in different growth stages and development situations were recorded in detail and realistically visualized, which affords historic characterizations and comprehensive information to guide current forest cultivation practices and silvicultural managements.

2. Material and Methods

2.1. Study Area and Field Measurements

The study area was located in the Hainan Province, China, between (109°43′E, 19°28′N) and (109°51′E, 19°38′N) with an area of 5000 ha. The main landform type is hilly plateau, and the average altitude is 188 m. The climate type is tropical savanna climate. Since the 1950s, rubber trees (Hevea brasiliensis) have been planted extensively for commercial purposes in this area. Presently, rubber plantations cover ca. 3000 ha. Due to the complicated ownership of property, rubber trees of each age group are mixed together and have not developed a clear spatial pattern if we observe from the sky.

Three sample plots of rubber trees were measured as shown in Figure 1d. The average tree ages in those sample plots were five years, ten years, and twenty years old, respectively. As shown in Figure 1a, each sample plot consisted of ten trees surrounded by five LiDAR scan positions. All rubber trees in sample plots were planted in a narrow-wide-row planting pattern (2 m by 4 m). We used Leica ScanStation C10 for scanning. The primary capabilities of Leica ScanStation C10 were listed as the following: accuracy of single measurement, 6 mm (position), 4 mm (distance); target acquisition, 2 mm std. deviation; range, 300 m at 90%; 134 m at 18% albedo (minimum range 0.1 m); scan rate, up to 50,000 points/sec, maximum instantaneous rate. Figure 1b shows the obtained point clouds from the five-year-old rubber tree sample plot. Figure 1c shows a photo of when a LiDAR scan was performed.

In addition to LiDAR scanning, conventional field measurements, collecting three tree variables, were performed. These measurements provided a ground truth reference for evaluating new methods developed in this study. First of all, diameter of breast height (DBH) was measured using plastic diameter tape for each tree. Then, the stem diameter at an arbitrary height along a tree bole was measured using diameter tape or a novel diameter clipper [50].

The third variable was the diameter of first-order branches. This work was labor-intensive and could not be applied to each tree in the three sample plots. Therefore, three adjacent trees in each sample plot were picked and then measured. The criterion for the selection of the sample trees was point cloud completeness. Therefore, these measurements were performed after the data preprocessing. We inspected the aligned point clouds of each sample plot to find trees with relatively minor defects. Then, crews in our survey team built and climbed scaffoldings and performed measurements of the branch base diameters for three trees in each sample plot.

2.2. Data Preprocessing

Point clouds of ten trees from each sample plot, for a total of 30 trees, were used to carry out VDT measurement and create forest scenario rendering. To prepare the qualified inputs, we first performed a data preprocessing work. As shown in Figure 2, this work consisted of two sections. The first section was the conventional data preprocessing procedure for the ground LiDAR scanning. The second section was the amelioration procedure for the new method developed in this study.

2.2.1. Conventional Data Preprocessing

The technical steps in data preprocessing included registration of point clouds from scanning stations, normalization of point cloud height, single tree segmentation, wood-leaf segmentation, and making slices of tree components. Depending on conditions, each slice could be identified as an ideal point cloud [7,51], or a non-ideal point cloud. For the non-ideal point clouds, a mathematical procedure was further applied (in Section 2.2.2). After that, all point clouds were regarded as ideal point clouds.

For registration, Cyclone (Leica Geosystem, Switzerland) software was used to align multi-station point cloud data. The root-mean-square error (RMSE) varied from 4.1 mm to 8.9 mm. In addition to registration, we applied manual deletion of abnormal objects in the aligned point cloud, such as people being scanned and rubber collectors attached to trees. After that, the denoising process was applied using various commonly used filters [3,52] to eliminate the isolated points with the aim of reserving compact tree point clouds. We used a progressive morphological filter to classify ground and non-ground points [53], and the ground points were removed. Consequently, we applied the previously developed individual tree segmentation method to separate point clouds to individual trees [2], and a wood-leaf segmentation algorithm was adopted here for bole and branch component extraction [52]. This single tree point cloud at the whole tree scale was used to verify the growth direction of the stem, and provided the projection coefficient for the following step [54]. A cross-section of the tree skeleton was sliced for stem diameter assessment. Slicing was performed every 0.5 m along the stems. There were thirty trees in total, ten trees in each age group with the acquired horizontal slices obtained. These prepared inputs for the VDT method could be regarded as the tree skeleton presented in the form of point clouds. From an informatics perspective, points are disordered data. The rest of the process of this work was designed to establish order out of disorder.

2.2.2. Amelioration of the Quality of the Wood Points with Gap Filling

A significant proportion of the raw data was made up of non-ideal point clouds, which caused errors in VDT measurement. The function of the amelioration was to make synthetic ideal point clouds that serve as qualifying inputs for VDT. The ideal point cloud is a concept we developed in our previous work [7]. It refers to a point cloud with a high point cloud density, and the gap between each point is generally even. If this condition is not satisfied, then a point cloud is a non-ideal point cloud.

The ideal situation of the collected point clouds should uniformly cover the entire surface of the tree’s woody facets. However, in the interlocked forest plot, mutual occlusion effects by vegetative elements make an occlusion-free position impossible to arrange, as the terrestrial laser scanners for acquainting the data for all the details of the tree body is impracticable, resulting in many tree branches and twigs not being fully hit by laser beams and the gaps emerging in the extracted slice contours regarding local point cloud deficiency. Hence, filling the gaps based on ellipse-fitting methods was deployed herein.

Figure 3 shows the workflow of the amelioration process. When the LiDAR scan failed to capture an ideal point cloud of a cross-section of stem (in Figure 3a), a non-ideal point cloud would be generated. In Figure 3b, we used purple points to represent the recorded points and white points to represent the missed points. To prevent the VDT detector from moving within the internal area in the woods, the missing area (white points) had to be compensated for. In Figure 3c, a fitting circle (red) was employed to enclose the cross-section area. In the next step, a synthetic ideal point cloud was formed through the steps shown in Figure 3d–g. The synthetic ideal point cloud consisted of (purple) points of raw data and (red) synthetic boundary lines. As shown in Figure 3g, a synthetic ideal point cloud preserves the original information from the LiDAR scan, which effectively maintains partial objectivity compared to a pure mathematical fitting as shown in Figure 3c. Finally, a VDT measurement was proceeded on this synthetic ideal point cloud, as shown in Figure 3h. A detailed description of the VDT measurement is Figure 4.

The key process of the amelioration was making the optimal ellipse fitting in Figure 3c, containing three steps: step 1, the determination of the long axis of the fitted ellipse; step 2, the determination of the short axis of the fitted ellipse; and step 3, the determination of the rotation angle of the fitted ellipse. Key mathematical approaches were listed as the following:

Step 1: determination of the long axis of the fitted ellipse.

The input was a cross-section of a tree skeleton. The 3D point cloud was converted to a two-dimensional (2D) point cloud by vertical projection onto the X-Y plane. Then, one or a few incomplete ellipse shapes can be seen in the X-Y plane. $P_j$ is used to represent the point set of each individual ellipse, i.e., a slice for each individual trunk or branch. In $P_j$, $p_i$ refers to a point in the point cloud set and $p_i\in P_j$. Using the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm [55], we finally separated the wood point clouds belonging to a fixed height interval into different components $P_j$.

As shown in Figure 4a, an unenclosed ellipse E (in Equation (3)) consisting of point clouds was discernible.

$$\frac{{\left(x-y\right)}^2}{a^2}+\frac{{\left(x\mathit{sin}\theta +y\mathit{cos}\theta \right)}^2}{b^2}=1$$

(3)

where x and y refer to the horizontal and vertical coordinates of the point in the 2D plane, respectively; a and b refer to the lengths of the semi-major and semi-minor axes of the ellipse separately.

The shading effect could be the reason for this. In order to determine the long axis of the fitted ellipse, a straight line $l_{p_i{p}_{i+1}}$ containing any two points $p_i\left(x_i,y_i\right)$ and $p_{i+1}\left(x_{i+1},y_{i+1}\right)$ was randomly selected. There was no order relationship between $p_i$ and $p_{i+1}$, which simply represent any two points in $P_j$. The Euclidean distance between $p_i$ and $p_{i+1}$ was calculated using Equation (4).

$$\left|p_i{p}_{i+1}\right|=\sqrt{{\left(x_i-x_{i+1}\right)}^2+{\left(y_i-y_{i+1}\right)}^2}$$

(4)

where $p_i$ and $p_{i+1}$ are two randomly selected projection points; x and y refer to the coordinates of the two points.

After iterative selection of any tow points belonging to $P_j$, the maximum value of $\left|p_i{p}_{i+1}\right|$ was determined as the potential length of the long axis of the predicted ellipse. It was marked as $l_{p_i{p}_{i+1}}$ ($\left|l_{p_i{p}_{i+1}}\right|=2a$ in standard pole-polar equation of ellipse). The midpoint of the long axis $l_{p_i{p}_{i+1}}$ was marked as $p_{mid}\left(x_{mid},y_{mid}\right)$, $x_{mid}=\frac{x_i+x_{i+1}}{2}$ $y_{mid}=\frac{y_i+y_{i+1}}{2}$. The slope of $l_{p_i{p}_{i+1}}$ was expressed as $\tau =\frac{y_{i+1}-y_i}{x_{i+1}-x_i}$. The equation $l_{p_i{p}_{i+1}}$ was transformed into Equation (5). Meanwhile, the $p_{mid}\left(x_{mid},y_{mid}\right)$ was taken as the centroid of the predicted ellipse.

$$y-y_i=\tau \left(x-x_i\right)$$

(5)

where $x_i$ and $y_i$ refer to the coordinates of $p_i$; $x and y$ refer to coordinates of the point on the straight line $l_{p_i{p}_{i+1}}$; $\tau$ refers to the slope of $l_{p_i{p}_{i+1}}$.

Step 2: determination of the short axis of the fitted ellipse.

The short axis $p_k{p}_{k+1}$ is perpendicular to the long axis $l_{p_i{p}_{i+1}}$ of the ellipse and intersected at center point $p_{mid}$ shown in Figure 4b. Then, the length of short axis was denoted as Equation (6).

$$\left|p_k{p}_{k+1}\right|=\sqrt{{\left(x_k-x_{k+1}\right)}^2+{\left(y_k-y_{k+1}\right)}^2}=2\mathrm{b}$$

(6)

where $p_k$ and $p_{k+1}$ refer to endpoints of the short axis; $x and y$ refer to coordinates of the endpoints; $2\mathrm{b}$ refers to the length of short axis of the fitted ellipse.

The predicted short axis could only have one endpoint, i.e., $p_{k+1}$ does not exist due to the vacancy of point clouds. In this condition, $p_{k+1}^{\prime }$ was used to represent the missing endpoint of the short axis according to the length of $\left|p_{k+1}^{\prime }{p}_{mid}\right|=\mathrm{b}$ and the direction perpendicular to the long axis $l_{p_i{p}_{i+1}}$.

Step 3: determination of the rotation angle of the fitted ellipse.

Before the process of the finding of the optimal value of $\theta$, we performed a coordinate transformation. The coordinate of $p_c$ was designed as the coordinate centroid in the 2D Euclidean space. Then, we gradually rotated the ellipse at a pace of $\frac{\pi }{180}$ counterclockwise starting from $-\frac{\pi }{2}$ to $\frac{\pi }{2}$, yielding the translated equation of the fitted ellipse as $E_r$ in Equation (7).

$$E_r: \frac{{\left(x\mathit{cos}\theta -y\mathit{sin}\theta \right)}^2}{a^2}+\frac{{\left(x\mathit{sin}\theta +y\mathit{cos}\theta \right)}^2}{b^2}=1$$

(7)

where x and y refer to the coordinates of the points on the ellipse; a and b refer to the lengths of the semi-major and semi-minor axes of the ellipse, respectively; $\theta$ refers to the angle between the long axis of the ellipse and the horizontal axis.

Then, a variable point on the elliptic boundary was defined as $q\left(x_q,y_q\right)$, and the distance between $p_i$ and $q$ denoted as $\left|p_{i}q\right|=\sqrt{{\left(x_i-x_q\right)}^2+{\left(y_i-y_q\right)}^2}$. Here, we utilized the Lagrange multiplier method to calculate the minima distance $d_i^{min}$ between each point $p_i\in P_j$ and $q$, and then the following Lagrange function was constructed as:

$$F\left(p_i,q,\mu \right)=f\left(p_i,q\right)+\mu ·\phi \left(p_i\right)$$

(8)

where $p_i$ refers to a random projection point; $q$ refers to a variable point on the elliptic boundary; $\mu$ refers to a variable called the “Lagrange multiplier”; $\phi$ refers to equation of an ellipse; $F\left(x\right)$ refers to the Lagrange function; f(x) refers to the Euclidean distance between $p_i$ and $q$.

Where $f\left(p_i,q\right)={\left(x_i-x_q\right)}^2+{\left(y_i-y_q\right)}^2$ and $\phi \left(p_i\right)=\frac{{\left(x_i\mathit{cos}\theta -y_i\mathit{sin}\theta \right)}^2}{a^2}+\frac{{\left(x_i\mathit{sin}\theta +y_i\mathit{cos}\theta \right)}^2}{b^2}$. Solving Equation (8) satisfies the following partial derivative equal to 0 under the conditions of a given value for $\theta$, i.e., $F_{x_i}^{\prime }=0;F_{y_i}^{\prime }=0;F_{\mu }^{\prime }=0$. Then, the distance $dis{t}_{sum}$ between all the projected point sets of any branch slice $P_j$ and the fitted rotation ellipse $E_r$ at $\theta$ angle can be expressed as follows: $dis{t}_{sum}={\displaystyle \sum }_{i=1}^N\left(d_i^{min}\right)$, where $N$ represents the total number of points for a branch slice. For all the points of a branch slice, pursuing the rotation angle $\theta$ in correspondence to the minimal value of $dis{t}_{sum}$ is the aim of deriving the optimized solution of the rotated ellipse equation, i.e., determination of the rotation angle $\theta$.

2.3. Analysis of Physical Scenario of Diameter Tapes

The successful development of a virtual measurement instrument (VMI) depends on a detailed analysis of the physical scenario of its corresponding real measurement instrument (RMI) [7]. In this study, the principle of the physical mechanism of diameter tape was elucidated, acting as a virtual diameter tape to determine the diameter, perimeter of each branch, and bole in pursuit of the spatial information collection of tree woody facets.

The physical scenario of diameter tapes when measuring DBH is shown in Figure 5a–d. Before the measurement was applied, the space around the stem was physically divided into two regions in Figure 5a. There was a basal area (internal area) and external area. The first encirclement of the stem was achieved by the diameter tape shortly after the start of the measurement process in Figure 5b. Accordingly, an enclosed area was produced by the diameter tape. As shown in Figure 5c, with the continuous shrinking of the diameter tape, the enclosed area was approached the basal area. This process was a non-discrete process. At each moment, the instantaneous closed area could be considered as a measure of the basal area. How to determine the best measure among an infinite number of measures? On the face of it, there was no doubt that the smallest value should be chosen. The objectivity of this choice was supported by the sample and clear termination condition of the measuring process. As shown in Figure 5d, the diameter tape (considered a rigid body) reached a steady state when it had no room for further contraction. This termination condition was independent of the outcome of the measuring process. The exact same physics was utilized each time to measure DBH using the diameter tape. In our opinion, this provides the basis for making the diameter tape a standard measurement tool for forest sample sites worldwide. The raw measurement datum was the tree parameter without the need of additional mathematical processes.

2.4. Design of Physical Scenario of Virtual Diameter Tape (VDT)

The first task in designing the physical scenario of the virtual diameter tape was to consider the role of point clouds from another perspective. As shown in Figure 5d, the outer contour line of the stem was detected by the diameter tape, as it did not allow external objects to penetrate the stem. Learned from this phenomenon, the points in the point clouds were considered solid objects in virtual space. At the single-point level, external objects in the virtual space cannot pass through over a single point. Accordingly, the set of point clouds would react to the (virtual) physical contact in a similar way as the stem would if it occurred in reality. Therefore, the points in the point cloud were no longer considered data points; they were components composed of virtual stems. With this foundation, the conventional mathematical estimation of the basal area was instead performed by computational virtual measurement on the virtual stem.

According to our previous knowledge [7], an exact duplication of the physical scenario of an RMI may exceed the abilities of current (desktop) computer resources. Consequently, additional predictive mathematical procedures may be introduced into the CVM process. These procedures are expected to have a negative impact on the mechanical objectivity of VMI. Therefore, redundant simulations of physical behavior, for example realistic forms of diameter tape, are excluded from this study.

As shown in Figure 5b, the key physical behavior of the diameter tape was that it could occupy an exclusive area regardless of its morphological and positional distinctions. In accordance with this, the VDT was designed as a spot or circle on a 2D plane, as shown in Figure 5f. Both the VDT and the stem disk filled an area in this plane entirely. After that, instead of the directional motion of the diameter tape (in Figure 5c), the motion pattern of the VDT was determined using random motion. First of all, a VDT was released outside the stem disk at Position (i) in Figure 5g. Then, it moved straight and made physical contact with the point cloud. At that moment, a physics engine, such as PhysX, was responsible for securing non-crossover between VDTs and points. Therefore, there was no chance of a VDT moving within the stem disc, as shown in Position (i’) in Figure 5g. Simultaneously, the direction of motion of the VDT was forced to change to a new random heading, such as at Position (ii) in Figure 5g. When the VDT had been functioning for a long enough period of time, as shown in Figure 5h in green, all areas that were not part of the stem disk could be accessed by the VDT as the same physical scenario in Figure 5d.

2.5. VDT Implementation

Unity (Unity Technologies, San Francisco, CA, USA) software was used to implement the VDT. Unity provided a basic 2D virtual space that could contain VDT and point clouds of stems. The simulation of basic physical laws (Newtonian mechanics) was provided by PhysX (NVIDIA Corporation, Santa Clara, CA, USA), which was integrated with Unity. VDT consists of the proper set of physical laws owing to the design of the physical scenario of VDT.

Compared to conventional mathematical-based procedures, the implementation of VDT was so simple that it nominally precluded the use of advanced knowledge in mathematics. The virtual space was supported by the default setting when creating a 2D project. In the created virtual space, points, VDT, and boundaries of virtual space were designated as Rigidbody. After that, the physical scenario was inputted in this virtual space. From a user’s perspective, VDT also has a similar operating logic to RMI, which initiated the measurement process by pushing a button and waiting for the results.

The pseudo-code of VDT is shown in Algorithm 1. The development environment was as the following: hardware: CPU: Intel i5-11400F at 2.6–4.4 GHz; RAM: 40 GB DDR4 2933 MHz; graphics card: NVIDIA GeForce GTX 1650S with RAM 4 GB; software: Unity 5.7.0, PhysX SDK 3.4 (for implementation of VDT).

Algorithm 1: VDT measurement process

while (termination_condiation is false)

VDT_Move(direction)

SaveFootprint();

if CollisionDetected();

direction.new(random);

end

end

End

2.6. Retrieval of DBH from Raw VDT Measurement Outcome

Figure 6a shows a raw VDT measurement outcome. It was not necessary to wait for full access of VDT detector in the entire area. When the size of the central area varied at a low rate, reached a potential termination condition for VDT measurement. Then, simple image segmentation methods were applied to the raw outcome [56]. As shown in Figure 6a, the perimeter, as well as the DBH of a cross-section of trees were extracted.

2.7. 3D Tree Skeleton Reconstruction and Actual Forest Scenario Rendering

We used Speedtree (IDV, South Carolina, CA, USA) as a tree modeling tool [57]. Compared to quantitative structure model (QSM) tree modeling methods [8], Speedtree was characterized by not applying algorithm logic to generate false predictions of tree components. The skeletons of trees were constructed based on the VDT measurement results. After that, actual forest scenario rendering was applied. Eventually, visual realistic tree models were completed with the help of VDT.

2.8. Reference Methods

We prepared three conventional mathematically based methods to compare with the VDT measurement. The methods were circle fitting [43], circle fitting (Hough) [46], and cylinder fitting [48]. Circle fitting (Hough) is short for circle fitting with Hough transformation. Firstly, we applied all three methods to the ideal point cloud dataset. Circle fitting (Hough) showed the best performance. Therefore, we further applied circle fitting (Hough) to the dataset of point clouds of medium and low completeness to compare with the VDT measurement. In this step, it was not necessary to apply the other two methods. The reason was that datasets of ideal point clouds and point clouds of medium and low completeness were subsets of the same dataset. Consistency of algorithmic performance was expected.

3. Results and Discussion

3.1. VDT Measurements on Ideal Point Clouds

Point clouds at DBH could be considered ideal point clouds in this work, as confirmed by our visual inspection. Therefore, the amelioration process was not applied. VDT had the best performance in this data set. Figure 7 shows the result of DBH extraction using four methods on 30 trees. Three of them, i.e., VDT, circle fitting, and circle fitting (Hough), show similar performance in terms of RMSE. In the case of the same sample size, the VDT method has a higher coefficient of determination value and lower RMSE in comparison to all other methods. The performance of cylinder fitting is among the worst.

3.2. VDT Measurements on Point Clouds of Medium Completeness

Point clouds of stem diameters along the stem curve consisted of both ideal and non-ideal point clouds. Under this condition, an amelioration step was applied. Then, VDT could effectively improve the accuracy compared to pure mathematical-based procedures. We believed that this was the power of the CVM procedure, even if it was only a partial CVM. Figure 8 shows the result of the extraction of stem diameters using two methods. One hundred and twenty stem diameters at different heights were produced using the circle fitting (Hough) method in Figure 8b. Based on this, VDT measurements were applied, as shown in Figure 8a.

3.3. VDT Measurements on Point Clouds of Low Completeness

Compared to point clouds in stem diameters, the point cloud completeness in the first order branches degraded dramatically. Figure 9 shows the result of extraction of the base diameter of the first order branches from nine trees in three sample plots. Due to the shading effect, a large percentage of slices were poor in completeness. In this case, the natural boundary that consisted of the original points took a small percentage. From the results, the VDT showed the adaptation of this poor-completeness data, and provided possibilities to make improvements. However, compared to conventional mathematical-based methods, VDT did not show improvements because the mathematical process was a dominant part of the process of estimating branch diameter. Accordingly, the partial CVM procedure had a limited contribution to the overall procedures.

Compared to three conventional mathematically based methods, VDT shows the best performance. From our point of view, the delicate mechanism of VDT made a significant contribution. Furthermore, the characteristics of the ideal point cloud had a favorable impact, widening the performance gap between VDT and conventional mathematical methods. The reason was that each cross-section of trees was not a perfect geometric shape. Ideal point clouds possessed various irregular spatial information on the cross-section surface. In this condition, irregular spatial information would be more likely to be ignored by the regular presuppositions of conventional mathematical methods. On the contrary, the mechanism of VDT was insensitive to specific shapes. However, there had been no previous studies of DBH retrieval using CVM due to two reasons. The first one was that CVM is a theory was recently proposed in 2021. The second reason was that it would not be easy to collect an ideal point cloud that allows implementations of CVM, e.g., VDT, to be executed. Owing to this fact, we designed the amelioration process to make the VDT method feasible to the current LiDAR data.

Table 1. Comparison between virtual diameter tape and circle fitting with Hough transformation.

Name	Virtual Diameter Tape	Circle Fitting (Hough)
Methodology	CVM	Mathematical modeling
Prior knowledge	No need	Cross-section area is round-shaped
Adoption of various shapes	Yes	No
Role of data	Objects in virtual space	Elements to be processed
Processing logic	Simulation of physical mechanism of diameter tape	Making prediction using presets
Raw result	Cross-section area	Model of cross-section area
Number of result (s)	Infinite	1
Final result	Perimeter (DBH)	Perimeter (DBH)
Need for validation	No (using ideal point clouds)/Yes	Yes
R2 (in IDL 1/Mid 2/Low 3)	0.97/0.96/0.88 4	0.94/0.94/0.85 4
RMSE(cm) (in IDL 1/Mid 2/Low 3)	1.02/1.30/0.56 4	1.29/1.56/0.55 4

In the table: ¹ IDL refers to the inputted point clouds are the ideal point cloud; ² Mid refers to the inputted point clouds are the point cloud of medium completeness; ³ Low refers to the inputted point clouds are the point cloud of low completeness; ⁴ the mean value from nine trees in three sample plots; CVM refers to computational virtual measurement; R² refers to determination coefficient; RMSE refers to root-mean-square deviation; DBH refers to diameter at breast height.

shows the comparison between the VDT method and the circle fitting (Hough) method on three different quality point clouds. The fundamental distinction of those two methods was in the methodology, which led to several differences. The first difference was the request for prior knowledge. VDT shared a similar mechanism on diameter tape and didn’t need to know the shapes of its measuring target. Meanwhile, mathematical approaches, such as the circle fitting (Hough) method, constrains the cross-section of trees to the established regular shape, which did not preserve the spontaneous morphological shape during the natural growth stage and resulted in inevitable deviations for the measurement of the woody horizontal cross-sections.

It is usual to achieve accurate estimation using conventional mathematical approaches for DBH. For example, Maas reported that DBH could be determined with a precision of about 1.8 cm in 2008 [58]; RMSE of DBH varied between 1.17 cm and 1.22 cm, reported by Liu et al. in 2018 [59]; the mean average deviation of DBH was ±2.90 cm, reported by Ritter [60] in 2017; the best RMSE of DBH in a work was 1.21 cm, as reported by Vandendaele et al. in 2022 [61]. In comparison with our method, Table 1 shows our method exceeds roughly one to three percent of the determination coefficient R² than the circle fitting (Hough) algorithm, as a representative of mathematical methods, with a decrement (approximately 0.27 cm in the ideal point cloud and 0.26 cm in point clouds of medium completeness) in RMSE for our method. Due to the diversity of forests and ecological environments, circle fitting algorithms always greatly depend on the parameter tuning operations and iteratively seeking optimized solutions. Conversely, the VDT method did not make these demands. Free from stem shape variations, VDT has the capability to measure any shape of branch cross-sections. We expect that the proposed method has strengthened universality with improved contour portrayal accuracy. This was also the initial ideological connotation when we conceived the CVM theory.

3.4. Rubber Tree Model Reconstruction Based on the Derived Growth Parameters

Our modeling results are based on the retrieved growth parameters and include the comparison with the original point clouds from various perspectives, as shown in Figure 10 and Figure 11. Statical information of tree models is shown in

Table 2. The growth variables retrieved from the modeling process and scanned points.

Age of Trees	Tree Height/m	Crown Width/m	Crown Volumes/m3	Average Length of Primary Branches of Rubber Tree/m	Average Included Angles between Stem and Branches/°	Tilt Angles of Stem/°
Five years	7.46 ± 0.61	3.02 ± 0.71 (N-S) 3.09 ± 0.97 (E-W)	69.63	1.840 ± 0.44	39.53 ± 4.84	3.92 ± 3.33
Ten years	9.18 ± 0.48	2.89 ± 0.95 (N-S) 4.21 ± 1.6 (E-W)	85.39	2.49 ± 0.68	41.11 ± 10.66	7.86 ± 7.67
Twenty years	14.47 ± 3.99	2.66 ± 1.55 (N-S) 5.22 ± 2.59 (E-W)	146.52	2.57 ± 2.33	45.91 ± 18.75	14.22 ± 13.9

. In addition to tree models, the environment was also added to the 3D scenes of three rubber tree sample plots. As shown in Figure 10 and Figure 11, morning and night scene effects were set to add ornamental completeness to the whole model. We regarded them as visual virtual sample plots and as a prelude to a physically realistic sample plot.

The partial modeling parameters used in the modeling process were acquired by measuring point clouds, which are listed in Table 2. It could be seen that the average lengths of primary branches of rubber trees increased with age of the trees, which is consistent with our own experience and related study [62]. Besides the average lengths, the range of average branch angle also varied with the age of the trees. Usually, trees accelerate their growth as they get older and bigger [63]. Some studies have overturned the assumption that old trees are less productive [64]. Meanwhile, height growth of old trees virtually stops, but bulks up like a bodybuilder. Twenty-year-old rubber tree bodies manifest a larger included angle between stem and branches, which indicates a bigger tree crown that provides more space occupation for foliage clump development. Moreover, as trees get older, the variations in the included angle also increase, which suggests that rubber trees are subject to intraspecific competition resulting in growth space exploration in pursuit of optimum physiological conditions.

In narrow-wide-row planting arrangements, wide rows with vacant space provide enough light illumination and growing space, which attracts branch stretching for optimizing photosynthesis conduction [65]. Hence, a reasonable explanation of rubber trees that continuously extended to the space on both lateral sides (where the wide row exists) with increasing age is reflected in Table 2, especially in the middle and upper layers of the forest skeleton. It can be clearly seen in Table 2 that the average tilt angle climbed up to 14.22° for 20-year-old rubber trees from 3.92° for 5-year-old rubber trees. This makes the 20-year-old rubber trees in Figure 10e,f have relatively more crooked bodies. Meanwhile, the computer graphics technology employed in this paper exactly shows the portrayal of the modeling results for rubber trees of various ages. Furthermore, tree height, branching angle, branch length, etc., can also be precisely depicted through this technology to achieve the best simulation effect of the models.

3.5. Sources of Error in CVM Procedure

The main error sources in this study were divided into two parts. The first part was in the non-CVM procedure, e.g., LiDAR scanning in a forest field [66] and the data pre-processing [67]. Those sources of error were comprehensively discussed in various studies [68,69,70]. Therefore, in this section, we focused on the second part of the error sources, which were presented in the CVM process.

The systematic error in the CVM process came from the deviation between the physics in reality and the simulation of physics. Usually, a third-party physics engine is employed, such as PhysX, in a CVM implementation to provide a simulation of basic physical laws. Those virtual physical laws are the basis of a virtual space. However, we could not find the specification of simulation accuracy provided by the manufacturer of the physics engine. Instead, Erez et al. made a comparison of simulation accuracy between five popular physic engines [71]. The result showed that the accuracy of the simulation of movement varied from 10⁰ to 10⁻¹³ (in different physical units, e.g., N·m·s and meter). Therefore, with an appropriate setting in a physic engine, the simulation accuracy could be orders of magnitude higher than the error in the data collected by measuring instruments in reality. Under these circumstances, we believed that the error in the simulation of basic physical laws could be neglected.

The random error was related to a specific physical scenario in an implementation CVM, e.g., the physical scenario of VDT in this study. In general, the diameter of the VDT detector was related to the average gap between points, and usually the diameter of the VDT detector was set to a number of times larger than the average gap. This setting prevented the VDT detector from moving inside the stem. Points in a point cloud not homogeneously distributed. Therefore, as shown in Figure 12a,b, the accessed area of the VDT detector was sensitive to the different gaps. To minimize these variations, a simple strategy could be applied, which was to increase the ratio between VDT detector diameter and the average gap. However, an unrestricted increment in this ratio would lead to the missing of the detection of area. Figure 12c,d simulates this situation.

3.6. Forward for the Realization of the Digital Twin of Trees

Compared to current methods that could provide ground truth reference, the advantage of VDT was exhibited with a higher degree of objectivity. To be specific, the need for external reference was no longer necessary. In other words, the auxiliary measurements performed by instruments in nature were replaced by computational virtual measurements (CVM) using VDT in virtual reality. Accordingly, ground truth reference for tree modeling could be provided by VDT. By reviewing the results of VDT in high-completeness point clouds with other conventional methods, the accuracy of VDT measurement was obviously higher than others when the same data set was processed. VDT also had better performance in hybrid-raw data. These facts confirmed that VDT was not only theoretically sound but also technically practical.

In our previous study [7,20], we proposed the theoretical description of virtual trees and CVM. Due to the current limitations from both instruments and the (civilian) computer environment, practical utilization of a complete CVM procedure on a full tree is not likely to become a reality anytime soon. In this study, we tried our best to integrate partial CVM in a conventional tree modeling process to provide key reference points without external references. This integration could be a small step on the evolutionary path from tree models to virtual trees or digital twins of trees by providing a higher level of mechanism objectivity. On this basis, we foresee a completely automated approach to transforming actual trees into virtual twins. Then, the digital twin can be subjected to simulated experiments. For the time being, we genuinely think that CVM could still serve as a future method. No matter what, we have faith that every 1000-mile-long journey begins with a single step.

3.7. Contribution of Digital Twin Technology to Rubber Tree Management

The virtual replication of rubber trees of different ages based on digital twin technology provided a new pattern for tree growth monitoring and virtual cultivation practices for tree plantation. Moreover, the presented work propagates the immersion, fidelity, interoperability, and development of silvicultural management for various rubber tree cultivars. We pursued a realistic representation of the study tree species in two aspects, i.e., phenotypic trait similarity and physical property consistency with guiding outcomes for created counterparts. In terms of phenotypic trait similarity, tree organ adhesion and partial crown overlapping somewhat exacerbate the meticulous portrayal of tree morphogenesis. With regard to physical property consistency, VDT included (partial) CVM procedures that could provide (partial) true values as key control points attached to visual models. In this way, a considerable percentage of real physical information was integrated into the visual models for enhancement of the fidelity of the reproduced forest scenarios. However, more physiological properties of wood materials of rubber trees, such as moisture, shrinkage, internal stresses, swelling, cracking, density, and sound/electrothermal conductivity, could not be effectively reflected by our models. This aspect could be improved further, as it would be especially beneficial to the exhibition of the interaction between the intermittence of natural wind and the yielded mechanical loading on rubber tree branches. Meanwhile, the requirement of rubber tree modeling employing digital twin techniques is recording the target tree morphologies spatiotemporally, seasonally, and meteorologically, rendering not only computer graphics-driven plant phenotypic trait characterization but also auxiliary research into the natural environment–plant biophysiological process. Furthermore, deep learning techniques [72] play a vital role in calculating actual forest growth parameters from various remote sensing data, which expedites interactive and immersive modeling procedures and enables users to readily covey their cognition and creation by integration into the mapped tree structures [73]. Since our laptop was not equipped with a graphics processing unit (GPU) for accelerating the creation of 3D forest models and rendering of complex scene lighting, the average tree modeling took approximately 0.3 Hours. In the future, with GPU assistance for graphical rendering, the time complexity and computational expenditure will decrease for better visualization and modification of real-time-generated branch and tree crown results.

4. Conclusions

A physical-based method, VDT, was successfully developed in this study. Technically, the VDT measures the diameter of the stems/branches as well as the geometry of their contour lines. VDT has been shown to have not only a theoretically better degree of objectivity but also higher accuracy when compared to conventional mathematical-based methods. In ideal point clouds, the RMSE of DBH was reported at 1.02 cm using VDT, while RMSEs were reported from 1.29 cm to 1.73 cm using conventional methods. In medium- and low-completeness point clouds, VDT shows good compatibility with the current point clouds, which indicates its compatibility with current LiDAR data. Based on that, real tree skeleton reconstruction and actual forest scenario rendering were further implemented. In addition to visual consistency, we believe that partial physical consistency was achieved.

This research envisioned a technical path for the realization of digital twins of trees. It featured the use of natural law to enhance the mechanical objectivity of the tree-modeling process. It featured the use of natural laws to improve the mechanism objectivity during the modeling of the tree. The current major limitation of VDT is the collection of ideal point clouds, which is only a small percentage of the raw data of LiDAR scanning. It would be expected the appearance of a new disruptive technology that can produce a high percentage of ideal point clouds.

Based on our experience, we can expect that a CVM procedure that includes various realistic physical simulations would yield a qualified digital twin of the tree in the future. Meanwhile, we anticipate the emergence of a field of measuring instruments using new mechanisms that will overcome the existing inadequacies found in LiDAR scanners. A joint effort between virtual and real measurement instruments can lead to the realization of digital twins of trees and forests.